# Pauli limited d-wave superconductors: Quantum breached pair phase and
thermal transitions

Madhuparna Karmakar  Indian Institute of Technology, Madras, Chennai-600036,
India<EMAIL_ADDRESS>
###### Abstract

We report a quantum phase transition in Pauli limited $d$-wave superconductors
and give the mean field estimates of the associated quantum critical point.
For a population imbalanced $d$-wave superconductor a stable ground state
phase viz. quantum breached pair phase has been identified which comprises of
spatial coexistence of gapless superconductivity and nonzero magnetization.
Based on the thermodynamic and quasiparticle indicators we for the first time
analyze this phase, discuss the thermal behavior of Pauli limited $d$-wave
superconductor, give accurate estimates of the thermal scales associated with
such systems and map out the pseudogap regime. Our work shows that while the
Pauli limited superconductors are known to exhibit exotic modulated
superconducting phase at large imbalance of fermion populations; in the regime
of weak imbalance an intriguing phase of competing orders is realized. We have
established that rather than the superconducting pairing field, it is the
average magnetization of the system that quantifies this quantum phase
transition. Given that the existing Pauli limited superconductors possess
unconventional pairing state symmetry of the superconducting order, our work
promises to open up new avenues in the experimental research of these
materials. We have also demonstrated an alternate scenario wherein the quantum
breached pair phase is a natural outcome for a $d$-wave superconductor with
unequal effective masses of the fermion species.

##  1 Introduction

Superconductivity, in competition or coexistence with magnetic correlations
has been a primary area of research in condensed matter physics, over the past
few decades. While magnetism is usually considered to be detrimental towards
the existence of superconducting order, there are examples of materials where
magnetic correlations reside proximate to or in coexistence with conventional
[1, 2, 3, 4, 5, 6, 7] or unconventional [8, 9, 10, 11, 12, 13, 14, 15, 16, 17]
superconducting orders. The fascinating phenomena of coexisting orders as
observed in these materials is dictated by the precise tuning of the energy
landscape, often carried out via external control parameters such as, doping,
pressure etc [8, 9, 10, 18, 19].

Figure 1: Color online: Magnetization-temperature ($m-T$) phase diagram of
population imbalanced $d$-wave superconductor. The solid curve correspond to
$T_{c}$, the low $m$ regime of second order transition is shown by the black
curve while the red curve at high $m$ correspond to the regime of first order
phase transition.

An equally intriguing but relatively less explored phenomenon of coexisting
superconducting and magnetic correlations is observed in Pauli limited
superconductors. The Pauli limited superconductors are characterized by their
loss of superconducting order via Pauli paramagnetic pair breaking effect,
which originates from the Zeeman splitting of the single electron energy
levels [20, 21]. When a superconducting pair undergoes Pauli paramagnetic pair
breaking it gives rise to unpaired fermions in the system. In the presence of
such unpaired fermions an uniform (zero momentum paired) superconducting state
is no longer stable, rather the superconducting pair acquires a finite
momentum which shows up as spatial modulation in the superconducting order.
The unpaired fermions coexist with such spatially modulated superconducting
state giving rise to what is now well known as the Fulde-Ferrell-Larkin-
Ovchinnikov (FFLO) superconductivity [22, 23]. Pauli limited superconductivity
has been realized both in solid state (e. g. CeCoIn5 [24, 25, 26, 27, 28, 29,
30, 31, 32, 33], $\kappa$-BEDT [34, 35, 36, 37, 38, 39, 40, 41], KFe2As2 [42,
43, 44] etc.) as well as in ultracold atomic gas [45, 46, 47, 48] systems. The
central requirement for realizing Pauli limited superconductivity is the
creation of an imbalance in the population of the fermionic species and thus a
Fermi surface mismatch, via an applied Zeeman field in solid state systems or
by loading different populations of the fermionic species in optical lattice
in ultracold atomic gas set ups.

While much attention has been paid to the FFLO phase over the past few years
[24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42,
43, 44], another exotic phase of the Pauli limited superconductors exhibiting
spatial coexistence of superconductivity and non zero magnetization, remains
almost unexplored. The said phase is known as the “Breached pair” (BP) and
corresponds to a situation where the imbalance in fermionic population is not
strong enough to give rise to a FFLO phase but is sufficient to give rise to a
Fermi surface mismatch [49, 50, 51]. The issue of the BP phase was first
raised decades ago by Sarma [52] 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. However,
he found that such a phase if it exists would be energetically unfavorable as
compared to the uniform superfluid (BCS) phase.

Figure 2: Color online: Ground state energy landscape. (a) Energy optimization
w. r. t. $\phi^{rel}$ at different $\mid U\mid/t$. Inset: $\phi^{rel}$
corresponding to minimum energy state at different $\mid U\mid/t$. (b) Energy
w. r. t population imbalance (Zeeman field) at different $\phi^{rel}$ for
$\mid U\mid=4t$ and $\mu=-0.2t$. The lowest energy state correspond to
$d_{x^{2}-y^{2}}$ pairing state symmetry with $\phi^{rel}=\pi$. (c)
Optimization of pairing field amplitude $\mid\Delta\mid$ at different $h/t$
for $\mid U\mid=4t$, $\mu=-0.2t$ and $\phi^{rel}=\pi$. (d) Optimized pairing
momentum (${\bf q}$) at different Zeeman field ($h/t$) at $\mid U\mid=4t$,
$\mu=-0.2t$ and $\phi^{rel}=\pi$.

Advent of ultracold atomic gas experiments brought renewed interest in the
physics of the BP phase. Through detailed analytic calculations it was
demonstrated by Liu and Wilczek [53, 54, 55] that while for an imbalanced
Fermi system the usual single band Hubbard model with on site interactions
between the fermion species does not allow the BP phase to be a stable ground
state, the same is indeed possible under non trivial circumstances such as,
(i) imbalance in fermion effective masses in presence of contact interaction
between them, (ii) momentum dependent interaction, (iii) same species
repulsion etc. [53, 54, 55]. It was further shown that for suitable choice of
parameters a mass imbalanced $s$-wave superconductor hosts the BP phase as a
stable ground state, which undergoes a second order quantum phase transition
to the uniform superconducting phase as a function of decreasing Zeeman field
[53]. This transition was said to belong to the class of topological Lifshitz
transitions in the sense that the Fermi surface topology changes across this
transition [56]. Note that this topological transition does not involve the
change in any topological invariant. The fate of the BP phase in an imbalanced
unconventional (non $s$-wave) superconductor however continues to remain an
open question till date, in spite of almost all the known Pauli limited
superconductors being unconventional in their pairing state symmetry.

Figure 3: Color online: Ground state phase diagram showing the thermodynamic
phases in the (a) chemical potential-Zeeman field ($\mu-h$) and (b) number
density- magnetization ($n-m$) plane. The order of phase transition from USC
to coexistence, from coexistence to FFLO and from FFLO to PPFL are second,
first and second, respectively, and are marked by the corresponding critical
fields and magnetization. In the $n-m$ plane the USC phase collapses to the
x-axis.

A sufficient volume of literature has already been devoted to theoretically
capture the physics of FFLO phase in the strong imbalance regime of Pauli
limited nodal ($d$-wave) superconductors [57, 58, 59, 60, 61, 62, 63, 64]. In
the present paper we focus on the regime of weak population imbalance in such
superconductors and for the first time demonstrate that the BP phase is a
stable ground state over a significant regime of population imbalance in Pauli
limited $d$-wave superconductors. We refer this ground state BP phase as the
“Quantum Breached Pair” (QBP) phase and establish that this phase comprises of
zero-momentum gapless superconductivity and nonzero magnetization, in
coexistence. A second order quantum phase transition takes place between the
QBP and unmagnetized $d$-wave superconductor (USC) phases, however, there is
no explicit symmetry breaking across this phase transition. We emphasize that
the superconducting pairing field is not a suitable order parameter when it
comes to quantifying the phase transition between QBP and USC phases, as the
pairing field remains unchanged across this transition. The phase transition
between the QBP and USC phases is rather characterized by an alternative order
parameter viz. magnetization ($m_{i}$). As we demonstrate in the following
sections, the ground state phases of the Pauli limited $d$-wave superconductor
are categorized based on the, (i) superconducting pairing field
($\Delta_{ij}$), (ii) single particle density of states at the Fermi level
($N(0)$) and (iii) magnetization ($m_{i}$) as, (a) USC ($\Delta_{ij}({\bf
q}=0)\neq 0$, $N(0)=0$, $m_{i}=0$), (b) QBP ($\Delta_{ij}({\bf q}=0)\neq 0$,
$N(0)\neq 0$, $m_{i}\neq 0$), (c) FFLO ($\Delta_{ij}({\bf q}\neq 0)\neq 0$,
$N(0)\neq 0$, $m_{i}\neq 0$) and (d) partially polarized Fermi liquid (PPFL)
($\Delta_{ij}=0$, $N(0)\neq 0$, $m_{i}\neq 0$), where, ${\bf q}$ corresponds
to the superconducting pairing momentum.

Figure 4: Color online: Zeeman field dependence of spin resolved fermionic
density of states ($N_{\sigma}(\omega)$), at $T=0$. The black (red) curves
correspond to $N_{\uparrow}(\omega)$ ($N_{\downarrow}(\omega)$), respectively.
Note the accumulation of spectral weight at the Fermi level ($\omega=0$) with
increasing field.

Our real space approach to the problem enables us to identify and map out
spatial coexistence of gapless $d$-wave superconducting pairing and nonzero
magnetization, characteristic to the QBP phase. Apart from identifying the
phase transitions at the ground state, the merit of this work rests in
tracking the thermal evolution of these phases using a non perturbative
numerical technique, mapping out the relevant thermal scales associated with
the system and demonstrating that the QBP phase smoothly crosses over to its
classical counterpart with increasing temperature. It must be noted that
within the premises of lattice fermion models the existing works addressing
the thermal physics of the Pauli limited $d$-wave superconductors are limited
to mean field theory (MFT) [63, 64] which are known to significantly
overestimate the thermal scales in the interaction regimes away from the weak
coupling limit [65, 49]. The unconventional $d$-wave superconductors are
established to be in the regime of intermediate coupling, which renders the
MFT unsuitable to address the thermal behavior of these systems. Our technique
takes into account the thermal (spatial) fluctuations of the pairing field at
all orders and not just the saddle point fluctuations. Consequently, it is
well suited to capture the thermal scales associated with the system
accurately and gets progressively accurate with increasing temperature. As
$T\rightarrow 0$, the thermal fluctuations die out and our approach becomes
akin to the MFT at the ground state. In this spirit the quantum phase
transition discussed in this work is basically the mean field estimate of the
same. While we do not expect any qualitative change in our ground state
results via inclusion of quantum fluctuations, quantitative shifts in the
ground state phase boundaries are possible.

While we discuss our results in the following sections, we highlight our
principal observations here: (a) We for the first time demonstrate that the
breached pair phase is a stable ground state of the Pauli limited $d$-wave
superconductor. We refer to this ground state breached pair phase as QBP and
show that it comprises of spatial coexistence of zero-momentum gapless
superconducting order and nonzero magnetization. (b) By tuning the applied
Zeeman field (or imbalance in the populations of the fermionic species) a
second order quantum phase transition can be realized between the USC and QBP
phases, quantified by the evolution of magnetization. In this work, we give
the mean field estimate of this quantum phase transition. (c) Based on the
thermodynamic and quasiparticle signatures we track the thermal evolution of
the ground state phases and map out the relevant thermal scales; we
demonstrate that short range superconducting pair correlations survive up to
temperatures $T\gg T_{c}$ and gives rise to the “pseudogap” regime. (d) We
show that for a $d$-wave superconductor with imbalance in the effective masses
of the fermion species, a QBP phase is realized in the absence of any
imbalance in population. The thermal phase diagram of such a system reveals
the existence of “species dependent” pseudogap scales which should be
experimentally detectable.

The rest of the paper is organized as follows. In section II we discuss the
model and the numerical technique used in this work, as well as the indicators
used to characterize the phases, section III comprises of the main results of
this work and their analysis. We discuss the mass imbalanced $d$-wave
superconductors in section IV, and follow it up by drawing the conclusions of
our work.

##  2 Model, method and indicators

In this section we sketch out the steps involved in studying the many body
quantum Hamiltonian, discuss our choice of the parameters and the relevant
indicators to characterize the phases of this system.

###  2.1 Model

Our starting Hamiltonian corresponds to a two-dimensional (2D) square lattice
with nearest neighbor attractive interaction between the fermions, in presence
of a Zeeman field, and reads as,

$\displaystyle H$ $\displaystyle=$ $\displaystyle-\sum_{\langle
ij\rangle,\sigma}t_{ij}c^{\dagger}_{i,\sigma}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}\sigma_{i}^{z}\hat{n}_{i,\sigma}$
(1)

where, $t_{ij}$ corresponds to the hopping parameter such that $t_{ij}=t=1$
for the nearest neighbors and is zero otherwise, $\mid U\mid>0$ corresponds to
the nearest neighbor attractive interaction between the fermions, the net
number density of the fermions in the system is maintained through the
chemical potential $\mu=(1/2)(\mu_{\uparrow}+\mu_{\downarrow})$, while
$h=(1/2)(\mu_{\uparrow}-\mu_{\downarrow})$ corresponds to the effective Zeeman
field.

###  2.2 Method

The partition function is basically 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)}$ (2)
$\displaystyle{\cal L}(\tau)={\cal L}_{0}(\tau)+{\cal L}_{U}(\tau)$ (3)
$\displaystyle{\cal L}_{0}(\tau)=\sum_{\langle
ij\rangle,\sigma}\\{\bar{\psi}_{i\sigma}((\partial_{\tau}-\mu)\delta_{ij}+t_{ij})\psi_{j\sigma}\\}$
(4) $\displaystyle{\cal L}_{U}(\tau)=-U\sum_{\langle
ij\rangle\sigma\sigma^{\prime}}\bar{\psi}_{i\sigma}\psi_{i\sigma}\bar{\psi}_{j\sigma^{\prime}}\psi_{j\sigma^{\prime}}$
(5)

where, $\beta$ is the inverse temperature. Our strategy is to decompose
$\hat{n}_{i}\hat{n}_{j}$ in terms of the bosonic auxiliary $d$-wave pairing
singlet $\Delta_{ij}(\tau)$ using Hubbard-Stratonovich transformation [66,
67]. $ij$ and $\tau$ refers to the spatial and imaginary time dependence of
the pairing field, respectively. In terms of the Matsubara frequency
$\Omega_{n}=2\pi nT$ the pairing field reads as, $\Delta_{ijn}$, where $T$ is
temperature. 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)}$ (6)
$\displaystyle{\cal L}(\tau)={\cal L}_{0}(\tau)+{\cal L}_{U}(\tau)+{\cal
L}_{cl}(\tau)$ (7) $\displaystyle{\cal L}_{0}(\tau)=\sum_{\langle
ij\rangle,\sigma}\\{\bar{\psi}_{i\sigma}((\partial_{\tau}-\mu)\delta_{ij}+t_{ij})\psi_{j\sigma}\\}$
(8) $\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.$
(9) $\displaystyle{\cal L}_{cl}(\tau)=4\sum_{i\neq
j}\frac{\mid\Delta_{ij}\mid^{2}}{\mid U\mid}$ (10)

Since the fermions are now quadratic, the $\int{\cal D}\psi{\cal
D}{\bar{\psi}}$ integral can be performed to generate the effective action for
the random background fields,

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

where, ${\cal G}$ is the electronic Green’s function in the $\\{\Delta\\}$
background. There are several options now, (i) 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. (ii) 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. (iii) 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$. (iv) Static path approximation (SPA)
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. 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$.
Consequently, it is akin to the MFT only at $T=0$ but captures the thermal
physics of the system accurately.

We choose the last option (SPA) as our numerical technique. The resulting
superconducting Hamiltonian reads as,

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

where, the last term of $H_{SC}$ corresponds to the classical stiffness cost
associated with the auxiliary field. Here, the $d$-wave singlet is defined as
$\Delta_{ij}=(c_{i\uparrow}c_{j\downarrow}+c_{j\uparrow}c_{i\downarrow})$. It
can further be expressed as $\Delta_{ij}=\mid\Delta_{ij}\mid
e^{i{\bf\phi}_{ij}}$, where
${\bf\phi}_{ij}\in\\{\phi_{ij}^{x},\phi_{ij}^{y}\\}$ is the direction
dependent phase of the complex pairing field and $\mid\Delta_{ij}\mid$ is the
pairing field amplitude, considered to be isotropic in the $xy$-plane. Note
that in the usual mean field approach to $H_{SC}$ the superconducting gap
function $\Delta_{ij}$ is assumed to be a real number, but here we retain the
degrees of freedom associated with the pairing field phases and amplitudes.
The system can thus be envisioned as free fermions moving on a random
background of $\Delta_{ij}$.

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

$\displaystyle P\\{\Delta_{ij}\\}\propto Tr_{cc^{\dagger}}e^{-\beta H_{SC}}$
(14)

which is connected to the free energy of the system. For large and random
$\Delta_{ij}$ the trace is taken numerically. We generate the random
background of $\\{\Delta_{ij}\\}$ by using Monte Carlo, diagonalizing $H_{SC}$
for each attempted update of $\Delta_{ij}$. The relevant fermionic correlators
are computed on the optimized configurations at different temperatures.
Evidently, the technique is computationally expensive. The computation cost is
cut down by using a traveling cluster approximation (TCA), wherein instead of
diagonalizing $H_{SC}$ for each attempted update of the auxiliary field, we
diagonalize a smaller cluster surrounding the update site. Both SPA and TCA
have been extensively bench marked and used for several quantum many body
systems. [68, 49, 69, 70].

Figure 5: Color online: Ground state phase diagram at $\mu=-0.2t$ and $\mid
U\mid=4t$ as determined based on the indicators, (i) $x$-component of average
superconducting phase correlation ($\phi^{xx}$), (ii) average magnetization
($m$) and (iii) inverse of DOS at the Fermi level ($1/N(0)$) (see text).

Variational ground state:-

Even for a selected interaction strength the ground state parameter space is
huge in terms of $\mu-h$. In order to get a handle on this parameter space we
first carry out a mean field variational calculation for the ground state. In
the spirit of MFT we stripe the pairing (auxiliary) field $\Delta_{ij}$ off
all fluctuations such that $\mid\Delta\mid$ is now a real number. For a
selected $\mu-h$ cross section we fix the relative phase between the $x$\- and
$y$-components of the superconducting phase $\phi^{rel}=\phi^{y}-\phi^{x}$,
such that $\phi^{rel}$ can take discrete values as, $0,\pi/4,\pi/2...\pi$ etc.
For each such choice of the relative phase, $\mid\Delta\mid$ is optimized so
as to obtain the lowest energy configuration for the selected $\mu-h$ cross
section. The process is carried out for different $\mu-h$ cross sections such
that one obtains the optimized $\mid\Delta\mid$ and $\phi^{rel}$ corresponding
to the minimum energy configuration. The optimized
$\Delta\in\\{\mid\Delta\mid,\phi^{rel}\\}$ configurations are then used to
compute the magnetization of the system for the corresponding $\mu-h$ cross
section. The ground state phase diagram is mapped out based on these
quantities in the $\mu-h$ plane.

Monte Carlo ground state:- The ground state at different $\mu-h$ cross
sections are verified by using Monte Carlo simulations and the results are
found to be in excellent agreement with those obtained via the variational
calculations. Within the framework of Monte Carlo protocol the system is
cooled down from a random high temperature phase and is allowed to attain the
lowest energy state at each temperature. Such unrestricted cooling allows us
to retain the spatial fluctuations of all orders in the pairing field, which
progressively dies out as the system approaches the ground state. For the
Monte Carlo simulations the (a) average superconducting phase correlation
($\phi^{xx}$, $\phi^{yy}$) and (b) the average magnetization ($m$) are used as
suitable indicators to map out the ground state.

Figure 6: Color online: Zeeman field dependence of spin resolved fermionic
occupation number ($n_{\sigma}({\bf k})$), at $T=0$. The $x$ and $y$-axes of
the panels correspond to $k_{x}$ and $k_{y}$, respectively and the intensity
of color shows the magnitude of $n_{\sigma}({\bf k})$. Note that the mismatch
in the Fermi surface is restricted at the $\\{\pm\pi/2,\pm\pi/2\\}$ ${\bf
k}$-points.

###  2.3 Parameters and indicators

In order to study the thermal behavior of the system we select a particular
cross section of the ground state phase diagram, such that the net chemical
potential is selected to be $\mu=-0.2t$ corresponding to a total fermonic
number density of $n\approx 0.94$. The interaction is selected to be $\mid
U\mid=4t$ corresponding to the intermediate interaction regime, which is known
to be suitable for realizing $d$-wave superconductivity [8, 71]. Note that
since the system belongs to the intermediate regime of interaction ($\mid
U\mid\geq t$) a perturbative approach (such as MFT) breaks down at finite
temperatures and the fluctuations dominate. One thus need to resort to non
perturbative techniques such as the one discussed in this article to take into
account the effects of fluctuations. Unlike the assumption of the MFT,
superconductivity in this intermediate regime of interaction is lost at high
temperatures via the loss of long range phase coherence of the superconducting
pairing field rather than the suppression of the pairing field amplitude. The
consequence are the pre-formed pairs at $T\neq 0$ leading to the pseudogap
phase in the $d$-wave superconductors. Our estimates based on the variational
as well as Monte Carlo calculations suggest a $d_{x^{2}-y^{2}}$ pairing state
symmetry with the relative phase between the pairing field components being
$\phi^{rel}=\pi$, in this parameter regime, in agreement with the existing
literature [71]. The results presented in this paper correspond to a system
size of $N=L\times L=24\times 24$, unless specified otherwise. The effect of
the finite system size is discussed towards the end of the paper. We have
analyzed our results based on the thermodynamic and quasiparticle indicators
defined as below,

Site resolved and averaged magnetization:-

$\displaystyle m_{i}=n_{i\uparrow}-n_{i\downarrow}$ $\displaystyle
m=\frac{1}{N}\sum_{i}\langle n_{i\uparrow}-n_{i\downarrow}\rangle$

Average phase correlation of pairing field:-

$\displaystyle\phi^{xx}=\frac{1}{N}\sum_{i\neq j}\langle
e^{i\phi_{i}^{x}}.e^{-i\phi_{j}^{x}}\rangle$
$\displaystyle\phi^{yy}=\frac{1}{N}\sum_{i\neq j}\langle
e^{i\phi_{i}^{y}}.e^{-i\phi_{j}^{y}}\rangle$

Mixed phase correlation of pairing field:-

$\displaystyle\phi^{xy}=\frac{1}{N}\sum_{i\neq j}\langle
e^{i\phi_{i}^{x}}.e^{-i\phi_{j}^{y}}\rangle$

Spin resolved single particle fermionic density of states (DOS):-

$\displaystyle N_{\uparrow}(\omega)=\frac{1}{N}\langle\sum_{i}\mid
u_{n}^{i}\mid^{2}\delta(\omega-E_{n})\rangle$ $\displaystyle
N_{\downarrow}(\omega)=\frac{1}{N}\langle\sum_{i}\mid
v_{n}^{i}\mid^{2}\delta(\omega+E_{n})\rangle$

Spectral lineshapes:-

$\displaystyle A_{\sigma}({\bf k},\omega)=-(1/\pi)ImG_{\sigma}({\bf
k},\omega)$ (15)

Fermion occupation number:-

$\displaystyle n_{\sigma}({\bf k})=\langle c_{{\bf k}\sigma}^{\dagger}c_{{\bf
k}\sigma}\rangle$ (16)

Distribution of pairing field amplitude:-

$\displaystyle P(\mid\Delta_{ij}\mid)=\langle\sum_{i\neq
j}\delta(\mid\Delta\mid-\mid\Delta_{ij}\mid)\rangle$ (17)

where, $i$ and $j$ correspond to two different sites on the lattice.
$\langle\cdots\rangle$ correspond to thermal average and $\sigma$ is the spin
label. $n_{i\sigma}$ are the number of the individual fermionic species, while
$u_{n}^{i}$ and $v_{n}^{i}$ are Bogoliubov eigenfunctions corresponding to the
eigenvalue $E_{n}$. $n_{\sigma}({\bf k})$ is the Fourier transform of the
single particle Green’s function. $G({\bf k},\omega)=lim_{\delta\rightarrow
0}G({\bf k},i\omega_{n})|_{i\omega_{n}\rightarrow\omega+i\delta}$, where
$G({\bf k},i\omega_{n})$ is the imaginary frequency transform of $\langle
c_{\bf k}(\tau)c^{\dagger}_{\bf k}(0)\rangle$.

Figure 7: Color online: Real space maps at $T=0$ showing the evolution of the
QBP phase over a range of Zeeman field. The $x$ and $y$-axes of each panel
correspond to the lattice indices, while the intensity of the color signifies
the magnitude of the corresponding indicator. The first row corresponds to the
amplitude of superconducting pairing field ($\mid\Delta_{ij}\mid$), middle row
shows the $x$-component of phase correlation of the superconducting pairing
field ($\cos(\phi_{x}^{0}-\phi_{x}^{i})$), and the bottom row shows the
magnetization ($m_{i}$). Note that the pairing field undergoes depletion at
isolated regions in space while the emergence of magnetization complements
this depletion.

##  3 Results

Fig. 1 constitutes one of the main results of this work, wherein we show the
magnetization-temperature ($m-T$) phase diagram of Pauli limited $d$-wave
superconductor. The highlight of this phase diagram is the $T=0$ weak
magnetization regime marked as the quantum breached pair (QBP) phase. Finite
temperature leads the QBP phase to undergo smooth crossover to its classical
counterpart and over a large part of the phase diagram the BP phase undergoes
second order transition to partially polarized Fermi liquid (PPFL). The strong
magnetization regime of the phase diagram belongs to the FFLO phase
characterized by finite momentum superconducting pairing and first order
thermal transition to PPFL phase. Sandwiched between the two is the unstable
(phase separated) regime where the BP phase undergoes weak first order thermal
transition to PPFL phase.

The phase diagram shows that the QBP phase undergoes a quantum phase
transition (QPT) to the USC phase, as a function of decreasing magnetization.
The mechanism which drives the QPT is the Pauli paramagnetic pair breaking and
generation of low energy gapless excitation, which hosts unpaired fermions.
For a population imbalanced $s$-wave superconductor the gapless
superconducting state is realized beyond a critical Zeeman field $h\geq
h_{c1}$ (or population imbalance). Concomitant to this phenomena the system
becomes unstable towards ${\bf q}=0$ superconducting pairing and the system
undergoes transition to the FFLO phase. On the contrary, for a nodal $d$-wave
superconductor, the critical field for paramagnetic pair breaking and
generation of gapless superconducting states is different from the one
($h_{c2}$, say) at which the system undergoes transition to FFLO state. Over
the regime $h_{c1}<h\leq h_{c2}$ the system thus hosts a coexistent QBP phase
comprising of “gapless” ${\bf q}=0$ superconducting state and non zero
magnetization. In what follows, we establish and analyze this phase diagram
based on the thermodynamic and quasiparticle indicators.

###  3.1 Ground state

We begin our discussion of the ground state by showing the energy landscape of
the system. Given the huge parameter space we are in the energy optimization
process (via variational calculation) involves multiple steps which are summed
up in Fig. 2 and are discussed below,

  * (a)

As the first step we consider the population balanced system ($h=0.0t$) at a
fixed $\mu=-0.2t$ and optimize the energy at different interactions $U/t$ over
$\mid\Delta\mid$ and $\phi^{rel}$ (in the absence of fluctuations both these
parameters are real numbers). The results obtained via this optimization is
shown in Fig. 2(a) and we note that for weaker interactions the relative phase
corresponding to the minimum energy is $\pi$ showing a $d_{x^{2}-y^{2}}$
pairing state symmetry of the superconducting state, in agreement with the
existing literature [71]. Increasing $U$ shifts the relative phase minima to a
lower value such that for $7t\leq U<9t$ a $d_{x^{2}-y^{2}}+id_{xy}$ pairing
state with $\phi^{rel}=2\pi/3$ is stabilized. At still stronger interactions
$\phi^{rel}$ shifts to even lower values. The results discussed in this
article corresponds to $\mid U\mid=4t$ and our energy minimization suggests a
$d_{x^{2}-y^{2}}$ pairing state as the stable ground state at this
interaction. For the minimized energy the corresponding value of
$\mid\Delta\mid$ (not shown in this figure) gives the pairing field amplitude.
The inset of Fig. 2(a) sums up the change in $\phi^{rel}$ w. r. t. $U/t$.

  * (b)

We next select the particular interaction $\mid U\mid=4t$, introduce
population imbalance via the Zeeman field $h/t$ and optimize the energy over
$|\Delta|(\cos({\bf q.r}))$ and $\phi^{rel}$, where ${\bf q}$ corresponds to
the pairing momentum; for the uniform $d$-wave superconductor ${\bf q}=0$.
Note that there are now three variational parameters as $|\Delta|$,
$\phi^{rel}$ and ${\bf q}$. This is to verify whether the pairing state
undergoes a change (in terms of $\phi^{rel}$) with the imbalance. We show our
results in Fig. 2(b) as the change in energy w. r. t. $h/t$ for selected
$\phi^{rel}$. At all values of $h/t$ the minimum energy configuration clearly
corresponds to the the relative phase $\phi^{rel}=\pi$, at this interaction,
suggesting that the pairing state symmetry of the superconducting state
remains unaltered with population imbalance. This narrows down the parameter
space to $\mid U\mid=4t$, $\mu=-0.2t$ and $\phi^{rel}=\pi$, and different
$h/t$. The remaining task is to determine how the pairing field amplitude
$\mid\Delta\mid$ varies with $h/t$, for this choice of the parameters. This
would enable us to determine the phase boundary between the superconducting
and non superconducting phases at the ground state. Additionally, one needs to
keep track of the pairing momentum ${\bf q}$ so as to be able to determine the
phase boundary between the uniform and non uniform (FFLO) superconducting
phases.

  * (c)

The minimization over $\mid\Delta\mid$ and ${\bf q}$ is carried out over
different choices of spatial modulations such as,
$\Delta_{ij}\sim\mid\Delta\mid\cos(qx_{i})$,
$\Delta_{ij}\sim\mid\Delta\mid\cos[q(x_{i}+y_{j})]$,
$\Delta_{ij}\sim\mid\Delta\mid[\cos(q.x_{i})+\cos(q.y_{j})]$ etc. Since we are
not discussing the FFLO phase in detail in this article, we show the energy
landscape for one such choice of $\Delta_{ij}$ at different $h/t$, in Fig.
2(c). Note that for weak imbalance the $\mid\Delta\mid$ corresponding to the
minimum energy state does not undergo significant change with $h/t$, larger
$h/t$ leads to the development of a weak minimum at smaller $\mid\Delta\mid$
value indicating that population imbalance leads to suppression in the pairing
field amplitude. This is expected because with increasing imbalance there are
now lesser number of fermions which can undergo pairing.

  * (d)

Finally in Fig. 2(d) we show how the pairing momentum ${\bf
q}=\sqrt{(}q_{x}^{2}+q_{y}^{2})$ varies with the population imbalance. For
$h\gtrsim 0.9t$ we note that the pairing momentum ${\bf q}$ picks up a non
zero value indicating an FFLO state. Fig. 2(c) and 2(d) together gives the
phase boundary between the uniform and FFLO superconducting states.

The energy minimization process is carried out over different $\mu-h$ cross
sections so as to map out the ground state phase diagram. By computing the
magnetization over the optimized configurations one can demarcate the regimes
of coexisting superconducting order and non zero magnetization. The ground
state obtained from variational calculations is reconfirmed by Monte Carlo
simulation, wherein fluctuations are taken into account in the pairing field
amplitude and phase, i. e. $\mid\Delta\mid\rightarrow\mid\Delta_{ij}\mid$ and
$\phi^{rel}=\phi_{ij}^{y}-\phi_{ij}^{x}$, respectively.

We next proceed to analyze the ground state thus obtained, in terms of
different indicators. Fig. 3(a) shows the ground state phase diagram of the
system in the $\mu-h$ plane, as determined through Monte Carlo simulation at a
selected interaction strength of $\mid U\mid=4t$. In the $\mu-h$ plane the
system hosts four different phases as, (i) unmagnetized $d$-wave
superconductor (USC), (ii) coexistence, (iii) Fulde-Ferrell-Larkin-Ovchinnikov
(FFLO) and (iv) partially polarized Fermi liquid (PPFL). The phases are
demarcated by critical magnetic fields $h_{c1}$, $h_{c2}$ and $h_{c3}$,
corresponding to a second order phase transition from the USC to coexistence
phase at $h_{c1}\sim 0.2t$, a first order transition of coexistence to FFLO at
$h_{c2}\sim 0.9t$ and a second order transition from FFLO to PPFL phase at
$h_{c3}\sim 1.25t$, respectively. We note that the order of phase transition
from the FFLO to PPFL phase is in agreement with the existing results [63].
The phase transition between the USC and coexistence phase, discussed for the
first time in this paper, is of second order, as expected [53].

The indicators used to demarcate these phases are the $x$-component of the
average phase correlation $\phi^{xx}$ (the $y$-component of the phase
correlation behaves identically) and the average magnetization $m$. These
indicators however does not give any information regarding the spectral
behavior of the underlying state. Consequently, based solely on these
indicators a distinction between the superconducting state being gapped or
gapless can not be made. The regime $h_{c1}<h\leq h_{c2}$ is thus shown as the
“coexistence” phase. While both the USC ($0<h\leq h_{c1}$) and coexistence
($h_{c1}<h\leq h_{c2}$) phases comprise of ${\bf q}=0$ $d$-wave pairing, the
distinction between the gapped and gapless superconductivity can be made via
the single particle DOS, which we discuss next. A (nodal) gapped
superconductor with finite magnetization would correspond to a phase separated
state, while a gapless superconductor with finite magnetization corresponds to
a QBP phase. While our calculations are carried out in the grand canonical
ensemble we recast the phase diagram in the $n-m$-plane in Fig. 3(b). The USC
phase now collapses to the x-axis, while the coexistence, FFLO and PPFL phases
are demarcated by $m_{c1}$ and $m_{c2}$ corresponding to the critical
magnetization. We choose a particular cross section of this phase diagram at
$\mu=-0.2t$ ($n\approx 0.94$) to understand the physics of the system both at
the ground state as well as at the finite temperatures. We have verified that
the qualitative behavior of the system is independent of the choice of $\mu$
(or $n$), except for very small filling. The Zeeman field regime being probed
is $h=[0:1.4]$ where $h=[0:0.9]$ correspond to the “notional” weak imbalance
regime, where the superconducting state has ${\bf q}=0$ pairing.

In order to characterize the underlying superconducting state we next show the
spin resolved single particle DOS at different Zeeman field, in Fig. 4. In the
absence of any population imbalance the DOS corresponding to the two spin
species are identical and centered around the Fermi level ($\omega=0$), giving
rise to a nodal gap. The $d_{x^{2}-y^{2}}$ character of the pairing field is
evident with $N(\omega)\propto\omega$ as $\omega\rightarrow 0$. Increasing
imbalance shifts the spin resolved DOS away from $\omega=0$ with them being
now centered around the shifted Fermi level $\omega=\pm h$. However, an
interesting behavior is observed at $\omega=0$ where there is now an
accumulation of finite spectral weight leading to a gapless superconducting
state. The behavior of the DOS is in dramatic contrast with that of imbalanced
$s$-wave superconductors where even though the DOS corresponding to the
different spin species are centered around the shifted Fermi levels, the
spectrum remains gapped and there are no low energy states at $\omega=0$,
unless the system hits the FFLO phase [49]. The coexistence phase shown in
Fig. 3 thus hosts a gapless ${\bf q}=0$ superconducting phase along with a
nonzero magnetization, characteristic to the QBP phase. The transition from
the gapped to gapless superconducting phases is the signature of the quantum
phase transition between the USC and QBP phases.

The information obtained from the thermodynamic and quasiparticle signatures
is summed up in Fig. 5 which shows the ground state phase diagram at
$\mu=-0.2t$ and $\mid U\mid=4t$. The phase diagram shows only the USC and QBP
phases and for the time being we ignore the ${\bf q}\neq 0$ pairing state.
There are three indicators based on which the phases are demarcated, along
with the $x$-component of the average phase correlation ($\phi^{xx}$) and the
average magnetization ($m$) we now show the inverse of the DOS at the Fermi
level ($1/N(0)$) which gives a measure of the energy gap ($E_{g}\sim 1/N(0)$).
As $N(0)\rightarrow\infty$, $E_{g}\rightarrow 0$, i. e. a large spectral
weight at the Fermi level correspond to a gapless superconducting state. Based
on this indicator we show that for $h\leq 0.2t$ the system in a gapped
superconductor with zero magnetization corresponding to the USC phase. In the
regime $0.2t<h\leq 0.85t$, $1/N(0)$ vanishes indicating a gapless phase,
simultaneous non zero $\phi^{xx}$ and $m$ qualifies the phase to be a QBP. For
$h\gtrsim 0.9t$, $\phi^{xx}$ collapses to zero while $m$ remains finite,
suggesting a non superconducting PPFL phase.

Having mapped out the ground state phase diagram we next move on to a deeper
analysis of these phases. We begin with the spin resolved momentum occupation
number ($n_{\sigma}({\bf k})$) as a function of increasing Zeeman field, shown
in Fig. 6. The momentum occupation number, which is one of the standard
indicators to map out the underlying Fermi surface of the superconducting
state, is the Fourier transform of the single particle Green’s function [72].
In the regime of balanced ($h=0t$) or very weakly imbalanced population the
momentum occupation number corresponding to the two spin species are
identical, i. e. they exhibit the usual Fermi type distribution corresponding
to the $d$-wave superconducting state. Increasing imbalance leads to mismatch
in this distribution and accumulation (or depletion) of weight at isolated
${\bf k}$-points, along the nodal directions ($\pm\pi/2,\pm\pi/2$), of the
Brillouin zone. The $n_{\sigma}({\bf k})$ is thus an important quantity which
gives evidence of the change in the “notional” Fermi surface topology with
population imbalance. In that spirit, $n_{\sigma}({\bf k})$ indicates a
Lifshitz like transition between the USC to QBP phases. Along the nodal
directions the Fermi surfaces changes from being identical to disparate with
increasing imbalance across this transition, which can be considered as a
signature of the Lifshitz like transition [53, 56].

Fig. 7 shows the real space signatures of the QBP phase as the spatial maps
for a typical Monte Carlo configuration of superconducting pairing field and
magnetization at selected Zeeman field cross sections. The rows in the figure
correspond to three different quantities viz. (i) amplitude of superconducting
pairing field ($\mid\Delta_{ij}\mid$), (ii) correlation of the $x$-component
of superconducting phase between a reference site and other sites of the
lattice ($\cos(\phi_{0}^{x}-\phi_{i}^{x})$) and (iii) site resolved
magnetization ($m_{i}$). At the lowest field shown in the figure the
superconducting state is robust with nearly uniform pairing field amplitude
and (quasi) long range phase coherence. On the other hand, magnetization is
very weak ($\sim 0.1$) in this regime. With increasing Zeeman field isolated
islands of depleted superconductivity emerges in the system, as characterized
by regions with suppressed pairing field amplitude and phase correlation. The
islands grow in size and begin to merge with each other as the field is
increased further. Increase in spatial inhomogeneity promotes spectral weight
accumulation at the Fermi level i. e. a gapless superconducting state.
Interestingly, at these regions of depleted superconductivity the
magnetization begins to gain weight as shown in the bottom row of Fig. 7. The
snapshots reveal that over a regime of Zeeman field gapless superconductivity
coexists with nonzero magnetization, characteristic to the QBP phase.

After establishing the ground state behavior of the system we now focus on the
thermal evolution of the various phases. At the onset, we note that the
thermal transition discussed in this paper are basically Berezinsky-
Kosterlitz-Thouless (BKT) transitions, corresponding to algebraic decay of
quasi long range order in two-dimensions.

###  3.2 Finite temperature

We show thermal evolution of the global indicators viz. average
superconducting phase correlations and average magnetization in Fig. 8. In
panels (a)-(c) we show the field dependence of the average phase correlation
components $\phi^{xx}$, $\phi^{yy}$ and $\phi^{xy}$ corresponding to the
$d$-wave superconducting pairing field, across the USC and QBP regime, as the
system evolves in temperature. Note that the superconducting state is fairly
robust in the field regime under consideration, with a nonzero (indicated by
the point of inflection of the curves) $T_{c}$. The system loses the (quasi)
long range phase coherence at $T_{c}$, which sets the thermal transition scale
of the system. In Fig. 8(d) we present the thermal evolution of average
magnetization at different Zeeman field. Increasing temperature leads to
thermal pair breaking and gives rise to unpaired fermions, consequently,
magnetization is larger at high temperatures. For $h\geq 0.2t$ the system is
in QBP phase as indicated by non zero magnetization at $T=0$.

Figure 8: Color online: Thermal evolution of average phase correlation (a)
$\phi^{xx}$, (b) $\phi^{yy}$, (c) $\phi^{xy}$ and (d) average magnetization
($m$), with Zeeman field. The point of inflection of the curves correspond to
the $T_{c}$. Figure 9: Color online: Temperature dependence of the single
particle DOS at the Fermi level corresponding to the the up-spin (a)-(d) and
down-spin (e)-(g) fermion species, at selected Zeeman fields.

The spin resolved DOS at different temperatures, at selected Zeeman fields are
presented next in Fig. 9. For $h>0$ the DOS at the ground state exhibits nodal
gap at the shifted Fermi levels, increasing temperature leads to piling up of
spectral weight and progressive closure of these gaps. Simultaneously, the
coherence peaks at the gap edges flatten out via large transfer of spectral
weight away from the (shifted) Fermi levels, indicating the loss of (quasi)
long range phase coherence. While the depletion of spectral weight at the
shifted Fermi level is observed even at high temperatures owing to the
correlation effects, we note that beyond a temperature, $T_{pg}$, say, the
behavior of the gap closure becomes non monotonic. The temperature range
$T_{c}<T<T_{pg}$ corresponds to the regime where short range superconducting
pair correlations survive in the system even after the loss of (quasi) long
range phase coherence and gives rise to the pseudogap phase. The figure
demonstrates how the pseudogap scale undergoes progressive suppression with
increasing imbalance in population.

Thermal evolution of the momentum resolved spectral lineshapes are presented
next in Fig. 10, for a selected Zeeman field of $h=0.6t$, which corresponds to
the QBP phase at the ground state. The figure shows the lineshape
$A_{\uparrow}({\bf k},\omega)$ for the momentum trajectory
$(0,0)\rightarrow(\pi,0)\rightarrow(\pi,\pi)\rightarrow(0,0)$, across the
Brillouin zone. We note that at the lowest temperature prominent gap opens up
along $(\pi/2,0)\rightarrow(\pi,0)\rightarrow(\pi,\pi/2)$. Progressive
increase in temperature suppresses the gap and finally leads to its closure at
$T\geq 1.25T_{c}$, where $T_{c}\sim 0.08t$. The behavior of
$A_{\downarrow}({\bf k},\omega)$ (not shown here) is similar to
$A_{\uparrow}({\bf k},\omega)$. Overall, in the nodal direction
$(\pm\pi/2,\pm\pi/2)$ the lineshape is characterized by a single peak which is
roughly immune to thermal evolution, while there is opening of gap in the
antinodal direction which undergoes progressive closure with temperature.

The real space signature is demonstrated next in Fig. 11 via the spatial maps.
Once again we use (i) $\mid\Delta_{ij}\mid$, (ii)
$\cos(\phi_{0}^{x}-\phi_{i}^{x})$, and (iii) $m_{i}$, as our indicators and
track them as they evolve in temperature at $h=0.6t$. Apart from the
complementary spatial realization of superconducting pairing field and
magnetization characteristic to the QBP phase at the lowest temperature, we
note enhancement in spatial fragmentation of $\mid\Delta_{ij}\mid$ with
increasing temperature. The observation is validated by the phase coherence
maps which demonstrates the increasing inhomogeneity and thus the loss of
(quasi) long range phase coherence due to thermal fluctuations. The
fragmentation and loss of phase coherence of superconducting pairing field is
accompanied by emergence of regions of large magnetization, a signature of the
“finite temperature BP phase” [49]. This regime of phase uncorrelated
superconducting islands (at $T>T_{c}$) is the visual realization of the
pseudogap phase, and can be thought to be a system of phase uncorrelated
Josephson junctions. As the isolated islands of superconducting pairing field
progressively shrinks with temperature the system loses its superconducting
order and undergoes transition to the PPFL phase.

Figure 10: Color online: Thermal evolution of the momentum resolved spectral
line shape (A${}_{\uparrow}({\bf k},\omega)$) across the
$(0,0)\rightarrow(\pi,0)\rightarrow(\pi,\pi)\rightarrow(0,0)$ trajectory in
the Brillouin zone, at $h=0.6t$. Figure 11: Color online: Real space maps
corresponding to the (a) pairing field amplitude $\mid\Delta_{ij}\mid$, (b)
pairing field phase coherence $\cos(\phi_{0}^{x}-\phi_{i}^{x})$ and (c)
magnetization $m_{i}$, as they evolve in temperature, at a selected Zeeman
field of $h=0.6t$. Increase in temperature leads to loss of phase coherence
and eventual spatial fragmentation of the superconducting state.

Note that both the QBP and the finite temperature pseudogap phases are
characterized by spatially inhomogeneous superconducting state and non zero
spectral weight at the Fermi level of the quasiparticle spectra. The key
distinction between these two phases is the presence (or absence) of (quasi)
long range phase coherence of the superconducting pairing field (i. e. the
superconducting phase stiffness). While the QBP phase is characterized by a
finite superconducting phase stiffness, i. e. a quasi long range phase
coherence between the fragmented superconducting islands, the said phase
coherence is lost in the finite temperature pseudogap phase, indicating the
loss of global superconducting order. Both the USC and QBP phases undergo
smooth cross over to the finite temperature BP phase and then undergoes
transition to the pseudogap phase with increasing temperature.

We now sum up the information gathered based on our analysis of thermodynamic
and quasiparticle indicators and revise the thermal phase diagram shown in
Fig. 1. The revised phase diagram in the $h-T$ plane is presented in Fig. 12.
Along with the thermodynamic phases discussed earlier, the figure shows the
pseudogap phase over a regime of temperature $T_{c}<T<T_{pg}$, and $T_{pg}$
undergoes progressive suppression with Zeeman field. The high field regime
correspond to the FFLO phase and as shown in the phase diagram, is
characterized by large suppression in $T_{c}$, which makes its experimental
realization, non trivial. Over a large part of the $h-T$ plane the system
loses its global superconducting order via a second order thermal phase
transition. The high field low temperature regime hosting the FFLO phase
undergoes a first order thermal phase transition to lose its superconducting
order. The order of these phase transitions are in agreement with the
experimental observations of Pauli limited superconductors such as CeCoIn5
[31].

Next, we take two temperature cross sections of this phase diagram to
highlight the quasiparticle behavior as the system transits from BP to
pseudogap to PPFL regimes, in Fig. 13. In panel (a) we show the up-spin DOS at
$T=0.05t$ as a function of increasing Zeeman field. For $h\lesssim 0.6t$ the
prominent coherence peaks signify (quasi) long range phase coherence while the
system is in the BP phase. At $h\sim 0.7t$ the system is at the verge of
transition to the pseudogap phase as suggested by a large accumulation of
spectral weight at the shifted Fermi level. The coherence peaks smear out
considerably and transfers large spectral weight away from the shifted Fermi
level. The DOS at $h=0.8t$ and $0.9t$ are representative of the pseudogap
regime. Note that the behavior of the DOS in this regime is different from
that observed at $h=0.7t$, owing to the finite temperature short range FFLO
fluctuations at higher Zeeman fields. $h=t$ corresponds to the PPFL phase at
$T=0.05t$ and we find the DOS to be akin to the free electron tight binding
spectra of a square lattice. All superconducting correlations have died out at
this field.

Panel (b) shows the distribution of the superconducting pairing field
amplitude across the $T-h$ plane mentioned in panel (a). We note that for weak
Zeeman fields the distribution remains unchanged from the one observed for a
$d$-wave superconductor, at $h=0t$. The mean amplitude of the pairing field
also remains roughly constant. Increasing field shifts most of the weight to
low amplitude $\mid\Delta_{ij}\mid$s, indicating suppression of the pairing by
Zeeman field. The distribution shows that as expected the FFLO phase comprises
of superconducting pairing with suppressed amplitudes.

A high temperature scan at $T=0.15t>T_{c0}$ is shown in panel (c), where
$T_{c0}$ correspond to the $T_{c}$ at $h=0t$. At all values of $h$ we find
that the magnitude of the coherence peak is reduced to almost half of what was
observed at $T=0.05t$, as there are now only short range pair correlations
surviving in the system without any (quasi) long range order. Consequently,
the coherence peaks have smeared out as is expected from the pseudogap phase.
The system transits to the PPFL phase at $h\sim 0.7t$, indicated by the
complete disappearance of the coherence peaks. The weak depletion in the DOS
at the shifted Fermi level as observed for $h\gtrsim 0.7t$ arises due to the
strong correlation regime we are in, and continues to survive even at higher
temperatures. The corresponding distribution of $\mid\Delta_{ij}\mid$ is shown
in panel (d). We notice that the distribution is significantly broader due to
thermal fluctuations and is roughly independent of the choice of $h$ since the
long range superconducting pair correlations have died out at this
temperature, as compared to its $T=0.05t$ counterpart.

Figure 12: Color online: Zeeman field-temperature ($h-T$) phase diagram of
population imbalanced $d$-wave superconductor, showing the BP, QBP, FFLO, PPFL
and pseudogap phases. The solid curve correspond to $T_{c}$, with the regime
of second order transition shown by black curve while the red curve shows the
regime of first order phase transition. The dashed curve correspond to
$T_{pg}$ marking the crossover from the pseudogap to the PPFL phase.

##  4 Discussion and conclusion

###  4.1 Imbalance in mass

In the earlier sections of this paper we have discussed about how a Pauli
limited $d$-wave superconductor gives rise to a quantum breached pair phase
with coexisting zero-momentum gapless superconducting order and nonzero
magnetization. We demonstrated that due to the shift in the Fermi level of the
individual fermion species there is a pile up of spectral weight at the
“unshifted” Fermi level as suggested by the spatial depletion of the
superconducting order. The regions of depleted superconductivity serves as
host to the unpaired fermions in the system which gives rise to nonzero
magnetization.

Figure 13: Color online: Single particle DOS ($N_{\uparrow}(\omega)$) ((a) and
(c)) and distribution of pairing field amplitude $P(\mid\Delta_{ij}\mid)$ ((b)
and (d)) at selected temperature cross sections $T=0.05t$ and $T=0.15t$,
respectively, for different Zeeman fields. $\mid\Delta\mid$ correspond to the
pairing field amplitude at the ground state.

In this section, we discuss an alternative scenario where similar physical
phenomena plays out. At this end we consider a Fermi-Fermi mixture with
fermion species having unequal effective masses. In the context of ultracold
atomic gases such Fermi-Fermi mixtures have already been realized as,
$Li^{6}$-$K^{40}$ mixture [73, 74, 75, 76], albeit with an isotropic $s$-wave
symmetry of the pairing state. While superfluidity is yet to be achieved in
such mixtures, the degenerate Fermi regime [73, 76], Fermi resonance between
$Li^{6}$-$K^{40}$ atoms [77, 75, 76] and formation of $Li^{6}$-$K^{40}$ hetero
molecules [74] is already a reality. Moreover, other Fermi-Fermi mixtures such
as, $Dy^{161}$, $Dy^{163}$, $Er^{167}$ etc. are expected to be experimentally
realizable in the near future [78, 79]. The theoretical efforts that have been
put in to understand such hetero molecules are primarily concentrated on
continuum models. Density functional theory combined with local density
approximation [80], functional renormalization group analysis [81], mean field
theory with Gaussian fluctuations [82, 83], T-matrix and extended T-matrix
approaches [84, 85, 86] have been employed to understand the physics of Fermi-
Fermi mixtures. Within the purview of lattice fermion models, quantum Monte
Carlo study on one dimensional mass imbalanced system [87], non perturbative
lattice Monte Carlo based analysis of two-dimensional system [88] and a recent
static path approximation based Monte Carlo study [89] are some of the efforts
worth mentioning.

In the present section we present a similar scenario with the pairing symmetry
being non local, a target that can be achieved in experiments pertaining to
ultracold atomic gases. 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. Below we demonstrate that a quantum
breached pair state and the associated phase transitions can be realized in
such systems even without an imbalance in the fermionic populations.

The Fermi-Fermi mixture comprising of unequal mass fermion species subjected
to a non local interaction can be depicted by the Hamiltonian,

$\displaystyle H_{SC}$ $\displaystyle=$ $\displaystyle-\sum_{\langle
ij\rangle,\alpha}t_{ij}^{\alpha}(c_{i\alpha}^{\dagger}c_{j\alpha}+h.c.)+\sum_{i\neq
j}\Delta_{ij}(c_{iL}^{\dagger}c_{jH}^{\dagger}+c_{jL}^{\dagger}c_{iH}^{\dagger})+h.c.$
(18) $\displaystyle-\mu\sum_{i\alpha}\hat{n}_{i\alpha}+4\sum_{i\neq
j}\frac{\mid\Delta_{ij}\mid^{2}}{\mid U\mid}$

where, $\alpha=L,H$ corresponds to the light and heavy fermion species,
respectively. Note that the imbalance in mass is imbibed in the hopping
parameter $t_{ij}^{\alpha}$, wherein, $t_{ij}^{\alpha}\sim 1/m_{\alpha}$;
$m_{H}$ and $m_{L}$ correspond to the masses of the heavy and light species,
respectively. Once again we set the hopping to be nearest neighbor and the
energy scales are measured in units of $t^{L}$, which is set to unity. As our
tuning parameter we define the ratio between the hopping parameters
corresponding to the fermion species as, $\eta=t_{H}/t_{L}=m_{L}/m_{H}$. Thus,
$\eta=1$ corresponds to the balanced (equal mass) limit of the system and
$\eta=0$ corresponds to maximum imbalance, wherein one of the fermion species
is completely localized.

Figure 14: Color online: (a) $x$-component of average pairing field phase
coherence $\phi^{xx}$ and (b) average magnetization ($m(T)$) for $d$-wave
superconductor with imbalance in the effective masses of the fermion species.
At large imbalance $\eta\lesssim 0.4$ the $T=0$ state shows coexisting
superconducting order and non zero magnetization.

We show the thermodynamic behavior corresponding to this system in Fig. 14, in
terms of the $x$-component of the average phase correlation $\phi^{xx}$ and
average magnetization $m$, at different $\eta$s. We note that though the phase
correlation gets progressively suppressed with increasing mass imbalance, over
a regime of imbalance $\eta\leq 0.4$ the system hosts nonzero magnetization
along with a finite phase correlation, at the ground state. In order to verify
whether this coexistence phase is a QBP phase we have examined the single
particle DOS and have found them to be gapless.

In the limit of mass balance ($\eta=1$) the DOS corresponding to the two
species exhibit nodal gap at the Fermi level corresponding to an uniform
$d$-wave superconductor. As the imbalance increases one of the species gets
semi localized, as a result the gain in the condensation energy through
pairing is less than that in case of the balanced mixture, i. e. not all the
fermions get paired up. There is now accumulation of finite spectral weight at
the Fermi level, the resulting gapless spectra accommodates the excess
(unpaired) fermions and give rise to a coexistence phase characteristic to the
QBP. For weaker interaction strength, gapless superconductivity is realized
for smaller imbalance ($\eta\rightarrow 1$ as $\mid U\mid\rightarrow 0$).

In Fig. 15 we show the composite phase diagrams for the proposed scenario of
Fermi-Fermi mixture, both at the ground state as well as at finite
temperature. The ground state phase diagram is shown in panel (a). In the
regime of $\eta\sim 0$ the system is a partially polarized Fermi liquid (PPFL)
with $m\neq 0$ and $\phi^{xx}=0$, with increasing $\eta$ the pairing field
phase correlation $\phi^{xx}$ picks up weight and progressively increases as
more fermions participate in the pairing. Simultaneously, the magnetization
gets suppressed $m\rightarrow 0$ (due to decrease in the number of unpaired
fermions) and eventually drops to zero at $\eta\sim 0.3$, beyond which the
system is an unmagnetized $d$-wave superconductor. The regime with $m\neq 0$
and $\phi^{xx}\neq 0$ correspond to the QBP phase, in the phase diagram.
$\eta\sim 0.3$ is thus the mean field estimate of the quantum critical point
for the quantum phase transition between the USC and QBP phases, at the
parameters under consideration. The average phase correlation ($\phi^{xx}$)
continues to increase and becomes saturated for $\eta>0.5$, to the magnitude
expected from the balanced system.

In panel (b) we present the corresponding finite temperature behavior. The
figure shows four distinct phases as, (i) superconductor (BP), (ii) PG-I,
(iii) PG-II and (iii) PPFL; along with three important thermal scales as, (i)
$T_{c}$, (ii) $T_{pg}^{I}$ and (iii) $T_{pg}^{II}$. Before we characterize the
phases individually we emphasize that unlike a system with imbalance in
fermion populations, there are two pseudogap regimes for systems with mass
imbalance, as shown in the figure. In one of the earlier works the author had
carried out detail analysis of mass imbalanced Fermi-Fermi mixtures with on-
site interaction [89], such as, $Li^{6}$-$K^{40}$. The analysis showed that
since the two fermion species are being subjected to different “scaled”
temperatures, the regime over which the short range correlations survive in
each of them are different. Consequently, while in the PG-I regime both the
species are pseudogapped, in the PG-II regime it is only the lighter species
which is pseudogapped while the heavier species is a partially polarized Fermi
liquid. We do not expect any qualitative change from this picture when the
pairing state symmetry is $d$-wave, as demonstrated in Fig. 15(b). As
expected, the BP regime involves finite temperature coexistence of $d$-wave
superconductivity and non zero magnetization, while the PPFL phase is highly
magnetized with vanishing superconducting correlations. The transitions from
superconductor to PG-I, crossover from PG-I to PG-II and from PG-II to PPFL
are marked by $T_{c}$, $T_{pg}^{I}$ and $T_{pg}^{II}$, respectively. As
$\eta=1.0$ corresponds to balanced limit, the $T_{pg}^{I}$ and $T_{pg}^{II}$
scales collapse to one at this point. The survival of the pseudogap regimes
upto $T\gg T_{c}$ ensures that signatures of mass imbalance in such systems
can be accessed through species resolved spectroscopic experiments, even when
the superconducting transition scales are strongly suppressed due to
imbalance.

Figure 15: Color online: (a) Ground state and (b) thermal phase diagram of
$d$-wave superconductor with unequal effective masses of the fermion species.
Thermal phase diagram shows species resolved pseudogap phases, with both the
species being pseudogaped in the PG-I regime and only the lighter species
being pseudogapped in the PG-II regime.

In conclusion, in this paper we have carried out a systematic study of Pauli
limited $d$-wave superconductors, within the purview of a lattice fermion
model. Unlike the existing body of literature, we focus our analysis on the
regime of population imbalance where it is not strong enough to give rise to
the exotic FFLO phase, but at the same time is sufficient to give rise to
Fermi surface mismatch and unpaired fermions. To the best of our knowledge, we
for the first time demonstrate that in $d$-wave Pauli limited superconductors,
imbalance in fermion populations give rise to a quantum phase transition from
an uniform superconductor to a quantum breached pair phase. We give a mean
field estimate of the quantum critical point of this transition. This quantum
phase transition is quantified by the average magnetization of the system,
rather than the superconducting pairing field. The merit of this work rests
not just in establishing the quantum breached pair phase of the Pauli limited
$d$-wave superconductors but also in accessing the thermal phases of such
systems. While the existing literature on the Pauli limited $d$-wave
superconductors seem to be restricted to studies based on the mean field
theory, we for the first time implement a non perturbative numerical technique
which enables us to access the thermal scales accurately, owing to its
inclusion of the spatial fluctuations of the superconducting pairing field.
Apart from the population imbalanced $d$-wave superconductors we have
discussed about an alternate scenario where a quantum breached pair phase can
be realized, viz. a mass imbalanced $d$-wave superconductor. We have mapped
out the ground state and finite temperature phase diagrams for such systems
and showed that while the ground state does host a quantum breached pair
phase, the finite temperature phase comprises of species selective regimes for
the survival of short range pair correlations, giving rise to two pseudogap
phases.

While much attention has been paid to understand the FFLO physics of the Pauli
limited $d$-wave superconductors, this is the first work which establishes the
existence of a quantum breached pair state in these systems over a significant
regime of imbalance. We have discussed several thermodynamic and quasiparticle
indicators which should be accessible to the existing experimental probes. We
believe that this work is likely to open up exciting new avenues for
experimental research to observe quantum breached pair phase in solid state
materials as well as in ultracold atomic gas setups.

##  5 Acknowledgements

The author duly acknowledges the use of HPC cluster at Harish Chandra Research
Institute, Prayagraj (Allahabad), India and Virgo HPC cluster at IIT, Madras,
India.

##  Appendix A Appendix

Figure A1: Color online: Thermal evolution of the $x$-component of average
superconducting phase correlation (a) $\phi^{xx}$ and (b) average
magnetization ($m$), with Zeeman field, for system size of $L=36$.

###  A.1 Finite size effect

The results presented in the main text correspond to a particular lattice size
of $L=24$. Any lattice calculation is however likely to be plagued by finite
size effect. In order to verify whether the quantum breached pair phase
discussed in this paper is an artifact of finite lattice sizes we have carried
out the simulations at different system sizes (upto $L=40$). In Fig. A1 we
show the thermal evolution of average of $x$-component of the superconducting
phase correlation and average magnetization at $L=36$ for different Zeeman
field. As is evident from the figure, there is indeed nonzero magnetization at
$T=0$ for $h\geq h_{c1}\sim 0.2t$ suggesting that the quantum breached pair
phase is robust and stable against finite lattice size effects.

## References

## References

  * [1] Fulde P and Keller J Superconductivity and Magnetism, edited by M. B. Maple and O. Fischer (Springer, Berlin), 1982.
  * [2] Muller K H and Narozhnyi V N 2001 Rep. Prog. Phys. 64 943 URL https://doi.org/10.1088
  * [3] Lynn J W, Skanthakumar S, Huang Q, Sinha S K, Hossain Z, Gupta L C, Nagarajan R and C G 1997 Phys. Rev. B 55 6584 URL https://link.aps.org/doi/10.1103/PhysRevB.55.6584
  * [4] Braun H F Rare Earth Transition Metal Borocarbides (Nitrides): Superconducting, Magnetic and Normal State Properties, edited by K. H. Muller and V. N. Narozhnyi, (Kluwer, Dordrecht), 2001.
  * [5] Schneider M, Fuchs G, Muller K H, Nenkov K, Behr G, Souptel D and Drechsler S L 2009 Phys. Rev. B 80 224522 URL https://link.aps.org/doi/10.1103/PhysRevB.80.224522
  * [6] Baba T, Yokoya T, Tsuda S, Kiss T, Shimojima T, Ishizaka K, Takeya H, Hirata K, Watanabe T, Nohara M, Takagi H, Nakai N, Machida K, Togashi T, Watanabe S, Wang X Y, Chen C T and Shin S 2008 Phys. Rev. Lett. 100 017003 URL https://link.aps.org/doi/10.1103/PhysRevLett.100.017003
  * [7] Naidyuk Y G, Kvitnitskaya O E, Tiutrina L V, Yanson I K, Behr G, Fuchs G, Drechsler S L, Nenkov K and Schultz L 2011 Phys. Rev. B 84 094516 URL https://link.aps.org/doi/10.1103/PhysRevB.84.094516
  * [8] Lee P A, Nagaosa N and Wen X G 2006 Rev. Mod. Phys. 78(1) 17–85 URL https://link.aps.org/doi/10.1103/RevModPhys.78.17
  * [9] Pfleiderer C 2009 Rev. Mod. Phys. 81(4) 1551–1624 URL https://link.aps.org/doi/10.1103/RevModPhys.81.1551
  * [10] Stewart G R 2011 Rev. Mod. Phys. 83(4) 1589–1652 URL https://link.aps.org/doi/10.1103/RevModPhys.83.1589
  * [11] Nicklas M, Stockert O, Park T, Habicht K, Kiefer K, Pham L D, Thompson J D, Fisk Z and Steglich F 2007 Phys. Rev. B 76(5) 052401 URL https://link.aps.org/doi/10.1103/PhysRevB.76.052401
  * [12] Pham L D, Park T, Maquilon S, Thompson J D and Fisk Z 2006 Phys. Rev. Lett. 97(5) 056404 URL https://link.aps.org/doi/10.1103/PhysRevLett.97.056404
  * [13] Johnston D C 2010 Advances in Physics 59 803–1061 URL https://doi.org/10.1080/00018732.2010.513480
  * [14] Hsu F C, Luo J Y, Yeh K W, Chen T K, Huang T W, Wu P M, Lee Y C, Huang Y L, Chu Y Y, Yan D C and Wu M K 2008 Proceedings of the National Academy of Sciences 105 14262–14264 ISSN 0027-8424 URL https://www.pnas.org/content/105/38/14262
  * [15] Fang M H, Pham H M, Qian B, Liu T J, Vehstedt E K, Liu Y, Spinu L and Mao Z Q 2008 Phys. Rev. B 78(22) 224503 URL https://link.aps.org/doi/10.1103/PhysRevB.78.224503
  * [16] Guo J, Jin S, Wang G, Wang S, Zhu K, Zhou T, He M and Chen X 2010 Phys. Rev. B 82(18) 180520 URL https://link.aps.org/doi/10.1103/PhysRevB.82.180520
  * [17] Fang M H, Wang H D, Dong C H, Li Z J, Feng C M, Chen J and Yuan H Q 2011 EPL (Europhysics Letters) 94 27009 URL https://doi.org/10.1209
  * [18] Dai P, Hu J and Dagotto E 2012 Nature Physics 8 709 URL https://doi.org/10.1038/nphys2438
  * [19] Wiesenmayer E, Luetkens H, Pascua G, Khasanov R, Amato A, Potts H, Banusch B, Klauss H H and Johrendt D 2011 Phys. Rev. Lett. 107(23) 237001 URL https://link.aps.org/doi/10.1103/PhysRevLett.107.237001
  * [20] Matsuda Y and Shimahara H 2007 Journal of the Physical Society of Japan 76 051005 URL https://doi.org/10.1143/JPSJ.76.051005
  * [21] D Saint James G S and Thomas E J Type II Superconductivity, Pergamon, New York, 1969
  * [22] Fulde P and Ferrell R A 1964 Phys. Rev. 135(3A) A550–A563 URL https://link.aps.org/doi/10.1103/PhysRev.135.A550
  * [23] Larkin A I and Ovchinnikov Y N 1965 Sov. Phys. JETP 20 762
  * [24] Bianchi A, Movshovich R, Capan C, Pagliuso P G and Sarrao J L 2003 Phys. Rev. Lett. 91(18) 187004 URL https://link.aps.org/doi/10.1103/PhysRevLett.91.187004
  * [25] Tayama T, Harita A, Sakakibara T, Haga Y, Shishido H, Settai R and Onuki Y 2002 Phys. Rev. B 65(18) 180504 URL https://link.aps.org/doi/10.1103/PhysRevB.65.180504
  * [26] Koutroulakis G, Mitrovic V F, Horvatic M, Berthier C, Lapertot G and Flouquet J 2008 Phys. Rev. Lett. 101(4) 047004 URL https://link.aps.org/doi/10.1103/PhysRevLett.101.047004
  * [27] Kenzelmann M, Strässle T, Niedermayer C, Sigrist M, Padmanabhan B, Zolliker M, Bianchi A D, Movshovich R, Bauer E D, Sarrao J L and Thompson J D 2008 Science 321 1652–1654 URL https://science.sciencemag.org/content/321/5896/1652
  * [28] Capan C, Bianchi A, Movshovich R, Christianson A D, Malinowski A, Hundley M F, Lacerda A, Pagliuso P G and Sarrao J L 2004 Phys. Rev. B 70(13) 134513 URL https://link.aps.org/doi/10.1103/PhysRevB.70.134513
  * [29] Martin C, Agosta C C, Tozer S W, Radovan H A, Palm E C, Murphy T P and Sarrao J L 2005 Phys. Rev. B 71(2) 020503 URL https://link.aps.org/doi/10.1103/PhysRevB.71.020503
  * [30] Gerber S, Bartkowiak M, Gavilano J L, Ressouche E, Egetenmeyer N, Niedermayer C, Bianchi A D, Movshovich R, Bauer E D, Thompson J D and Kenzelmann M 2013 Nature Physics 10 126 URL https://doi.org/10.1038/nphys2833
  * [31] Kumagai K, Saitoh M, Oyaizu T, Furukawa Y, Takashima S, Nohara M, Takagi H and Matsuda Y 2006 Phys. Rev. Lett. 97(22) 227002 URL https://link.aps.org/doi/10.1103/PhysRevLett.97.227002
  * [32] Kim D Y, Lin S Z, Weickert F, Kenzelmann M, Bauer E D, Ronning F, Thompson J D and Movshovich R 2016 Phys. Rev. X 6(4) 041059 URL https://link.aps.org/doi/10.1103/PhysRevX.6.041059
  * [33] Shi-Zeng Lin Duk Y Kim E D B F R J D T and Movshovich R arXiv:1902.04797v1
  * [34] Beyer R and Wosnitza J 2013 Low Temperature Physics 39 225–231 URL https://doi.org/10.1063/1.4794996
  * [35] Lortz R, Wang Y, Demuer A, Böttger P H M, Bergk B, Zwicknagl G, Nakazawa Y and Wosnitza J 2007 Phys. Rev. Lett. 99(18) 187002 URL https://link.aps.org/doi/10.1103/PhysRevLett.99.187002
  * [36] Mayaffre H, Kramer S, Horvatić M, Berthier C, Miyagawa K, Kanoda K and Mitrovic V F 2014 Nature Physics 10 928 URL https://doi.org/10.1038/nphys3121
  * [37] Wright J A, Green E, Kuhns P, Reyes A, Brooks J, Schlueter J, Kato R, Yamamoto H, Kobayashi M and Brown S E 2011 Phys. Rev. Lett. 107(8) 087002 URL https://link.aps.org/doi/10.1103/PhysRevLett.107.087002
  * [38] Coniglio W A, Winter L E, Cho K, Agosta C C, Fravel B and Montgomery L K 2011 Phys. Rev. B 83(22) 224507 URL https://link.aps.org/doi/10.1103/PhysRevB.83.224507
  * [39] Bergk B, Demuer A, Sheikin I, Wang Y, Wosnitza J, Nakazawa Y and Lortz R 2011 Phys. Rev. B 83(6) 064506 URL https://link.aps.org/doi/10.1103/PhysRevB.83.064506
  * [40] Agosta C C, Jin J, Coniglio W A, Smith B E, Cho K, Stroe I, Martin C, Tozer S W, Murphy T P, Palm E C, Schlueter J A and Kurmoo M 2012 Phys. Rev. B 85(21) 214514 URL https://link.aps.org/doi/10.1103/PhysRevB.85.214514
  * [41] Cho K, Smith B E, Coniglio W A, Winter L E, Agosta C C and Schlueter J A 2009 Phys. Rev. B 79(22) 220507 URL https://link.aps.org/doi/10.1103/PhysRevB.79.220507
  * [42] Zocco D A, Grube K, Eilers F, Wolf T and Löhneysen H v 2013 Phys. Rev. Lett. 111(5) 057007 URL https://link.aps.org/doi/10.1103/PhysRevLett.111.057007
  * [43] Cho K, Kim H, Tanatar M A, Song Y J, Kwon Y S, Coniglio W A, Agosta C C, Gurevich A and Prozorov R 2011 Phys. Rev. B 83(6) 060502 URL https://link.aps.org/doi/10.1103/PhysRevB.83.060502
  * [44] Khim S, Lee B, Kim J W, Choi E S, Stewart G R and Kim K H 2011 Phys. Rev. B 84(10) 104502 URL https://link.aps.org/doi/10.1103/PhysRevB.84.104502
  * [45] Shin Y i, Schunck C H, Schirotzek A and Ketterle W 2008 Nature 451 689 URL https://doi.org/10.1038/nature06473
  * [46] Schunck C H, Shin Y, Schirotzek A, Zwierlein M W and Ketterle W 2007 Science 316 867–870 URL https://science.sciencemag.org/content/316/5826/867
  * [47] Shin Y, Zwierlein M W, Schunck C H, Schirotzek A and Ketterle W 2006 Phys. Rev. Lett. 97(3) 030401 URL https://link.aps.org/doi/10.1103/PhysRevLett.97.030401
  * [48] Liao Y a, Rittner A S C, Paprotta T, Li W, Partridge G B, Hulet R G, Baur S K and Mueller E J 2010 Nature 467 567 URL https://doi.org/10.1038/nature09393
  * [49] Karmakar M and Majumdar P 2016 Phys. Rev. A 93(5) 053609 URL https://link.aps.org/doi/10.1103/PhysRevA.93.053609
  * [50] Karmakar M and Majumdar P 2016 The European Physical Journal D 70 220 URL https://doi.org/10.1140/epjd/e2016-70250-2
  * [51] Koponen T, Kinnunen J, Martikainen J P, Jensen L M and Törmä P 2006 New Journal of Physics 8 179–179 URL https://doi.org/10.1088%2F1367-2630%2F8%2F9%2F179
  * [52] Sarma G 1963 Journal of Physics and Chemistry of Solids 24 1029 – 1032 ISSN 0022-3697 URL http://www.sciencedirect.com/science/article/pii/0022369763900076
  * [53] Liu W V, Wilczek F and Zoller P 2004 Phys. Rev. A 70(3) 033603 URL https://link.aps.org/doi/10.1103/PhysRevA.70.033603
  * [54] Forbes M M, Gubankova E, Liu W V and Wilczek F 2005 Phys. Rev. Lett. 94(1) 017001 URL https://link.aps.org/doi/10.1103/PhysRevLett.94.017001
  * [55] Liu W V and Wilczek F 2003 Phys. Rev. Lett. 90(4) 047002 URL https://link.aps.org/doi/10.1103/PhysRevLett.90.047002
  * [56] Lifshitz I M 1960 Sov. Phys. JETP 11 1130
  * [57] Vorontsov A B, Sauls J A and Graf M J 2005 Phys. Rev. B 72(18) 184501 URL https://link.aps.org/doi/10.1103/PhysRevB.72.184501
  * [58] Vorontsov A B and Graf M J 2006 Phys. Rev. B 74(17) 172504 URL https://link.aps.org/doi/10.1103/PhysRevB.74.172504
  * [59] Kaczmarczyk J, Sadzikowski M and Spałek J 2011 Phys. Rev. B 84(9) 094525 URL https://link.aps.org/doi/10.1103/PhysRevB.84.094525
  * [60] Beaird R, Vorontsov A B and Vekhter I 2010 Phys. Rev. B 81(22) 224501 URL https://link.aps.org/doi/10.1103/PhysRevB.81.224501
  * [61] Yanase Y 2008 Journal of the Physical Society of Japan 77 063705 URL https://doi.org/10.1143/JPSJ.77.063705
  * [62] Vorontsov A B, Vekhter I and Graf M J 2008 Phys. Rev. B 78(18) 180505 URL https://link.aps.org/doi/10.1103/PhysRevB.78.180505
  * [63] Zhou T and Ting C S 2009 Phys. Rev. B 80(22) 224515 URL https://link.aps.org/doi/10.1103/PhysRevB.80.224515
  * [64] Cui Q, Hu C R, Wei J Y T and Yang K 2012 Phys. Rev. B 85(1) 014503 URL https://link.aps.org/doi/10.1103/PhysRevB.85.014503
  * [65] Randeria M and Taylor E 2014 Annual Review of Condensed Matter Physics 5 209–232 URL https://doi.org/10.1146/annurev-conmatphys-031113-133829
  * [66] Hubbard J 1959 Phys. Rev. Lett. 3(2) 77–78 URL https://link.aps.org/doi/10.1103/PhysRevLett.3.77
  * [67] Schulz H J 1990 Phys. Rev. Lett. 65(19) 2462–2465 URL https://link.aps.org/doi/10.1103/PhysRevLett.65.2462
  * [68] Tarat S and Majumdar P 2015 Eur. Phys. Jour. B 88 68 URL https://doi.org/10.1140/epjb/e2015-50284-6
  * [69] Swain N, Tiwari R and Majumdar P 2016 Phys. Rev. B 94(15) 155119 URL https://link.aps.org/doi/10.1103/PhysRevB.94.155119
  * [70] Patel N D, Mukherjee A, Kaushal N, Moreo A and Dagotto E 2017 Phys. Rev. Lett. 119(8) 086601 URL https://link.aps.org/doi/10.1103/PhysRevLett.119.086601
  * [71] Mayr M, Alvarez G, Sen C and Dagotto E 2005 Phys. Rev. Lett. 94(21) 217001 URL https://link.aps.org/doi/10.1103/PhysRevLett.94.217001
  * [72] Sensarma R, Randeria M and Trivedi N 2007 Phys. Rev. Lett. 98(2) 027004 URL https://link.aps.org/doi/10.1103/PhysRevLett.98.027004
  * [73] Taglieber M, Voigt A C, Aoki T, Hänsch T W and Dieckmann K 2008 Phys. Rev. Lett. 100(1) 010401 URL https://link.aps.org/doi/10.1103/PhysRevLett.100.010401
  * [74] Voigt A C, Taglieber M, Costa L, Aoki T, Wieser W, Hänsch T W and Dieckmann K 2009 Phys. Rev. Lett. 102(2) 020405 URL https://link.aps.org/doi/10.1103/PhysRevLett.102.020405
  * [75] Costa L, Brachmann J, Voigt A C, Hahn C, Taglieber M, Hänsch T W and Dieckmann K 2010 Phys. Rev. Lett. 105(12) 123201 URL https://link.aps.org/doi/10.1103/PhysRevLett.105.123201
  * [76] Naik D, Trenkwalder A, Kohstall C, Spiegelhalder F M, Zaccanti M, Hendl G, Schreck F, Grimm R, Hanna T M and Julienne P S 2011 The European Physical Journal D 65 55–65 ISSN 1434-6079 URL https://doi.org/10.1140/epjd/e2010-10591-2
  * [77] Wille E, Spiegelhalder F M, Kerner G, Naik D, Trenkwalder A, Hendl G, Schreck F, Grimm R, Tiecke T G, Walraven J T M, Kokkelmans S J J M F, Tiesinga E and Julienne P S 2008 Phys. Rev. Lett. 100(5) 053201 URL https://link.aps.org/doi/10.1103/PhysRevLett.100.053201
  * [78] Lu M, Burdick N Q and Lev B L 2012 Phys. Rev. Lett. 108(21) 215301 URL https://link.aps.org/doi/10.1103/PhysRevLett.108.215301
  * [79] Frisch A, Aikawa K, Mark M, Ferlaino F, Berseneva E and Kotochigova S 2013 Phys. Rev. A 88(3) 032508 URL https://link.aps.org/doi/10.1103/PhysRevA.88.032508
  * [80] Braun J, Drut J E, Jahn T, Pospiech M and Roscher D 2014 Phys. Rev. A 89(5) 053613 URL https://link.aps.org/doi/10.1103/PhysRevA.89.053613
  * [81] Roscher D, Braun J and Drut J E 2015 Phys. Rev. A 91(5) 053611 URL https://link.aps.org/doi/10.1103/PhysRevA.82.013624
  * [82] Baarsma J E, Gubbels K B and Stoof H T C 2010 Phys. Rev. A 82(1) 013624 URL https://link.aps.org/doi/10.1103/PhysRevA.82.013624
  * [83] Gubbels K B, Baarsma J E and Stoof H T C 2009 Phys. Rev. Lett. 103(19) 195301 URL https://link.aps.org/doi/10.1103/PhysRevLett.103.195301
  * [84] Hanai R, Kashimura T, Watanabe R, Inotani D and Ohashi Y 2013 Phys. Rev. A 88(5) 053621 URL https://link.aps.org/doi/10.1103/PhysRevA.88.053621
  * [85] Hanai R and Ohashi Y 2014 Phys. Rev. A 90(4) 043622 URL https://link.aps.org/doi/10.1103/PhysRevA.90.043622
  * [86] Guo H, Chien C C, Chen Q, He Y and Levin K 2009 Phys. Rev. A 80(1) 011601 URL https://link.aps.org/doi/10.1103/PhysRevA.80.011601
  * [87] Batrouni G G, Wolak M J, Hubert F and Rousseau V G 2009 EPL (Europhysics Letters) 86 47006 URL https://doi.org/10.1209
  * [88] Braun J, Drut J E and Roscher D 2015 Phys. Rev. Lett. 114(5) 050404 URL https://link.aps.org/doi/10.1103/PhysRevLett.114.050404
  * [89] Karmakar M 2018 Phys. Rev. A 97(3) 033617 URL https://link.aps.org/doi/10.1103/PhysRevA.97.033617