10.0pt

# Imbalanced Fermi systems and exotic Superconducting Phases

Madhuparna Karmakar  Department of Physics, Indian Institute of technology,
Madras, Chennai-600036, India.

Abstract

Quantum many body phases, in competition or coexistence has always been a
fascinating area of research in condensed matter physics. One of the most
widely investigated scenario of phase competition is that between
superconductivity and magnetism, which in spite of several decades of research
continues to bring up new surprises till date. While superconductivity and
magnetism are naively considered to be antagonists, there are several
circumstances where these many body phases can be found to be cooperating,
coexisting and even promoting each other. The results of such cooperative
behavior are found to show up as exotic quantum phases.

In this chapter we discuss about such an exotic superconducting phase, which
comes to life when a particular class of superconducting material viz. “Pauli
limited superconductors” are subjected to an imbalance between the population
of the fermionic species undergoing pairing. The consequence of such an
imbalance is the emergence of a non-trivial spatially modulated
superconducting phase, popularly known as the Fulde-Ferrell-Larkin-Ovchinnikov
(FFLO) phase. Characterized by a finite momentum pairing between the fermionic
species the FFLO phase hosts a non zero magnetization, leading to an
unconventional coexistent phase of superconductivity and magnetism.

Starting with an introduction to the Pauli limited superconducting systems, in
this article we will review the solid state and ultra cold atomic gas setups
which hosts the FFLO superconducting phase, along with the relevant
experimental diagnostics and their reported observations. This would be
followed up by the theoretical attempts made to understand the physics of
these novel superconducting systems, and the progress made therein. Finally, a
non perturbative numerical technique would be discussed, which is suitable to
capture the behavior of such many body systems, and edges over the existing
ones in several aspects. Being a vibrant area of strong correlation physics,
the Pauli limited superconductors hold several exciting promises for future
investigations. We will touch upon them towards the end of the chapter.

##  0.1 Introduction

Superconductivity is one of the most fascinating phenomena of quantum many
body systems. Right from its discovery in mercury (Hg) by Kamerlingh Onnes in
1911 [1], superconductivity has continued to intrigue the condensed matter
physics community over the past more than a century. The basic idea that
underlies any superconducting material is the formation of a macroscopic
condensate with long range order, which manifests itself via properties like
perfect diamagnetism, absence of electrical resistivity etc. below a
transition temperature (Tc). It was established by Bardeen, Cooper and
Schrieffer (BCS) in their pioneering work that the microscopic mechanism that
gives rise to the superconducting behavior in a system is the formation of
“Cooper pairs”, wherein, two fermions with equal and opposite momenta (${\bf
K}_{\uparrow}$, ${-\bf K}_{\downarrow}$) pairs up by overcoming the repulsion
between them [2].

The BCS theory was originally devised to capture the physics of “conventional”
superconductors, with low transition temperature such as, lead, mercury etc.,
and it was shown that in these materials the lattice phonons served as the
attractive “glue” to pair up the fermions [1, 2]. While the BCS theory is
still the most celebrated microscopic theory in the history of
superconductivity, the discovery of high temperature superconducting cuprates
brought forth the limitations of the BCS theory [3]. The transition
temperatures reported for the cuprates and subsequently for several other
categories of superconducting materials were well above what could possibly be
achieved by a phonon mediated coupling between the fermions [4, 5, 6, 7]. In
other words, there was some alternate mechanism which served as the attractive
“glue” between the fermions in these “unconventional” superconductors with
high transition temperature.

Inspite of the change in the mechanism of superconductivity, some of the
aspects however remained unaltered between the conventional and unconventional
superconductors. For example the superconductivity was still dictated by the
pairing between the fermions having equal and opposite momenta, i. e. the net
momentum of the Cooper pair was zero. Keeping it consistent with the existing
literature, we would refer to such superconducting pairing as the “zero-
momentum” superconductivity, so as to distinguish them from the “finite-
momentum” superconductivity, which would be introduced in the later sections
of this chapter.

It is obvious that once this pairing between the fermions is broken by some
perturbation, superconductivity is destroyed. One of the most prominent
superconducting pair breaking perturbation is the magnetic field. It is now
well established that an orbital magnetic field when applied to a type-II
superconductor leads to the formation of superconducting vortices beyond a
lower critical field, H$\geq$H${}_{c1}^{orb}$ [8, 9]. As the magnetic field is
increased the number of vortices penetrating the superconductor progressively
increases and gives rise to the Abrikosov vortex lattice [10]. As large number
of vortices crowd the superconducting system the vortex cores begin to overlap
each other at an upper critical field, H${}_{c2}^{orb}=\Phi_{0}/2\pi\xi^{2}$,
where $\Phi_{0}=\pi\hbar c/|e|$ is the flux quantum, and $\xi$ is the
superconducting coherence length. The upper critical field (H${}_{c2}^{orb}$)
thus sets the strength of the orbital magnetic field required to destroy the
superconductivity of a material [1, 8, 9].

Superconductivity can also be destroyed by magnetic field via an alternate
mechanism known as the Pauli paramagnetic pair breaking effect, originating
from the Zeeman splitting of the single electron energy levels. When a Zeeman
field is applied to a normal metal the electrons undergo polarization because
the Zeeman effect splits the Fermi surfaces corresponding to the up
($\uparrow$) and down ($\downarrow$) fermion species. This is known as Pauli
paramagnetism. On the other hand, a superconductor comprising of Cooper pairs
is not readily polarized and one needs to break the pair in order to polarize
the system. This pair breaking takes place when the Pauli paramagnetic energy,
$E_{p}=(1/2)\chi_{n}H^{2}$ equals the superconducting condensation energy
$E_{c}=(1/2)N(0)\Delta^{2}$. Here, $\chi_{n}=(1/2)(g\mu_{B})^{2}N(0)$ is the
spin susceptibility of the normal state, $g$ is the spectroscopic splitting
factor of an electron, $\mu_{B}$ is the Bohr magneton, $\Delta$ is the
superconducting energy gap and $N(0)$ is the density of states at the Fermi
level. The upper critical field at which the Pauli paramagnetic pair breaking
effect takes place in a superconductor is estimated to be
H${}_{c2}^{P}=\sqrt{2}\Delta/g\mu_{B}$ and is known as the Chandrasekhar-
Clogston limit [11, 12].

In general, the upper critical magnetic field of a superconductor is
determined based on both the orbital and Pauli paramagnetic pair breaking
effects. The relative importance of these effects is determined in terms of
the Maki parameter $\alpha$, defined as [9],

$\displaystyle\alpha$ $\displaystyle=$
$\displaystyle\sqrt{2}\frac{H_{c2}^{orb}}{H_{c2}^{P}}\sim\Delta/\epsilon_{F}$
(1)

which is the ratio of the two upper critical magnetic fields at zero
temperature (T=0). Here, $\epsilon_{F}$ is the Fermi energy. For most of the
superconductors $\alpha<1$, and thus the system loses superconductivity via
the overlap of superconducting vortices much before the applied magnetic field
reaches H${}_{c2}^{P}$. There are however situations when $\alpha>1$ and thus
H${}_{c2}^{P}<$ H${}_{c2}^{orb}$, viz, (i) the Fermi energy is strongly
suppressed, (ii) the superconducting system is two dimensional and is
subjected to an in-plane magnetic field, such that the orbital effects can be
neglected, and (iii) the system is charge neutral. The class of
superconductors with $\alpha>1$ and which loses its superconductivity by Pauli
paramagnetic pair breaking effect are known as the Pauli limited
superconductors. Prominent examples include heavy fermion superconductor
CeCoIn5 [13], layered two-dimensional organic superconductor $\kappa$-(BEDT-
TTF)2Cu(NCS)2 [14], iron superconductor KFe2As2 [15] and ultracold Fermi gas
with imbalance in the populations of fermionic species [16].

###  0.1.1 Pauli limited superconductivity

Figure 1: Color online: Schematic phase diagram showing (a) the
Chandrashekhar-Clogston limit of Pauli limited superconductor, (b) spatially
modulated superconducting (FFLO) phase at high magnetic fields.

A Pauli limited superconductor in the absence of any orbital effect is known
to demonstrate non trivial behavior in the in-plane magnetic field-temperature
($H-T$) phase diagram [9, 17], as shown by the schematic in Figure 1(a). In
the vicinity of Tc, dHc2/dT $\propto$ $\sqrt{T_{c}/(T_{c}-T)}$ and the initial
slope of H${}_{c2}^{P}$ is infinite, unlike the orbital limiting case for
which the slope is finite. In the absence of Pauli limiting effect the loss of
superconductivity at H${}_{c2}^{orb}$ is via a second order phase transition.
For a Pauli limited superconductor, however, the order of thermal phase
transition changes below a tricritical temperature, T∗, say. While the
transition is of second order for T$>$T∗, a first order transition is realized
for T$<$T∗. In Figure 1(a) we represent the second and first order phase
transitions by solid and dashed curves, respectively. The second order
transition line (dotted) at T$<$T∗ correspond to a metastable state. The low
temperature high magnetic field regime of this phase diagram thus correspond
to a first order thermal transition regime, and the associated field at T=0 is
H${}_{c2}^{P}$.

This Chandrasekhar-Clogston limit H${}_{c2}^{P}$ was considered to be the
upper critical field of Pauli limited superconductivity until Fulde-Ferrell
(FF) [18] and Larkin-Ovchinnikov (LO) [19] pointed out that this critical
field can be enhanced further, by taking into account the possibility of an
inhomogeneous superconducting state in the high magnetic field regime. The
Pauli paramagnetic pair breaking effect is reduced by a new pairing state
(${\bf K}_{\uparrow}$, ${-{\bf K}+{\bf Q}}_{\downarrow}$) between the Zeeman
split Fermi surfaces, where ${\bf Q}\neq 0$ is the centre of mass momentum of
the Cooper pair. Thus, rather than a zero momentum superconducting state, the
system hosts a finite momentum superconducting state in the regime of high
magnetic field and low temperatures. We represent this possibility of finite
momentum pairing by the schematic shown in Figure 1(b). Below the tricritical
point (T∗) the transition is always first order, as was demonstrated later for
continuum two and three-dimensional systems [20, 21]. This exotic finite
momentum superconducting phase is known as the Fulde-Ferrell- Larkin-
Ovchinnikov (FFLO) state. The schematic presented in Figure 2 shows the finite
momentum pairing between the fermions in the FFLO state in comparison with the
zero momentum pairing in the BCS state. As is evident from the figure an FFLO
state comprises of a mismatch in the size of the Fermi surfaces corresponding
to the up and down fermion species. Such a Fermi surface mismatch gives rise
to unequal populations of the fermionic species i. e. a population imbalance.

Figure 2: Color online: Schematic diagram showing the Fermi surface mismatch
in the FFLO phase in comparison with the Fermi surface in the BCS state.

The symmetry breaking arising due to the finite momentum ${\bf Q}$ gives rise
to spatial modulations in the superconducting order parameter. There had been
several proposals regarding the form of this spatially modulated order
parameter with the first one suggested by Fulde and Ferrell being [18],

$\displaystyle\Delta({\bf r})$ $\displaystyle=$ $\displaystyle|\Delta|e^{i{\bf
Q}.{\bf r}}$ (2)

Here, the superconducting amplitude is homogeneous but the phase undergoes
periodic modulations in the real space. In a similar spirit Larkin and
Ovchinnikov suggested that a linear combination of two plane wave solutions
corresponding to degenerate superconducting states as [19],

$\displaystyle\Delta({\bf r})$ $\displaystyle=$
$\displaystyle|\Delta|(e^{i{\bf Q}.{\bf r}}+e^{-i{\bf Q}.{\bf
r}})=2|\Delta|\cos({\bf Q}.{\bf r})$ (3)

can give rise to a spatially modulated superconducting state, such that the
amplitude rather than the phase of the superconducting order parameter is
modulated. In general, depending upon the symmetry of the system more than two
plane waves can undergo linear combination to give rise to a spatially
modulated superconducting order parameter [19, 20, 21, 22, 23, 24, 25, 26,
27]. It has been found that for a two-dimensional thin film the linear
combination of two and three plane waves give rise to square and hexagonal
states, respectively [22, 27]. The spatial modulations give rise to nodes in
the superconducting order parameter, which serve as hosts to the unpaired
fermions arising out of the Pauli paramagnetic pair breaking. The FFLO phase
can thus be envisaged as a coexistent phase wherein finite momentum
superconductivity and finite magnetization (arising out of the unpaired
fermions) compliment each other in real space.

For a long time, the FFLO superconducting state though theoretically proposed
early on, failed to garner the interest it deserved. One of the principal
reasons for the same was its inaccessibility to the existing experimental
probes. The FFLO state in two-dimensional materials is highly susceptible to
impurities [28, 29], consequently the sample needs to be very clean in order
to observe the FFLO state. Even for a clean sample the applied magnetic field
should be applied precisely parallel to the plane because any small
perpendicular component of the magnetic field would give rise to Abrikosov
vortices and destroy the FFLO state [30]. Secondly, the transition temperature
of the FFLO state is strongly suppressed and the state gets rapidly destroyed
by thermal fluctuations even in three dimensions, which makes its experimental
detection difficult [31, 32]. It was however suggested that crystal anisotropy
of real materials can stabilize the FFLO phase against thermal fluctuations
[31].

Though the FFLO phase continues to be the most sought after superconducting
phase in systems with Fermi surface mismatch, there are other candidate phases
which can be stabilized in these systems, viz. Breached Pair (BP) [33, 34, 35]
and deformed Fermi surface superfluid phases [36, 37]. Unlike FFLO, both the
BP and deformed Fermi surface phases comprise of zero momentum pairing between
the fermions. The breached pair phase is characterized by the coexistence of
finite magnetization and “gapless” zero momentum superconductivity. The issue
of BP phase was first raised by Sarma in his seminal work where he discussed
the possibility of self consistent mean field solution with gapless mode in
$s$-wave superconductors, in presence of an applied magnetic field [33]. He
however found that such a phase even if it exists would be energetically
unfavorable as compared to the BCS state. Later Liu and Wilczek demonstrated
that while a contact interaction between the species in a population
imbalanced system does not allow the BP phase to be a stable ground state, the
same can indeed be achieved under non trivial circumstances such as, imbalance
in fermion effective masses with contact interaction between them, momentum
dependent interaction, same species repulsion etc [34, 35, 38, 39]. Another
possible mechanism to achieve stable ground state in an imbalanced Fermi
system is through deformed Fermi surfaces which spontaneously breaks the
rotational symmetry [36, 37].

##  0.2 Experimental signatures

The renewed interest on imbalanced superconducting systems is attributed to
two key experimental breakthroughs viz. (i) the advent of ultracold atomic
gases and (ii) improvisation of experiments on solid state materials. Before
going into the details of the experimental observations we highlight the
candidate systems amenable to the realization of imbalanced superconductivity
and superfluidity. In general, the term superconductivity is used in the
context of solid state systems containing charged electrons, while for the
charge neutral ultracold atomic gases, superfluidity is the correct choice of
nomenclature to describe the state with long range macroscopic phase
coherence. In the remaining sections of this chapter we will follow this
protocol to refer to the solid state and neutral cold atomic gases.

We mentioned in the earlier sections that in order to realize FFLO phase the
superconducting material should have large Maki parameter ($\alpha>$1). This
requirement is stringent and is achievable only under specific conditions.

  * •

In a three dimensional superconductor if H${}_{c2}^{orb}$ is significantly
larger than H${}_{c2}^{P}$ then it would lead to a large Maki parameter. In
heavy fermion superconductors such as CeCoIn5 the large mass of the electrons
give rise to low electron velocity and thus a high H${}_{c2}^{orb}$. The
resulting high Maki parameter $\alpha\sim 5$ makes it a suitable candidate to
realize FFLO superconductivity.

  * •

In layered two-dimensional materials the orbital effects are strongly
suppressed when the magnetic field is applied, which makes a large Maki
parameter plausible in such systems. A candidate material is quasi two-
dimensional organic superconductor $\kappa$-(BEDT-TTF)2Cu(NCS)2, with Maki
parameter of $\alpha\sim 8$.

  * •

If the system is a neutral Fermi gas, i. e. comprising of neutral atoms then
the applied magnetic field can couple only to the spin degrees of freedom.
This makes the ultracold atomic gases a promising candidate to realize FFLO
phase. In these systems the superfluid pairing takes place between the
fermions present in two hyperfine energy states of the neutral atoms or it can
take place between the fermions in the energy states of the atoms belonging to
two different atomic species or to two different isotopes of the same atomic
species.

We now discuss the experimental observations made on these candidate systems,
which suggests that these systems host exotic phases such as, FFLO and BP, in
presence of applied magnetic field or population imbalance of fermionic
species.

###  0.2.1 Heavy fermion superconductor:- CeCoIn5

One of the earliest experiments which suggested the possibility of FFLO phase
in CeCoIn5 (T${}_{c}\sim$2.3K) were the specific heat measurements reported by
Bianchi et al. [13]. They demonstrated that below a critical temperature
T${}_{0}\sim$ 0.3Tc and in presence of high magnetic field (H $>$ 10T) applied
either in-plane or out of plane of the material, the thermal transition of
superconducting to normal state changes from second to first order. They
identified a specific heat anomaly in this regime of first order transition
and attributed it to an underlying spatially modulated superconducting state.
Similar observations were earlier made by Tayama et al. [40] in their dc
magnetization measurements on CeCoIn5, which showed a first order phase
transition, indicating strong Pauli paramagnetic suppression of the
superconducting pairing in the regime of high magnetic field and low
temperatures (T${}_{0}\sim$ 0.3Tc). The possibility of FFLO phase in CeCoIn5
was further established by thermal conductivity [41] as well as by resistivity
and penetration depth measurements [42]. In presence of an in-plane applied
magnetic field the thermal conductivity data showed a discontinuous jump
indicating a change from second to first order phase transition below a
critical temperature of T${}_{0}\sim$1K. Similar signatures were detected by
Martin et al. in their resistivity measurements [42].

Based on the nuclear magnetic resonance (NMR) measurements carried out on
single crystal of CeCoIn5, Koutroulakis et al. [43] mapped out the (in-plane)
magnetic field-temperature phase diagram of this system. Based on the results
of NMR shift they identified three different regimes in the thermal phase
diagram. In the regime of low magnetic field (H$\lesssim$ 9.2T) a long range
ordered (LRO) homogeneous superconducting phase was identified. In the high
magnetic field regime (H$\geq$ 10.2T), in-commensurate spin density wave (IC-
SDW) order was found to coexist with FFLO at low temperatures. The
intermediate magnetic field regime (9.2T $\leq$ H $<$ 10.2T) was inferred to
comprise of FFLO phase surviving with antiferromagnetic (AFM) correlations.
This phase was found to undergo first order thermal transition. Similar
inferences were made earlier by Kumagai et al. [44] via NMR Knight shift
measurements on CeCoIn5. They observed that the FFLO superconductivity is
realizable for both in-plane and out of plane magnetic fields and the
corresponding phase diagrams are qualitatively different from each other.

By carrying out high magnetic field neutron scattering experiments on CeCoIn5
Gerber et al. [45] identified a “Q-phase” in the high magnetic field low
temperature regime, corresponding to a finite momentum pairing. They further
confirmed that the finite momentum pairing coexists with an SDW order and
undergoes first order thermal phase transition at H$\sim 9.8$T. In a
relatively recent work Kim et al. [46] suggested that though there is indeed a
Q-phase in the high field low temperature regime of the CeCoIn5 phase diagram
in presence of an in-plane applied magnetic field, it does not correspond to a
FFLO phase. It was inferred based on their thermal conductivity measurements
that the Q-phase comprises of intertwinned SDW, $d$-wave superconducting and
inhomogeneous $p$-wave pair density wave (PDW) orders. Evidently, the true
nature of the Q-phase seems to be unsettled. A very recent thermal
conductivity measurement by the same group suggests that the characteristic of
the Q-phase is intimately connected to the direction of the applied magnetic
field [47]. For a magnetic field applied along the ab-plane of the crystal the
Q-phase is likely to host a SDW order. On the other hand for a field applied
perpendicular to the ab-plane a FFLO order is realizable in the high field low
temperature Q-phase of CeCoIn5 [47].

###  0.2.2 Layered organic superconductor:-
$\kappa$-(BEDT-TTF)2Cu(NCS)2

In two-dimensional layered organic superconductor $\kappa$-BEDT (T${}_{c}\sim$
9.5K), the FFLO phase was inferred from specific heat measurements in presence
of applied in-plane magnetic field [14]. Unlike CeCoIn5 it was observed that
in $\kappa$-BEDT the thermal transitions from both the zero momentum and FFLO
superconducting phases are of first order. A Zeeman field driven phase
transition was detected in $\kappa$-BEDT via NMR measurements by Wright et al.
[48]. They suggested that this transition and the associated discontinuity in
the derivative of magnetization, at H$\sim$213$\pm$3kOe amounts to a possible
indication of spatially inhomogeneous superconducting phase. Similar inference
was made by Bergk et al. [49] and Agosta et al. [50] based on magnetic torque
and penetration depth measurements, respectively. The high magnetic field (H
$\geq$ 21T) regime was suggested to host FFLO superconducting phase. The most
clear evidence of FFLO phase in $\kappa$-BEDT was given by Mayaffre et al.
through NMR measurements [51]. They observed that the NMR spin-lattice
relaxation rate gets significantly enhanced in the high magnetic field regime,
and attributed this observation to the Andreev bound states arising out of the
unpaired quasiparticles residing in the nodes of the spatially modulated FFLO
state.

In addition, the presence of FFLO superconducting phase has been recently
proposed in iron superconductors such as KFe2As2, LiFeAs etc. [15, 52, 53],
which are suggested to be two-band superconductors.

###  0.2.3 Cold atomic gases and Fermi superfluids

Over the past decade ultracold Fermi gases have proved to be the testbed of
research on imbalanced superfluid systems, owing to the flexibility of tuning
the interaction parameters (via Feshbach resonance), optical lattice and trap
geometries, that these systems provide. Unlike the solid state superconducting
systems where a magnetic field is applied, population imbalance in ultracold
atomic gases is created by loading unequal populations of fermionic species in
the optical lattice. Alternatively, equal populations of different atomic
species with unequal effective masses can be loaded in the optical lattice to
create mass imbalanced Fermi gases. Clear evidence of superfluid pair
correlations have been found in case of both balanced [54, 55] as well as
imbalanced [16] Fermi gases. Most of these experiments are carried out at or
close to unitarity.

The concept of unitarity is intimately related to the $s$-wave scattering
length which can be controlled via the Feshbach resonance. In three-
dimensional (3D) continuum Fermi gas the scattering length
$a_{3D}\rightarrow\infty$ as $g\rightarrow g_{c}$. Here, $g_{c}$ is finite in
3D and defines the interaction strength at which the first two-body bound
state formation takes place in vacuum. The dimensionless coupling
$1/k_{F}a_{3D}=0$ (where, $k_{F}$ is the Fermi momentum) at $g_{c}$ and also
corresponds to the regime of maximum Tc of the 3D Fermi gas. For two-
dimensional (2D) Fermi gases unitarity corresponds to the coupling strength at
which $ln(k_{F}a_{2D})\rightarrow 0$ ($a_{2D}$ is $s$-wave scattering length
in 2D), which once again corresponds to the interaction for maximum Tc. For a
lattice the notion of $k_{F}$ is not very meaningful particularly if one is
away from the regime of small filling. In the context of lattice fermion
models, the crossover regime or unitarity corresponds to the interaction
strength at which Tc is maximum, where neither the fermionic nor bosonic
description of the system holds true completely [54, 56, 57, 58]. We will
revisit the issue of unitarity in the later sections of this chapter, when we
discuss the BCS-BEC crossover picture. We now discuss some of the important
experiments on imbalanced Fermi gas in ultracold cold atomic gas setups, and
their principal observations.

Wolfgang Ketterle and his group mapped out the thermal phase diagram of
strongly interacting imbalanced spin mixture created between two hyperfine
states of 6Li atoms [59]. Based on radio frequency (rf) spectroscopy
measurements they identified a gap in the single particle excitation spectra
which indicated that the minority species is paired up. It was shown that in
presence of large imbalance, at high temperature the spectra shows only the
atomic peak which was considered to be the evidence of unpaired fermions in
the system. As the temperature is lowered, in addition to the atomic a pairing
peak also appeared in the spectra indicating a coexistent phase with
superfluid pairs and unpaired fermions. Further decreasing the temperature
left behind only the pairing peak. The work showed high temperature
coexistence of paired and unpaired fermions but did not find any indication of
superfluidity at low temperatures in the regime of large imbalance. Note that
at large imbalance there would not be enough minority spins to pair up with
the majority species. Consequently, signatures of pair correlations (if any)
would be strongly suppressed.

Shin et al. [60] reported the experimental observation of phase separation
between superfluid and normal regions in a strongly interacting Fermi gas of
6Li with population imbalance. They used a special phase contrast imaging
technique and 3D image reconstruction to determine the density difference
between the two components of the gas in an optical trap. A shell structure
was observed, in which the superfluid region containing equal densities of the
components is surrounded by a normal region of unequal densities. Furthermore,
using phase contrast imaging they characterized the normal to superfluid phase
transition.

A quantum phase transition was detected in the unitary imbalanced Fermi gas of
6Li at a critical polarization (magnetization) of $\sim 36\%$ by Shin et al.
[16]. Further, they detected the change in the order of thermal phase
transition as a function of polarization. In the regime of low temperatures
and increasing polarization a first order phase transition was detected via a
jump in the polarization. On the other hand the transition at high
temperatures and weak polarization was found to be of second order. These two
orders of phase transition are connected at a tricritical point. While they
carried out the experiment in presence of trapping potential, using
tomographic reconstruction of local Fermi temperatures and spin polarization
they determined the phase diagram of the homogeneous system free from any
inhomogeneity arising due to the trapping potential. The resulting phase
diagram comprised of a superfluid phase in the regime of weak polarization
characterized by spatial coexistence of superfluid order and non zero
polarization at finite temperatures, akin to the BP phase. This regime
undergoes second order thermal phase transition to a normal phase. The large
polarization regime is a phase separated (unstable) regime which undergoes a
first order thermal phase transition. At still larger polarizations the system
is in the normal phase. Note that signatures of FFLO phase was not detected in
this experiment, however Shin et al. did not rule out the possibility of such
a phase at large polarizations [16].

In a different experimental set up, containing arrays of one-dimensional tubes
constructed using a two- dimensional optical lattice, Liao et al. [61]
investigated the phases of an imbalanced Fermi gas of 6Li. Due to nesting of
Fermi surfaces the FFLO phase is expected to have a larger regime of stability
in one dimension as compared to the higher dimensions. Liao et al. [61] mapped
out the phase diagram of this system in polarization-temperature plane and
demonstrated that, at low imbalance a partially polarized phase is realized at
the centre of the trap, which radially spreads out with increasing imbalance.
This partially polarized region extends to the edge of the cloud at a critical
imbalance. Further increase in imbalance leads the edges of the cloud to
become completely polarized. Liao et al. found their experimental phase
diagram to be in agreement with theoretical predictions and suggested that
even though signatures of FFLO phases was not visible in their experiment it
was certainly a possibility [61].

Experimental investigations of mass imbalanced Fermi-Fermi mixtures are
relatively fewer. A mass imbalanced Fermi-Fermi mixture is achievable for
example in a 6Li-40K mixture. While superfluidity in such a system is yet to
be attained in experiments, the Fermi degenerate regime [62, 63] as well as
the Feshbach resonance between 6Li and 40K atoms [63, 64, 65] and formation of
6Li-40K heteromolecules [66] are already a reality. Furthermore, experimental
realization of mixtures of other fermion species (such as, 161Dy, 163Dy,
167Er) are expected in future [67, 68]. It has been reported that for a
double-degenerate 6Li-40K mixture the Fermi temperatures are TLi=390 nK and
TK=135 nK, for Li and K species, respectively [69]. In comparison, for a
balanced Fermi gas of 6Li, the Fermi temperature is known to be TF=1.0$\mu$K
[70] with the corresponding Tc scale being T${}_{c}\sim$ 0.15TF [71]. In other
words, imbalance in mass strongly suppresses the superfluid Tc of the Fermi-
Fermi mixtures, making their experimental detection non trivial. However,
short range superfluid pair correlations are likely to survive at temperatures
significantly higher than the superfluid Tc and accessible to experimental
probes.

##  0.3 Theoretical formalism

The experimental observations across different classes of solid state
materials and ultracold atomic gases lead to a flurry of activities in the
theoretical front. The theoretical studies can broadly be classified in two
categories as, (i) the lattice fermion models and (ii) the continuum models.
In this chapter we will mainly focus on (i). For a detail discussion on
continuum models with imbalance the readers are referred to [72, 73, 74].
Before getting into the details of the results obtained from different lattice
fermion based calculations, let us first set up the theoretical framework
which can be used to investigate imbalanced Fermi superfluids and
superconductors. For this we resort to the simplest tool of mean field theory.
We will use the framework of BCS theory and Bogoliubov-Valatin transformation
to set up the mean field formalism. The generic Hamiltonian of the two
component Fermi gas with an on-site interaction between the fermion species is
used to model the system. The BCS mean field approximation is used such that
the order parameters are striped of all fluctuations. We take into account the
imbalance between the population of the fermion species and allow for the
possibilities of exotic superfluid phases such as, FFLO and BP, in addition to
the homogeneous BCS phase. Bogoliubov-Valatin transformations [8] are used to
diagonalize the Hamiltonian and the resulting eigenvalues and eigenfunctions
are used to compute different observables relevant from the point of view of
experiments. Once the basic formalism is in place we will discuss the
modifications required to address different families of imbalanced
superconductors and superfluids.

###  0.3.1 Mean field theory

We consider the Hamiltonian of a two component Fermi gas comprising of up
($\uparrow$) and down ($\downarrow$) fermion species, with an attractive
interaction between them. For the sake of simplicity we do not take into
account any external potential acting on the fermions. The Hamiltonian
corresponding to such a system reads as,

$\displaystyle\hat{H}$ $\displaystyle=$ $\displaystyle\int d{\bf
r}\sum_{\sigma}\hat{\psi}_{\sigma}^{\dagger}({\bf r})(\hat{K}_{\sigma}({\bf
r})-\mu_{\sigma})\hat{\psi}_{\sigma}({\bf r})$ (4) $\displaystyle+\int\int
d{\bf r}d{\bf r}^{\prime}V_{\uparrow\downarrow}({\bf r},{\bf
r}^{\prime})\hat{\psi}_{\uparrow}^{\dagger}({\bf
r})\hat{\psi}_{\downarrow}^{\dagger}({\bf
r}^{\prime})\hat{\psi}_{\downarrow}({\bf
r}^{\prime})\hat{\psi}_{\uparrow}({\bf r})$

where, $\sigma=\uparrow,\downarrow$ correspond to the two fermion species.
$V_{\uparrow\downarrow}({\bf r},{\bf r}^{\prime})$ is the interaction between
the species, with ${\bf r}$ and ${\bf r}^{\prime}$ being the position
coordinates, corresponding to a two dimensional system. The kinetic energy
contribution to the Hamiltonian is given as, $\hat{K}_{\sigma}({\bf
r})=-\frac{\hbar^{2}}{2m_{\sigma}}\frac{\partial^{2}}{\partial{\bf r}^{2}}$
and $\mu_{\sigma}$ is the species dependent chemical potential. Note that the
effective mass ($m_{\sigma}$) in the kinetic energy contribution is dependent
on the species. For a system with only imbalance in fermionic populations
$m_{\uparrow}=m_{\downarrow}=m$ since the mass of the fermions of both the
species are the same. On the other hand, for a system comprising of two
different isotopes of the same atom or of two different atomic species,
$m_{\uparrow}\neq m_{\downarrow}$. Such systems can be realized in ultracold
atomic gases and are called “mass imbalanced” Fermi-Fermi mixtures. For the
present discussion we consider equal effective masses for both the species and
set $m_{\uparrow}=m_{\downarrow}=m$.

We simplify the model by approximating the interaction to be a contact
interaction as,

$\displaystyle V_{\uparrow\downarrow}({\bf r},{\bf r}^{\prime})$
$\displaystyle=$ $\displaystyle V_{0}\delta({\bf r}-{\bf r}^{\prime})$ (5)

This approximation is particularly good for ultracold atomic gases in which
the interaction is short ranged [75]. With this approximation the Hamiltonian
takes the form,

$\displaystyle\hat{H}$ $\displaystyle=$ $\displaystyle\int d{\bf
r}\sum_{\sigma}\hat{\psi}_{\sigma}^{\dagger}({\bf r})(\hat{K}_{\sigma}({\bf
r})-\mu_{\sigma})\hat{\psi}_{\sigma}({\bf r})$ (6) $\displaystyle+V_{0}\int
d{\bf r}\hat{\psi}_{\uparrow}^{\dagger}({\bf
r})\hat{\psi}_{\downarrow}^{\dagger}({\bf r})\hat{\psi}_{\downarrow}({\bf
r})\hat{\psi}_{\uparrow}({\bf r})$

We consider an underlying lattice potential and for that expand the continuum
Fermi operators in terms of Wannier functions [76] as,

$\displaystyle\hat{\psi}_{\sigma}^{\dagger}({\bf r})$ $\displaystyle=$
$\displaystyle\sum_{{\bf n},i}\hat{c}_{\sigma,{\bf n},i}^{\dagger}w_{\bf
n}^{*}({\bf r}-{\bf r}_{i})$ (7)

where, $\hat{c}_{\sigma,{\bf n},i}^{\dagger}$ creates a $\sigma$ particle in
the lattice site $i$ in the energy band ${\bf n}$, and $w_{\bf n}^{*}({\bf
r}-{\bf r}_{i})$ are the Wannier functions. The Wannier functions can be
approximated by harmonic oscillator states in the tight-binding limit. At low
temperatures and weak interactions only the lowest energy band ${\bf n}=0$
would be occupied by the fermions, we can thus omit the band index henceforth.
Using the above expansion the continuum Hamiltonian is re-written as,

$\displaystyle\hat{H}$ $\displaystyle=$
$\displaystyle-\sum_{\sigma}\sum_{\langle
ij\rangle}t_{\sigma,ij}\hat{c}_{\sigma,i}^{\dagger}\hat{c}_{\sigma,j}+\sum_{\sigma,i}(\epsilon_{\sigma}-\mu_{\sigma})\hat{c}_{\sigma,i}^{\dagger}\hat{c}_{\sigma,i}$
(8)
$\displaystyle+U\sum_{i}\hat{c}_{\uparrow,i}^{\dagger}\hat{c}_{\downarrow,i}^{\dagger}\hat{c}_{\downarrow,i}\hat{c}_{\uparrow,i}$

Here, $\langle ij\rangle$ corresponds to the summation over nearest neighbors
of a lattice and $t_{\sigma,ij}$ is the hopping amplitude between the nearest
neighbors; $\epsilon_{\sigma}$ is the onsite energy and $U$ is the on-site
interaction, and are expressed as [77],

$\displaystyle t_{\sigma,ij}$ $\displaystyle=$ $\displaystyle-\int d{\bf
r}w^{*}({\bf r}-{\bf r}_{i})\hat{K}_{\sigma}({\bf r})w({\bf r}-{\bf r}_{j})$
(9) $\displaystyle\epsilon_{\sigma}$ $\displaystyle=$ $\displaystyle\int d{\bf
r}w^{*}({\bf r}-{\bf r}_{i})\hat{K}_{\sigma}({\bf r})w({\bf r}-{\bf r}_{i})$
(10) $\displaystyle U$ $\displaystyle=$ $\displaystyle V_{0}\int d{\bf
r}|w({\bf r}-{\bf r}_{i})|^{4}$ (11)

Equation (8) corresponds to the standard single band Hubbard model for a two
component Fermi gas. In order to take into account the imbalance in the
population of fermionic species, we set $\mu_{\uparrow}\neq\mu_{\downarrow}$.
We consider the system to be mass balanced and the hopping amplitude to be a
constant, such that $t_{\uparrow,ij}=t_{\downarrow,ij}=t_{ij}=t$. We also
consider $\epsilon_{\uparrow}=\epsilon_{\downarrow}=\epsilon=0$, i. e. the on-
site energy is taken to be a constant and set to zero. The resulting
Hamiltonian thus takes the form,

$\displaystyle\hat{H}$ $\displaystyle=$
$\displaystyle-t\sum_{\sigma}\sum_{\langle
ij\rangle}\hat{c}_{\sigma,i}^{\dagger}\hat{c}_{\sigma,j}-\sum_{\sigma,i}\mu_{\sigma}\hat{c}_{\sigma,i}^{\dagger}\hat{c}_{\sigma,i}+U\sum_{i}\hat{c}_{\uparrow,i}^{\dagger}\hat{c}_{\downarrow,i}^{\dagger}\hat{c}_{\downarrow,i}\hat{c}_{\uparrow,i}$

Here, $U<0$ corresponding to the attractive interaction between the fermions.

#### BCS theory of uniform superconductors

Equation (12) comprises of a four-fermion interaction term which is not
readily solvable. In order to make the Hamiltonian tractable we use the BCS
mean field approximation. Within this approximation we introduce a pairing
field ($\hat{\Delta}_{i}$) to describe the Cooper pairs. We do not go into the
details of the BCS theory here and guide the interested readers to [1, 8, 78,
79]. The pairing field thus introduced serves as the superconducting “order
parameter” and is related to the expectation value of
$\hat{\psi}_{\downarrow}\hat{\psi}_{\uparrow}$ as,

$\displaystyle\langle\hat{\Delta}_{i}\rangle\equiv\Delta_{i}$ $\displaystyle=$
$\displaystyle V_{0}\langle\hat{\psi}_{\downarrow}({\bf
r})\hat{\psi}_{\uparrow}({\bf
r})\rangle=U\langle\hat{c}_{i\downarrow}\hat{c}_{i\uparrow}\rangle$ (13)

The standard BCS theory describes the transition from a normal Fermi gas to a
zero momentum uniform superconducting state, and the pairing field is
considered to be a constant, $\Delta_{i}=\Delta$. In order to capture the
inhomogeneous superconducting state such as FFLO or BP, one needs to retain
the spatial dependence of the pairing field.

Note that the pairing field defined so far is on-site and does not take into
account any directional dependence. In the language of superconducting pairing
state symmetry, such an on-site pairing correspond to a $s$-wave
superconducting state. For the sake of simplicity we will discuss the
theoretical formalism in detail for a $s$-wave superconducting state and will
outline the modifications of this formalism for a non $s$-wave pairing state,
at the end of this section.

The BCS approximation in addition to the pairing field introduces Hartree
field corresponding to the fermionic number densities, which is defined as,
$n_{\sigma}=\langle\hat{c}_{\sigma}^{\dagger}\hat{c}_{\sigma}\rangle$.
Including the Hartree field in the formalism renormalizes the chemical
potential as,

$\displaystyle(-\mu_{\sigma}+Un_{-\sigma})\hat{c}_{\sigma}^{\dagger}\hat{c}_{\sigma}$
$\displaystyle\equiv$
$\displaystyle-\tilde{\mu}_{\sigma}\hat{c}_{\sigma}^{\dagger}\hat{c}_{\sigma}$
(14)

Here, $\tilde{\mu}_{\sigma}$ is the renormalized chemical potential and the
densities are considered to be uniform. In terms of the average pairing and
the Hartree fields the mean field Hamiltonian describing a $s$-wave
superconducting pairing of the fermions in a two component Fermi gas, reads
as,

$\displaystyle\hat{H}$ $\displaystyle=$
$\displaystyle-t\sum_{\sigma}\sum_{\langle
ij\rangle}\hat{c}_{\sigma,i}^{\dagger}\hat{c}_{\sigma,j}-\sum_{\sigma}\sum_{i}\tilde{\mu}_{\sigma}\hat{c}_{\sigma,i}^{\dagger}\hat{c}_{\sigma,i}$
(15)
$\displaystyle+\sum_{i}(\hat{c}_{\uparrow,i}^{\dagger}\hat{c}_{\downarrow,i}^{\dagger}\Delta_{i}+\Delta_{i}^{*}\hat{c}_{\downarrow,i}\hat{c}_{\uparrow,i}-\frac{|\Delta_{i}|^{2}}{U})$

In the Hamiltonian given by Equation (15) the pairing field is replaced by its
approximated average value (see Equation 13), but we do not assume any
specific form of this average.

Note that the approximation of uniform density should be used carefully. The
approximation breaks down in case of superconductors with spatial
inhomogeneities, such as those arising in a FFLO of BP phase or in disordered
superconductors. In such cases the correct tool of choice is Bogoliubov-de-
Gennes (BdG) formalism which takes into account the spatial inhomogeneities of
the pairing field amplitude ($|\Delta_{i}|$) as well as of the number
densities [1, 8].

#### Mean field theory of FFLO superconductivity

The FFLO phase comprises of spatially modulated superconducting state
coexistent with non zero magnetization. The corresponding Cooper pair centre
of mass momentum ${\bf Q}$ is determined by minimizing the energy of the
system. While the FF state corresponds to a superconducting state with
modulation in the pairing field phase ($\Delta({\bf r})=\Delta e^{i{\bf
Q}.{\bf r}}$), the LO ($\Delta({\bf r})=\Delta\cos({\bf Q}.{\bf r})$) state on
the other hand contains modulation in the pairing field amplitude.

There is now a consensus that at the ground state of two component population
imbalanced Fermi gas the LO phase is energetically favorable as compared to
the FF phase, both in case of lattice as well as continuum systems [21, 80,
81, 82]. An exception however is the population imbalanced system with spin-
orbit coupling, where the FF phase is found to have lower energy as compared
to that of the LO phase. Strictly speaking, the regime of stability of the FF
and LO phases should depend upon the system under consideration, however for
many systems they are found to overlap.

In this section we consider the simpler FF ansatz to discuss the mean field
theory of the spatially modulated superconducting phase. The theory can
however be easily generalized to describe the LO superconducting phase. Using
the ansatz for the FF state ($\Delta_{i}=\Delta_{i}e^{i{\bf Q}.{\bf
r}}=\Delta_{i}e^{i(Q_{x}x+Q_{y}y)}$) the superconducting Hamiltonian for the
population imbalanced two component Fermi gas in two-dimensions is given as,

$\displaystyle\hat{H}$ $\displaystyle=$
$\displaystyle-t\sum_{\sigma}\sum_{\langle
ij\rangle}\hat{c}_{\sigma,i}^{\dagger}\hat{c}_{\sigma,j}-\sum_{\sigma}\sum_{i}\tilde{\mu}_{\sigma}\hat{c}_{\sigma,i}^{\dagger}\hat{c}_{\sigma,i}$
$\displaystyle+\sum_{i}(\Delta_{i}e^{i(Q_{x}.x+Q_{y}.y)}\hat{c}_{\uparrow,i}^{\dagger}\hat{c}_{\downarrow,i}^{\dagger}+\Delta_{i}e^{-i(Q_{x}.x+Q_{y}.y)}\hat{c}_{\downarrow,i}\hat{c}_{\uparrow,i}-\frac{\Delta_{i}^{2}}{U})$

Here, $Q_{x}$ and $Q_{y}$ correspond to the components of the Cooper pair
centre of mass momentum ${\bf Q}$, in two dimensions. In order to determine
the superconducting gap and number equations this Hamiltonian needs to be
diagonalized. Since the Hamiltonian is not diagonal in the electronic basis we
carry out a basis transformation using Bogoliubov-Valatin transformations
defined as [8],

$\displaystyle\hat{c}_{\uparrow,i}$ $\displaystyle=$
$\displaystyle\sum_{n}(u_{n,\uparrow}^{i}\gamma_{n,\uparrow}-v_{n,\uparrow}^{*i}\gamma_{n,\downarrow}^{\dagger})$
(17) $\displaystyle\hat{c}_{\downarrow,i}$ $\displaystyle=$
$\displaystyle\sum_{n}(u_{n,\downarrow}^{i}\gamma_{n,\downarrow}+v_{n,\downarrow}^{*i}\gamma_{n,\uparrow}^{\dagger})$
(18)

Here, $n$ is a complete set of states, $u_{n}$ and $v_{n}$ are the
eigenfunctions of the Hamiltonian $\hat{H}$, which satisfies the condition
$u_{n}^{i2}+v_{n}^{i2}=1$ [8]. The new fermionic operators $\gamma_{n}$ are
quasiparticle operators which satisfies the anti commutation relation
$\gamma_{n}^{\dagger}\gamma_{n}+\gamma_{n}\gamma_{n}^{\dagger}=1$ and are
chosen in a way such that,

$\displaystyle\hat{H}$ $\displaystyle=$ $\displaystyle
E_{0}+\sum_{n,\sigma}\epsilon_{n}\gamma_{n,\sigma}^{\dagger}\gamma_{n,\sigma}$
(19)

Thus, we have,

$\displaystyle[\hat{H},\gamma_{n,\sigma}]$ $\displaystyle=$
$\displaystyle-\epsilon_{n}\gamma_{n,\sigma}$ (20)
$\displaystyle[\hat{H},\gamma_{n,\sigma}^{\dagger}]$ $\displaystyle=$
$\displaystyle\epsilon_{n}\gamma_{n,\sigma}^{\dagger}$ (21)

where, $\epsilon_{n}$ corresponds to the eigenvalues of the Hamiltonian
$\hat{H}$. Using the above expressions one can determine the coefficients
$u_{n}^{i}$ and $v_{n}^{i}$ by solving the matrix eigenvalue equation. For
this we define the matrices A and $\Psi$, such that,

$\displaystyle A$ $\displaystyle=$
$\displaystyle\begin{bmatrix}-t-\mu_{\uparrow}&\hat{\Delta}_{i}\\\
\hat{\Delta}_{i}^{*}&t+\mu_{\downarrow}\end{bmatrix}$ (22)

and

$\displaystyle\Psi$ $\displaystyle=$
$\displaystyle\begin{bmatrix}\\{u_{n}\\}\\\ \\{v_{n}\\}\end{bmatrix}$ (23)

The resulting eigenvalue equation is then given as,

$\displaystyle A\Psi$ $\displaystyle=$ $\displaystyle\epsilon_{n}\Psi$ (24)

In order to solve the eigenvalue equations, one first needs to determine the
superconducting pairing field ($\Delta_{i}$) as well as the chemical
potentials $\mu_{\uparrow}$ and $\mu_{\downarrow}$. For this the set of
coupled equations corresponding to $\Delta_{i}$ and fermionic number densities
$n_{i\sigma}$ is to be solved self consistently. In terms of the BdG
eigenfunctions and eigenvalues the self consistent mean field equations are
given as,

$\displaystyle\Delta_{i}$ $\displaystyle=$ $\displaystyle
U\sum_{n}\langle\hat{c}_{i,\downarrow}\hat{c}_{i,\uparrow}\rangle$ (25)
$\displaystyle=$ $\displaystyle
U\sum_{n}\\{v_{n,\downarrow}^{*i}u_{n,\uparrow}^{i}f(\epsilon_{n})+u_{n,\downarrow}^{*i}v_{n,\uparrow}^{i}f(\epsilon_{-n})\\}$
$\displaystyle n_{i\uparrow}$ $\displaystyle=$
$\displaystyle\sum_{n}\langle\hat{c}_{i,\uparrow}^{\dagger}\hat{c}_{i,\uparrow}\rangle$
(26) $\displaystyle=$
$\displaystyle\sum_{n}\\{|u_{n,\uparrow}^{i}|^{2}f(\epsilon_{n})+|v_{n,\uparrow}^{i}|^{2}f(\epsilon_{-n})\\}$
$\displaystyle n_{i\downarrow}$ $\displaystyle=$
$\displaystyle\sum_{n}\langle\hat{c}_{i,\downarrow}^{\dagger}\hat{c}_{i,\downarrow}\rangle$
(27) $\displaystyle=$
$\displaystyle\sum_{n}\\{|u_{n,\downarrow}^{i}|^{2}[1-f(\epsilon_{n})]+|v_{n,\downarrow}^{i}|^{2}[1-f(\epsilon_{-n})]\\}$

where, $f(\epsilon_{n})$ and $f(\epsilon_{-n})$ are the Fermi functions
corresponding to the positive and negative eigenvalues, respectively. Starting
from an initial guess these equations are iterated to obtain self consistent
values of the pairing field and the number densities.

For a lattice of dimension $L\times L$ the dimension of the matrix is
$2N\times 2N$, where, $N=L^{2}$. The net chemical potential of a population
imbalanced Fermi gas is defined as,
$\mu=\frac{\mu_{\uparrow}+\mu_{\downarrow}}{2}$. The difference between the
chemical potentials act as a Zeeman field defined as,
$h=\frac{\mu_{\uparrow}-\mu_{\downarrow}}{2}$. In terms of $\mu$ and $h$ we
can thus write the spin resolved chemical potentials as,
$\mu_{\uparrow}=\mu+h$ and $\mu_{\downarrow}=\mu-h$. While working in grand
canonical ensemble we work at a fixed chemical potential $\mu$ and tune the
Zeeman field $h$ such that, the system undergoes phase transition between an
uniform superconducting phase and inhomogeneous FFLO or BP phases, at a fixed
interaction $U$.

By examining the real space self consistent mean field solutions of the
pairing field and fermionic number densities one can analyze the underlying
state. A non zero $\Delta_{i}$ corresponds to a superconducting state, while a
state with $\Delta_{i}=0$ is non superconducting. In the same spirit,
$m_{i}=(n_{i\uparrow}-n_{i\downarrow})\neq 0$ corresponds to a state with
finite magnetization. Before analyzing the self consistent solution it is
however essential to determine whether the particular solution is stable, in
other words whether it corresponds to the minimum energy state for the given
choice of parameters.

We thus next calculate the free energy of the system and for that we go back
to Equation 19, which gives the Hamiltonian in the quasiparticle basis. Here
$E_{0}$ correspond to the ground state energy and in order to determine
$E_{0}$ we note that,

$\displaystyle\langle\Psi_{G}|\hat{H}|\Psi_{G}\rangle$ $\displaystyle=$
$\displaystyle E_{0}$ (28)

where, $|\Psi_{G}\rangle$ correspond to the ground state. The free energy of
the system reads as,

$\displaystyle F$ $\displaystyle=$ $\displaystyle
E_{0}+\sum_{\epsilon_{n}}(-\epsilon_{-n}-\frac{1}{\beta}[\ln(1+e^{-\beta\epsilon_{n}})+\ln(1+e^{-\beta\epsilon_{-n}})]$
(29) $\displaystyle-\frac{\Delta_{i}^{2}}{U})$

The global minima of this free energy functional gives the stable ground state
of the system. Based on the pairing field amplitude and average magnetization
($\langle m_{i}\rangle$) the ground state of two component imbalanced Fermi
gas can be broadly classified as, (i) uniform (BCS) superconducting phase
($\Delta_{i}({\bf Q}=0)\neq 0$, $\langle m_{i}\rangle=0$), (ii) (partially)
polarized Fermi liquid (PPFL) ($\Delta_{i}=0$, $\langle m_{i}\rangle\neq 0$),
(iii) FFLO phase ($\Delta_{i}({\bf Q}\neq 0)\neq 0$, $\langle m_{i}\rangle\neq
0$), (iv) BP phase ($\Delta_{i}({\bf Q}=0)\neq 0$, $\langle m_{i}\rangle\neq
0$). Note that even though both BCS and BP phases are characterized by ${\bf
Q}=0$ pairing state, the BP phase comprises of gapless excitations unlike the
BCS superconducting state [33, 34, 35, 38, 39]. These gapless excitations are
essential to accommodate the unpaired fermions which in turn gives rise to
finite magnetization.

###  0.3.2 Beyond the mean field approximation

So far we have taken into account only the inhomogeneity in the amplitude of
the pairing field (BdG formalism). The complex superconducting pairing field
defined as, $\Delta_{i}=|\Delta_{i}|e^{i\theta_{i}}$ however comprises of
amplitude and phase, both of which in principle can fluctuate and become
inhomogeneous. The fluctuations in the pairing field amplitude and phase
becomes progressively important with increasing temperature where thermal
fluctuations can serve as pair breaking agent. In order to understand the
importance of fluctuations on the superconductivity of a system let us take a
detour to briefly summarise the BCS-BEC crossover in Fermi gases [73, 83, 84,
85].

#### BCS-BEC crossover

Figure 3 shows the schematic phase diagram of the BCS-BEC crossover of the two
component Fermi gas comprising of equal populations of up and down fermion
species, in the interaction-temperature (U-T) plane. The phase diagram
comprises of four distinct regimes viz. superconductor, (bosonic) insulator,
Fermi liquid, non Fermi liquid and two thermal scales as Tc and Tpg. While in
a solid state system tuning the interaction strength is non trivial in
ultracold atomic gases, Feshbach resonance is used to control the $s$-wave
scattering length and in turn the interspecies interaction [75, 83, 86]. The
ground state of the BCS-BEC crossover can be captured by the BCS theory with
suitably modified coefficients, as shown by Leggett [87]. The two thermal
scales that we show in this phase diagram can naively be understood to be the
one determined through the experiments on ultracold atomic gases (Tc) and the
one determined through BCS-meanfield approximation (Tpg). Evidently, they do
not agree with each other except for the weak coupling regime (U$\ll$t) i. e.
deep inside the BCS regime.

Figure 3: Color online: Schematic phase diagram showing the BCS-BEC crossover
picture in lattice fermion model.

According to the BCS mean-field theory, superconductivity of a system vanishes
due to the suppression of the superconducting pairing field amplitude. This
obeys the thermal scale T${}_{BCS}\sim e^{-1/U}$ and the corresponding
superconducting pairing field amplitude varies as $\Delta\sim$ TBCS. Thus,
TBCS increases linearly with increasing $\Delta$ and in turn with $U$ which
governs the magnitude of $\Delta$ [1, 2, 8]. As shown in Figure 3, Tpg bears
out this linear relation with respect to $U$. The agreement of Tpg with the
experimentally determined Tc in the weak coupling regime suggests that the BCS
theory is well suited to capture the finite temperature behavior of the system
in the weak coupling regime. The high temperature regime which can be captured
by the BCS theory is marked as the Fermi liquid.

Unlike the BCS prediction the experimental Tc however behaves non
monotonically and in the strong coupling regime T${}_{c}\sim 1/U$, as shown in
Figure 3. It is obvious that the BCS theory ceases to be applicable in the
strong coupling regime where the pair of fermions tightly bound together and
undergoes Bose-Einstein condensation, unlike the weak coupling regime where
the Cooper pairs are characterized by large coherence length [85, 88, 89, 90].
An efficient way to understand the BEC regime of this crossover picture is by
mapping the fermionic attractive Hubbard model on to an equivalent spin model.
It is now established that the attractive Hubbard model yields an effective
field theoretic description in terms of an $SO(3)$ non linear sigma model. In
the strong coupling regime this model maps on to the Heisenberg model in an
uniform field. Except for the regime of weak magnetic field (i. e. low filling
of the attractive Hubbard model) the Heisenberg model reduces to the quantum
XY model [91, 92]. In the regime of strong coupling, the system undergoes
transition from the superfluid state to a bosonic insulator with increasing
temperature. The bosonic insulator is characterized by phase uncorrelated
large pairing field amplitudes and a hard single particle spectral gap at the
Fermi level.

The regime of intermediate coupling is the most complex part of this phase
diagram and none of the perturbative schemes are suitable to capture its
behavior. In this regime marked as the “non Fermi liquid” in Figure 3, neither
a completely fermionic nor a bosonic description holds true because the phase
comprises of coexisting unpaired fermions and preformed pairs (bosons). The
intuitive basis of understanding the BCS-BEC crossover was provided by the
earlier works of Leggett [87] and Nozeires and Schmitt-Rink [93]. Later
powerful semi analytic [91, 94, 95, 96, 97] and sophisticated numerical
schemes such as determinant quantum Monte Carlo (DQMC) [56, 57, 98, 99, 100,
101, 102, 103, 104, 105], dynamical mean field theory (DMFT) [106, 107, 108,
109, 110, 111] and static path approximation (SPA) [112] were utilized with
success to investigate the physics of the intermediate coupling regime of the
BCS-BEC crossover picture. This non Fermi liquid regime of the BCS-BEC
crossover is now known as the pseudogap phase and leaves distinct signatures
in the quasiparticle spectrum of the system. The quasiparticle spectra of the
pseudogap regime (T${}_{c}<$T$<$Tpg) is characterized by renormalized phase
coherence peaks arising out of the existing short range pair correlations in
this region. Further, at the Fermi level the quasiparticle spectra is no
longer hard gapped (unlike T $<$ Tc) and contains finite spectral weight due
to unpaired fermions arising out of thermal pair breaking.

The two thermal scales are thus defined as the Tc, at which the
superconducting state loses its long range macroscopic phase coherence and
Tpg, at which short range pair correlations are lost. The pseudogap regime
with short range pair correlations but no long range order, is characterized
by spatial fragmentation of the state into isolated islands akin to
uncorrelated Josephson junctions. The finite temperature behavior of the BCS-
BEC crossover shows that in the regime of intermediate (and strong) coupling
it is the thermally induced phase fluctuations which leads to the loss of long
range phase coherence and of the global superconductivity, unlike the weak
coupling regime where superconductivity is destroyed via the suppression of
pairing field amplitude. One thus needs to take into account the effect of
phase fluctuations in order to correctly capture the behavior of
superconductors or superfluids with interactions away from the weak coupling
regime, at finite temperatures.

#### BCS-BEC crossover and FFLO systems

Having established the importance of fluctuations in superconducting systems
we now once again take a look at the candidate systems where FFLO or BP phases
are reported. The key candidates as discussed earlier are, three dimensional
heavy fermion superconductor CeCoIn5, two-dimensional layered organic
superconductor $\kappa$-BEDT, both with large Maki parameter of $\alpha>1$ and
neutral ultracold atomic gases.

Apart from the large Maki parameter there is another property which serves as
a common link between these apparently very distinct systems. An important
characteristic of any superconductor is the ratio between the superconducting
gap at zero temperature and Tc. The BCS estimate of the same is
$2\Delta_{0}/k_{B}T_{c}\sim 3.5$ [1, 8], which corresponds to the weak
coupling regime of Figure 3. The BCS estimate sets the scale, based on which
the interaction regime of other superconducting systems are determined. Based
on the point contact [113] and scanning tunnelling spectroscopy [114] used to
determine the magnitude of the gap of any superconducting material, it has
been observed that for CeCoIn5 ($\alpha\sim 5$) $2\Delta_{0}/k_{B}T_{c}\sim
7.73$ [114], while for $\kappa$-BEDT ($\alpha\sim 8$) it is
$2\Delta_{0}/k_{B}T_{c}\sim 4.4$ [115]. The iron superconductor is conjectured
to be a two gap superconductor with the ratio for the smaller gap being
$2\Delta_{0}/k_{B}T_{c}\sim 1.4$ and that of the larger gap being
$2\Delta_{0}/k_{B}T_{c}\sim 4.6$ [116]. These experimental observations
indicate that all the candidate Pauli limited solid state superconductors are
situated away from the weak coupling regime. In case of the ultracold atomic
gases the interaction can be tuned externally and most of the experiments are
carried out at unitarity, which roughly corresponds to the regime of maximum
Tc in the BCS-BEC crossover, shown in Figure 3.

The experimental observations thus bring forth a common link between the
various Pauli limited superconducting systems, viz. all these systems belong
to the regime of intermediate coupling. While at one hand this sets up the
stage for complex physical properties to play out, on the other hand it calls
for a non perturbative approach to address these systems. It is thus
imperative to take into account the effect of fluctuations while addressing
the finite temperature behavior of the imbalanced Fermi superconductors and
superfluids. Consequently, a mean field approach to the problem will not
suffice.

We discussed above that apart from the semi analytic and T-matrix approaches
[91, 94, 95, 96, 97] there are powerful numerical techniques which takes into
account the effect of fluctuations in quantum many body systems [56, 98, 99,
100, 101, 102, 103, 104, 105, 57, 106, 107, 108, 109, 110, 111]. We now
discuss a path integral formalism based on which the different techniques and
their regimes of validity can be understood.

#### Path integral formalism for imbalanced Fermi systems

We begin by rewriting the two-dimensional attractive Hubbard model with
population imbalance, introduced in Equation (12), which can be used to model
the solid state as well as the ultracold atomic gas systems. We will discuss
about the system specific modifications in the later sections and for the sake
of simplicity discuss here the theoretical formalism in detail for a square
lattice with on-site attractive interaction $|U|$ which gives rise to $s$-wave
superconductivity. The corresponding Hamiltonian reads as,

$\displaystyle\hat{H}$ $\displaystyle=$ $\displaystyle\sum_{\langle
ij\rangle,\sigma}(t_{ij}-\mu\delta_{ij})\hat{c}_{i\sigma}^{\dagger}\hat{c}_{j\sigma}-h\sum_{i}\sigma_{iz}-|U|\sum_{i}\hat{n}_{i\uparrow}\hat{n}_{i\downarrow}$

where, $t_{ij}=-t$ only for the nearest neighbor hopping and is zero
otherwise; $\sigma_{iz}$ is the Pauli matrix. The hopping amplitude $t$ is set
to unity and is chosen to be the reference energy scale for the problem. $\mu$
is the net chemical potential and the Zeeman field $h$ is applied along the
$\hat{z}$ direction; $|U|>0$.

The Hubbard partition function is a functional integral over the Grassman
fields $\psi_{i\sigma}(\tau)$ and $\bar{\psi}_{i\sigma}(\tau)$,

$\displaystyle Z=\int{\cal D}\psi{\cal
D}{\bar{\psi}}e^{-\int^{\beta}_{0}d\tau{\cal L}(\tau)}$ (31)
$\displaystyle{\cal L}(\tau)={\cal L}_{0}(\tau)+{\cal L}_{U}(\tau)$ (32)
$\displaystyle{\cal L}_{0}(\tau)=\sum_{\langle
ij\rangle,\sigma}\\{\bar{\psi}_{i\sigma}((\partial_{\tau}-\mu_{\sigma})\delta_{ij}+t_{ij})\psi_{j\sigma}\\}$
(33) $\displaystyle{\cal
L}_{U}(\tau)=-|U|\sum_{i\sigma\sigma^{\prime}}\bar{\psi}_{i\sigma}\psi_{i\sigma}\bar{\psi}_{i\sigma^{\prime}}\psi_{i\sigma^{\prime}}$
(34)

where, $\beta$ is the inverse temperature, $\mu_{\sigma}$ takes into account
the field $h$. Our strategy is to decompose the interaction in terms of the
bosonic auxiliary $s$-wave pairing field
$\Delta_{i}(\tau)=|\Delta_{i}(\tau)|e^{i\theta_{i}(\tau)}$, using Hubbard-
Stratonovich transformation [117, 118]. Here $|\Delta_{i}(\tau)|$ and
$\theta_{i}(\tau)$ corresponds to the amplitude and phase of the pairing
field, respectively, and encodes spatial ($i$) and imaginary time ($\tau$)
dependence. This leads to,

$\displaystyle Z=\int{\cal D}\psi{\cal D}{\bar{\psi}}{\cal D}\Delta{\cal
D}{\Delta}^{*}e^{-\int^{\beta}_{0}d\tau{\cal L}(\tau)}$ (35)
$\displaystyle{\cal L}(\tau)={\cal L}_{0}(\tau)+{\cal L}_{U}(\tau)+{\cal
L}_{cl}(\tau)$ (36) $\displaystyle{\cal L}_{0}(\tau)=\sum_{\langle
ij\rangle,\sigma}\\{\bar{\psi}_{i\sigma}(\tau)((\partial_{\tau}-\mu_{\sigma})\delta_{ij}+t_{ij})\psi_{j\sigma}(\tau)\\}$
(37) $\displaystyle{\cal
L}_{U}(\tau)=\sum_{i}[\Delta_{i}(\tau)\bar{\psi}_{i\uparrow}(\tau)\bar{\psi}_{i\downarrow}(\tau)+h.c.]$
(38) $\displaystyle{\cal
L}_{cl}(\tau)=\sum_{i}\frac{|\Delta_{i}(\tau)|^{2}}{|U|}$ (39)

The $\psi$ integral is now quadratic, but an additional integration over the
field $\Delta_{i}(\tau)$ has been introduced. The weight factor for the
$\Delta$ configurations can be determined by integrating out the $\psi$,
$\bar{\psi}$ and using these weighted configurations, one goes back and
compute the fermionic properties. Formally,

$\displaystyle Z=\int{\cal D}\Delta{\cal
D}{\Delta}^{*}e^{-S_{eff}\\{\Delta,\Delta^{*}\\}}$ (40) $\displaystyle
S_{eff}=\ln Det[{\cal
G}^{-1}\\{\Delta,\Delta^{*}\\}]+\int^{\beta}_{0}d\tau{\cal L}_{cl}(\tau)$ (41)

where, ${\cal G}$ is the electronic Green’s function in the $\\{\Delta\\}$
background. The weight factor for an arbitrary space-time configuration
$\Delta_{i}(\tau)$ involves computation of the fermionic determinant in that
background. If we write the auxiliary field $\Delta_{i}(\tau)$ in terms of its
Matsubara modes as, $\Delta_{i}(\Omega_{n})$ then the different options are
readily recognized,

  * •

Quantum Monte Carlo retains the full “$i,\Omega_{n}$” dependence of $\Delta$
computing $\ln Det[{\cal G}^{-1}\\{\Delta\\}]$ iteratively for importance
sampling. The approach is valid at all T, but does not readily yield real
frequency spectra.

  * •

Mean field theory (MFT) is time independent, neglects the phase fluctuations
completely but can handle spatial inhomogeneity in amplitude of the pairing
field. Thus, $\Delta_{i}(i\Omega_{n})\rightarrow\mid\Delta_{i}\mid$. When the
mean field order parameter vanishes at high temperature the theory
trivializes.

  * •

Dynamical mean field theory (DMFT) retains the full dynamics but keeps
$\Delta$ at effectively one site, i.e,
$\Delta_{i}(\Omega_{n})\rightarrow\Delta(\Omega_{n})$. This is exact when
dimensionality $D\rightarrow\infty$.

  * •

Static path approximation (SPA) [119, 120, 121] approach retains the full
spatial dependence in $\Delta$ but keeps only the $\Omega_{n}=0$ mode, i.e,
$\Delta_{i}(\Omega_{n})\rightarrow\Delta_{i}$. It thus includes classical
fluctuations of arbitrary magnitude but no quantum ($\Omega_{n}\neq 0$)
fluctuations. One may consider different temperature regimes. (1) $T=0$: Since
classical fluctuations die off at $T=0$, SPA reduces to standard Bogoliubov-
de-Gennes (BdG) MFT. (2) At $T\neq 0$ we consider not just the saddle-point
configuration but all configurations following the weight $e^{-S_{eff}}$
above. These involve the classical amplitude and phase fluctuations of the
order parameter, and the BdG equations are solved in all these configurations
to compute the thermally averaged properties. This approach suppresses the
order much quicker than in MFT. (3) High T : Since the $\Omega_{n}=0$ mode
dominates the exact partition function, the SPA approach becomes exact as
$T\rightarrow\infty$. It is thus akin to the MFT only at $T=0$ but captures
the thermal behavior of the system accurately. Within this approach the system
can be envisaged as free fermions moving in a random background of
$\Delta_{i}$.

#### Static path approximation (SPA)

Evidently, all the numerical approaches have different set of merits and
caveats, and consequently different regimes of applicability. We now discuss
the numerical implementation of the SPA technique in some detail [122, 123,
124, 125, 126]. The merits of this technique rests on its applicability, in
principle to systems of arbitrary dimensions and lattice sizes, with moderate
computation cost and its access to real frequency properties of the system.
The technique takes into account spatial fluctuations and agreement with the
results of numerically exact DQMC suggests that it can capture the thermal
scales and high temperature phases accurately [112]. Using the Hubbard
Stratonovich decomposition to decompose the interaction term in Equation (29)
and then dropping the imaginary time dependence of the bosonic auxiliary field
($\Delta_{i}(\tau)$) the effective Hamiltonian reads as,

$\displaystyle H_{eff}$ $\displaystyle=$ $\displaystyle-t\sum_{\langle
ij\rangle,\sigma}\hat{c}_{i\sigma}^{\dagger}\hat{c}_{j\sigma}+\sum_{i}\Delta_{i}\hat{c}_{i\uparrow}^{\dagger}\hat{c}_{i\downarrow}^{\dagger}+h.c.-\mu\sum_{i,\sigma}\hat{n}_{i\sigma}$
(42)
$\displaystyle-h\sum_{i}\sigma_{i}^{z}+\sum_{i}\frac{\mid\Delta_{i}\mid^{2}}{\mid
U\mid}$

where, the last term $H_{cl}=\frac{|\Delta_{i}|^{2}}{|U|}$ is the stiffness
cost associated with the now “classical” auxiliary field. In the standard mean
field approach the superconducting pairing field $\Delta_{i}$ is assumed to be
a real number, but here we retain the degrees of freedom associated with the
pairing field phase and amplitude.

The background field $\Delta_{i}$ obeys the Boltzmann distribution,

$\displaystyle P\\{\Delta_{i}\\}\propto Tr_{cc^{\dagger}}e^{-\beta H_{eff}}$
(43)

which is connected to the free energy of the system. For large and random
$\Delta_{i}$ the trace is taken numerically. We generate the random background
of $\\{\Delta_{i}\\}$ by using Monte Carlo, diagonalizing $H_{eff}$ for each
attempted update of $\Delta_{i}$. The relevant fermionic correlators are
computed on the optimized configurations at different temperatures.

Evidently, the technique is numerically expensive and involves large
computation cost in diagonalizing the matrix to compute the free energy for
each attempted update of the auxiliary field. For an $N\times N$ matrix the
computation cost scales as $\sim{\cal O}(N^{3})$ per update and $\sim{\cal
O}(N^{4})$ per lattice sweep. This imposes constraint on the lattice size that
can be accessed within reasonable computation time. In order to access larger
system sizes a cluster based update scheme could be utilized wherein instead
of diagonalizing for the entire lattice at each attempted update, a smaller
cluster of size $N_{c}=L_{c}\times L_{c}$ surrounding the update site is
diagonalized [123, 124, 126]. This reduces the cost per update to $\sim{\cal
O}(N_{c}^{3})$ and to $\sim N\times{\cal O}(N_{c}^{3})$ per lattice sweep,
which scales linearly with the system size. This allows one to access system
sizes of upto $\sim 40\times 40$, which is beyond the reach of DQMC. Access to
larger system sizes aid in to explore lower interaction regimes where the
coherence length of the Cooper pairs are large. Moreover, in order to
faithfully capture spatially inhomogeneous orders, access to lager lattice
sizes is essential.

We have now set up a theoretical framework to investigate the behavior of
spatially inhomogeneous superconducting and superfluid phases. Different
numerical tools to access these systems both at the level of ground state as
well as at finite temperature through mean field and non perturbative
techniques are discussed. In order to understand the theoretical formalism we
considered that the interaction between the fermion species is onsite and the
resulting pairing state symmetry is $s$-wave. This is a valid approximation as
far as ultracold atomic gas set ups are concerned and most of the experiments
carried out on imbalanced Fermi systems in ultracold atomic gases are for
$s$-wave superfluids. Furthermore, here a single channel decomposition of the
interaction term was carried out. In principle, along with the pairing channel
the interaction can be decomposed in density and spin channels as well. This
however would add to the computation cost and more importantly would not lead
to any qualitative difference in the behavior of the present system.

The behavior of the imbalanced Fermi superfluid (or superconductor) is
analyzed based on different thermodynamic and quasiparticle indicators, some
of which are,

  * •

Magnetization:- $m_{i}=n_{i\uparrow}-n_{i\downarrow}$.

  * •

Pairing field structure factor:- $S_{\Delta}({\bf
q})=\frac{1}{N^{2}}\sum_{i,j}\langle\Delta_{i}\Delta_{j}^{*}\rangle e^{i{\bf
q}.({\bf r}_{i}-{\bf r}_{j})}$

  * •

Spin resolved single particle density of states (DOS):-
$N_{\uparrow}(\omega)=\frac{1}{N}\langle\sum_{i}|u_{n\uparrow}^{i}|^{2}\delta(\omega-\epsilon_{n})\rangle$,
$N_{\downarrow}(\omega)=\frac{1}{N}\langle\sum_{i}|v_{n\downarrow}^{i}|^{2}\delta(\omega-\epsilon_{n})\rangle$

  * •

Spectral functions and lineshapes:- $A_{\sigma}({\bf
k},\omega)=(-1/\pi)ImG_{\sigma}({\bf k},\omega)$.

  * •

Low energy spectral weight distribution:- $A({\bf k},0)$.

Here, $G({\bf k},\omega)=lim_{\delta\rightarrow 0}G({\bf
k},i\omega_{n})|_{i\omega_{n}\rightarrow\omega+i\delta}$, with $G({\bf
k},\omega)$ being the imaginary frequency transform of $\langle c_{\bf
k}(\tau)c_{\bf k}^{\dagger}(0)\rangle$.

Note that the static path approximation approach, as mentioned above does not
take into account the quantum fluctuations at the ground state. While in
continuum this could be a poor approximation, on a lattice it is reasonable.
The low energy fluctuations in a FFLO like state arises from, (i) the phase
symmetry of the $U(1)$ order parameter, (ii) the translational symmetry
breaking and (iii) rotational symmetry breaking. As a result, in continuum a
long range order can not survive even at T=0 in two dimensions. The mean field
theory, which predicts such order is therefore invalid.

On a lattice, the spatial symmetry is already broken and thus the
translational and rotational modes are gapped out. The only low energy
excitations thus left are the XY type, in presence of which the two-
dimensional system can have long range order at T=0 and undergoes Berezinskii-
Kosterlitz-Thousless (BKT) transition at the finite temperatures. The issue of
fluctuations thus reduces to how well the $U(1)$ phase fluctuations and the
corresponding Tc are captured by this approach.

#### Non $s$-wave pairing

The superconducting pairing state symmetry in solid state materials depend on
the material parameters and the coupling constants involved. As discussed
earlier, the all the known Pauli limited superconductors belong to the regime
of intermediate coupling. This in turn naturally categorises them to be
superconductors with short coherence length, whose pairing state symmetry can
no longer be $s$-wave. Indeed the pairing state symmetry of CeCoIn5 and
$\kappa$-BEDT are established to be $d$-wave [43, 51]. Unlike $s$-wave, the
$d$-wave pairing state requires an intersite interaction between the fermions.

The Hubbard Hamiltonian can be suitably modified to take into account the
nearest neighbor interaction which gives rise to $d$-wave pairing state
symmetry as,

$\displaystyle\hat{H}$ $\displaystyle=$ $\displaystyle\sum_{\langle
ij\rangle,\sigma}t_{ij}\hat{c}^{\dagger}_{i,\sigma}\hat{c}_{j,\sigma}-U\sum_{\langle
ij\rangle}\hat{n}_{i}\hat{n}_{j}-\mu\sum_{i,\sigma}\hat{n}_{i,\sigma}-h\sum_{i}\sigma_{i}^{z}$

where, the attractive interaction $|U|>0$ is now between the nearest
neighbors, on a square lattice. After carrying out the Hubbard-Stratonovich
transformation to decompose the four-fermion interaction, the corresponding
partition function reads as,

$\displaystyle Z=\int{\cal D}\psi{\cal D}{\bar{\psi}}{\cal D}\Delta{\cal
D}{\Delta}^{*}e^{-\int^{\beta}_{0}d\tau{\cal L}(\tau)}$ (45)
$\displaystyle{\cal L}(\tau)={\cal L}_{0}(\tau)+{\cal L}_{U}(\tau)+{\cal
L}_{cl}(\tau)$ (46) $\displaystyle{\cal L}_{0}(\tau)=\sum_{\langle
ij\rangle,\sigma}\\{\bar{\psi}_{i\sigma}((\partial_{\tau}-\mu_{\sigma})\delta_{ij}+t_{ij})\psi_{j\sigma}\\}$
(47) $\displaystyle{\cal L}_{U}(\tau)=-\sum_{i\neq
j}\Delta_{ij}(\bar{\psi}_{i\uparrow}\bar{\psi}_{j\downarrow}+\bar{\psi}_{j\uparrow}\bar{\psi}_{i\downarrow})+h.c.$
(48) $\displaystyle{\cal L}_{cl}(\tau)=4\sum_{i\neq
j}\frac{\mid\Delta_{ij}\mid^{2}}{\mid U\mid}$ (49)

The bosonic auxiliary field $\Delta_{ij}$ encodes the direction dependence of
the $d$-wave superconducting order parameter.

Apart from the superconducting pairing state symmetry there is another aspect
which needs to be taken into account while constructing a model Hamiltonian
for the solid state Pauli limited superconductors, viz the lattice structure.
The heavy fermion CeCoIn5 comprises of a three dimensional cubic lattice
structure [51], while on the other hand the quasi two-dimensional
$\kappa$-BEDT is a material with triangular lattice structure [115]. KFe2As2
contains an added level of complexity which needs to be taken into account
while setting up a theoretical framework, this material is conjectured to be a
multiband superconductor with two superconducting gaps (and thus the pairing
fields) of different magnitudes [116].

We conclude this section by stating that a realistic theoretical model of
solid state Pauli limited superconductors would be significantly complicated
owing to the various interaction parameters and complex band structure that
needs to be taken into account. However, a relatively simplistic model which
takes into account the dominant coupling scales and energy bands is amenable.
Such a theoretical framework should be able to capture the experimental
observations on these materials with qualitative and in some cases with
quantitative accuracy. It should however be kept in mind that solving even a
minimal quantum many body Hamiltonian is a demanding and computationally
extensive task. Consequently, there still remains several open areas to be
explored in the playground of Pauli limited superconductors as well as of
Fermi superfluids in ultracold atomic gas setups.

##  0.4 Mean field and beyond mean field studies

This section discusses some of the studies carried out on lattice fermion
models to understand population imbalanced Fermi gases, Pauli limited
superconductors and mass imbalanced Fermi-Fermi mixtures. Most of the T=0
investigations are carried out within the framework of mean field theory and
there is a large volume of literature discussing the same. We do not attempt
to be comprehensive and would discuss only some of the important results and
the inferences drawn based on them.

Regarding the beyond mean field attempts, broadly the observations made by
three different numerical approaches are covered in this section, viz. (i)
dynamical mean field theory (DMFT) [127, 128, 129], (ii) determinant quantum
Monte Carlo (DQMC) [130, 131] and (iii) static path approximation (SPA) [123,
125].

The thermal transitions pertaining to two-dimensional systems, discussed in
this section are Berezinskii-Kosterlitz-Thouless (BKT) transitions
corresponding to algebraic decay of quasi long range order in two dimensions
[132, 133, 134].

###  0.4.1 Population imbalanced $s$-wave superfluids

Within the framework of lattice mean field theory the two-component Fermi
mixture with population imbalance in an optical lattice was studied by Koponen
et al. [135]. The authors found that both the FF and BP phases correspond to
stable ground state of the system, with the FF state corresponding to the
global minimum in a fixed chemical potential calculation, and the BP phase
corresponding to the minimum energy state for a fixed fermionic number density
calculation. The same authors further investigated the finite temperature
behavior of this system, wherein they accessed the finite temperature physics
by using the mean field theory. The resulting temperature-magnetization phase
diagram showed that at fixed number density the FF phase constitutes a large
part of the phase diagram along with the BP and phase separated regimes [136].

Using an elaborate variational mean field scheme based on BdG method with
Hartree corrections, Loh and Trivedi [137] investigated the LO phase in a
cubic lattice. Their mean field scheme comprised of six variational parameters
corresponding to the complex valued pairing field, chemical potential and
three Zeeman fields. The study found that for a cubic lattice the LO phase is
a stable ground state over a significant parameter regime. The authors further
discussed the quasiparticle behavior of the LO state and showed that the state
is characterized by additional van-Hove singularities in the single particle
density of states (DOS). Based on the momentum resolved fermionic occupation
number ($n_{\sigma}({\bf k})=\langle\hat{c}_{{\bf
k},\sigma}^{\dagger}\hat{c}_{{\bf k},\sigma}\rangle$) they inferred that the
LO phase breaks the four fold symmetry but preserves the inversion symmetry,
unlike the FF phase.

Figure 4: Color online: Comparison of Tc scales obtained from the mean-field
calculation (upper curve) and static path approximation (SPA) Monte Carlo
technique (lower curve). The solid and dashed curves in the mean field result
corresponds to second and first order phase transitions, respectively. In the
SPA results BP-II represents a breached pair state that undergoes a second-
order transition to the partially polarized Fermi liquid (PPFL), while BP-I
undergoes a first-order thermal transition to the PPFL. Beyond BP-I the system
exhibits FFLO order up to a critical field. Figure adapted with permission
from [123], Copyright (2016) by the American Physical Society.

An extensive numerical study of the ground state was carried out by Chiesa and
Zhang as well as by Rosenberg et al on a two-dimensional square lattice, using
Hartree-Fock-Bogoliubov-de-Gennes (HFBDG) theory [138, 139]. It was observed
that an uniaxial LO phase covered a large part of the phase diagram in the
magnetization-number density ($m-n$) plane, at an intermediate interaction
strength. In the regime of small magnetization and large number density,
diagonal stripe orders were found for intermediate interaction, corresponding
to a “supersolid” phase, in which charge and pairing orders coexisted. Apart
from these, mean field based comparative study with LO and FF ansatz showed
that for both two and three dimensional systems the LO state correspond to
lower energy phase as compared to the FF [82]. The observation was in
agreement with the preceding attempts along this direction [137, 140, 141].

Mean field investigations are carried out in one dimensional systems as well,
for example by Sun and Bolech, who studied the system of coupled one
dimensional tubes [142]. It was found that the stability of the FFLO state was
dependent on the magnitude of pair tunnelling between the tubes. In addition
to the FFLO phase, an exotic $\eta$-phase was realized in the mean field study
of this one dimensional system. In the $\eta$-phase the FFLO wave vector was
found to occur at the edge of the first Brillouin zone.

The behavior of the LO state across dimensional crossover between one to three
dimensions was investigated using DMFT by Törmä and co-authors [128, 129]. It
was shown that at the ground state the LO phase is stable over large regime of
parameters, extending well into the 3D-like regime, unlike the predictions of
mean field [143] and effective field theory [144]. When extended to finite
temperatures by using QMC solver in single site DMFT impurity problem, the
same system showed that the transition temperature of the LO phases (TFFLO) is
suppressed to upto one third of that of uniform superfluid at zero imbalance
[128]. Effect of non local fluctuations were taken into account in this
framework by using cluster-DMFT and it was demonstrated that CDMFT correctly
captured the behavior of the system in terms of decreasing robustness of long
range phase coherence with reducing dimensionality of the system [129].

In their DQMC study of imbalanced $s$-wave superfluid Wolak et al. [131]
observed signatures of finite momentum pairing, suggesting instability towards
the FFLO phase. At high temperature indications of a magnetized paired phase
was observed. The work however failed to explore the low temperature regime (T
$<$ 0.1t, with $t$ being the hopping scale), owing to the fermionic sign
problem of DQMC. Note that the access to low temperature regime is important
in case of FFLO phase, for which the transition temperature is expected to be
strongly suppressed, due to the population imbalance. This work therefore
failed to give an estimate of the Tc scale associated with the FFLO phase. The
authors predicted the zero imbalance Tc to be $\sim 0.1$t but suggested that
FFLO correlations survive up to a scale T $>>$ Tc.

Figure 5: Color online: Thermal evolution of the various indicators in the BP
phase. Starting from the top row, we show maps of superconducting amplitude
($|\Delta_{i}|$), phase correlation [($\cos{\theta_{0}-\theta_{i}}$),
$\theta_{0}$ is a reference site], pairing structure factor ($S_{\Delta}({\bf
q})$), magnetization $m_{i}$, magnetic structure factor ($S_{m}({\bf q})$),
and fermionic number density ($n_{i}$). The temperature, along the rows, is
marked at the bottom of the figure. Reprinted figure with permission from
[123], Copyright (2016) by the American Physical Society.

Gukelberger et al. [130] used diagrammatic QMC to map out the thermal phase
diagram of imbalanced $s$-wave superfluid in the magnetization-temperature
plane at intermediate lattice filling, based on their observation of the
divergence of the pair susceptibility. Interestingly, they found that the weak
magnetization low temperature regime of the phase diagram comprises of a zero
momentum pairing state with finite magnetization. As discussed in the earlier
sections, such a coexistence of superfluidity and magnetization corresponds to
the BP phase. The intermediate magnetization regime was found to be a FFLO
phase, followed by the high magnetization regime corresponding to a Fermi
liquid. The lowest temperature that was probed in this study was T$\sim 0.02$t
which in principle can be higher than the TFFLO. Thus, it remains unclear
whether the BP phase shown in the phase diagram by Gukelberker et al. [130]
corresponds to a true ground state or a finite temperature phase. It must
however be noted that earlier mean field studies of similar system, carried
out at fixed lattice filling suggested the BP phase to be a possible stable
ground state in the regime of weak magnetization [82].

A comprehensive analysis of this system in grand canonical ensemble (fixed
chemical potential) was carried out using SPA, close to the half filling and
at an intermediate interaction [123, 125]. Apart from mapping out the ground
state and finite temperature phase diagrams and capturing the thermal scales
accurately, the work discussed real space and quasiparticle signatures of the
LO phase as well as their thermal evolution. The ground state of the system
was found to comprise of uniform superfluid (USF), LO and partially polarized
Fermi liquid (PPFL) phases, as the function of increasing magnetic field. It
was demonstrated, as in Figure 4, that inclusion of thermal fluctuations
significantly suppresses the Tc corresponding to both uniform and LO
superfluid phases. Comparison with the results of mean field calculations
showed that the mean field theory overestimates the Tc by a factor of $\sim
4$. In the LO phase the Tc was about $20\%$ of its mean field estimate and in
the absolute term only $2\%$ of the hopping scale. The Tc corresponding to the
two component balanced Fermi gas in two dimensions is estimated to be
T${}_{c}\sim$ 0.15t, based on DQMC [57] as well as SPA [112, 123]
calculations. The maximum Tc of the LO phase is thus $\sim$0.03t, which is far
lower than the lowest temperature accessed by Wolak et al. [131] and very
close to that accessed by Gukelberger et al. [130].

Figure 6: Color online: Spatial maps characterizing the thermal evolution of
the LO state, through (a) pairing field amplitude, (b) phase correlation, (c)
magnetization and (d) fermionic number density. Reprinted figure with
permission from [123], Copyright (2016) by the American Physical Society.

The work showed that in the regime of weak magnetic field the system loses
uniform superfluid order to PPFL via a second order thermal transition, while
for the strong magnetic field LO phase, the thermal transition is of first
order. Moreover, a first order thermal transition regime was identified near
the quantum critical point between the uniform superfluid and FFLO phases. In
this regime, the uniform (${\bf Q}$=0) superfluid phase undergoes first order
thermal transition to the PPFL phase. Over the entire regime of magnetic field
(magnetization) signatures of superfluid pairing were observed at the high
temperatures as short range pairing correlations, which survives upto
temperatures T $\sim$ 1.5Tc and gives rise to the pseudogap phase. Akin to the
Tc the pseudogap scale is strongly suppressed by the magnetic field.

Based on the real space maps the low magnetic field high temperature phase was
identified to be a BP phase with spatial coexistence of zero momentum
superfluid correlations and non zero magnetization, as shown in Figure 5
[123]. In Figure 5, along with the pairing field amplitude ($|\Delta_{i}|$),
phase coherence ($\cos(\theta_{0}-\theta_{i})$) and magnetization ($m_{i}$),
the momentum resolved pairing field ($S_{\Delta}({\bf q})$) and magnetic
($S_{m}({\bf q})$) structure factors, as well as spatially resolved fermionic
number density ($n_{i}$) are shown. The lowest temperature state corresponds
to a phase cohered uniform superfluid with vanishingly small magnetization.
Increase in temperature leads to suppression of long range phase coherence and
fragmentation of the superfluid state. Concomitantly, magnetization begins to
develop spatially in the system, in the form of isolated islands, which are
roughly complementary to the suppression of the pairing field. The real space
observation is further attested by the momentum resolved structure factor
corresponding to the pairing field and magnetization.

Figure 7: Color online: Thermal evolution of the momentum and spin resolved
fermionic occupation number $n_{\sigma}({\bf k})$ in the LO phase. Reprinted
figure with permission from [123], Copyright (2016) by the American Physical
Society.

At the ground state the LO phase was characterized by uniaxial stripes with
the exact pairing momentum being dictated by the magnitude of the Zeeman field
for a fixed interaction and chemical potential. The thermal disordering of the
LO phase is demonstrated via the real space maps in Figure 6 [123]. The maps
show that the thermal scales are strongly suppressed in the LO phase and the
long range order is lost even at T$\sim$0.005t, even though short range
correlations survive upto high temperatures. The superfluid order and
magnetization are spatially complementary to each other and the minimum of one
corresponds to the maximum of the other.

Distribution of spin and momentum resolved fermionic occupation number which
maps out the Fermi surface shows that in the LO phase the Fermi surface is
two-fold symmetric, rather than the four-fold symmetry of the uniform
superconducting state. Signatures of LO phase was observed upto T $\sim$ 2Tc,
at still higher temperatures the signatures of pairing correlations are lost
and the distribution was found to be akin to the PPFL phase, as shown in
Figure 7.

Figure 8: Color online: Thermal evolution of spectral function
($A_{\uparrow}({\bf k},\omega)$) (top row) and low energy spectral weight
distribution ($A({\bf k},0)$) (bottom row), in the LO phase. The spectral
function is plotted along the $(0,0)-(\pi,\pi)$ trajectory across the
Brillouin zone. The $x$-axis corresponds to the momentum $k_{x}=k_{y}$ and the
$y$-axis corresponds to the energy eigenvalues. The low energy spectral weight
distribution in shown in the $k_{x}-k_{y}$-plane.

The quasiparticle properties such as single particle DOS, spectral functions
and distribution of low energy spectral weight carried intriguing imprints of
the underlying modulated superfluid phase [123, 125]. The single particle DOS
contained additional van-Hove singularities, in agreement with the mean field
predictions [137]. This was the first and so far the only work on population
imbalanced $s$-wave superfluids which mapped out the spectral function
($A({\bf k},\omega)$) of the LO phase, within a lattice fermion model. Owing
to the multiple scatterings that the fermion undergoes, the dispersion spectra
of the LO phase was found to contain six branches, in contrast to the two
branches of the BCS (uniform) superconductors. Each fermion species gives rise
to three dispersive branches, separated by soft gaps, located at ($\omega=0$)
and away ($\omega=\pm h$), from the Fermi level. In agreement with such a non
trivial dispersion spectra the corresponding low energy spectral weight
distribution ($A({\bf k},0)$) was found to be highly intriguing. It was
observed that in spite of a $s$-wave pairing state symmetry a nodal superfluid
gap structure emerges out of the LO correlations. Figure 8 shows the thermal
evolution of the spectral function (for one of the fermionic species) and low
energy spectral weight distribution of the LO phase. Thermal fluctuations lead
to the smearing out of the nodal architecture of the low energy spectral
weight distribution and at high temperatures the four-fold symmetry of the
$s$-wave superfluid gap is restored.

Figure 9: Color online: The imbalance-temperature phase diagram inferred from
measurements on the ${}^{6}Li$ cold atomic gas at unitarity (left) [16],
compared with the results obtained at an intermediate coupling (peak Tc) of 2D
Hubbard model, using SPA (right). The normalization of the x axis is same in
both panels, while the y axis have different reference scales. Figure adapted
with permission from [123], Copyright (2016) by the American Physical Society.

This work drew comparison between the numerically computed [123] thermal phase
diagram with that obtained experimentally in ultracold atomic gas set up [16].
As shown in Figure 9, the agreement between the two, both in terms of the
order of phase transitions as well as thermal scales, is fairly good.
Numerical calculations suggest that the large magnetization regime hosts the
LO phase. Though the signature of spatially modulated LO phase was not
observed in the experiment, its possibility has not been ruled out either. In
a related work the same authors determined an anomalous pseudogap phase with
non trivial spectral signatures, near the finite temperature BP-LO boundary at
low magnetizations and a re-entrant pseudogap phase at large magnetizations
[125].

###  0.4.2 Pauli limited $d$-wave superconductors

Studies carried out on the FFLO phase of non $s$-wave superconductors are
relatively fewer, inspite of all the solid state Pauli limited superconductors
being $d$-wave. Vorontsov et al. used self consistent theory to study the
effect of in-plane Zeeman field on a two-dimensional superconducting system
with $d_{x^{2}-y^{2}}$ pairing state symmetry [145, 146]. They solved
Eilenberger equations in presence of magnetic field for the superconducting
order parameter with spatial dependence, to map out the phase diagram in the
temperature-magnetic field (T-H) plane. Their study showed the LO phase to be
more stable as compared to the FF, over a large part of the parameter space.
The single particle DOS as determined from their calculations showed van-Hove
singularities akin to Andreev bound states, which are characteristic to FFLO
state.

Figure 10: Color online: (Left) Magnetization-temperature (m-T) and (Right)
Zeeman field-temperature (h-T) phase diagram of population imbalanced $d$-wave
superconductor. The solid curves correspond to Tc, with second and first order
phase transitions being represented by black and red curves, respectively. The
h-T phase diagram shows the thermal scale Tpg in terms of dotted curve and the
associated pseudogap regime. Reprinted figure with permission from [126],
Copyright (2020) IOP Publishing Ltd. Printed in the UK.

Within the framework of Ginzburg-Landau theory Beaird et al. studied Pauli
limited superconductors with $d$-wave pairing state symmetry [147]. They
included thermally induced classical magnetic fluctuations in their
theoretical framework to explain the experimentally observed phase diagram of
CeCoIn5. Their work showed that the inclusion of magnetic fluctuations affect
the thermodynamic phase diagram significantly such that, a regime of first
order thermal transition from the normal state emerges near the quantum
critical point where the system undergoes transition between the uniform
$d$-wave superconducting state and a spatially modulated LO phase. Beaird et
al. [147] showed that the regime of first order thermal transition straddles
both the uniform and LO superconducting phases in agreement with the
experimental observations of the heavy fermion Pauli limited superconductor
CeCoIn5 [44].

The effect of non magnetic impurities on Pauli limited $d$-wave superconductor
was studied by Vorontsov et al. using quasiclassical Eilenberger equations in
which the effect of impurities was included via self consistent T-matrix
approximation [148]. It was inferred that contrary to the naive expectations,
the FFLO phase is robust against impurities in the limit of strong disorder.

Figure 11: Color online: Real space maps at T=0 showing the evolution of the
QBP phase over a range of Zeeman field in terms of (i) the amplitude of
superconducting pairing field, (ii) phase correlation of the superconducting
pairing field and (iii) the magnetization. The color signifies the magnitude
of the corresponding indicator. Reprinted figure with permission from [126],
Copyright (2020) IOP Publishing Ltd. Printed in the UK.

Using the BCS weak coupling model Yang and Sondhi made certain intriguing
predictions about $d$-wave superconductor subjected to an in-plane Zeeman
field [149]. They showed that owing to the nodal architecture of the $d$-wave
pairing state the Zeeman field would destroy superconductivity in parts of the
Fermi surface and give rise to pockets of normal electrons. The authors
suggested that this phenomenon should be observable in the quasiparticle
signatures of the system.

Using the BdG mean field theory on a two-dimensional lattice model Zhou and
Ting mapped out the thermal phase diagram of Pauli limited $d$-wave
superconductor as a function of in-plane Zeeman field [150]. Apart from the
uniform superconducting phase in the low field regime they observed one and
two-dimensionally modulated FFLO phases in the regime of high field and low
temperatures. Since the approach does not take into account the effects of
thermal fluctuations the stability of the modulated phases remain unknown. The
single particle DOS at the ground state as obtained by the BdG mean field
theory showed signatures of LO state in terms of Andreev bound states.

The beyond mean field studies of $d$-wave Pauli limited superconductors are
much lesser than that of $s$-wave. In a series of papers Ikeda and co-authors
studied the problem of FFLO phase in $d$-wave superconductor in presence of
antiferromagnetic (AFM) fluctuations, relevant for materials such as CeCoIn5
[151, 152, 153, 154, 155]. Within the framework of repulsive Hubbard model
which gives rise to $d$-wave superconductivity as well as AFM order, the
effect of paramagnetic pair breaking was taken into account through an in-
plane Zeeman field. Using fluctuation exchange approximation (FLEX) the
authors showed that the high magnetic field low temperature regime comprises
of coexisting AFM and $d$-wave superconducting orders, as well as a FFLO
phase. The authors claimed their results to be in agreement with the
experimental observations on CeCoIn5. A Ginzburg Landau theory based study was
carried out by Adachi and Ikeda on two-band Pauli limited superconductors,
which is expected to provide useful insights to the behavior of Pauli limited
iron superconductors such as KFe2As2 [156].

Figure 12: Color online: Real space maps as the system evolves in temperature,
at a selected Zeeman field in the QBP phase. Reprinted figure with permission
from [126], Copyright (2020) IOP Publishing Ltd. Printed in the UK.

The static path approximation (SPA) technique was used at a fixed chemical
potential and intermediate interactions to address the behavior of two-
dimensional population imbalanced $d$-wave superconductor, in detail [126].
The corresponding phase diagram as shown in Figure 10 was mapped out in the
magnetization-temperature (m-T) as well as in the field-temperature (h-T)
plane. It was demonstrated that the low field (magnetization) regime comprises
of a “quantum breached pair” (QBP) phase, characterized by spatial coexistence
of “gapless” $d$-wave superconductivity and non zero magnetization. This was
the first numerical realization of the coexistent phase predicted by Yang and
Sondhi, in a lattice fermion model [149]. The QBP phase can be understood as
the T=0 counterpart of the finite temperature BP phase realized in the
population imbalanced $s$-wave superconductors/superfluids. It was observed
that increase in field (magnetization) leads to a Lifshitz transition from the
uniform $d$-wave superconductor to the QBP phase, quantified by the average
magnetization of the system. Both the uniform $d$-wave superconductor and QBP
phases undergo second order thermal phase transition to the PPFL phase. The
high magnetic field (magnetization) regime hosts FFLO state at low
temperatures, which undergoes first order thermal phase transition to PPFL.
Across the quantum critical point corresponding to the transition between the
QBP and FFLO phases the thermal phase transition was found to be of first
order, in agreement with the earlier predictions [147]. Inclusion of
fluctuations showed significant suppression of the superconducting transition
temperature, both for the uniform superconducting as well as FFLO regimes.

Based on the real space signatures shown in Figure 11 and Figure 12 it, was
inferred that both field and temperature leads to spatial inhomogeneity in the
superconducting state, which shows up in the quasiparticle signatures as the
“pseudogap” behavior. The phase diagram thus comprises of an additional
thermal scale Tpg as shown in Figure 10, where short range superconducting
pair correlations die out, long after the loss of long range superconducting
phase coherence at Tc [126].

###  0.4.3 Mass imbalanced $s$-wave superfluids

In ultracold atomic gases pairing between equal populations of two atomic
species with unequal atomic masses (Fermi-Fermi mixture) can give rise to mass
imbalanced Fermi superfluid. Experimentally a mass imbalanced Fermi-Fermi
mixture can be achieved in 6Li-40K mixtures. The theoretical studies are
primarily restricted to continuum models, using both mean field and beyond
mean field techniques [157, 158, 159, 160, 161, 162, 163]. The cumulative
outcome of these works brought forth two important observations, (i) imbalance
in mass promotes the instability of the system towards a supersolid phase
accompanied by a Lifshitz transition, (ii) the mass imbalanced Fermi-Fermi
mixture contains more than one pseudogap scales at high temperature.

Works on lattice fermion models of mass imbalanced Fermi-Fermi mixtures are
very limited. Wang et al. [164] studied this problem in one dimension using
time-evolving block decimation technique and found that increasing mass
imbalance ratio shrinks the FFLO regime of the phase diagram and at large
imbalance the FFLO phase disappears.

Figure 13: Color online: Ground-state dispersion spectra at different mass
imbalance ratio $\eta$. The last two panels show the $\eta$ dependence of the
density of states (DOS) for the (a) light- and (b) heavy-fermion species.
Reprinted figure with permission from [124], Copyright (2018) by the American
Physical Society.

Using density matrix renormalization group (DMRG) and quantum Monte Carlo
(QMC) techniques Dalmonte et al. [165] showed that mass imbalance favors
paired phases as compared to population imbalance, and concluded that in mass
imbalanced systems it might be possible to detect FFLO in experimentally
realizable parameter regimes. A mean field study on two-dimensional square
lattice by Pahl and Koinov revealed that the ground state phase diagram of
this system comprises of normal, FF and BP regimes and that a majority of
heavy species promotes the FF phase [166]. A non perturbative lattice Monte
Carlo study was carried out in two-dimensions, which revealed that a mean
field approach to the problem grossly overestimates the ground state energy
and therefore makes the stability of the phases questionable [167].

The only beyond mean field study of mass imbalanced Fermi-Fermi mixture in
two-dimensions, at and above the ground state was carried out using the SPA
[124]. Within the framework of attractive Hubbard model the imbalance in mass
was taken into account through imbalance in the hopping amplitude of the two
species as, $t_{H}/t_{L}=m_{L}/m_{H}=\eta$. Here, $H$ and $L$ corresponds to
the heavy and light fermion species, respectively and $\eta$ serves as the
tuning parameter to control the ratio of imbalance. The work demonstrated that
imbalance in mass significantly alters the dispersion spectra of the fermionic
species and gives rise to subgap and supergap features in the species resolved
single particle DOS, as shown in Figure 13.

Figure 14: Color online: Mass imbalance-temperature ($\eta$-T) phase diagram
at fixed population imbalance. Along with the superfluid regime, the figure
shows the pseudogap regimes based on the species-resolved DOS. T${}_{max}^{L}$
and T${}_{max}^{H}$ correspond to the pseudogap scales. PG-I corresponds to
the regime where both the fermionic species are pseudogapped and in the PG-II
regime only the lighter species is pseudogapped. Reprinted figure with
permission from [124], Copyright (2018) by the American Physical Society.

At a fixed imbalance of population, the phase diagram comprises of LO, BP and
unstable (phase separated) regimes in the imbalance-temperature plane and a
larger imbalance in mass reduces the regime of stability of these phases. The
thermal phase diagram in the imbalance-temperature ($\eta-T$)-plane presented
in Figure 14 shows the intriguing feature of two pseudogap phases as PG-I and
PG-II, in agreement with the observation of continuum models [158, 161]. Owing
to the mismatch in mass and thus of the hopping amplitudes the two species are
subjected different “scaled” temperatures. Consequently, while in the
pseudogap regime PG-I both the fermion species are pseudogapped, in the
pseudogap regime PG-II it is only the lighter species which is pseudogapped
while the heavier species is a polarized Fermi liquid. At $\eta=1$
corresponding to zero imbalance in mass, both the pseudogap scales collapses
to one. The phase diagram shows that, while akin to the population imbalanced
system the Tc is strongly suppressed by imbalance in mass, the short range
pair correlations survive over a large temperature regime of T $\sim$ 2Tc
[124].

Figure 15: Color online: Interaction-temperature (U-T) phase diagram at
$\eta$=0.15 corresponding to the ${}^{6}Li-{}^{40}K$ mixture. The figure shows
the Tc scale along with the pseudogap scales for this mixture at a fixed
population imbalance. Reprinted figure with permission from [124], Copyright
(2018) by the American Physical Society.

The work further discussed the BCS-BEC crossover picture for mass imbalanced
Fermi-Fermi mixture at a fixed population imbalance and established the
presence of two pseudogap regimes across the crossover. The BCS-BEC crossover
revealed a crucial observation. It showed that a critical interaction strength
(Uc) is required to realize superfluid order in mass imbalanced Fermi-Fermi
mixture in presence of a population imbalance. This observation is in striking
contrast with those of two component balanced Fermi gases where superfluid
phase is realized for any arbitrary strength of attractive interaction. Figure
15 presents the BCS-BEC crossover of the Fermi-Fermi mixture, which shows that
the superfluid Tc collapses to zero at U$\sim 3t_{L}$ for the parameter regime
mentioned in the figure caption.

A similar study was carried using SPA by Karmakar [126], for Fermi-Fermi
mixtures with non local $d$-wave interaction between the fermion species. In
the context of solid state systems, materials with unequal masses of the
fermions can be envisaged as different fermion species belonging to different
electronic bands. The thermal phase diagram of this system was found to be
similar to the one obtained for a mass imbalanced $s$-wave Fermi gas, with two
finite temperature pseudogap regimes even for a population balanced system.
The ground state calculations showed that over a large regime of mass
imbalance the QBP phase is stablized with coexisting gapless $d$-wave
superfluidity and non zero magnetization, which undergoes Lifshitz transition
to an uniform $d$-wave superfluid phase, with decreasing imbalance in mass.

##  0.5 Imbalance and spin-orbit coupling in Fermi systems

Over the past few years there has been immense progress both in the
experimental and theoretical front to understand the spin-orbit coupled
ultracold atomic gases. Since the ultracold atomic gases are one of the
prospective systems to realize exotic superfluid phases such as, FFLO and BP,
there has been obvious interest to investigate the effect of spin-orbit
coupling on imbalanced Fermi gases. In this section we will summarize some of
the major developments in this area and the observations therein. We will also
touch upon some of the open problems in this area towards the end of this
section. For further details on spin orbit coupled Fermi gases the readers are
encouraged to consult the excellent reviews by Galitski and Spielman [168] and
Dalibard et al. [169]. The subject is also discussed in a recent review by
Kinnunen et al. [170].

Spin orbit coupling has remained a subject of interest since the early days of
quantum mechanics and can be understood as follows. If an electron is allowed
to move in an electric field and we move the electron’s rest frame, this will
lead to the generation of a magnetic field by the associated Lorentz
transformation. This magnetic field will couple itself to the electron spin
leading to spin-orbit coupling. The exact form of the spin-orbit coupling
depends upon the direction of the electric field and the electron velocity;
one of the commonly used form of the spin-orbit coupling is the Rashba spin-
orbit coupling (RSOC), given as,

$\displaystyle H_{RSOC}$ $\displaystyle=$
$\displaystyle\alpha(k_{y}\hat{\sigma}_{x}-k_{x}\hat{\sigma}_{y})$ (50)

where, $\alpha$ is the strength of the spin-orbit coupling interaction,
$\hat{\sigma}_{x}$ and $\hat{\sigma}_{y}$ are the Pauli matrices. In the
absence of inversion symmetry the Rashba spin-orbit coupling takes the form of
Dresselhaus spin-orbit coupling and is given as,

$\displaystyle H_{DSOC}$ $\displaystyle=$
$\displaystyle\alpha(k_{y}\hat{\sigma}_{x}+k_{x}\hat{\sigma}_{y})$ (51)

In real space the Rashba spin-orbit coupling is given as,

$\displaystyle H_{RSOC}$ $\displaystyle=$ $\displaystyle-\sum_{\langle
ij\rangle,\sigma}t_{ij,\alpha}\hat{c}_{i,\sigma}^{\dagger}\\{{\bf\sigma}\cdot{\bf\alpha}_{ij}\\}_{\sigma,\sigma^{\prime}}\hat{c}_{j,\sigma^{\prime}}$
(52)

where, ${\bf\sigma}$ is the vector of Pauli matrices and
${\bf\alpha}_{ij}=\alpha\hat{z}\times{\bf r}_{ij}$ is the spin-orbit vector on
each bond of a two-dimensional square lattice. $t_{ij,\alpha}$ correspond to
the strength of the spin-orbit coupling.

In ultracold atomic gas experiments spin-orbit coupling is generated using
synthetic electric field [169] and pseudospins of the atomic species. The
basic principle that is followed in such experiments is to couple different
internal energy states of atoms through a two-photon Raman transition. The
atom first absorbs the photon from one laser beam then emits the photon to a
second laser beam, by stimulated emission. The process conserves momentum by
imparting some of the momentum to the atom. The Doppler shift between the
laser beams and the atom gives rise to the momentum dependent spin-orbit
coupling [171]. Spin-orbit coupling in ultracold atomic gases was
experimentally realised for bosons [172] as well as for fermions [173, 174]. A
detail discussion on the experimental realization of spin-orbit coupling in
ultracold atomic gases can be found in [169, 175, 176, 177].

Theoretically, the effect of spin-orbit coupling on the FFLO superfluid state
has largely been investigated using mean field theory. This has been carried
out either by determining the ground state energy using the FF ansatz in
presence of spin-orbit coupling [178], or by solving the self consistent BdG
equations for arbitrary choice of spatially inhomogeneous superfluid order
parameter, so as to determine the global minima of the energy landscape [179,
180, 181, 182]. Obviously, the second approach is computationally more
demanding than the first one but allows for a larger family of possible ground
state solutions.

Based on BdG mean field theory in two-dimensions, Xu et al. showed that in a
continuum system in presence of spin-orbit coupling the FF state is
energetically favoured over the LO state, particularly at strong spin-orbit
couplings [179]. This observation was reconfirmed by several other works. A
continuum system in three-dimensions with spin-orbit coupling and in presence
of in-plane Zeeman field was studied by Zheng et al. [178]. They showed that
the FF state is stable over a large parameter regime. In a similar spirit, Hu
and Liu pointed out that at sufficiently large spin-orbit coupling the
superfluid state is always an FF state [183]. They worked out the finite
temperature phase diagram at a broad Feshbach resonance and predicted a
superfluid transition temperature of T${}_{c}\sim$ 0.2TF, where, TF is the
Fermi temperature. Similar estimates were made by these authors when they
extended their model by including out of plane Zeeman field [184].

Beyond mean field corrections were incorporated for the continuum model of
three dimensional spin-orbit coupled Fermi gas by Liu [185]. Once again the
superfluid transition temperature was found to be upto $\sim$ 0.2TF. Similar
conclusions were drawn by Dong et al. [186] and Zhou et al. [187] for spin-
orbit coupled three-dimensional Fermi gases.

Within the framework of lattice fermion models, Iskin investigated the spin-
orbit coupled Fermi gas in presence of an optical trap in two dimensions,
using BdG mean field theory [188]. He showed that in presence of spin-orbit
coupling and in-plane Zeeman field the FF state is stabilized, in agreement
with the observations of continuum model. By investigating the finite
temperature behavior Iskin suggested that since the modulations of the
superfluid order parameter are restricted near the edge of the trap, it would
be difficult to detect them at high temperatures [188].

The topological aspect of spin-orbit coupled Fermi gases was investigated by
Tewari et al. [189], who showed that spin-orbit coupling along with Zeeman
field can be used to drive a phase transition between a topologically trivial
and a topologically non-trivial superfluid. Similar observations were reported
for imbalanced Fermi gases in presence of both spin-orbit coupling and Zeeman
field [190, 191]. In their studies carried out on two dimensional continuum
system Zhang et al. [190] and Qu et al. [191] identified a topological FF
state in addition to gapped FF and gapless nodal FF states. They further
confirmed the existence of chiral edge states of the topological FF state, on
a two-dimensional square lattice [191]. Observation of gapless topological FF
states was also reported by Cao et al. [192]. They further confirmed the
robustness of the Majorna modes by adding disorder potential to the system.
Presence of topological FF phase has also been confirmed in one-dimensional
lattices [193].

Mean field theory has been pushed to investigate the finite temperature
behavior of spin-orbit coupled imbalanced Fermi gases, as well [194, 195]. It
was inferred that together with in-plane and out of plane Zeeman field, spin
orbit coupling gives rise to BKT transition in two dimensions, both for gapped
as well as gapless FF phases [194, 195].

In a recent work, Wang et al. investigated the spin-orbit coupled Fermi gas
within the framework of continuum model in two and three dimensions [196].
They took into account the effect of pairing field fluctuations and studied
the system using two schemes, viz. (i) T-matrix approximation and (ii)
fluctuation exchange approximation (FLEX). Based on their observations from
both these approaches Wang et al. inferred that a finite temperature
realization of FFLO phase is not possible in this system. They remarked that
the FFLO phase is unstable against pair fluctuations in continuum in both two
and three dimensions [196]. Within the purview of lattice fermion models, an
auxiliary field QMC study was carried out for balanced two-component Fermi gas
with spin-orbit coupling at the ground state, by Rosenberg et al [197]. It was
shown that such a system comprises of exotic phases such as, supersolid phase
containing both singlet and triplet pairings. Rosenberg et al. showed the
presence of edge currents arising out of the spin-orbit coupling. So far no
such investigation has been carried out in the context of spin-orbit coupled
imbalanced Fermi gases.

As is evident from the above discussions, the theoretical works on spin-orbit
coupled imbalanced Fermi gases are largely the ones based on mean field
theory. This leaves a void in our understanding as far as the effects of
fluctuations are concerned. For example, the exact thermal scales of the spin-
orbit coupled imbalanced Fermi gas in an optical lattice is so far unknown.
The presence of lattice breaks the translational symmetry and can potentially
stabilize phases which are unstable in continuum. It is worth investigating
whether a modulated superfluid phase can be stabilized in a lattice with spin
orbit coupling in presence of pair fluctuations. Moreover, spin-orbit coupling
allows for the possibility of realizing triplet superfluid pairing. Most of
the works we discussed here do not take this possibility into account. Even if
a triplet pairing state does not lead to a long range order at the ground
state, it is certainly possible that triplet pair correlations are present at
high temperatures. In the same spirit, with the energy of the FF and LO phases
being comparable, even though the ground state hosts FF superfluid phase in
presence of spin-orbit coupling, competing LO correlations can be present at
finite temperatures. Signatures of the presence of such competing correlations
would be borne out in the quasiparticle properties of the system.
Investigations of such effects require one to include fluctuations beyond the
mean field theory, and are certainly worth pursuing.

##  0.6 Conclusions and outlook

The last few years have witnessed a great deal of activity in the theoretical
as well as experimental front to unveil the mysteries of the exotic
superconducting phases like FFLO and BP. Though the possibility of
superconducting phases with finite centre of mass momentum of the Cooper pairs
was predicted decades ago, it was the advent of ultracold atomic gases, and
significant advancement of material science experiments, that has paved the
way to access and understand such novel quantum many body phases.

Amongst the solid state materials, so far three candidate materials have stood
out as possible hosts of FFLO phase viz., heavy fermion superconductor
CeCoIn5, organic superconductor $\kappa$-(BEDT-TTF)2Cu(NCS)2 and iron
superconductor KFe2As2. With the continuous development that is going on in
material science, the hope is that other materials would join the class of
Pauli limited superconductors in not so distant future. The ultracold atomic
gas systems with the unprecedented control of the tuning parameters, that it
provides, have taken the centre stage to engineer and analyze exotic quantum
phases, such as FFLO. While a clear spatial signature of the modulated
superfluid phase is yet to be realized in these systems, coexistence of
superfluid pair correlations and unpaired fermions are indeed observed,
suggesting a possible BP phase.

The theoretical attempts are not far behind either, and the problem of
imbalanced Fermi systems have been addressed by mean field as well as non mean
field approaches. We now have a broad understanding of not just the ground
state but also of how the effects of fluctuations play out and dictate the
thermal scales associated with these systems.

The world of quantum many body physics is however a fast evolving one and
newer experimental breakthroughs are not infrequent. For example, experimental
realization of spin-orbit coupling in fermionic and bosonic systems in
ultracold atomic gas setups have been reported recently. In the theoretical
front however, we have not moved beyond the realms of mean field theory so
far, except occasionally. It is obvious that even for the population balanced
scenario, these spin-orbit coupled bosonic and fermionic systems are
theoretically not well understood at this moment and a concerted effort is
required to achieve the same. In the same spirit, the effect of disorder adds
another interesting dimension to these systems. While disorder is largely
expected to be detrimental for the FFLO phase, it can possibly aid in to
stabilize the BP phase, at the ground state of a $s$-wave superfluid.

This chapter has covered some of the important theoretical and experimental
advances and observations pertaining to imbalanced Fermi systems, and is no
way comprehensive. The subject is not yet well understood and promises a
plethora of avenues and areas which needs to be explored inquisitively.

## References

  * [1] Michael Tinkham.  Introduction to $\mathrm{S}$uperconductivity, $\mathrm{M}$c$\mathrm{G}$raw-$\mathrm{H}$ill $\mathrm{B}$ook $\mathrm{C}$o., $\mathrm{N}$ew $\mathrm{Y}$ork, 1975.
  * [2] J. Bardeen, L. N. Cooper, and J. R. Schrieffer.  Theory of $\mathrm{S}$uperconductivity.  Phys. Rev., 108:1175–1204, 1957.
  * [3] J. G. Bednorz and K. A. Müller.  Possible high $\mathrm{T}$c superconductivity in the $\mathrm{B}$a-$\mathrm{L}$a$\mathrm{Cu}$-$\mathrm{0}$ system.  Zeitschrift für Physik B Condensed Matter, 64:189–193, 1986.
  * [4] Patrick A. Lee, Naoto Nagaosa, and Xiao-Gang Wen.  Doping a $\mathrm{M}$ott insulator: $\mathrm{P}$hysics of high-temperature superconductivity.  Rev. Mod. Phys., 78:17–85, 2006.
  * [5] Christian Pfleiderer.  Superconducting phases of $f$-electron compounds.  Rev. Mod. Phys., 81:1551–1624, 2009.
  * [6] G. R. Stewart.  Superconductivity in iron compounds.  Rev. Mod. Phys., 83:1589–1652, 2011.
  * [7] David C. Johnston.  The puzzle of high temperature superconductivity in layered iron pnictides and chalcogenides.  Advances in Physics, 59(6):803–1061, 2010.
  * [8] P. G. de Gennes.  Superconductivity Of $\mathrm{M}$etals $\mathrm{A}$nd $\mathrm{A}$lloys, $\mathrm{T}$aylor and $\mathrm{F}$rancis, $\mathrm{A}$bingdon, 1999.
  * [9] G. Sarma D. Saint James and E. J. Thomas.  Type II $\mathrm{S}$uperconductivity, $\mathrm{P}$ergamon, $\mathrm{N}$ew $\mathrm{Y}$ork, 1969.
  * [10] A. A. Abrikosov.  On the $\mathrm{M}$agnetic $\mathrm{P}$roperties of $\mathrm{S}$uperconductors of the $\mathrm{S}$econd $\mathrm{G}$roup.  Sov. Phys. JETP, 5:1174, 1957.
  * [11] B. S. Chandrasekhar.  A note on the maximum critical field of high$-$field superconductors.  Appl. Phys. Lett., 1:7, 1962.
  * [12] A. M. Clogston.  Upper $\mathrm{L}$imit for the $\mathrm{C}$ritical $\mathrm{F}$ield in $\mathrm{H}$ard $\mathrm{S}$uperconductors.  Phys. Rev. Lett., 9:266–267, 1962.
  * [13] A. Bianchi, R. Movshovich, C. Capan, P. G. Pagliuso, and J. L. Sarrao.  Possible $\mathrm{F}$ulde-$\mathrm{F}$errell-$\mathrm{L}$arkin-$\mathrm{O}$vchinnikov $\mathrm{S}$uperconducting $\mathrm{S}$tate in $\mathrm{CeCo}\mathrm{In}_{5}$.  Phys. Rev. Lett., 91:187004, 2003.
  * [14] R. Lortz, Y. Wang, A. Demuer, P. H. M. Böttger, B. Bergk, G. Zwicknagl, Y. Nakazawa, and J. Wosnitza.  Calorimetric $\mathrm{E}$vidence for a $\mathrm{F}$ulde-$\mathrm{F}$errell-$\mathrm{L}$arkin- $\mathrm{O}$vchinnikov $\mathrm{S}$uperconducting $\mathrm{S}$tate in the $\mathrm{L}$ayered $\mathrm{O}$rganic $\mathrm{S}$uperconductor $\kappa\mathrm{\text{$-$}}(\mathrm{BEDT}\mathrm{\text{$-$}}\mathrm{TTF}{)}_{2}\mathrm{Cu}(\mathrm{NCS}{)}_{2}$.  Phys. Rev. Lett., 99:187002, 2007.
  * [15] D. A. Zocco, K. Grube, F. Eilers, T. Wolf, and H. v. Löhneysen.  Pauli-$\mathrm{L}$imited $\mathrm{M}$ultiband $\mathrm{S}$uperconductivity in $\mathrm{K}\mathrm{Fe}_{2}\mathrm{As}_{2}$.  Phys. Rev. Lett., 111:057007, 2013.
  * [16] Yong-il Shin, Christian H. Schunck, Andre Schirotzek, and Wolfgang Ketterle.  Phase diagram of a two-component $\mathrm{F}$ermi gas with resonant interactions.  Nature, 451:689, 2008.
  * [17] J. B. Ketterson and S. N. Song.  Superconductivity, $\mathrm{C}$ambrigde $\mathrm{U}$niversity $\mathrm{P}$ress, $\mathrm{C}$ambridge, 1999.
  * [18] Peter Fulde and Richard A. Ferrell.  Superconductivity in a $\mathrm{S}$trong $\mathrm{S}$pin-$\mathrm{E}$xchange $\mathrm{F}$ield.  Phys. Rev., 135:A550–A563, 1964.
  * [19] A. I. Larkin and Yu. N. Ovchinnikov.  Nonuniform state of superconductors.  Sov. Phys. JETP, 20:762, 1965.
  * [20] C Mora and R Combescot.  Nature of the $\mathrm{F}$ulde-$\mathrm{F}$errell-$\mathrm{L}$arkin-$\mathrm{O}$vchinnikov phases at low temperature in 2 dimensions.  Europhysics Letters (EPL), 66(6):833–839, 2004.
  * [21] C. Mora and R. Combescot.  Transition to $\mathrm{F}$ulde-$\mathrm{F}$errell-$\mathrm{L}$arkin-$\mathrm{O}$vchinnikov phases in three dimensions: $\mathrm{A}$ quasiclassical investigation at low temperature with $\mathrm{F}$ourier expansion.  Phys. Rev. B, 71:214504, 2005.
  * [22] Hiroshi Shimahara.  Structure of the $\mathrm{F}$ulde-$\mathrm{F}$errell-$\mathrm{L}$arkin-$\mathrm{O}$vchinnikov $\mathrm{S}$tate in $\mathrm{T}$wo-$\mathrm{D}$imensional $\mathrm{S}$uperconductors.  Journal of the Physical Society of Japan, 67(3):736–739, 1998.
  * [23] Jeffrey A. Bowers and Krishna Rajagopal.  Crystallography of color superconductivity.  Phys. Rev. D, 66:065002, 2002.
  * [24] Cornelis J. van der Beek, Marcin Konczykowski, Luc Fruchter, René Brusetti, Thierry Klein, Jacques Marcus, and Christophe Marcenat.  Thermodynamics of the vortex liquid in heavy-ion-irradiated superconductors.  Phys. Rev. B, 72:214504, 2005.
  * [25] Qian Wang, H.-Y. Chen, C.-R. Hu, and C. S. Ting.  Local $\mathrm{T}$unneling $\mathrm{S}$pectroscopy as a $\mathrm{S}$ignature of the $\mathrm{F}$ulde-$\mathrm{F}$errell-$\mathrm{L}$arkin-$\mathrm{O}$vchinnikov $\mathrm{S}$tate in $s$-and $d$-$\mathrm{d}$wave $\mathrm{S}$uperconductors.  Phys. Rev. Lett., 96:117006, 2006.
  * [26] R. Combescot and G. Tonini.  Cascade of two-dimensional $\mathrm{F}$ulde-$\mathrm{F}$errell-$\mathrm{L}$arkin-$\mathrm{O}$vchinnikov phases with anisotropy.  Phys. Rev. B, 72:094513, 2005.
  * [27] Yuji Matsuda and Hiroshi Shimahara.  Fulde, $\mathrm{F}$errell, $\mathrm{L}$arkin, $\mathrm{O}$vchinnikov $\mathrm{S}$tate in $\mathrm{H}$eavy $\mathrm{F}$ermion $\mathrm{S}$uperconductors.  Journal of the Physical Society of Japan, 76(5):051005, 2007.
  * [28] L. G. Aslamazov.  Influence of $\mathrm{I}$mpurities on the $\mathrm{E}$xistence of an $\mathrm{I}$nhomogeneous $\mathrm{S}$tate in a $\mathrm{F}$erromagnetic $\mathrm{S}$uperconductor.  Sov. Phys. JETP, 28:773, 1969.
  * [29] Satoshi Takada.  Superconductivity in a $\mathrm{M}$olecular $\mathrm{F}$ield. ii: $\mathrm{S}$tability of $\mathrm{F}$ulde $-$$\mathrm{F}$errel $\mathrm{P}$hase.  Prog. Theor. Phys., 43:27, 1970.