# Perturbative computation of thermal characteristics of the Stoner phase
transition

Oskar Grocholski$~\\!{}^{1}$ and Piotr H. Chankowski$~\\!{}^{2}$111Emails:
<EMAIL_ADDRESS><EMAIL_ADDRESS>
${}^{1}\phantom{a}$Deutsches Electronen-Synchrotron DESY
Notkestr. 85, 22607 Hamburg, Germany
2Faculty of Physics, University of Warsaw,
Pasteura 5, 02-093 Warsaw, Poland

###### Abstract

We apply the thermal (imaginary time) perturbative expansion to the relevant
effective field theory to compute characteristics of the phase transition to
the ordered state which can occur at low temperatures in the gas of
(nonrelativistic) spin $1/2$ fermions interacting through a short range spin
independent repulsive binary interaction potential. We show how to obtain a
systematic expansion of the system’s free energy depending on the densities
$n_{+}$ and $n_{-}$ of spin up and spin down fermions. In this paper we
truncate this expansion at the second order and determine, by numerically
minimizing the free energy, the equilibrium proportions of $n_{+}$ and $n_{-}$
(that is, the system’s polarization) as functions of the temperature, the
system’s overall density $n=n_{+}+n_{-}$ and the strength of the interaction.

Keywords: Diluted gas of interacting fermions, phase transitions,
thermodynamic potentials, effective field theory, scattering lengths

##  1 Introduction

There is a qualitative argument that in the gas of spin 1/2 fermions
interacting through a short range repulsive spin-independent binary potential
a phase transition to the ordered state should occur, if the interaction
strength and/or the system’s density is sufficiently large. Indeed, at zero
temperature, when the entropy factor does not intervene, the configuration of
the system in which there are more fermions in one spin state than in the
other one may be energetically favoured because, due to the Pauli exclusion
principle the $s$-wave interaction of fermions in the same spin state is
impossible and the resulting decrease of the interaction energy may be greater
than the associated increase of the kinetic energy (increase of the Fermi
energy of the more populated spin state).

Theoretical investigation of this phenomenon, called the Stoner transition,
taking into account its temperature dependence, requires the full machinery of
statistical mechanics. The standard textbook treatment of the problem [1, 2],
equivalent to the so-called mean field approach or the Hartree-Fock
approximation, employs the pseudo-potential method which allows to determine
in the first order approximation the Hamiltonian spectrum and to compute the
Canonical Ensemble partition function of the system. In this approximation the
phase transition is continuous (with the divergent magnetic susceptibility
characterized by the critical exponent $\gamma=1$ and a finite discontinuity
of the heat capacity) and at low temperatures (where the Sommerfeld expansion
can be used to obtain analytical expression for the relevant chemical
potentials) it occurs when [3, 1, 2]

$\displaystyle k_{\rm F}a_{0}\geq{\pi\over 2}\left[1+{\pi^{2}\over
12}\left({k_{\rm B}T\over\varepsilon_{\rm F}}\right)^{2}+\dots\right],$

where the (overall) Fermi wave vector and energy

$\displaystyle k_{\rm F}=\left(3\pi^{2}~{}\\!{N\over
V}\right)^{1/3},\phantom{aaaa}\varepsilon_{\rm F}={\hbar^{2}k_{\rm F}^{2}\over
2m_{f}}~{}\\!,$ (1)

characterize the density of the system and $a_{0}>0$ is the $s$-wave
scattering length characterizing the strength of the (repulsive) interaction.
The continuous character of the Stoner transition obtained in this
approximation is, however, accidental - it is due to a numerical coincidence
specific for a (three-dimensional) system of spin $s=1/2$ fermions only (in
the same approximation the transition is of first order if $s>1/2$ and/or
$D\neq 3$). In fact, computing the system’s free energy beyond the mean field
approximation, using the ordinary second order perturbative expansion, it was
found [4] that at low temperatures it is of the first order, just as had been
suggested in [5] on the basis of the generic presence of non-analytic terms
(resulting from the coupling of the order parameter to the gap-less modes) in
the free energy which cause the transition to have the first order character.

The character of the considered transition (its dependence on the parameter
$k_{\rm F}a_{0}$) can be most easily investigated at zero temperature because
then the problem reduces to the computation of the ground-state energy density
$E_{\Omega}/V$ of the system of fermions interacting through a binary spin-
independent repulsive potential as a function of the system’s density $n=N/V$
and its polarization

$\displaystyle P=(N_{+}-N_{-})/N~{}\\!.$ (2)

Such a computation is most easily performed using the modern effective field
theory approach, the application of which to this problem has been pioneered
in [6]. In this approach the underlying spatially nonlocal field-theory
interaction (see e.g. [7] for the exposition of the relevant formalism of the
second quantization), resulting from the ordinary potential two-body
interaction, is replaced by an infinite set of local (contact) effective
interactions

$\displaystyle\hat{V}_{\rm
int}=C_{0}\\!\int\\!d^{3}\mathbf{x}~{}\\!\psi_{+}^{\dagger}\psi_{+}\psi_{-}^{\dagger}\psi_{-}+\hat{V}_{\rm
int}^{(C_{2})}+V_{\rm int}^{(C_{2}^{\prime})}+\dots$ (3)

$\psi_{\pm}(\mathbf{x})$ are here the usual field operators of spin up and
down fermions; the terms $V_{\rm int}^{(C_{2})}$, $V_{\rm
int}^{(C_{2}^{\prime})}$ (which will be not needed in this work) represent
local operators of lower length dimension with four fermionic fields and two
spatial derivatives and the ellipsis stands for other local operators (with
more derivatives and/or field operators) of yet lower length dimension (see
[6]). The amount of work needed to obtain the systematic expansion of the
ground state energy in powers (which can be modified by logarithms) of $k_{\rm
F}R$, where $R$ is the characteristic length scale of the underlying two-body
spin independent interaction potential, is in this way greatly reduced. This
is because in this approach the coupling constants, like $C_{0}$ in (3), of
the effective local interactions are directly determined in terms of the
scattering lengths $a_{\ell}$ and effective radii $r_{\ell}$,
$\ell=0,1,\dots,$ (which are assumed to be of order $\sim R$) parametrizing
the low energy expansion in powers of the relative momentum
$\hbar|\mathbf{k}|$ of the elastic scattering amplitude of two fermions. The
simplifications brought in by the effective field theory method allowed to
easily reproduce [8] and generalize to arbitrary repulsive potentials and
arbitrary spins $s$ [9] the old result of Kanno [10] who computed the order
$(k_{\rm F}R)^{2}$ correction to the energy density using the specific hard
sphere interaction of spin $s=1/2$ fermions. The first order character of the
phase transition at $T=0$ is then clearly seen in the form of the energy
density obtained in this approximation plotted as a function of the order
parameter $P$: starting from some value of $k_{\rm F}a_{0}$ the energy density
develops the second minimum well away from the one at $P=0$ and at $k_{\rm
F}a_{0}=1.054$ (for $s=1/2$) this second minimum becomes deeper than that at
$P=0$.

However the analysis of the dependence on the order parameter of the system’s
energy density which includes the complete order $(k_{\rm F}R)^{3}$
corrections obtained recently in [11, 12] using the same effective field
theory approach shows that, independently of the value $s$ of the fermion
spin, they have the effect of erasing the first order character of the Stoner
transition, making it almost indistinguishable from the continuous one. This
is reflected in the fact that the height of the hill separating the minimum at
$P\neq 0$ from the one at $P=0$ is greatly reduced (for higher spins also the
position of the nontrivial minimum of $E_{\Omega}/V$ as a function of the
relevant order parameter is strongly shifted towards its small values)
compared to the situation without these corrections. Moreover, there are
claims [13] based on a resummation of an infinite subclass of Feynman diagrams
contributing to the ground-state energy density that the transition (at $T=0$)
is indeed continuous. Although it is not obvious that the contributions taken
into account in this resummation are really the dominant ones [12], the
results it leads to seem to agree well, as far as the critical value of
$k_{\rm F}a_{0}$ is concerned, with the numerical quantum Monte Carlo
simulations [14].

In view of this situation it is desirable to investigate how the higher order
corrections influence the character of the Stoner phase transition at nonzero
temperatures. With this goal in mind in this paper we formulate a systematic
perturbative expansion of the thermodynamic potentials of the system in
question applying the standard imaginary time formalism [7] within the
effective field theory. We show that the expansion of the free energy is in
this approach particularly simple being given by the same connected vacuum
Feynman diagrams which give nonzero contributions to the energy density
expressed in terms of the chemical potentials of the noninteracting system. In
the numerical analysis we restrict ourselves in this paper only to the second
order contributions reproducing the results obtained in [4], but with more
labour the computations can be extended to higher orders as well.

##  2 Perturbative expansion of the thermodynamic potential
$\Omega(T,V,\mu_{+},\mu_{-})$

The natural equilibrium statistical physics formalism in which to treat the
problem of the gas of fermions the interactions of which preserve their spins,
and therefore the numbers $N_{\sigma}$ of particles with the spin projection
$\sigma$, is the Grand Canonical Ensemble with separate chemical potentials
$\mu_{\sigma}$ associated with the individual spin projections. One is
therefore interested in the statistical operator (as usually, $\beta\equiv
1/k_{\rm B}T$)

$\displaystyle\hat{\rho}={1\over\Xi_{\rm
stat}}~{}\\!e^{-\beta\hat{K}}~{}\\!,\phantom{aaa}{\rm
in~{}which}\phantom{aaa}\hat{K}=\hat{H}_{0}-\sum_{\sigma}\mu_{\sigma}\hat{N}_{\sigma}+\hat{V}_{\rm
int}\equiv\hat{K}_{0}+\hat{V}_{\rm int}~{}\\!,$ (4)

and in computing the statistical sum (we specify the notation to the case of
spin $1/2$ fermions, so that $\sigma=+,-$)

$\displaystyle\Xi_{\rm stat}(T,V,\mu_{+},\mu_{-})={\rm
Tr}\\!\left(e^{-\beta\hat{K}}\right),$ (5)

from which all the necessary thermodynamic potentials can in principle be
obtained by performing the standard steps. The free part $\hat{K}_{0}$ of the
operator $\hat{K}=\hat{K}_{0}+\hat{V}_{\rm int}$, where $\hat{V}_{\rm int}$
will be taken in the form (3), is

$\displaystyle\hat{K}_{0}=\sum_{\mathbf{p},\sigma}\left(\varepsilon_{\mathbf{p}}-\mu_{\sigma}\right)a_{\mathbf{p},\sigma}^{\dagger}a_{\mathbf{p},\sigma}=\sum_{\sigma}\int\\!{d^{3}\mathbf{p}\over(2\pi)^{3}}\left(\varepsilon_{\mathbf{p}}-\mu_{\sigma}\right)a^{\dagger}_{\sigma}(\mathbf{p})~{}\\!a_{\sigma}(\mathbf{p})~{}\\!,$
(6)

with $\varepsilon_{\mathbf{p}}\equiv\hbar^{2}\mathbf{p}^{2}/2m_{f}$, in the
normalizations in the finite volume $V$ and in an infinite space,
respectively.

To compute perturbatively the statistical sum $\Xi_{\rm
stat}(T,V,\mu_{+},\mu_{-})$ one introduces [7] the (imaginary time)
interaction picture evolution operator

$\displaystyle{\cal
U}_{I}(\tau_{2},\tau_{1})=e^{\tau_{2}\hat{K}_{0}}~{}\\!e^{-(\tau_{2}-\tau_{1})\hat{K}}~{}\\!e^{-\tau_{1}\hat{K}_{0}}~{}\\!,$
(7)

which satisfies the differential equation

$\displaystyle{d\over d\tau_{2}}~{}\\!{\cal
U}_{I}(\tau_{2},\tau_{1})=-V^{I}_{\rm int}(\tau_{2})~{}\\!{\cal
U}_{I}(\tau_{2},\tau_{1})~{}\\!,$

($V^{I}_{\rm int}(\tau_{2})=e^{\tau_{2}\hat{K}_{0}}V_{\rm
int}e^{-\tau_{2}\hat{K}_{0}}$) with the “initial” condition ${\cal
U}_{I}(\tau,\tau)=\hat{1}$ and which formally can be written in the form

$\displaystyle{\cal U}_{I}(\tau_{2},\tau_{1})={\rm
T}_{\tau}\exp\\!\left\\{-\int_{\tau_{1}}^{\tau_{2}}\\!d\tau~{}\\!V^{I}_{\rm
int}(\tau)\right\\},$

in which ${\rm T}_{\tau}$ is the symbol of the “chronological” ordering. Since
$e^{-\beta\hat{K}}=e^{-\beta\hat{K}_{0}}~{}\\!{\cal U}_{I}(\beta,0)$, the
statistical sum can be represented as

$\displaystyle\Xi_{\rm stat}={\rm Tr}\\!\left(e^{-\beta\hat{K}_{0}}~{}\\!{\cal
U}_{I}(\beta,0)\right)\equiv\Xi_{\rm stat}^{(0)}{\rm
Tr}\\!\left(\hat{\rho}^{(0)}~{}\\!{\cal U}_{I}(\beta,0)\right),$ (8)

where $\hat{\rho}^{(0)}$ and $\Xi_{\rm stat}^{(0)}$ are the statistical
operator and the statistical sum of the noninteracting system, respectively.
The perturbative expansion of $\Xi_{\rm stat}$ is then given by the series

$\displaystyle\Xi_{\rm stat}=\Xi_{\rm
stat}^{(0)}\sum_{n=0}^{\infty}{(-1)^{n}\over
n!}\\!\int_{0}^{\beta}\\!d\tau_{n}\dots\int_{0}^{\beta}\\!d\tau_{1}{\rm
Tr}\left(\hat{\rho}^{(0)}~{}\\!{\rm T}_{\tau}[V^{I}_{\rm int}(\tau_{n})\dots
V^{I}_{\rm int}(\tau_{1})]\right).$ (9)

The corresponding expansion of the potential
$\Omega(T,V,\mu_{+},\mu_{-})=-{1\over\beta}\ln\Xi_{\rm
stat}(T,V,\mu_{+},\mu_{-})$ is

$\displaystyle\Omega=\Omega^{(0)}-{1\over\beta}\ln\\!\left\\{\sum_{n=0}^{\infty}{(-1)^{n}\over
n!}\\!\int_{0}^{\beta}\\!d\tau_{n}\dots\int_{0}^{\beta}\\!d\tau_{1}{\rm
Tr}\\!\left(\hat{\rho}^{(0)}~{}\\!{\rm T}_{\tau}[V^{I}_{\rm
int}(\tau_{n})\dots V^{I}_{\rm int}(\tau_{1})]\right)\right\\},$ (10)

its first term $\Omega^{(0)}$ being the textbook expression [1]
($\varepsilon_{\mathbf{k}}=\hbar^{2}\mathbf{k}^{2}/2m_{f}$)

$\displaystyle\Omega^{(0)}(T,V,\mu_{\sigma})=-{1\over\beta}~{}\\!\sum_{\sigma}V\\!\int\\!{d^{3}\mathbf{k}\over(2\pi)^{3}}\ln\\!\left(1+e^{-\beta(\varepsilon_{\mathbf{k}}-\mu_{\sigma})}\right).$
(11)

Or, since the logarithm picks up connected contributions only,

$\displaystyle\Omega=\Omega^{(0)}-{1\over\beta}\sum_{n=0}^{\infty}{(-1)^{n}\over
n!}\\!\int_{0}^{\beta}\\!d\tau_{n}\dots\int_{0}^{\beta}\\!d\tau_{1}{\rm
Tr}\\!\left(\hat{\rho}^{(0)}~{}\\!{\rm T}_{\tau}[V^{I}_{\rm
int}(\tau_{n})\dots V^{I}_{\rm int}(\tau_{1})]\right)^{\rm con}.$ (12)

In this form the expression for $\Omega$ is just the thermal analog of the
expansion of the formula222Here $T$ denotes time and not the temperature.

$\displaystyle E_{\Omega}=E_{\Omega_{0}}-\lim_{T\rightarrow\infty}{\hbar\over
iT}~{}\\!\langle\Omega_{0}|{\rm
T}_{t}\exp\\!\left(-{i\over\hbar}\\!\int_{-T/2}^{T/2}\\!dt~{}\\!V^{I}_{\rm
int}(t)\right)|\Omega_{0}\rangle^{\rm con}~{}\\!,$ (13)

used in [6, 8, 11, 12] for computing the ground state energy $E_{\Omega}$ of
the system. It is clear that the correspondence between the two formalisms is
$\beta\leftrightarrow iT/\hbar$ (it transforms the $K$-picture operators into
the Heisenberg picture ones and vice versa). The formula (13) for the ground
state energy is thus obtained from the thermal expansion (12) by taking the
limit $\beta\rightarrow\infty$ and simultaneously adjusting the chemical
potential $\mu_{\sigma}$ so that there are $N_{\sigma}$ particles with the
spin projection $\sigma$ (see below).

Evaluation of the successive terms of the expansion (12) reduces, owing to the
thermal analog of the Wick formula (see [7]), to drawing all possible
connected Feynman diagrams with a given numbers of different interaction
vertices arising from $\hat{V}_{\rm int}$ joined by the oriented lines and
integrating over the positions $\mathbf{x}$ and “times” $\tau$ ascribed to
these vertices the corresponding products of free thermal propagators

$\displaystyle-{\cal
G}^{(0)}_{\sigma_{2},\sigma_{1}}(\tau_{2}-\tau_{1};\mathbf{x}_{2}-\mathbf{x}_{1})={1\over\beta}\sum_{n}\int\\!{d^{3}\mathbf{k}\over(2\pi)^{3}}~{}\\!e^{-i\omega^{F}_{n}(\tau_{2}-\tau_{1})}~{}\\!e^{i\mathbf{k}\cdot(\mathbf{x}_{2}-\mathbf{x}_{1})}\left(-\tilde{\cal
G}^{(0)}_{\sigma_{2},\sigma_{1}}(\omega^{F}_{n},\mathbf{k})\right),$

the Fourier transforms $-\tilde{\cal G}^{(0)}_{\sigma_{2},\sigma_{1}}$ of
which have the form [7] (the definition of $\omega^{F}_{n}$ as well as of
$\omega^{B}_{n}$ are given in (A.1) and (A.2))

$\displaystyle-\tilde{\cal
G}^{(0)}_{\sigma_{2},\sigma_{1}}(\omega^{F}_{n},\mathbf{k})={-\delta_{\sigma_{2},\sigma_{1}}\over
i\omega^{F}_{n}-(\varepsilon_{\mathbf{k}}-\mu_{\sigma_{1}})}~{}\\!,$

associated with the (oriented) lines connecting the vertices of the diagram.
The resulting Feynman rules in the “momentum” space are almost identical with
the ordinary ones except that the integrations over frequencies are replaced
by summations over the (fermionic) Matsubara frequencies
$\omega^{F}_{n}=(\pi/\beta)(2n+1)$, $n\in\mathbb{Z}$. In this way one obtains
the expansion of the potential $\Omega(T,V,\mu_{+},\mu_{-})$ the successive
terms of which depend on the chemical potentials $\mu_{+}$ and $\mu_{-}$ which
must be adjusted in successive orders of the expansion to yield through the
relations

$\displaystyle
N_{\pm}=-\left(\partial\Omega/\partial\mu_{\pm}\right)_{T,V}~{}\\!,$ (14)

the prescribed densities $n_{+}=N_{+}/V$ and $n_{-}=N_{-}/V$ of particles with
the spin projections up and down.

Figure 1: The first order correction $\Omega^{(1)}$ to the thermodynamic
potential $\Omega(T,V,\mu_{+},\mu_{-})$. Solid and dashed lines represent
fermions with opposite spin projections.

It will be instructive to recover first, using this formalism, the textbook
results [1, 2] of the mean field approximation. The first correction
$\Omega^{(1)}$ to the Grand potential is given by the single diagram shown in
Figure 1. The corresponding expression reads

$\displaystyle\Omega^{(1)}={1\over\beta}~{}\\!C_{0}\\!\int_{0}^{\beta}\\!d\tau\\!\int\\!d^{3}\mathbf{x}~{}\\!{\rm
Tr}\\!\left(\hat{\rho}^{(0)}{\rm T}_{\tau}[\hat{\psi}^{\dagger
I}_{+}\hat{\psi}^{I}_{+}\hat{\psi}^{\dagger
I}_{-}\hat{\psi}^{I}_{-}]\right)=C_{0}~{}\\!V~{}\\!{\cal
G}^{(0)}_{++}(0,\mathbf{0})~{}\\!{\cal G}^{(0)}_{--}(0,\mathbf{0})~{}\\!.$
(15)

Using the summation formula (A.1) one obtains

$\displaystyle{\cal
G}^{(0)}_{\pm\pm}(0,\mathbf{0})=\int\\!{d^{3}\mathbf{k}\over(2\pi)^{3}}\left[1+e^{\beta(\varepsilon_{\mathbf{k}}-\mu_{\pm})}\right]^{-1}.\phantom{aaaaaaaa}$
(16)

As will be shown in the next section, to the first order in the coupling
$C_{0}$ the free energy $F(T,V,N_{+},N_{-})$ is given by

$\displaystyle
F(T,V,N_{+},N_{-})=\Omega^{(0)}(T,V,\mu^{(0)}_{+},\mu^{(0)}_{-})+N_{+}\mu^{(0)}_{+}+N_{-}\mu^{(0)}_{-}$
$\displaystyle+~{}\\!\Omega^{(1)}(T,V,\mu^{(0)}_{+},\mu^{(0)}_{-})+\dots,\phantom{aaaaaaaaaaa}\\!$
(17)

where $\mu^{(0)}_{\pm}$ are the zeroth order chemical potentials determined by
the conditions analogous to (14) but with $\Omega$ replaced by $\Omega^{(0)}$
given by (11). It is convenient to define the function

$\displaystyle f(\nu)\equiv{3\over
2}\\!\int_{0}^{\infty}\\!d\xi~{}\\!{\xi^{1/2}\over
1+e^{\xi-\nu}}\equiv{3\sqrt{\pi}\over 4}~{}\\!f_{3/2}(\nu)~{}\\!,$ (18)

and to rewrite these conditions in the form

$\displaystyle f(\nu_{\pm})=\left({\varepsilon^{(0)}_{\rm F}(n_{\pm})\over
k_{\rm B}T}\right)^{3/2}~{}\\!,$ (19)

in which $\nu_{\pm}\equiv\mu_{\pm}/k_{\rm B}T$ and $\varepsilon^{(0)}_{\rm
F}(n)=(6\pi^{2}n)^{2/3}\hbar^{2}/2m_{f}$ is the Fermi energy of the system of
$N=nV$ spin $0$ noninteracting fermions enclosed in the volume $V$. The
function $f(\nu)$, which is a decent monotonically growing function of $\nu$
mapping $\mathbb{R}$ onto $\mathbb{R}_{+}$, has the inverse, so after writing
$n_{\pm}$ as $(n/2)(1\pm P)$ the solutions take the form333Inverting the
appropriate expansions of the integral (18) given e.g. in [1] it is
straightforward to find that asymptotically  $\displaystyle
f^{-1}(x)=\left\\{\matrix{\ln\\!\left(\sqrt{2}-\sqrt{2-(4x/3)\sqrt{8/\pi}}\right)+\dots,\phantom{aa}x\ll
1\cr
x^{2/3}\left[1-(\pi^{2}/12)~{}\\!x^{-4/3}-(\pi^{4}/80)~{}\\!x^{-8/3}-(1511\pi^{6}/207360)~{}\\!x^{-4}+\dots\right],\phantom{aa}x\gg
1}\right.$

$\displaystyle{\mu^{(0)}_{\pm}\over k_{\rm B}T}=f^{-1}\\!\left((1\pm
P)\left({\varepsilon_{\rm F}(n)\over k_{\rm B}T}\right)^{3/2}\right),$ (20)

in which $\varepsilon_{\rm F}(n)$ the system’s overall Fermi energy (1).

Expressed in terms of the zeroth order chemical potentials $\mu_{\pm}^{(0)}$,
the first order correction (15) can be simply written as
$\Omega^{(1)}(T,V,\mu^{(0)}_{+},\mu^{(0)}_{-})=C_{0}V(N_{+}/V)(N_{-}/V)$, i.e.
it is independent (when expressed in terms of the particle densities) of the
temperature. Minimization with respect to $N_{+}$ of $F(T,V,N_{+},N-N_{+})$
truncated to the first order in the coupling $C_{0}$ (at fixed $N$) then leads
to the equilibrium condition

$\displaystyle\mu^{(0)}_{+}(N_{+})-\mu^{(0)}_{-}(N-N_{+})+{C_{0}\over
V}~{}\\!(N-2N_{+})=0~{}\\!,$

which, because $N-2N_{+}=-NP$ and (to this order)
$C_{0}=(4\pi\hbar^{2}/m_{f})a_{0}$, can be rewritten in the familiar form [1]

$\displaystyle\mu^{(0)}_{+}(N_{+})-\mu^{(0)}_{-}(N_{-})={8\over
3\pi}~{}\\!\varepsilon_{\rm F}~{}\\!(k_{\rm F}a_{0})~{}\\!P~{}\\!.$ (21)

This leads to the continuous phase transition.

The effect of the external magnetic field ${\cal H}$ can be also taken into
account by simply including the interaction with it in the free part of the
Hamiltonian, i.e. by replacing $\mu_{\pm}$ in (6) by
$\tilde{\mu}_{\pm}=\mu_{\pm}\pm{\cal H}$ (the magnetic moment has been here
included in ${\cal H}$ which has therefore the dimension of energy). Since
ultimately the free energy will be cast in the form in which its dependence on
$N_{\pm}$ and ${\cal H}$ enters only through $\tilde{\mu}^{(0)}_{\pm}$ which
should be determined from the conditions

$\displaystyle{\tilde{\mu}^{(0)}_{\pm}\over k_{\rm
B}T}={\mu^{(0)}_{\pm}\pm{\cal H}\over k_{\rm B}T}=f^{-1}\\!\left((1\pm
P)\left({\varepsilon_{\rm F}(n)\over k_{\rm B}T}\right)^{3/2}\right),$ (22)

this prescription remains valid to all orders of the expansion. In particular,
in the first order approximation the equilibrium condition, written in the
convenient dimensionless variables

$\displaystyle t\equiv{T\over T_{\rm F}}\equiv{k_{\rm B}T\over\varepsilon_{\rm
F}}~{}\\!,\phantom{aaa}h\equiv{{\cal H}\over\varepsilon_{\rm
F}}~{}\\!,\phantom{aaa}\delta_{\pm}\equiv{\mu^{(0)}_{\pm}\over\varepsilon_{\rm
F}}~{}\\!,$ (23)

takes the form

$\displaystyle{8\over 3\pi}~{}\\!(k_{\rm
F}a_{0})~{}\\!P+2h=t\left[f^{-1}\\!\left({1+P\over
t^{3/2}}\right)-f^{-1}\\!\left({1-P\over t^{3/2}}\right)\right].$ (24)

If the asymptotic expansion of $f^{-1}(x)$ for $x\gg 1$ is used, this
reproduces the equilibrium condition derived in [1].

For further applications it will be convenient to write down explicitly the
formula (17) (including the external magnetic field ${\cal H}$) expressing it
through the introduced dimensionless variables (23) and the polarization (2):

$\displaystyle{6\pi^{2}\over k_{\rm F}^{3}}~{}\\!{F\over\varepsilon_{\rm
F}V}=-{3\sqrt{\pi}\over
4}~{}\\!t^{5/2}\left[f_{5/2}(\tilde{\nu}_{+})+f_{5/2}(\tilde{\nu}_{-})\right]\phantom{aaaaaaaaaaaaaaaaaaaaaaaaa}$
$\displaystyle+~{}\\!(1+P)(\tilde{\delta}_{+}-h)+(1-P)(\tilde{\delta}_{-}+h)+(k_{\rm
F}a_{0})~{}\\!{4\over 3\pi}~{}\\!(1-P^{2})+\dots$ (25)

Here444By the appropriate change of variables and the integration by parts
$\Omega^{(0)}$ given by (11) is written in terms of the standard integral (26)
with $p=5/2$ [1].

$\displaystyle
f_{p}(\nu)={1\over\Gamma(p)}\int_{0}^{\infty}\\!{d\xi~{}\\!\xi^{p-1}\over
1+e^{\xi-\nu}}~{}\\!,$ (26)

and $\tilde{\nu}_{\pm}$ (and $\tilde{\delta}_{\pm}\equiv t\tilde{\nu}_{\pm}$)
are given by (22). In the limit $T\rightarrow 0$ ($t\rightarrow 0$) in which
$\tilde{\nu}_{\pm}\gg 1$, $f_{5/2}(\nu)=(4/3\sqrt{\pi})(2/5)\nu^{5/2}+\dots,$
while (c.f. the expansion of the function $f^{-1}(x)$ given in the footnote)
$\tilde{\nu}_{\pm}=(1\pm P)^{2/3}+\dots$ and the right hand side of (25) tends
to

$\displaystyle-{2\over
5}\left[(1+P)^{5/3}+(1-P)^{5/3}\right]+(1+P)\left[(1+P)^{2/3}-h\right]\phantom{aaaaaaaaaaaaaaaaaa}~{}\\!$
$\displaystyle+(1-P)\left[(1-P)^{2/3}+h\right]+(k_{\rm F}a_{0})~{}\\!{4\over
3\pi}~{}\\!(1-P^{2})~{}\\!,$

reproducing, of course, the well known formula for the ground-state energy
given (for ${\cal H}=0$) e.g. in [12].

##  3 Expansion of the free energy

From the thermodynamic point of view much more convenient to work with than
the potential $\Omega$ is the free energy $F=\Omega+\mu_{+}N_{+}+\mu_{-}N_{-}$
which canonically depends on $T$, $V$ and the particle numbers $N_{\pm}$. It
turns out that the expansion of this potential is also simpler. We will derive
it here up to the third order following the method outlined in [15]. To make
the notation more transparent we will denote the chemical potentials as

$\displaystyle\mu_{+}\equiv
x=x_{0}+x_{1}+x_{2}+\dots,\phantom{aaa}\mu_{-}\equiv
y=y_{0}+y_{1}+y_{2}+\dots,$ (27)

where the successive terms $x_{n}$, $y_{n}$ correspond to the successive terms
$\Omega^{(n)}$ of the expansion of the potential $\Omega$. Introducing the
notation $\Omega^{(n)}_{x}$, $\Omega^{(n)}_{y}$, $\Omega^{(n)}_{xx}$, etc. for
the first, second, etc. derivatives of $\Omega^{(n)}$ with respect to their
chemical potential arguments and expanding the right hand side of the relation
($N_{x}\equiv N_{+}$, $N_{y}\equiv N_{-}$)

$\displaystyle
F=\Omega^{(0)}(x_{0}+x_{1}+x_{2}+x_{3}+\dots,~{}\\!y_{0}+y_{1}+y_{2}+y_{3}+\dots)\phantom{aaaaaaaaa}~{}\\!$
$\displaystyle+~{}\\!\Omega^{(1)}(x_{0}+x_{1}+x_{2}+\dots,~{}\\!y_{0}+y_{1}+y_{2}+\dots)\phantom{aaaaaaaaaaaaaaaaa}~{}$
$\displaystyle+~{}\\!\Omega^{(2)}(x_{0}+x_{1}+\dots,~{}\\!y_{0}+y_{1}+\dots)+\Omega^{(3)}(x_{0}+\dots,y_{0}+\dots)+\dots$
$\displaystyle+~{}\\!(x_{0}+x_{1}+x_{2}+x_{3}+\dots)N_{x}+(y_{0}+y_{1}+y_{2}+y_{3}+\dots)N_{y}~{}\\!,\phantom{aaaa}~{}\\!$

one obtains, using the zeroth order relations $\Omega^{(0)}_{x}=-N_{x}$ and
$\Omega^{(0)}_{y}=-N_{y}$ and the fact that
$\Omega^{(0)}(x_{0},y_{0})=\Omega^{\rm free}(x_{0})+\Omega^{\rm free}(y_{0})$
(cf. the formula (11)), i.e. that $\Omega^{(0)}_{xy}=0$,

$\displaystyle
F=\left(\Omega^{(0)}+x_{0}N_{x}+y_{0}N_{y}\right)+\left(\Omega^{(1)}\right)\phantom{aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa}$
$\displaystyle+\left(\Omega^{(2)}+x_{1}~{}\\!\Omega^{(1)}_{x}+y_{1}~{}\\!\Omega^{(1)}_{y}+{1\over
2}~{}\\!x_{1}^{2}~{}\\!\Omega^{(0)}_{xx}+{1\over
2}~{}\\!y_{1}^{2}~{}\\!\Omega^{(0)}_{yy}\right)\phantom{aaaaaaaaaaaaaaaaaaaaaa}$
$\displaystyle+\left(\Omega^{(3)}+x_{1}~{}\\!\Omega^{(2)}_{x}+y_{1}~{}\\!\Omega^{(2)}_{y}+{1\over
2}~{}\\!x_{1}^{2}~{}\\!\Omega^{(1)}_{xx}+{1\over
2}~{}\\!y_{1}^{2}~{}\\!\Omega^{(1)}_{yy}+x_{1}~{}\\!y_{1}~{}\\!\Omega^{(1)}_{xy}\right.\phantom{aaaaaaaaaaaaa}$
(28)
$\displaystyle\left.+~{}\\!x_{2}~{}\\!\Omega^{(1)}_{x}+y_{2}~{}\\!\Omega^{(1)}_{y}+x_{1}~{}\\!x_{2}~{}\\!\Omega^{(0)}_{xx}+y_{1}~{}\\!y_{2}~{}\\!\Omega^{(0)}_{yy}+{1\over
6}~{}\\!x_{1}^{3}~{}\\!\Omega^{(0)}_{xxx}+{1\over
6}~{}\\!y_{1}^{3}~{}\\!\Omega^{(0)}_{yyy}\right)+\dots,$

all functions being now evaluated at $x_{0}$ and $y_{0}$ (at
$\tilde{x}_{0}=\mu_{+}^{(0)}+{\cal H}$ and $\tilde{y}_{0}=\mu_{-}^{(0)}-{\cal
H}$ if there is an external magnetic field). The terms in the successive
brackets are the successive terms of the expansion of the free energy. The
first order correction $F^{(1)}$ used in the preceding section is indeed given
by $\Omega^{(1)}(x_{0},y_{0})$ (by
$\Omega^{(1)}(\tilde{x}_{0},\tilde{y}_{0})$). Furthermore, expanding around
$x_{0}$ and $y_{0}$ (or $\tilde{x}_{0}$ and $\tilde{y}_{0}$) the right hand
side of the relation (14) which determines the chemical potential $x$

$\displaystyle-
N_{x}=\Omega^{(0)}_{x}+(x_{1}+x_{2})~{}\\!\Omega^{(0)}_{xx}+{1\over
2}~{}\\!x_{1}^{2}~{}\\!\Omega^{(0)}_{xxx}+\Omega^{(1)}_{x}+x_{1}~{}\\!\Omega^{(1)}_{xx}+y_{1}~{}\\!\Omega^{(1)}_{xy}+\Omega^{(2)}_{x}+\dots,$

and the other similar relation for $y$, and taking into account that $x_{0}$
and $y_{0}$ are such that $-N_{x}=\Omega^{(0)}_{x}$,
$-N_{y}=\Omega^{(0)}_{y}$, one obtains

$\displaystyle
x_{1}=-{\Omega^{(1)}_{x}\over\Omega^{(0)}_{xx}}~{}\\!,\phantom{aaaaaaaaaaaaaaaaaaaaaaaaaaaaaa}~{}\\!$
$\displaystyle
x_{2}=-{\Omega^{(2)}_{x}\over\Omega^{(0)}_{xx}}+{\Omega^{(1)}_{xx}\Omega^{(1)}_{x}\over[\Omega^{(0)}_{xx}]^{2}}+{\Omega^{(1)}_{xy}\Omega^{(1)}_{y}\over\Omega^{(0)}_{xx}\Omega^{(0)}_{yy}}-{\Omega^{(0)}_{xxx}[\Omega^{(1)}_{x}]^{2}\over
2[\Omega^{(0)}_{xx}]^{3}}~{}\\!.$ (29)

$y_{1}$ and $y_{2}$ are given by the analogous formulae. Inserting the
corrections to the chemical potentials determined in this way into the
formulae for $F^{(2)}$ and $F^{(3)}$ one finds that (again all functions are
evaluated at $x_{0}$ and $y_{0}$ or at $\tilde{x}_{0}$ and $\tilde{y}_{0}$)

$\displaystyle F^{(2)}=\Omega^{(2)}-{[\Omega^{(1)}_{x}]^{2}\over
2~{}\\!\Omega^{(0)}_{xx}}-{[\Omega^{(1)}_{y}]^{2}\over
2~{}\\!\Omega^{(0)}_{yy}}~{}\\!.$ (30)

and (the formulae for $x_{1}$ and $y_{1}$ immediately imply that the first
four terms in the last line of the formula (28) sum up to zero) that

$\displaystyle
F^{(3)}=\Omega^{(3)}-{\Omega^{(2)}_{x}\Omega^{(1)}_{x}\over\Omega^{(0)}_{xx}}-{\Omega^{(2)}_{y}\Omega^{(1)}_{y}\over\Omega^{(0)}_{yy}}\phantom{aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa}~{}$
(31) $\displaystyle+~{}\\!{\Omega^{(1)}_{xx}[\Omega^{(1)}_{x}]^{2}\over
2~{}\\![\Omega^{(0)}_{xx}]^{2}}+{\Omega^{(1)}_{yy}[\Omega^{(1)}_{y}]^{2}\over
2~{}\\![\Omega^{(0)}_{yy}]^{2}}+{\Omega^{(1)}_{xy}\Omega^{(1)}_{x}\Omega^{(1)}_{y}\over\Omega^{(0)}_{xx}\Omega^{(0)}_{yy}}-{\Omega^{(0)}_{xxx}[\Omega^{(1)}_{x}]^{3}\over
6~{}\\![\Omega^{(0)}_{xx}]^{3}}-{\Omega^{(0)}_{yyy}[\Omega^{(1)}_{y}]^{3}\over
6~{}\\![\Omega^{(0)}_{yy}]^{3}}~{}\\!.$

It will be seen that the extra terms in (30) precisely cancel the
contributions to $\Omega^{(2)}$ of those diagrams which do not contribute to
the expansion of the formula (13) for the ground state energy density. The
analogous cancellation of the extra terms in (31) and the in $\Omega^{(3)}$ is
demonstrated in Appendix B.

$b)$$c)$ Figure 2: The order $C^{2}_{0}$ contributions $\Omega^{(2)b}$ and
$\Omega^{(2)c}$.

##  4 Computation of $F^{(2)}$

Diagrams contributing to $\Omega^{(2)}$ are shown in Figures 2 and 3 (the left
one). It is straightforward to check that the contributions $\Omega^{(2)b}$
and $\Omega^{(2)c}$ of the ones of Figure 2 cancel against the last two terms
in the formula (30). Indeed, with the help of the summation rules collected in
Appendix A and taking into account that these contributions are evaluated at
$x_{0}$ and $y_{0}$ one easily obtains ($\Omega^{(2)c}$ is given by an
analogous formula)

$\displaystyle\Omega^{(2)b}={C_{0}^{2}V\over 2}\left({N_{-}\over
V}\right)^{2}\int\\!{d^{3}\mathbf{p}\over(2\pi)^{3}}\left[{d\over
da}~{}\\!{1\over 1+e^{\beta
a}}\right]_{a=\varepsilon_{\mathbf{p}}-x_{0}}=-{1\over
2}~{}\\!C_{0}^{2}V\beta~{}\\!(n_{x}-n_{xx})n_{y}^{2}~{}\\!,$ (32)

where the second form of $\Omega^{(2)b}$ is given in the notation introduced
in Appendix B. With the help of the formulae (B.1), (B.2) it is immediately
seen that it is canceled by the second term of (30). Thus

$\displaystyle\Omega^{(2)b}+\Omega^{(2)c}-{[\Omega^{(1)}_{x}]^{2}\over
2~{}\\!\Omega^{(0)}_{xx}}-{[\Omega^{(1)}_{y}]^{2}\over
2~{}\\!\Omega^{(0)}_{yy}}=0~{}\\!.$

Hence, $F^{(2)}=\Omega^{(2)a}$ evaluated at $x_{0}$ and $y_{0}$ (or at
$\tilde{x}_{0}$ and $\tilde{y}_{0}$).

The integrals and sums corresponding to the left diagram of Figure 3 giving
$\Omega^{(2)a}$ can be written in three different forms (corresponding to
three different routings of the internal momenta and frequencies) of which two
can be composed of two “elementary” blocks $A$ and $B$ shown in Figure 3
right:

$\displaystyle\Omega^{(2)a}=-{1\over
2}~{}\\!C_{0}^{2}V~{}\\!{1\over\beta}\sum_{l\in\mathbb{Z}}\\!\int\\!{d^{3}\mathbf{q}\over(2\pi)^{3}}~{}\\![A(\omega_{l}^{B},\mathbf{q})]^{2}=-{1\over
2}~{}\\!C_{0}^{2}V~{}\\!{1\over\beta}\sum_{l\in\mathbb{Z}}\\!\int\\!{d^{3}\mathbf{q}\over(2\pi)^{3}}~{}\\![B(\omega_{l}^{B},\mathbf{q})]^{2}.$

where (here
$n_{\pm}(\mathbf{p})\equiv[1+\exp\\{\beta(\varepsilon_{\mathbf{p}}-\mu^{(0)}_{\pm})\\}]^{-1}$,
$\mu^{(0)}_{+}\equiv x_{0}$, $\mu^{(0)}_{-}\equiv y_{0}$)

$\displaystyle
A(\omega_{l+1}^{B},\mathbf{q})={1\over\beta}\sum_{n\in\mathbb{Z}}\\!\int\\!{d^{3}\mathbf{k}\over(2\pi)^{3}}~{}\\!{1\over
i\omega_{n}^{F}-(\varepsilon_{\mathbf{k}}-x_{0})}~{}\\!{1\over
i\omega_{l-n}^{F}-(\varepsilon_{\mathbf{q}-\mathbf{k}}-y_{0})}$
$\displaystyle=\int\\!{d^{3}\mathbf{k}\over(2\pi)^{3}}~{}\\!{n_{+}(\mathbf{k})+n_{-}(\mathbf{q}-\mathbf{k})-1\over
i\omega_{l+1}^{B}-(\varepsilon_{\mathbf{k}}-x_{0}+\varepsilon_{\mathbf{q}-\mathbf{k}}-y_{0})}~{}\\!.\phantom{aaaaaaaaa}$
(33)

and

$\displaystyle
B(\omega_{l}^{B},\mathbf{q})={1\over\beta}\sum_{n\in\mathbb{Z}}\\!\int\\!{d^{3}\mathbf{k}\over(2\pi)^{3}}~{}\\!{1\over
i\omega_{n}^{F}-(\varepsilon_{\mathbf{k}+\mathbf{q}}-x_{0})}~{}\\!{1\over
i\omega_{n+l}^{F}-(\varepsilon_{\mathbf{k}}-y_{0})}$
$\displaystyle=\int\\!{d^{3}\mathbf{k}\over(2\pi)^{3}}~{}\\!{n_{+}(\mathbf{k}+\mathbf{q})-n_{-}(\mathbf{k})\over
i\omega_{l}^{B}-(\varepsilon_{\mathbf{k}}-y_{0}-\varepsilon_{\mathbf{k}+\mathbf{q}}+x_{0})}~{}\\!.\phantom{aaaaaaaaaa}~{}$
(34)

(The contributions $\Omega^{(3)a}$ and $\Omega^{(3)b}$ of the left and right
diagrams shown in Figure 7 can be written analogously with
$[A(\omega_{l}^{B},\mathbf{q})]^{3}$ and $[B(\omega_{l}^{B},\mathbf{q})]^{3}$,
respectively [11, 12]).

$A(\omega_{l+1}^{B},\mathbf{q})=$${}^{q-k,~{}l-n}$${}_{k,~{}n}$$B(\omega_{l}^{B},\mathbf{q})=$${}^{k,~{}n+l}$${}_{k+q,~{}n}$
Figure 3: The order $C_{0}^{2}$ diagram contributing to the correction
$\Omega^{(2)}$ and two “elementary” one-loop diagrams out of which the second
order and the third order corrections with the $C_{0}$ couplings can be
constructed. Solid and dashed lines denote propagators of fermions with the
spin projections $+$ and $-$, respectively.

With the help of the sum rule (A.5) the sum over $l$ of two $A$-blocks can be
done and gives (the symbol $\int_{\mathbf{k}}$ stands for the integral over
the measure $d^{3}\mathbf{k}/(2\pi)^{3}$)

$\displaystyle\int_{\mathbf{k}}\\!\int_{\mathbf{p}}{[n_{+}(\mathbf{k})+n_{-}(\mathbf{q}-\mathbf{k})-1][n_{+}(\mathbf{p})+n_{-}(\mathbf{q}-\mathbf{p})-1]\over\varepsilon_{\mathbf{k}}+\varepsilon_{\mathbf{q}-\mathbf{k}}-\varepsilon_{\mathbf{p}}-\varepsilon_{\mathbf{q}-\mathbf{p}}}\phantom{aaaaaaaaaaaaaaaaaa}$
$\displaystyle\times\left({1\over
1-e^{\beta(\varepsilon_{\mathbf{k}}-x_{0})}e^{\beta(\varepsilon_{\mathbf{q}-\mathbf{k}}-y_{0})}}-{1\over
1-e^{\beta(\varepsilon_{\mathbf{p}}-x_{0})}e^{\beta(\varepsilon_{\mathbf{q}-\mathbf{p}}-y_{0})}}\right).$

The identity

$\displaystyle
n_{+}(\mathbf{k})+n_{-}(\mathbf{q}-\mathbf{k})-1=n_{+}(\mathbf{k})~{}\\!n_{-}(\mathbf{q}-\mathbf{k})\left[1-e^{\beta(\varepsilon_{\mathbf{k}}-x_{0})}e^{\beta(\varepsilon_{\mathbf{q}-\mathbf{k}}-y_{0})}\right].$
(35)

and the fact that the two terms in the bracket above gives equal contributions
allows then to write

$\displaystyle\Omega^{(2)a}=C_{0}^{2}V\\!\int_{\mathbf{q}}\\!\int_{\mathbf{p}}\\!\int_{\mathbf{k}}\\!{n_{+}(\mathbf{k})~{}\\!n_{-}(\mathbf{q}-\mathbf{k})[1-n_{+}(\mathbf{p})-n_{-}(\mathbf{q}-\mathbf{p})]\over\varepsilon_{\mathbf{k}}+\varepsilon_{\mathbf{q}-\mathbf{k}}-\varepsilon_{\mathbf{p}}-\varepsilon_{\mathbf{q}-\mathbf{p}}}~{}\\!.$
(36)

It is interesting to notice that because the integral of the quartic product
$n_{+}(\mathbf{k})~{}\\!n_{-}(\mathbf{q}-\mathbf{k})n_{+}(\mathbf{p})~{}\\!n_{-}(\mathbf{q}-\mathbf{p})$
vanishes (the numerator is even with respect to the interchange
$\mathbf{k}\leftrightarrow\mathbf{p}$ while the denominator is odd), the
expression for $\Omega^{(2)a}$ can be written (after the change
$\mathbf{p}=-\mathbf{u}+\mathbf{s}$, $\mathbf{k}=-\mathbf{t}+\mathbf{s}$,
$\mathbf{q}=2\mathbf{s}$ of the integration variables) in the form completely
analogous to the expression giving $E_{\Omega}/V$ (see [8]), the only
modification being the change in the prefactor and the replacement of
$\theta(k-|\mathbf{v}|)$ and $\theta(|\mathbf{v}|-k)$ by $n(|\mathbf{v}|)$ and
$1-n(|\mathbf{v}|)$, respectively. (Curiously enough, we have found that this
simple analogy does not work for the diagrams of Figure 7).

It is straightforward to see that the expression (36) is divergent, the
divergence arising from the unity in the square bracket in the numerator. In
the variables $\mathbf{s}$, $\mathbf{t}$ and $\mathbf{u}$ the integral over
$\mathbf{u}$ is the one evaluated with the cutoff $\Lambda$ in [8] and using
this result and changing once more the variables to
$\mathbf{k}=\mathbf{t}-\mathbf{s}$, $\mathbf{p}=\mathbf{t}+\mathbf{s}$, after
adding the contribution $\Omega^{(1)}$ and expressing $C_{0}$ in terms of the
scattering lengths $a_{0}$

$\displaystyle C_{0}(\Lambda)={4\pi\hbar^{2}\over
m_{f}}~{}\\!a_{0}\left(1+{2\over\pi}~{}\\!a_{0}\Lambda+\dots\right),$ (37)

[16, 11, 12], one arrives at the finite (to the second order) result

$\displaystyle\Omega^{(1)}+\Omega^{(2)a}={4\pi\hbar^{2}\over
m_{f}}~{}\\!a_{0}\left(1+{2\over\pi}~{}\\!\Lambda
a_{0}+\dots\right)V~{}\\!{N_{-}\over V}~{}\\!{N_{+}\over
V}\phantom{aaaaaaaaaaaaaaaaaaaaaaaaaaa}$ $\displaystyle-~{}\\!{\Lambda\over
2\pi^{2}}\left({4\pi\hbar^{2}\over
m_{f}}~{}\\!a_{0}\right)^{2}{m_{f}\over\hbar^{2}}~{}\\!V~{}\\!{N_{-}\over
V}~{}\\!{N_{+}\over V}+\Omega^{(2)a}_{\rm finite}={4\pi\hbar^{2}\over
m_{f}}~{}\\!a_{0}~{}\\!V~{}\\!{N_{-}\over V}~{}\\!{N_{+}\over
V}+\Omega^{(2)a}_{\rm finite}~{}\\!.$

The finite part of $\Omega^{(2)a}$,

$\displaystyle\Omega^{(2)a}_{\rm
fin}=-C_{0}^{2}V\\!\int_{\mathbf{q}}\int_{\mathbf{k}}n_{+}(\mathbf{k})~{}\\!n_{-}(\mathbf{q}-\mathbf{k})\int_{\mathbf{p}}{n_{+}(\mathbf{p})+n_{-}(\mathbf{q}-\mathbf{p})\over\varepsilon_{\mathbf{k}}+\varepsilon_{\mathbf{q}-\mathbf{k}}-\varepsilon_{\mathbf{p}}-\varepsilon_{\mathbf{q}-\mathbf{p}}}~{}\\!,$

upon setting first $\mathbf{k}=\mathbf{k}_{1}$,
$\mathbf{q}-\mathbf{k}=\mathbf{k}_{2}$ and then replacing in the term with
$n_{-}(\mathbf{k}_{1}+\mathbf{k}_{2}-\mathbf{p})$ the variable
$\mathbf{k}_{1}+\mathbf{k}_{2}-\mathbf{p}$ by $\mathbf{p}^{\prime}$ (upon
which
$\varepsilon_{\mathbf{k}_{1}+\mathbf{k}_{2}-\mathbf{p}}\rightarrow\varepsilon_{\mathbf{p}^{\prime}}$
but at the same time
$\varepsilon_{\mathbf{p}}\rightarrow\varepsilon_{\mathbf{k}_{1}+\mathbf{k}_{2}-\mathbf{p}^{\prime}}$)
can be cast in the convenient symmetric form

$\displaystyle\Omega^{(2)a}_{\rm
fin}=F^{(2)}=-C_{0}^{2}V\\!\int_{\mathbf{k}_{1}}\int_{\mathbf{k}_{2}}n_{+}(\mathbf{k}_{1})~{}\\!n_{-}(\mathbf{k}_{2})\int_{\mathbf{p}}{n_{+}(\mathbf{p})+n_{-}(\mathbf{p})\over\varepsilon_{\mathbf{k}_{1}}+\varepsilon_{\mathbf{k}_{2}}-\varepsilon_{\mathbf{p}}-\varepsilon_{\mathbf{k}_{1}+\mathbf{k}_{2}-\mathbf{p}}}~{}\\!.$
(38)

The expression (38) is very similar555Recall the standard rule
$\sum_{\mathbf{k}}\rightarrow V\int d^{3}\mathbf{k}/(2\pi)^{3}\equiv
V\int_{\mathbf{k}}$ for passing from the box normalization used in [2, 4] to
the continuum one used here. to the formula (5) used in [4] as the second
order contribution to the system’s internal energy density $u$, except that
the latter has an extra factor of 2. The foundation of the formula for
$f=u-Ts$ (which, apart from this factor of 2, is equivalent to our one) used
in [4] is, however, somewhat unclear: to obtain their second order correction
to the energy density $u$ these authors took the expression (15) given in
Section 11.4 of [2] which is obtained by simply using the finite temperature
distributions in place of the zero temperature ones in the ordinary second
order correction to the ground state energy of the system and have taken the
entropy density $s$ as given by the zeroth order textbook formula. In contrast
our expression (38) results from a systematic, well founded expansion and the
coefficient in (38) is unambiguously fixed by the cancellation of the
divergence.

After integrating over the cosine of the angle between $\mathbf{p}$ and
$\mathbf{k}_{1}+\mathbf{k}_{2}$ one can write the resulting expression in the
form

$\displaystyle
F^{(2)}=-V~{}\\!{C_{0}^{2}m_{f}\over(2\pi)^{2}\hbar^{2}}\\!\int_{\mathbf{k}_{1}}\int_{\mathbf{k}_{2}}{n_{+}(\mathbf{k}_{1})~{}\\!n_{-}(\mathbf{k}_{2})\over|\mathbf{k}_{1}+\mathbf{k}_{2}|}\int_{0}^{\infty}\\!dp~{}\\!p\left[n_{+}(\mathbf{p})+n_{-}(\mathbf{p})\right]\phantom{aaaaaaaaa}$
$\displaystyle\times\ln\\!\left|{(p-\Delta_{+})(p-\Delta_{-})\over(p+\Delta_{+})(p+\Delta_{-})}\right|,$
(39)

in which

$\displaystyle\Delta_{\pm}\equiv{1\over
2}~{}\\!|\mathbf{k}_{1}+\mathbf{k}_{2}|\pm{1\over
2}~{}\\!|\mathbf{k}_{1}-\mathbf{k}_{2}|~{}\\!.$

It is clear that the singularity at $|\mathbf{k}_{1}+\mathbf{k}_{2}|=0$ in
(39) is spurious: if $\mathbf{k}_{1}+\mathbf{k}_{2}=\mathbf{0}$ then
$\Delta_{-}=-\Delta_{+}$ and the innermost integral vanishes.

##  5 Numerical evaluation

The most difficult part of the computation is the accurate and efficient
numerical evaluation of the multiple integrals in the expression (39).
Rescaling the momentum integration variables $\mathbf{k}_{1}=k_{\rm
F}\mathbf{v}_{1}$, etc. and inserting $C_{0}=(4\pi\hbar/m_{f})a_{0}$ one can
write the second order contribution to the right hand side of (25) as

$\displaystyle{6\pi^{2}\over k_{\rm F}^{3}}~{}\\!{F^{(2)}\over\varepsilon_{\rm
F}V}=-(k_{\rm
F}a_{0})^{2}~{}\\!{6\over\pi^{2}}\int_{0}^{\infty}\\!dv_{1}~{}\\!v_{1}^{2}~{}\\!n(v_{1},\nu_{+},t)\\!\int_{0}^{\infty}\\!dv_{2}~{}\\!v_{2}^{2}~{}\\!n(v_{2},\nu_{-},t)$
$\displaystyle\times\sum_{\sigma=\pm}\sum_{\sigma^{\prime}=\pm}\\!\int_{-1}^{1}\\!d\xi~{}\\!{I(\Delta_{\sigma},\nu_{\sigma^{\prime}},t)\over\sqrt{v_{1}^{2}+v_{2}^{2}+2\xi
v_{1}v_{2}}}~{}\\!,\phantom{aaa}$ (40)

where $\nu_{\pm}=\mu^{(0)}_{\pm}/k_{\rm B}T\equiv\delta_{\pm}/t$
($\delta_{\pm}=\mu^{(0)}_{\pm}/\varepsilon_{\rm F}$),

$\displaystyle n(v,\nu,t)=\left[1+\exp\\!\left({v^{2}\over
t}-\nu\right)\right]^{-1},$ $\displaystyle
I(\Delta,\nu,t)=\int_{0}^{\infty}\\!du~{}\\!u~{}\\!n(u,\nu,t)~{}\\!\ln\\!\left|{u-\Delta\over
u+\Delta}\right|.$

and

$\displaystyle\Delta_{\pm}(v_{1},v_{2},\xi)={1\over
2}\sqrt{v_{1}^{2}+v_{2}^{2}+2\xi v_{1}v_{2}}\pm{1\over
2}\sqrt{v_{1}^{2}+v_{2}^{2}-2\xi v_{1}v_{2}}~{}\\!.$

The trick allowing to realize the numerical computation is to make first, for
fixed values of $t$ (temperature) and $P$ (the system’s polarization) which
together, through (20) determine $\nu_{+}$ and $\nu_{-}$, an interpolation of
the functions $I(|\Delta|,\nu_{+},t)$ and $I(|\Delta|,\nu_{-},t)$ (because,
obviously, $I(-|\Delta|,\nu_{\pm},t)=-I(|\Delta|,\nu_{\pm},t)$) in the
variable $w=1/(1+|\Delta|)$ (to interpolate on the compact interval $[0,1]$)
and then performing numerically the integrations over $v_{1}$, $v_{2}$ and
$\xi$ using these interpolations. In the actual code written in the Python
programming language the functions $I(|\Delta|,\nu_{\pm},t)$ are evaluated
with the help of the adaptive integration routine (scipy.integrate.quad; the
integration domain is split ted into three subdomains to accurately handle the
logarithmic singularity - in the relevant regions near $w=1/(1+\Delta)\equiv
w_{0}$ we substitute $r^{3}=|w-w_{0}|$ so that the integrand behaves like
$r^{2}\ln(r)$ and can be treated using the quadrature methods - and its sharp
falloff, especially for small temperatures $t$, near $u^{2}=t\nu$ of the
distribution $n(u,\nu,t)$) and then interpolated using the cubic spline
interpolation routines of Python. The remaining triple integral over $v_{1}$,
$v_{2}$ and $\xi$ are again performed with the help of the Clenshaw-Curtis
quadrature in the variables $w_{1,2}=1/(1+v_{1,2})$ (again to have a compact
integration domain and again splitting it into subdomains to better handle the
regions $v_{1}^{2}\approx t\nu_{+}$ and $v_{2}^{2}\approx t\nu_{-}$); the
spurious singularity at $|\mathbf{v}_{1}+\mathbf{v}_{2}|=0$ is taken care of
by simply taking somewhat different numbers for the $v_{1}$ and $v_{2}$ grids.

To check the correctness of the code we have first compared its results for
$t\rightarrow 0$ (replacing the distributions $n(u,\nu,t)$ by the Heaviside
theta functions) with the second order correction $E^{(2)}_{\Omega}$ to the
system’s ground-state energy which as a function of $P$ is known analytically
[10, 9] (the function $J_{K}(x,y)$ is given e.g. by the formula (4) in [12]):

$\displaystyle{6\pi^{2}\over k_{\rm
F}^{3}}~{}\\!{E^{(2)}_{\Omega}\over\varepsilon_{\rm F}V}=(k_{\rm
F}a_{0})^{2}~{}\\!{6\over
5\pi^{2}}~{}\\!J_{K}((1+P)^{1/3},~{}\\!(1-P)^{1/3})~{}\\!.$

At $P=0$ (equal densities of spin up and spin down fermions) $J_{K}=4(11-\ln
4)/21$ and the right hand side of the above formula (setting in all these
comparisons $k_{\rm F}a_{0}=1$) equals $0.22264482$ while the Python code for
the right hand side of (40) gives the value $0.22264522$. For $P=0.5$ the code
gives $0.17184256$ to be compared with $0.17184207$ while at $P=0.9$ the
numbers to be compared are $0.046470057$ and $0.046470077$ (at $P=1$ both are
zero reflecting the impossibility of the $s$-wave interactions of two fermions
in the same spin state). For nonzero temperatures the results obtained using
the Clenshaw-Curtis quadrature have been compared with the ones obtained using
the more accurate (but more time consuming) adaptive integration routine. The
comparison shows that the relative uncertainty $\Delta_{F}$ (the difference of
the results of the two methods divided by their mean) is typically of order
$10^{-5}$, varying rather irregularity with $P$ and increasing somewhat with
$t$; in our further estimates we set $\Delta_{F}=10^{-5}$ for
$t\stackrel{{\scriptstyle<}}{{{}_{\sim}}}0.1$, $\Delta_{F}=1.5\times 10^{-5}$
for $0.1<t\leq 0.2$ and $\Delta_{F}=2\times 10^{-5}$ for $0.2<t$. While this
accuracy superficially looks quite satisfactory, it is, nevertheless, barely
sufficient: for values of the parameters ($t$ and/or $k_{\rm F}a_{0}$) at
which spontaneous ordering appears there is a very delicate cancellation
between different contributions to $F$ and the (relative) error of order
$10^{-5}$ in $F^{(2)}$ can, and in some cases indeed does, lead do the
appearances of very shallow fake minimum near $P=0$.

##  6 Results

Figure 4: Plots of the differences $F(P)-F(0)$ in units of $(k_{\rm
F}^{3}/6\pi^{2})(\hbar^{2}k_{\rm F}^{2}/2m_{f})$ of the system of spin $1/2$
fermions as a function of the order parameter $P$ for two representative
values of the temperature $t\equiv T/T_{\rm F}$ as obtained in the second
order of the perturbative expansion. Left: $t=0.1$; the successive curves
(from below) correspond to $k_{\rm F}a_{0}=1.0718$ (the lowest, blue, line),
$k_{\rm F}a_{0}=1.0723$ (yellow), $k_{\rm F}a_{0}=1.0728$ (green), $1.0733$
(red) and $1.0738$ (the highest, blue, line). Right: $t=0.15$; the successive
curves (from below) correspond to $k_{\rm F}a_{0}=1.0978$ (the lowest, blue,
line), $k_{\rm F}a_{0}=1.983$ (yellow), $k_{\rm F}a_{0}=1.0988$ (green),
$1.0993$ (red) and $1.0998$ (the highest, blue, line).

For a fixed value of the temperature, the system’s free energy $F$ as a
function of the polarization (and of the parameter $k_{\rm F}a_{0}$ in which,
in the approximation to which our analysis is restricted, it is a polynomial
of the second order) can be efficiently obtained by evaluating numerically the
integrals in (40) for several values of $P$ and constructing the cubic spline
interpolation. The resulting free energy differences, $F(P)-F(0)$, are plotted
in Figure 4 as functions of the polarization $P$ for two temperatures: $t=0.1$
and $0.15$ and several values of $k_{\rm F}a_{0}$ (obtained by constructing
the interpolation based on 11 points in $P$ only).

In view of the mentioned uncertainty in the computation of $F^{(2)}$ the
critical value of $k_{\rm F}a_{0}$ and the value of the polarization $P$ at
the transition must be determined by requiring that the value of $F$ at a
minimum developing away from $P=0$ differs from the one at $P=0$ at least by
$\Delta_{F}F^{(2)}(0)$. In this way one can properly handle the mentioned fake
minima close to $P=0$ one of which can be observed in the right panel of Fig.
4 (for $t=0$ that such a minimum is indeed produced by the inaccuracies of the
numerical code can be substantiated by comparing with the analytically known
dependence of the ground-state energy on $P$). The actual procedure which has
been adopted to determine the polarization and its uncertainty is as follows.
For a fixed value of the parameter $k_{\rm F}a_{0}$ (which is successively
increased from 0 in steps $\Delta(k_{\rm F}a_{0})=0.001$) the values of $F$ on
a preliminary grid of $P$-values $P_{n}=n\Delta P$ with $n=0,\dots,n_{\rm
max}=32$ are obtained. If the minimal value of $F$ occurs for $n_{\rm
min}=n_{\rm max}$ the polarization is taken as the maximal ($P=1$); if $n_{\rm
min}=0$ or $|F(P_{n_{\rm min}})-F(0)|\leq\Delta_{F}F^{(2)}(0)$ the
polarization is taken as vanishing ($P=0$). If $n_{\rm min}\neq 0,n_{\rm max}$
and $|F(P_{n_{\rm min}})-F(0)|>\Delta_{F}F^{(2)}(0)$ the polarization is taken
as truly nonvanishing. If it is nonvanishing for the first time (as far as the
increasing values of $k_{\rm F}a_{0}$ are concerned), one determines $n_{\rm
down}$ and $n_{\rm up}$ such that $|F(P_{n_{\rm min}})-F(P_{n_{\rm
down/up}})|<2\Delta_{F}F^{(2)}(0)$ (of course, $n_{\rm down}=0$ and/or $n_{\rm
up}=n_{\rm max}$ if these criteria cannot not be fulfilled for intermediate
values of $n$) and finds the values of $F$ on a finer grid of $P$-values with
$|P_{j+1}-P_{j}|=0.0001$ and $P_{n_{\rm down}}\leq P_{j}\leq P_{n_{\rm up}}$.
If the minimum of $F$ found on the finer grid occurs for $P_{j_{\rm
min}}<0.02$, it is assumed that it is a numerical artifact and the
polarization is taken as vanishing. In the opposite case the polarization is
taken to be nonvanishing and that value of $k_{\rm F}a_{0}$ is recorded as the
critical one (for the considered temperature). In this case on the finer grid
one seeks a range $(P_{j_{\rm down}},~{}\\!P_{j_{\rm up}})$ of $P$ around
$P_{j_{\rm min}}$ in which $|F(P_{j_{\rm
min}})-F(P_{j})|>\Delta_{F}F^{(2)}(0)$ for $P_{j_{\rm down}}\leq P_{j}\leq
P_{j_{\rm up}}$; if such a range cannot be found the transition is classified
as continuous (the polarization at the considered temperature is assumed to
increase continuously from zero as $k_{\rm F}a_{0}$ is increased) while if a
nontrivial range is obtained, the transition is classified as first order and
$P_{j_{\rm min}}-P_{j_{\rm down}}$ and $P_{j_{\rm up}}-P_{j_{\rm min}}$ are
taken as the uncertainties of the determination of the polarization right at
the transition. For values of $k_{\rm F}a_{0}$ higher than the critical one
(determined as described above for the considered temperature) $F$ is
evaluated on a finer grid of points $P_{j}$ with $P_{n_{\rm min}-1}\leq
P_{j}\leq P_{n_{\rm min}+1}$ and $P_{j+1}-P_{j}=0.001$ and the corresponding
polarization is determined as the position of the minimum of $F$ on this finer
grid. In this way one finds that $(k_{\rm F}a_{0})_{\rm cr}=1.05409$ at $t=0$
(which perfectly agrees with the known value obtained by computing the
system’s ground-state energy [11] and with [4]), $(k_{\rm F}a_{0})_{\rm
cr}=1.05858$ at $t=0.05$, $(k_{\rm F}a_{0})_{\rm cr}=1.07282$ at $t=0.1$,
$(k_{\rm F}a_{0})_{\rm cr}=1.09881$ at $t=0.15$ and $(k_{\rm F}a_{0})_{\rm
cr}=1.13845$ at $t=0.2$. The corresponding values of the polarization right at
the transition point are $P_{\rm cr}=0.575_{-0.019}^{+0.017}$ (again in
agreement with the value found in [11]) $0.558_{-0.017}^{+0.017}$,
$0.477_{-0.021}^{+0.019}$, $0.325_{-0.048}^{+0.035}$ and
$0.197_{-0.096}^{+0.045}$.

Figure 5: Polarization $P=(N_{+}-N_{-})/N$ of the system of spin $1/2$
fermions with a short range repulsive interaction obtained from the free
energy $F$ computed up to the second order of the perturbative expansion. In
the left panel as a function of the “gas parameter” $k_{\rm F}a_{0}$ for
several values of the temperature (counting from the left): $t\equiv
T/T_{F}=0$, $0.1$, $0.15$ and $0.2$ ($T_{\rm F}\equiv\varepsilon_{\rm
F}/k_{\rm B}$). Marked are also uncertainties of the value of $P$ right at the
transition points. In the right panel as a function of the temperature for
several fixed values of $k_{\rm F}a_{0}$.

The dependence of the polarization as a function of the “gas parameter”
$k_{\rm F}a_{0}$ is shown, for a few values of the temperature $t$, in the
left panel of Figure 5. This is essentially the same plot as the one presented
in [4] (the agreement with the critical values of the gas parameters at
successive temperatures that can be read off from the plot there seems to be
quite good) except that in Figure 5 marked are also the uncertainties in the
determination (following from the procedure just described) of the
polarization right at the transition.

Owing to the efficiency of our numerical code (stemming basically from the
trick with the interpolations) the procedure of finding the polarization of
the system described above can be applied also at fixed values of $k_{\rm
F}a_{0}$ (replacing the grid in $k_{\rm F}a_{0}$ by a one in $t$). The
resulting polarization of the system as a function of the temperature for
several fixed values of the gas parameter is shown in the right panel of
Figure 5. Knowing the polarization as a function of the other parameters it is
possible to construct the free energy $F(T,V,N)\equiv F(T,V,N,P(T,N/V))$ for
several values of $k_{\rm F}a_{0}$ and to determine also other thermodynamic
characteristics of the system. For example, using the grid in $t$ the second
derivative of the free energy $F(T,V,N))$ with respect to the temperature can
in principle be obtained yielding the system’s heat capacity. The result of
such an exercise is shown in Figure 6 for two values of the “gas parameter”.
It is shows that the discontinuity of the heat capacity at the transition
point grows with the value of $k_{\rm F}a_{0}$ (i.e. also with the increasing
temperature, if $k_{\rm F}a_{0}$ is varied). However for higher values of
$k_{\rm F}a_{0}$ the numerical inaccuracies do not allow for a reliable
computation. Indeed, as the transition at higher temperatures becomes
continuous a divergence of the heat capacity probably starts to build up
making the numerical computation of the second derivative of the free energy
unstable for $t\stackrel{{\scriptstyle>}}{{{}_{\sim}}}0.12$. Similarly, it is
in principle possible to determine the system’s polarization taking into
account an infinitesimally weak external magnetic field (this as explained
influences only the determination of the zero-th order chemical potentials
$\tilde{\mu}^{(0)}_{\pm}$ from the conditions (22)) and to compute the
system’s magnetic susceptibility $\chi_{T}$ by constructing the derivative of
the polarization with respect to ${\cal H}$. While such a computation seems to
indicate that at least at low temperatures, at which the transition is (in the
approximation to which our computation is restricted) is first order, the
susceptibility also has a finite discontinuity at the transition point, it is
not sufficiently stable numerically to yield reliable values of $\chi_{T}$ and
it is probably more practical to obtain it by computing the (connected) two
point correlation function from the formula
$\chi_{T}=(\beta/V)\tilde{G}^{(2)}_{\rm con}(\mathbf{0})$. We do not attempt
this here.

Figure 6: Heat capacity (in units of $Nk_{\rm B}$) of the system of spin $1/2$
fermions with a short range repulsive interaction as a function of the
temperature for two different fixed values of the parameter $k_{\rm F}a_{0}$
obtained from the free energy $F$ computed up to the second order.

##  7 Conclusions

We have developed a systematic perturbative expansion of the grand
thermodynamic potential $\Omega$ and of the free energy $F$ of the system of
(nonrelativistic) interacting spin 1/2 fermions. We have applied this
expansion within the effective field theory in which the underlying repulsive
spin-independent binary interaction of fermions is replaced by an infinite
number of contact interaction terms and which allows to directly express
computed quantities in terms of the scattering lengths and effective radii
which characterize the underlying interaction potential. We have shown (up to
the third order but the result seems to be valid in general) that to the
expansion of the free energy effectively contribute only those Feynman
diagrams which give nonvanishing contributions to the ground-state energy of
the system evaluated at zeroth order chemical potentials (associated with spin
up and spin down fermions).

Our numerical analysis has been restricted here to the first nontrivial order
of the perturbative expansion (i.e. the first one going beyond the textbook
mean field approximation) in which the results are still universal, i.e. on
the form of the underlying interaction depend only through the $s$-wave
scattering length $a_{0}$ (in the next order the results start to depend also
on the $p$-wave scattering length $a_{1}$ and the effective radius $r_{0}$).
We have devised a method for efficient numerical evaluation of the requisite
nested integrals and used it to compute the system’s polarization and its
value right at the transition point paying attention to the uncertainty of the
determination of the latter quantity which is crucial is assessing the
character of the transition. For low temperatures,
$T\stackrel{{\scriptstyle<}}{{{}_{\sim}}}0.1~{}\\!T_{\rm F}$, we have also
managed to determine the system’s heat capacity encountering, however, some
problems with the accuracy of numerical evaluation of the derivatives of the
free energy which seem to prevent obtaining (at least without substantial
improvements in the method) reliable values of the heat capacity for higher
temperatures as well as determining the system’s magnetic susceptibility.

Of course, since the perturbative computation of the system’s ground-state
energy agrees with the results obtained (for specific forms of the underlying
interaction) using the Quantum Monte Carlo approach only for $k_{\rm
F}a_{0}\stackrel{{\scriptstyle<}}{{{}_{\sim}}}0.5$, the results presented here
cannot be taken very seriously. Moreover it is now known that already the
inclusion of the third order corrections to the system’s ground-state energy
(free energy at zero temperature) significantly weaken the first order
character of the transition (at zero temperature) to the ordered state. For
these reasons our effort summarized here should be treated rather as a
preliminary step taken towards extending the computation to a higher order and
towards a possible implementation of a resummation of some class of the
contributions to the free energy in the spirit of the approach of [13]. Such a
resummation can probably also allow to overcome the limitation, inherent in
the effective field theory approach, to sufficiently small temperatures only:
as this approach relies on the clean separation of the scales ($R\ll k_{\rm
F}^{-1}$, where $R$ is the characteristic length of the underlying
interaction) it cannot be applied, at least if restricted to a finite order of
the perturbative expansion, when $k_{\rm B}T$ becomes comparable with the
energy scale set by $\varepsilon_{\rm F}$. We plan to return to these issues
in the forthcoming paper.

## Appendix A

The following summation formulae hold [7] (the limit $\eta\rightarrow 0^{+}$
is implicit):

$\displaystyle{1\over\beta}\sum_{n\in\mathbb{Z}}{e^{i\eta\omega^{F}_{n}}\over
i\omega^{F}_{n}-x}={1\over 1+e^{\beta
x}}~{}\\!,\phantom{aaa}\omega^{F}_{n}\equiv{\pi\over\beta}~{}\\!(2n+1)~{}\\!,$
(A.1)
$\displaystyle{1\over\beta}\sum_{n\in\mathbb{Z}}{e^{i\eta\omega^{B}_{n}}\over
i\omega^{B}_{n}-x}={1\over 1-e^{\beta
x}}~{}\\!,\phantom{aaa}\omega^{B}_{n}\equiv{2\pi\over\beta}~{}\\!n~{}\\!,$
(A.2)

and, by decomposing into simple fractions,

$\displaystyle{1\over\beta}\sum_{n\in\mathbb{Z}}{1\over
i\omega^{F}_{n}-x}~{}\\!{1\over
i\omega^{F}_{n+l}-y}={1\over\beta}\sum_{n\in\mathbb{Z}}{1\over
i\omega^{F}_{n}-x}~{}\\!{1\over
i\omega^{F}_{n}-(y-i\omega_{l}^{B})}\phantom{aaaaaaa}$ $\displaystyle={1\over
i\omega^{B}_{l}-(y-x)}\left({1\over 1+e^{\beta x}}-{1\over 1+e^{\beta
y}}\right),\phantom{aa}$ (A.3)
$\displaystyle{1\over\beta}\sum_{n\in\mathbb{Z}}{1\over
i\omega^{F}_{n}-x}~{}\\!{1\over
i\omega^{F}_{l-n}-y}=-{1\over\beta}\sum_{n\in\mathbb{Z}}{1\over
i\omega^{F}_{n}-x}~{}\\!{1\over
i\omega^{F}_{n}-(i\omega^{B}_{l+1}-y)}\phantom{aaaaa}$ $\displaystyle={1\over
i\omega^{B}_{l+1}-(y+x)}\left({1\over 1+e^{\beta x}}-{1\over 1+e^{-\beta
y}}\right).$ (A.4)

Similarly,

$\displaystyle{1\over\beta}\sum_{l\in\mathbb{Z}}{1\over
i\omega^{B}_{l}-x}~{}\\!{1\over i\omega^{B}_{l}-y}={1\over x-y}\left({1\over
1-e^{\beta x}}-{1\over 1-e^{\beta y}}\right).$ (A.5)

Useful can be also the formula

$\displaystyle{1\over x-a_{1}}~{}\\!{1\over x-a_{2}}~{}\\!\dots~{}\\!{1\over
x-a_{n}}=\sum_{l=1}^{n}\left(\prod_{k\neq l}^{n}{1\over
a_{l}-a_{k}}\right){1\over x-a_{l}}~{}\\!.$ (A.6)

## Appendix B

Here we demonstrate the cancellation of the additional terms in the formula
(31) for $F^{(3)}$ against the contributions of diagrams which vanish at zero
temperature. Analogous cancellation of the contribution of the diagrams shown
in Figure 2 against the additional terms in (30) has been checked in the main
text. It will be convenient to introduce first the following notation:

$\displaystyle n_{x}\phantom{a}~{}\\!\equiv\int_{\mathbf{k}}\\!{1\over
1+e^{\beta(\varepsilon_{\mathbf{k}}-x_{0})}}~{}\\!,$ $\displaystyle
n_{xx}~{}\equiv\int_{\mathbf{k}}\\!{1\over[1+e^{\beta(\varepsilon_{\mathbf{k}}-x_{0})}]^{2}}~{}\\!,$
$\displaystyle
n_{xxx}\equiv\int_{\mathbf{k}}\\!{1\over[1+e^{\beta(\varepsilon_{\mathbf{k}}-x_{0})}]^{3}}~{}\\!,$

etc. From (11) one immediately obtains (all functions are evaluated at $x_{0}$
and $y_{0}$)

$\displaystyle\Omega^{(0)}_{x}=-V~{}\\!n_{x}~{}\\!,$
$\displaystyle\Omega^{(0)}_{xx}=-V\beta\left(n_{x}-n_{xx}\right),$ (B.1)
$\displaystyle\Omega^{(0)}_{xxx}=-V\beta^{2}\left(n_{x}-3n_{xx}+2n_{xxx}\right).$

Analogously can be written the derivatives of $\Omega^{(0)}$ with respect to
$y$. The necessary derivatives of $\Omega^{(1)}=C_{0}V~{}\\!n_{x}~{}\\!n_{y}$
take the form

$\displaystyle\Omega^{(1)}_{x}=C_{0}V\beta\left(n_{x}-n_{xx}\right)n_{y}~{}\\!,$
$\displaystyle\Omega^{(1)}_{y}=C_{0}V\beta~{}\\!n_{x}\left(n_{y}-n_{yy}\right),$
$\displaystyle\Omega^{(1)}_{xx}=C_{0}V\beta^{2}\left(n_{x}-3n_{xx}+2n_{xxx}\right)n_{y}~{}\\!,$
(B.2)
$\displaystyle\Omega^{(1)}_{yy}=C_{0}V\beta^{2}~{}\\!n_{x}\left(n_{y}-3n_{yy}+2n_{yyy}\right),$
$\displaystyle\Omega^{(1)}_{xy}=C_{0}V\beta^{2}\left(n_{x}-n_{xx}\right)\\!\left(n_{y}-n_{yy}\right).$

To $\Omega^{(3)}$, in addition to the two “mercedes-type” diagrams shown in
Figure 7 (the contributions $\Omega^{(3)a}$ and $\Omega^{(3)b}$), contribute
also the “mitsubishi-type” diagrams of Figure 8 (the contributions
$\Omega^{(3)c}$ and $\Omega^{(3)d}$), the two diagrams of Figure 9 (the
contributions $\Omega^{(3)e}$ and $\Omega^{(3)f}$) and the single “audi-type”
diagram of Figure 9 ($\Omega^{(3)g}$).

Figure 7: The particle-particle and the particle-hole diagrams (the “mercedes-
like” diagrams) contributing in the order $C^{3}_{0}$ to $\Omega^{(3)}$.

The computation of $\Omega^{(3)c}$ and $\Omega^{(3)d}$ is straightforward (it
is analogous to that of $\Omega^{(2)b}$ and $\Omega^{(2)c}$) and yields

$\displaystyle\Omega^{(3)c}+\Omega^{(3)d}={1\over
6}~{}\\!C_{0}^{3}V\beta^{2}\left[\left(n_{x}-3n_{xx}+2n_{xxx}\right)n_{y}^{3}+n_{x}^{3}\left(n_{y}-3n_{yy}+2n_{yy}\right)\right].$

One then readily sees that in (31) this is cancelled by the last two terms:

$\displaystyle\Omega^{(3)a}+\Omega^{(3)b}-{\Omega^{(0)}_{xxx}[\Omega^{(1)}_{x}]^{3}\over
6~{}\\![\Omega^{(0)}_{xx}]^{3}}-{\Omega^{(0)}_{yyy}[\Omega^{(1)}_{y}]^{3}\over
6~{}\\![\Omega^{(0)}_{yy}]^{3}}=0~{}\\!.$

One has now to consider the terms
$-\Omega^{(2)}_{x}\Omega^{(1)}_{x}/\Omega^{(0)}_{xx}$ and
$-\Omega^{(2)}_{y}\Omega^{(1)}_{y}/\Omega^{(0)}_{yy}$ in (31). $\Omega^{(2)}$
is given by three diagrams shown in Figure 3 left and 2. It is convenient to
write the contribution the first one, $\Omega^{(2)a},$ in the form

$\displaystyle\Omega^{(2)a}=-{C_{0}^{2}V\over
2}~{}\\!{1\over\beta}\sum_{l\in\mathbb{Z}}\\!\int_{\mathbf{q}}\\!\int_{\mathbf{k}}\\!\int_{\mathbf{p}}{1\over\beta}\sum_{n\in\mathbb{Z}}{1\over
i\omega^{F}_{n}-(\varepsilon_{\mathbf{k}}-x)}~{}\\!{1\over
i\omega^{F}_{n+l}-(\varepsilon_{\mathbf{k}+\mathbf{q}}-x)}\phantom{a}$
$\displaystyle\times{1\over\beta}\sum_{m\in\mathbb{Z}}{1\over
i\omega^{F}_{m}-(\varepsilon_{\mathbf{p}}-y)}{1\over
i\omega^{F}_{m-l}-(\varepsilon_{\mathbf{p}-\mathbf{q}}-y)}~{}\\!.$ (B.3)

Differentiating it with respect to $x$ one obtains the expression which is a
sum of the two terms

$\displaystyle\Omega^{(2)a}_{x}={C_{0}^{2}V\over
2}~{}\\!{1\over\beta}\sum_{l\in\mathbb{Z}}\\!\int_{\mathbf{q}}\\!\int_{\mathbf{k}}\\!\int_{\mathbf{p}}{1\over\beta}\sum_{m\in\mathbb{Z}}{1\over
i\omega^{F}_{m}-(\varepsilon_{\mathbf{p}}-y)}{1\over
i\omega^{F}_{m-l}-(\varepsilon_{\mathbf{p}-\mathbf{q}}-y)}\phantom{aaaaa}$
$\displaystyle\times{1\over\beta}\sum_{n\in\mathbb{Z}}\left\\{{1\over[i\omega^{F}_{n}-(\varepsilon_{\mathbf{k}}-x)]^{2}}~{}\\!{1\over
i\omega^{F}_{n+l}-(\varepsilon_{\mathbf{k}+\mathbf{q}}-x)}\right.\phantom{a}~{}$
$\displaystyle+\left.{1\over
i\omega^{F}_{n}-(\varepsilon_{\mathbf{k}}-x)}~{}\\!{1\over[i\omega^{F}_{n+l}-(\varepsilon_{\mathbf{k}+\mathbf{q}}-x)]^{2}}\right\\}.$

These two terms are equal - to see this, it suffices to set in the second term
$\mathbf{q}=-\mathbf{q}^{\prime}$,
$\mathbf{k}=\mathbf{k}^{\prime}+\mathbf{q}^{\prime}$,
$\mathbf{p}=\mathbf{p}^{\prime}-\mathbf{q}^{\prime}$ and $l=-l^{\prime}$,
$n=n^{\prime}+l^{\prime}$ $m=m^{\prime}-l^{\prime}$. Thus, multiplying this by
$-\Omega^{(1)}_{x}/\Omega^{(0)}_{xx}$ which simply equals $C_{0}~{}\\!n_{y}$
one readily sees (taking into account that ${\cal
G}_{--}^{(0)}(0,\mathbf{0})=n_{y}$) that this precisely cancels
$\Omega^{(3)e}$. Thus, in $F^{(3)}$

$\displaystyle\Omega^{(3)e}+\Omega^{(3)f}-{\Omega^{(1)}_{x}\over\Omega^{(0)}_{xx}}~{}\\!\Omega^{(2)a}_{x}-{\Omega^{(1)}_{y}\over\Omega^{(0)}_{yy}}~{}\\!\Omega^{(2)a}_{y}=0~{}\\!,$

Figure 8: Two order $C^{3}_{0}$, “mitsubishi-type” diagrams contributing g to
$\Omega^{(3)}$.

After these cancellations one is left with

$\displaystyle
F^{(3)}=\Omega^{(3)a}+\Omega^{(3)b}+\Omega^{(3)g}+C_{0}~{}\\!n_{y}\left(\Omega^{(2)b}_{x}+\Omega^{(2)c}_{x}\right)+C_{0}~{}\\!n_{x}\left(\Omega^{(2)b}_{y}+\Omega^{(2)c}_{y}\right)$
$\displaystyle+~{}\\!{\Omega^{(1)}_{xx}[\Omega^{(1)}_{x}]^{2}\over
2~{}\\![\Omega^{(0)}_{xx}]^{2}}+{\Omega^{(1)}_{yy}[\Omega^{(1)}_{y}]^{2}\over
2~{}\\![\Omega^{(0)}_{yy}]^{2}}+{\Omega^{(1)}_{xy}\Omega^{(1)}_{x}\Omega^{(1)}_{y}\over\Omega^{(0)}_{xx}\Omega^{(0)}_{yy}}~{}\\!.\phantom{aaaaaaaaaa}$

Explicitly the last line of $F^{(3)}$ reads

$\displaystyle{1\over
2}~{}\\!C_{0}^{3}V\beta^{2}\left\\{\left(n_{x}-3n_{xx}+2n_{xxx}\right)n_{y}^{3}\right.\phantom{aaaaaaaaaaaaaaaaaaaaaaaaaaa}~{}$
$\displaystyle\left.+n_{x}^{3}\left(n_{y}-3n_{yy}+2n_{yyy}\right)+2~{}\\!n_{x}\left(n_{x}-n_{xx}\right)\\!\left(n_{y}-n_{yy}\right)n_{y}\right\\},$

while the contribution $\Omega^{(3)g}$ of the “audi-type” diagram of Figure 8
can be written as

$\displaystyle\Omega^{(3)g}=C_{0}^{3}V\beta^{2}~{}\\!n_{x}\left(n_{x}-n_{xx}\right)\\!\left(n_{y}-n_{yy}\right)n_{y}~{}\\!.$

Finally,

$\displaystyle\Omega^{(2)b}_{x}+\Omega^{(2)c}_{x}=-{1\over
2}~{}\\!C_{0}^{2}V\beta^{2}\left[\left(n_{x}-3n_{xx}+2n_{xxx}\right)n_{y}^{2}+2~{}\\!n_{x}\left(n_{x}-n_{xx}\right)\\!\left(n_{y}-n_{yy}\right)\right],$

(the sum $\Omega^{(2)b}_{y}+\Omega^{(2)c}_{y}$ is given by an analogous
expression) and after a straightforward algebra all the extra terms cancel
out, so that eventually, $F^{(3)}=\Omega^{(3)a}+\Omega^{(3)b}$ that is, it
given solely by the “mercedes-like” diagrams evaluated at the zeroth order
chemical potentials $x_{0}$, $y_{0}$ (or $\tilde{x}_{0}$, $\tilde{y}_{0}$, if
there is an external magnetic field). The diagrams canceled by the extra terms
in the formulae (30) and (31) are precisely those (see e.g. [6]) which vanish
at zero temperature, that is do not contribute to the expansion of the formula
(13) for the ground state energy.

Figure 9: The contributions to $\Omega^{(3)e}$ and $\Omega^{(3)f}$.

One can also simplify the formula (29) for the second order correction $x_{2}$
to the $\mu_{+}$ chemical potential (and the analogous formula for $y_{2}$).
After a straightforward algebra one obtains

$\displaystyle x_{2}=-{\Omega^{(2)a}_{x}\over\Omega^{(0)}_{xx}}~{}\\!,$

all other terms neatly canceling. Of course, $x_{1}=C_{0}~{}\\!n_{-}$,
$y_{1}=C_{0}~{}\\!n_{+}$ but it is perhaps more instructive to write666This
form clearly shows, since the cancellation of the divergences in the sum
$\Omega^{(1)}+\Omega^{(2)a}$ has already been demonstrated, that the computed
perturbatively chemical potentials $x$ and $y$ are to this order finite, after
the cutoff dependence of the coupling $C_{0}$ is taken into account. The
argument obviously generalizes to all orders: if the free energy $F$ is made
finite by the renormalization of the couplings, so must be the chemical
potentials.

$\displaystyle
x=x_{0}-{1\over\Omega^{(0)}_{xx}}\left(\Omega^{(1)}_{x}+\Omega^{(2)a}_{x}+\dots\right).$

The cancellation found here is obviously necessary for the equivalence of two
ways of determining the system’s polarization. Indeed, only if this
cancellation holds is the minimization of the free energy written in the form

$\displaystyle
F=\Omega^{(0)}(x_{0},y_{0})+(x_{0}~{}\\!n_{x}+y_{0}~{}\\!n_{y})~{}\\!V+\Omega^{(1)}(x_{0},y_{0})+\Omega^{(2)a}(x_{0},y_{0})+\dots,$

with respect to $n_{x}$ (keeping $n_{y}=n-n_{x}$), which (taking into account
that $\Omega^{(0)}_{x}+n_{x}V=0$, $\Omega^{(0)}_{y}+n_{y}V=0$) amounts to

$\displaystyle(x_{0}-y_{0})~{}\\!V=-\left[\Omega^{(1)}_{x}+\Omega^{(2)a}_{x}+\dots\right]{\partial
x_{0}\over\partial
n_{x}}+\left[\Omega^{(1)}_{y}+\Omega^{(2)a}_{y}+\dots\right]{\partial
y_{0}\over\partial n_{y}}~{}\\!.$

equivalent to the condition $\mu_{+}=\mu_{-}$ written in the form
$x_{0}+x_{1}+x_{2}+\dots=y_{0}+y_{1}+y_{2}+\dots$, that is

$\displaystyle
x_{0}-y_{0}={\Omega^{(1)}_{x}+\Omega^{(2)a}_{x}+\dots\over\Omega^{(0)}_{xx}}-{\Omega^{(1)}_{y}+\Omega^{(2)a}_{y}+\dots\over\Omega^{(0)}_{yy}}~{}\\!.$

The equivalence follows from noticing that since $n_{x}=-\Omega^{(0)}_{x}/V$,
$n_{y}=-\Omega^{(0)}_{y}/V$, the derivatives of $x_{0}$ and $y_{0}$ are
precisely equal to

$\displaystyle{\partial x_{0}\over\partial
n_{x}}=-{V\over\Omega^{(0)}_{xx}}~{}\\!,\phantom{aaa}{\partial
y_{0}\over\partial n_{y}}=-{V\over\Omega^{(0)}_{yy}}~{}\\!.$

From this argument it immediately follows that
$x_{3}=-(\Omega^{(3)a}_{x}+\Omega^{(3)b}_{x})/\Omega^{(0)}_{xx}$ and
$y_{3}=-(\Omega^{(3)a}_{y}+\Omega^{(3)b}_{y})/\Omega^{(0)}_{yy}$.

Figure 10: The “audi”-type, order $C^{3}_{0}$ contribution to $\Omega^{(3)}$.

Restricted to the first order the left hand side of the equality
$x_{0}+x_{1}=y_{0}+y_{1}$ reads

$\displaystyle x_{0}+x_{1}+\dots=k_{\rm
B}T~{}\\!f^{-1}\\!\left(\left({\varepsilon^{(0)}_{\rm F}(n_{+})\over k_{\rm
B}T}\right)^{3/2}\right)+C_{0}~{}\\!n_{-}+\dots$

The right hand side is given by the analogous formula. If one sets here
$n_{\pm}=(N/2V)(1\pm P)$, this reproduces the condition (21).

## References

  * [1] K. Huang, Statistical Mechanics, John Willey and Sons, Inc., New York 1963.
  * [2] R.K. Pathria, Statistical Mechanics, Pergamon Press, Oxford, 1972.
  * [3] E. Stoner, Phil. Mag. 15 (1933), 1018.
  * [4] R.A. Duine and A.H. MacDonald, Phys. Rev. Lett. 95, 230403 (2005).
  * [5] D. Belitz, T.R. Kirkpatrick and T. Vojta, Phys. Rev. Lett. 82, 4707 (1999); D. Belitz and T.R. Kirkpatrick, Phys. Rev. Lett. 89, 247202 (2002).
  * [6] H.-W. Hammer, R. J. Furnstahl, Nucl. Phys. A678, 277 (2000); arXiv:0004043 [nucl-th].
  * [7] A.L. Fetter and J.D. Walecka, Quantum Theory of Many Particle Systems, McGraw Hill, 1971.
  * [8] P.H. Chankowski and J. Wojtkiewicz, Phys. Rev. B104 144425 (2021) arXiv:2108.00793 [cond-mat.quant-gas].
  * [9] J. Pera, J. Casulleras and J. Boronat, arXiv:2205.13837 [cond-mat.quant-gas], SciPost. Phys. 14, 038 (2023).
  * [10] S. Kanno, Prog. Theor. Phys. 44, 813 (1970).
  * [11] P.H. Chankowski and J. Wojtkiewicz, Acta. Phys. Pol. B53 (2022), 9-A4, arXiv:2303.05076 [cond-mat.quant-gas].
  * [12] P.H. Chankowski, J. Wojtkiewicz and S. Augustynowicz, Phys. Rev. A107 (2023), 063311, arXiv:2206.09921 [cond-mat.quant-gas].
  * [13] L. He and X.-G. Huang, Phys. Rev. A85, 043624 (2012), arXiv:1106.1345; L. He, Ann. of Phys. 351, 477 (2014), arXiv:1405.3338, Phys. Rev. A90, 053633, arXiv:1405.5242; L. He, X.-J. Liu, X.-G. Huang and H. Hui, Phys. Rev. A93, 063629 (2016), arXiv:1412.2412.
  * [14] S. Pilati, G. Bertaina, S. Giorgini and M. Troyer, Phys. Rev. Lett. 105, 030405 (2010), arXiv:1004/1169 [cond-mat.quant-gas].
  * [15] R.J. Furnstahl, H.-W. Hammer and S.J. Puglia Ann. of Phys. 322, 2703 (2007), arXiv:0612086 [nucl-th].
  * [16] C. Wellenhofer, C. Drischler and A. Schwenk, Phys. Lett. B 802 135247 (2020) arXiv:1812.08444 [nucl-th]; Phys. Rev. C104, 014003 (2021) arXiv:2102.05966 [cond-mat.quant-gas].