# Coexistence of density wave and superfluid order in a dipolar Fermi gas

###### Abstract

We analyse the coexistence of superfluid and density wave (stripe) order in a quasi-two-dimensional gas of dipolar fermions aligned by an external field. Remarkably, the anisotropic nature of the dipolar interaction allows for such a coexistence in a large region of the zero temperature phase diagram. In this region, the repulsive part of the interaction drives the stripe formation and the attractive part induces the pairing, resulting in a supersolid with -wave Cooper pairs aligned along the stripes. From a momentum space perspective, the stability of the supersolid phase is due to the fact that the stripe order renders the Fermi surface only partially gapped, leaving gapless regions that are most important for -wave pairing. We finally demonstrate how this supersolid phase can be detected in time-of-flight experiments.

###### pacs:

03.75.Ss, 67.85.Lm, 67.10.Db, 67.80.kbSince the prediction of superfluidity in solid helium several decades ago Andreev (); Leggett (); Chester (); Boninsegni (), the intriguing possibility of coexisting diagonal (density) and off-diagonal (superfluid) order forming a supersolid has been subject to intense investigations. However, the supersolid phase has not been observed unequivocally as interpretations of state-of-the-art experiments in helium are still debated Kim (); Balibar (). The recent experiments on cold dipolar gases Ospelkaus (); Gunton (); Heo (); Wu (); Lev (); Schulze (); Tung (); Repp (), may finally allow an observation of this conceptually important phase. Supersolidity has been predicted to exist for dipolar bosons in an optical lattice Capogrosso (); Pollet (); Cinti (), dipolar bosons with three-body forces Lu (), and spinor Bose condensates with spin-orbit coupling Li1 (); Li2 (). For Fermions however, relevant studies are fewer and limited to the case of an optical lattice Gadsbolle (); Zeng (). In this paper, we expand the scope of study concerning supersolidity of dipolar Fermi gases and show that a supersolid phase is in fact the ground state in a large region of the phase diagram of a two-dimensional (2D) Fermi gas of dipoles aligned by an external field.

We consider spinless Fermions with a dipole moment at zero temperature, confined in the plane by a harmonic trapping potential along the -direction. We take (), where is the Fermi energy of a 2D non-interacting gas with areal density and . In this limit, the Fermions are “frozen” in the harmonic oscillator ground state in the direction and the system is effectively 2D. The dipole moments are aligned by an external field into a direction which is perpendicular to the -axis and forms an angle with respect to the -axis as illustrated in Fig. 1. The interaction between the dipoles is where for electric dipoles, and is the angle between the relative displacement vector of the two dipoles and the dipole moment . For , the Fourier transform of the effective 2D interaction is given by (up to an irrelevant constant term) Fischer ()

(1) |

Here and , where is the trapping length in -direction and is the polar angle of . We note that saturates in the limit of large . This of course is not physical since any true molecular potential has a strong repulsive core which effectively provides a high momentum cut-off for the potential. Such a cut-off can also be introduced by considering the two-body scattering problem as we will discuss below. With a cut-off in mind, we can use the fact that the 2D limit is equivalent to , and make the approximation in (1).

The strength of this interaction is measured by the ratio of the typical interaction and kinetic energy . In addition to , the system is characterised by the dipole tilting angle , which controls the degree of anisotropy of the interaction in the -plane. For weak to moderate interaction strengths and small tilting angles, the system is well described by Landau Fermi liquid theory Baranov2 (); Lu (). For , the repulsion between two dipoles in the plane is strongest when their relative displacement vector is along the direction and weakest along the direction. The anisotropy is predicted to give rise to stripe formation for interaction strengths beyond a critical value Block1 (); Baranov2 (); Babadi (); Parish (); Block2 (). In this phase the dipoles form stripes parallel to the -axis to minimise the repulsion, corresponding to a density modulation with wave vector as illustrated in Fig. 1. As increases the dipolar interaction eventually becomes partially attractive. A Fermi liquid to -wave superfluid phase transition is predicted to occur for tilting angles greater than a critical angle , due to the attractive part of the dipolar interaction Bruun ().

The zero temperature mean-field phase diagram shown in Fig. 2 summarises the discussion given above. In this phase diagram the critical coupling strength for stripe formation is obtained from a Hartree-Fock (HF) calculation Block1 (); Block2 () and the normal Fermi liquid to superfluid transition critical angle is determined by BCS theory (see below). Although each of these three phases have been studied extensively, an interesting question remains unaddressed: what is the nature of the ground state in the region of the phase diagram where the superfluid and stripe phases overlap? In this paper, we provide an answer to this question. We demonstrate that the system in the density wave phase eventually becomes unstable towards pairing as the tilting angle increases and the interaction becomes more attractive. Importantly, the resulting superfluid order does not exclude the stripe order, thus making the system a supersolid as understood in the sense mentioned above.

We use self-consistent HF theory to analyse the stripe phase, as it is the only theory so far that allows us to determine quantitatively properties of this phase. In the stripe phase, the translational symmetry is broken in the direction, and , where annihilates a dipole with momentum . This yields a modulated density where is the stripe order parameter. As described in detail in Ref. Block2 (), the resulting mean-field Hamiltonian can be diagonalised by an unitary transformation , where is restricted to the first Brillouin zone (BZ) , are the reciprocal lattice vectors, and is a band index. The diagonalised Hamiltonian takes the form , where annihilates a quasiparticle with energy . The most salient outcome of the HF analysis is that the Fermi surface contains gapped regions around , as well as gapless regions around and , see Fig. 4. The gapped regions are a manifestation of the stripe order, and the gap magnitude is roughly proportional to .

The key fact for the present purpose is that the stripe order still leaves gapless regions on the Fermi surface, which opens up the intriguing possibility of superfluid pairing. To explore this, we use BCS theory with the Hamiltonian . Here,

describes pairing between the time-reversed states, where is the interaction between the quasiparticles SM (). To derive a gap equation that is amenable to a partial wave expansion, we switch to the “extended zone scheme”, whereby a single particle state in the ’th band in the first BZ is mapped onto a state in the ’th BZ in the standard way SM (), where the vector is now unrestricted. The effective pairing interaction shall be denoted by . Pairing between time-reversed quasiparticles gives rise to the gap parameter , which satisfies the gap equation

(2) |

Here and , where the chemical potential is approximated by the value in the stripe phase. The Cauchy principal value term in (2) renders the gap equation well defined with no need for a high momentum cut-off. Such a term can be introduced by renormalizing the gap equation in terms of scattering amplitude of two dipoles in a vacuum Levinsen (); Cooper (); Baranov2 (). In the absence of experimental data for dipole-dipole scattering in 2D, one can simply regard it as a specific procedure to provide a cut-off.

To solve the gap equation, we expand the gap parameter as where restricts the summation to odd indices, since for spinless Fermions. A more general expansion contains both and terms. However, the terms are favoured by the attractive part of the potential and the gap parameter given by the previous expression maximises the pairing SM (). Using the expansion in (2) we obtain a system of equations

(3) |

where

(4) |

and

(5) |

Equations (3)-(4) with (18) are the fundamental equations to be solved numerically, and they form the basis of the results presented in the rest of this paper.

Before we turn to fully numerical solutions of (3), it is important to understand under what conditions the gap equation admits a solution. To do so, we shall first examine the Fourier components . As shown in the Supplementary Material SM (), these Fourier components calculated numerically from (18) differ very little from those between the bare particles, which are obtained by replacing in (18) by the bare interaction . The latter components, denoted by , can be determined analytically and obey the selection rule only if . In addition, the lowest component is in general dominant over the higher components. We find that the Fourier components possess all the above properties to a very good approximation. The agreement between and holds even deep into the stripe phase, which seems initially surprising since a large stripe amplitude gives rise to extended gapped regions around the Fermi surface. The reason is that the quasiparticle interaction is altered from only in the gapped regions centred at , which are precisely the regions of integration in (18) suppressed by the factor. In light of earlier work on the superfluid transition Bruun (), the fact that the quasiparticle pairing interaction is approximately the same as that between the bare particles strongly suggests that pairing in the stripe phase is possible, provided that the Fermi surface is not fully gapped.

From these results it can be shown that the dominant component of the gap parameter is in fact the first harmonic. Thus the simplest approximation is the momentum-independent ansatz . Since the integrand in (3) is peaked around the partially gapped Fermi surface, a good estimate for when the gap equation admits a finite solution simply follows from the requirement that the effective -wave interaction in the vicinity of the Fermi surface is attractive, i.e.,

(6) |

where is the Fermi wave vector at . The critical angle is therefore , which is the same as that obtained for a normal Fermi liquid to superfluid transition at the same level of approximation Bruun ().

With an estimate of the critical angle, we now solve the gap equation self-consistently including higher harmonics and retaining the full momentum dependence of . The quasiparticle energies and the effective interactions are calculated from the HF theory for the stripe phase and are then used as input to the gap equation (3)-(4). This approach assumes that the pairing has a negligible effect on the stripes, which we will demonstrate is correct. As an example of the calculations, we show in Fig. 3 (left) the amplitudes of the first three partial wave components of the gap parameter for and . For these parameters, the system is deep in the stripe phase with in the absence of pairing. We see that the pairing occurs dominantly in the -wave channel , but also has a noticeable -wave () component; all the higher partial wave components are completely negligible. This feature is in fact typical of solutions to the gap equation. With the solutions of the gap equation we can further determine the redistribution of the quasiparticles and hence, the change of the stripe amplitude as a result of pairing. We find very small relative changes in the stripe amplitude in all our calculations. This demonstrates that our approach is consistent and that the stripe and superfluid orders can indeed coexist forming a type of supersolid. Figure 3 (right) depicts the gap parameter at the tip of the Fermi surface, , as a function of for various values of . For negative but small effective -wave interaction, the behaviour of is well described by the weak pairing approximation , which follows from the ansatz mentioned earlier.

To understand the coexistence of stripe and superfluid orders, we examine the bare particle pair correlation function which is plotted in Fig. 4 (left) for and . It clearly shows that pairing is concentrated in the gapless regions of the underlying Fermi surface for the stripe phase. Consequently, it does not affect the particle distribution in the gapped region which is responsible for the stripe formation. We analyse this further by determining the pair wave function in real space , where is the field operator of the dipoles. In Fig. 4 (right) we show for a Cooper pair with the centre of mass at the origin of the coordinates. The -wave nature of the pairing is clearly visible with strongly peaked along the -axis, where the dipole-dipole interaction is most attractive. In addition, we plot in Fig. 5 the relative probability density of finding the centre-of-mass of a Cooper pair at a specific location. This is given by , which depends only on the -coordinate of the Cooper pair centre-of-mass due to the translational symmetry in the direction. We see that the probability density varies in phase with that of the density of the dipoles, such that the -wave pairing has a maximum on the stripes and a minimum in between. Thus, from the real space perspective, the stripe and superfluid orders coexist due to the anisotropy of the dipolar interaction. Namely, the repulsive part induces the stripe formation while the attractive part induces pairing, resulting in Cooper pairs with -wave symmetry along the stripes as illustrated in Fig. 1. In the limit of strong interaction, the density between the stripes presumably vanishes and the stripes well separate. This raises the interesting possibility of realising an array of 1D -wave superconductors which have topological properties Alicea (). The study of this strong coupling limit is beyond the scope of the present paper.

We argued earlier that is a good approximation for the boundary separating the stripe and the supersolid phases for . We now obtain a more accurate result by varying the tilting angle and determining the critical angle below which the gap equation ceases to admit a finite solution. Such calculations can also be carried out for the normal Fermi liquid to superfluid transition. The overall phase boundary obtained this way is shown in Fig. 2. We see that our initial estimate is in fact remarkably accurate and the phase boundary has a rather weak dependence on the interaction strength . This suggests that the onset of pairing is primarily determined by the degree of anisotropy of the dipolar potential. The identification of the supersolid region in the phase diagram bounded by this boundary, and the collapse line, and our elucidation of the nature of this phase, are the main results of this paper.

Finally, we discuss how the supersolid phase can be detected in time-of-flight (TOF) experiments, which have been used to probe a variety of phases and correlations Greiner (); Folling (); Rom (); Spielman (). As shown in Ref. Block2 (), the density wave order can be detected by measuring the correlation function in TOF experiments. Similarly the superfluid order can be detected by a measurement of the pair correlation function , which can then be compared to theoretical results such as that shown in Fig. 4 (left). However, we need to bear in mind that in a standard experiment, the imaging system introduces a smoothening of the absorption images in the -plane, which can be modelled by convolution of the absorption density with a Gaussian Folling (); BruunAnti (). This reduces the magnitude of the correlation peak considerably from the theoretical maximum value of approximately shown in Fig. 4. Nevertheless, we expect that the TOF experiments can be used to detect the supersolid phase.

In conclusion, we demonstrate that a 2D gas of fermionic dipoles aligned by an external field allow for a coexistence of stripe and superfluid order in a large region of the zero temperature phase diagram. This occurs as a result of the anisotropic nature of the dipolar interaction, where the repulsive part drives the stripe formation, and the attractive part induces the formation of -wave Cooper-pairs along the stripes. In momentum space, the existence of the supersolid phase can be understood from the fact that the stripe order renders the Fermi surface partially gapped, leaving gapless the regions most important for -wave pairing. We finally discuss how the supersolid phase can be detected in TOF experiments. Our results point to several interesting future research directions. This includes realising an array of 1D topological superconductors in the limit of strong interaction, and investigating parallels to the high cuprates, where the co-existence of charge-density-wave order and superconductivity was recently observed Hayward ().

###### Acknowledgements.

GMB would like to acknowledge the support of the Hartmann Foundation via grant A21352 and the Villum Foundation via grant VKR023163.## References

- (1) A. F. Andreev and I. M. Lifshitz, Sov. Phys. JETP 29, 1107 (1969).
- (2) A. J. Leggett, Phys. Rev. Lett. 25, 1543 (1970).
- (3) G. V. Chester, Phys. Rev. A 2, 256 (1970).
- (4) M. Boninsegni and N. V. Prokof’ev, Rev. Mod. Phys. 84, 759 (2012).
- (5) E. Kim and M. H. W. Chan, Nature (London) 427, 225 (2004).
- (6) S. Balibar, Nature (London) 464, 176 (2010).
- (7) K.-K. Ni, S. Ospelkaus, D. Wang, G. Quéméner, B. Neyenhuis, M. H. G. de Miranda, J. L. Bohn, J. Ye, and D. S. Jin, Nature (London) 464, 1324 (2010).
- (8) B. Deh, W. Gunton, B. G. Klappauf, Z. Li, M. Semczuk, J. Van Dongen, and K. W. Madison, Phys. Rev. A 82, 020701 (2010).
- (9) M.-S. Heo, T. T. Wang, C. A. Christensen, T. M. Rvachov, D. A. Cotta, J.-H. Choi, Y.-R. Lee, and W. Ketterle, Phys. Rev. A 86, 021602 (2012).
- (10) Mingwu Lu, Nathaniel Q. Burdick, and Benjamin L. Lev, Phys. Rev. Lett. 108, 215301 (2012).
- (11) C.-H. Wu, J. W. Park, P. Ahmadi, S. Will, and M. W. Zwierlein, Phys. Rev. Lett. 109, 085301 (2012).
- (12) T. A. Schulze, I. I. Temelkov, M. W. Gempel, T. Hartmann, H. Knöckel, S. Ospelkaus, and E. Tiemann, Phys. Rev. A 88, 023401 (2013).
- (13) S.-K. Tung, C. Parker, J. Johansen, C. Chin, Y. Wang, and P. S. Julienne, Phys. Rev. A 87, 010702 (2013).
- (14) M. Repp, R. Pires, J. Ulmanis, R. Heck, E. D. Kuhnle, M. Weidemüller, and E. Tiemann, Phys. Rev. A 87, 010701 (2013).
- (15) B. Capogrosso-Sansone, C. Trefzger, M. Lewenstein, P. Zoller, and G. Pupillo, Phys. Rev. Lett. 104, 125301 (2010).
- (16) L. Pollet, J. D. Picon, H. P. Büchler, and M. Troyer, Phys. Rev. Lett. 104, 125302 (2010).
- (17) F. Cinti, P. Jain, M. Boninsegni, A. Micheli, P. Zoller, and G. Pupillo, Phys. Rev. Lett. 105, 135301 (2010).
- (18) Zhen-Kai Lu, D. S. Petrov, and G. V. Shlyapnikov, arXiv:1409.7737.
- (19) Y. Li, L. P. Pitaevskii, and S. Stringari, Phys. Rev. Lett. 108, 225301 (2012).
- (20) Y. Li, G. I. Martone, L. P. Pitaevskii, and S. Stringari, Phys. Rev. Lett. 110 235302 (2013).
- (21) Anne-Louise Gadsbolle and G. M. Bruun, Phys. Rev. A 85, 021604(R) (2012).
- (22) T. S. Zeng and L. Yin, Phys. Rev. B 89, 174511 (2014).
- (23) U. R. Fischer, Phys. Rev. A 73, 031602 (2006).
- (24) Y. Yamaguchi, T. Sogo, T. Ito and T. Miyakawa, Phys. Rev. A 82, 013643 (2010).
- (25) L. M. Sieberer and M. A. Baranov, Phys. Rev. A 84 063633 (2011).
- (26) J. K. Block, N. T. Zinner, and G. M. Bruun, New Journal of Physics 14, 105006 (2012).
- (27) M. Babadi and E. Demler, Phys. Rev. B 84, 235124 (2011).
- (28) M. M. Parish and F. M. Marchetti, Phys. Rev. Lett. 108, 145304 (2012).
- (29) J. K. Block and G. M. Bruun, Phys. Rev. B90, 155102 (2014)
- (30) G. M. Bruun and E. Taylor, Phys. Rev. Lett. 101, 245301 (2008).
- (31) See Supplemental Material for a description of the “extended zone scheme”, pairing symmetry and a comparison of the bare and quasiparticle pairing interactions.
- (32) J. Levinsen, N. R. Cooper, and G. V. Shlyapnikov, Phys. Rev. A 84, 013603 (2011).
- (33) N. R. Cooper and G. V. Shlyapnikov, Phys. Rev. Lett. 103, 155302 (2009).
- (34) J. Alicea, Rep. Prog. Phys. 75, 076501 (2012).
- (35) M. Greiner, C. A. Regal, J. T. Stewart, and D. S. Jin, Phys. Rev. Lett. 94, 110401 (2005).
- (36) S. Fölling, F. Gerbier, A. Widera, O. Mandel, T. Gericke, and I. Bloch, Nature (London) 434, 481 (2005).
- (37) T. Rom, T. Best, D. van Oosten, U. Schneider, S. Fölling, B. Paredes, and I. Bloch, Nature (London) 444, 733 (2006).
- (38) I. B. Spielman, W. D. Phillips, and J. V. Porto, Phys. Rev. Lett. 98, 080404 (2007).
- (39) T. Sogo, L. He, T. Miyakawa, S. Yi, H. Lu, and H. Pu, New J. Phys. 11, 055017 (2009).
- (40) E. Altman, E. Demler, and M. D. Lukin, Phys. Rev. A 70, 013603 (2004).
- (41) G. M. Bruun, O. F. Syljuåsen, K. G. L. Pedersen, B. M. Andersen, E. Demler, and A. S. Sørensen, Phys. Rev. A 80, 033622 (2009).
- (42) L. E. Hayward, D. G. Hawthorn, R. E. Melko, and S. Sachdev, Science 343, 1336 (2014).

## Appendix A Supplemental Material

### a.1 1: Hartree-Fock theory on the stripe phase and the extended zone scheme

The mean-field Hamiltonian used to describe this phase is given by Block2 ()

(7) |

where , is the single particle Hartree-Fock energy

(8) |

and is a real off-diagonal element defined by

(9) |

The inclusion of the second term in Eq. (7) accounts for the possibility of formation of the density wave along the direction. The Hamiltonian in Eq. (7) resembles (although is not identical to) that of non-interacting particles in a potential periodic in the direction with periodicity . Consequently the quasiparticle eigenlevels of Eq. (7) exhibit a band-like structure along the direction of the wave vector, where is the band index and is restricted to the first Brillouin zone. The corresponding quasiparticle wave function can be expressed as , where is the reciprocal lattice vector and the expansion coefficients are determined the Schrödinger equation

(10) |

Equation (10) is analogous to the Schrödinger equation of a particle in a periodic lattice, where plays the role of the Fourier components of a “periodic potential”. Unlike a true periodic potential, however, depends explicitly on due to the inclusion of the exchange interaction. In terms of the quasiparticle wave function basis, the Hamiltonian in (7) can now be brought into a diagonalised form , where is the annihilation operator of the quasiparticle. The quasiparticle occupation number in the ground state is specified by the chemical potential of the density wave phase, which in turn is determined by the density of the gas as

(11) |

The quasiparticle energy and the expansion coefficients are implicit functions of the Hartree-Fock elements and . Therefore these quantities as well as the chemical potential are determined self-consistently through Eqs. (8)-(11). The reader is referred to Ref. Block2 () for a detailed account of their numerical calculation.

It turns out that the effects of the off-diagonal terms in Eq. (7) to the quasiparticle dispersion are only perturbative, due to the fact the magnitudes of are generally small compared to the Fermi energy Block2 (). Consequently the quasiparticle dispersion does not in fact deviate significantly from the usual parabolic form except in regions close to the Brillouin zone boundaries where band gaps open up. It is thus meaningful to use the “extended zone scheme” Ashcroft () instead of the “reduced zone scheme” in labelling the quasiparticle energy levels. More specifically, each of the physical quantities associated with the single particle state can be labelled by a single wave vector in the -th Brillouin zone , which is defined as

(12) |

for and

(13) |

for . Likewise, the physical quantities associated with the time-reversal state can be labelled by the vector . The effective pairing interaction, given by

(14) |

shall be denoted by . These correspondences are made clear if one considers the limit of vanishing off-diagonal Hartree-Fock elements . In this limit the Bloch state simply approaches the plane wave state and the effective pairing interaction approaches the bare interaction .

### a.2 2: The pairing symmetry

Here we provide a rationale for the choice of a gap parameter with even parity with respect to as expressed in the main Letter. The gap equation (4) in the main Letter can be written as

(15) |

where is the anti-symmetrized interaction matrix

(16) |

It can be shown that the interaction matrix has the following expansion

(17) |

where restricts the summation to odd indices,

(18) |

and

(19) |

The sine and cosine terms in the expansion (17) are not coupled and, as a consequence, the gap equation (15) admits solutions with either even or odd parity with respect to the variable. For even solutions only the first part of the potential in (17) contributes to the integral in Eq. (15) and for odd solutions only the second part does. In order to see which part of the potential favours Cooper pairing, we express and in terms of the interaction potential in real space . This can be done by approximating by the bare interaction matrix

(20) |

in Eqs. (18) and (19). Using the expansion , where is the Bessel function of the first kind, and performing the integrals with respect to and , we find

(21) |

and

(22) |

From these expressions we see that mostly samples the attractive part of the potential in the real space (the sliver around the axis) while mostly samples the repulsive part. Let us take the the first diagonal elements and for example. These are in fact the most dominant matrix elements for and respectively. As the attractive sliver of the potential expands from the axis with an increasing tilting angle , can potentially become negative due to the fact the attractive part of potential is more significantly weighted in the integral. The matrix element , on the other hand, remains positive for all tilting angles. This analysis motivates us to look for solutions to the gap equation with even parity with respect to .

### a.3 3: Comparisons between and

The Fourier components (which are defined by Eq. (18) with replaced by ) for the bare interaction can be evaluated analytically. Using Eq. (2) of the main Letter in Eq. (18) we find that the only non-vanishing matrix elements are those whose indices differ by 0 or . That is, the matrix has the following tridiagonal structure

(23) |

For odd we find

(24) |

(25) |

and

(26) |

where

(27) |

with . The integral can generally be expressed in terms of the complete elliptic integrals. As a few examples, for are shown below as

(28) |

(29) |

and

(30) |

where and are the complete elliptic integrals of the first and second kind respectively. As a useful gauge of the relative importance of the matrix elements for various , we consider the special case of . In this case the expressions in (24)-(26) simplify and we find

(31) |

and

(32) |

We see that the magnitudes of these matrix elements decreases rapidly as . The Fourier components for the quasiparticles are calculated numerically using Eq. (18). We find that the formation of density wave has minimal effects on the effective paring interaction, namely agrees very well with . (see Figs. 6-7 as an example) for a wide range of and .

## References

- (43) J. K. Block and G. M. Bruun, Phys. Rev. B 90, 155102 (2014)
- (44) N. W. Ashcroft and N. D. Mermin, Solid State Physics, (Holt, Rinehart, and Winston, New York, 1976).