# Communal pairing in spin-imbalanced Fermi gases

D.C.W. Foo1 T.M. Whitehead1 G.J. Conduit1

(1 _Cavendish Laboratory, J.J. Thomson Avenue, Cambridge, CB3 0HE, United
Kingdom_
)

###### Abstract

A spin-imbalanced Fermi gas with an attractive contact interaction forms a
superconducting state whose underlying components are superpositions of Cooper
pairs that share minority-spin fermions. This superconducting state includes
correlations between all available fermions, making it energetically favorable
to the Fulde–Ferrell–Larkin–Ovchinnikov superconducting state. The ratio of
the number of up- and down-spin fermions in the instability is set by the
ratio of the up- and down-spin density of states in momentum at the Fermi
surfaces, to fully utilize the accessible fermions. We present analytical and
complementary Diffusion Monte Carlo results for the state.

For over a century, the phenomenon of superconductivity has captured the
attention of theorists, who have provided fundamental revelations about its
underlying principles. Bardeen, Cooper, and Schrieffer (BCS) [1, 2] gave
seminal insights into the mechanism of superconductivity in systems with spin-
balanced Fermi surfaces, showing that superconductivity may be understood as
the collective behavior of coherent Cooper pairs of fermions. Fulde, Ferrell,
Larkin, and Ovchinnikov (FFLO) [3, 4] extended this result, demonstrating that
even in systems with spin-imbalanced Fermi surfaces Cooper pairs may still
form the basis of a superconducting state. However, in spin-imbalanced
systems, the density of states in momentum at the Fermi surface of the
majority-spin fermions is greater than that of the minority-spin fermions, so
the number of Cooper pairs that can exist is limited by the number of
minority-spin fermions at their Fermi surface. This leaves many of the
majority-spin fermions unpaired and so uncorrelated, wasting their potential
for contributing correlation energy to the system.

In a few-fermion context, an instability [5] containing more majority than
minority-spin fermions maximizes the binding energy captured in spin-
imbalanced systems by taking advantage of the correlations between all
available momentum states. Such an instability involves non-exclusive pairing
between several majority-spin fermions and one minority-spin fermion in an
ensemble that we call a communal state. An example state is shown in Fig. 1,
with three majority- (up-)spin fermions each paired with the same minority-
(down-)spin fermion. This inspires us to merge Cooper pairs that share a
minority-spin fermion to construct a communal superconducting state that
correlates all available momentum states on the Fermi surfaces. We show that
this superconducting state with fermions shared between pairs is energetically
favorable over the exclusive Cooper pair-based FFLO superconductivity in spin-
imbalanced systems.

Figure 1: (Color online) Idealized representation of the spin-imbalanced
system showing Fermi surfaces for the down- (light blue circle) and up-spin
(light red arc) species, with occupiable momentum states extending over a
momentum scale set by the Debye frequency, forming annuli. The intensity of
color in the annuli indicate the approximate extent of the superconducting
correlations. Also shown are the momenta of
$(N_{\uparrow},N_{\downarrow})=(3,1)$ up- and down-spin fermions with
corresponding $q$-vectors $\mathbf{q}_{\sigma i}$. The angular spread of the
up-spin fermion momenta is exaggerated for clarity.

Current experimental developments enable the study of exotic superconducting
phases in solid-state spin-imbalanced Fermi gases [6, 7, 8, 9, 10], and spin-
orbit coupling may give rise to inhomogeneous superconductivity [11, 12]; but
no single experiment has provided unambiguous evidence for the existence of
FFLO superconductivity, leaving the true nature of the ground state an open
question. However, the recent development of uniform trapping potentials for
ultracold atomic gases [13] promises unparalleled experimental accuracy,
presenting an ideal opportunity to revisit the structure of the
superconducting ground state of spin-imbalanced Fermi gases.

In this Letter we examine the ratio of number of majority- to minority-spin
fermions in communal states underlying the superconducting state of a spin-
imbalanced Fermi gas, compare the energetics of communal superconductivity
favorably to that of FFLO superconductivity, and present analytical and
complementary Diffusion Monte Carlo results for the state. We also discuss
unique experimental consequences of the proposed communal superconductor.

To explore communal superconductivity we examine a two-spin fermionic system
with an attractive contact interaction. The quantum partition function,
$\mathcal{Z}=\int\mathcal{D}(\psi,\bar{\psi})\mathrm{e}^{-S[\psi,\bar{\psi}]}$,
depends on the BCS action

$\displaystyle S[\psi,\bar{\psi}]=\\!$
$\displaystyle\sum_{\omega,\mathbf{k},\sigma}\\!\bar{\psi}_{\mathbf{k},\sigma}\\!(-i\omega\\!+\\!\xi_{\mathbf{k},\sigma})\psi_{\mathbf{k},\sigma}$
$\displaystyle-g\\!\\!\\!\sum_{\omega,\mathbf{k},\mathbf{k}^{\prime},\mathbf{q}}\\!\\!\\!\bar{\psi}_{\mathbf{k},\uparrow}\bar{\psi}_{\mathbf{q}-\mathbf{k},\downarrow}\psi_{\mathbf{q}-\mathbf{k}^{\prime},\downarrow}\psi_{\mathbf{k}^{\prime},\uparrow},$

where $\psi_{\mathbf{k},\sigma}$ and $\bar{\psi}_{\mathbf{k},\sigma}$ are a
fermion field and its Grassmann conjugate, for momentum $\mathbf{k}$ and spin
species $\sigma\in\\{\uparrow,\downarrow\\}$,
$\xi_{\mathbf{k},\sigma}\equiv\epsilon_{\mathbf{k},\sigma}-\mu_{\sigma}$ where
$\epsilon_{\mathbf{k},\sigma}$ and $\mu_{\sigma}$ are the species-dependent
dispersion and chemical potential respectively, $g>0$ is the strength of the
attractive contact interaction, and $\omega$ is a fermionic Matsubara
frequency. In this expression the momenta $\mathbf{q}$, referred to henceforth
as $q$-vectors, give the net momenta of coupled fermions. Our strategy is to
build on the original BCS and FFLO theories that are directly applicable to
the solid state, and so here adopt a Debye frequency cutoff on the sums over
$\mathbf{k}$, however similar results are obtained in cold atom gases provided
proper regularization is carried out.

We perform a Hubbard-Stratonovich decoupling in the Cooper channel, using a
concise matrix formalism to express the action as

$\displaystyle S[\psi,\\!\Delta]=$
$\displaystyle\sum_{\omega,\mathbf{k}}\begin{pmatrix}\bar{\boldsymbol{\psi}}_{\uparrow}\\\
\boldsymbol{\psi}_{\downarrow}\end{pmatrix}^{\mathrm{T}}\\!\begin{pmatrix}\boldsymbol{\mathcal{G}}^{-1}_{\uparrow}&&-\boldsymbol{\Delta}\\\
-\boldsymbol{\Delta}^{\dagger}&&\boldsymbol{\mathcal{G}}^{-1}_{\downarrow}\end{pmatrix}\\!\begin{pmatrix}\boldsymbol{\psi}_{\uparrow}\\\
\bar{\boldsymbol{\psi}}_{\downarrow}\end{pmatrix}\\!$
$\displaystyle+\\!\sum_{\omega}\\!\frac{\mathrm{Tr}(\boldsymbol{\Delta}^{\dagger}\\!\boldsymbol{\Delta})}{g},$
(1)

where the vectors $\boldsymbol{\psi}_{\sigma}=(\psi_{(\mathbf{q}_{\sigma
1}+\varsigma_{\sigma}\mathbf{k}),\sigma},\psi_{(\mathbf{q}_{\sigma
2}+\varsigma_{\sigma}\mathbf{k}),\sigma},\ldots)^{\mathrm{T}}$, with
$\varsigma_{\uparrow}=+1$ and $\varsigma_{\downarrow}=-1$, the Grassman
conjugates $\bar{\boldsymbol{\psi}}_{\sigma}$ are similar, the matrices
$\boldsymbol{\mathcal{G}}^{-1}_{\sigma}=\mathrm{diag}(\mathcal{G}^{-1}_{(\mathbf{q}_{\sigma
1}+\varsigma_{\sigma}\mathbf{k}),\sigma},\mathcal{G}^{-1}_{(\mathbf{q}_{\sigma
2}+\varsigma_{\sigma}\mathbf{k}),\sigma},\ldots)$, for
$\mathcal{G}^{-1}_{\mathbf{p},\sigma}=-i\omega+\varsigma_{\sigma}\xi_{\mathbf{p},\sigma}$,
and

$\displaystyle\boldsymbol{\Delta}=\left(\begin{matrix}\Delta_{\mathbf{q}_{\uparrow
1}+\mathbf{q}_{\downarrow 1}}&&\Delta_{\mathbf{q}_{\uparrow
1}+\mathbf{q}_{\downarrow 2}}&&\cdots\\\ \Delta_{\mathbf{q}_{\uparrow
2}+\mathbf{q}_{\downarrow 1}}&&\Delta_{\mathbf{q}_{\uparrow
2}+\mathbf{q}_{\downarrow 2}}&&\cdots\\\
\vdots&&\vdots&&\ddots\end{matrix}\right),$

where the $\mathbf{q}_{\sigma i}$ run over all the $q$-vectors of species
$\sigma$. We label the number of fermions in the underlying instability, and
hence the number of $q$-vectors, per species by $N_{\sigma}$: therefore the
$\boldsymbol{\mathcal{G}}^{-1}_{\sigma}$ are $N_{\sigma}\times N_{\sigma}$
matrices and $\boldsymbol{\Delta}$ is an $N_{\downarrow}\times N_{\uparrow}$
matrix. We shall find that in spin-imbalanced systems $N_{\uparrow}\neq
N_{\downarrow}$, so that $\boldsymbol{\Delta}$ is rectangular rather than
square, and different numbers of fermions from each species are involved in
the underlying instability. In the system represented by Fig. 1, where there
are three up-spin and one down-spin fermions involved in the underlying
instability, $\boldsymbol{\Delta}$ would be a $1\times 3$ matrix. We focus our
analysis on the Cooper channel as recent work [14] shows that in the BCS
limit, screening [15] and pairing mechanisms decouple so the reduction in
critical temperature due to particle-hole interactions for both our communal
state and for FFLO will be the same. We note for completeness that decoupling
through the magnetic channel was also considered but had no consequence.

The elements of the $\boldsymbol{\Delta}$ matrix gap the dispersion. For the
Fulde–Ferrell (FF) state [3] (also referred to as single-plane-wave
superconductivity) $\boldsymbol{\Delta}$ has only a single entry, and for
crystalline FFLO superconductivity it is diagonal, see [16, 17, 18] and
references therein. The non-diagonal form here allows communal
superconductivity, as in common with the few-fermion analysis [5] multiple
majority-spin fermions share a minority-spin fermion. We focus on
superconductivity where the sharing majority-spin fermions have nearly aligned
$q$-vectors, comparable to the FF state. For simplicity of analysis we assume
that none of the $\mathbf{q}_{\uparrow i}+\mathbf{q}_{\downarrow j}$ pairs of
$q$-vectors in $\boldsymbol{\Delta}$ are degenerate. Following Ref. [5] the
$\mathbf{q}_{\sigma i}$ vectors are taken to be not equal to each other so
that each $\psi_{\mathbf{p},\sigma}$ appears only once in Eq. (1). This
corresponds to assigning states on the Fermi surfaces into non-overlapping
communal states of equal angular width.

With this expression for the action, working in the mean-field approximation,
we can carry out a Ginzburg-Landau expansion of the regularized thermodynamic
potential to obtain

$\displaystyle\Omega=T\sum_{\omega,\mathbf{k}}\sum_{n=1}^{\infty}\frac{1}{n}\mathrm{Tr}(\boldsymbol{\mathcal{G}}_{\uparrow}\boldsymbol{\Delta}\boldsymbol{\mathcal{G}}_{\downarrow}\boldsymbol{\Delta}^{\dagger})^{n}+\frac{\mathrm{Tr}(\boldsymbol{\Delta}^{\dagger}\boldsymbol{\Delta})}{g},$

where $T$ is the temperature. To make progress with this expression, we
symmetrize the coupling amplitudes, $\Delta_{\mathbf{q}}=\Delta$. Near the
second-order transition to the normal state we may neglect high-order terms in
$\Delta$ and truncate the expression for the thermodynamic potential to

$\displaystyle\Omega=\alpha\Delta^{2}+\frac{1}{2}\beta\Delta^{4}+\ldots,$

where

$\displaystyle\alpha$
$\displaystyle=\sum_{\mathbf{q}_{\uparrow},\mathbf{q}_{\downarrow}}\left(\frac{1}{g}+T\sum_{\omega,\mathbf{k}}\mathcal{G}_{\mathbf{q}_{\uparrow}+\mathbf{k},\uparrow}\mathcal{G}_{\mathbf{q}_{\downarrow}-\mathbf{k},\downarrow}\right),$
$\displaystyle\beta$
$\displaystyle=\\!\\!\sum_{\begin{subarray}{c}\mathbf{q}_{\uparrow
1},\mathbf{q}_{\downarrow 1},\\\ \mathbf{q}_{\uparrow
2}\,\mathbf{q}_{\downarrow 2}\end{subarray}}J(\mathbf{q}_{\uparrow
1},\mathbf{q}_{\downarrow 1},\mathbf{q}_{\uparrow 2},\mathbf{q}_{\downarrow
2}),$ (2)

with

$\displaystyle J(\mathbf{q}_{1},\mathbf{q}_{2},$
$\displaystyle\mathbf{q}_{3},\mathbf{q}_{4})\\!$
$\displaystyle=T\sum_{\omega,\mathbf{k}}\mathcal{G}_{\mathbf{q}_{1}+\mathbf{k},\uparrow}\mathcal{G}_{\mathbf{q}_{2}-\mathbf{k},\downarrow}\mathcal{G}_{\mathbf{q}_{3}+\mathbf{k},\uparrow}\mathcal{G}_{\mathbf{q}_{4}-\mathbf{k},\downarrow}.$
(3)

To evaluate these expressions, we specialize to the case of small Debye
frequency, found for many conventional superconductors [19, 20, 21]. In this
limit, the vectors $\mathbf{q}_{\sigma i}$ are expected to be approximately
parallel to maximise the number of contributing occupiable momentum states.
Approximately parallel but unequal $\mathbf{q}_{\sigma i}$ vectors provide a
natural tiling of the Fermi surfaces into non-overlapping communal states.
This enables us to factorize out combinatorial factors, giving

$\displaystyle\alpha=N_{\uparrow}N_{\downarrow}\Bigg{(}\frac{1}{g}$
$\displaystyle+T\sum_{\omega,\mathbf{k}}\mathcal{G}_{\mathbf{q}+\mathbf{k},\uparrow}\mathcal{G}_{\mathbf{q}-\mathbf{k},\downarrow}\Bigg{)},$
$\displaystyle\beta=N_{\uparrow}N_{\downarrow}\>\\!\big{[}J_{0}$
$\displaystyle+(N_{\uparrow}-1)J_{\uparrow}+(N_{\downarrow}-1)J_{\downarrow}$
$\displaystyle+(N_{\uparrow}-1)(N_{\downarrow}-1)J_{\uparrow\downarrow}\big{]},$
(4)

where $J_{0}=J(\mathbf{q},\mathbf{q},\mathbf{q},\mathbf{q})$,
$J_{\uparrow}=J(\mathbf{q}+\delta\mathbf{q}_{\uparrow},\mathbf{q},\mathbf{q}-\delta\mathbf{q}_{\uparrow},\mathbf{q})$,
$J_{\downarrow}=J(\mathbf{q},\mathbf{q}+\delta\mathbf{q}_{\downarrow},\mathbf{q},\mathbf{q}-\delta\mathbf{q}_{\downarrow})$,
and
$J_{\uparrow\downarrow}=J(\mathbf{q}+\delta\mathbf{q}_{\uparrow},\mathbf{q}+\delta\mathbf{q}_{\downarrow},\mathbf{q}-\delta\mathbf{q}_{\uparrow},\mathbf{q}-\delta\mathbf{q}_{\downarrow})$.
Here $\mathbf{q}$ is taken to represent the average $q$-vector for the
fermions, symmetrized between species, and $\delta\mathbf{q}_{\sigma}$ is half
the average separation between $q$-vectors for species $\sigma$, which in the
small Debye frequency limit is orthogonal to the vector $\mathbf{q}$. We
follow the prescription of Ref. [5] that the angular widths of the regions of
Fermi surface involved in the communal superconducting state are the same
between species, so the arc lengths $\delta\mathbf{q}_{\sigma}$ are
proportional to the Fermi momenta and
$|\delta\mathbf{q}_{\uparrow}|/|\delta\mathbf{q}_{\downarrow}|=k_{\uparrow\mathrm{F}}/k_{\downarrow\mathrm{F}}$,
where $k_{\sigma\mathrm{F}}$ is the Fermi momentum of species $\sigma$. For a
free dispersion $J_{\uparrow\downarrow}$ may be evaluated at zero temperature
as

$\displaystyle
J_{\uparrow\downarrow}=\frac{\mathcal{N}_{d}Q^{2}}{(Q^{4}\\!-\\!k_{\perp}^{2}(\delta\mathbf{q}_{\uparrow}\\!+\\!\delta\mathbf{q}_{\downarrow})^{2})(Q^{4}\\!-\\!k_{\perp}^{2}(\delta\mathbf{q}_{\uparrow}\\!-\\!\delta\mathbf{q}_{\downarrow})^{2})},$
(5)

where $\mathcal{N}_{d}$ is a dimension $d\in\\{2,3\\}$ dependent normalization
factor, $Q^{2}\equiv
2k_{\mathrm{F}}k_{\mathrm{D}}+\tfrac{1}{2}|\delta\mathbf{q}_{\uparrow}|^{2}+\tfrac{1}{2}|\delta\mathbf{q}_{\downarrow}|^{2}+k_{\mathrm{D}}^{2}+k_{\perp}^{2}$,
$k_{\mathrm{F}}=(k_{\uparrow\mathrm{F}}+k_{\downarrow\mathrm{F}})/2$ is the
average Fermi momentum, $k_{\mathrm{D}}$ is the Debye frequency and
$k_{\perp}$ is the average extent of $\mathbf{k}$ in the direction
perpendicular to $\mathbf{q}$ such that the fermions are within the Debye
frequency of the Fermi energy. Similar expressions for $J_{\uparrow}$,
$J_{\downarrow}$ and $J_{0}$ can be found by taking
$\delta\mathbf{q}_{\downarrow}=0$, $\delta\mathbf{q}_{\uparrow}=0$, or both.
For $k_{\mathrm{D}}\ll k_{\mathrm{F}}$,
$k_{\perp}\sim\sqrt{k_{\mathrm{F}}k_{\mathrm{D}}}$ and Eq. (5) confirms that
for a single instability $\beta\geq 0$ for realistic values of the
$\mathbf{q}_{i}$ [17], justifying the truncation in Eq. (Communal pairing in
spin-imbalanced Fermi gases).

To identify the optimal ratio of number of fermions involved in the communal
superconductor, we express Eq. (Communal pairing in spin-imbalanced Fermi
gases) as a function of $N_{\uparrow}/N_{\downarrow}$ and
$N_{\uparrow}N_{\downarrow}$, and then perform a saddle point analysis to
optimize $\Omega$ with respect to $N_{\uparrow}/N_{\downarrow}$,
$N_{\uparrow}N_{\downarrow}$, and $\Delta$ simultaneously. We note from Ref.
[5] that singly-excited state fluctuations in $N_{\uparrow}/N_{\downarrow}$
have a linear energy dependence, so we focus on this saddle point analysis.
This gives the expected ratio of number of fermions involved in the underlying
instability as

$\displaystyle\frac{N_{\uparrow}}{N_{\downarrow}}=\frac{J_{\uparrow\downarrow}-J_{\downarrow}}{J_{\uparrow\downarrow}-J_{\uparrow}}=\left(\frac{|\delta\mathbf{q}_{\uparrow}|}{|\delta\mathbf{q}_{\downarrow}|}\right)^{2}=\left(\frac{\nu_{\uparrow}}{\nu_{\downarrow}}\right)^{2/(d-1)},$
(6)

where the second equality was obtained from Eq. (5) and $\nu_{\sigma}$ is the
density of states in momentum at the Fermi surface of species $\sigma$.

In the case of cold atoms in the absence of a Debye frequency, regularization
of the divergent momentum summation in $\alpha$ can be done by using
scattering theory to replace the weak interaction coupling parameter $g$ with
a formal expression involving the scattering length and another summation with
the same UV divergence characteristics, resulting in a convergent expression
[22]. The momentum summation in $\beta$, and indeed in all higher terms of the
expansion of $\Omega$, do not exhibit UV divergences, as can be seen in how
the expression for $J_{\uparrow\downarrow}$ in Eq. (5) is finite for any value
of $k_{\mathrm{D}}$, indicating that the high $\mathbf{k}$ contributions to
the summation are not dominant. Indeed, provided the interaction is weak, the
optimal placement of the $\mathbf{q}$ are unchanged and since the dominant
contributions to the $J$ come from around the Fermi surface, the result of Eq.
(6) remains unchanged, that is making the shift from solid-state to ultracold
atoms does not change any of the significant physics.

This result confirms that the superconducting state is indeed communal, with
pairs sharing minority-spin fermions to take advantage of all available
correlations in spin-imbalanced systems. Eq. (6) also aligns with our
heuristic expectation that the instability involves more fermions of the
species with the larger density of states in momentum at its Fermi surface, as
was also found in the few-fermion case [5]. For spin-balanced systems,
$\nu_{\uparrow}=\nu_{\downarrow}$ and so $N_{\uparrow}/N_{\downarrow}=1$,
recovering the BCS theory result, whilst in the polaron limit of a single
minority-spin impurity in a full Fermi sea of majority-spin fermions, the
single minority-spin fermion couples with all the majority-spin fermions at
their Fermi surface, in agreement with results from the literature [23, 24,
25].

Our conclusions do not contradict the well-known result that FFLO
superconductivity is the ground state in one dimension, obtained separately by
density matrix renormalization group [26] and time-evolving block decimation
methods [27]. In a one-dimensional system $\nu_{\uparrow}=\nu_{\downarrow}$
regardless of the spin-imbalance, removing our initial physical motivation for
communal pairing. Furthermore, in one dimension $\delta\mathbf{q}_{\sigma}$
cannot be orthogonal to $\mathbf{q}$, and so if one attempts to form a
communal state, it necessarily involves fermions of a particular spin species
with different energies, invalidating the symmetrization of the gap amplitudes
$\Delta_{\mathbf{q}}$ so there is no evidence that communal pairing can be
energetically favourable to the FFLO state found in previous studies of one-
dimensional systems.

The same optimization procedure that gave Eq. (6) for
$N_{\uparrow}/N_{\downarrow}$ also provides an expression for
$N_{\uparrow}N_{\downarrow}$, as

$\displaystyle
N_{\uparrow}N_{\downarrow}=\frac{\left(J_{0}-J_{\uparrow}-J_{\downarrow}+J_{\uparrow\downarrow}\right)^{2}}{\left(J_{\uparrow\downarrow}-J_{\uparrow}\right)\left(J_{\uparrow\downarrow}-J_{\downarrow}\right)}.$
(7)

For reasonable values of the $|\delta\mathbf{q}_{\sigma}|$ and $|\mathbf{q}|$
in spin-imbalanced systems this expression gives values of
$N_{\uparrow}N_{\downarrow}>1$, confirming that the communal superconductor is
indeed made up of multiple fermions of at least one spin species. Excessively
high $N_{\uparrow}N_{\downarrow}$ is energetically penalized by the highest
term in the expansion of the thermodynamic potential, which goes as
$(N_{\uparrow}N_{\downarrow}\Delta^{2})^{n}$, and so we expect communal
superconductivity to have both $N_{\sigma}$ being reasonably small integers.
In the spin-balanced limit Eq. (7) collapses to the BCS result
$N_{\uparrow}N_{\downarrow}=1$.

The analysis may be adapted to the number conserving canonical ensemble by
constructing the Helmholtz free energy
$F=\Omega+\mu_{\uparrow}n_{\uparrow}+\mu_{\downarrow}n_{\downarrow}$, where
$n_{\sigma}=-\partial\Omega/\partial\mu_{\sigma}$ is the total particle number
of spin species $\sigma$, so
$F=(1-\mu_{\uparrow}\partial/\partial\mu_{\uparrow}-\mu_{\downarrow}\partial/\partial\mu_{\downarrow})\Omega$.
As the chemical potential only appears in the propagator, it suffices to note
that
$\partial\mathcal{G}_{\mathbf{k},\sigma}/\partial\mu_{\nu}=\varsigma_{\sigma}\mathcal{G}^{2}_{\mathbf{k},\sigma}\delta_{\sigma,\nu}$
and so the net result on $N_{\uparrow}N_{\downarrow}$ and
$N_{\uparrow}/N_{\downarrow}$ of moving from the grand canonical ensemble to
the canonical ensemble is to shift the functions $J$ as

$\displaystyle J\to J-$ $\displaystyle
T\sum_{\omega,\mathbf{k}}\mathcal{G}_{\mathbf{q}_{1}+\mathbf{k},\uparrow}\mathcal{G}_{\mathbf{q}_{2}-\mathbf{k},\downarrow}\mathcal{G}_{\mathbf{q}_{3}+\mathbf{k},\uparrow}\mathcal{G}_{\mathbf{q}_{4}-\mathbf{k},\downarrow}$
$\displaystyle\times(\mu_{\uparrow}(\mathcal{G}_{\mathbf{q}_{1}+\mathbf{k},\uparrow}+\mathcal{G}_{\mathbf{q}_{3}+\mathbf{k},\uparrow})$
$\displaystyle\quad-\mu_{\downarrow}(\mathcal{G}_{\mathbf{q}_{2}-\mathbf{k},\downarrow}+\mathcal{G}_{\mathbf{q}_{4}-\mathbf{k},\downarrow})),$

leading to the same conclusions as in the grand canonical ensemble, namely
that in the canonical ensemble the superconducting state is indeed communal
with $N_{\uparrow}N_{\downarrow}>1$ and
$N_{\uparrow}/N_{\downarrow}=(\nu_{\uparrow}/\nu_{\downarrow})^{(2/(d-1))}$.

Now that we have shown that the communal superconductor is energetically
favourable over single-plane-wave FFLO superconductivity in a spin-imbalanced
Fermi gas, we need to confirm that we have not compromised the stability of
the superconducting state. We validate this by examining the phase boundaries
between the communal superconductor and three competitor phases.

With increasing spin-imbalance, BCS superconductivity becomes unstable against
FFLO superconductivity at the Chandrasekhar-Clogston limit [28, 29]. Although
communal superconductivity is energetically favourable over FFLO
superconductivity, BCS superconductivity still has a large density-of-states
advantage over communal superconductivity, and to a first approximation the
phase boundary between BCS superconductivity and communal superconductivity
will remain at the same Chandrasekhar-Clogston value.

The phase boundary between communal superconductivity and the normal state
will also remain the same as between single-plane-wave FFLO superconductivity
and the normal state. In both cases the second-order phase transition occurs
when $\alpha=0$, and this condition is identical between FFLO and communal
superconductivity, up to an irrelevant multiplicative factor of
$N_{\uparrow}N_{\downarrow}$ in Eq. (4), and so the phase boundary is also
identical.

Stability against phase separation can be expressed as the positive-
definiteness of the total particle number susceptibility matrix [30, 31]. This
condition includes the possibility of separation into two superconducting
phases, with ratios of number of fermions differing from that predicted in Eq.
(6), and may be expanded following Eq. (Communal pairing in spin-imbalanced
Fermi gases), to leading order in $\Delta$ giving

$\displaystyle\alpha\frac{\partial^{2}\alpha}{\partial
q^{2}}>2\left(\frac{\partial\alpha}{\partial q}\right)^{2}.$

This is the same as the equivalent expression for FFLO superconductivity, up
to a factor of $(N_{\uparrow}N_{\downarrow})^{2}$ that cancels between the two
sides of the inequality. This indicates that the line of stability against
phase separation is the same for communal superconductivity as for FFLO
superconductivity to leading order.

Although the discussion above focuses on nearly aligned $q$-vectors,
comparable to the FF state, it is known that the Larkin–Ovchinnikov (LO) state
[4] built from two plane-waves can be energetically favorable to single-plane-
wave superconductivity. Therefore, we now follow the prescription of Larkin
and Ovchinnikov and consider a communal superconducting state out of two
instabilities on opposite sides of the Fermi surfaces. The only differences in
the theory of communal superconductivity for one and two instabilities are a
multiplicative factor of 2 in Eq. (Communal pairing in spin-imbalanced Fermi
gases) and additional terms in the expression for $\beta$ in Eq. (4).
Similarly to the single instability case the optimal instability contains more
up- than down-spin fermions, and so the communal superconducting state is also
energetically favorable over the LO state.

Figure 2: (Color online) Left: Plot of condensate fraction in pair momentum
space in units of $1/2\sqrt{3}\pi r_{s}$. The positions of expected peaks
corresponding to a communal state ($r_{s}q=1/\sqrt{3}\pi$), FFLO pairing
($r_{s}q=1/2\pi$), and BCS pairing ($r_{s}q=0$) are shown with green stars,
black triangles and a magenta circle respectively. Communal type peaks are
prominent with no obvious FFLO peaks and a mild BCS peak. The k-space grid is
drawn in black. Right: Diagrams of pairing orbitals in the communal and FFLO
phases. The colored points indicate the filled Fermi areas of the down-spin
(blue) and up-spin (red) species in the non-interacting limit. The colored
clouds indicate which states contribute to the instabilities with intensity
indicating the strength of contribution. Black lines indicate a possible
pairing arrangement.

We supplement the preceding analysis with numerical evidence obtained from a
quantum Monte Carlo study of a finite spin-imbalanced 2D homogeneous fermion
gas with attractive interactions using the casino program [32]. To minimize
finite size effects [33, 34], we place the fermions in a rhomboidal box with
vertex angle $60^{\circ}$ so that the discretized momentum points form a
triangular lattice. This allows for the densest tiling of discrete momentum
points in 2D giving the closest to circular Fermi surfaces. We set up a system
of 61 spin-up and 19 spin-down fermions and work in atomic units with the
average inter-fermion separation $r_{s}=1$. Such particle numbers are
consistent with those used in other DMC studies [35, 36], giving us confidence
that the results obtained should be at least qualitatively related to the
analytics done in the thermodynamic limit. We note for completeness that
qualitatively similar results were obtained for systems with both a smaller
and larger number of fermions, and in systems on a square simulation cell.

As the Fermi surfaces form hexagons and have a 2:1 ratio of fermions, we
expect a $(N_{\uparrow},N_{\downarrow})=(2,1)$ instability. An
ultratransferable pseudopotential [37] with scattering length $a=6.5991$ and
zero effective range was introduced so that the BCS coherence length of an
equivalent spin-balanced system was approximately equal to the simulation cell
size. We note here for completeness that while such a pseudopotential is
indeed meant to emulate the scattering properties of a contact interaction, it
nevertheless has a finite extent in space that introduces a natural momentum
cutoff. The UTP therefore does not have infinite range in momentum space and
so the assumption of small Debye frequency used in the analytical derivation
is applicable to the simulation.

Following previous work [38], we employ a Slater-Jastrow trial wavefunction of
the form
$\Psi_{\mathrm{T}}=\mathrm{e}^{-J(\mathsf{r_{\uparrow}},\mathsf{r_{\downarrow}})}\det[\phi(\mathbf{s}_{i,j})].$
The pairing orbital is

$\displaystyle\phi(\mathbf{s}_{i,j})=$
$\displaystyle\sum_{l=1}^{4}a_{l}\cos(\mathbf{k}_{l}\cdot\mathbf{s}_{i,j})$
$\displaystyle+\Theta(N_{\downarrow}-i)(1-\tfrac{s}{r_{s}})^{3}\Theta(1-\tfrac{s}{r_{s}})\sum_{m=0}^{2}b_{m}r^{m},$

where
$\mathbf{s}_{i,j}\equiv\mathbf{r}_{\uparrow,i}-\Theta(N_{\downarrow}-i)\mathbf{r}_{\downarrow,j}$,
$s\equiv\lvert\mathbf{s}_{i,j}\rvert$, $\mathbf{r}_{\sigma,i}$ is the position
vector of the $i-$th fermion of spin species $\sigma$, $\mathbf{k}_{l}$ is the
l$th$ shortest reciprocal-space vector, $\Theta$ is the Heaviside step
function, and the $\\{a_{l}\\},\\{b_{m}\\}$ are optimisable parameters. The
Jastrow factor is given by

$\displaystyle
J(\mathsf{r_{\uparrow}},\mathsf{r_{\downarrow}})=\sum_{i,j}\Big{[}$
$\displaystyle\sum_{k}u_{k}r_{i,j}^{k-1}(1-\frac{r_{i,j}}{r_{s}})^{3}\Theta(1-\frac{r_{i,j}}{r_{s}})$
$\displaystyle+$
$\displaystyle\sum_{m}p_{m}\cos(\mathbf{G}_{m}\cdot\mathbf{r}_{i,j})+\nu\
\mathrm{terms}\Big{]}$

where $\mathsf{r_{\sigma}}$ denotes the set of position vectors for all
fermions of spin species $\sigma$,
$\mathbf{r}_{i,j}\equiv\mathbf{r}_{\uparrow,i}-\mathbf{r}_{\downarrow,j}$,
$r\equiv\lvert\mathbf{r}_{i,j}\rvert$, the $\mathbf{G}_{m}$ are reciprocal-
space vectors through which anisotropy may be introduced, and the
$\\{u_{k}\\},\\{p_{m}\\}$ are optimisable parameters.
$J(\mathsf{r_{\uparrow}},\mathsf{r_{\downarrow}})$ is thus a function of all
opposite-spin fermion separations containing a short range isotropic $u$ term,
anisotropic $p$ terms [39] and a $\nu$ term [40] that reflects the simulation
cell symmetry and whose form is omitted for brevity. While the trial
wavefunction was originally used to capture pairing between electrons and
holes in a bilayer, it nevertheless has three attractive properties that
warrant usage in this context; namely that it is a pairing wavefunction and
that it deforms continuously into the accepted form for a balanced
superconductor and the non-interacting state. Furthermore, the nodal surface
can vary and be optimized through the $b_{m}$ parameters.

The trial wavefunction was optimised using Variational Monte Carlo [41] then
equilibrated using Diffusion Monte Carlo before the condensate fraction in
momentum space [42] was accumulated. In casino, the condensate fraction is
defined as $f_{\mathbf{q}}\equiv\sum_{\mathbf{k}}(\langle
c^{\dagger}_{\mathbf{k},\uparrow}c^{\dagger}_{\mathbf{q}-\mathbf{k},\downarrow}c_{\mathbf{q}-\mathbf{k},\downarrow}c_{\mathbf{k},\uparrow}\rangle-
n_{\mathbf{k},\uparrow}n_{\mathbf{q}-\mathbf{k},\downarrow})$, where
$c^{\dagger}_{\mathbf{k},\sigma}(c_{\mathbf{k},\sigma})$ is the creation
(annihilation) operator for a fermion of momentum $\mathbf{k}$ and spin
$\sigma$ and $n_{\mathbf{k},\sigma}\equiv\langle
c^{\dagger}_{\mathbf{k},\sigma}c_{\mathbf{k},\sigma}\rangle$ is the momentum
density. The results display strong anisotropy of the condensate fraction at
$r_{s}q=1/\sqrt{3}\pi$, regardless of whether in the Jastrow factor the
$p_{m}$ are selected to permit anisotropy, constrained to ensure isotropy, or
set to zero. The state observed is therefore not an artefect of (an)isotropy
of the trial wavefunction, consistent with the communal state being the robust
ground state. Furthermore, we confirm the sharing of minority-spin fermions
between pairs by noting that the anisotropy is consistent with the pairing
scheme shown in the upper right of Fig. 2; the black lines indicate the up-
down pairs present in a single $(N_{\uparrow},N_{\downarrow})=(2,1)$
instability which when combined with its time-reversed partner gives 4 pairing
momenta, corresponding to 4 peaks in the condensate fraction. We only indicate
the peak positions for plane wave FFLO; crystalline FFLO would have peaks at
all $q$-vectors of equal magnitude to plane-wave FFLO, none of which coincide
with the 4 observed peaks of the communal pairing state.

The lack of a strong condensate fraction at $r_{s}q=1/2\pi$ indicates that the
FFLO phase is not the major contributor to the ground state. Finally we note
for completeness that while there appears to be a BCS contribution to the
condensate fraction, the uncertainty in the condensate fraction at $r_{s}q=0$
was twice as large as at other points and so that peak is half as significant
as the others. Similar results were seen in systems with a 3:1 and 4:1 ratio
at the Fermi surface, affirming the pairing result for
$N_{\uparrow}/N_{\downarrow}$, and for systems with the same 2:1 ratio at the
Fermi surfaces but different total numbers of fermions, namely 19 spin-up and
7 spin-down fermions, and 127 spin-up and 37 spin-down fermions. The results
obtained for systems with different numbers of fermions give us confidence
that the observed state is not an artefact of finite size effects.

We note for completeness that the condensate fraction was also gathered for
pairs of the same spin species in an attempt to simultaneously search for
induced p-wave superfluidity [43], however the values of intra-spin condensate
fraction were more than 10 orders of magnitude smaller than those for the
inter-spin condensate fraction and were indistinguishable from zero at all
pair momenta.

Having seen that communal superconductivity is energetically favorable over
FFLO superconductivity, we now consider its possible experimental
consequences. We focus on two effects where the communal nature of the
underlying instability should be directly observable.

Multiple phase transitions: With increasing spin-imbalance, that is increasing
$\nu_{\uparrow}/\nu_{\downarrow}$, Eq. (6) predicts that the ratio
$N_{\uparrow}/N_{\downarrow}$ should increase. Starting from the BCS state
with $\nu_{\uparrow}/\nu_{\downarrow}=1$ the system should progress through
several communal superconducting phases with increasing values of
$N_{\uparrow}/N_{\downarrow}$, where both $N_{\sigma}$ are integers, giving a
series of different superconducting states. The product
$N_{\uparrow}N_{\downarrow}$ governed by Eq. (7) restricts $N_{\uparrow}$ and
$N_{\downarrow}$ to small integers and so restricts the number of transitions
observed. Each transition is expected to be second order, and so the communal
superconducting phase would be characterized by a series of discontinuities in
the heat capacity and the compressibility, which is directly observable in
ultracold atomic gases [44], as the spin-imbalance is changed. No such phase
transitions are expected within the FFLO phase. Plots of the specific heat
$C_{\mathrm{V}}$ against $\nu_{\uparrow}/\nu_{\downarrow}$ at low temperature
are shown in Fig. 3 for the communal and FFLO phases. The plot used the same
parameters as our DMC study. A finite discontinuity in the specific heat is
present in both curves at the Chandrasekhar-Clogston limit at
$\nu_{\uparrow}/\nu_{\downarrow}\approx 1.2$, with further discontinuities
manifesting in the communal phase as $N_{\uparrow}$ and $N_{\downarrow}$
change. With increasing imbalance the heat capacity jumps upward as an
additional up-spin fermion enters the strongly correlated state, increasing
its binding energy. The behavior of the communal state is in marked contrast
to the FFLO phase that has near constant specific heat, giving a significant
observable difference between the two phases.

Figure 3: Plot of specific heat $C_{\mathrm{V}}$ against
$\nu_{\uparrow}/\nu_{\downarrow}$ for the energetically favourable communal
(solid) and not energetically favourable FFLO (dotted) phases. Communal phases
are labelled by the pair $(N_{\uparrow},N_{\downarrow})$. The BCS phase
denotes where the Chandrasekhar-Clogston limit is breached, indicated by the
vertical dashed line. The two curves coincide in the BCS and (1,1) phases.

Superconducting order parameter: In real space the order parameter will
exhibit a beat pattern due to the interference between similar $q$-vectors,
which could allow the identification of the particular $q$-vectors in the
superconductor. The order parameter and its spread in momentum could be
determined in an ultracold atomic gas experiment through density-density
correlations measured from time-of-flight experiments [45]. While the width of
the narrow peaks of either the BCS or FFLO phase can be resolved in absorption
imaging experiments to better than $0.1k_{\mathrm{F}}$ [46], the spread of
$q$-vectors from communal pairing seen in Fig. 2 and indicated by exact
diagonalisation studies [5] is on the order of $k_{\mathrm{F}}$ and should be
readily distinguished. This can be seen in Fig. 2 where alongside the clear
peaks, the condensate fraction is generally nonzero for momenta $r_{s}q\leq
1/\sqrt{3\pi}$ and only falls to zero on the border of the region shown. In
contrast, FFLO and crystalline FFLO theories predict sharp peaks in the
condensate fraction, as in spin-balanced BCS theory, at fixed magnitude of the
pairing momenta.

Andreev reflection: The elementary excitations above the proposed ground state
are predicted to be well-described by the few fermion analysis [5]. This
should have novel consequences especially concerning Andreev reflection
experiments as the strong correlations between a group of fermions held in a
communal state should result in multiple retroreflected holes for a single
incident fermion, in sharp contrast to the single hole per fermion expected in
normal FFLO theory.

Beyond these experimental signatures, communal superconductivity also
introduces the notion that the building block of a superconductor may involve
the sharing of fermions between Cooper pairs. In particular, there can be
fluctuations in the number of fermions in the underlying instability, which
could lead to the renormalization of the properties of a spin-balanced
superconductor. The analysis is also generalizable to systems with multiple
underlying instabilities, more akin to crystalline FFLO superconductivity
[16], and we use that system as a guide for the likely modifications when
communal superconductivity is built from several instabilities.

We have introduced the idea of a communal superconductor, whose underlying
instability is composed of multiple pairs with shared fermions, to enable the
use of all available inter-fermion correlations. We have shown that communal
superconductivity is energetically favorable over FFLO superconductivity in
spin-imbalanced Fermi gases, both analytically and with complementary DMC
results, and that a communal superconductor has clear experimental signatures.

Data used for this Letter are available online [47]. The authors thank Adam
Nahum, Johannes Hofmann, Johannes Knolle, Jens Paaske, Robin Reuvers, and
Victor Jouffrey for useful discussions, and acknowledge the financial support
of the NUS, the EPSRC, and the Royal Society.

## References

  * [1] J. Bardeen, L.N. Cooper, and J.R. Schrieffer, Phys. Rev. 106, 162 (1957).
  * [2] J. Bardeen, L.N. Cooper, and J.R. Schrieffer, Phys. Rev. 108, 1175 (1957).
  * [3] P. Fulde and R.A. Ferrell, Phys. Rev. 135, A550 (1964).
  * [4] A.I. Larkin, Y.N. Ovchinnikov, Sov. Phys. JETP 20, 762 (1965).
  * [5] T.M. Whitehead and G.J. Conduit, Phys. Rev. B 97, 014502 (2018).
  * [6] A. Bianchi, R. Movshovich, C. Capan, P.G. Pagliuso, and J.L. Sarrao, Phys. Rev. Lett. 91, 187004 (2003).
  * [7] H. Mayaffre, S. Krämer, M. Horvatić, C. Berthier, K. Miyagawa, K. Kanoda, and V.F. Mitrović, Nature Physics 10, 928 (2014).
  * [8] R. Prozorov, M.D. Vannette, S.A. Law, S.L. Bud’ko, and P.C. Canfield, Phys. Rev. B 77, 100503(R) (2008).
  * [9] P.C. Canfield, S.L. Bud’ko, and B.K. Cho, Physica C 262, 249 (1996).
  * [10] M.H. Hamidian, S.D. Edkins, S.H. Joo, A. Kostin, H. Eisaki, S. Uchida, M.J. Lawler, E.-A. Kim, A.P. Mackenzie, K. Fujita, J. Lee, and J.C.S. Davis, Nature 532, 343 (2016).
  * [11] W. Zhang and W. Yi, Nat. Commun. 4, 2711 (2013).
  * [12] S.-T. Lo, S.-W. Lin, Y.-T. Wang, S.-D. Lin, and C.-T. Liang, Sci. Rep. 4, 5438 (2014).
  * [13] B. Mukherjee, Z. Yan, P.B. Patel, Z. Hadzibabic, T. Yefsah, J. Struck, and M.W. Zwierlein, Phys. Rev. Lett. 118, 123401 (2017).
  * [14] L. Pisani, A. Perali, P. Pieri and G. Calvanese Strinati, Phys. Rev. B 97, 014528 (2018)
  * [15] L.P. Gor’kov and T.K. Melik-Barkhudarov, JETP 13, 1018 (1964)
  * [16] J.A. Bowers and K. Rajagopal, Phys. Rev. D 66, 065002 (2002).
  * [17] R. Casalbuoni and G. Nardulli, Rev. Mod. Phys. 76, 263 (2004).
  * [18] C. Mora and R. Combescot, EPL 66, 833 (2004).
  * [19] J.A. Kok and W.H. Keesom, Physica 4, 835 (1937).
  * [20] C.A. Reynolds, B. Serin, and L.B. Nesbitt, Phys. Rev. 84, 691 (1951).
  * [21] M. Horowitz, A.A. Silvidi, S.F. Malaker, and J.G. Daunt, Phys. Rev. 88, 1182 (1952).
  * [22] W. Ketterle and M. W. Zwierlein, arXiv:0801.2500 (2008).
  * [23] F. Chevy, Phys. Rev. A 74, 063628 (2006).
  * [24] R. Combescot and S. Giraud, Phys. Rev. Lett. 101, 050404 (2008).
  * [25] P. Massignan, M. Zaccanti, and G.M. Bruun, Rep. Prog. Phys. 77, 033401 (2014).
  * [26] F. Heidrich-Meisner, G. Orso, and A. E. Feiguin, Phys. Rev. A 81, 053602 (2010).
  * [27] H. Lu, L. O. Baksmaty, C. J. Bolech and H. Pu, Phys. Rev. Lett. 108, 225302 (2012).
  * [28] B.S. Chandrasekhar, Appl. Phys. Lett. 1, 7 (1962).
  * [29] A.M. Clogston, Phys. Rev. Lett. 9, 266 (1962).
  * [30] Q. Chen, Y. He, C.-C. Chien, and K. Levin, Phys. Rev. A 74, 063603 (2006).
  * [31] J. Wang, Y. Che, L. Zhang, and Q. Chen, Sci. Rep. 7, 39783 (2017).
  * [32] R. Needs, M. Towler, N. Drummond and P.L. Ríos, casino User’s Guide Version 2.13 (2015).
  * [33] N.D. Drummond, R.J. Needs, A. Sorouri and W.M.C. Foulkes, Phys. Rev. B 78, 125106 (2008).
  * [34] D. Frenkel, Eur. Phys. J. Plus (2013) 128: 10.
  * [35] A. Galea, H. Dawkins, S. Gandolfi, and A. Gezerlis, Phys. Rev. A 93, 023602 (2016).
  * [36] G. E. Astrakharchik, J. Boronat, J. Casulleras, and S. Giorgini, Phys. Rev. Lett 93, 200404 (2004).
  * [37] T.M. Whitehead, L.M. Schonenberg, N. Kongsuwan, R.J. Needs and G.J. Conduit, Phys. Rev. A 93, 042702 (2016).
  * [38] R. Maezono, P.L. Ríos, T. Ogawa and R.J. Needs, Phys. Rev. Lett. 110, 216407 (2013).
  * [39] N.D. Drummond, M.D. Towler and R.J. Needs, Phys. Rev. B 70, 235119 (2004)
  * [40] T.M. Whitehead, M.H. Michael and G.J. Conduit, Phys. Rev. B 94, 035157 (2016).
  * [41] P.R.C. Kent, R.J. Needs and G. Rajagopal, Phys. Rev. B 59, 12344 (1999)
  * [42] L.M. Schonenberg and G.J. Conduit, Phys. Rev. A 95, 013633 (2017)
  * [43] A. Bulgac, M. M. Forbes and A. Schwenk, Phys. Rev. Lett. 97, 020402 (2006).
  * [44] M.J.H. Ku, A.T. Sommer, L.W. Cheuk, and M.W. Zwierlein, Science 335, 563 (2012).
  * [45] E. Altman, E. Demler, and M.D. Lukin, Phys. Rev. A 70, 013603 (2004).
  * [46] M. G. Ries, A. N. Wenz, G. Zürn, L. Bayha, I. Boettcher, D. Kedar, P. A. Murthy, M. Neidig, T. Lompe, and S. Jochim, Phys. Rev. Lett. 114, 230401 (2015).
  * [47] D.C.W. Foo, T.M. Whitehead, and G.J. Conduit, Cambridge University DSpace repository.


# Electron liquid state in the symmetric Anderson lattice

Igor N. Karnaukhov  G.V. Kurdyumov Institute for Metal Physics, 36 Vernadsky
Boulevard, 03142 Kiev, Ukraine

###### Abstract

Using mean field approach, we provide analytical and numerical solution of the
symmetric Anderson lattice for arbitrary dimension at half filling. The
symmetric Anderson lattice is equivalent to the Kondo lattice, which makes it
possible to study the behavior of an electron liquid in the Kondo lattice. We
have shown that, due to hybridization (through an effective field due to
localized electrons) of electrons with different spins and momenta k and
$\textbf{k}+\overrightarrow{\pi}$, the gap in the electron spectrum opens at
half filling. Such hybridization breaks the conservation of the total magnetic
momentum of electrons, the spontaneous symmetry is broken. The state of
electron liquid is characterized by a large Fermi surface. A gap in the
spectrum is calculated depending on the magnitude of the on-site Coulomb
repulsion and value of s-d hybridization for the chain, as well as for square
and cubic lattices. Anomalous behavior of the heat capacity at low
temperatures in the gapped state, which is realized in the symmetric Anderson
lattice, was also found.

## Introduction

The processes of electron scattering with spin flip play a dominant role in
the Kondo problem, they are determined by the exchange s-d interaction [1, 2].
They should also be explicitly taken into account when considering the
behavior of the Kondo lattice, since at half-filling the insulator state is
also determined by this interaction. The one-particle modifications of the
Kondo lattice (see for example [3]), exact solvable many-particle models [4,
5], which do not taken into account the electron scattering with spin flip
unable to describe the Kondo insulator state. Numerical calculations of the
Anderson and Kondo lattices take into account clusters with a small number of
particles, which significantly affects the adequacy of the results obtained
[6, 7, 8]. This leads to the fact that the mechanism of the formation of the
large Fermi surface and the gap in the electronic spectrum in the Kondo
lattice at half filling is still unclear. In contrast to the single-impurity
Anderson and Kondo models in their lattice version, the processes of electron
scattering with spin flip cannot be calculated exactly.

In the weak s-d hybridization limit the Hamiltonian of the Anderson lattice
reduces to that of the Kondo lattice. It should also be taken into account,
that the local density of localized electrons must be equal to unity. In the
weak s-d hybridization limit, the local density of d-electrons in the
symmetric Anderson lattice is equal to unity. d-electrons determine the local
spin-$\frac{1}{2}$ at the sites of the lattice, only in this case we can talk
about the spin-$\frac{1}{2}$ Kondo lattice. Thus the symmetric Anderson model
is similar to the Kondo lattice. In the strong coupling limit of the Hubbard
model, when $U\gg t$ ( $U$ and $t$ are on-site repulsion and hopping
integral), it is shown that the Hubbard model and the Kondo lattice model
become identical [9]. It gives possibility to use the formalism, proposed for
calculation of the Mott transition in the Hubbard model [10], for solution of
the Kondo lattice problem.

In the paper, we consider the solution of the symmetric Anderson model, using
a mean field approach. An effective $\lambda$-field connects the states of
d-electrons with the different spins and momenta k,
$\textbf{k}+\overrightarrow{\pi}$. Due to the s-d hybridization, the states of
s-electrons with different spins and momenta k,
$\textbf{k}+\overrightarrow{\pi}$ also hybridize, the gap opens in the
electron spectrum at half filling. The value of the gap is determined by the
magnitude of the $\lambda$-field (which in turn depends on the on-site
repulsion) and the value of the s-d hybridization. We shall show that breaking
spontaneous symmetry in the Anderson lattice Hamiltonian makes it possible to
take into account the scattering of conduction electrons on d-electrons with a
spin flip. They determine the ground state of symmetric Anderson and Kondo
lattices, the value of the gap in the electron spectrum, cause formation of an
insulator state at half-filling. The electron spectrum is symmetric about the
zero energy, it corresponds to the symmetric Anderson model for a nontrivial
solution of the $\lambda$-field. The gap in the electron spectrum has a many-
particle nature.

## Model

The Hamiltonian of the Anderson lattice is the sum of two terms, the first of
which is determined by energy of the bands of s- and d-electrons and
hybridization between them, the second takes into account the on-site
repulsion of d-electrons ${\cal H}={\cal H}_{0}+{\cal H}_{int}$

$\displaystyle{\cal
H}_{0}=-\sum_{j=1}^{N-1}\sum_{\sigma=\uparrow,\downarrow}(c^{\dagger}_{j,\sigma}c_{j+1,\sigma}+c^{\dagger}_{j+1,\sigma}c_{j,\sigma})+\epsilon_{g}\sum_{j=1}^{N}\sum_{\sigma=\uparrow,\downarrow}n_{j,\sigma}+v\sum_{j=1}^{N}\sum_{\sigma=\uparrow,\downarrow}(c^{\dagger}_{j,\sigma}d_{j,\sigma}+d^{\dagger}_{j,\sigma}c_{j,\sigma}),$
$\displaystyle{\cal H}_{int}=U\sum_{j=1}^{N}n_{j,\uparrow}n_{j,\downarrow},$
(1)

where $c^{\dagger}_{j,\sigma},c_{j,\sigma}$ and
$d^{\dagger}_{j,\sigma},d_{j,\sigma}$ ($\sigma=\uparrow,\downarrow)$ are the
fermion operators determined on a lattice site $j$, $U$ is the value of the
on-site Hubbard interaction determined by the density operator
$n_{j,\sigma}=d^{\dagger}_{j,\sigma}d_{j,\sigma}$, the band width of
c-fermions is determined by the hopping integral equal to one, the energy of
flat band of d-fermions equal to $\epsilon_{g}$, $v$ defines the
hybridizations of s- and d-electrons, N is the total number of atoms.

First of all, we will focus on the consideration of the model for single
Anderson impurity, exact solution of which has been obtained by Wiegmann [1].
At $\epsilon_{g}=-\frac{U}{2}$ and $U\gg\Gamma$ ($\Gamma=v^{2}$) trickly one
electron is localized on the impurity, so we can talk about an impurity with
spin-$\frac{1}{2}$. This case corresponds to the symmetric Anderson model.
According to [2] the behavior of impurity in the symmetric Anderson model is
equivalent to the spin-$\frac{1}{2}$ Kondo impurity. In this context, in a
non-magnetic (paramagnetic) state the behavior of an electron liquid in the
symmetric Anderson lattice is similar to that in the spin-$\frac{1}{2}$ Kondo
lattice. Should be notes also, when the hybridization is small the Hamiltonian
(1) can be mapped into the Kondo lattice model with the exchange imtegral
$J\simeq-\frac{2\Gamma U}{\epsilon_{g}(\epsilon_{g}+U)}$.

Figure 1: (Color online) Electron spectrum of the chain a) and low energy
excitations b) (where $\Delta=0.004$) as function of wave vector calculated at
$\epsilon_{g}=-1$ (or $U=2$), $v=0.1$, $\lambda=0.2$.  Figure 2: (Color
online) Electron spectrum of the square lattice a) and low energy excitations
b) (where $\Delta=0.004$) as function of the wave vector calculated at
$\epsilon_{g}=-1$ (or $U=2$), $v=0.1$, $\lambda=0.2$.  Figure 3: (Color
online) The gap $\Delta$ in the electron spectrum as function of $v$ and
$\lambda$ calculated at $\epsilon_{g}=-1$ (or $U=2$), in the coordinates $v$,
$\lambda$ (the value of $\Delta$ does not depend on the dimension of the
lattice).

## The ground-state

$\lambda$-field breaks the spontaneous symmetry [11, 12], the magnetic moments
of conduction and d-electrons are not conserved, the total number of electrons
is conserved [10]. Only for $\textbf{q}=\overrightarrow{\pi}$, the electron
spectrum corresponds to the symmetric Anderson lattice. Thus, we are talking
about the lattice with the doubled period, as in the case of the Peierls
phase-transition. In this case the electron spectrum is symmetric with respect
to zero and has the following form for the chain (see in Fig 1a) and square
lattice (see in Fig 2a). Four branches of the spectrum $\pm
E_{\gamma}(\textbf{k)}$ ($\gamma=1,2$) are determined by the following
expression (see section Methods))

$E_{\gamma}(\textbf{k})=\frac{1}{\sqrt{2}}\sqrt{\alpha(\textbf{k})+(-1)^{\gamma}\sqrt{\alpha^{2}(\textbf{k})-4\beta^{2}(\textbf{k})}},$
(2)

where
$\alpha(\textbf{k})=\varepsilon^{2}(\textbf{k})+\epsilon_{g}^{2}+\lambda^{2}+2\Gamma$,
$\beta^{2}(\textbf{k})=(\lambda^{2}+\epsilon_{g}^{2})\varepsilon^{2}(\textbf{k})-2\epsilon_{g}\Gamma\varepsilon(\textbf{k})+\Gamma^{2}$.

For $v\neq 0$, $\lambda\neq 0$ and an arbitrary dimension of the lattice the
spectrum is gapped at half-filling. The gap is equal to zero at $v=0$ or
$\lambda=0$ (see in Fig 3), $v=\lambda=0$ corresponds to the atomic limit of
the symmetric Anderson lattice with a bare on-site repulsion.

The dependence of the gap on $v$ and $\lambda$ has quite universal in nature,
its value does not depend on the dimension of the model (see in Fig 3). As
shown below, the value of $\lambda$ is the solution of Eq (3), at given $U$
its value is determined by $v$, depends on dimension of the lattice. In the
chain the low energy spectrum is determined by two special points (see in Fig
1b). In the square lattice, they form a square in the k-plane (see in Fig 2b).
The k-surface, formed by the points on the k-plane of the low-energy spectrum,
is constructed as follows: each projection of the wave vector corresponds to
two values of its other projections. The $\gamma$-branches of the spectrum are
hybridized, the electron density is defined as $\sum_{\gamma=1,2}\int
d\textbf{k}=n_{s}+1$, here $n_{s}$ is the density of s-electrons ($n_{s}=1$ at
half-filling). Thus, we can talk about a large Fermi surface in the Kondo
lattice also, the behavior of which is similar to the symmetric Anderson
lattice.

In the ground state the equation for $\lambda$ follows from the form of the
action (6) and the spectrum (2)

$\frac{\lambda}{U}=\frac{1}{4}\sum_{\gamma=1,2}\int
d\textbf{k}\frac{1}{E_{\gamma}(\textbf{k})}\frac{\partial
E_{\gamma}^{2}(\textbf{k})}{\partial\lambda},$ (3)

where the value of the on-site repulsion corresponds to a condition
$U=-2\epsilon_{g}$. The energy of the quasi-particle excitations depend on
$\epsilon_{g}$ (2), thus (3) is self-consistent equation. As a result, the
state of the symmetric Anderson lattice can be realized for $\lambda\neq 0$.

As we noted above the gap in the spectrum depends both $v$ and $\lambda$ (in
the Anderson model the scattering matrix depends on $v$ and $U$ [1]).
Numerical solutions for $\lambda$ as function of $v$ are presented in Fig 4
for different dimension of the lattice. Results of calculations are obtained
at $\epsilon_{g}=-1$, $U=2$ (see in Fig 4a)) and $\epsilon_{g}=-1/2$, $U=1$
(see in Fig 4b)). Trivial solution $\lambda=0$ and $v=0$ corresponds to the
symmetric Anderson model with $\epsilon_{g}=-\frac{U}{2}$ in which s- and
d-electrons are not coupled. Nontrivial solution for $\lambda$ and $v$
determines state of electron liquid in the symmetric Anderson model. Solution
for $\lambda$ takes place at finite values of $v$, the $\lambda$-value
increases with increasing $v$ (see in Fig 4). Such at $v=0.45$ the gaps in the
spectrum have the following values for U=2 (U=1): $\Delta=0.015$
($\Delta=0.093$) in the chain, $\Delta=0.018$ ($\Delta=0.074$) in square and
$\Delta=0.005$ ($\Delta=0.088$) in cubic lattices. The value of the gap
decreases with increasing $U$.

The stability of gapped ground state also depends on the magnitude of the on-
site repulsion. Comparing actions (6) at T = 0K for the gapped (paramagnetic)
state and the phase state with the maximum magnetic moment
$M=M_{s}+M_{d}=\frac{N}{2}+\frac{N}{2}$, we determine the region of stable
insulator state. The minimum value of the on-site Hubbard repulsion $U_{c}$,
at which the gapped state is realized, is $U_{c}=0.999205$ (1D),
$U_{c}=0.99933$ (2D), $U_{c}=0.99914$ (3D) at $v=\frac{1}{2}$ and
$U_{c}=0.96987$ (1D), $U_{c}=0.98969$ (2D), $U_{c}=0.97796$ (3D) at $v=1$. The
value of $U_{c}$ weakly depends on the dimension of the lattice and the
magnitude of the s-d hydridization $v$. The non-trivial solution for
$\lambda$, shown in Fig. 4, corresponds to stable insulator state at $U=2$,
the gapped state at $U=1$, is unstable. Calculating the actions $\delta
S=\frac{S(\lambda)}{\beta}-\frac{S(0)}{\beta}$ at T = 0K, which correspond to
the states with the gap (for $\lambda\neq 0$) and gapless with a maximum
magnetic moment (for $\lambda=0$), we determine the stability of the insulator
state at $U=2$. The numerical calculations are shown in Fig 5 for different
dimensions of the lattice. The region of stability of the insulator state is
highlighted in color, the solution for $\lambda$ obtained for $U=2$ (see Fig.
4a) is in the colored region. Note, that $\delta S=0$ at $\lambda=0$ and
$\delta S\sim\lambda^{2}$ at $\lambda\to\infty$, thus intermediate values of
$0<\lambda<0.1$ correspond to the stable insulator state (see in Fig 5).

Taking into account the spectrum of quasi particle excitations (2) in the
action (6), we determine the low temperature behavior of the specific heat
(its electronic part) for 3D system (cubic lattice). In the insulator state,
the value of $C$ exponentially depends on the value of the gap. Numerical
calculations of the temperature dependence $\frac{C}{T}$ are shown in Fig. 6.
The calculations are made for $U=2$ and $v=0.43$, $\lambda=0.0075$,
$\Delta=0.0023$ (see Fig. 6a) and $v=0.5$, $\lambda=0.07$, $\Delta=0.0279$
(see Fig.6b). Note that for these parameters, the gaps in the electron
spectrum differ by an order of magnitude. The experimentally measured heat
capacity $\frac{C}{T}$ $SmB_{6}$ [14] shows the anomalous behavior in the low
temperature range $2K<T<10K$, namely at T = 4K there is a minimum. According
to numerical calculations (see insets in Fig. 6), a local minimum of
$\frac{C}{T}$ is realized at low temperatures. The temperature range and the
minimum value depend on the value of the gap.

Figure 4: (Color online) $\lambda$ as a function of $v$ calculated at
$\epsilon_{g}=-1$, $U=2$ a) and $\epsilon_{g}=-0.5$, $U=1$ b) for the chain,
square and cubic lattices.  Figure 5: (Color online) $\delta S\times 10^{5}$
as a function of $v$ and $\lambda$ calculated at $\epsilon_{g}=-1$ ($U=2$) for
the chain a), square b) and cubic c) lattices, the region of stability of the
gapped state is highlighted in color.

Figure 6: (Color online) The temperature dependence of the specific heat
$\frac{C}{T}$ calculated for cubic lattice at $U=2$ : $v=0.43$,
$\lambda=0.0075$, $\Delta=0.0023$ a) and $v=0.5$, $\lambda=0.07$,
$\Delta=0.0279$. Inserts illustrate the regions of abnormal behavior of the
specific heat at low temperatures.

## Methods

### Canonical functional of the symmetric Anderson model

Let us introduce the operator
$\chi_{j}^{\dagger}=d^{\dagger}_{j,\uparrow}d_{j,\downarrow}$ and redefine the
term ${\cal H}_{int}$ (1) is the following form ${\cal
H}_{int}=-U\sum_{j}\chi_{j}^{\dagger}\chi_{j}$ [10]. The Hubbard-Stratonovich
transformation converts the interacting problem into a non-interacting one in
a stochastic external field (hereinafter we will define it as the
$\lambda$-field). We define the interaction term, taking into account the
action ${S}_{0}$

$\displaystyle
S=S_{0}+\sum_{j}\frac{\lambda^{\ast}_{j}\lambda_{j}}{U}+\sum_{j}(\lambda_{j}\chi_{j}+H.c.)$
(4)

The canonical functional is determined as

${\cal Z}=\int{\cal D}[\lambda]\int{\cal D}[\chi\dagger,\chi]e^{-S},$

where the action
$S=\frac{1}{U}\sum_{j}\lambda^{\ast}_{j}\lambda_{j}+\int_{0}^{\beta}d\tau\sum_{\textbf{k}}\Psi_{\textbf{k}}^{\dagger}(\tau)[\partial_{\tau}+{\cal
H}_{eff}(\textbf{k})]\Psi_{\textbf{k}}(\tau)$, $\Psi_{\textbf{k}}(\tau)$ is
the wave function, k is the wave vector of an electron. We expect that
$\lambda_{j}$ does not depend on $\tau$, since the translation invariance is
conserved in an electron liquid state.

At the on-site hybridization between d-states of electrons with different
spins and due to translation invariance of the model Hamiltonian, only the
phases of $\lambda_{j}$ depends on $j$, a namely $\lambda_{j}=\exp(i\textbf{q
}\textbf{j})\lambda$ [10], where q is an unknown wave vector. The phase can
also fluctuate at the lattice site; subsequent averaging over local
fluctuations restores translational invariance [13]. The task is reduced to
moving fermions in a static inhomogeneous $\lambda$-field that sets the form
of ${\cal H}_{eff}(\textbf{k})$

${\cal
H}_{eff}(\textbf{k})=\left(\begin{array}[]{cccc}-\varepsilon(\textbf{k})&v&0&0\\\
v&\epsilon_{g}&\lambda&0\\\ 0&\lambda^{*}&-\epsilon_{g}&v\\\
0&0&v&-\varepsilon(\textbf{k}+\textbf{q})\end{array}\right)$ (5)

The spectrum of non-interacting s-electrons is given by
$\varepsilon(\textbf{k})=-2\sum_{i=1}^{D}\cos k_{i}$, here D is dimension of
the model.

We can integrate out fermions to obtain the following action $S$ per an atom

$\displaystyle\frac{S(\lambda)}{\beta}=-\frac{{T}}{{N}}\sum_{\textbf{k}}\sum_{n}\sum_{\gamma=1}^{4}\ln[-i\omega_{n}+E_{\gamma}(\textbf{k},\textbf{q})]+\frac{|\lambda|^{2}}{{U}},$
(6)

where $\omega_{n}=T(2n+1)\pi$ are the Matsubara frequencies, k, q are the
momenta of electrons, four quasiparticle excitations
$E_{\gamma}(\textbf{k},\textbf{q})$ ($\gamma=1,...,4$) determine the electron
states in the $\lambda$-field. In the saddle point approximation the canonical
functional ${\cal Z}$ will be dominated by the minimal action $S$ (6), that
satisfies the following equation $\partial S/\partial\lambda=0$. In the Kondo
problem, the processes of electron scattering from spin flip dominate, in our
case, the $\lambda$-field connects the states of s-electrons with opposite
spins and different momenta k and $\textbf{k}+\textbf{q}$. Due to the on-site
repulsion between electrons in the Hubbard model[10], an effective field
connects states of electrons with different spins and different momenta k,
$\textbf{k}+\overrightarrow{\pi}$. This state is stable at half-filling, a gap
opens in the electron spectrum [10]. In this case, the same mechanism of
coupling between s-electrons is realized indirectly through an intermediate
subsystem of d-electrons. The gap in the electron spectrum will naturally be
smaller.

## Conclusion

Using mean field approximation we have considered the solution of the
symmetric Anderson lattice at half-filling for different dimensions of the
lattice. It is shown that an effective field, which binds the states of
s-electrons with different spins and momenta, leads to the gap in the electron
spectrum. The electron spectrum is mirror symmetric (with respect to zero
energy), has the type of the Majorana spectrum. The states of s- and
d-electrons are hybridized, therefore the Fermi surface is determined by the
total density of electrons. In the symmetric Anderson lattice model the gaped
state of electron liquid corresponds to breaking spontaneous symmetry. This
made it possible to take into account the processes of spin flip scattering of
conduction electrons on localized electrons. This is important because it is
these processes that lead to the appearance of the Abrikosov - Suhl resonance
in the Kondo problem and an insulator state in the Kondo lattice. The gapped
state is formed at finite values of the s-d hybridization and on-site
repulsion. Numerical calculations of the low-temperature behavior of the heat
capacity are carried out in the gapped state. The results explained the
presence of a local minimum of the heat capacity $\frac{C}{T}$ at low
temperatures observed in $SMB_{6}$. In the symmetric Anderson lattice, a local
spin-$\frac{1}{2}$ is realized at lattice site, therefore its behavior is
similar to the spin-$\frac{1}{2}$ Kondo lattice. The proposed approach allows
us to describe the symmetric Anderson and Kondo lattices in one formalism. One
can speak of a large Fermi surface and the gapped state of an electron liquid
in the Kondo lattice at half filling.

## References

  * [1] Wiegmann P.B. Towards an exact solution of the Anderson model. _Phys.Lett.A_ 80, 163 (1980)
  * [2] Tsvelick A.M. & Wiegmann P.B. Exact results in the theory of magnetic alloys. _Advances in Physics_ 32, 453 (1983)
  * [3] Yang W.-W., Zhao J., Luo H.-G., & Zhong Y. Exactly solvable Kondo lattice model in the anisotropic limit. _Phys. Rev. B_ 100, 045148 (2019)
  * [4] Karnaukhov I.N. Exactly solvable model of a one-dimensional Kondo lattice. _Phys.Rev.B_ 56, R4313 (1997)
  * [5] Karnaukhov I.N. & Slieptsov I.O. Exactly solvable 2D topological Kondo lattice model. _Europhysics Letters_ 109 57005 (2015)
  * [6] Guerrero M. & Yu C.C. Kondo insulators modeled by the one-dimensional Anderson lattice: A numerical-renormalization-group study. _Phys. Rev. B_ 51 10301 (1995)
  * [7] Sykora S. & Becker K.W. Heavy fermion properties of the Kondo Lattice model. _Scientific Reports_ 3 2691 (2013)
  * [8] Wang Y.-O., Lin H.-Q., & Gubernatis J.E. Zero temperaturenumerical studies of multiband lattice models of strongle correlated electrons. _Commun.Comput.Phys._ 1 575 (2006)
  * [9] Ivantsov I., Ferraz A. & Kochetov E. Hubbard model vs. Kondo model: Strong coupling limit, arXiv:2010.12255 [cond-mat.str-el]
  * [10] Karnaukhov I.N. Mott transition in two-band fermion model with on-site Coulomb repulsion. _Annals of Phys._ 433 168308 (2020)
  * [11] Karnaukhov I.N. Spontaneous breaking of time-reversal symmetry in topological insulators. _Phys.Lett.A_ 381 1967 (2017)
  * [12] Karnaukhov I.N. Spontaneous breaking of timereversal symmetry in topological superconductors. _Scientific Reports_ 7 7008 (2017)
  * [13] Golosov D,I. Onsite density matrix description of the extended Falicov-Kimball model at finite temperatures. _Phys.Rev.B_ 101 165130 (2020)
  * [14] Phelan W.A., Koohpayeh S.M., Cottingham P., Freeland J.M., Leiner J.C., Broholm C.L., & McQueen T.M. Correlation between bulk thermodynamic measurements and the low-temperature-resistance plateau in SmB6. _Physical Review X_ , 4 031012 (2014).

## Author contributions statement

I.K. is an author of the manuscript

## Additional information

The author declares no competing financial interests.