The Fulde-Ferrell-Larkin-Ovchinnikov state for ultracold fermions in lattice and harmonic potentials: a review

TL;DR

This review discusses the theoretical understanding and experimental prospects of the FFLO state in ultracold fermionic gases, emphasizing systems like optical lattices and spin-orbit coupled setups, and addressing key open questions.

Contribution

It provides a comprehensive overview of the theoretical basis, potential experimental realizations, and open questions regarding the FFLO state in ultracold quantum gases.

Findings

Identifies promising systems for FFLO observation, such as optical lattices and spin-orbit coupled gases.

Analyzes theoretical bounds and experimental limitations affecting FFLO detection.

Summarizes current ideas and open questions about FFLO stability and properties.

Abstract

We review the concepts and the present state of theoretical studies of spin-imbalanced superfluidity, in particular the elusive Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) state, in the context of ultracold quantum gases. The comprehensive presentation of the theoretical basis for the FFLO state that we provide is useful also for research on the interplay between magnetism and superconductivity in other physical systems. We focus on settings that have been predicted to be favourable for the FFLO state, such as optical lattices in various dimensions and spin-orbit coupled systems. These are also the most likely systems for near-future experimental observation of the FFLO state. Theoretical bounds, such as Bloch's and Luttinger's theorems, and experimentally important limitations, such as finite-size effects and trapping potentials, are considered. In addition, we provide a comprehensive review of the various ideas presented for the observation of the FFLO state. We conclude our review with an analysis of the open questions related to the FFLO state, such as its stability, superfluid density, collective modes and extending the FFLO superfluid concept to new types of lattice systems.

Equations

$$\hat{H}=\int d\mathbf{r}\sum _{\sigma =\uparrow ,\downarrow }{\hat{\psi }}_{\sigma }^{†}\left(\mathbf{r}\right)\left({\hat{K}}_{\sigma }(\mathbf{r})-{\mu }_{\sigma }^0\right){\hat{\psi }}_{\sigma }\left(\mathbf{r}\right)$$

$$+\int \int d\mathbf{r}d{\mathbf{r}}^{\prime }{V}_{\uparrow \downarrow }\left(\mathbf{r},{\mathbf{r}}^{\prime }\right){\hat{\psi }}_{\uparrow }^{†}\left(\mathbf{r}\right){\hat{\psi }}_{\downarrow }^{†}\left({\mathbf{r}}^{\prime }\right){\hat{\psi }}_{\downarrow }\left({\mathbf{r}}^{\prime }\right){\hat{\psi }}_{\uparrow }\left(\mathbf{r}\right),$$

$$V_{\uparrow \downarrow }\left(\mathbf{r},{\mathbf{r}}^{\prime }\right)=V_0\delta \left(\mathbf{r}-{\mathbf{r}}^{\prime }\right),$$

$$\hat{H}=$$

$${\hat{\psi }}_{\sigma }^{†}(\mathbf{r})=\sum _{\mathbf{n},\mathbf{i}}{\hat{c}}_{\sigma \mathbf{n},\mathbf{i}}^{†}{w}_{\mathbf{n}}^{\ast }(\mathbf{r}-{\mathbf{r}}_{\mathbf{i}}),$$

$$\hat{H}=-$$

$$t_{\sigma \mathbf{ij}}$$

$${\epsilon }_{\sigma }$$

$$U$$

$$\hat{H}=-\sum _{\sigma }\sum _{⟨\mathbf{r},{\mathbf{r}}^{\prime }⟩}{t}_{\mathbf{ij}}{\hat{c}}_{\sigma \mathbf{r}}^{†}{\hat{c}}_{\sigma {\mathbf{r}}^{\prime }}-\sum _{\sigma \mathbf{r}}{\mu }_{\sigma }^0{\hat{c}}_{\sigma \mathbf{r}}^{†}{\hat{c}}_{\sigma \mathbf{r}}$$

$$+U\sum _{\mathbf{r}}{\hat{c}}_{\uparrow \mathbf{r}}^{†}{\hat{c}}_{\downarrow \mathbf{r}}^{†}{\hat{c}}_{\downarrow \mathbf{r}}{\hat{c}}_{\uparrow \mathbf{r}},$$

$$⟨\hat{\Delta }\left(\mathbf{r}\right)⟩\equiv \Delta (\mathbf{r})=V_0⟨{\hat{\psi }}_{\downarrow }\left(\mathbf{r}\right){\hat{\psi }}_{\uparrow }\left(\mathbf{r}\right)⟩,$$

$$\Delta \left(\mathbf{r}\right)=U⟨{\hat{c}}_{\downarrow \mathbf{r}}{\hat{c}}_{\uparrow \mathbf{r}}⟩.$$

$$(-{\mu }_{\sigma }^0+U{n}_{-\sigma }){\hat{c}}_{\sigma }^{†}{\hat{c}}_{\sigma }\equiv -{\mu }_{\sigma }{\hat{c}}_{\sigma }^{†}{\hat{c}}_{\sigma },$$

$$\hat{H}=-\sum _{\sigma }\sum _{⟨\mathbf{r},{\mathbf{r}}^{\prime }⟩}{t}_{\mathbf{ij}}{\hat{c}}_{\sigma \mathbf{r}}^{†}{\hat{c}}_{\sigma {\mathbf{r}}^{\prime }}-\sum _{\sigma }\sum _{\mathbf{r}}{\mu }_{\sigma }{\hat{c}}_{\sigma \mathbf{r}}^{†}{\hat{c}}_{\sigma \mathbf{r}}$$

$$+\sum _{\mathbf{r}}\left({\hat{c}}_{\uparrow \mathbf{r}}^{†}{\hat{c}}_{\downarrow \mathbf{r}}^{†}\Delta (\mathbf{r})+{\Delta }^{\ast }(\mathbf{r}){\hat{c}}_{\downarrow \mathbf{r}}{\hat{c}}_{\uparrow \mathbf{r}}-\frac{{\Delta }^{\ast }(\mathbf{r})\Delta (\mathbf{r})}{U}\right),$$

$$\Delta (\mathbf{r})=\Delta e^{i\mathbf{q}\cdot \mathbf{r}},$$

$$\Delta (\mathbf{r})=\Delta \cos (\mathbf{q}\cdot \mathbf{r}).$$

$$\hat{H}=$$

$${\hat{c}}_{\sigma \mathbf{r}}$$

$$\sum _{\mathbf{r}}{\hat{c}}_{\sigma \mathbf{r}}^{†}{\hat{c}}_{\sigma \mathbf{r}}$$

$$\sum _{\mathbf{r}}$$

$$\sum _{\sigma }$$

$$\hat{H}=\sum _{\mathbf{k}}$$

$$\xi_{\sigma{\bf k}}=\epsilon_{\bf k}-\mu_{\sigma}=\sum_{\alpha\in\{x,y,z\}}2t_{\alpha}\big{[}1-\cos(k_{\alpha}d_{\alpha})\big{]}-\mu_{\sigma},$$

$$2\sum _{\sigma ,\alpha ,\mathbf{k}}\left(t_{\alpha }{\hat{c}}_{\sigma \mathbf{k}}^{†}{\hat{c}}_{\sigma \mathbf{k}}\right)=2\sum _{\sigma ,\alpha }{t}_{\alpha }{N}_{\sigma }.$$

$$\hat{H}=\sum _{\mathbf{k}}$$

$$+\sum _{\mathbf{k}}\left({\xi }_{\downarrow \mathbf{q}-\mathbf{k}}-\frac{{\Delta }^2}{U}\right),$$

$$\displaystyle\sum_{\bf k}\bigg{\{}{\hat{c}}_{\uparrow{\bf k}}^{\dagger}{\hat{c}}_{\downarrow{\bf q-k}}^{\dagger}\Delta$$

$${\hat{H}}_{\mathrm{LO}}=$$

$$\displaystyle\times\left(\begin{array}[]{ccc}\xi_{\uparrow{\bf k}}&\Delta&\Delta\\ \Delta&-\xi_{\downarrow{\bf q-k}}&0\\ \Delta&0&-\xi_{\downarrow{\bf-q-k}}\\ \end{array}\right)\left(\begin{array}[]{c}{\hat{c}}_{\uparrow{\bf k}}\\ {\hat{c}}_{\downarrow{\bf q-k}}^{\dagger}\\ {\hat{c}}_{\downarrow{\bf-q-k}}^{\dagger}\end{array}\right)$$

$$\displaystyle\left(\begin{array}[]{c}\hat{\gamma}_{+,{\bf k,q}}\\ \hat{\gamma}_{-,{\bf k,q}}^{\dagger}\end{array}\right)=\left(\begin{array}[]{cc}u_{\bf k,q}&-v_{\bf k,q}\\ v_{\bf k,q}&u_{\bf k,q}\end{array}\right)\left(\begin{array}[]{c}{\hat{c}}_{\uparrow{\bf k}}\\ {\hat{c}}_{\downarrow{\bf q-k}}^{\dagger}\end{array}\right)$$

$$\displaystyle\left(\begin{array}[]{c}{\hat{c}}_{\uparrow{\bf k}}\\ {\hat{c}}_{\downarrow{\bf q-k}}^{\dagger}\end{array}\right)=\left(\begin{array}[]{cc}u_{\bf k,q}&v_{\bf k,q}\\ -v_{\bf k,q}&u_{\bf k,q}\end{array}\right)\left(\begin{array}[]{c}\hat{\gamma}_{+,{\bf k,q}}\\ \hat{\gamma}_{-,{\bf k,q}}^{\dagger}\end{array}\right).$$

$$E_{\pm ,\mathbf{k},\mathbf{q}}$$

Full text

\ioptwocol

1. 1 Introduction

2. 1.1 Superconductivity, superfluidity and magnetism

3. 1.2 Chandrasekhar-Clogston limit and exotic superfluidity

4. 1.3 FFLO studies with ultracold atomic gases

5. 1.4 Structure of this review

6. 2 Mean-field theory of the FFLO state in a lattice

7. 2.1 The Hamiltonian for a two-component Fermi gas in a lattice

8. 2.2 The BCS mean-field approximation

9. 2.3 The FF and LO mean-field approximation

10. 2.4 Bogoliubov transformation in the FFLO case

11. 2.5 The quasiparticle dispersions

12. 2.6 Self-consistent gap and number equations

13. 2.7 Thermodynamic potential and Helmholtz free energy

14. 2.7.1 Traps and Ensembles

15. 2.8 Ginzburg-Landau expansion

16. 2.9 The FFLO state and the Bloch’s theorem

17. 2.10 Luttinger’s theorem and Fermi surfaces

18. 2.11 Connection to the repulsive Hubbard model: the FFLO state and striped phases vs. d-wave superconductors

19. 3 Mean-field studies on FFLO in 2D and 3D optical lattices

20. 4 Mass-imbalanced Fermi gases in lattices

21. 5 FFLO state in the Green’s function formalism

22. 6 Collective modes for the FFLO state

23. 7 Beyond mean-field studies of the FFLO state in 2D and 3D lattices

24. 7.1 Dynamical mean-field theory (DMFT) for studies of the FFLO state

25. 7.2 The FFLO state in the 1D-3D crossover

26. 7.3 Quantum Monte Carlo studies of the FFLO state in 2D square lattices

27. 8 FFLO in trapped gases and phase separation

28. 8.1 Local density approximation

29. 8.2 Bogoliubov-de Gennes equations in spherical traps

30. 8.3 The FFLO state and proximity effects in trapped 3D gases

31. 8.4 FFLO in elongated traps and dimensional crossover

32. 8.5 BdG in practice

33. 9 FFLO and spin-orbit coupling

34. 9.1 Spin-orbit coupling in the context of cold atomic gases

35. 9.2 Spin-orbit coupling in fermionic systems and FFLO

36. 9.3 Topological states in systems with spin-orbit coupling

37. 9.4 Summary of the status of FFLO in the presence of spin-orbit coupling

38. 10 Detection of FFLO state

39. 10.1 Ramps and quenches

40. 10.2 Density profiles

41. 10.3 Interferometric methods

42. 10.4 Collective modes

43. 10.5 Noise correlations

44. 10.6 Spectroscopies

45. 11 Conclusions and outlook

1 Introduction

1.1 Superconductivity, superfluidity and magnetism

Superfluidity appears in range of systems from ultracold atomic gases to neutron stars, and from electrons in metals to quark-gluon plasma. The microscopic Bardeen-Cooper-Schrieffer (BCS) [1] theory explains superfluidity in fermionic systems as the formation and condensation of Cooper pairs of indistinguishable fermions, for instance particles with different internal spin states 111Throughout this review, by ”superfluid” and ”superfluidity” we refer to superfluidity of paired fermions, unless stated otherwise.. Magnetism and superfluidity of charged particles (i.e., superconductivity) are, to a large extent, mutually incompatible: magnetic fields cause eddy currents in superconductors, which either results in expulsion of the magnetic field from the superconductor known as the Meissner effect, or the suppression of superconductivity in the regions of non-zero magnetic field. If the magnetic field becomes sufficiently strong, the free energy associated with the superfluid state may exceed the free energy of the normal state and the material loses superconductivity altogether. This is because Cooper pairing favours equal amounts of particles in the two spin states while Pauli paramagnetism implies that the spins of the charge carriers prefer to align with the magnetic field. Still, it has been theorized that superfluidity and magnetism can coexist, and there has been tremendous interest in trying to realize such systems experimentally.

Notable theory proposals of ground states that combine superfluidity and magnetism, such as the Fulde-Ferrell-Larkin-Ocvhinnikov (FFLO) state [2, 3, 4], assume a finite Zeeman energy, but no orbital effects that are typically caused by magnetic fields. The search for such exotic forms of superconductivity and superfluidity thus calls for systems where orbital effects are absent. By considering two-dimensional layered superconductors for example, one can suppress orbital effects against in-plane magnetic fields, as circulating eddy currents around magnetic field lines are blocked by geometric constraints. This is the approach used in many solid-state experiments to realize magnetized superconductors [5, 6, 7, 8]. Another possibility to avoid orbital effects is to study superfluidity in neutral systems. Experiments on ultracold neutral atoms allow independent knobs for orbital effects (via fluid rotation) and the Zeeman effect (via pseudo-spin species imbalance). This provides a serious advantage of ultracold quantum gases (UQG) over solid state fermionic systems for exploring the competition between paired superfluidity and Zeeman energy. Here, we review the latest advances in the theoretical and numerical study of FFLO superfluids in neutral atomic gases. We focus on lattice systems but cover also recent developments on the effects of trap geometry and spin-orbit coupling.

1.2 Chandrasekhar-Clogston limit and exotic superfluidity

In spin-singlet superconductors, electrons are paired with opposite spins. Switching the electron spin would thus result in breaking a Cooper pair, imposing energy cost equal to the binding energy of the pair. This is known as the Chandrasekhar-Clogston limit [9, 10], which tells us that the Zeeman energy due to the magnetic field $h$ must exceed Cooper pair binding energy $\Delta$ before Pauli paramagnetism will play a role. A closer inspection yields the bound $h=\Delta/\sqrt{2}$. If the magnetic field is less than this critical value, the system will remain non-magnetised and a superfluid. Beyond this limit, the magnetic field provides sufficient energy to switch the spins, break the pairs, and suppress superfluidity. However this argument is limited to simple BCS-type superfluids. The salient question is whether some more exotic superfluid phases, such as breached pair [11, 12, 13], FFLO [2, 4] or deformed Fermi surface [14, 15] superfluids can exceed the Chandrasekhar-Clogston limit and realize magnetized superfluidity. Furthermore, the limit can be broken in spin-orbit coupled superfluids where spin-orbit coupling can stabilize the superfluid and enable different topological phases [16].

Breached pair superfluids retain the homogeneous nature of the BCS superfluid but allow the coexistence of spin-imbalance by opening an unoccupied breach in the momentum distribution. This results in gapless excitations, showing that not all atoms are paired. Deformation of Fermi surfaces, associated with spontaneous breaking of the rotational symmetry, has also been predicted as a possible mechanism for spin-imbalanced superfluidity with zero-momentum pairs [14, 15]. In the FFLO state, Cooper pairs are condensed in one or more finite-center-of-mass-momentum states. In a Fulde-Ferrell (FF) [2] superfluid, all pairs have the same momentum ${\bf q}$ while in the Larkin-Ovchinnikov (LO) [4] one pairs involve two opposite momenta $\pm{\bf q}$. Gapless excitations exist in FFLO superfluids like in the breached pair one, but in the FFLO state the order parameters become spatially varying, with the period given by the inverse of the FFLO wave vector ${\bf q}$. As condensates typically, the FFLO state is associated with breaking of the $U(1)$ symmetry for the phase of the order parameter $\Delta$, but in addition, other symmetries may be broken. The essential difference between FF- and LO-type states is that the FF states break the time-reversal symmetry but retain the translational one while LO states are time-reversal-symmetric but there is no translational invariance. The symmetry breaking properties are determined by having a complex (FF) or real (LO) order parameter, and by its spatial dependence and symmetry; the number of FFLO vectors $\mathbf{q}$ involved in the pairing can in general be more than one or two. The Goldstone modes that are related to breaking of continuous symmetries are thus distinct for different types of FFLO states.

1.3 FFLO studies with ultracold atomic gases

In UQG, the pseudospin states of the atoms are determined in the preparation of the sample, and these spin degrees of freedom are effectively frozen. The system is thus somewhat different from the Pauli paramagnetism in metals, where electrons will realign the magnetic moments with the external magnetic field. Instead, magnetism in ultracold gases is itinerant: atoms are mobile and free to rearrange into different magnetic domains. Any exotic magnetized superfluidity would appear as a coexistence region involving both superfluidity and spin imbalance, whereas absence of such exotic phases would be observable as either phase separation into a spin-balanced superfluid and spin-imbalanced normal fluid, or as an absence of superfluidity altogether.

A phase with spin-imbalance and vanishing Cooper pair momentum ${\bf q}=0$ is known to be unstable at low temperatures, with instability in favor of the FFLO state. On the other hand the FFLO state is unstable towards the aforementioned phase separation in a large region of parameter space, in particular in three-dimensional (3D) systems. Indeed early UQG experiments on spin-imbalanced superfluids [17, 18] showed phase separation with a spin-balanced superfluid core surrounded by a spin-imbalanced normal state gas in a harmonic trap.

There are already reviews regarding the FFLO state in continuum 3D UQG [19, 20, 21]. Reduced dimensionality has been predicted to favour FFLO. This was exemplified by promising experimental results in references [22, 23], in which a spin-imbalanced fermionic gas was prepared in an array of one-dimensional (1D) tubes. Density profiles were in good agreement with theories and suggestive of an FFLO state, but clear observation of FFLO-type pairing was not obtained. These experiments were done in continuum and are thus outside the scope of this review, but they have been covered in a review [24]. In low dimensions, the long-range order of the FFLO state is hampered by quantum and thermal fluctuations. Even in 3D the nature of the FFLO state is quite peculiar, with the FF state having vanishing superfluid stiffness in transverse (to the FFLO momentum) direction and in the LO state the superfluid density being highly anisotropic [19]. The small parameter regime where the FFLO state is predicted in continuum systems, together with vulnerability to phase separation and fluctuations, motivates the search for systems and conditions that could improve the stability of this elusive state of matter.

The focus of this review is on the FFLO state in 2D and 3D optical lattices. We mention some works of FFLO in 1D lattices but do not provide a comprehensive review of the 1D case, for that see [25, 24, 26]. Lattice systems improve the stability of the FFLO state due to the shape of Fermi surfaces that can increase the available phase space for pair formation for finite-momentum pairs. This effect is called Fermi surface nesting, and is closely connected with existence of van Hove singularities that provide high density of states. The nesting effects are further enhanced in lattices with a lower dimension. Low-dimensional optical lattices are thus promising for experimental realization of the FFLO state. Two-dimensional (2D) lattices are particularly attractive with the advent of quantum gas microscopes for fermionic atoms [27, 28, 29, 30, 31, 32], which allow observation of individual atoms and probing of non-local correlations. Mott [33, 34] and short-range antiferromagnetic [35, 36, 37, 38] correlations have been observed for repulsively interacting fermions in optical lattices, and recently the use of quantum gas microscopes has enabled observing antiferromagnetic correlations ranging over the whole sample [39]. These are closely connected with the FFLO state with attractive interactions through particle-hole mapping: the striped phase in the doped repulsive Hubbard model corresponds to the FFLO state in the attractive Hubbard model. This striped phase may compete with d-wave superfluidity, which would be important for understanding the origin of high temperature superconductivity. Quantum gas microscopes provide a promising platform for studying these phenomena and for the actual observation of the FFLO state.

An interesting new candidate for realizing FFLO superfluids is spin-orbit coupled atomic gases. UQG systems with spin-orbit coupling have recently been created for bosons [40, 41, 42] and for fermions [43, 44]. These systems are realized by utilizing Raman lasers that couple the spin- and spatial degrees of freedom of the atoms, effectively providing a ${\bf L\cdot S}$-type spin-orbit coupling (SOC). This is a neutral atom analogy of the spin-orbit coupling in solid state physics or of the electrons orbiting the nucleus within the atom. In SOC systems, the spatial degrees of freedom become mixed with the spin degrees of freedom, making the concept of magnetization less well defined. Like optical lattices, these systems also exhibit large parameter regions with FFLO-type superfluids, but here finite-momentum pairing can already be realized with spin-balanced systems. The SOC introduces preferred momenta in the dispersion, moreover, some forms of SOC may break the rotational symmetry already at the Hamiltonian level. Furthermore, unlike in non-spin-orbit coupled systems, even the single-plane wave FF state has been shown to be realizable in the ground state. In sum, the FFLO state in SOC systems is expected to be quite different from its conventional counterpart in continuum and lattices.

1.4 Structure of this review

We start this review by providing in section 2 a detailed pedagogic description of the mean-field theory applied to spin-imbalanced gases leading to FFLO states and other spin-imbalanced superfluid states. We describe the use of thermodynamic potentials and Ginzbug-Landau theory for determining the stability of the various phases and nature of the phase transitions. Important theorems, such as the Bloch’s theorem for the absence of current in the ground state and the Luttinger’s theorem connecting Fermi surfaces with particle densities, are discussed. The mapping to the repulsive Hubbard model and the corresponding connection of the FFLO state to the striped phase and d-wave superconductors is explained.

We continue with mean-field studies in section 3 by providing a review of theoretical mean-field works on the FFLO state in optical lattices and explaining the important role played by Fermi surface nesting and van Hove singularities. Also works on mass imbalanced gases in optical lattices are briefly reviewed in section 4.

As a step toward a quantum field theoretical treatment, section 5 describes the FFLO physics using Green’s function formalism. The formalism is applied in section 6 for studying collective modes. In section 7 we review the studies done with more advanced numerical methods, namely dynamical mean-field theory (DMFT) and quantum Monte Carlo (QMC) studies. These studies are important in determining the actual stability of the FFLO state against fluctuations neglected in mean-field studies. Especially the role of dimensional crossover is considered, as low dimensional systems provide larger parameter ranges for the FFLO state while higher dimensions give the necessary stability for the long-range order. Despite the focus on lattice systems, we present also a review on trapped continuum systems in section 8. We concentrate on advances missing from the previous review [21], particularly new studies done in non-spherical traps.

We give an analysis of the literature on the FFLO state in spin-orbit coupled systems in section 9. In section 10, we provide a comprehensive review of the various detection schemes that have been proposed for observing the FFLO state in continuum systems, optical lattices, various dimensionalities, and spin-orbit coupled systems. This is a stand-alone part of the review and can be accessed without first reading the rest of the review. Finally, in section 11 we conclude with an overview of the salient points of the latest research on FFLO physics in UQG, and of the key open questions.

2 Mean-field theory of the FFLO state in a lattice

Here, we present the basic mean-field theory description of the FFLO state using the approach of BCS theory and the Bogoliubov transformation. We start from a generic Hamiltonian for two component fermions, which describes an ultracold Fermi gas well. A lattice potential is included and continuous position variables are used to begin with, discrete spatial variables are introduced later. The reader can easily convert the following lattice description back to free space, if desired. We perform the BCS mean-field approximation, including the possibility of a population imbalance between the two spin components and exotic pairing, such as the FFLO and or breached pair (BP) states. The Hamiltonian is expressed in the (lattice) momentum basis and a Bogoliubov transformation is applied to diagonalize the Hamiltonian. The grand potential and the Helmholtz free energy can then be determined in terms of the eigenenergies. We discuss those quasiparticle eigenenergies in detail since they provide intuitive understanding about the various exotic superfluids. Finally, the fundamental relations between the FFLO state and Bloch’s and Luttinger’s theorems are described.

2.1 The Hamiltonian for a two-component Fermi gas in a lattice

We present here the Hamiltonian for a two-component Fermi system, including a lattice potential, and show how to obtain the lattice Hamiltonian from the continuous one. For a more detailed derivation and information on the approximations used, see for instance [45, 46, 47, 48, 49]. We consider a system of two types of fermions, which we label spin up ($\uparrow$) and spin down ($\downarrow$). The main focus are systems of ultracold gases, where this labeling is a pseudospin, which means atoms in two different internal atomic states, such as hyperfine states, or different atomic species or isotopes. The interaction between the two components of the gas is described by the potential $V_{\uparrow\downarrow}\left(\mathbf{r},\mathbf{r}^{\prime}\right)$, and the components experience an external lattice potential $V_{\mathrm{ext},\sigma}(\mathbf{r})$, which could include also trapping potential. We consider three dimensional systems, but transferring these results to other dimensions is straightforward. The system is described by the Hamiltonian

[TABLE]

where ${\hat{K}_{\sigma}}({\mathbf{r}})=-\frac{\hbar^{2}\nabla^{2}}{2m_{\sigma}}+V_{\mathrm{ext},\sigma}({\bf r})$, $m_{\sigma}$ is the mass and $\mu^{0}_{\sigma}$ the chemical potential for a particle in state $\sigma$. A lattice potential for ultracold atoms is created by interfering laser beams, where the periodicity of the lattice is half the laser wavelength and the lattice depth depends on the intensity of the laser [50, 46]. Nowadays, a lattice can also be created using holographic masks [51] and such a lattice was used for a single site resolution of individual atoms. For other experiments on quantum gas microscopes, see references [52, 53, 29, 28, 32]. For ultracold dilute gases, where the interactions are short range and temperatures low, the interaction potential can be approximated by a contact interaction, see for instance references [46, 48, 54, 49],

[TABLE]

where $V_{0}<0$ for attractive interactions. The Hamiltonian then reads

[TABLE]

In the presence of a lattice potential, it is convenient to expand the Fermi operators in the continuum as

[TABLE]

where ${\hat{c}}^{\dagger}_{\sigma{\bf n,i}}$ creates a $\sigma$ particle in the lattice at site ${\bf i}$ in an energy band ${\bf n}$ and where $w^{*}_{\bf n}({\bf r-r_{i}})$ are the Wannier functions [55]. In the tight-binding limit, the Wannier functions can be approximated by harmonic oscillator states on each site. At low temperatures and for sufficiently small interaction energies, which are criteria met by ultracold atoms in optical lattices, particles only occupy the lowest band with ${\bf n}=0$ and consequently the band index can be omitted. Note, however, that if one wishes to consider very strong interactions created by Feshbach resonances, a proper description involves higher bands [56]. Using the above expansion, the continuum Hamiltonian can be rewritten into the standard Hubbard Hamiltonian

[TABLE]

where the summation indices $\langle{\bf i,j}\rangle$ are over nearest neighbour indices and nearest neighbour tunneling (hopping) energy $t_{\sigma{\bf ij}}$, the on-site energy $\varepsilon_{\sigma}$, and the on-site interaction energy $U$ are [57]

[TABLE]

In the equations (5)–(8), we have allowed for spin dependent and anisotropic hopping $t_{\sigma{\bf ij}}$. Anisotropic lattices can be created by using different laser intensities and wavelengths for different spatial directions. A spin dependence can result both from different masses for the two spin components or from a spin-dependent lattice. In both cases, the onsite energies are also spin dependent $\varepsilon_{\uparrow}\neq\varepsilon_{\downarrow}$. A spin dependent lattice can be created by applying a magnetic field gradient modulated in time [58].

In the rest of this section, we assume the masses of the two spin components to be the same and the lattice potential to be spin independent, so that $t_{\uparrow{\bf ij}}=t_{\downarrow{\bf ij}}\equiv t_{{\bf ij}}$ and $\varepsilon_{\uparrow}=\varepsilon_{\downarrow}\equiv\varepsilon$. The latter is in this case a constant energy shift and can be taken equal to zero. The Hamiltonian now reads

[TABLE]

where we use the same notation as in the continuum for the position ${\bf r}$, which is here a discrete variable corresponding to the lattice sites, and here $U<0$ for attractive interactions.

2.2 The BCS mean-field approximation

To study the transition from a normal state to a superfluid consisting of paired fermions, the so-called Cooper pairs, we introduce a field describing the pairs. This is done using the BCS mean-field approximation, which also introduces the Hartree fields and replaces, through the use of Wick’s theorem, some of the operators in the four-operator interaction term by complex numbers. More details on BCS theory in general and in cold gases can be found in numerous places, for instance references [59, 60, 54, 49].

The order parameter of the BCS theory is the pairing field $\hat{\Delta}({\bf r})$, which is on average related to the expectation value of $\hat{\psi}_{\downarrow}\hat{\psi}_{\uparrow}$

[TABLE]

and describes the Cooper pairs. As we see in section 2.4, the pairing field is related to the energy gap in the excitation spectrum of the BCS theory. The usual BCS theory describes transitions from a normal gas of fermionic particles to a homogeneous superfluid consisting of Cooper pairs, which is described by a mean-field ansatz independent of position $\Delta({\bf r})=\Delta$. However to describe exotic superfluidity, the spatial dependence is a key feature.

In the lattice model, the pairing field is

[TABLE]

We include the Hartree fields by redefining the chemical potentials. Denoting $\left\langle{\hat{c}}_{\sigma}^{\dagger}{\hat{c}}_{\sigma}\right\rangle=n_{\sigma}$, we get

[TABLE]

where it is assumed that the densities are uniform. In the mean-field approximation, the pairing field is approximated by its average value, $\hat{\Delta}({\bf r})=\Delta({\bf r})$, but we do not yet assume any specific form for this average. The Hamiltonian now reads

[TABLE]

where we also included the Hartree shifts. Note that the assumption of uniform densities is not always strictly correct because a spatial modulation of the order parameter amplitude $|\Delta({\bf r})|$ is typically accompanied by a modulation of density in the FFLO state. To take this into account one should keep the spatial dependence of the Hartree terms and solve the problem self-consistently for instance by a Bogoliubov-de Gennes (BdG) approach. In practice this is quite difficult, and fortunately (for the theorist, not for the experimentalist), the modulations of the density are quite small so one can justify the uniform density approximation for the Hartree shifts in the Hamiltonian (13). One should keep in mind, however, that if such Hamiltonian then leads to modulated density of the ground state, the treatment is not fully self-consistent.

2.3 The FF and LO mean-field approximation

The possibility of pairing with unequal chemical potentials and with a spatially dependent order parameter is included in the general pairing order parameter introduced in equation (11). In the FF state, $\uparrow$ and $\downarrow$ particles in different momentum states form a Cooper pair, such that the Cooper pair carries a net momentum ${\bf q}$. This is in contrast to the BCS superfluid, where the Cooper pairs are formed between particles in momentum states ${\bf k}$ and ${\bf-k}$. When describing the FF superfluid state, the Cooper pair momentum ${\bf q}$ is a parameter determined by minimizing the system’s energy, and its absolute value turns out to scale with the population imbalance between the two spin components $|\mathbf{q}|\propto|k_{F\uparrow}-k_{F\downarrow}|$ 222Note that the relation $|\mathbf{q}|\propto|k_{F\uparrow}-k_{F\downarrow}|$ is given by mean-field theory, and some beyond mean-field studies on finite lattices have found FFLO wavevectors deviating from the mean-field results [61].. The FF order parameter is

[TABLE]

where $\Delta\in\mathbb{R}$ is a real number. This plane wave form of the order parameter was first considered by Fulde and Ferrell [2]. Larkin and Ovchinnikov [4] around the same time considered an ansatz containing both ${\bf q}$ and $-{\bf q}$, so that the order parameter is of the form

[TABLE]

In the FF phase, the order parameter is complex as equation (14) shows, and time-reversal symmetry is broken. The term FF phase (or FF state) can be used to denote also time-reversal symmetry broken states with several pairing wave vectors ${\bf q}$. Likewise, the essential feature of the LO state is the preservation of the time-reversal symmetry; two opposite FFLO wave vectors ${\bf q}$ and $-{\bf q}$ is the minimal setting for this, but there could be more.

According to many studies, both for homogeneous and for lattice systems, the cosine wave LO ansatz gives a lower energy than the FF plane wave, see for instance references [62, 63, 64, 65, 66]. Fortunately in many different systems, the use of FF and LO ansatzes predict the same parameter regime for the existence of the density-imbalanced superfluid. This means that one can use the simpler FF ansatz to explore different systems and parameter regimes, keeping in mind that the true ground state could be the LO state. However even when the parameter regimes of the FF and LO state overlap, the order of the phase transitions predicted by using an FF or LO ansatz can be different, for instance compare references [67, 68] and reference [66].

Since the LO state is in general more stable than the FF state, loosely speaking we can also say that the FFLO state is described by the FF ansatz in equation (14). However, there are also systems where the above is simply not true, for example for Fermi gases with spin orbit coupling. In these systems, it is actually found that the FF state is energetically more favorable than the LO state, see section 9. Throughout this review as well as often in the literature, the term FFLO superfluid is used even when only the FF ansatz is actually applied. In general, the FFLO order parameter can be composed of any number of Fourier components with different ${\bf q}$s, leading to arbitrarily complicated spatial dependence $\Delta({\bf r})$. This is discussed further in section 8 where the Bogoliubov-deGennes formalism is introduced.

Let us now proceed here with the FF ansatz (14), so that the Hamiltonian becomes

[TABLE]

It is convenient to write the Hamiltonian in the (lattice-)momentum representation, that is, represent the operators in the plane wave basis. This means Fourier transforming equation (2.3) by using

[TABLE]

where $M$ is the number of lattice sites and the momentum ${\bf k}$ runs through the reciprocal lattice. The density terms transform as

[TABLE]

since $\frac{1}{M}\sum_{{\bf r}}e^{i({\bf k}^{\prime}-{\bf k})\cdot{\bf r}}=\delta_{{\bf k},{\bf k}^{\prime}}$ for the reciprocal lattice vectors ${\bf k},{\bf k}^{\prime}$. Similarly, the pairing terms become

[TABLE]

We now assume a cubic lattice, where the nearest neighbor hopping term gives rise to a cosine dispersion, for example, the Fourier transform for the $x$-direction is

[TABLE]

where $d_{x}$ is the lattice spacing in the $\hat{x}$-direction.

Finally, we arrive at the mean-field FFLO Hubbard Hamiltonian in momentum space,

[TABLE]

where the single particle energies are given by

[TABLE]

with $t_{\alpha}$ the tunneling coefficient for hopping in the direction $\alpha$. In the case that all particles in the system occupy the lowest energy levels $\xi_{\sigma{\bf k}}$, the highest occupied energy state is the Fermi energy $E_{\mathrm{F}}$, and the corresponding momentum state is the Fermi momentum $k_{\mathrm{F}}$. At zero temperature, the chemical potential is equal to the Fermi energy for a non-interacting system. To get the dispersions to correspond to those of free particles in the limit of small $k$, the following terms have been added to the Hamiltonian:

[TABLE]

2.4 Bogoliubov transformation in the FFLO case

In this subsection, we diagonalize the Hamiltonian in equation (2.3) by a Bogoliubov transformation to be able to determine the number and gap equations, and to calculate the thermodynamic potential. By performing this transformation, from the particle basis to the quasiparticle basis, we find the quasiparticle energies. Let us start by rewriting the Hamiltonian equation (2.3) using matrix formalism

[TABLE]

where we have used that ${\hat{c}}_{\downarrow{\bf k}}^{\dagger}{\hat{c}}_{\downarrow{\bf k}}=1-{\hat{c}}_{\downarrow{\bf k}}{\hat{c}}^{\dagger}_{\downarrow{\bf k}}$. The essential point in this representation is that the Hamiltonian can be expressed as a sum of independent $2\times 2$ matrices for each lattice momentum ${\bf k}$ that can be diagonalized separately. In other words, if the above summation over momentum ${\bf k}$ is written explicitly, the Hamiltonian can be expressed with a (large) block-diagonal matrix, which allows for rewriting the Hamiltonian as a sum of $2\times 2$ matrices.

This simplicity results from the plane wave ansatz used for the FF state. If a Cooper pair ansatz is considered including more plane waves, it is no longer possible to rewrite the Hamiltonian using $2\times 2$ matrices. This is already the case for the LO ansatz in equation (15), containing two plane waves. Namely, the pairing part of the Hamiltonian for the LO Cooper pairs reads

[TABLE]

which contains more terms than the FF pairing Hamiltonian equation (2.3), since in the LO case an $\uparrow$ particle in momentum state ${\bf k}$ couples to both a $\downarrow$ fermion with momentum ${\bf q-k}$ as well as to one with ${\bf-q-k}$. For rewriting the complete LO Hamiltonian using matrix multiplication, it means that the resulting matrix is not block-diagonal and cannot be split into smaller independent matrices. Instead, even with the summation over momentum states ${\bf k}$, the LO matrix is limited only by the number of lattice sites and therefore diverges for infinite lattices. We truncate this matrix here to show a small part

[TABLE]

where the factor 2/3 is added to compensate for the fact that two spin components are here described by a $3\times 3$ matrix. In the above matrix (33), all couplings from the ${\hat{c}}_{\uparrow{\bf k}}$ field to the $\downarrow$ fields are included, but not all couplings are in turn included for the $\downarrow$ fermions. For instance according to equation (31), ${\hat{c}}_{\downarrow{\bf q-k}}$ couples to both ${\hat{c}}^{\dagger}_{\uparrow{\bf k}}$ and ${\hat{c}}^{\dagger}_{\uparrow{\bf k-2q}}$. The latter coupling is not included in the above matrix.

To determine thermodynamic quantities for the LO state, one needs to diagonalize the Hamiltonian using a Bogoliubov transformation. This is done in the same way as explained in the following for the FF state, except that it can only be done numerically. In practice, thermodynamic quantities are calculated for a certain matrix dimension, and consequently the matrix is enlarged by including more coupling terms until convergence is reached. When describing a continuous, second order phase transition from the normal state to the superfluid phase the order parameter $\Delta$ is small, and a $6\times 6$ or $10\times 10$ matrix could be large enough (for LO phase breaking translational symmetry in one out direction only). However when describing the superfluid phase or a first order phase transition, much bigger matrix sizes are needed. For more details, see references [65, 66].

Let us now turn back to the FF Hamiltonian and apply a Bogoliubov transformation for the matrices in (30), by introducing the new basis operators $\hat{\gamma}_{+}$ and $\hat{\gamma}_{-}$

[TABLE]

and we require $u_{\bf k,q}^{2}+v_{\bf k,q}^{2}=1$ to ensure the fermionic anticommutation relations for the $\hat{\gamma}$ operators. The inverse of the above transformation is

[TABLE]

By inserting this transformation in equation (30) and demanding that the resulting Hamiltonian is diagonal, we arrive at the quasiparticle energies

[TABLE]

For comparison, the quasiparticle energies for the standard BCS theory with ${\bf q}=0$ read

[TABLE]

When the chemical potentials are equal for the two spin components, $\mu_{\uparrow}=\mu_{\downarrow}$, then $E_{+}=E_{-}$ and thier minimum value is $\Delta$. This shows that the quasiparticle excitations have a minimum energy, that is, the energy gap, the magnitude of which is $\Delta$.

The Bogoliubov coefficients are found by diagonalizing the Hamiltonian:

[TABLE]

[TABLE]

and we arrive at

[TABLE]

where in the last step we used the anticommutation relation $\gamma_{-{\bf k}}\gamma_{-{\bf k}}^{\dagger}=1-\gamma_{-{\bf k}}^{\dagger}\gamma_{-{\bf k}}$.

The operators $\hat{\gamma}_{\pm}$ are the quasiparticle operators that obey fermionic anticommutation relations, like the original particle operators. Within the BCS theory, quasiparticles are non-interacting. To derive some basic properties of the superfluids such as collective modes, quasiparticle interactions have to be introduced via the generalized random phase approximation (GRPA), for instance, as discussed in section 6. For the moment, we consider non-interacting quasiparticles, which means that the system is described as an ideal Fermi gas with dispersions given by $E_{+}$ and $E_{-}$ in equation (55).

2.5 The quasiparticle dispersions

From the diagonalized Hamiltonian, equation (2.4), the mean-field thermodynamic potential can be calculated. Consequently the phase diagram of the two-component Fermi gas can be determined from the thermodynamic potential by minimizing it with respect to the gap $\Delta$ and the FFLO wavevector ${\bf q}$. The thermodynamic potential depends on the quasiparticle dispersions $E_{\pm}$, and before turning to the calculation of the thermodynamic potential, we therefore take a closer look at the quasiparticle dispersions.

The quasiparticles describe the creation of excitations in the BCS state. Explicitly, the operator $\hat{\gamma}_{+}$ is a superposition of creating an up ($\uparrow$) particle in some momentum state and at the opposite momentum creating a hole for the down ($\downarrow$) component. Similarly, the quasiparticle $\hat{\gamma}_{-}$ creates a $\downarrow$ particle and an $\uparrow$ hole. The BCS state has either both $\uparrow$ and $\downarrow$ components occupying a certain $({\bf k},-{\bf k})$ pair of states, or no particles occupying these states. Therefore an excitation, that is a single $\uparrow$ particle at state $\vec{k}$ but no $\downarrow$ particle at $-{\bf k}$, can be created either by destroying the $\downarrow$ particle in the case where initially there were particles in the $({\bf k},-{\bf k})$ states, or by creating an $\uparrow$ particle in the case where the states $({\bf k},-{\bf k})$ were empty. This is why the quasiparticle operators are superpositions of particles and holes.

When the system is in the normal state, there are no Cooper pairs, $\Delta=0$. The quasiparticle dispersions are then equal to the particle dispersions, $E_{+,{\bf k}}=\xi_{\uparrow{\bf k}}$ and $E_{-,{\bf k,q}}=\xi_{\downarrow{\bf q-k}}$, for momenta above the Fermi levels. For $k>k_{\mathrm{F}}$, one can then associate $E_{+,{\bf k}}$ with an $\uparrow$ particle and $E_{-,{\bf k}}$ with a $\downarrow$ one, see the left panel in figure 1. There also the hole branch is shown, which is $-E_{\pm}$. As soon as the Fermi levels for the two spin components are unequal, the association of the quasiparticles with the particles becomes more complicated due to the square root term in equation (55).

In the two-component Fermi mixture, a phase transition can occur to a superfluid, consisting of condensed Cooper pairs. For an equal number of particles, $\mu_{\uparrow}=\mu_{\downarrow}$, this transition in the square lattice is toward a BCS superfluid with ${\bf q}=0$. The quasiparticles for this case are presented in the right panel in figure 1, where $\Delta\neq 0$. In the paired case, an excitation gap is formed. Here the two eigenenergies are degenerate, meaning that $E_{+}=E_{-}$.

If there is a chemical potential imbalance in the Fermi mixture, $\mu_{\uparrow}\neq\mu_{\downarrow}$, a phase transition to a superfluid phase with ${\bf q}=0$ could in principle occur, see figure 2. With respect to the degenerate case in figure 1, the quasiparticle dispersions are shifted by the difference in chemical potentials, $h=(\mu_{\uparrow}-\mu_{\downarrow})/2$. In the left figure, $h<\Delta$ and both excitations are gapped, corresponding to the BCS superfluid. In the right figure, the chemical potential imbalance is even bigger and $h>\Delta$. Now $E_{+}$ crosses zero, and there is a range of momenta for which gapless excitations exist in $E_{+}$. This solution of the gap equation is has the peculiar feature that even when most particles (momentum states) are associated with the pairing gap $\Delta$, there are momentum values where gapless excitations are possible. In contrast, in the BCS state all eigenvalues are positive, and there are no quasiparticles at zero temperature. It is important to note that, for fixed chemical potentials, weak interactions, continuum systems or simple square lattices, this non-BCS solution of the gap equation is not a state of matter since it is a maximum, not minimum of the grand potential. Sarma explained this in his original work [11]. This solution of the gap equation was later suggested as ”breached pair (BP) state” [12], although the authors soon noted that one needs to have fixed particle numbers, not fixed chemical potentials, in order for it to be a state of matter (minimum of energy) [13]. Later it was shown that the same argument about fixed particle numbers versus chemical potentials holds in square lattices [69]. The terms Sarma state and BP state are nowadays used for referring to this type of solution of the gap equation, which under specific conditions (e.g., fixed particle numbers or complex lattice geometry) can also correspond to a stable state of matter. In the following, we mostly use the terms Sarma phase or Sarma state.

The phase transition from the normal state is can also be toward an FFLO superfluid phase, where the Cooper pairs carry a net momentum ${\bf q}$. The quasiparticle dispersions for that case are shown in figure 3. It turns out that the wave vector ${\bf q}$ found by minimizing the thermodynamic potential is proportional to the population imbalance $|{\bf q}|\simeq|k_{F\uparrow}-k_{F\downarrow}|$. In other words, the bigger the difference between the chemical potentials or particle numbers, the bigger $|{\bf q}|$. Also in the FFLO case, there are gapless excitations, and quasiparticles can exist in the superfluid state. The difference with the Sarma state is that the two quasiparticle branches $E_{+}$ and $E_{-}$ are not only shifted with respect to each other as in figure 2, but they also become asymmetric with respect to $k$ and $-k$ due to the nonzero ${\bf q}$, see figure 3.

2.6 Self-consistent gap and number equations

It is now possible to derive a set of coupled equations consisting of the gap and number equations, from which the gap parameter $\Delta$ and the chemical potentials, $\mu_{\uparrow}$ and $\mu_{\downarrow}$, can be solved. Let us start with deriving the number equations. The total number of particles in one spin component is given by

[TABLE]

Using the inverse of the Bogoliubov transformation equation (54), we obtain for the $\uparrow$ component

[TABLE]

where the fact that the expectation values of the form $\langle\hat{\gamma}_{+}\hat{\gamma}_{-}\rangle$ are zero for non-interacting (quasi)particles [59, 48, 49] was used. The chemical potential of the $\uparrow$ component depends on the number of $\downarrow$ particles via the Hartree shifts, see equation (12). A similar equation holds for the number of down ($\downarrow$) particles:

[TABLE]

where the chemical potentials contain the Hartree shifts.

Using the original definition of the order parameter, $U\langle{\hat{c}}_{\downarrow{\bf r}}{\hat{c}}_{\uparrow{\bf r}}\rangle=\Delta e^{i{\bf q\cdot r}}$, it is possible to derive the so-called gap equation. The BCS gap equation, where ${\bf q}=0$, is

[TABLE]

After dividing the above equation by $\Delta$ and subsequently using equation (58), the gap equation reads

[TABLE]

where the prefactor on the right side is positive, since $U<0$ for attractive interactions.

The summation on the right side of the above gap equation equation (70) diverges at high momenta ${\bf k}$. This divergence directly results from the approximation of the interatomic interaction by the unphysical contact interaction and is absent if the true interatomic potential is used. In the case we consider here, the single band lattice model, this divergence does not cause any issues because the single band approximation provides a natural cutoff for the above summation, namely the band width (assuming that interactions are so weak that the single band approximation is valid). In the homogeneous Fermi gas, without a lattice, the ultraviolet divergence can be removed, for example, by introducing a cutoff and renormalizing the coupling constant, see for instance references [46, 48, 49].

This coupled set of equations, consisting of the number equations (67) and (68), and the gap equation (70), is equivalent to the standard BCS-Leggett theory in a lattice (although there often the chemical potentials of the components are assumed equal). In the population balanced case, where $N_{\uparrow}=N_{\downarrow}$, it is possible to choose $\mu_{\uparrow}=\mu_{\downarrow}$ by hand and eliminate one of the number equations. If the interactions are weak, that is $|U|$ is small, the chemical potential can be approximated with the Fermi energy. However when the interaction strength increases, the chemical potential has to be solved from the number equation to get the correct results. The interactions become large when exploring the BCS-BEC crossover [70, 54, 71] and approaching the BEC regime from the BCS one, then the use of the number equation becomes mandatory. Therefore the combined gap and number equations are sometimes called crossover equations. It is essential to remember that the gap equations have also solutions that are not minima of energy and are thus unphysical. It is always essential to examine the physical nature of the solutions, which can be done for instance by considering the thermodynamic potential.

2.7 Thermodynamic potential and Helmholtz free energy

To determine the ground state of a two-component Fermi mixture in a lattice for certain system parameters, such as densities and temperatures, one needs to minimize the total energy of the system. Depending on the ensemble considered, this total energy is either the grand thermodynamic potential or the (Helmholtz) free energy. The minimum of the grand thermodynamic potential $\Omega$ corresponds to thermodynamic equilibrium in the grand-canonical ensemble, where the chemical potentials are fixed and the particle numbers are fluctuating. The thermodynamic potential is the Legendre transform of the free energy $F$. In the canonical ensemble, where the particle numbers are fixed, thermodynamic equilibrium is reached by minimizing the free energy. Here we show a short derivation of the thermodynamic potential $\Omega$. For a more elaborate derivation and more background on ensembles and Legendre transforms, see for instance references [59, 60, 70, 72, 48].

To derive the thermodynamic potential, from which subsequently other thermodynamic quantities can be obtained, we first need the partition function for the grand-canonical ensemble, which can be calculated from the diagonalized Hamiltonian,

[TABLE]

where the sum goes through the quasiparticle basis (i.e., the basis where the quasiparticle occupation numbers are good quantum numbers) of the Hilbert space and $\beta=1/k_{\mathrm{B}}T$, the inverse thermal energy with $k_{\mathrm{B}}$ Boltzmann’s constant. The Hamiltonian in equation (2.4) is of the form

[TABLE]

where $C=\sum_{{\bf k}}\left(\xi_{\downarrow-{\bf k+q}}-E_{-,{\bf k,q}}-\Delta^{2}/U\right)$. As a constant, $C$ is factored out of the expectation value and $Z$ reads

[TABLE]

where $\hat{n}_{\sigma,{\bf k}}=\hat{\gamma}_{\sigma,{\bf k}}^{\dagger}\hat{\gamma}_{\sigma,{\bf k}}$, with $\sigma=\pm$. Since all commutators of the type $[\hat{n}_{+,{\bf k}},\hat{n}_{\pm,{\bf k^{\prime}}}]$ are zero for all ${\bf k,k^{\prime}}$, the exponential of the sum can be written as the product of exponentials,

[TABLE]

For fermionic quasiparticles, the possible occupation numbers for a given state are [math] and $1$, leading to

[TABLE]

The thermodynamic grand potential can now be obtained via

[TABLE]

and reads in this case

[TABLE]

In thermodynamic equilibrium, the order parameter takes the value that minimizes the thermodynamic potential. Here the order parameter is 4-dimensional because $\Delta$ is the amplitude of the order parameter and ${\bf q}$ gives the magnitude and direction of the wave vector in the plane wave form. To find the equilibrium state of the system and determine the phase diagram, one thus has to minimize the thermodynamic potential. When the global minimum of the thermodynamic potential is located at $\Delta=0$ and ${\bf q}=0$, it means there are no Cooper pairs and the system is in the normal phase. The system is in the BCS state when the particle densities are equal $n_{\uparrow}=n_{\downarrow}$ and the global minimum of $\Omega$ is attained at $\Delta\neq 0$ and ${\bf q}=0$. The Sarma phase would form the ground state of the system if the energy was minimized with a nonzero number of Cooper pairs carrying no momentum, $\Delta\neq 0$ and ${\bf q}=0$, in the presence of a population imbalance $n_{\uparrow}\neq n_{\downarrow}$. Finally when the global minimum of the thermodynamic potential is located at $\Delta\neq 0$ and ${\bf q}\neq 0$, the system is in the FFLO state.

At zero temperature, the BCS state can exist as long as the chemical potential difference between components does not exceed the so-called Chandrasekhar-Clogston limit. At weak coupling, this implies that the chemical potential difference is smaller than $\Delta/\sqrt{2}$ [73]. At finite temperature, a sufficiently small density difference may coexist with the superfluid since thermal quasiparticles can accommodate some imbalance even in the BCS state, see figure 2a) where one of the quasiparticle branches has a smaller gap and can thus have larger thermal population. In the Sarma state, imbalance can be accommodated both at zero and finite temperature because some of the quasiparticle energies take negative values as seen from figure  2b). Equation 56 shows that existence of negative quasiparticle energies requires $\Delta$ to be smaller than $(\mu_{\uparrow}-\mu_{\downarrow})/2$; this condition should be valid for the Sarma/BP state. The finite temperature BCS state with imbalance and the Sarma/BP state at finite temperature both display uniform superfluid order parameter and spin-density imbalance, therefore these terms are sometimes used rather loosely and interchangeably in the literature. However, the absence and presence of a Fermi surface is a fundamental difference between these states.

Alternatively, the above mean-field thermodynamic potential can be calculated using a path-integral formalism and subsequently making a saddle-point approximation. The criteria obtained for phase transitions are the same, see references [74, 75] and the references therein.

2.7.1 Traps and Ensembles

It is important to choose correctly whether to minimize the thermodynamic potential or the free energy, based on which system parameters are fixed. In our cases of interest, the fixed parameters are temperature, and depending on the situation, either the chemical potentials $\mu_{\uparrow}$ and $\mu_{\downarrow}$ or the total particle numbers, $N_{\uparrow}$ and $N_{\downarrow}$, as explained below.

In ultracold gas experiments, there is typically no particle bath. Rather the numbers of particles in the total system, $N_{\uparrow}$ and $N_{\downarrow}$, are constant in the experiment. Furthermore since the time scale of spin relaxation is much larger than the time scale of the experiment, the atoms are prevented from converting their spin. However in a harmonic trap, the particles can distribute themselves non-uniformly, and therefore the local particle number is not fixed by external constraints. In this situation, the atom distribution assumes a form where the chemical potentials are constant throughout the trap and the grand potential $\Omega(\Delta,{\bf q},\mu_{\uparrow},\mu_{\downarrow})$ is the relevant energy to be minimized. This is usually the case in studies of trapped ultracold gases, also in the case of optical lattices since the lattice potential is typically superimposed by a harmonic trap.

In principle it is possible to create potentials that are nearly flat, such as higher order traps approaching box potentials or new types of traps created by masks or holographic techniques [76, 77, 78, 79]. In these cases, the density of the system is homogeneous, and now the relevant thermodynamic quantity to minimize is the Helmholtz free energy $F$, which can be calculated from the grand potential with

[TABLE]

where in the lattice case $f_{\uparrow}$ and $f_{\downarrow}$ are the numbers of atoms per lattice site, that is the filling fractions or filling factors. Because of the single band model, it is not possible to have more than one atom of one species at a single lattice site, thus the filling factors are bounded between [math] and $1$. When minimizing the free energy $F$, it is important to take care that the chemical potentials satisfy the number equations (67) and (68). Numerically this means that for every $(\Delta,{\bf q})$ -pair, before the value of $F$ can be calculated, the number equations must be solved. Since in the interacting case $\mu_{\uparrow}$ and $\mu_{\downarrow}$ are not independent, and the value of $\mu_{\uparrow}$ affects the value of $\mu_{\downarrow}$ and vice versa, these equations have to be solved iteratively until both chemical potentials have converged.

2.8 Ginzburg-Landau expansion

As explained in section 2.7, the phase diagram of a two-component Fermi mixture can be determined from the thermodynamic potential $\Omega$ (or free energy $F$) by minimizing it with respect to the order parameter. However if one is interested in the phase transitions of the system like the onset of superfluidity, it can be enough to study the change in behavior of the thermodynamic potential, which is explained in this subsection by doing a Ginzburg-Landau expansion of the thermodynamic potential. For further reading, see for instance references [59, 60, 21].

Near a continuous phase transition, the BCS order parameter $\Delta({\bf r})$ is small, and the thermodynamic potential in equation (77) can be expanded as

[TABLE]

where the dots denote the higher orders in $|\Delta|^{2}$ and in the gradients $|\nabla\Delta|^{2}$. The expansion contains only even terms in $\Delta$, because the thermodynamic potential is even in $\Delta$, see equation (77) and figure 4. We only include the smallest possible gradient term that is even in $\Delta$, because here we are only interested in the onset of possible inhomogeneous FFLO phases.

The Landau coefficients in the equation (79) all depend on the system parameters, such as chemical potentials and temperature. If all these Landau coefficients are positive, the minimum of the thermodynamic potential is located at $\Delta({\bf r})$ equal to zero, and the system is in the normal state. A phase transition to a superfluid state has occurred when the position of the global minimum is located at a nonzero order parameter $\Delta({\bf r})$, which describes a condensate of Cooper pairs.

If the coefficient $\gamma$ is positive, it costs energy to have a spatially varying superfluid, such as an FFLO state. A homogeneous phase is then energetically more favorable, and we can therefore restrict ourselves to a pairing field $\Delta$ independent of position. In the opposite case, where $\gamma<0$, the system can gain energy when the order parameter varies in space. The plane wave FF and standing wave LO state are the two simplest examples of the infinitely many possibilities of position-dependent order parameters.

Now in the case where $\gamma$ is positive and the system is homogeneous, both continuous and discontinuous phase transitions from the normal to the superfluid phase can occur. If in the above Landau thermodynamic potential $\alpha$ is negative and all other coefficients are positive, the minimum of the thermodynamic potential is attained at some nonzero $\Delta$ and the Fermi gas is in the superfluid state. Thus as $\alpha$ changes sign, a continuous phase transition takes place, which can be determined by solving $\alpha=0$.

It is also possible to have a first-order phase transition, when in the thermodynamic potential density all coefficients but the fourth-order coefficient $\beta$ are positive. The thermodynamic potential now has two minima, one located at $\Delta$ equal to zero and the other at a nonzero value of the order parameter. For temperatures above $T_{c}$, the global minimum of the thermodynamics potential is located at $\Delta$ equal to zero and the Fermi gas is in the normal state. At the critical temperature, the two minima are equal, and for even lower temperatures, the minimum located at a nonzero order parameter $\Delta$ is the global minimum and the system is in the superfluid state. The criterion for a first-order phase transition thus is $\Omega(0)=\Omega(\Delta)$. In contrast to a second-order phase transition, the location of the global minimum of the thermodynamic potential now changes discontinuously, see figure 4b, which is associated with phase separation. In contrast to the continuous transition, for determining a first-order phase transition, the full thermodynamic potential is needed.

For population imbalanced Fermi gases in optical lattices, both continuous and discontinuous phase transitions occur, as we discuss in section 3. Therefore there must exist a point in the phase diagram where the character of the phase transition changes, that is, a multicritical point, and a sketch of this scenario is given in the left panel in figure 5. If the phase transition changes from being second order to a first order transition, the multicritical point is a tricritical point, which can be found by setting

[TABLE]

Also inhomogeneous phases can occur for the population imbalanced Fermi gas in optical lattices. Therefore there must exist a point in the phase diagram where the phase transition changes from a transition from a normal to a homogeneous superfluid phase into a transition from the normal to an inhomogeneous superfluid phase, which is called a Lifshitz point, see figure 5 right panel. The criterion for the Lifshitz point is

[TABLE]

Superfluidity at nonzero momentum can be established in many ways, of which FF and LO are only two examples, and due to this variety of possibilities it is difficult to predict which kind of inhomogeneous superfluidity will actually occur in the system. However they all have to emerge from the Lifshitz point, making it worthwhile to know about its existence and position. Depending on the strength of the interaction, the Lifshitz point and tricritical point can overlap, $\alpha=\beta=\gamma=0$. In that case, it is necessary to study the full thermodynamic potential to determine what the ground state of the system below this multicritical point actually is.

For the two-component Fermi gas in continuum, many studies exist that investigate the vicinity of the Lifshitz point, for instance references [80, 5, 62, 81]. For Fermi gases in optical lattices, we are not aware of any such studies.

Now that we have some understanding about the different coefficients in equation (79), we can start to calculate them. In mean-field theory, the quadratic Landau coefficient $\alpha$ is explicitly given by

[TABLE]

and can straightforwardly be obtained by taking a derivative of the full thermodynamic potential in equation (77). Note that when $\alpha$ is set to zero, the above equation is exactly the gap equation in equation (70) with ${\bf q}=0$. An alternative way to determine the above coefficient, which is rather simple and elegant, is by using Feynman diagrams [82]. Also the other coefficients $\gamma$ and $\beta$ can be determined using a diagrammatic language. Namely the quadratic coefficient $\alpha$ corresponds to the amplitude of the so-called ladder diagram where the incoming and outgoing bosonic pairing fields $\hat{\Delta}$ carry zero momentum, see figure 6a. The external momentum of the diagram is here taken equal to zero, since the Cooper pairs in the BCS state do not have a net momentum. Physically, $\alpha$ can be interpreted as being proportional to the chemical potential of the Cooper pairs. The fourth-order coefficient $\beta$ has a diagrammatic representation with four external bosonic fields carrying zero momentum, see figure 6b, or can be calculated by taking derivatives from the full thermodynamic potential. It is given by

[TABLE]

In the FFLO phase, the Cooper pairs carry a net momentum. So for the transition to the inhomogeneous superfluid phase, we consider the ladder diagram where the bosonic pairing fields carry nonzero momentum $\mathbf{q}$. The expression for the amplitude of this diagram reads

[TABLE]

For the onset of an FFLO phase, it is enough to study the dependence of this quadratic coefficient on system parameters. To find the location of the Lifshitz point, the first point where FFLO occurs, $\alpha({\bf q})$ is expanded in powers of ${\bf q}$

[TABLE]

which contains only even powers, because $\alpha({\bf q})$ is even in ${\bf q}$, see figure 7, and where the dots denote higher order powers in $\mathbf{q}^{2}$. If all coefficients $a_{i}$ are positive, $\alpha(\mathbf{q})$ has its global minimum at ${\bf q}=0$, and the corresponding phase transition is from the normal to a homogeneous superfluid phase, see the left panel in figure 7. Whereas, if $a_{1}$ is negative, $\alpha(\mathbf{q})$ has its global minimum at a nonzero momentum ${\bf q}$, right panel in figure 7. In the latter case, the phase transition is from a normal to an inhomogeneous superfluid phase. Comparing this with the Landau expansion of the thermodynamic potential in equation (79), we see that in the expansion of $\alpha(\mathbf{q})$, $a_{0}$ can be identified with the quadratic coefficient $\alpha$ and $a_{1}$ can be identified with the gradient coefficient $\gamma$. Thus from the ladder diagram with external momentum, an expression for the $\gamma$ coefficient can be found, namely

[TABLE]

Physically, $\gamma$ can be interpreted as being proportional to the inverse effective mass of the Cooper pairs.

The above expression for $\alpha({\bf q})$ can actually be obtained from the FF thermodynamic potential, by differentiating it with respect to $|\Delta|^{2}$ and consequently taking $\Delta=0$, but not taking ${\bf q}$ to zero, and is equal to the gap equation in equation (70). But importantly (and as seen from the above explanation using Feynman diagrams) the expression for $\alpha({\bf q})$ is an expression for a position dependent order parameter containing a single wave vector ${\bf q}$. Namely in the above calculation, the only assumption made was that the bosonic field carries a net momentum, whereas no assumption was made on the dependence on position of the order parameter.

Also importantly, the Ginzburg-Landau expansion can only be applied close to a continuous phase transition. In practice, one thus always needs the full thermodynamic potential to check if this is the case. The Ginzburg-Landau expansion is therefore expected to apply at the transition from normal Fermi liquid to the FFLO state, but not for the transition from FFLO to BCS which happens for smaller values of imbalance (chemical potential difference).

2.9 The FFLO state and the Bloch’s theorem

The Bloch’s theorem states that there is no net current in the ground state. Since the FFLO state is characterized by a non-zero momentum of the Cooper pairs, it is important to confirm whether it fulfils the Bloch’s theorem. We use the following definition for the net momentum contained in the FFLO state:

[TABLE]

Periodic boundary conditions were assumed. Now one can apply the inverse of the Bogoliubov transformation equation (54) to obtain

[TABLE]

At zero temperature this becomes

[TABLE]

The first term corresponds to paired particles and the momentum $\mathbf{q}$ they carry. The terms $\mathbf{k}$ and $-\mathbf{k}$ cancel in the paired case. Note that the case where both eigenenergies would be negative cannot occur, as seen from the definition (55). The other two terms in equation (88) describe the momenta of the unpaired particles. In general, the three contributions cancel each other to fulfill Bloch’s theorem. The last term corresponds to the unpaired majority particles, thus one can immediately see that an area of unpaired particles should form to momenta opposite the direction of the FFLO wave vector $\mathbf{q}$ to cancel the first term. The case where only majority particles are unpaired corresponds to the left panel of figure 3 where the dispersions cross zero (gapless excitations, i.e., free particles) only for one sign of momenta. But for some parameter regimes, there can be also free (gapless) minority particles, corresponding to the second term in the above equation (88). These are visible in the eigenvalues of figure 3, the right panel, as crossings of zero at both positive and negative momenta. The formation of unpaired regions is illustrated schematically in figure 8, and via numerical data in figure 11.

One can easily calculate the net momentum equation (2.9) or equation (88) from a numerical solution of the problem and check whether it is zero. The zero temperature case can be approached analytically. Bloch’s theorem was central in the original work of Fulde and Ferrell [2]: it was used to select the ground state from a larger class of possible solutions. A further analytical treatment was done in early classic work by Takada and Izuyama [83]. They derived analytical formulas for the areas in momentum space where the unpaired particles exist, and they expressed the system energy in terms of these formulas. In the calculation, an approximation was made by neglecting terms of the order $q^{2}/k_{F}$, that is, it was assumed that the FFLO wave vector’s absolute value is small compared to that of the Fermi wave vector. The authors then showed that the net current defined analogously to equation 86 becomes proportional, up to constants, to the derivative of the total energy with respect to the FFLO wave vector:

[TABLE]

Bloch’s theorem is therefore equivalent to the demand that the FFLO solution is an extremum of energy. Fulfilling Bloch’s theorem is thus a necessary condition for the existence of an FFLO solution that can minimize the system energy. However it is not sufficient because one further needs to show that the extremum is actually a minimum. The connection between the current being zero and minimizing the energy is also shown in detail in [72].

A similar analysis can be performed in the case where the order parameter ansatz is more complicated than the FF single $\mathbf{q}$ case. When the order parameter ansatz contains two opposing momenta, as with the LO ansatz, it is quite intuitive to assume that there is no net current. In general, one should always confirm whether the solution for a more complex FFLO state fulfills Bloch’s theorem.

The FF state breaks time-reversal symmetry (TRS) because it is characterized by a single wavevector $\mathbf{q}$ with a definite direction. It might seem counter intuitive that nevertheless net current is absent. However note that there actually exists spin current in the FF state: the pairs with wavevector $\mathbf{q}$ have the same amount of up and down particles while the unpaired particles that carry momenta that cancel $\mathbf{q}$ are mostly or only majority particles, leading to net spin current.

2.10 Luttinger’s theorem and Fermi surfaces

Luttinger’s theorem states that a volume enclosed by a Fermi surface is directly proportional to the particle density. This is easy to see in the case of non-interacting, zero temperature Fermi gas. In case of the BCS state, the Fermi surface disappears. This is consistent with Luttinger’s theorem because in the BCS state, the particle number (density) fluctuates, only its mean value is fixed. This is connected to the breaking of the $U(1)$ symmetry and appearance of a condensate phase. Now in FFLO and Sarma/BP states, there are gapless excitations, that is, Fermi surface(s). One may ask, does this somehow contradict Luttinger’s theorem, since these are superfluids and the particle number fluctuates? But remember that it is the number of pairs that fluctuates in the superfluid: actually the difference $n_{\uparrow}-n_{\downarrow}$ remains constant. The volumes enclosed by Fermi surfaces in the FFLO and Sarma states are related to particle density difference.

Sachdev and Yang [84] considered a model where two fermion species can form a boson (which can be understood as a molecule or a Cooper pair) and rigorously showed that the particle number difference of the two species is given by

[TABLE]

where $\Omega_{\uparrow}$ and $\Omega_{\downarrow}$ are the volumes of the two Fermi surfaces and $d$ the dimension. This holds also for vanishing Fermi surface volumes, that is, there can be two, one or zero Fermi surfaces depending on the phase of the system. In case of a periodic order parameter, such as in FFLO state, the Fermi surface volumes are well-defined modulo the Brillouin zone volume. Quantum phase transitions are associated with the change in the number of Fermi surfaces when the system parameters are varied. In a recent article [85], the validity of the Luttinger’s theorem was proven for any approximate theory based on a conserving approximation in case of spin-density imbalance.

2.11 Connection to the repulsive Hubbard model: the FFLO state and striped phases vs. d-wave superconductors

In bipartite lattices, it is possible to make simple particle-hole transformations that connect phenomena that occur for attractive and repulsive interactions [86]. This provides an interesting connection between the FFLO state in lattices and important open questions in high-$T_{c}$ superconductivity. From experiments, the order parameter in many high-$T_{c}$ materials is known to have d-wave symmetry. The repulsive Hubbard model is viewed as a minimal model that could perhaps describe high-temperature superconductors, such as cuprates. Analytical solutions of the Hubbard model are not available, and numerical methods like quantum Monte Carlo are plagued by the numerical sign problem in the regime where superconducting phases would occur, namely away from half filling. Thus the true ground state over the whole parameter space is not known.

At half filling, where quantum Monte Carlo can give reliable results, numerical simulations predict anti-ferromagnetic (AFM) order [87, 88]. Away from half filling, several competing phases have been predicted. For instance, approximative numerical approaches predict a d-wave superfluid where the pairing involves correlations between particles on several neighbouring sites. Such pairing cannot be described by a mean-field theory within the basic Hubbard model that only has on-site interactions. A d-wave superfluid is of course exciting, considering the observed d-wave symmetry of the order parameter in many real materials. The d-wave superfluid may compete or co-exist with the AFM order in the Hubbard model [87, 89], but not only with that. There are suggestions that also a so-called striped phase could exist [90]. In the striped phase, the amount of doping varies spatially, on top of AFM order, see figure 9. Such a phase can be described within mean-field theory. It is not known which of these phases are true ground states at each area of the parameter space, and whether there are regimes where some or all of them co-exist. Co-existence of AFM and superconducting order has been studied, see for instance [91, 92] and references therein, as well as the interplay between stripes and d-wave superconductivity ([93, 94] and references therein). Recently cluster DMFT calculations that included the possibility of AFM order, d-wave superconductivity and the striped phase simultaneously were performed. The results show coexistence of the d-wave superconductor and striped phases for a large doping region and a smaller area of a uniform d-wave superfluid [95]. Existence of the striped phase for one value of doping was also recently shown by a study combining several different numerical approaches [96], and for a Lieb-lattice geometry [97].

As we explain below, the striped phase is connected to the FFLO state via a particle-hole mapping [98, 99, 100, 101]. Therefore, studies of the FFLO state in optical lattices will provide new knowledge on the coexistence of superfluidity and magnetization, and profoundly important information about the Hubbard model in general and its relevance in explaining high-$T_{c}$ superconductivity.

Let us consider the Hubbard model in equation (9), but with the terms $\sum_{\sigma{\bf r}}\left(\frac{\mu^{0}_{\sigma}}{2}+\frac{U}{4}\right)$ and $-\sum_{\sigma{\bf r}}\frac{U}{2}\hat{c}_{\sigma{\bf r}}^{\dagger}\hat{c}_{\sigma{\bf r}}$ added to the Hamiltonian. The first one is a constant, and the second means a trivial shift of the chemical potentials. This produces the so-called particle-hole symmetric Hubbard Hamiltonian

[TABLE]

This Hamiltonian possesses certain symmetries with respect to the repulsive and attractive interactions in bipartite lattices. We first consider the case of having the same chemical potential and particle numbers for $\uparrow$ and $\downarrow$ particles. There exists a unitary transformation $U_{1}$ (the label one is just to make a difference to the interaction energy $U$ of the Hubbard model) that connects the two at half filling, namely

[TABLE]

where $\epsilon(i)=1$ for one sublattice of the bipartite lattice and $\epsilon(i)=-1$ for the other. With this transformation

[TABLE]

which means that the Hamiltonians for both signs of $U$ must have the same energy spectrum and the same form of the ground state. Therefore at half filling, if one knows the ground state for $U$, the ground state for $-U$ can be straightforwardly obtained. It turns out that the usual BCS superfluid for attractive interactions maps to a Néél-ordered anti-ferromagnet on the repulsive side. If strong attractive interactions are considered, on-site pairs are formed and they Bose-condense; there exists a BCS-BEC crossover in the lattice. The BEC of on-site pairs corresponds to a Mott antiferromagnet at strong interactions. The existence of on-site pairs at higher temperatures, before one reaches the superfluid below the critical temperature, corresponds to having a Mott state before reaching a Mott antiferromagnet. The identical behaviour at negative and positive interaction $U$ has been seen, for instance, in expansion experiments of ultracold Fermi gases [102, 103]. Importantly in an experimental situation where the symmetry of the system may be broken by some unwanted or desired features, the mapping may not be exact.

Away from half filling and for differing chemical potentials and particle densities, there is still a mapping, but it is more complicated. Now it turns out that the negative and positive $U$ Hubbard models map to each other in such a way that spin imbalance, that is, different chemical potentials $\mu_{\uparrow}\neq\mu_{\downarrow}$ in one model maps to chemical potential away from half filling in the model with opposite sign of $U$. To be precise, with the notation $\mu\equiv(\mu_{\uparrow}+\mu_{\downarrow})/2$, $\Delta\mu\equiv(\mu_{\uparrow}-\mu_{\downarrow})/2$ the mapping is the following:

[TABLE]

Since the total density $\frac{1}{M}\sum_{{\bf r}}\langle\hat{c}_{\uparrow{\bf r}}^{\dagger}\hat{c}_{\uparrow{\bf r}}\rangle+\langle\hat{c}_{\downarrow{\bf r}}^{\dagger}\hat{c}_{\downarrow{\bf r}}\rangle=n$ ($M$ is the number of the lattice sites) and density difference $\frac{1}{M}\sum_{{\bf r}}\langle\hat{c}_{\uparrow{\bf r}}^{\dagger}\hat{c}_{\uparrow{\bf r}}\rangle-\langle\hat{c}_{\downarrow{\bf r}}^{\dagger}\hat{c}_{\downarrow{\bf r}}\rangle=\Delta n$ transform as $\frac{1}{M}\sum_{{\bf r}}\langle 1-\hat{d}_{\uparrow{\bf r}}^{\dagger}\hat{d}_{\uparrow{\bf r}}\rangle\pm\langle\hat{d}_{\downarrow{\bf r}}^{\dagger}\hat{d}_{\downarrow{\bf r}}\rangle$ one straightforwardly obtains relations for the total densities in the mapping:

[TABLE]

Likewise the polarization $P=(n_{\uparrow}-n_{\downarrow})/(n_{\uparrow}+n_{\downarrow})=\Delta n/n$ maps as $P\longleftrightarrow(1-\Delta n)/(1-n)$. At half filling, $n_{\uparrow}+n_{\downarrow}=n=1$, $\mu=0$ and the relations can be further simplified. For instance assuming an attractive interaction system with finite polarization $P$, the situation where the FFLO state may occur, then at half filling for the attractive system (($\mu_{\uparrow\mathrm{attr}}+\mu_{\downarrow\mathrm{attr}})/2=0$, $n_{\mathrm{attr}}=1$) we have $\Delta\mu_{\mathrm{attr}}=-\mu_{\mathrm{rep}}$ and $1-P_{\mathrm{attr}}=n_{\mathrm{rep}}$. The essence of the particle-hole transformation is essentially unchanged when using the standard Hubbard Hamiltonian of the form in equation (9), but then one has to take into account extra chemical potential shifts and constants appearing in the transformation when comparing actual numbers. The particle-hole symmetric Hamiltonian (2.11) in contrast leads to the simple relations between chemical potentials and total densities as given in equations (96)–(99).

This means that the FFLO state in the attractive Hubbard model corresponds to the striped phase that has been predicted for the doped repulsive Hubbard model [61]. Therefore, observing the FFLO state would address open questions, discussed above, related to the repulsive Hubbard model as well, with potentially great relevance to high temperature superconductivity.

In ultracold quantum gas (UQG) experiments, however, one has to be cautious with the connection between the repulsive and attractive Hubbard models because the lattice potential is typically superimposed by another trapping potential, for instance a harmonic one. This has important consequencies to the particle-hole mapping and to the possible equivalence of phases predicted for the attractive and repulsive models. This is thoroughly analysed by Ho, Cazalilla and Giamarchi [104]. The basic ingredients of the Hubbard Hamiltonian are the hopping and the interaction terms as seen from equation (2.11); let us write here a generic form to the rest of the Hamiltonian and call it $\hat{H}_{\mathrm{ext}}$ (external potential), also including the chemical potentials:

[TABLE]

where $V_{i}$ is the energy shift caused by a spatially varying external trapping potential and $h_{i}$ is a Zeeman field. If the particle number polarization is determined in the preparation of the gas and is uniform, as is usually the case in quantum gas experiments, $h_{i}=h$ does not depend on position and thus does not disturb the particle-hole mapping. But the term $\sum_{i}(V_{i}-\mu)$ does depend on position and has important consequencies. In the repulsive case, it might be difficult to obtain the desired doping uniformly since it is energetically favorable to have the center of the trap fully filled and the edges with lower density. On the attractive side with population imbalance, phase separation into a BCS-paired trap centre and polarized edges might be favored analogously to trapped continuum gases (see section 8). All this depends on parameters though, and one may hope to find parameter regimes where the trapping effects are negligible, or that the trap potential changes so slowly that particular phases can be found in different trap areas (c.f. discussion in section 7.2). Also the newly achieved box potentials [76, 78, 105, 79] offer interesting solutions for uniform gases and thus also straightforward applicability of the particle-hole mapping.

3 Mean-field studies on FFLO in 2D and 3D optical lattices

The mean-field theory of the FFLO state in two-component Fermi mixtures in optical lattices was considered by Koponen et al. [69, 106]. The authors found that both the Sarma and the FF state can form the ground state of the system. The Sarma state was found to be stable for fixed densities while the FF state is the energy minimum for fixed chemical potentials; this result is similar to the continuum case where fixed densities can stabilize the Sarma state [13, 107]. The authors of [69] show that the parameter regime for which the system is in the FF state is relatively large compared to the continuum due to the nesting effect. Namely in the lattice, the overlap of the Fermi levels for the two components can become substantial by displacing the minority Fermi sea with the FF wavevector. This effect is illustrated in figure 10. Furthermore the authors present the momentum distributions for the FF superfluid phase, see figure 11 where this displacement can be seen. It is shown that the system has no total momentum because the momentum distributions of the individual components cancel the effect of the FF wavevector ${\bf q}$.

In the papers following this first study by the same authors [67, 68], the complete phase diagrams were determined in one-, two- and three-dimensional lattices. In a Letter [67] from 2007, the three-dimensional case was considered and the phase diagram as a function of temperature and polarization at fixed filling was shown to contain an FF region. Koponen et al. [67] also computed the phase diagram as a function of polarization and filling at zero temperature and showed it has has a substantial FFLO region, see figure 12. Importantly, they found that even after considering phase separated states (phase separated normal phase and a BCS phase), FFLO states cover a large area of the phase diagram. The broad range of parameters for which the FFLO state was the ground state is in striking contrast to three dimensional homogeneous systems, where the FFLO exists only in an exceedingly narrow region of the parameter space [2, 4, 108, 109]. In reference [68], Koponen et al. extended a mean field exploration of the FFLO state in optical lattices to different dimensions and discussed the role of the Hartree corrections. In all dimensions, nesting and the Van Hove singularity plays an important role, explaining the large regions of the parameter space where FFLO is the ground state. Momentum distributions and density-density correlations were also computed.

By using a more elaborate variational mean-field approach based on the Bogoliubov-de Gennes (BdG) method and including Hartree corrections, Loh et al. [110] studied the stability and observable signatures of the LO phase in a cubic lattice. In their approach, six variational parameters were self-consistently computed at each site: complex valued order parameter, chemical potential, and three Zeeman fields. In a cubic lattice, the LO phase is stable in a much broader range of parameters than expected in the continuum. Furthermore, they computed that the LO phase could exist in a lattice if entropy per particle is below $\sim 0.5\,k_{B}$. If the lattice was anisotropic, the stability of the LO phase would be enhanced further. While the pairing energy is small compared to the overall energy, compared to the FF phases that break the time-reversal symmetry, the LO phases could have a pairing energy that exceeds the pairing energy of the FF phases even by a factor of $50$.

A mean-field Hartree-Fock-Bogoliubov theory was applied by Chiesa and Zhang as well as by Rosenberg et al. [111, 112] to study a polarized Fermi gas in a three-dimensional square lattice. Without targeting states of specific form and via extensive numerical work, they found that the LO state with a single ${\bf q}$ is energetically favored and intriguingly that the favored direction of the LO state ordering can change as the interaction strength is varied. In particular, a ”diagonal” order with ${\bf q}\propto(1,1,1)$ can become favored over ${\bf q}\propto(0,0,1)$- order as interactions become stronger and the system becomes more strongly filled. An example of the relevant phase diagram is shown in figure 13 for a relatively large interaction strength. Diagonal order is indicated by green circles and appears for high densities and relatively small polarizations.

In [66], an explicit comparison was made between the FF and the LO ansatz within mean-field theory in a two-dimensional square lattice and, in agreement with earlier studies [72, 110, 113], the LO state is energetically more favorable than the FF ansatz. Moreover in contrast to results in a three-dimensional system [112], it was concluded that it is energetically more favorable to have the LO wave vector ${\bf q}$ along one of the lattice vectors, as compared to a LO ${\bf q}$ that makes a 45o angle with the lattice vectors. In the LO state the translational symmetry is broken, which is reflected in the particle densities. An important comparison between between the mean-field phase diagram and a phase diagram obtained by quantum Monte Carlo calculations [114] was also made in [66]. Interestingly the transition lines agree very well if the mean-field temperatures and polarizations are simply scaled by a constant, see figure 14.

A somewhat different problem, intermediate between continuum and lattice results, was studied by Devreese et al. [115, 74]. They considered a three dimensional system where a one-dimensional optical lattice was imposed along the $z$-direction. By studying the mean-field theory as a function of interaction strength and lattice parameters, they found that the FFLO state can persist to higher spin imbalance when the lattice period was reduced.

Reduced dimensionality is known to enhance possibilities for FFLO states. Using a mean-field theory, with a possibility of also pair tunneling, Sun and Bolech [116] mapped the behavior of a system with coupled one-dimensional tubes. When barriers separating tubes were high, the FFLO state appeared in the center of the tubes. As barriers were lowered, tunneling rate increased, and the central FFLO state gave way to a central BCS-like state. Mean field studies in a 1D lattice found in addition to the FFLO state the possibility of an $\eta$ phase, that is, a phase where the FFLO wave vector is at the first Brillouin zone edge [117].

4 Mass-imbalanced Fermi gases in lattices

It is also possible to have mixtures of fermions with different masses, and such mixtures have indeed been experimentally realized in ultracold atom experiments [118, 119, 120, 121, 122]. Mass imbalance causes some important physical changes. First, having equal chemical potentials does not mean equal densities in the mass-imbalanced limit. Second, the dispersions for the two components are different, meaning that already in the population balanced case there is a Fermi level mismatch.

Qualitatively it is easy to see that mass imbalance and chemical potential imbalance can give rise to similar effects. For an ideal gas at zero temperature in a continuum, we have

[TABLE]

and consequently $k_{F,\sigma}=\sqrt{2m\mu_{\sigma}}/\hbar$. Now if the chemical potentials differ but the masses are the same, Fermi surfaces are separated by

[TABLE]

On the other hand, if the masses differ but the chemical potentials are the same

[TABLE]

Mismatch in Fermi surfaces can be achieved either by changing the chemical potentials or by changing the particle masses. Similar arguments also apply in the lattice. For example for tight binding dispersion in a one-dimensional lattice, the connection between chemical potential and $k_{F,\sigma}$ is given by $2J_{\sigma}\cos(k_{F,\sigma})=\mu_{\sigma}$ and inverse of the hopping strength $1/J_{\sigma}$ plays an equivalent role as a particle mass. However the analogue in Fermi momentum mismatch does not necessarily mean that the pairing physics is the same in both cases.

BCS quasiparticle dispersions in a continuum system with unequal masses [123] are

[TABLE]

where $m_{\pm}=1\pm m_{\uparrow}/m_{\downarrow}$, $\tilde{\mu}_{\pm}=(\mu_{\uparrow}\pm\mu_{\downarrow})/2$, and $\epsilon_{k}=\hbar^{2}k^{2}/(2m_{\uparrow})$. If masses are the same $m=m_{\downarrow}=m_{\uparrow}$ while chemical potentials are not, we have

[TABLE]

On the other hand, we could consider the case with unequal masses, but with equal chemical potentials $\mu_{\uparrow}=\mu_{\downarrow}=\mu$. In the weak coupling limit, the wavevectors for which $\epsilon_{k}m_{+}/2\approx\mu$ are the most relevant. Then quasiparticle dispersions can be approximated as

[TABLE]

where

[TABLE]

and $M_{r}=(1/m_{\uparrow}+1/m_{\downarrow})^{-1}$ is the reduced mass. We then see that equations (104) and (105) are the same as long as we make the mappings

[TABLE]

[TABLE]

and

[TABLE]

This demonstrates that in this limit the phenomena expected in mass imbalanced mixtures are similar to those in mixtures with unequal chemical potentials. However this argument is only valid in the weak coupling limit (and in the low filling limit in lattices). In the strong coupling regime, the physics in these two cases is expected to be different. Note that the above argument did not consider the possibility of a FFLO state, and more complex physics results when both mass and chemical potential imbalance are present.

To physically realize an unequal mass mixture in a lattice, one could put unequal mass fermions in a lattice, but effective mass imbalance can also be introduced via different hopping parameters $J_{\sigma,\alpha}$ for the two different spin components. This is because the hopping parameters are inversely proportional to the (effective) masses of the particles, namely

[TABLE]

As an example, a very promising prospect in this regard is the recent experiment presented in [58], where a versatile method for creating state-dependent optical lattices by applying a magnetic field gradient modulated in time was demonstrated. Thus different hopping parameters can be engineered for different internal spin states.

Mass imbalanced mixtures in optical lattices in the context of FFLO states were considered by Wang et al. [124]. Using a time-evolving block decimation (TEBD) technique they find that with increasing mass imbalance, the FFLO region in the phase diagram decreases in size and eventually disappears for a large enough mass ratio. The one-dimensional system was also considered by Dalmonte et al. [125] by applying density matrix renormalization group and quantum Monte Carlo techniques to study population and mass imbalanced Fermi mixtures in a trap. They found that mass imbalance helps to stabilize paired phases against spin-imbalance and concluded that unequal mass mixtures are promising systems to observe FFLO superfluidity in an experimentally realistic system.

Pahl and Koinov [126] considered 6Li-40K mixtures in square optical lattices using a mean-field theory. They found a normal, Sarma and FF phase to be present in the polarization versus a temperature phase diagram. The FF region occurs for non-zero polarization, and the region is larger for a majority of heavier 40K atoms. While the literature on FFLO states in mass imbalanced systems in lattices is not very broad, based on the results in continuum [127], there are reasons to believe that mass imbalance might enhance FFLO type pairing also in a lattice, motivating further studies on this topic.

5 FFLO state in the Green’s function formalism

Here we briefly present how superfluids, such as the FFLO state, can be described using Green’s functions. This is convenient for many purposes, such as obtaining the collective modes and determining the ground state by beyond-mean-field approaches. We first define the Green’s functions in the Matsubara formalism, which allows consideration of finite temperatures. The possibility of pairing is taken into account by using the Nambu formalism where the Green’s function becomes a matrix with the off-diagonal elements accounting to pairing correlations (sometimes called anomalous Green’s functions), while the diagonals are the usual Green’s functions. Forms of the Green’s functions in the Hartree-Fock-Gor’kov mean-field approximations are given and are shown to lead to eigenenergies and coherence factors of the same form as obtained in section 2. After introducing the formalism, we proceed in section 6 to describe interactions between quasiparticles and thereby obtain collective modes of the system. In section 7 beyond mean-field studies of the FFLO state are discussed.

The single particle Green’s function is defined as follows in the finite temperature Matsubara Green’s function formalism [59, 128, 129]

[TABLE]

Here, we use the notation 1 for the variables $\boldsymbol{r}_{1}\tau_{1}\sigma_{1}$, denoting the position coordinate, the imaginary time of the Matsubara formalism, and the (pseudo-)spin, respectively, while $T_{\tau}$ is the time ordering operator and the brackets denote the thermodynamic average defined for the operator ${\hat{O}}$ as

[TABLE]

where $\beta=\frac{1}{k_{\mathrm{B}}T}$ is the inverse thermal energy with $k_{\mathrm{B}}$ the Boltzmann constant and $T$ the temperature. We describe the system in Nambu formalism with an extended spin index $\sigma\in\{1,2,3,4\}$ where ${\hat{c}}_{1}={\hat{c}}_{\uparrow}$, ${\hat{c}}_{2}={\hat{c}}_{\downarrow}$, ${\hat{c}}_{3}={\hat{c}}^{\dagger}_{1}={\hat{c}}^{\dagger}_{\uparrow}$ and ${\hat{c}}_{4}={\hat{c}}^{\dagger}_{2}={\hat{c}}^{\dagger}_{\downarrow}$. Compared to equation (108), we moved the spin index from the argument of ${\hat{c}}$ to its subscript. Thus the anomalous Green’s functions with two creation or two annihilation operators are also defined by equation (108). These are essential for describing pairing correlations on the mean field level. The total Green’s function can be also written in matrix form by taking the spin indices to correspond to the rows and columns of a matrix, as done explicitly in equation (117) and the discussion following it.

The single particle Green’s function has the following equation of motion

[TABLE]

The integral sign is a shorthand notation for summation over position and spin as well as integration over time, while variables of summation and integration are indicated by the overbar. Here $G_{0}$ is the non-interacting Green’s function, which is usually easy to determine exactly. The self energy $\Sigma$, for a general two-body interaction, is defined via

[TABLE]

where $V(1,2)$ is the inter-particle interaction potential between two particles and can generally depend on the particle spatial coordinates and spin as well as on time. Here the notations $\tau^{+}$ and $\tau^{-}$ specify the time ordering $\tau^{+}>\tau>\tau^{-}$ and imply taking the limit $\tau^{\pm}\rightarrow\tau$, while the two particle Green’s function $G_{2}$ is defined as

[TABLE]

Here we consider on-site interactions described by the potential $V(1,2)=U_{\sigma_{1}\sigma_{2}}\delta(\boldsymbol{r}_{1}-\boldsymbol{r}_{2})\delta(\tau_{1}-\tau_{2})$, with $U_{34}=U_{12}=U$ for the extended index. In the Hartree-Fock-Gor’kov mean-field approximation, the self energy – now written explicitly in matrix form – is

[TABLE]

where $\delta_{11^{\prime}}=\delta(\boldsymbol{r}_{1}-\boldsymbol{r}_{1}^{\prime})\delta(\tau_{1}-\tau_{1}^{\prime})$. Here the single particle Green’s functions appearing in $\Sigma$ have the variables $(\boldsymbol{r}_{1}\tau_{1},\boldsymbol{r}_{1}\tau_{1}^{+})$, and the spin indices $\sigma$, $\sigma^{\prime}$ now appear in these Green’s functions as $G_{\sigma\sigma^{\prime}}(\boldsymbol{r}_{1}\tau_{1},\boldsymbol{r}_{1}\tau_{1}^{+})=-\left\langle T_{\tau}\left({\hat{c}}_{\sigma}(\boldsymbol{r}_{1}\tau_{1}){\hat{c}}_{\sigma^{\prime}}^{\dagger}(\boldsymbol{r}_{1}\tau_{1}^{+})\right)\right\rangle$.

The Hartree terms can be included in the chemical potentials on the diagonal of the self-energy matrix in a way similar to equation (12). The Fock-exchange terms $G_{12}$, $G_{21}$, $G_{34}$ and $G_{43}$ are zero in a Fermi gas where spin-flips do not occur. Finally we introduce the key element of the mean field FFLO theory. We assume that the pairing fields of the self-energy have an oscillating structure

[TABLE]

where the FFLO pairing vector is denoted by $\boldsymbol{q}$. The choice $\boldsymbol{q}$ instead of $2\boldsymbol{q}$, which is sometimes used, corresponds to the ansatz for the FF pairing field as in equation (14).

Although the choice of the spin indexing $\sigma\in\{1,2,3,4\}$ leads to four by four matrices, as used in equations (117)–(122), there is redundancy in the labeling because indices $3$ and $4$ simply correspond to Hermitian conjugates of $1$ and $2$. Therefore the essential physics of the system is described by a two by two matrix with usual single particle Green’s functions for the two spins in the diagonal and the anomalous Green’s functions in the off-diagonals; this is the standard form of the Nambu Green’s function. After taking a proper Fourier transform and inversion of equation (108) [130], one can write the Nambu Green’s function as 333Note that in [130, 131], the convention $\Delta\exp(2i{\bf q\cdot x})$ was used instead of $\Delta\exp(i{\bf q\cdot x})$ chosen in this review, another difference in notation is that $v(\boldsymbol{p})$ here corresponds to $-v(\boldsymbol{p})$ in [130].

[TABLE]

In the spin-block $\sigma_{1},\sigma_{2}\in\{2,3\}$, a similar calculation yields the expression

[TABLE]

Here appear the quasiparticle energies $E_{\pm}$

[TABLE]

where $\xi^{\pm}_{\boldsymbol{p},\boldsymbol{q}}=\left(\xi_{1}(\boldsymbol{p})\pm\xi_{2}(\boldsymbol{q}-\boldsymbol{p})\right)/2$, the free particle energies correspond to those defined in equation (22) as $\xi_{1}=\xi_{\uparrow}$ and $\xi_{2}=\xi_{\downarrow}$, and the coherence factors $u$ and $v$

[TABLE]

The quasiparticle energies and the coherence factors have the same form as when derived within the simple BCS formalism, for comparison of the quasiparticle energies see equations (55) and (57) in section 2.

6 Collective modes for the FFLO state

Superfluidity can be explained within the Landau symmetry breaking paradigm. The $U(1)$ symmetry that allows freedom to choose the overall phase factor of a quantum state becomes broken: the superfluid wave function obtains a well-defined phase [59, 60]. In the case of a superfluid of Cooper-paired fermions, this is the phase $\varphi$ of the order parameter $e^{i\varphi}\Delta$, usually chosen to be zero for convenience. As an important consequence, a low-energy collective mode is formed related to this symmetry breaking: it is possible to distort the phase $\varphi$ by a small gradient with low energy cost. Physically this corresponds to a collective movement of all condensed particles, for instance the Cooper pairs. This is an example of the more general Goldstone theorem which states that spontaneous breaking of a continuous symmetry leads to the existence of a Nambu-Goldstone mode, which is a bosonic excitation of the ground state and has a gapless dispersion at the long wavelength limit.

In the case of superfluids, the Nambu-Goldstone mode is called the Anderson-Bogoliubov phonon and is one of the defining characteristics of superfluidity. However it exists as a low-energy excitation only for neutral superfluids, such as those in ultracold quantum gases (UQG). For charged systems like superconductors, the situation is different. There the order parameter phase (i.e., the movement of the charged particles) becomes coupled to the vector potential. As a result, the collective mode becomes gapped, that is, massive; this is an example of the Higgs mechanism. The Anderson-Bogoliubov phonon naturally affects the system’s thermodynamics. BCS theory in its simplest form does not describe such collective excitations, nevertheless it is able to predict some basic quantities like the critical temperature. However one has to be cautious of phases of matter predicted by minimizing the energy within a simple mean-field theory because collective modes, related to order parameter and density fluctuations, may destabilize the state. This can be especially critical for the FFLO case where fluctuations may become strong [109, 81]; in their original article [2], Fulde and Ferrell already showed that the FF state does not have supercurrent for small values of the vector potential in the direction transverse to the FF wavevector, hinting to instability to fluctuations. Extending mean-field theories, collective excitations can be captured by introducing interactions between quasiparticles, as described below.

The existence and properties of low-energy collective modes becomes an intriguing question in the case of the FFLO state because not only the $U(1)$ symmetry is broken as the superfluid order parameter forms, but also the translational and rotational symmetries may be broken due to the spatial dependence of the order parameter. The breaking of rotational symmetry is apparent from the FF form of the order parameter $\Delta e^{i\mathbf{q}\cdot\mathbf{r}}$, but since the spatial dependence there is of plane wave type, one needs an LO-type state (e.g., $\Delta cos(\mathbf{q}\cdot\mathbf{r})$) to have the translational symmetry broken as well. In LO-type states also the densities and order parameter amplitudes depend on position. One can ask whether and what kind of collective modes appear as a consequence of such symmetry breaking. Moreover the FFLO state hosts unpaired particles in addition to Cooper pairs, that is, a Fermi surface exists, which could affect the collective modes.

The collective mode spectrum of the FFLO state has been considered in the context of superconductivity in quantum chromodynamics (see section IV of the review [5]). Effective Lagrangians for the Goldstone bosons were derived, and it was shown that indeed Nambu-Goldstone modes exist. The system is quite complex since quarks have both color and flavor degrees of freedom, in addition to different masses, and the condensates can be of vector form. When two-component, attractively interacting ultracold Fermi gases emerged as a potential (charge neutral) system to experimentally realize the FFLO state, it became relevant to reconsider the question of FFLO states and collective modes in such systems. The purpose of this section is to review theory work on FFLO state collective modes in the ultracold gases context. We first briefly present one method that can be used for calculating the collective modes, namely linear response with the generalized random phase approximation (GRPA). We proceed to discuss some selected publications where the question of collective modes of the FFLO state in UQG has been approached.

Calculation of the collective mode spectrum can be done within the linear response theory, assuming thermodynamic equilibrium. The change in the expected value of a relevant observable ${\hat{O}}$ in response to a time-dependent perturbation described by a interaction-picture operator ${\hat{H}}^{\prime}(t)$ is given by the Kubo formula [59, 128, 129]

[TABLE]

The perturbation ${\hat{H}}^{\prime}(t)$ is assumed vanishing for $t<0$, and $\Theta(t)$ is the Heaviside step function. Here both ${{\hat{H}}}^{\prime}(t)$ and ${\hat{O}}(t)$ are single-particle operators of the generic form

[TABLE]

where $\phi$ now characterizes the strength and spatial, temporal and spin-dependence of the perturbing field. One should understand $\tilde{\phi}$ as a generic way of defining the observable of interest: it is usually independent of time, selects the spin components of interest, and defines whether the observable is dependent on a spatial coordinate or whether all spatial coordinates are integrated over.

In the case of single particle operators as in equation (145), the integrand of the Kubo formula, that is the integrand of $\delta\left\langle{\hat{O}}(t)\right\rangle$ is a retarded two-body correlation function, or linear response function, multiplied by the strength of the perturbation $\phi$. The definition of the retarded linear response function is

[TABLE]

One can now obtain collective modes as a response to, for instance a density perturbation, from the linear response function. However its exact evaluation is in practice an impossible task for an interacting many-body system in dimensions higher than one, for realistic system sizes. Approximations are needed, but with too crude approximations, the collective modes can be missed. Basically one has to include the existence of quasiparticles as well as interactions between them to arrive at the relevant collective modes. Random phase approximation (RPA) [132] was introduced originally to describe collective effects in an electron gas and was later generalized (GRPA) to describe collective modes in superconductors [133]. The GRPA form of the linear response function is usually sufficient to find the low-energy collective modes of superfluids.

The Kadanoff-Baym formalism [134, 135] is a method to rigorously derive the linear response function (for example the GRPA response) so that it is conserving and thermodynamically consistent regarding the approximations made in the many-body problem. A conserving approximation is such that it obeys the conservation laws arising from the symmetries of the physical system. Thermodynamical consistency means that the partition function is unique: direct evaluation of observables and their calculation by differentiating the partition function give the same result. The GRPA linear response function derived from the Kadanoff-Baym formalism is (compared to equation (146), the spin indices are now moved from the subscript to the arguments)

[TABLE]

where the subindex [math] signifies the evaluation in absence of perturbation, $\phi=0$. In mean-field approximations of the Hartree-Fock-Bogoliubov type, the functional derivative $\left(\frac{\delta\Sigma(\bar{3},\bar{4})}{\delta G({5},{6})}\right)_{\phi=0}$ is a combination of delta functions, as one can conclude from equation (117), and does not produce any additional Green’s functions to the equation.

Edge and Cooper studied a 1D two-component Fermi gas with attractive contact interactions using a mean-field BdG approach (c.f. section 8 of this review) and the GRPA formalism [136, 137]. The BdG approach can describe a spatially oscillating order parameter, for instance of cosine-type, going beyond the simple FF ansatz $\Delta e^{i\mathbf{q}\cdot\mathbf{r}}$. In this case, both translational and rotational symmetries and the densities and absolute value of the order parameter amplitude are spatially varying. Edge and Cooper studied the response of the system to a density perturbation periodic in time and space, thus giving the momentum $k$ and frequency $\omega$ to the system. They found two collective mode branches in the density response at low frequency and wavelength. One of them appears at $k=0$, and the other one at $k^{*}=2(k_{F\uparrow}-k_{F\downarrow})$, see figure 15. Here the factor of two compared with the FFLO wave vector $q=k_{F\uparrow}-k_{F\downarrow}$ appears because the densities oscillate in space with a period only half of the cosine-form order parameter oscillation, because they follow the absolute value. The appearance of two gapless modes is due to the FFLO state breaking both the $U(1)$ and translational symmetries. The period of the density and order parameter absolute value oscillations give $k^{*}$. The two modes involve spatial oscillations of both density and spin-density. This is a consequence of the breaking of time reversal symmetry by the imbalance. A detailed analysis of how the equivalents of such modes appear in a trapped gas is presented in references [136, 137].

While the mean-field analysis of refs. [136, 137] can provide qualitative insight to the behaviour of collective modes, it cannot be an accurate description of the 1D system. In 1D no continuous symmetries are spontaneously broken, and only power-law decaying correlations exist. However the qualitative results match well with the knowledge obtained by exact methods in 1D [138, 139, 140, 141, 26]. A 1D interacting two-component (spin-half) system is known to form a Luttinger liquid whose elementary excitations at zero spin-polarization are collective density fluctuations that carry only charge or spin. When transitioning from the unpolarized 1D superfluid to the 1D FFLO equivalent, the spin gap is closed and a second gapless sound mode emerges. This can be viewed as a Luttinger liquid of the excess fermions with a Fermi wavevector $k_{F\uparrow}-k_{F\downarrow}$. With spin polarization, the excitations become of mixed spin and density type, reflecting the breaking of time-reversal symmetry by the spin-polarization. For more information about the 1D FFLO state, we recommend [25, 24, 26].

Collective modes in 2D and 3D lattice systems of two-component Fermions were studied by Heikkinen and Törmä [131], using the FF ansatz. The density response was derived based on the Kadanoff-Baym approach to GRPA, as given by equation (147). The FF ansatz of the form $\Delta e^{i\mathbf{q}\cdot\mathbf{r}}$, with a fixed FF wavevector $\mathbf{q}$, makes the system asymmetric with respect to directions. In a lattice system as discussed in section 3, $\mathbf{q}$ aligned with one of the lattice axes, in most cases, minimizes the energy. It was found in [131] that this anisotropy of the order parameter leads to a collective mode dispersion relation that is different in the directions parallel and perpendicular to $\mathbf{q}$, see figure 16. This means that the speed of sound is anisotropic as well. A rigorous derivation is presented in the Ph.D. thesis of M.O.J. Heikkinen [130], section 4.2.1. which shows that the mode obtained by this approach indeed corresponds to the Nambu-Goldstone mode. Specifically, the response function is shown to diverge at zero frequency and momentum, which proves that the mode is gapless. The proof explicitly shows that the mode is gapless because the BCS gap equation is consistently included in the formalism calculating the collective mode spectrum.

Heikkinen and Törmä also considered a quasi-1D lattice to compare with the work of Edge and Cooper [136, 137]. They found two branches around $k=0$, similar to [136, 137], see figure 17. One branch corresponds to the collective mode in the higher dimensional case, while the other comes from the quasi-1D nature: quasiparticle excitations gather into a narrow stripe. The existence of quasiparticles, even at low energy, is due to the gapless nature of the FF state, however in higher dimensions, they are more dispersed and do not form a stripe. The results of [131] differ at higher momenta from those in [136, 137]: there is no gapless sound mode around $k^{*}=2(k_{F\uparrow}-k_{F\downarrow})$ because in the FF case, $|\Delta|$ is constant and the excess fermions are uniformely distributed. The backbending of the collective mode dispersion in 2D and 3D lattices at large momenta has also been interpreted as a roton-like behaviour [142, 143, 144, 145].

Radzihovsky [109, 146] derived a low-energy Landau theory for the LO state in a two-component Fermi gas and found two Goldstone modes, one corresponding to superfluidity and one that he refers to as a smectic phonon. Drawing analogy to the above discussion, the former one is related to the Anderson-Bogoliubov phonon and the latter to the existence of the excess fermions. However, Radzihovsky emphasizes the liquid-crystalline character and softness of the smectic phonon, which may in some cases lead to vanishing FFLO order parameter, or peculiar topological defects such as composite half-integer vortex-dislocation defects. In a lattice, the rotational symmetry is broken and such features are not expected to appear, except perhaps for special conditions such as low filling where the dispersion felt by the particles is almost rotationally symmetric.

The Anderson-Bogoliubov phonon can be understood as fluctuation of the phase of the order parameter. Fluctuation of the order parameter amplitude is referred to as the pair vibration mode, or the Higgs mode. This has been theoretically considered in the ultracold Fermi gas context in refs. [147, 148, 149, 150, 151] for trapped gases. The Leggett mode is a collective mode associated with pairing fluctuations of multiple bands, and a harmonic trap analogue of such mode was found in [148]. We are not aware of any works on Higgs or Leggett modes in the FFLO state, at least not in the ultracold gases context.

7 Beyond mean-field studies of the FFLO state in 2D and 3D lattices

The mean-field approximation is known to be typically good in the case of weak interactions and high dimensions when the system is not likely to be strongly correlated. Pairing in the FFLO state is obviously stronger if the interaction strength is increased. Moreover, there is evidence that reduced dimensionality helps to stabilize the FFLO state. Thus both from the conceptual and practical (e.g., higher critical temperatures) point of view, it is of interest to explore the FFLO state in strongly interacting Fermi gases and in reduced dimensions. Correspondingly because quantum and thermal fluctuations are typically stronger in low dimensions and for strong interactions, one has to harness beyond mean-field methods for the theoretical description of the state. Even in the weakly interacting limit in 3D, the subtle nature of the FFLO order parameter (c.f. the discussion on fluctuations in section 6), deserves a closer look with beyond mean-field methods to test its robustness against both quantum and thermal fluctuations.

There are a handful of beyond mean-field studies of the FFLO state in lattices, in 2D, 3D, and quasi-1D systems [152, 153, 154, 155, 156, 114, 157].444In 1D, of course, the situation is competely different since exact methods such as Bethe ansatz can be applied; we do not review the works that focus on 1D FFLO state [24]. Of these, three references [153, 155, 156] are based on dynamical mean-field theory (DMFT), and the other works [154, 114, 157] on quantum Monte Carlo (QMC). The early inspirational work [152] presents an in-principle exact solution of the problem using the algebraic techniques of the Richardson-Gaudin model, but in practice they resort to solving the exact equations for a small (6x6) lattice, also utilizing QMC for the search of candidate ground states. The obtained zero-temperature phase diagram does not easily compare to later mean-field and beyond-mean-field work, possibly due to the small system size. They predict the possibility of a breached FFLO state.

To be able to describe how the results obtained by DMFT go beyond mean field thery, we first sketch the basic idea of how to apply the DMFT method to study fermion superfluidity and especially FFLO (for a general review on the DMFT method, see [158]). We then review results on the FFLO state in the 1D-3D crossover using the method in section 7.2. In section 7.3, we discuss results on FFLO physics in 2D lattices obtained utilizing QMC.

7.1 Dynamical mean-field theory (DMFT) for studies of the FFLO state

DMFT [158] is a non-perturbative method where the many-body problem is self-consistently mapped to a quantum impurity problem, as described in figure 18. It is assumed that the self-energy of the system is local (here $j$ and $l$ label lattice sites, and $\omega_{n}$ is the Matsubara frequency)

[TABLE]

This assumption can be relaxed by using a cluster version of DMFT where the self-energy can be non-local within the cluster, as discussed later. The DMFT method as described in figure 18(a) applies to spatially homogeneous lattice systems, meaning the absence of external potentials, such as a harmonic trap potential typical in ultracold gases. This assumption of homogeneity leads also to a spatically homogeneous order parameter. To study the FFLO state where the order parameter is non-uniform in space, one needs a real-space version of DMFT [159, 160, 161, 153]. The idea of real-space DMFT is depicted in figure 18(b). We briefly present how fermionic superfluidity and the FFLO state can be described by the real-space DMFT method.

To describe superfluid symmetry breaking in an interacting fermion system, the anomalous pairing correlations need to be incorporated in to DMFT. This can be done using the Nambu formalism. The Nambu Green’s function is defined as a two by two matrix

[TABLE]

This is equivalent to the $\sigma=\{1,4\}$ matrix block of Section 5. The indices $j,l$ label the lattice sites. Here for the anomalous Green’s function, we have $F^{\dagger}_{jl}(\tau)=[F_{jl}(\tau)]^{\dagger}$, and $[{\hat{O}}(\tau)]^{\dagger}={\hat{O}}^{\dagger}(-\tau)$ for the Matsubara time evolution. The lattice dispersion in the non-interacting case enters the non-interacting Green’s function, which in the Nambu formalism is of the form

[TABLE]

where $G_{\sigma,jl}^{0}(i\omega_{n})=1/\left[i\omega_{n}+t_{jl}+\left(\mu_{\sigma}-V_{j,\sigma}\right)\delta_{jl}\right]$, $t_{jl}$ is the hopping energy and $V_{j,\sigma}$ the single-particle potential that can correspond to the trap potential, for instance. The order parameter of the superfluid phase at a given site is of the usual form

[TABLE]

where $U<0$ for attractive interactions.

The Dyson equation gives the relation between the interacting and non-interacting Green’s functions and the self-energy, and is of the form (in general, the Dyson equation can have terms corresponding to spin-flip processes and spin-dependent interactions, but this form is sufficient to describe the two-component ultracold Fermi gases typically considered for the FFLO state):

[TABLE]

In DMFT, the self-energy is approximated to be site-diagonal, $\boldsymbol{\Sigma}_{jl}(i\omega_{n})\approx\delta_{jl}\boldsymbol{\Sigma}_{j}(i\omega_{n})$. Direct matrix inversion is used for solving the Dyson equation in real-space DMFT. Now all this defines the full lattice problem. In DMFT, one divides the problem into a local impurity problem and a bath that represents the rest of the lattice. We now proceed to first discuss the impurity problem and then how to connect it to the bath.

The quantum impurity problem of lattice site $j$ can be defined as the problem of solving the full local Green’s function $\boldsymbol{\mathcal{G}}_{j}$ for the interaction strength $U$ and a local bath Green’s function $\boldsymbol{\mathcal{G}}_{\mathrm{Bath},j}$. This means summation of all connected diagrams of $\boldsymbol{\mathcal{G}}_{\mathrm{Bath},j}$ and $U$. Once the full local Green’s function $\boldsymbol{\mathcal{G}}_{j}$ is obtained, one can obtain the local self-energy $\boldsymbol{\Sigma}_{j}(i\omega_{n})$ from the Dyson equation

[TABLE]

Next, the quantum impurity problem (158) should be connected to the full lattice problem (157). This is done by forming the bath of the impurity problem via removing the self-energy of site $j$ from the local Green’s function $\boldsymbol{G}_{jj}(i\omega_{n})$. This means that the bath Green’s function is given by the equation

[TABLE]

The right hand side of this equation specifically contains the inverse of the local component of $\boldsymbol{G}$ and not the $jj$-element of $\boldsymbol{G}^{-1}$. This bath Green’s function includes the effect of interactions on all other lattice sites except the site $j$.

The solution of the full problem is found iteratively, for instance by making an initial guess for the self-energy at all sites. Combined with the fact that the non-interacting lattice Green’s function is known, this allows us to calculate the full lattice Green’s function and the bath Green’s function. With the bath Green’s function thus determined, one can solve the impurity problems at each site, and thereby obtain improved self-energies. The iteration is repeated until convergence is found.

The quantum impurity problem at the heart of the DMFT method, basically obtaining $\boldsymbol{\mathcal{G}}_{j}$, can be solved by various means. For instance, exact diagonalization (ED) utilizing an Anderson impurity problem, QMC, or density matrix renormalization group (DMRG) methods can be used. Typically ED and DMRG are suited for zero or low temperature descriptions, while QMC is most efficient at intermediate temperatures. The use of a QMC solver in DMFT thus captures, locally, both quantum and thermal fluctuations.

DMFT goes beyond the mean-field approximation because it solves the local lattice problem exactly. To describe interactions in the lattice system, all local diagrams are included. In some problems, however, non-local quantum correlations may become significant. These correlations can be addressed by a cluster version of DMFT, described in figure 18(c). There are many variants of cluster DMFT, the most widespread being cellular dynamical mean field theory (CDMFT) and dynamical cluster approximation (DCA). In CDMFT, the idea is to define the cluster on the real-space lattice, while DCA can be characterized as a coarse graining in the momentum space. The cluster selection of the CDMFT method breaks the translation invariance of the original lattice, whereas DCA enforces the translation invariance.

We now discuss the results obtained in beyond mean-field studies of the FFLO state using the DMFT and QMC methods. The first work using DMFT to study the FFLO state by Kim et al. [153] is discussed in section 8 since it specifically addresses issues related to the trapping geometry. Another FFLO study in the context of ultracold gases going beyond mean field was done by Wolak et al. [154] using determinant QMC for 2D square lattice. This work, as discussed in section 7.3 with a recent related work by Gukelberger et al. [114]. Next, we discuss DMFT studies that consider the FFLO state, particularly in the context of 1D-3D crossover.

7.2 The FFLO state in the 1D-3D crossover

In 1D, exact methods can be used to show [24] that at zero temperature, the 1D equivalent of the FFLO state exists for any interaction and any polarization $P=(N_{\uparrow}-N_{\downarrow})/(N_{\uparrow}+N_{\downarrow})=(n_{\uparrow}-n_{\downarrow})/(n_{\uparrow}+n_{\downarrow})$, where $N_{\sigma}$ and $n_{\sigma}$ are the particle numbers and densities of the two spin components. The 1D equivalent of the FFLO state is extremely stable due to perfect nesting: the Fermi surfaces in 1D are just points, see the discussion in section 3. On the other hand, no long range order exists in 1D, and the polynomially decaying FFLO correlations are highly vulnerable to thermal fluctuations. In higher dimensions, the parameter regime (polarization, interaction, chemical potentials, etc.) for a stable FFLO state is extremely small, although lattice geometries provide a bigger parameter window, as explained in section 3. Thermal fluctuations, in contrast, are known to have less effect on quantum states the higher the dimension is, and true long range order is possible in 3D. This has inspired the hypothesis that an intermediate dimensionality between 1D and 3D might optimally combine the large parameter regime of FFLO in 1D with the stabilization towards finite temperature typical for higher dimensions. Mean-field [162, 116] and effective field theory [139] calculations for coupled 1D tubes have suggested that this could be possible. The question has also been considered by DMRG calculations on a Hubbard-ladder [163].

The behaviour of the FFLO state in 1D-3D crossover in lattices was considered with the DMFT method by Kim and Törmä [164] and Heikkinen et al. [155, 156]. In [164], the attractive Hubbard model combined with a trap in one direction was approached by DMFT with an exact diagonalization impurity solver, allowing access to zero and low temperatures. The interpolation between 1D and 3D was done by keeping tunneling parallel to one fixed direction (the direction of the trap potential), $t_{\parallel}$, while tuning the tunneling in the perpendicular direction, $t_{\perp}$. The interaction strength $U$ was chosen for each value of the transverse hopping $t_{\perp}$ such that it corresponded to the lattice equivalent of the unitarity limit [165]. For small $t_{\perp}\leq 0.3t_{\parallel}$, it was found that the shell structure within the trap was of the 1D type: the FFLO state in the middle of the trap surrounded by balanced superfluid areas. For larger coupling between the chains, the 3D-like behaviour of balanced superfluid in the central trap regions surrounded by FFLO areas was found, in contrast. This 1D-3D crossover behaviour occurred for small polarizations, while for larger polarizations, the FFLO state exists everywhere in the trap, and was replaced by the normal state for polarizations beyond the critical value. A clear difference to the previous mean-field predictions was found: the DMFT results predicted the FFLO state to be stable for a wide range of $t_{\perp}$ values throughout the crossover, extending well in to the 3D-like regime. This contrasts with having an optimal region near the transition from 1D to 3D, as suggested in articles presenting mean field [162] and effective field theory [139] studies. There can be several reasons for this difference, first that the work [164] includes beyond mean-field fluctuations. Second it considers coupled 1D lattice systems while the earlier mean-field works studied coupled 1D tubes (continuum systems), thus nesting effects are present in [164] in all dimensions.

Heikkinen et al. [155] extended these studies to finite temperatures by using a QMC solver in the DMFT impurity problem. The configuration and the parameters were similar to reference [164] described above. The total particle number was held approximately constant while the polarization was varied. The main results are summarized in figure 19. Again the FFLO state is found to be stable over a large range of parameters throughout the crossover. The critical temperature for the FFLO state is found to be about one-third of the critical temperature of the balanced case, see figure 19 (a), (d), (g).

Both references [164] and [155] showed that at the turning point between 1D-like and 3D-like behavior, at $t_{\perp}=0.3t_{\parallel}$, the FFLO order parameter demonstrates an intriguing behavior: it has a remarkably constant oscillation amplitude despite the overall trapping potential, see figure 20. This is because in the 1D limit, the FFLO oscillations occur at the edges, and in the 3D limit at the center of the trap. Thus in the intermediate regime both tendencies balance to produce a flat profile [162]. There is thus a sweet spot in the intermediate regime where the FFLO order parameter in a trapped system resembles the one expected for a homogeneous system. Moreover many suggested experimental probes of the FFLO state give more clear signatures in case of homogeneous FFLO oscillations compared to a typical trapped system FFLO order parameter profile.

Since fluctuations at the level of single sites (as captured by DMFT) already change the picture in comparison to mean-field results, it is essential to ask what might non-local fluctuations do to the stability of the FFLO state. This question was considered in [156] using a cluster version of DMFT. The choice of the cluster is illustrated in figure 21: clusters consisting of 1D chains of 36-42 sites were solved exactly, and the chains were connected in two perpendicular directions, and usual DMFT was applied to the colletion of the chains. The self-energy is of the form

[TABLE]

that is, it is block diagonal in the interchain index $l$, while all the indices $i$ and $i^{\prime}$ run through all the self-energy terms, both local and non-local, within the 1D chain. Periodic boundary conditions were applied in the 1D cluster but otherwise arbitrary spatial dependence of the order parameter was allowed. This approach is different from the DCA version of DMFT where translational invariance is imposed, but it also deviates from the standard CDMFT with fixed cluster size. In this way, a spontaneous breaking of the translational invariance as expected for the FFLO state was allowed, but it was not imposed by construction as in CDMFT.

As the main result, the FFLO state was found to be stable even in the presence of non-local quantum fluctuations. The critical temperatures for both the FFLO and BCS phases were found to be smaller than those predicted by mean-field theory, see figure 22. There was a clear difference between the cluster and single-site, real-space DMFT results: the cluster version captured the increasing fluctuations when the 1D limit was approached, while the critical temperature in the single-site version remained rather constant when the interchain tunneling was decreased to realize quasi-1D systems. Notably, even the BCS critical temperature was strongly affected by non-local quantum fluctuations at the quasi-1D limit. This indicates that the cluster DMFT approach used was capable of describing the decreasing robustness of long-range coherence when going towards the 1D limit and the related drop in critical temperature.

Hulet, Mueller and coworkers [22] have experimentally studied an imbalanced Fermi gas in a continuum 1D system and found consistency with the FFLO state, although no smoking gun signature. We do not discuss this important experiment here because it belongs to the topic of a 1D FFLO phase, which was excluded from this review. The same experimental group recently studied the 1D-3D crossover [23], note that the difference to the crossover studied in this section is that their 1D system is a continuum tube, not a 1D lattice system as here. The crossover was driven by tuning the tunneling strength between the 1D tubes. They observed that indeed there is a critical tunneling energy separating the 1D and 3D regimes: in the 1D regime, the polarized superfluid appeared in the middle of the tube, while in the 3D case it was found at the edges, consistent with our discussion of what happens in a corresponding lattice 1D-3D crossover. Direct signatures of the FFLO state were not observed.

7.3 Quantum Monte Carlo studies of the FFLO state in 2D square lattices

As discussed in section 3, in 2D the nesting effects near van Hove singularities are predicted to be particularly strong. However, quasi-long range order should be possible in 2D, in contrast to 1D, via BKT type mechanisms. This makes 2D lattices a promising setting for the FFLO state but also imposes demands on the theoretical treatment since fluctuations should be important. In 2012, Wolak et al. [154] used determinant QMC to tackle this question and found for low fillings that the pair correlator indeed had a feature at finite momentum, suggesting instability towards the FFLO phase. Simultaneously they calculated normalized double occupancy that describes pairing. At higher temperatures, the peak in the pair correlator reduces to zero instead of the finite value, while the normalized double occupancy stays constant. This was interpreted as a polarized paired phase (PPP), without the broken translational symmetry typical for the FFLO state. However, when the authors calculated the pairing susceptibility, they found no $s$-wave superfluidity in the polarized case for the parameter regime they were able to study. The fermion sign problem limited their study to low densities and temperatures at $0.1-0.2t$, where $t$ is the hopping, or higher. Thus the existence of an FFLO superfluid remained to be confirmed.

Gukelberger et al. [114] used unbiased diagrammatic QMC simulations and determined the FFLO superfluid phase boundaries from the divergence of the pair susceptibility. For attractive interactions of strength of about half the band width, they found a reasonably large parameter area in the polarization - temperature phase diagram where the FFLO state is predicted. A zero momentum pairing phase is found, and the high temperature - high polarization phase is identified as a Fermi liquid, see figure 23. The critical temperatures for the FFLO phase were found to be optimally about half of the critical temperature of the balanced case, in accordance with the DMFT studies [155, 156]. For very large polarizations, a triplet superfluid was found, but at an exponentially small temperature.

Interestingly the phase diagram matches the mean-field result with simple scaling of the temperatures and polarizations, as was found by Baarsma and Törmä [66]. In that work, the two-component Fermi gas in a 2D square lattice was considered using the mean-field ansatz, which describes the LO state, see section 3. As expected, there is a clear quantitative difference between the mean-field and the QMC results, the latter producing smaller critical temperatures and polarizations by a factor of about four. However, when the temperatures are rescaled by the scaling factor $T^{\prime}$ between the two critical temperatures at $P=0$, and the polarizations similarly by the ratio of polarizations at a low temperature $P^{\prime}$, the two phase diagrams match remarkably well for all values of $T$ and $P$, see figure 24.

In [157] a hybrid of mean-field theory and the QMC approach that aims to capture thermal fluctuations of the order parameter was used. Critical temperatures were found to be about 4-5 times smaller than predicted by mean field theory, similar to the difference in the QMC and mean-field calculations in [66] and [114] discussed above. However the temperature-magnetization (magnetization divided by the density equals polarization) phase diagram obtained by their hybrid method predicts the LO region to be much smaller and in a different place than the exact QMC results of [114], although direct comparison is not possible due to the near half-filling and quarter filling densities used in [157] and [114], respectively. Similarly to [154], the authors of [157] find a finite momentum signature in the pairing structure factor well above $T_{c}$.

Note that since the focus of this review is on the FFLO state in lattices, we do not discuss here beyond mean-field theory studies treating continuum imbalanced Fermi gases. Likewise, we will not give a comprehensive review of beyond mean-field treatments of trapped imbalanced gases, but a few examples of such studies are discussed in section 8.

8 FFLO in trapped gases and phase separation

The effect of trapping potential and its relation with FFLO has already been reviewed in [21], so here we give only a brief description and then concentrate on topics not covered by their review. Particularly we consider here the effect of trap elongation, beyond mean-field theories, and we give a more detailed analysis of the Bogoliubov-de Gennes (BdG) method. Most of the discussion in this section involves gases without optical lattices, although a lattice can be used for producing highly elongated potentials, and the BdG method can be used also in lattices [166].

Ultracold atom gas experiments are always conducted in some kind of external trapping potential. This has profound effects on the behaviour of the gas as the translational symmetry is broken and atom gas density becomes position dependent. Furthermore it has particularly dramatic effects on spin-imbalanced gases, as the position dependence of the trap naturally drives a phase separation in which, in 3D, the edges of the trap become fully polarized while the center of the trap hosts both atomic species.

Fortunately, the external trapping potential can be incorporated in the mean-field theory. The external trapping potential in the Hamiltonian (1) is described by the term

[TABLE]

In practice the trapping potential is often assumed to be spherically symmetric $V({\bf r})=V(r)$, or at least cylindrically symmetric. The mean-field approximation proceeds as outlined in section 2, with the main differences being that the order parameter becomes position dependent $\Delta({\bf r})$ with a priori unknown position dependence, and that the Hartree energy shift, which follows the shape of the density profile $n_{\sigma}({\bf r})$, can no longer be simply included in the chemical potential. The mean-field Hamiltonian in a trap becomes

[TABLE]

where the local single-particle Hamiltonian is defined as

[TABLE]

In practice the Hartree shift, the $V_{0}n_{-\sigma}$-term in the single-particle Hamiltonian, is usually completely neglected when studying trapped gases as it is ill-behaved in the strongly interacting limit. For strong attractive interactions $|k_{\mathrm{F}}a|>\frac{4}{3\pi}$, the Hartree shift predicts an unphysical collapse of the gas. This instability is driven by the uncontrolled increase of the Hartree shift when the density of the gas increases. This critical bound is easily obtained by considering the total energy cost of adding a particle: the kinetic energy equals the Fermi energy $E_{\mathrm{F}}=\frac{\hbar^{2}k_{\mathrm{F}}^{2}}{2m}$ and the (Hartree) interaction energy gain is $\frac{4\pi\hbar^{2}a}{m}n$, where $n=k_{\mathrm{F}}^{3}/(6\pi^{2})$ is the density of one of the spin components. Together these energies yield

[TABLE]

As seen from the equation, the interaction energy increases faster than the kinetic energy cost as function of the Fermi momentum $k_{\mathrm{F}}$, and hence the density $n$. In particular when the factor $1+\frac{4}{3\pi}k_{\mathrm{F}}a$ becomes negative, the total energy cost of adding a particle becomes negative. It is thus energetically favorable for the gas to collapse when attractive interaction becomes strong: $\left|k_{\mathrm{F}}a\right|>\frac{3\pi}{4}$. Such collapse is completely unphysical but it plagues mean-field theories and the BdG method in particular. The problem can be solved by using correct (i.e., determined by Monte Carlo methods or experimentally) energy functional for the Hartree energy shift [167], but often the Hartree energy shift is simply neglected in mean-field studies.

The mean-field Hamiltonian equation (163) can be diagonalized formally in the same way as in the homogeneous case, leading eventually to the BdG equations:

[TABLE]

where $i$ is the eigenstate index. The physical interpretation of each of the elements in these BdG equations are the same as in the homogeneous case discussed in section 2.4. However, the apparent simplicity of these equations is misleading. The position dependence of both the order parameter $\Delta({\bf r})$ and the differential operators in the single-particle Hamiltonian $K_{\sigma}({\bf r})$ require one to solve either coupled differential equations or turn each of the elements into matrices. Notice that this is also the case of a homogeneous LO-state, which does not allow separation into different momentum states ${\bf k}$, but instead involves a coupling between all states, see section 2.4.

Before analyzing the BdG equations in more detail, we first briefly review the use of the local density approximation for solving equation (166). We also shortly review other more advanced theories used for describing FFLO physics in trapped systems, not discussed in earlier reviews.

8.1 Local density approximation

In the local density approximation (LDA), the gas is assumed to be locally homogeneous, with the external potential $V_{\mathrm{ext}}({\bf r})$ providing only a position-dependent energy shift. This assumption is expected to be valid for sufficiently large systems in which all physical quantities are only slowly varying with position. With such an assumption, the position coordinate ${\bf r}$ can be neglected in BdG equations (166), becoming only a parameter that determines the local chemical potential

[TABLE]

and the BdG equations become identical to the ones in the homogeneous case

[TABLE]

where $\epsilon_{k}=\frac{\hbar^{2}{\bf k}^{2}}{2m}$. However this needs to be solved separately for each different value of the (local) chemical potential $\tilde{\mu}_{\sigma}$. Thus the LDA to the BdG equations means solving the homogeneous mean-field theory locally for each spatial point ${\bf r}$, with chemical potentials that are shifted by the value of the external potential $V_{\mathrm{ext}}({\bf r})$.

Assuming that the two spin components feel the same trapping potential, the chemical potential difference $\delta\mu=\tilde{\mu}_{\uparrow}-\tilde{\mu}_{\downarrow}$ does not depend on the position, whereas the average chemical potential $\mu_{\mathrm{avg}}({\bf r})=(\tilde{\mu}_{{\bf r},\uparrow}+\tilde{\mu}_{{\bf r},\downarrow})/2$ does. To solve the mean-field BCS theory for a trapped gas in the LDA, it is thus sufficient to solve the homogeneous order parameter $\Delta$ and the densities $n_{\sigma}$ for a range of average chemical potentials. Once one has the functions $\Delta(\mu_{\mathrm{avg}})$, $n_{\sigma}(\mu_{\mathrm{avg}})$ mapped for a fixed chemical potential difference $\delta\mu$ and temperature $T$, the local densities $n_{\sigma,{\bf r}}$ and order parameter $\Delta_{\bf r}$ can immediately be determined as

[TABLE]

and

[TABLE]

The trapped gas thus effectively spans a range of parameters in the homogeneous gas phase diagram, and from the phase diagram one can immediately read the corresponding layered phase structure of the gas, see figure 25. In practice, one needs to iterate the chemical potential difference $\delta\mu$ if one wants to solve the model for a specific number of atoms.

As described by the LDA, the inhomogeneous potential drives the system to phase separation [72, 168]. In 3D, at the edge one has a single-component gas in a normal state, whereas the center of the trap involves an equal density BCS-superfluid. This general property is seen in experiments [17, 169], but also in practically all theories using LDA, both mean-field [21] and beyond mean-field theories [170, 171].

Indeed LDA is very versatile, as it is not limited to mean-field theories, but it can be applied to any theory for homogeneous gases. The homogeneous mean-field results suggest that the FFLO phase involves only a sliver of parameter space between phase separated and normal states in the homogeneous gas phase diagram [108]. This translates into a narrow shell of FFLO state in a trapped gas for specific parameters for atom numbers, interaction strength and temperature. More involved beyond mean-field theories have similarly been applied to trapped gases using the LDA [170, 171]. These works considered the role played by pseudogap physics at temperatures above the critical temperature, but they did not consider the FFLO state.

LDA involves sharp phase boundaries, and the zero-temperature mean-field theory in particular produces sharp density changes at the superfluid-normal state interfaces as seen in figure 26. Such rapid variations violate the assumptions of slowly varying potentials in LDA. Improvements to LDA involve gradient corrections that describe the surface tension or domain wall energies [173, 174, 175, 176, 21]. When an energy cost is associated with the superfluid-normal interfaces, rapid variations of the order parameter $\Delta({\bf r})$ and densities $n_{\sigma}({\bf r})$ are suppressed. Mean-field theories showed that this should lead to deformation of the superfluid core to minimize the size of the superfluid-normal interfaces. These effects were hoped to solve discrepancies observed in the experiments by Rice [18] and MIT groups [17, 169]. The surface tension was too weak [176] to provide an explanation. Instead it was suggested [177] and later experimentally confirmed [178] that the features in the Rice experiment were due to the evaporative cooling from highly elongated trap, resulting in a depolarized nonequilibrium state at the center of the trap. However the concept of surface tension is useful for understanding deformations of the density profiles predicted by BdG equations in elongated traps, which will be discussed below.

8.2 Bogoliubov-de Gennes equations in spherical traps

The BdG equations (166) can be solved also without the LDA. However this requires either the solution of coupled differential equations or expressing the different elements in matrix form. The first approach using differential equations was reviewed by Gubbels and Stoof [21], so here we review the matrix formalism, which is the approach used by a number of authors. We focus on the case of spin-imbalanced Fermi gas in a spherically symmetric harmonic trap [179, 180, 181], but the formalism itself also works with minor modifications in non-spherically symmetric traps [153, 182, 183]. Also, the early works on FFLO by Mizushima et al. involved BdG equations in cylindrically symmetric system with periodic boundary conditions [184].

The standard approach is to expand single-particle operators in the eigenbasis of the trapping potential. This results in a diagonal form for the single-particle part of the mean-field Hamiltonian. Another possible choice is to use an evolving quasiparticle basis, in which the eigenbasis of the previous iteration of the gap equation $\left(u_{i}({\bf r}),v_{i}({\bf r})\right)^{T}$ is used as the basis for the next iteration.

In the case of a spherically symmetric harmonic trap, $V(r)=\frac{1}{2}m\omega^{2}r^{2}$, where $\omega$ is the trapping frequency, the single-particle creation and annihilation operators are expressed in the basis of 3D harmonic oscillator eigenstates

[TABLE]

where the quantum numbers $n,l,m$ are: the radial quantum number $n$ counting the nodes in the radial function, the orbital angular momentum $l$, and the projected angular momentum $m$. The angular part of the eigenstate $Y_{lm}(\hat{r})$ is the spherical harmonic with the radial unit vector $\hat{r}=(\theta,\varphi)$. The radial part of the wavefunction is

[TABLE]

where $\bar{r}=r/r_{\mathrm{osc}}$, the harmonic oscillator length $r_{\mathrm{osc}}=[\hbar/(m\omega)]^{1/2},$ and $L_{n}^{l+1/2}(\bar{r}^{2})$ is the associated Laguerre polynomial. The corresponding single-particle eigenenergies are $\varepsilon_{nl}=\hbar\omega\left(2n+l+3/2\right)$, with different $m$-states being degenerate. While such analytical formulas for the single-particle eigenstates are useful, the BdG method is also quite applicable with numerically solved eigenstates.

The spherical symmetry, which is assumed to also hold for the order parameter $\Delta(r)$, allows performing the angular integrations. This makes the mean-field Hamiltonian in equation (163) diagonal in $m$-quantum numbers and the resulting radial 1D Hamiltonian $(2l+1)$-fold degenerate. Dropping the redundant $m$-quantum number, the Hamiltonian becomes

[TABLE]

Here the Hartree interaction is described by the elements

[TABLE]

and the pairing field is described by

[TABLE]

We note that the $\sigma$ dependence of the Hartree term is due to the population imbalance which implies that the corrections are different for the two components. The density of $\sigma$ atoms is

[TABLE]

and the order parameter is

[TABLE]

The additional factor $(2l+1)$ comes from the degeneracy of the $m$-states.

The nature of the four matrix entries in the BdG equations (166) can now be explained. The upper left corner element, $\hat{K}_{\uparrow}({\bf r})-\mu_{\uparrow}$ becomes a matrix with diagonal elements $\epsilon_{nl}+J_{nn\downarrow}^{l}-\mu_{\uparrow}$ describing single-particle energies and Hartree shifts, and off-diagonal elements $J_{nn^{\prime}\downarrow}^{l}$ provide transitions between the $n$-levels due to a position dependent Hartree shift. Similarly, the upper right corner element, $\Delta({\bf r})$, becomes a matrix with elements $F_{nn^{\prime}}^{l}$. The diagonal elements of these blocks, $F_{nn}^{l}$, describe intrashell pairing in a manner analogous to the zero center-of-mass momentum ${\bf k,-k}$ pairing in a homogeneous system. Similarly, the off-diagonal elements of the block, $F_{nn^{\prime}}^{l}$ describe pairing between different harmonic trap levels $n$ and $n^{\prime}$ (the $l$-quantum number being the same for both). Solving the BdG equations involves diagonalizing the resulting matrix, yielding the eigenenergies $E_{i}$ and the eigenstates $(u_{i},v_{i})$. Obviously one needs to impose a cutoff to keep the dimensionality of the matrices finite, but even more importantly since the gap equation is ultraviolet divergent. The regularization of the gap equation, and the associated renormalization of the bare interaction $U$ is generally done using LDA. However it is not clear how renormalization should be done when studying 1D-3D crossover, as LDA based schemes are either for 1D (in which case regularization/renormalization is not needed) or for 3D. For the present purposes, it is enough to simply assume some large energy cutoff $\epsilon_{\mathrm{c}}$. With the energy spectrum of the 3D harmonic oscillator, the energy cutoff implies an angular momentum $l$-dependent cutoff for the radial quantum number $n<N_{l}=\epsilon_{\mathrm{c}}/(2\hbar\omega)-l/2$.

In spherically symmetric systems, the matrix equation has a block structure, since the Hamiltonian can be expressed as a sum of mutually commuting terms for different $l$-quantum numbers. One can thus express it as $H_{\mathrm{MF}}=\sum_{l}H_{\mathrm{MF}}^{l}$, and write a separate set of BdG equations for each $H_{\mathrm{MF}}^{l}$ block. Diagonalizing these yields eigenenergies $E_{jl}$ and the corresponding eigenstates $(u_{j}^{l},v_{j}^{l})^{T}$, where the $u_{j}^{l}$ and $v_{j}^{l}$ are now vectors in the initial single-particle basis and $j$-index enumerates the different eigenstates of the sub-Hamiltonian $H_{\mathrm{MF}}^{l}$. Notice that the $j$-index also includes the $+/-$-quasiparticle branches introduced in section 2.4, and consequently the sign convention for the quasiparticle energies $E_{jl}$ is altered from the one in equation (55). Here half of the eigenstates have negative energies, and the zero-temperature ground state corresponds to the state in which all quasiparticle states with negative energy states are occupied and others empty, whereas in section 2.4, the corresponding quasiparticle states describe particle and hole excitations.

The atom densities can now be written as

[TABLE]

where the Fermi distribution $n_{\mathrm{F}}(E)=1/(1+e^{E/k_{\mathrm{B}}T})$ and $u_{jn}^{l}$-factors describe the weight of the spin-$\uparrow$ component in each state $nl$. Likewise, the density of atoms in $\downarrow$ state is

[TABLE]

where the sign in the Fermi function is changed because the $\downarrow$ component of the eigenstates correspond to holes as compared to the particles in the $\uparrow$ components. The order parameter is given by

[TABLE]

BdG equations are ideally suited for describing FFLO-type physics, since they do not in principle restrict the spatial dependence of the order parameter. However in practice, most of the studies have assumed that the order parameter is real and spherically symmetric $\Delta({\bf r})=\Delta(r)$. These assumptions neglect the possibility of FF-type superfluids in which the phase of the order parameter is changing, and the angular FFLO-state in which the order parameter varies not only radially but also around the trap. However, later BdG studies for gases in elongated (non-spherical) traps seem to confirm these early assumptions. On the other hand, angular FFLO has been studied in toroidal traps [185] and in those systems even the complex FF-type states can be stabilized when the system is rotated [186]. Angular FFLO seems to also be present in 2D lattices at high filling [166].

8.3 The FFLO state and proximity effects in trapped 3D gases

The BdG theory yields density and order parameter profiles that are in good qualitative agreement over a wide range of parameters with experimental results but also with the LDA [17, 169, 180, 181, 187, 172, 21]. Figure 27 shows how increasing the system size, or atom numbers, produces profiles that are increasingly similar to the LDA results. In particular, BdG equations provide the correct general features of a balanced BCS-type superfluid in the center of the trap and fully polarized edges.

However between the balanced superfluid in the center and unpaired normal fluid at the edge, BdG equations provide a narrow region with an oscillating order parameter that has been associated with the FFLO-state. Despite the similarity with the usual FFLO state, the oscillations are present whenever one has a superfluid-polarized normal state interface [181]. On the other hand, this interface region becomes relatively narrower when the size of the system is increased. The interface region and the associated density and order parameter modulations can be understood as a proximity effect [187, 172, 21]. The proximity effect is caused at the interface of the two phases by the correlations inside the superfluid penetrating the partially polarized normal gas. The effect has been studied in interfaces of superconductors and ferromagnets in solid state systems [188]. Indeed as the center of the trap is a balanced BCS-type superfluid, and the edge is a polarized normal gas, one does expect superfluid correlations to reach into the normal gas. This manifests as oscillations of the order parameter, with the period of the oscillations determined by the difference in the local Fermi momenta of the two atomic species, in this sense similar to the FFLO state. Similar order parameter oscillations have also been predicted for dipolar gases, in which the $p$-wave interaction channel itself provides explicit intershell coupling [189, 190].

The penetration depth has been shown to be independent of the system size [172]. On one hand, this implies that the region with order parameter oscillations cannot be considered to form a separate shell of an FFLO-type superfluid phase, as the relative size of the region shrinks when the number of atoms is increased. On the other hand, it is not a finite size effect or an artifact of the BdG method, but rather an interface effect that can be understood through the surface tension discussed above.

8.4 FFLO in elongated traps and dimensional crossover

In the case of elongated traps, as is the case of actual experiments, the external trapping potential becomes

[TABLE]

where $\omega_{\rho}$ and $\omega_{z}$ are the trapping frequencies in the radial and longitudinal directions, respectively. For $\omega_{z}<\omega_{\rho}$, the trap is elongated in a cigar-like fashion, with the trap aspect ratio defined as $\lambda=\omega_{\rho}/\omega_{z}$.

The elongation has trivial effects when using LDA, since the local value of the potential $V_{\mathrm{ext}}({\bf r})$ uniquely determines the state in point ${\bf r}$ and the state is the same along an equipotential surface. This also implies that all phase boundaries and density profiles in the trap follow the shape of the trap, and particularly, they have the same aspect ratios as the trap itself. As discussed above, gradient corrections to the LDA change the picture by introducing energy costs to domain walls. These surface tension effects are most naturally included in the BdG method and, as will be seen, the BdG method predicts significant deformation for the minority component atoms for highly elongated traps. Interestingly more advanced theories, such as asymmetric superfluid local density approximation (ASLDA) density functional theory study by Pei et al. [191] (see also [167]) and dynamical mean-field theory (DMFT) in [153], did not find similar elongation of the minority component density. With these caveats, we can analyze the effect of elongation on BdG mean-field theory.

Unfortunately breaking the spherical symmetry increases the numerical complexity of the BdG method, as the block diagonal form of the Hamiltonian is lost. Actually the early works on the FFLO-state [184, 192] involved BdG equations applied to traps with cylindrical symmetry but with periodic boundary conditions in the longitudinal direction. Similarly the case of extreme elongation, in which the system is predominantly 1D, has been considered in the review by Guan, Batchelor and Lee [24]. In that limit, the spin-imbalance results in density difference at the center of the trap instead of depletion of the minority component at the edges [193]. Consequently the FFLO-type state is strongly enhanced in 1D traps, the case which has also been studied experimentally [22]. However already moderate elongation of the trap has important effects on the density and order parameter profiles. In particular, a DMFT study involving a trapped gas in a lattice in [153] predicts a polarization window for elongated traps in which the FFLO state can be stabilized, see figure 28.

The BdG method applied to elongated traps yields significant deformation in the minority component density profile, evident in figure 29. Baksmaty et al. [183] studied the effect in detail, showing that it is enhanced in highly elongated systems. At the highest elongation they studied, the radial density distributions of the minority and majority components were equal, reflecting that the minority component was, at the largest polarizations, more 2D rather than cigar-like as the majority component. The cause for the deformation is in the surface tension at the superfluid-normal gas interface, which prefers the radial density profiles of the two components to match at the expense of longitudinal direction. BdG studies for elongated traps also predict order parameter oscillations, like in spherically symmetric traps. With increasing elongation, these oscillations arranged themselves along the longitudinal axis, at the same time becoming more prominent [183]. While these interface effects are expected to also be present in larger systems, the role they play in the actual aspect ratios of the atom gases must decrease when the system size is increased.

More advanced theories, such as the ASLDA study by Pei et al. [191] and the DMFT study in [153] did not exhibit similar deformation of the minority component density profile, see Figs. 28 and 30. This would suggest that surface tension may not play a significant role after all. The DMFT work was done for 210 atoms, but the ASLDA study involved significantly larger systems. These findings underline the importance of normal state correlations and spin-dependent Hartree energy shifts. The ASLDA theory particularly is, despite the name, essentially a BdG theory but with proper Monte Carlo-based variational coefficients for the Hartree and superfluid self-energies. The lack of a proper normal state description in the BdG method is also believed to be the main reason for the anomalously high critical polarization $P_{\mathrm{c}}$ predicted by the theory.

8.5 BdG in practice

Spin-imbalanced gases are particularly difficult to solve with BdG equations. The main reason is that the nodes or zeroes of the order parameter $\Delta({\bf r})$ propagate very slowly in an iterative process. This can be easily understood from the LDA point of view: since zero order parameter $\Delta=0$ is always a solution to the gap equation, any zeroes in the order parameter profiles remains. Since BdG involves non-localized single-particle eigenstates, these nodes are still somewhat mobile. However the larger the number of atoms, the more localized the high energy eigenstates are, and the bigger the problem. Compared with balanced gases, solving spin-imbalanced BdG profiles can take up to a hundred times more iterations for convergence, and even then one has no guarantee of finding the ground state. More advanced root solving methods, such as the Broyden method, do improve the convergence and have been used by some groups [153, 183]. Also the initial state for the gap profile is a good approximation of the final solution as big changes in the profiles are often associated with appearance of zeroes in the gap profile during the iteration. However the rigidity of the order parameter zeroes in BdG equations has also been used for explaining the observations of the Rice experiment [194, 178]. In references [191, 182], initial states with different oscillating order parameter ansatzes were considered, and given sufficiently large atom numbers, the iteration was found to converge to different solutions, see figure 31. Ground state solutions were the ones with the lowest number of order parameter zeroes, but the excited states obtained through this procedure were argued to describe the findings of the Rice experiment as a signature of a metastable state. Notice that the work in [191] involved the ASLDA theory and thus goes beyond the usual mean-field BdG method.

9 FFLO and spin-orbit coupling

9.1 Spin-orbit coupling in the context of cold atomic gases

Recently there has been rapid progress both on the theory and experiments of spin-orbit coupled ultracold quantum gases (UQG). Review articles by Galitski and Spielman [42] and Dalibard et al. [41] summarize the status well, but let us also review the most salient features of these system here. After this quick introduction, we proceed to review the growing body of literature on spin-orbit coupled fermions in the context of FFLO states.

Spin-orbit coupling (SOC) is of great interest since it plays a crucial role in many condensed matter phenomena of current interest such as Majorana fermions [195] and topological insulators [196]. SOC is often the key ingredient giving rise to a non-zero topological index for insulators. This in turn implies closing the gap and hence conducting states on the surface of the insulator where two topologically different insulating states are in contact. Somewhat analoguous phenomena can appear in (gapped) superconductors, where an energy gap exists for fermionic excitations. For a topologically non-trivial superconductor, this gap can close on the surface.

Since the early days of quantum mechanics, in solid state systems, spin-orbit coupling has been a distinct possibility. If we consider an electron moving in an electric field and move to the electrons rest frame, the associated Lorentz transformation will generate a magnetic field. This magnetic field can then couple with electron spin giving rise to a spin-orbit coupling that will depend on the electron velocity. A detailed form of this coupling will be related to the directions of electron velocity and electric field, but for the special case of a motion in the xy-plane $\hat{{\bf e}}_{x}$ and electric field $E_{0}\hat{{\bf e}}_{z}$ along the z-axis, we find magnetic field

[TABLE]

where $m$ is the particle mass and $c$ is the speed of light. This magnetic field then interacts with the magnetic moment $\hat{{\bf\mu}}=\gamma\hat{{\bf\sigma}}$, where $\gamma$ is the gyromagnetic ratio and $\hat{{\bf\sigma}}$ is the spin operator. As a result a spin-orbit coupling

[TABLE]

appears, where $\alpha$ is a constant determining the strength of the coupling. Sometimes the coupling coefficient $\alpha$ is expressed in terms of the characteristic velocity $v_{R}=\alpha/\hbar$. The above coupling is called Rashba-coupling, and the system with it lacks mirror symmetry. In the absence of inversion symmetry, the spin-orbit coupling is of Dresselhaus type, namely

[TABLE]

In atomic physics, the $\sim\hat{L}\cdot\hat{S}$ term (with $\hat{L}$ the orbital angular momentum and $\hat{S}$ the spin angular momentum) in the Schrödinger equation has its origin in the non-relativistic limit of the Dirac equation. Momentum and spin dependent contributions to the system energy are called spin-orbit couplings for this reason. Notice that unlike the Zeeman interaction for example, the above examples of SOC preserve time reversal symmetry since under time reversal both spin and momentum change sign. Consequently the Kramers theorem ensures that each eigenstate is at least doubly degenerate, and one expects a degenerate ground state in a system with spin-orbit coupling (in the absence of Zeeman terms).

Electric fields inside solid state systems can be very large $\sim 10^{12}\,V/m$ (unit is volts per meter) and comparable to fields inside atoms. Such large fields can give rise to substantial spin-orbit couplings, which are clearly unfeasible in a laboratory setting. Fortunately SOCs do not have to be created using electric fields and real spins, but can also be generated using synthetic fields [41] and pseudo-spins of atomic systems. Indeed spin-orbit coupling has been recently demonstrated in UQG both for bosons [40] as well as for fermions [197, 44]. In these experiments, demonstrations relied on synthetic spin-orbit coupling generated with laser beams coupling to multilevel atoms in specific ways. These were based on earlier theoretical proposals [198, 199]. The basic idea is that different atomic internal states are coupled via a two-photon Raman transition. The atom absorbs a photon from one laser beam, and then via stimulated emission, emits it to the second beam. To conserve momentum, this process can impart momentum on the atom, and its strength depends on the Doppler shifts between the atoms and laser beams, which creates the characteristic momentum dependence of spin-orbit coupling [200].

An example of a spin-orbit coupling that has been realized is an equal mixture of Rashba and Dresselhaus couplings. In other words, the couplings have been of the type

[TABLE]

so that the full non-interacting Hamiltonian in the experiment by Lin et al. [40] is

[TABLE]

where $\hat{I}$ is the identity matrix and $\hat{\sigma}$ are the Pauli matrices. The parameters $\Omega_{z}=-g\mu_{B}B_{z}$ and $\delta=-g\mu_{B}B_{y}$ of the Hamiltonian are related to experimentally controlled magnetic field components $B_{z}$ and $B_{y}$. The strength of the SOC $\alpha=E_{L}/k_{L}$ can be expressed in terms of the recoil energy $E_{L}=\hbar^{2}k_{L}^{2}/2m$ of the laser beam photons, where $k_{L}=\sqrt{2}\pi/\lambda$ and $\lambda$ is the wavelength of the Raman beam photons. This $2\times 2$ Hamiltonian can be easily diagonalized to find the eigenstates and eigenenergies. These states form the so-called helicity basis [197]. In the above Hamiltonian, the term with $\sigma_{z}$ corresponds to out-of-plane Zeeman field while the (constant) term with $\sigma_{y}$ is an in-plane Zeeman field.

Figure 32 demonstrates the eigenstates and the role of different terms in the Hamiltonian (186) along the x-axis. With Zeeman fields and SOC coupling zero, we simply have two degenerate parabolic dispersions with their minimum at $k_{x}=0$. As we turn on the spin-orbit coupling, the two parabolic dispersions move in opposite directions so that their (degenerate) minima occur at non-zero momentum. The dispersions cross at $k=0$. As the out-of-plane Zeeman field $\Omega$ is turned on, the band crossing disappears and a gap opens between bands at $k=0$. However the ground state remains degenerate. Finally if we turn on the in-plane Zeeman field $\delta$, the ground state degeneracy is lifted. When all terms are non-zero, we have an avoided crossing between bands with a non-degenerate ground state at $k\neq 0$.

The above Hamiltonian is somewhat special but does contain many features that reoccur in various spin-orbit coupled systems. A different SOC was experimentally realized by Huang et al. [201] in a gas of fermionic potassium atoms. In their case, the spin-orbit coupling was of Dresselhaus type $H_{SO}=\lambda_{x}k_{x}\sigma_{x}+\lambda_{y}k_{y}\sigma_{y}$. Another experiment by the same group [202] studied a system with spin-orbit coupling $H_{SO}=-\alpha k_{y}\sigma_{x}+(\beta_{x}k_{x}-\beta_{y}k_{y})\sigma_{y}$. Finally very recently, Wu et al. [203] experimentally realized a 2D SOC of type $H_{SO}=h(q_{x}^{2}+q_{z}^{2})\sigma_{z}+\lambda_{so}(q_{x}\sigma_{y}+q_{z}\sigma_{x})$ in a cloud of bosonic ${}^{87}Rb$ atoms. For a more thorough theoretical discussion on how different spin-orbit couplings might be realized experimentally, see references [198, 199, 204, 41].

The SOC in equation (186) only depends on the velocity along x-axis. If we instead assume a Rashba SOC, the single particle Hamiltonian is

[TABLE]

For (pseudo) spin-$\frac{1}{2}$ particle, this can be expressed as a matrix

[TABLE]

which we can diagonalize to find the eigenenergies

[TABLE]

In this case, it is clear that the ground state has a large degeneracy since all states on the ring defined by $\sqrt{k_{x}^{2}+k_{y}^{2}}=\alpha m/\hbar^{2}$ have the same energy. In the presence of interactions, it has been suggested that this degeneracy might lead to interesting effects such as quantum many body Schrödinger cat-states [205].

Creating such SOCs as Rashba and Dresselhaus couplings analoguous with those appearing in solid state systems is useful since it allows us to explore the associated physics in a very clean and tunable environment. However the SOCs that have been experimentally realized until now only capture a small segment of theoretical proposals [206, 207], and it is expected that synthetic SOCs that have no counterpart in solid state systems will be demonstrated in ultracold atomic systems. This will open a way to probe novel physics with unexpected properties.

9.2 Spin-orbit coupling in fermionic systems and FFLO

For fermions, spin-orbit coupling provides many interesting possibilities. By changing the structure of Fermi surfaces, SOC can reduce the effective dimensionality of the system and thus enhance pairing, and for example, make superconductivity possible at higher temperatures. By enabling topologically non-trivial band structures, it could imply topologically non-trivial superconductors as well. For identical fermions, $s$-wave interactions, which usually dominate physics in ultracold quantum gases, vanish. For a mixture with two different fermionic species, $s$-wave interactions are relevant and usually dominate over $p$- and $d$-wave interactions which could be of interest regarding topological states. When $s$-wave interactions are small, attempts to experimentally realize some interesting quantum states in, for example, a $p$-wave interacting fermionic gas are hampered by large losses in the region where $p$-wave interactions would be large [208]. Spin-orbit coupling provides an interesting way around these problems: a system with usual $s$-wave interactions can give rise to $p$-wave interactions in the basis of pseudo-spin eigenstates [209].

In addition to these possibilities, SOC in a Fermi gas can change the topology of the Fermi surfaces and give rise to Lifshitz transitions [210] already for non-interacting systems. For low densities, fermions accumulate to the lower helicity branch and start filling the states around each minima. The Fermi surfaces are then disconnected. For higher densities, the surfaces eventually merge, and for even higher densities, a separate closed Fermi surface appears inside the larger one since fermions start to fill the higher helicity branch as well. We demonstrate this in figure 33 reproduced from the article by Wang et al. [197] who studied spin-orbit interactions in a harmonically trapped Fermi gas. Changes in the Fermi surfaces can affect phase diagrams dramatically since they can, for example, influence nesting properties and densities of states in an important way.

Let us explore how FFLO states in systems with spin-orbit coupling are approached theoretically. Following the notation of Zhang and Yi [211] for a system of two-component fermions interacting via contact interaction with strength $U$, the Hamiltonian is given by

[TABLE]

where for attractive interactions $U<0$. Futhermore, $\xi_{k,\sigma}=\epsilon_{k,\sigma}-\mu$, $\epsilon_{k,\sigma}=\hbar k^{2}/2m_{\sigma}$, and $\mu$ is the chemical potential. If the system is in a lattice, this single particle dispersion in the absence of SOC would be changed to lattice dispersion. The Hamiltonian $\hat{H}_{SOC}$ describes the spin-orbit couplings, and in a 2D system with Rashba coupling, this term would be

[TABLE]

The last term $\hat{H}_{\mathrm{Z}}$ is due to Zeeman interactions and has a form

[TABLE]

Often simplifications are assumed for the Zeeman fields $h_{\beta}$ such that, for example, $h_{y}=0$ and only the out-of-plane Zeeman field $h_{z}$ and in-plane Zeeman field $h_{x}$ contribute.

The Hamiltonian in (9.2) is one typical example, and variations exist. For example, the spin-orbit coupling could be expressed on a different basis so that an equivalent Hamiltonian might superficially look different [212]. Also one might consider SOC that is a mixture of Rashba and Dresselhaus couplings with equal Rashba and Dresselhaus (ERD) spin-orbit couplings [213, 214, 215] being especially relevant to most experiments done so far. In lattices the momentum summations would be restricted to the first Brillouin zone and renormalization of contact interaction is unnecessary due to natural cut-off provided by the lattice.

Dimensionality can also be varied from 1D to 3D and sometimes SOC might be 2D even if the system itself is 3D [216]. The type of SOC is essential for symmetry breaking properties of the potential FFLO state: some forms of SOC explicitly break the rotational symmetry, for instance, which therefore cannot be broken spontaneously. However even in the presence of SOC, the density is typically uniform in the normal state, and the FFLO state can break the translational symmetry. The existence of a preferred momentum as given by SOC, may also influence the collective mode spectrum of the FFLO state.

Without interactions, the Hamiltonian (9.2) could be diagonalized into a helicity basis straight away, but interactions necessitate additional approximations. There are two dominant approaches. In the first approach, one assumes a Fulde-Ferrell type order parameter (note: we mostly follow the notation from reference [211] and therefore notation here deviates slightly from (14)

[TABLE]

and then derives a usual BCS type mean-field theory around it. This implies a mean-field Hamiltonian [211]

[TABLE]

where

[TABLE]

and $\mathcal{I}$ is a 4x4 unit matrix. Here a basis of $(a_{{\bf k},\uparrow},a_{{\bf Q-k},\uparrow}^{\dagger},a_{{\bf Q-k},\downarrow},a_{{\bf k},\downarrow}^{\dagger})^{T}$ was used and $\Lambda_{\bf k}=\alpha(k_{x}+ik_{y})-h_{x}$. From this mean-field Hamiltonian, the grand potential

[TABLE]

can be computed and then minimized to find the ground state for a set of parameters, see section 2.7.

In the second related but more general and numerically more demanding approach, one allows arbitrary spatial variation of the order parameter $\Delta({\bf r})$ and then solves the BdG equations for the quasiparticle amplitudes of the mean-field theory. Since we have already explained the BdG equations in section 8, here we simply guide the reader to references [217, 218, 219, 113] for descriptions of the SOC case.

In the absence of SOC, limitation to the FF-ansatz can sometimes be too restrictive since LO states typically have lower energy. However in the presence of SOC, the assumption of the FF-ansatz is often better justified. As explained previously, SOC can create ground state degeneracy while (in-plane) Zeeman fields break this degeneracy, and one minimum away from ${\bf k}=0$ becomes energetically favorable. This makes it more natural to assume a single plane wave ansatz for the order parameter as well. This has been studied by Xu et al. [113] who explored the competition of LO and FF states in a lattice in 2D optical lattices and found that FF states dominate once relatively large spin-orbit coupling is present. We demonstrate this with figure 34.

For polarized Fermi gases without SOC, it is known that population imbalance tends to drive the system towards phase separation. When such considerations were extended to systems with SOC, it was found by Iskin and Subaşi [220] and by Yi and Guo [221] in continuum that SOC can counteract the tendency towards phase separation and stabilize uniform superfluid phases. Furthermore in the presence of SOC, nodes of the quasiparticle dispersions might behave in a complex way. However these studies did not consider FFLO phases. A later article by Iskin [217] solved BdG equations in the presence of Rashba SOC in a 2D system in a continuum and found that SOC implies polarized superfluid states replacing superfluids with oscillating order parameters. However the chosen polar coordinate system imposed rotational symmetry on the solutions, and for example, FF type solutions were then excluded by construction.

As explained in earlier sections, usually in the absence of SOC one can conclude that FFLO states occupy only a narrow region of the parameter space in free space; although a considerable one in lattices. When researchers started to explore FFLO in systems with SOC, they noticed that this conclusion was no longer valid. Zheng et al. [216] studied 3D systems with Rashba SOC and an in-plane Zeeman field in continuum. They assumed an FF ansatz for the order parameter and found a large region of the parameter space where the FFLO state might appear. Since the system was 3D, this phase was expected to have long-range superfluid order as well. Zheng et al. explained the sudden dominance of the FFLO state by observing that the Zeeman field with SOC distort the Fermi surfaces and shift them in such a way that they are no longer centered around the origin. This suppresses the usual BCS state and gives rise to FFLO states.

Hu and Liu studied a closely related system in reference [212] and outlined a finite temperature phase diagram at a broad Feshbach resonance. They predicted that for sufficiently large SOC, the superfluid state is always an FF state. They also estimated a promisingly high critical temperature of $T_{c}\sim 0.2\,T_{F}$ for the transition, where $T_{F}$ is the Fermi temperature. In reference [222], Liu and Hu extended their model to also include an out-of-plane Zeeman field and again found large regions of the parameter space occupied by the FFLO state. We show their predicted phase diagram in figure 35. They again predicted that the critical temperature could be a substantial fraction of the Fermi temperature.

The instability of the normal state to FFLO with Rashba SOC was then elaborated by Liu [223], who used the Thouless criterion and beyond mean-field corrections to compute (among other things) the critical temperature. With Rashba SOC coupling $\alpha k_{\mathrm{F}}/E_{\mathrm{F}}=1$ (where $E_{\mathrm{F}}$ is the Fermi energy), in-plane Zeeman field $h_{x}=0.5E_{F}$, and in the unitarity limit the critical temperature for the FF superfluid was found to be quite high, $T_{c}=0.2\,T_{F}$.

Around the same time, similar conclusions were found by Dong et al. [224] and Zhou et al. [225]. In these studies, the model was somewhat different in that they focused on the symmetric 3D spin-orbit coupling in continuum $v_{x}k_{x}\hat{\sigma}_{x}+v_{y}k_{y}\hat{\sigma}_{y}+v_{z}k_{z}\hat{\sigma}_{z}$. Despite this difference, by assuming FF type order parameter Dong et al. [224] found that the FF state dominated the phase diagram both at zero as well as at non-zero temperatures. In addition to solving the mean-field phase diagram, Dong et al. also solved the two-body problem for the dimer bound state and found dimers with finite momentum. The Zhou et al. [225] study complements these results by exploring the instability towards the FFLO state using a small $Q$ expansion of the grand potential. Furthermore it classifies FFLO states according to the presence of a gap and structure of the nodal surfaces in momentum space.

In an experiment, a trapping potential would be present. A harmonically trapped system of two-component fermions in an optical lattice with SOC was studied by Iskin [218]. In this study, BdG equations were solved in a 2D lattice model for $150$ fermions. Iskin found that with Rashba SOC and in-plane Zeeman field, the order parameter was of the FF type, without a nodal structure in the order parameter expected from the LO state. We show an example of this in figure 36.

9.3 Topological states in systems with spin-orbit coupling

Topological states of matter have become an intense research area in recent years [226, 196]. In the context of spin-orbit coupling, Tewari et al. [227] pointed out that SOC with Zeeman fields in a fermionic system can be used to drive a topological phase transition into a topologically non-trivial superfluid. This happens when the Zeeman splitting becomes so large that the chemical potential is in the gap (opened by the Zeeman splitting) of the different dispersion branches. The question then arose whether or not the FFLO phases that appear in systems with spin-orbit couplings might also have non-trivial topological properties. This was answered affirmatively in several recent studies, that were published close to each other in time [211, 228, 229, 230]. The key feature making this possible was the presence of Zeeman fields in addition to spin-orbit coupling.

Let us first clarify what is meant by the term “topological” here. Earlier we pointed out that in the presence of the spin-orbit coupling, the topology of the Fermi surface can change for non-interacting particles. Similar arguments can be extended to the intertacting case by considering the locations on momentum space where quasi-particle dispersions change signs. Such an approach has been used to classify different phases as being topologically distinct [220, 221, 217, 214]. Historically transitions between phases with different Fermi surface topologies have sometimes been called topological transitions, but due to the present-day use of the term topological, it is better to refer to them as Lifshitz transitions.

Changes in the topological structure of the nodes of the dispersions should not be confused with classifying states according to their topological properties by computing their topological invariants, such as the Chern number. For each k-vector, diagonalization of the BdG equation gives rise to four eigenenergies, and we can label these eigenenergies in terms of quasi-particle and quasi-hole index $\eta=\{+,-\}$ and helicity index $\nu=\{1,2\}$ as $E_{{\bf k},\nu}^{\eta}$. If we then denote the associated wavefunctions with the shorthand notation $|n\rangle$ $n=(\eta,\nu)$, we can calculate the Berry curvature, see for instance [211]

[TABLE]

We then get the Chern number for each helicity band of quasi-particles and holes by integrating over the Berry curvature over the Brillouin zone

[TABLE]

The Chern number for the superfuid is then the sum of the Chern numbers of the quasi-hole (i.e. $\eta=-$) bands $C=\gamma_{1}^{-}+\gamma_{2}^{-}$. For topologically trivial phases, the Chern number vanishes when it is non-zero for topological phases.

Note that here we compute the Chern number for the interacting state. However one could also compute the Chern number for single particle bands and find that some bands are topologically non-trivial. This adds an additional layer of potential confusion, since having a topologically non-trivial single particle band does not imply that the realized state would also be topologically non-trivial. However it may lead to other interesting consequences, such as superfluidity in a flat band [231, 232, 233, 234, 235].

Zhang et al. [211] and Qu et al. [228] studied a 2D Fermi gas with Rashba SOC. They identified a gapped FF state with a center of mass momentum along the x-axis as well as a gapless nodal FF state (with two disconnected gapless contours in momentum space), and finally a topological FF state. A topologically trivial gapped FF state appears for large SOC or, alternatively a small out-of-plane Zeeman field. At the phase boundary between nFF and gapped FF state, quasi-hole and quasiparticle spectra touch the Fermi surface at different points. This is in contrast with the boundary between topological FF state and the gapped FF state where the gaps close at the same point. By solving the BdG equations of a tight-binding model in a 2D strip geometry, Qu et al. [228] also confirmed the existence of chiral edge states of topological FF states using a square lattice model.

Cao et al. [229] also found the possibility of gapped and gapless topological FF superfluid in a 2D Fermi gas with SOC in continuum. For an example of their full phase diagram, see the figure 37. By adding a disorder potential along the y-axis, they confirmed the robustness of the Majorana edge modes. In conventional systems without SOC, transverse superfluid stiffness of the FF states vanishes [81, 236]. Interestingly Cao et al. [229] found a positive definite superfluid density tensor when SOC was present. Furthermore they also extended the study beyond Rashba coupling by considering SOC of the type

[TABLE]

and varying the angle parameter $\psi$. If $\psi=\pi/4$, SOC is of Rashba type while $\psi=\pi/2$ corresponds to the experimentally realized equal Rasha and Dresselhaus SOC. Importantly Cao et al. confirmed that while topological states disappear in the limit of equal Rashba and Dresselhaus (ERD) coupling, the existence of topological states does not require pure Rashba type SOC.

Closely related to these studies is the paper by Xu et al. [230] who also identified the possibility of realizing (massless chiral) Weyl fermions in these systems. Dirac fermions are massive and described by four (or more) components, but in the zero mass limit, the Dirac equation allows solutions with two components known as Weyl fermions. By increasing the out-of-plane Zeeman field $h_{z}$, one can move from a gapped FF phase to a topological FF phase which is gapped except at Weyl nodes around which the dispersions are linear and can be described by the chiral Weyl equation. For relatively large values of $h_{x}$, a phase where excitations are gapless also outside Weyl nodes was found and designated as a gapless topological FF state. A quasiparticle excitation gap closes at the transition and identifies the points where Weyl fermions might appear. Xu et al. [230] noted that properties of these Weyl fermions could be controlled with Zeeman fields and suggested they could be observed by studying the changes in the speed of sound around the Weyl nodes.

SOC with the experimentally realized equal Rashba and Dresselhaus SOC was studied by Wu et al. [214, 215] in 2D systems and by Iskin and Subaşi [213] in a 3D system. Although Cao et al. [229] did not find topologically non-trivial states (meaning non-zero Chern number) for ERD, Wu et al. [214, 215] pointed out highly non-trivial nodal structures in momentum space that could be used to further classify, for example, FFLO states. Similar conclusions were found by Iskin and Subaşi [213]. Since SOC can generate effective $p$-wave interactions, the symmetry of the realized paired states is of interest. This question was studied by Chan and Gong [237] who found, for example, triplet pairing in the limit with vanishing in-plane Zeeman fields.

Dimensionality influences the appearance of topologically non-trivial states. Liu and Hu [238] solved BdG equations in a 1D system. In the absence of SOC, one expects either a BCS superfluid or an FFLO superfluid (or their 1D equivalents without long range order), but SOC with Zeeman field could drive the system into a superfluid with zero energy Majorana fermions [195] at the edge. This is called a topological superfluid since a zero energy edge mode in real space is taken as an indicator of an interface with two topologically distinct regions with the surrounding vacuum being obviously topologically trivial. Physically the appearance of Majorana fermions can be understood from the analogy with chiral $p$-wave superfluids. As explained earlier, SOC mixes states in such a way as to create effective $p$-wave interactions in a system with underlying $s$-wave interactions. A Zeeman field can then open a gap between two branches so that, for sufficiently large Zeeman fields, all atoms occupy the lowest branch, thus creating a possibility for $p$-wave pairing.

Liu and Hu [238] found that the FFLO state appears for small SOC while topological superfluid appears at large values of SOC. Chen [239] also studied the mean-field theory with different SOCs including ERD, using a 1D lattice model. In the system of [239], instead of the earlier mentioned Chern number, one could compute the $\mathbb{Z}_{2}$ number which is $1$ for topologically trivial systems and $-1$ for topologically non-trivial ones. With ERD, Chen found the possibility of topologically non-trivial FFLO state with $\mathbb{Z}_{2}=-1$ by tuning the band filling. These states could be further classified based on whether or not they were also gapped. Chen also solved the BdG equations in a harmonic trap and intepreted that the topologically non-trivial state would manifest itself there as zero energy Majorana fermions at the edge of the system.

Dimensional crossover and the competition between FFLO and Majorana states was studied further by Seo et al. [219] by tuning the dispersions from multi-dimensional to 1D using a lattice potential in the $xy$-plane. In agreement with Li and Hu  [238], they found FFLO states for smaller values of the (Rashba) SOC and (topological) uniform superfluid state with Majorana fermions at the edge for higher values of SOC. It was also found that lowering the dimensionality from 3D to 1D enables a crossover from the 3D uniform superfluid phase into a FFLO phase or to the topological uniform superfuid state. In a bilayer system with SOC, FFLO can appear without any population or chemical potential imbalance, and can transit to topological superconductivity as the interaction increases [240].

Qu et al. [241] used a 1D lattice model with SOC and a large out-of-plane Zeeman field to study the competition between FFLO superfluids and superfluids with Majorana fermions. They computed the phase diagram for the model and, among other things, found that with increasing strength of SOC superfluids with Majorana fermions occupy ever larger areas of the parameter space. This result is consistent with references [238, 219] which did not consider a lattice model.

As dimensionality is reduced, fluctuations become more pronounced and might destroy the ordered states predicted by simple mean-field theories. In particular in a 2D superfluid, Berezinsky - Kosterlitz - Thouless (BKT) transition is expected to appear. As temperature is lowered, exponentially decreasing correlations of the disordered state above the critical temperature are replaced with just algebraicly decaying correlations below the BKT transition temperature. This phase with quasi long-range order is a superfluid, and there is a universal jump in the superfluid density at the BKT critical temperature. BKT physics has recently been studied in the context of spin-orbit coupled Fermi gases for FF states [242, 243]. In the absence of SOC, the transverse superfluid density for the FF state vanishes, and zero critical temperature is expected in that case [236, 244]. When SOC with in- and out-of-plane Zeeman fields are introduced, Xu et al. and Cao et al. found non-zero BKT transition temperature for both gapped as well as gapless FF superfluids since transverse superfluid density could be non-zero. This suggests that observations of FFLO phases in 2D systems with spin-orbit couplings might be possible at non-zero temperatures.

All articles on SOC and FFLO that we have considered so-far rely on mean-field descriptions. An interesting deviation from this pattern is the work by Zvyagin and Schlottmann [245] who considered an exactly solvable 1D fermionic model with spin-orbit interactions and strongly attractive onsite interactions. Using Bethe ansatz they obtained the critical exponents for superfluidity and density waves and found that the exponent for the superfluid was the smallest. This suggested an instability to a FFLO phase in the presence of weak interchain coupling.

9.4 Summary of the status of FFLO in the presence of spin-orbit coupling

As is clear the literature on FFLO states in the presence of SOC has expanded rapidly. All studies so far have been mean-field theories and mostly focus on Rashba or equal Rashba and Dresselhaus spin-orbit couplings. Studies have revealed that compared to systems without SOC, FFLO states become favorable in the presence of SOC. Furthermore, in contrast with systems without SOC, FF states are expected to dominate over the LO ones. In 1D systems, it has also been pointed out that states with Majorana fermions at the edge of the system compete with the FFLO states. As SOC strength increases, eventually FFLO gives way to Majorana states. In higher dimensions, Majorana states can also appear if there are vortices in superfluids with spin-orbit couplings [246, 247, 248]. But the most remarkable insight in higher dimensions is the prediction of topological FFLO states in systems with Rashba SOC and both in- and out-of-plane Zeeman fields. In such topologically non-trivial states, it is possible that the gap closes even in the bulk.

It will be interesting to see how accurate the predictions of mean-field theories are. Beyond mean-field effects might very well lead to pronounced changes in, for example, lower dimensional systems as well as in optical lattices with relatively large filling fractions. Orbital physics in optical lattices can also enable new possibilities. A recent proposal to realize FFLO states in unpolarized Fermi gas by Zheng et al. [249] relied on moving optical lattices which couples the $s$- and $p$-orbital bands with each other. This gives rise to a single particle Hamiltonian equivalent to the presence of spin-orbit coupling, although in a lattice model, the off-diagonals also contain terms that are higher than lowest order in $k$.

10 Detection of FFLO state

The challenge in realizing the FFLO state is accompanied by the difficulty of actually detecting it. The FFLO state is predicted to occupy only a rather small part of the homogeneous phase diagram in 3D. Thus, in spatially inhomogeneous trapped gases, the state will be present in only a small part of the gas, except possibly in 1D and spin-orbit coupled systems where the FFLO state involves a larger parameter range. The effect of trapping potential has been considered in section 8. In particular, the size of the FFLO phase domain in the cloud should be larger than the FFLO wavelength $1/q_{\mathrm{FF}}$. Atomic gas experiments with more box-like trapping potentials [76, 77, 78, 79] are believed to improve the situation [250], as the size of the FFLO domain in the cloud can then be increased and the experimental signature of the state improved.

Here we provide a comprehensive review of the various schemes that have been proposed for detection of the FFLO state. Indeed, most of the theoretical works on the topic propose some way of probing the state, although only a few analyze the problem of detection in detail. We roughly divide them into six different categories: ramps and quenches, density profiles, interferometric methods, collective modes, noise correlations, and spectroscopies. Many of the methods have been proposed and analyzed in 1D systems, but we will consider them all on equal footing since many of the methods are applicable in any dimension. Likewise we consider methods proposed both for continuum and lattice cases, and in spin-orbit coupled systems; many of the methods are applicable for all systems.

10.1 Ramps and quenches

The now standard methods for detecting balanced BCS-type pairing in ultracold gases were already proposed early on also for detecting FFLO correlations [252]. One method is an interaction strength quench, in which an external magnetic field is swept across the BCS-BEC crossover. This process projects the Cooper pairs in the BCS side onto bound molecules in the BEC side, and the momentum distribution of the molecules is then probed by releasing the gas from the trap and following the free expansion of the molecular cloud [253] or by using momentum-resolved Raman spectroscopy [254]. Finding a molecular condensate on the BEC side is a signature of an initial condensate of Cooper pairs on the BCS side. In the present context, since the Cooper pairs in the FFLO state will have a finite center-of-mass momentum, the momentum distribution of the corresponding molecules would show the same center-of-mass motion [72, 255, 256, 124, 257, 258, 259, 260, 261, 125, 262]. While the projection to Feshbach molecules can be used for probing the FFLO state in lattices also [259], the presence of an optical lattice allows for other schemes. Indeed quenches and the expansions can be done in various ways, but the common denominator when considering detection of the FFLO state is that the pairs in the FFLO state are projected onto some more stable configuration (bound molecules, repulsively bound pairs, doublons) and the momentum distribution of these is then probed, see Fig. 38. Futhermore an interesting variant was proposed in [263] in which triplet correlations between atoms residing in the nodes of the oscillating LO-order parameter are mapped onto $p$-wave molecules. In 3D systems, the pair formation dynamics in the unitary regime were found to be very rapid [264]. However it is not clear whether this imposes constraints for seeing FFLO correlations in molecular momentum distributions, since formation of molecules is actually a desired property. However molecules are strongly interacting in the unitary regime, and a quantitative understanding of how FFLO correlations survive the BCS-BEC ramp in 3D systems is missing.

The stability of the correlations has been studied in 1D systems. Ref. [251] considered a 1D system in which interactions are quenched from the attractive to the repulsive side, and the initial FFLO correlations are imprinted on repulsively bound pairs. However signatures of FFLO correlations were found to vanish rapidly in the interaction strength sweep, showing the importance of a sufficiently rapid sweep. Similar decay of the FFLO correlations in the pair momentum distribution were found also in Refs. [266, 267, 26]. However also interestingly the work suggests that features in the spin-density are more robust and offer a possible way for detecting the FFLO state in the expansion of a 1D gas. Since the loss of FFLO correlations was caused by the unavoidable collisions when expanding along 1D lattice, the problem may be less prominent in higher dimensions. On the other hand, Ref. [265] considered similar 1D setting, and the separation of the doublons and single atoms was beneficial, as the expansion velocity of unpaired atoms was found to be connected with the FFLO momentum $q$. In this case, the atoms remained interacting during the expansion, that is, only the trap potential was quenched, see Fig. 39. Consequently the FFLO momentum $q$ was directly reflected in the expansion velocities of the unpaired fermions because in 1D they form a Fermi sea of non-interacting quasiparticles with Fermi momentum $k_{\mathrm{F}\uparrow}-k_{\mathrm{F}\downarrow}$. The use of time-of-flight expansion for measuring momentum distributions and detecting the FFLO state have also been considered in spin-orbit coupled systems [213, 216, 268].

Alternatively one could quench to a non-interacting regime, in which case the release of the gas will allow the measurement of the momentum distributions of the two atomic components separately. While this will not yield direct access to pair-distribution function, many works suggest that FFLO state will have observable effects in single-particle momentum distributions as well [69, 67, 64, 124, 110].

10.2 Density profiles

The FFLO state, particularly the LO state, has also been shown to influence the density distributions of the two atomic components [62, 256, 269, 270, 110, 271, 272, 261, 111, 156, 66]. While the total density profile is largely unaffected by the FFLO state, the local polarization appears to be a viable option for the detection [110, 66]. The state might thus be observable in species selective absorption imaging [184], or in phase contrast images [273], that directly probe the density difference between the two species. Unfortunately in 3D, even if the FFLO phase domain would be large enough to exhibit actual density modulations, any such off-centered density modulations are easily lost in the columnar absorption image, which probes the integrated density through the cloud (along the direction of the imaging laser).

Since the FFLO momentum tends to align itself in the longitudinal direction in elongated traps, one would expect density modulations, and FFLO correlations, to be most easily observed in quasi 1D-systems [162]. Also, quantum gas microscope setups for fermions [27, 28, 29, 30, 31, 32] with 2D lattices can avoid the problems with columnar integration altogether. At sufficiently high densities, 2D systems are predicted to exhibit a special angular FFLO state [166] in which the modulations are positioned in an angular fashion instead of radially. Still the density variations predicted even in 1D systems are quite minor, see for example Fig. 40 and Refs. [274, 153]. Minor modulations in the densities may be better observed through spectroscopic means that specifically couple to density modulations of a given periodicity, such as Bragg spectroscopy [136, 65, 275].

A quite different approach would be to utilize the connection between the trapping potential, local density and pressure for determination of the equation of state [276, 277, 278]. Equation of state, or particularly the heat capacity [279], is influenced by the FFLO state, and one would expect jumps in the specific heat at various phase transitions [72, 280, 281], although precise thermometry of the gas is still one of the outstanding experimental problems. From local densities one can calculate also local entanglement entropy, which was studied in Ref. [282] and shown to be a viable candidate for detecting the normal state to FFLO state transition.

10.3 Interferometric methods

Interference patterns between two atomic clouds with FFLO-correlations have also been considered by a number of groups [284, 274, 283]. Fig. 41 shows one such setup, in which two atomic clouds are allowed to expand, and upon overlapping an interference pattern is seen. The appearing interference fringes probe the phase difference between the sources, and in the case of an FFLO state with an oscillating order parameter, the fringes become staggered [283]. To be able to probe the interference between pairing fields in the first place, the method needs to be combined with the magnetic field sweep into the BEC side, in a similar fashion as considered in section 10.1.

An interesting alternative is the suggestion of using a Josephson current for detecting an FFLO state [285, 113]. The Josephson current is proportional to the local order parameter, and a local probe can thus identify position dependent order parameter fields. This can be realized with two clouds coupled via a small dimple potential channel [285]. Since the scheme couples directly to fermionic pair correlations, no BCS-BEC ramps are needed.

For interference studies, one needs a scheme for producing two (or more) clouds. One possibility is to prepare these in a double well potential or in an optical lattice. Alternatively a single atomic cloud could be split in two by lasers, or Bragg beams could be used for extracting parts of the cloud for the interference protocol [286, 274].

10.4 Collective modes

Superfluidity is a collective phenomenon, and it has profound effects on collective excitations. A natural question is then whether collective modes or excitations can be used for identifying the superfluid FFLO state.

Most studies of the collective modes in FFLO state consider uniform lattice or continuum systems and apply mean-field theories with the random-phase approximation [136, 137, 287, 131, 143, 144, 145], although 1D systems have also been studied with exact numerical methods [288, 289]. The breaking of the rotational symmetry in the FFLO state has been shown to result in anisotropic speed of sound [131]. This is discussed in detail in section 6. Probing the anisotropic speed of sound has also been considered in the context of spin-orbit coupled systems as a possible way of detecting the FFLO state [230, 242, 249]. Furthermore as one might expect, the FFLO state influences the nature of vortices in the superfluid state [109]. These may offer a way for detecting the FFLO state [290, 291], although the low visibility of vortices in the BCS and unitary regimes can prove challenging [292].

However more research would be needed for understanding the effect of FFLO state in other collective modes, such as the often used trap modes (breathing modes, dipole modes, etc.), and the interplay between the second sound and unpaired atoms in the FFLO state would deserve more analysis.

10.5 Noise correlations

Simple density probes suffer from rather weak signals (small density perturbations) and strong noise (large background density). Therefore higher order correlations, such as density-density and current-current correlations, have been considered as potential probes for the FFLO state. The use of shot noise correlations, in which the atom gas in the FFLO state is first quenched to a weakly interacting state and allowed to expand, and the correlated appearance of fermions at momentum states $-k$ and $k+q$ for variable $k$ but fixed $q$ was already proposed by Yang [252].

The use of noise correlations as a probing scheme in general was reviewed in Refs. [46, 294], but here we briefly review the ideas that have been presented for the detection of the FFLO state itself. Ref. [293] shows a detailed analysis of the various density and spin correlations in a 1D lattice and their applicability for detecting the FFLO state. While practically all higher order correlators show signatures of the FFLO state [274, 64], the $\langle n_{\uparrow}n_{\downarrow}\rangle$ correlator appears to provide particularly clear signatures, with peaks separated by the FFLO momentum, shown in Fig. 42. These results seem to also hold in higher dimensions [295, 296]. The use of noise correlations in spin-orbit coupled systems has also been proposed [230, 268, 249, 26].

In practice, the noise correlations can be studied either in-situ or in expanding clouds. In a freely expanding cloud, the observed spatial correlations yield information regarding the initial momentum correlations [297, 68]. Such time-of-flight probes need to be combined with interaction quenches to molecular states [124] or noninteracting states. Still one would expect the expansion dynamics to be problematic for the conservation of the correlations. Alternatively, the various density-density correlations can also be studied in-situ. A spatially resolved quantum polarization spectroscopy protocol was suggested [298, 299] as a way to probe spin-spin correlations, and density-density correlations might be used to probe using quantum gas microscopes [27, 28, 29, 30, 31]. Furthermore although the ground state FFLO state does not have a total mass current (see section 2.9), it does contain mutually cancelling superfluid and quasiparticle currents in the ground state. Thus current-current correlations, if one could come up with an idea on how to probe them, could provide a good signature of the FFLO state.

10.6 Spectroscopies

Since the FFLO state involves a spatially varying order parameter, an obvious way for detecting the state would be to use probes that have been used for probing the order parameter and excitation gap $\Delta$ in spin-balanced systems [300]. The spectral function can be probed by numerous methods [301, 302]: radio-frequency spectroscopy, Raman spectroscopy, Bragg spectroscopy, and lattice modulation spectroscopy, to name just a few commonly used methods. Radio-frequency spectroscopy is already a standard method for probing the single-particle excitation spectrum, although possible interactions in the final state make interpretation of results complicated. Bragg spectroscopy can be used for probing elementary particle-hole excitations, but in the case of long-wavelength Bragg spectroscopy, also collective excitations. A particularly appealing aspect of Bragg spectroscopy is the periodic coupling to atomic density fields, and thus a possibility of probing the weak density modulations in the FFLO state [136, 65, 275]. Lattice modulation spectroscopy for observing the FFLO state was considered in Ref. [303]. In this method, a periodic variation of the depth of the lattice potential excites the system, creating doublon excitations. The width of the double occupancy spectrum was shown to provide a measure of the FFLO momentum. In similar fashion, other spectroscopies probe different parts of the system’s spectral function, see for example Fig. 43.

Inhomogeneous order parameter results in Andreev mid-gap states [304, 110] and unpaired atoms yield gapless excitations [305, 155]. These provide single-particle excitations at energies below the superfluid gap or even at negative energies [306]. However a major problem in deducing the presence of the FFLO state from spectroscopic probes is that the inhomogeneity due to the trapping potential already results in low-energy excitations in the spectral function [180], as discussed in section 8. To avoid these problems, local probes such as methods utilizing lasers (Bragg spectroscopy, Raman spectroscopy) may be most useful, and they have indeed been used for probing the essentially homogeneous order parameter $\Delta$ [307]. Also methods acting on the whole cloud, such as radio-frequency spectroscopy, can be used as long as the response is considered only in a limited region of the gas, such as in the experiments done in JILA [308]. Again box traps will also help make the spectroscopic signature clearer.

An interesting addition to usual spectroscopies is the use of impurities. Spin-selective scattering potentials [309] and impurities [310, 311, 312] are influenced by the local pairing gap, allowing a very precise spatially resolved spectroscopy of the FFLO order parameter.

Use of spatially- and momentum-resolved spectroscopies have also been considered in spin-orbit coupled systems [228, 211, 225, 230, 238, 215, 230], and were studied in more detail in [222]. In spin-orbit coupled systems, the connection between spin-imbalance and the FFLO state is not clear, since the spin-orbit coupling itself is sufficient for producing a finite FFLO momentum. Spin-orbit coupled systems are analyzed in more detail in section 9. Fig. 44 shows momentum-resolved radio-frequency spectrum of a spin-orbit coupled system. A particularly interesting proposal is to probe the edge modes in a topological FFLO phase using spatially resolved spectroscopy [239]. In the case of both spatially- and momentum-resolved spectroscopy, the topological FFLO state is expected to manifest as gapless excitation at finite momentum $k$ [237].

11 Conclusions and outlook

In this review, we have presented the theory of the FFLO state in lattice and trap geometries, and reviewed literature on theory and simulations of the lattice FFLO state in the context of ultracold quantum gases (UQG). A review of recent literature on the FFLO state in the presence of spin-orbit coupling was also given, as well as a comprehensive description of possible methods to observe this elusive state of matter in quantum gas experiments.

As section 10 shows, there are plenty of options for experimental probes that may reveal salient features of the FFLO state. Promising recent developments that may be crucial for unambiguous observation of the FFLO state are the quantum gas microscopes [313, 52, 314], in particular their fermionic versions [27, 28, 29, 30, 31, 32]. In those systems, one has unprecedented access to small spatial variations in densities as well as spin and density correlations that all carry specific signatures of the FFLO state. The development has been similarly impressive in realizing spin-orbit coupled systems in UQG [40, 41], and it will be interesting to extend such studies to the case of imbalanced fermion systems. One more promising development is the realization of box potentials and homogeneous gases [76, 78, 105, 79, 315, 39] instead of the usual harmonically trapped gases. This can make possible FFLO signatures more clear, counteract the tendencies to phase separation, and realize the fixed particle number situation instead of fixed chemical potential, which may lead even to the stability of the Sarma state [13, 69]. In the future, there will probably also be many artificial solid state systems where the orbital magnetism and Pauli paramagnetic effects can be tailored in a way amenable to the FFLO state. In general, examples of such designer matter range from microscopic superconducting structures realizing artificial atoms [316] to nanostructured [317] and thin-film materials [318] to atomic scale artificial lattices [319].

The FFLO state in lattices poses some major open questions to be approached by future theory work, furthermore, several conceptually new directions of research can be envisioned. To connect with near future UQG experiments, further work on finite temperature properties is needed. Interestingly, some beyond-mean-field calculations [154, 157] find a finite momentum signature in the pairing structure factor well above the critical temperature for superfluidity $T_{c}$. It would be important to study in more detail whether and when such precursors of the FFLO state exist, since they may be observable in the UQG experiments even at presently available temperatures. This relates to the fundamental but almost unstudied question of the nature of the normal state right above the critical temperature systems susceptible to FFLO correlations, especially for strong interactions. One can also search for precursors of the FFLO state at the limits of imbalance: polaron systems in which there are only individual atoms of the minority spin state, and the case of vanishing polarization where the FFLO state is connected with pinned solitons and soliton trains that can arise either in the ground or excited state [305].

Another key question is the stability of the FFLO state at temperatures finite but well below $T_{c}$. In continuum systems, the stability of the FFLO state towards fluctuations of the superfluid phase is a subtle issue. Already the original works by Fulde and Ferrell [2], and Larkin and Ovchinnikov [4], point out the fundamental difference of the supercurrents in the directions parallel and transverse to the FFLO wave vector $\mathbf{q}$. The question of stability was studied in [320] where the Ginzburg-Landau theory in the context of type II superconductors predicted that any spatial dependence of the order parameter that varies in one dimension only (such as the (double)planewave FF(LO) ansatz) is unstable, however if the order parameter structure is essentially 2D such as triangular, stability is established. A more elaborate analysis based on mean-field approximation with quantum fluctuations, and taking into account anisotroptic superfluid density [236], showed that the FF state in a 2D continuum system is unstable due to vanishing superfluid density in the direction transverse to the FF wavevector. However the fate of the LO state remained an open question since the microscopic BCS-based approach does not easily accommodate the LO ansatz (c.f. the discussion in section 2.4). Radzihovsky and coworkers have extensively studied the stability of the continuum FFLO state and its collective modes in various dimensions using a Ginzburg-Landau type approach [109, 19, 81, 321]. They showed that the spontaneously broken rotational symmetry leads to a soft Goldstone mode, which limits the superfluid density in the transverse (to the FFLO momentum) direction. Indeed in 3D, the FF state was argued to have vanishing transverse superfluid density, and even in the LO state it was suppressed. Still, FF and LO states were argued to be energetically stable at zero temperature for dimensions $d>1$.

All the above theory knowledge on FFLO stability concerns continuum systems: much less is known about the effect of fluctuations in the lattice case. The Mermin-Wagner theorem, denying true long-range order in 1D and 2D systems is valid there as well, as is the possibility of the Berezinskii-Kosterlitz-Thouless (BKT) mechanism to produce quasi-long-range order and superfluidity in 2D. However due to the breaking of the rotational symmetry and continuous translation invariance as well as nesting effects in lattices, the role of fluctuations may differ from the continuum case. The existing beyond-mean-field studies of the density-imbalanced Fermi gases in lattices indicate stability of the FFLO state against fluctuations: a stable FFLO state was found by DMFT calculations that take into account local fluctuations [164, 155], and by cluster DMFT which also includes non-local fluctuations [156], as well as by QMC methods  [114]. Furthermore, recent cluster DMFT [95] and other numerical [96] studies found the particle-hole equivalent of the FFLO state [61], namely the striped phase, to persist in the doped repulsive Hubbard model. There is thus mounting theoretical and computational evidence for the general existence of the FFLO state in lattices, but which form of the order parameter is the stable one, and at exactly what temperatures, can be still regarded as essentially open questions. This calls for the evaluation of the superfluid density, even in the mean-field case, for estimating the BKT temperature of superfluidity. Furthermore, beyond-mean-field studies of unexplored parameter regimes and with less approximative treatment of fluctuations are needed.

Understanding FFLO superfluidity beyond the mere ground-state phase diagrams also demands studies of collective modes. In spin-orbit coupled systems, it would be interesting to study which symmetries are spontaneously broken by trying to identify the corresponding Goldstone (Anderson-Bogoliubov) modes. For instance the equal Rashba and Dresselhaus SOC breaks the rotational symmetry already at the Hamiltonian level and the dispersion has minima at finite momenta, which could affect the collective mode spectrum of the FFLO state. Further the Higgs mode related to order parameter amplitude fluctuations, and the Leggett mode appearing in multiband systems, are to our knowledge completely unexplored in the FFLO context. Considering the extended (nearest-neighbor interacting) Hubbard model [322], $p$-wave interactions [323], dipolar interactions [324, 325, 326], or dark-state control of Feshbach resonances [327] can also bring new features to FFLO physics.

Observation of the FFLO state, after half-a-century of theory work and experimental effort, would be a major breakthrough. It would profoundly affect the way we understand both superconductivity and magnetism. Realization of the FFLO state with quantum gases in optical lattices would also be a landmark in determining the phase diagram of the Hubbard model, the paradigm model of strongly correlated electrons [328]. Observation of the FFLO phase would suggest, assuming that the system is particle-hole symmetric (see section 2.11), that in the repulsive Hubbard model the striped phase is stable in the corresponding parameter regimes; this would give information about the possible existence of the d-wave superfluid and thereby the role of the Hubbard model in explaining high temperature superconductivity.

Observation and understanding of the classic FFLO state is a grand goal in itself, however, even more is expected when extending the basic ideas of the FFLO physics in simple square lattices to new contexts. In [329, 330] the idea of mixed geometry pairing was introduced: the $\uparrow$ and $\downarrow$ fermionic species were residing in different lattice geometries, which resulted in stable Sarma/BP states, among other things. Rich new physics is expected from multiband systems, especially those that can host flat, that is, dispersionless bands. In a flat band, interaction effects dominate over kinetic energy, which can lead to predictions of room temperature superconductivity [331, 332], for example. Remarkably even when the dispersion is zero, supercurrent can still flow: it was recently found [231, 232, 233, 234, 235] that superfluidity in flat bands is guaranteed whenever the band has non-zero Berry curvature, and a direct relation between flat band superfluidity and the quantum metric was discoverd. The quantum metric [333] describes the distance between quantum states and is the real part of the quantum geometric tensor whose imaginary part is the Berry curvature. This makes one obviously curious about the fate of the FFLO state in a lattice supporting flat bands. Flat bands can be achieved, for instance, by lattice geometry, as in case of the Lieb and kagome lattices, or by (effective) magnetic fluxes; both approaches have been already experimentally demonstrated in UQG [334, 335, 58, 336]. Another avenue for rich new physics would be to consider topological superconductivity [337, 226] for spin-density imbalance, a topic already touched by some of the SOC FFLO studies. Ultimately one may find a new density-imbalanced superfluid whose excitation spectra, broken symmetries and collective modes, and topological invariants differ so much from the FFLO and Sarma state paradigms that it is essentially a new state of matter.

This work was supported by the Academy of Finland through its Centres of Excellence Programme (Projects No. 284621, 13272490 and 303351) and by the European Research Council (ERC-2013-AdG-340748- CODE). We thank Leo Radzihovsky, Dong-Hee Kim, Aleksi Julku, Pramod Kumar, Long Liang, Marek Tylutki, and Tuomas Vanhala for careful reading of the manuscript and their valuable comments.