# Stabilising Fulde-Ferrel-Larkin-Ovchinnikov Superfluidity with long-range Interactions in a mixed dimensional Bose-Fermi System

###### Abstract

We analyse the stability of inhomogeneous superfluid phases in a system consisting of identical fermions confined in two layers that are immersed in a Bose-Einstein condensate (BEC). The fermions in the two layers interact via an induced interaction mediated by the BEC, which gives rise to pairing. We present zero temperature phase diagrams varying the chemical potential difference between the two layers and the range of the induced interaction, and show that there is a large region where an inhomogeneous superfluid phase is the ground state. This region grows with increasing range of the induced interaction and it can be much larger than for a corresponding system with a short range interaction. The range of the interaction is controlled by the healing length of the BEC, which makes the present system a powerful tunable platform to stabilise inhomogeneous superfluid phases. We furthermore analyse the melting of the superfluid phases in the layers via phase fluctuations as described by the Berezinskii-Kosterlitz-Thouless mechanism and show that the normal, homogeneous and inhomogeneous superfluid phases meet in a tricritical point. The superfluid density of the inhomogeneous superfluid phase is reduced by inherent gapless excitations, and we demonstrate that this leads to a significant suppression of the critical temperature as compared to the homogeneous superfluid phase.

## I Introduction

The interplay between population imbalance and superfluid pairing has been subject to intense study ever since Fulde and Ferrell (FF) as well as Larkin and Ovchinnikov (LO) predicted that they can co-exist Fulde and Ferrell (1964); Larkin and Ovchinnikov (1965). In condensed matter systems, an external magnetic field leads to a population imbalance between the two electron spin projections, which in general is at odds with superconductivity. FFLO however realised that the superconductor can accommodate some population imbalance at the price of giving the Cooper pairs a non-zero center-of-mass (COM) momentum, thereby forming a spatially inhomogeneous but periodic order parameter with no vortices. The fate of superfluid pairing in the presence of population imbalance is a fundamental question relevant for many systems in nature including cold atoms Chevy and Mora (2010); Radzihovsky and Sheehy (2010); Kinnunen et al. (2018), superconductors Matsuda and Shimahara (2007), and quark matter Casalbuoni and Nardulli (2004); Alford et al. (2008); Anglani et al. (2014). Nevertheless, an unambiguous observation of a FFLO phase is still lacking. A major problem for electronic superconductors is that orbital effects due to the magnetic field lead to the formation of vortices and eventually destroy pairing before any FFLO physics can be observed. One strategy to avoid this problem is to explore low dimensional systems, where orbital effects are suppressed due to the confinement. Indeed, results consistent with a FFLO phase have been reported for quasi-2D organic and heavy fermion superconductors Beyer and Wosnitza (2013); Mayaffre et al. (2014); Bianchi et al. (2003). Theoretical studies have furthermore concluded that the FFLO phase is favored in 2D as compared to 3D Shimahara (1994); Burkhardt and Rainer (1994); Pieri et al. (2007).

Quantum degenerate atomic gases are well suited to investigate FFLO physics, because they do not suffer from orbital effects as the atoms are neutral. In addition, it is relatively straightforward to make low dimensional systems, and signatures of FFLO physics has indeed been observed in a one-dimensional (1D) atomic Fermi gas Liao et al. (2010). There has been a number of investigations of the FFLO phase for 2D atomic gases with a short range interaction Gukelberger et al. (2016); Baarsma and Törmä (2016); Conduit et al. (2008); Radzihovsky and Vishwanath (2009); Wolak et al. (2012); Yin et al. (2014); Sheehy (2015). Recently, it was argued that long range interactions further increases the region of stability of the FFLO phase for a 2D gas of dipolar atoms as compared to a short range interaction Lee et al. (2017).

Here, we investigate how to stabilise FFLO superfluidity using a mixed dimensional system with a tuneable range interaction. The system consists of fermions confined in two layers immersed in a BEC. The induced interaction between the layers mediated by the BEC gives rise to pairing, and we analyse the stability of the corresponding superfluid phases. We show that the zero temperature phase diagram as a function of the chemical potential difference in the two layers and the range of the interaction has a large region, where an inhomogeneous superfluid phase is the ground state. This region grows with increasing range of the interaction and becomes much wider than for a zero range interaction. The interaction range can be tuned by varying the healing length of the BEC. We furthermore investigate the melting of the 2D superfluid phases via the Berezinskii-Kosterlitz-Thouless mechanism. The normal phase, the homogeneous and inhomogeneous superfluid phases are shown to meet in a tricritical point in the phase diagram, which determines the maximum critical temperature of the inhomogeneous superfluid. This maximum temperature is however significantly suppressed compared to the homogeneous superfluid phase, due to inherent gapless excitations, which decrease the superfluid density.

## Ii The Bilayer System

We consider the system illustrated in Fig. 1. Two layers contain fermions of a single species with mass , and they are immersed in a 3D weakly interacting BEC consisting of bosons with mass and density . The distance between the layers is and the surface density of fermions in each layer is with denoting the two layers. When used in equations, they are taken to mean , respectively. Occasionally we will also use the notation to mean the opposite layer from . The boson-fermion interaction is short range and characterised by the strength with the reduced mass, and the effective 2D-3D scattering length Nishida and Tan (2008). Throughout, this interaction strength is taken to be weak in the sense , where is the Fermi momentum of layer .

We treat the BEC using zero temperature Bogoliubov theory. This is a good approximation since the critical temperature of the superfluid phases of the fermions is much smaller than the critical temperature of the BEC. The bosonic degrees of freedom can then be integrated out, which yields an effective interaction between the fermions. In the static limit, this interaction is on the Yukawa form, and we end up with an effective Hamiltonian for the fermions in the two layers on the form Heiselberg et al. (2000); Wu and Bruun (2016); Midtgaard et al. (2017)

(1) |

Here, creates a fermion in layer with 2D momentum , the dispersion in each layer is with the chemical potentials, and is the system volume. We define the average chemical potential and the “magnetic field” , so that we can write . The Yukawa interaction is

(2) |

where is the healing length of the BEC with the boson-boson scattering length. The healing length determines the range of the induced interaction as can be seen by Fourier transforming Eq. (2) by to real space giving with the 3D distance between the fermions. It follows that the range of the interaction can be tuned by varying the density or the scattering length of the surrounding BEC, which turns out to be a key property for the following. We note that retardation effects can be neglected, when the speed of sound in the BEC is much larger than the Fermi velocity in the two planes Wu and Bruun (2016).

## Iii Pairing and Green’s functions

The attractive interaction given by Eq. (2) can lead to Cooper pairing within each layer (intralayer pairing), and between the two layers (interlayer pairing). We recently analysed the competition between intra- and interlayer pairing for the case of equal density in each layer, i.e. for Midtgaard et al. (2017). For a layer distance large compared to the range of the interaction, we found that the ground state is characterised by intralayer -wave pairing, whereas -wave interlayer pairing is stable for smaller . In addition, we identified a crossover phase for intermediate where both types of pairing co-exist.

In this paper, we focus on interlayer -wave pairing corresponding to . We shall investigate the case of a non-zero field giving rise to a population difference between the two layers, and the possibility of FFLO interlayer pairing with non-zero COM momentum. Such interlayer pairing with COM momentum is characterised by the anomalous averages , which leads us to define the corresponding pairing field

(3) |

We have dropped the subscripts on the induced interaction , as it here and in the following refers to the interlayer interaction only, i.e. . The pairing field obeys the Fermi anti-symmetry . In real space, it is on the form

(4) |

where and are the COM and relative coordinates respectively, and is the field operator for particles in layer . We note that the pairing field is not translationally symmetric in the FFLO phase, i.e. . The interaction becomes in the mean-field BCS approximation

(5) |

where we include only the pairing channel as we focus on the superfluid instability. The Hartree-Fock terms will in general lead to small effects in the weak coupling regime, which mostly can be accounted for by a renormalisation of the chemical potentials .

All results presented here can be obtained by a direct diagonalisation of the mean-field BCS Hamiltonian using a standard Bogoliubov transformation. We however use Green’s functions to analyse the FFLO states, since this formalism naturally allows us to go beyond mean-field theory to include effects such as retardation, if needed in the future. The normal and anomalous Green’s functions for the superfluid phases are defined in the standard way as

(6) |

where denotes imaginary-time ordering. Using , the Gor’kov equations for these Green’s functions are straightforwardly derived. They read

(7) | ||||

(8) | ||||

(9) |

where we have Fourier transformed to Matsubara frequency space with . The self-consistent gap equation is then from (3) and (6)

(10) |

## Iv Fulde-Ferrel Superconductivity

The values for which is non-zero determines the structure of the order parameter in the superfluid phase. It has been shown that can form very complicated 2D structures corresponding to for many s in Eq. (4) Shimahara (1998); Mora and Combescot (2004). However, the FFLO phase exhibits a second order transition to the normal phase in 2D at an upper critical field Shimahara (1994); Burkhardt and Rainer (1994), and at this transition it is sufficient to consider the case of only for a single vector. The reason is that any linear combination of the ’s, which are unstable towards pairing, is degenerate at , since it is only the non-linear part of the gap equation that determines the optimal combination that minimises the energy.

In the following, we therefore consider the case only for a single . This will give the correct upper critical field for the second order transition between the FFLO and the normal phase. We also expect that it will give a fairly precise value for the lower critical field determining the first order transition between the FFLO and the superfluid phase, since the energy difference between the FFLO phases with various spatial structures is small Shimahara (1998); Mora and Combescot (2004). Our scheme recovers the usual homogenous BCS pairing for , and it corresponds to a plane wave Fulde-Ferrel (FF) type of pairing when , as can be seen from Eq. (4). Since we only have one vector, we can simplify the notation for the gap as . The Gor’kov equations (7)-(9) are then easily solved giving

(11) |

Using Eq. (11) in Eq. (10), we obtain the gap equation

(12) |

where , with , and the Fermi distribution function is . In Eq. (12) and the rest of this paper, we have further simplified the notation by defining and . We solve Eq. (12) along with the number equation

(13) |

Note that in order to compare with previous results in the literature, we keep the total number of particles, , fixed, and not the number of particles in each plane. To solve the gap equation (12), we perform a partial wave expansion of the induced interaction , where is the azimuthal angle of . In the numerics, we keep the two leading terms, corresponding to -wave singlet and -wave triplet pairing respectively. When solving the above equations, we find three different phases: The homogeneous BCS phase with , the FF state with , and the normal phase with no superfluid pairing. To determine which of these phases is the ground state, we compare their energy . We have from mean-field theory

(14) |

## V Results for zero temperature

We now solve the Gor’kov equations for , varying the interaction range and density imbalance through the healing length and the field , respectively. The resulting phase diagrams are shown in Figs. 2 and 3 for two different values of the layer distance, and . The Fermi momentum is defined as . We set the effective scattering length to . This small value ensures that we stay in the valid range of mean-field theory for all parameter sets shown in the figures. As is standard in the literature, we measure the field in units of , which is the pairing field at the Fermi surface for , that is, with no density imbalance.

Consider first the case of the layer distance . Fig. 2 clearly shows that there is a large region in the phase diagram, where the FF phase is stable. We find that the phase transition between the BCS and the FF phase is first order at the lower critical field , whereas it is second order for the transition between the FF phase and the normal phase at the upper critical field . This is in agreement with the results for a short range interaction Shimahara (1994); Burkhardt and Rainer (1994). Moreover, the range of values of for which the FF phase is the ground state increases with the interaction range . This shows that a long range interaction stabilises the FF phase. The reason is that the relative strength of the -wave component compared to the -wave component of the interaction increases with increasing range, which favors FF pairing. To illustrate this important point further, we plot as horizontal lines in Fig. 2 the critical fields for a short range interaction Fulde and Ferrell (1964); Lee et al. (2017), and . While the FF region approaches that of a short range interaction for decreasing , it becomes much larger with increasing range .

Consider next the case of zero layer distance shown in Fig. 3. While a large still stabilises the FF phase, the effect here is much less pronounced as compared to the case . Increasing leads to a smaller increase in the range of values of for which the FF phase is the ground state than for . The reason is that the short range divergence of the Yukawa interaction between the fermions in the two layers given by Eq. (2) is cut off at for a finite layer distance . This makes the -wave part of the interaction stronger compared to the -wave part. As a result, the FF phase where there is pairing in both the - and -wave channels is favored compared to the pure -wave BCS state for a non-zero layer distance. Note that the reason the superfluid region (BCS and FF) seems larger for compared to even though the strength of the interaction obviously is smaller for a non-zero layer distance, is that we measure in units of , which is also smaller. Had we used the unit instead for instance, the superfluid region of the phase would be larger.

## Vi Berezinskii-Kosterlitz-Thouless Melting

Since the fermions are confined in 2D layers, phase fluctuations of the order parameter are significant and will eventually melt the superfluid through the Berezinskii-Kosterlitz-Thouless (BKT) mechanism at a critical temperature . To describe this, we first calculate the superfluid stiffness, or equivalently the superfluid density, by imposing a linear phase twist on the order parameter and calculating the corresponding energy cost to second order in the twist. For a given vector , the real space pairing becomes using Eq. (4)

(15) |

where and is the imposed spatially linear phase twist. From Eq. (15), it is clear that the direction of the phase twist relative to the COM of the Cooper pairs is important: When is parallel to , the phase twist corresponds to adding/removing COM momentum to the Cooper pairs which compresses/expands the wavelength of the plane wave pairing field ; when is perpendicular to , the phase twist corresponds to a small rotation of the COM momentum to the Cooper pairs which rotates the plane wave pairing field. These two effects are illustrated in Fig. 4(a).

The phase twist gives a free energy cost of the form

(16) |

where and are spatial derivatives parallel or perpendicular to the COM momentum of the Cooper pairs, corresponding to and respectively. The associated superfluid stiffness constants are and . In the second line of Eq. (16), we have rescaled the spatial coordinate perpendicular to by the factor to obtain an isotropic -model with the effective stiffness constant Wu et al. (2016). Alternatively, defining the superfluid densities parallel and perpendicular to as and , we can write the free energy cost as

(17) |

Here is the superfluid velocity parallel to and likewise for . A long but straightforward calculation of the second order energy shift due to the phase twist gives

(18) |

for the superfluid density along , where . Here, is the total surface density of fermions coming from the two layers. An equivalent formula holds for . Equation (18) is the 2D version of the usual expression for the superfluid density allowing for the effects of the spatial anisotropy of the FF state Leggett (1975); Landau et al. (1980). We can now determine the critical temperature of the superfluid using the BKT condition

(19) |

where we have defined the effective superfluid density as .

In Fig. 4(b), we plot the superfluid densities as a function of for layer distance and boson-fermion interaction strength . We see that , reflecting that the energy cost related to compressing/expanding the COM momentum is higher than that related to rotating it as expected. Note that both superfluid densities are smaller than the total density even for . This is due to the inherent presence of gapless quasiparticle states in the FF phase. These gapless excitations, which are shown in the inset of Fig. 4(b), reduce the superfluid density. We also plot the effective superfluid density as well as the line . It follows from Eq. (19) that the superfluid phase melts when crosses this line.

Having set up the theory for BKT melting, we can now analyse the BCS and FF phases at non-zero temperatures. In Fig. 5, we present an example of a phase diagram for and . We see that the critical temperature of the FF phase increases with decreasing field . The highest critical temperature is obtained just above the lowest critical field , where the FF, BCS, and normal phase meet in a tricritical point at and . This critical temperature is well below the theoretical maximum of obtained by setting in Eq. (19). The reason is that the gapless excitations in the FF phase decrease the superfluid density below , as we discussed above in connection with Fig. 4. Indeed, we note that the value of the two components of the superfluid density in Fig. 4 at are almost unchanged from their value at . This is a general result: While the flexibility of the bilayer system allows us to optimize the induced interaction to favour FF pairing, the gapless nature of this state prevents us from reaching critical temperatures close to .

The critical temperature of the BCS phase is much higher as can be seen from Fig. 5. Due to its fully gapped spectrum, can in fact relatively easily be tuned to be close to maximum by varying the layer distance , the interaction range , even while keeping the boson-fermion interaction strength weak.

## Vii Conclusion

We analysed a mixed-dimensional system consisting of two layers of identical fermions immersed in a BEC. This system was shown to support superfluid pairing due to an attractive induced interaction between the two layers mediated by the BEC. When the densities of the two layers are different, the resulting superfluid phase is inhomogeneous. Using a plane wave FF ansatz to describe this phase, we demonstrated that it is stabilised by the nonzero range of the induced interaction. Importantly, the FF phase can occupy much larger regions of the phase diagram as compared to the case of a short range interaction. The range of the induced interaction can be tuned by varying the BEC healing length, which makes the present system promising for realising FFLO physics experimentally. We furthermore analysed the melting of the superfluid phases due to phase fluctuations using BKT theory, and demonstrated that the normal, homogenous and inhomogeneous superfluid phases meet in a tricritical point in the phase diagram. The superfluid density of the FF phase was shown to be suppressed by intrinsic gapless excitations, and this leads to a significant reduction in the critical temperature compared to the homogeneous superfluid, which can be tuned to have a critical temperature close to the maximum .

###### Acknowledgements.

We wish to acknowledge the support of the Danish Council of Independent Research – Natural Sciences via Grant No. DFF - 4002-00336.## References

- Fulde and Ferrell (1964) P. Fulde and R. A. Ferrell, Phys. Rev. 135, A550 (1964).
- Larkin and Ovchinnikov (1965) A. Larkin and I. Ovchinnikov, Soviet Physics-JETP 20, 762 (1965).
- Chevy and Mora (2010) F. Chevy and C. Mora, Reports on Progress in Physics 73, 112401 (2010).
- Radzihovsky and Sheehy (2010) L. Radzihovsky and D. E. Sheehy, Reports on Progress in Physics 73, 076501 (2010).
- Kinnunen et al. (2018) J. J. Kinnunen, J. E. Baarsma, J.-P. Martikainen, and P. Törmä, Reports on Progress in Physics 81, 046401 (2018).
- Matsuda and Shimahara (2007) Y. Matsuda and H. Shimahara, Journal of the Physical Society of Japan 76, 051005 (2007), https://doi.org/10.1143/JPSJ.76.051005 .
- Casalbuoni and Nardulli (2004) R. Casalbuoni and G. Nardulli, Rev. Mod. Phys. 76, 263 (2004).
- Alford et al. (2008) M. G. Alford, A. Schmitt, K. Rajagopal, and T. Schäfer, Rev. Mod. Phys. 80, 1455 (2008).
- Anglani et al. (2014) R. Anglani, R. Casalbuoni, M. Ciminale, N. Ippolito, R. Gatto, M. Mannarelli, and M. Ruggieri, Rev. Mod. Phys. 86, 509 (2014).
- Beyer and Wosnitza (2013) R. Beyer and J. Wosnitza, Low Temperature Physics 39, 225 (2013), https://doi.org/10.1063/1.4794996 .
- Mayaffre et al. (2014) H. Mayaffre, S. Krämer, M. Horvatić, C. Berthier, K. Miyagawa, K. Kanoda, and V. F. Mitrović, Nature Physics 10, 928 EP (2014).
- Bianchi et al. (2003) A. Bianchi, R. Movshovich, C. Capan, P. G. Pagliuso, and J. L. Sarrao, Phys. Rev. Lett. 91, 187004 (2003).
- Shimahara (1994) H. Shimahara, Phys. Rev. B 50, 12760 (1994).
- Burkhardt and Rainer (1994) H. Burkhardt and D. Rainer, Annalen der Physik 506, 181 (1994), https://onlinelibrary.wiley.com/doi/pdf/10.1002/andp.19945060305 .
- Pieri et al. (2007) P. Pieri, D. Neilson, and G. C. Strinati, Phys. Rev. B 75, 113301 (2007).
- Liao et al. (2010) Y.-a. Liao, A. S. C. Rittner, T. Paprotta, W. Li, G. B. Partridge, R. G. Hulet, S. K. Baur, and E. J. Mueller, Nature 467, 567 EP (2010).
- Gukelberger et al. (2016) J. Gukelberger, S. Lienert, E. Kozik, L. Pollet, and M. Troyer, Phys. Rev. B 94, 075157 (2016).
- Baarsma and Törmä (2016) J. E. Baarsma and P. Törmä, Journal of Modern Optics 63, 1795 (2016), https://doi.org/10.1080/09500340.2015.1128009 .
- Conduit et al. (2008) G. J. Conduit, P. H. Conlon, and B. D. Simons, Phys. Rev. A 77, 053617 (2008).
- Radzihovsky and Vishwanath (2009) L. Radzihovsky and A. Vishwanath, Phys. Rev. Lett. 103, 010404 (2009).
- Wolak et al. (2012) M. J. Wolak, B. Grémaud, R. T. Scalettar, and G. G. Batrouni, Phys. Rev. A 86, 023630 (2012).
- Yin et al. (2014) S. Yin, J.-P. Martikainen, and P. Törmä, Phys. Rev. B 89, 014507 (2014).
- Sheehy (2015) D. E. Sheehy, Phys. Rev. A 92, 053631 (2015).
- Lee et al. (2017) H. Lee, S. I. Matveenko, D.-W. Wang, and G. V. Shlyapnikov, Phys. Rev. A 96, 061602 (2017).
- Nishida and Tan (2008) Y. Nishida and S. Tan, Phys. Rev. Lett. 101, 170401 (2008).
- Heiselberg et al. (2000) H. Heiselberg, C. J. Pethick, H. Smith, and L. Viverit, Phys. Rev. Lett. 85, 2418 (2000).
- Wu and Bruun (2016) Z. Wu and G. M. Bruun, Phys. Rev. Lett. 117, 245302 (2016).
- Midtgaard et al. (2017) J. M. Midtgaard, Z. Wu, and G. M. Bruun, Phys. Rev. A 96, 033605 (2017).
- Shimahara (1998) H. Shimahara, Journal of the Physical Society of Japan 67, 736 (1998), https://doi.org/10.1143/JPSJ.67.736 .
- Mora and Combescot (2004) C. Mora and R. Combescot, EPL (Europhysics Letters) 66, 833 (2004).
- Wu et al. (2016) Z. Wu, J. K. Block, and G. M. Bruun, Scientific Reports 6, 19038 EP (2016).
- Leggett (1975) A. J. Leggett, Rev. Mod. Phys. 47, 331 (1975).
- Landau et al. (1980) L. Landau, E. Lifshitz, and L. Pitaevskij, Statistical Physics: Part 2 : Theory of Condensed State, Landau and Lifshitz Course of theoretical physics (Oxford, 1980).