# Pixel Relationships-based Regularizer for Retinal Vessel Image Segmentation

Lukman Hakim1, Takio Kurita2∗ 1 Department of Information Engineering,
Graduate School of Engineering, Hiroshima University, Higashihiroshima, Japan
2 Graduate School of Advanced Science and Engineering, Hiroshima University,
Higashihiroshima, Japan  ∗E-mail<EMAIL_ADDRESS>
###### Abstract

The task of image segmentation is to classify each pixel in the image based on
the appropriate label. Various deep learning approaches have been proposed for
image segmentation that offers high accuracy and deep architecture. However,
the deep learning technique uses a pixel-wise loss function for the training
process. Using pixel-wise loss neglected the pixel neighbor relationships in
the network learning process. The neighboring relationship of the pixels is
essential information in the image. Utilizing neighboring pixel information
provides an advantage over using only pixel-to-pixel information. This study
presents regularizers to give the pixel neighbor relationship information to
the learning process. The regularizers are constructed by the graph theory
approach and topology approach: By graph theory approach, graph Laplacian is
used to utilize the smoothness of segmented images based on output images and
ground-truth images. By topology approach, Euler characteristic is used to
identify and minimize the number of isolated objects on segmented images.
Experiments show that our scheme successfully captures pixel neighbor
relations and improves the performance of the convolutional neural network
better than the baseline without a regularization term.

###### Index Terms:

neighboring, pixel, deep learning, graph, topology.

##  I Introduction

Image segmentation is an important task in the field of computer vision. The
task of image segmentation is to classify each pixel in the image based on the
appropriate label. Image segmentation applications have been widely used, such
as medical imagery (for example, Retinal Vessel Segmentation[3], Tumor
Segmentation[4], Breast Cancer Detection[5]), Autonomous Vehicles (Example:
Vehicle Detection[6], Pedestrian Detection[7]), Counting Objects[8], and Farm
Industry[9]. Several Image segmentation techniques have been developed, such
as thresholding, k-means clustering, and graph cut methods.

Recently, many deep learning techniques have been developed for image
segmentation that offers high accuracy and deep architecture[1][2]. However,
the deep learning technique uses a pixel-wise loss function for the training
process. Using pixel-wise loss neglected the pixel neighbor relationships in
the network learning process. The neighboring relationship of the pixels is
essential information in the image. Utilizing neighboring pixel information
provides an advantage over using only pixel-to-pixel information.

Figure 1: Illustration of graph based smoothing regularizer (GBS)
architecture. Figure 2: Representation of graph between two regions, a)
Ground-truth image, b) graph generated from background region and c) graph
generated from foreground region.

This study presents regularizers to provide the pixel neighbor relationship
information to the learning process. We have proposed three regularizers are
constructed by the graph theory approach and topology approach. Firstly, we
introduced a Graph-Based Smoothing Regularizer (GBS). The GBS considers the
graph laplacian from the foreground and background regions and then combines
it with the CNN baseline loss function. The combination of regularizers allows
the network to learn the pixel relationship efficiently. The illustration of
GBS is shown in Fig. 1.

Figure 3: Architecture of Graph Laplacian regularizer based on differences of
neighboring pixels.

Secondly, we introduced Graph Laplacian Regularization based on the
Differences of Neighboring Pixels (GLRDN) by constructing graph laplacian from
prediction and ground-truth images. A graph uses pixels as vertices and edges
defined by the ”differences” of neighboring pixels instead of similarities
between pixels. A graph uses pixels as vertices and edges defined by the
”differences” of neighboring pixels instead of similarities between pixels.
The basic idea is, if pair-wise pixels belonging a similar class, the
differences are small. Otherwise, the differences are significant if pair-wise
pixels are belonging a different class. The illustration of GLRDN is shown as
Fig 3.

Figure 4: Simple Illustration of Euler Characteristic Calculation.

Thirdly, we proposed a regularizer based on Euler Characteristic (EC) in the
segmented image. EC is a topological property of a shape in an image that can
identify the number of objects in the image. This EC-based regularizer is used
to identify the number of objects in the segmentation results. Then this
information is added to the objective function so that the network can
minimize the number of isolated objects in the segmented image.

To evaluate the effectiveness of the proposed methods, we implemented the
regularizers on retinal blood vessels segmentation in fundus images and
compared them with the baseline CNN without regularizers.

##  II Methods

This section explains three proposed regularizers we have done: GBS, GLRDN,
and EC Regularizer.

###  II-A Graph Based Smoothing (GBS)

The proposed graph-based smoothing regularizer is based on the relationship
between pixels[10]. The graph approach is used to obtain adjacency graphs and
diagonal matrices from the relationship information between pixels. The graph
can be transformed into a matrix by counting the edges between the two
adjacent vertices. The two vertices can be said to be adjacent if an edge
connects it. Then from this matrix, the Graph Laplacian is defined.

The GBS defines two different graphs for the foreground region and background
region. As shown in Fig. 2, let us consider the set of training samples
$X=\left(x_{m},t_{m}\right)|m=1,\ldots,M$ where $x_{m}$ is a $m^{th}$ input
image and $t_{m}$ is the $m^{th}$ target image, and $M$ is number of training
samples. In deep convolutional neural network, the network is trained to
predict the output image ${y}_{m}$ from the $m^{th}$ input image ${x}_{m}$.
For each edge of foreground $\left(j_{F},k_{F}\right)\in E_{F}$ and background
$\left(j_{B},k_{B}\right)\in E_{B}$ of the graph, the similarity
$\beta_{\left(j_{F},k_{F}\right)}$ and $\beta_{\left(j_{B},k_{B}\right)}$ is
defined as

$\beta_{(j_{F},k_{F})}=1-|t_{j_{F}}-t_{k_{F}}|$ (1)

$\beta_{(j_{B},k_{B})}=1-|t_{j_{B}}-t_{k_{B}}|$ (2)

We introduced the regularization term for smoothing $S$ based on foreground
region $F$ and background region $B$ as

$\sum_{(j_{F},k_{F})\in
G_{F}}\beta_{j_{F},k_{F}}{(y_{j_{F}}-y_{k_{F}})^{2}}=y^{T}(D_{F}-A_{F})y=y^{T}L_{F}y$
(3) $\sum_{(j_{B},k_{B})\in
G_{B}}\beta_{j_{B},k_{B}}{(y_{j_{B}}-y_{k_{B}})^{2}}=y^{T}(D_{B}-A_{B})y=y^{T}L_{B}y$
(4)

where, $L_{F}$ and $L_{B}$ is Graph Laplacian for foreground region and
background region, respectively. The adjacency and diagonal matrices is
defined as following

$\binom{A_{F}=\beta_{(j_{F},k_{F})},D_{F}=\sum_{j_{F}=1}^{N}\beta_{j_{F},k_{F}}}{A_{B}=\beta_{(j_{B},k_{B})},D_{B}=\sum_{j_{B}=1}^{N}\beta_{j_{B},k_{B}}}$
(5)

Graph based smoothing regularizer $R_{gbs}$ can be written as

$R_{gbs}=y^{T}(L_{F}+L_{B})y\\\ =y^{T}L_{G}y$ (6)

In these experiments, we use the Binary Cross Entropy (BCE) loss as objective
function. The BCE loss is given as below:

$E_{bce}=\sum_{i}^{M}t_{i}\log(y_{i})+(1-t_{i})\log(1-y_{i})$ (7)

where $t$ and $y$ is groundtruth and output, respectively.

The objective function $O$ applied in this study is the summation of the
binary cross entropy of each label with the regularization term using graph
based smoothing, which is defined as

$O_{1}=E_{bce}+\lambda R_{gbs}$ (8)

Parameter $\lambda$ is used to control the effect of the regularizer.

$O_{2}=E_{bce}+\lambda R_{glrdn}$ (9)

###  II-B Graph Laplacian based on Differences of Neighboring Pixels (GLRDN)

To define the GLRDN, we consider to utilize the differences of the differences
of neighboring pixels between the target image and estimated image[12]. Let
assume target image ${t}_{m}$ and the estimated images ${y}_{m}$, the GLRDN
define as

$\displaystyle R_{glrdn}(\bm{t}_{m},\bm{y}_{m})$ $\displaystyle=\sum_{(i,j)\in
E}\\{(t_{mi}-t_{mj})-(y_{mi}-y_{mj})\\}^{2}$ $\displaystyle=\sum_{(i,j)\in
E}(\Delta t_{m}^{ij}-\Delta y_{m}^{ij})^{2}$
$\displaystyle=(\Delta\bm{t}_{m}-\Delta\bm{y}_{m})^{T}(\Delta\bm{t}_{m}-\Delta\bm{y}_{m})$
$\displaystyle=(B\bm{t}_{m}-B\bm{y}_{m})^{T}(B\bm{t}_{m}-B\bm{y}_{m})$
$\displaystyle=(\bm{t}_{m}-\bm{y}_{m})^{T}B^{T}B(\bm{t}_{m}-\bm{y}_{m})$
$\displaystyle=(\bm{t}_{m}-\bm{y}_{m})^{T}L(\bm{t}_{m}-\bm{y}_{m})$ (10)

where $B$ is the incident matrix. It is noticed that $B^{T}B$ is equal to the
graph Laplacian matrix $L$. Graph Laplacian matrix assumes that the
differences of neighbor pixel between target images and estimated images,
denoted as $(\bm{t}_{m}-\bm{y}_{m})$ is smooth with respect to concerning
corresponding graph G. In particular, it enforce the value of
$(\bm{t}_{m}-\bm{y}_{m})^{T}L(\bm{t}_{m}-\bm{y}_{m})$ should be small. In the
training process, GLRDN is use with BCE objective function

###  II-C Euler Characteristic-based Regularizer

The Euler characteristic is a property of an object’s shape that is not
affected by topological changes such as scale transformations and rotations.
EC can define the number of objects in two-dimensional space provided there
are no holes in the object. EC can be obtained by calculation of the relations
between vertices (P), edges (S), and faces (F ) :

$EC=F-S+P$ (11)

Figure 5: Comparison results of blood vessel segmentation between our proposed
regularizer and without proposed regularizer on the DRIVE (top row) and STARE
(bottom row) datasets.

For a binary image, the EC can be constructed in the following steps: (a)
define pixels as imaginary vertices, (b) calculate the number of edges based
on 8x8 neighboring pixels, (c) count the number of vertices, (d) count the
number of edges, (e) count the number of faces, (f) calculate the EC using Eq.
(11).

We demonstrate the calculation of Euler Characteristic of simple binary image
as shown in Fig. 4. In this study, we calculated the EC from the CNN output.
Since the EC of an object equals the number of objects in the image, we use
this information to determine the number of isolated objects in the CNN
output. Then, information about the number of isolated objects is combined
with BCE loss as a regularizer in the CNN training process To minimize the
number of isolated objects.

For the EC calculation to be invariant, it is necessary to calculate the
vertices, sides, and faces in the opposite direction. Then, calculated the
average of the EC in two different directions as follows:

$R_{ec}=\frac{{EC}_{1}+{EC}_{2}}{2}$ (12)

In this study, EC estimate from the segmentation output of network[11]. We
combine the objective function with regularization term as:

$O_{3}=E_{bce}+\lambda R_{ec}$ (13)

where $\lambda$ is denoted the parameter to adjust the effect of EC
regularizer.

##  III Experiments

We measure the performance of our proposed method using Sensitivity,
Specificity, Accuracy, and Area Under curve on DRIVE, STARE, and CHASEDB1
datasets. To increase the number of datasets, we divide the image into patches
of 48 x 48 pixels to get 4,750 patches for each image. In this way, we get a
sufficient number of datasets to perform the training step. The first, we
compare our proposed method with baseline U-Net architecture with and without
regularizer term. For the experiment, we train the network with 100 epochs and
initialize the learning rate to 0.001. To achieve convergence, Every time we
reach 25 epochs, we reduce the learning rate by 10 times. For the training
sample, we set the batch to 32 and used Adam as the optimizer.

TABLE I: Comparison of architectural performance without and with our proposed GBS on the DRIVE datasets. | Methods | Sn | Sp | Acc | AUC
---|---|---|---|---|---
| Baseline | 0.6707 | 0.9867 | 0.9465 | 0.9652
| Baseline+GBS | 0.7064 | 0.9897 | 0.9536 | 0.9794
TABLE II: Comparison of architectural performance without and with our proposed GLRDN on the DRIVE datasets. | Methods | Sn | Sp | Acc | AUC
---|---|---|---|---|---
| Baseline | 0.7429 | 0.9840 | 0.9544 | 0.9686
| Baseline+GLRDN | 0.7914 | 0.9791 | 0.9561 | 0.9740
TABLE III: Comparison of architectural performance without and with our proposed Euler Characteristic-based Regularizer on the DRIVE datasets. | Methods | Sn | Sp | Acc | AUC
---|---|---|---|---|---
| Baseline | 0.7583 | 0.9826 | 0.9551 | 0.9691
| Baseline+EC | 0.8463 | 0.9759 | 0.9600 | 0.9824
Figure 6: Receiver Operating Characteristics (ROC) curve analysis of blood
vessel segmentation between our proposed regularizer and without proposed
regularizer on the different datasets. (a)DRIVE (b)STARE and (c)CHASEDB1.

##  IV Results

We show some segmentation results from the experiment in the figure 5. The
figure shows a very complicated retinal vessel consisting of large and small
vessels crossing each other. The segmentation results from the U-Net baseline
appear poor in detecting complex vessels and small vessels. However, the
proposed regularizer can improve the segmentation results on small vessels
with more detail and better.We also present the results of our quantitatively
proposed regularizer on the DRIVE dataset in the I, II, and III tables.
Without a regularizer, the baseline network loses a lot of small vessels and
creates many isolated objects. These tables indicate that our regularizer is
able to segment small vessels and minimize the number of isolated objects on
the vessel than without using a regularizer. Figure 6 shows the performance of
our proposed regularizer compared to the baseline network without regularizer
on different datasets using ROC curve analysis. In the figure, the blue line
shows the ROC curve of the baseline network without the regularizer. The
orange line shows the ROC curve of the baseline network merged with our
regularizer. The two ROC curve lines show that our proposed regularizer can
improve the performance of the baseline network. This indicates that our
proposed performance regularizer can provide better performance enhancement on
retinal blood vessel segmentation. To get better visualization, we zoomed in
on the ROC curve at the top left.

##  V Conclusions

This study takes advantage of the neighboring pixel information in the image
as additional information on the CNN. This information is constructed as a
regularization term in the objective function. The method we propose is able
to improve vessel segmentation, especially on small vessels, better than
without using a regularizer. In addition, isolated objects are also getting
minimal. The calculation of regularizer on a pixel-by-pixel basis is a
limitation of this approach. In the future, we plan to expand the work by
using sparse graphs to reduce computational costs.

## References

  * [1] O. Ronneberger, P. Fischer, and T. Brox, “U-Net: Convolutional Networks for Biomedical Image Segmentation,” in Lecture Notes in Computer Science, Springer International Publishing, 2015, pp. 234–241.
  * [2] V. Badrinarayanan, A. Kendall, and R. Cipolla, “SegNet: A Deep Convolutional Encoder-Decoder Architecture for Image Segmentation,” IEEE Transactions on Pattern Analysis and Machine Intelligence, no. 12, pp. 2481–2495, Dec. 2017, doi: 10.1109/tpami.2016.2644615.
  * [3] H. Fu, Y. Xu, S. Lin, D. W. Kee Wong, and J. Liu, “DeepVessel: Retinal Vessel Segmentation via Deep Learning and Conditional Random Field,” in Medical Image Computing and Computer-Assisted Intervention – MICCAI 2016, Springer International Publishing, 2016, pp. 132–139.
  * [4] M. Havaei et al., “Brain tumor segmentation with Deep Neural Networks,” Medical Image Analysis, pp. 18–31, Jan. 2017, doi: 10.1016/j.media.2016.05.004.
  * [5] D. Q. Zeebaree, H. Haron, A. M. Abdulazeez, and D. A. Zebari, “Machine learning and Region Growing for Breast Cancer Segmentation,” 2019 International Conference on Advanced Science and Engineering (ICOASE), Apr. 2019, doi: 10.1109/icoase.2019.8723832.
  * [6] B. Chougula, A. Tigadi, P. Manage, and S. Kulkarni, “Road segmentation for autonomous vehicle: A review,” 2020 3rd International Conference on Intelligent Sustainable Systems (ICISS), Dec. 2020, doi: 10.1109/iciss49785.2020.9316090.
  * [7] M. Ullah, A. Mohammed, and F. Alaya Cheikh, “PedNet: A Spatio-Temporal Deep Convolutional Neural Network for Pedestrian Segmentation,” Journal of Imaging, no. 9, p. 107, Sep. 2018, doi: 10.3390/jimaging4090107.
  * [8] H. Cholakkal, G. Sun, F. Shahbaz Khan, and L. Shao, “Object Counting and Instance Segmentation With Image-Level Supervision,” 2019 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), Jun. 2019, doi: 10.1109/cvpr.2019.01268.
  * [9] C. Wang, Y. Tang, X. Zou, W. SiTu, and W. Feng, “A robust fruit image segmentation algorithm against varying illumination for vision system of fruit harvesting robot,” Optik, pp. 626–631, Feb. 2017, doi: 10.1016/j.ijleo.2016.11.177.
  * [10] L. Hakim, N. Yudistira, M. Kavitha, and T. Kurita, “U-Net with Graph Based Smoothing Regularizer for Small Vessel Segmentation on Fundus Image,” in Communications in Computer and Information Science, Springer International Publishing, 2019, pp. 515–522.
  * [11] L. Hakim, M. S. Kavitha, N. Yudistira, and T. Kurita, “Regularizer based on Euler characteristic for retinal blood vessel segmentation,” Pattern Recognition Letters, pp. 83–90, Sep. 2021, doi: 10.1016/j.patrec.2021.05.023.
  * [12] L. Hakim, H. Zheng, and T. Kurita, “Improvement for Single Image Super-resolution and Image Segmentation by Graph Laplacian Regularizer Based on Differences of Neighboring Pixels,” International Journal of Intelligent Engineering and Systems, no. 1, Feb. 2022, doi: 10.22266/ijies2022.0228.10.


# Itinerant ferromagnetism in dilute SU(N) Fermi gases

Jordi Pera  Joaquim Casulleras  Jordi Boronat  Departament de Física, Campus
Nord B4-B5, Universitat Politècnica de Catalunya, E-08034 Barcelona, Spain

(August 27, 2024)

###### Abstract

We present exact analytic results for the energy of a SU(N) repulsive Fermi
gas as a function of the spin-channel occupation at second order in the gas
parameter. This is an extension of previous results that now incorporates the
degree of polarization of the system. Therefore, the magnetic properties of
the gas can be obtained, free from numerical uncertainties. For spin 1/2 we
find that second-order corrections change the itinerant ferromagnetic
transition from continuous to first-order. Instead, for spin larger than 1/2
the phase transition is always of first-order type. The transition critical
density reduces when the spin increases, making the phase transition more
accessible to experiments with ultracold dilute Fermi gases. Estimations for
Fermi gases of Yb and Sr with spin 5/2 and 9/2, respectively, are reported.

##  I Introduction

At low density an electron gas is paramagnetic rather than ferromagnetic
because its energy is essentially kinetic. At a same density, the Fermi energy
of a polarized (ferromagnetic) phase is larger than the one of the unpolarized
(paramagnetic) one [1]. However, increasing the density produces an increase
of the potential energy due to the interaction of electrons with different
spin. At some point, this increment of potential energy can exceed the kinetic
energy gap and the gas will become ferromagnetic. This is the well-known
Stoner scenario of itinerant ferromagnetism [2] where a quantum magnetic
transition is predicted in the absence of any crystal pattern and any magnetic
field.

However, the observation of itinerant ferromagnetism in different materials
has been extremely elusive because an increase in the density can produce non-
Fermi liquids [3] or crystallization [4] before the expected ferromagnetic
transition. The high tunability of trapped cold Fermi gases has offered a new
platform to study itinerant ferromagnetism. However, even in this case, the
observation of this transition has turned elusive because the repulsive branch
is metastable with respect to the formation of spin up-spin down dimers. A
first pioneering observation of the ferromagnetic transition [5] was later
revised and concluded that pair formation precludes the achievement of the
relatively large gas parameters required for observing this transition [6].
More recently, it has been claimed that the ferromagnetic state is effectively
obtained after the observation of spin-domain inmiscibility in Fermi 6Li
around a gas parameter $x=k_{F}a_{0}\simeq 1$, with $k_{F}$ the Fermi momentum
and $a_{0}$ the s-wave scattering length [7]. The theoretical study of
repulsive Fermi gases, to determine the gas parameter at which the
ferromagnetic transition appears, has been intense [8, 9, 10, 11, 12, 13, 14,
15]. Itinerant ferromagnetism is predicted to happen around $x\simeq 1$, with
the most reliable results derived from quantum Monte Carlo (QMC) simulations
[9]. At this value of the gas parameter, the Stoner model is not
quantitatively accurate and the use of microscopic approaches as QMC is the
best tool in spite of the sign problem.

The Stoner model is derived using first-order perturbation theory on the gas
parameter, that is it is the Hartree-Fock model for the repulsive Fermi gas
[2]. For particles with spin $1/2$, this model predicts a continuous
transition from a paramagnetic phase to a ferromagnetic one at $x=\pi/2$. The
transition point obtained with QMC is significantly smaller [9], pointing to
the relevance of interactions which Stoner model incorporates only at first
order. The first correction to the Hartree-Fock energy was obtained by Huang,
Lee, and Yang using pseudopotentials for the hard-sphere Fermi gas [16, 17].
The same result was also derived by Galitskii [18], using Green functions, and
Abrikosov and Khalatnikov by means of Landau theory of Fermi liquids [19].
This analytical result for the second-order perturbation theory gives the
total energy of the Fermi gas as a single function of the gas parameter.
However, one cannot use this equation when the number of particles in each
spin is different or, in other words, the polarization does not enter in it as
a variable.

In recent years, the experimental production of SU(N) fermions has renewed the
theoretical interest in their study. Cazalilla et al. [20] discussed that
Fermi gases made of alkaline atoms with two electrons in the external shell,
such as 173Yb, present an SU(N) emergent symmetry. They also claimed that for
$s>1/2$ the ferromagnetic transition had to be of first order because the
mathematical structure of SU(N$>$2) is significantly different than the one of
SU(2). In 2014, Cazalilla and Rey [21] reviewed the progress made with ultra-
cold alkaline-earth Fermi gases. Ref. [22] was one of the first observations
of SU(N) symmetric interactions in alkaline-earth-metal atoms. Interaction
effects in SU(N) Fermi gases as a function of N were studied in Ref. [23]: for
weak interactions it was predicted that inter-particle collisions are enhanced
by N. Collective excitations in SU(N) Fermi gases with tunable spin proved to
be instrumental to investigate collective properties of large spin systems
[24]. On the other hand, Ref. [25] studied the prethermalization of these
systems starting from different initial conditions finding that, under some
conditions, the imbalanced initial state could be stabilized for a certain
time. Recently, Ref. [26] performed a thorough study of the thermodynamics of
deeply degenerate SU(N) Fermi gases using 87Sr for which N can be tuned up to
10. For temperatures above the super-exchange energy, the behavior of the
thermodynamic quantities was found to be universal with respect to N [27].

In the present work, we solve analytically the second-order perturbation
energy including the polarization or, more precisely, the dependence on the
relative concentration of particles with different $z$-spin component. Our
results are derived for a generic spin and thus, we can apply it to SU(N)
fermions. This generalization allows, for instance, the study of dilute Fermi
gases of Ytterbium [28], with spin 5/2, and Strontium [29], with spin 9/2,
that have been already produced in experiments. With the analytic result for
the energy, as a function of the gas parameter and occupation of spin
channels, we perform a Landau analysis of the ferromagnetic transition for any
spin. We find that, for spin 1/2, the phase transition turns to be first-order
with respect to the polarization instead of the continuous character of it
derived in the Stoner model. Interestingly, the critical density for itinerant
ferromagnetism is observed to decrease monotonously when the spin increases,
opening new possibilities for experimental realizations in cold Fermi gases.

##  II Methodology

We study a repulsive Fermi gas at zero temperature with spin $S$ and spin
degeneracy $\nu=2S+1$. In the dilute gas regime, only particles with different
$z$-spin component interact via a central potential $V(r)$ ($s$-wave
scattering). The number of particles in each spin channel is
$N_{\lambda}=C_{\lambda}N/\nu$, with $N$ the total number of particles and
$C_{\lambda}$ being the fraction of $\lambda$ particles (normalized to be one
if the system is unpolarized, $N_{\lambda}=N/\nu$, $\forall\lambda$). The
Fermi momentum of each species is $k_{F,\lambda}=k_{F}C_{\lambda}^{1/3}$, with
$k_{F}=(6\pi^{2}n/\nu)^{1/3}$. The kinetic energy is readily obtained,

$\frac{T}{N}=\frac{3}{5}\epsilon_{F}\frac{1}{\nu}\sum_{\lambda}C_{\lambda}^{5/3},$
(1)

with $\epsilon_{F}=\hbar^{2}k_{F}^{2}/(2m)$ the Fermi energy. The lowest-order
contributions to the potential energy are given by [30]

$\displaystyle
V=\frac{\hbar^{2}\Omega}{2m}\sum_{\lambda_{1},\lambda_{2}}\int\frac{d\textbf{l}}{(2\pi)^{3}}n_{l}\int\frac{d\textbf{k}}{(2\pi)^{3}}n_{k}$
(2)
$\displaystyle\times\\{K(\textbf{k},\textbf{l};\textbf{k},\textbf{l})-\delta_{\lambda_{1},\lambda_{2}}K(\textbf{k},\textbf{l};\textbf{l},\textbf{k})\\}\
,$

with $\Omega$ the volume, $n_{l}$ and $n_{k}$ the momentum distributions of
the free Fermi gas. Up to second order in the $s$-wave scattering length
$a_{0}$, the scattering $K$ matrix is given by [30, 31]

$K(r,R)=4\pi a_{0}+(4\pi a_{0})^{2}\,I(r,R)+O(a_{0}^{3})\ ,$ (3)

with $\textbf{r}=(\textbf{k}-\textbf{l})/2$ and
$\textbf{R}=\textbf{k}+\textbf{l}$, the relative and total momentum in the
center of mass frame, respectively. The function $I(r,R)$ is defined as

$\displaystyle I(r,R)=\frac{1}{(2\pi)^{3}}$ $\displaystyle\times\int
2\,d\textbf{q}d\textbf{q}^{\prime}\frac{1-(1-n_{q})(1-n_{q^{\prime}})}{q^{2}+q^{\prime
2}-k^{2}-l^{2}}\delta(\textbf{q}+\textbf{q}^{\prime}-\textbf{k}-\textbf{l})\
.$

Considering only the first term in the expansion of the $K$-matrix (3), one
gets for the potential energy (2) the well-know Hartree-Fock energy [2],

$\left(\frac{V}{N}\right)_{1}=\frac{2\epsilon_{F}}{3\pi}\left[\frac{1}{\nu}\sum_{\lambda_{1},\lambda_{2}}C_{\lambda_{1}}C_{\lambda_{2}}(1-\delta_{\lambda_{1},\lambda_{2}})\right]\,x\
,$ (5)

with $x\equiv k_{F}a_{0}$ the gas parameter of the Fermi gas.

The second order term in the gas parameter $x$ is due to the second term of
the $K$-matrix. This second order term reads

$\left(\frac{V}{N}\right)_{2}=\frac{\epsilon_{F}}{k_{F}^{7}}\left[\frac{1}{\nu}\sum_{\lambda_{1},\lambda_{2}}I_{2}(k_{F,\lambda_{1}},k_{F,\lambda_{2}})(1-\delta_{\lambda_{1},\lambda_{2}})\right]x^{2}\
,$ (6)

with

$\displaystyle I_{2}(C_{\lambda_{1}},C_{\lambda_{2}})=\frac{3}{16\pi^{5}}\int
d\textbf{l}\,n_{l}\int d\textbf{k}\,n_{k}\int
2\,d\textbf{q}d\textbf{q}^{\prime}$ (7)
$\displaystyle\times\frac{1-(1-n_{q})(1-n_{q^{\prime}})}{q^{2}+q^{\prime
2}-k^{2}-l^{2}}\delta(\textbf{q}+\textbf{q}^{\prime}-\textbf{k}-\textbf{l})\
.$

The calculation of $I_{2}(C_{\lambda_{1}},C_{\lambda_{2}})$ (7) is rather
involved. For the unpolarized phase, i.e., when all the spin states are
equally populated, the integral in Eq. (7) was made in the fifties of the past
century [17, 16, 18, 19]. On the other hand, when the gas has a finite
polarization that integral becomes more cumbersome. In previous work, it was
solved partially but with a final numerical integration [32]. We have been
able to integrate Eq. (7) and found an analytical expression for it (See App.
C). Our result is the following,

$I_{2}(C_{\lambda_{1}},C_{\lambda_{2}})=\frac{4k_{F}^{7}}{35\pi^{2}}C_{\lambda_{1}}C_{\lambda_{2}}\frac{C_{\lambda_{1}}^{1/3}+C_{\lambda_{2}}^{1/3}}{2}\,F(y)\
,$ (8)

with

$F(y)=\frac{1}{4}\big{(}15y^{2}-19y+52-19y^{-1}+15y^{-2}\big{)}+\frac{7}{8}y^{-2}\big{(}y-1\big{)}^{4}\big{(}y+3+y^{-1}\big{)}\ln{\bigg{|}\frac{1-y}{1+y}\bigg{|}}-\frac{2y^{4}}{1+y}\ln{\bigg{|}1+\frac{1}{y}\bigg{|}}-\frac{2y^{-4}}{1+y^{-1}}\ln{\bigg{|}1+y\bigg{|}}$
(9)

and $y\equiv(C_{\lambda_{1}}/C_{\lambda_{2}})^{1/3}$. Our result reproduces
the formula derived by Kanno [33] for the specific case of $s=1/2$ and hard
spheres.

Assembling it all together, we can write the energy per particle up to second
order in $x$, and for any occupation of the $\nu$ available spin states, as

$\displaystyle\frac{E}{N}$ $\displaystyle=$
$\displaystyle\frac{3\epsilon_{F}}{5\nu}\left\\{\sum_{\lambda}C_{\lambda}^{5/3}+\frac{10}{9\pi}\left[\sum_{\lambda_{1},\lambda_{2}}C_{\lambda_{1}}C_{\lambda_{2}}(1-\delta_{\lambda_{1},\lambda_{2}})\right]\,x\right.$
(10) $\displaystyle+$
$\displaystyle\left.\frac{5}{3k_{F}^{7}}\left[\sum_{\lambda_{1},\lambda_{2}}I_{2}(C_{\lambda_{1}},C_{\lambda_{2}})(1-\delta_{\lambda_{1},\lambda_{2}})\right]\,x^{2}\right\\}\
,$

with $I_{2}(C_{\lambda_{1}},C_{\lambda_{2}})$ given by Eqs. (8,9). Eq. (10) is
a perturbative expansion in $x$ and works fine for low values of $x$, $x<1$.
If the gas is unpolarized, i.e., $C_{\lambda}=1,\,\forall\lambda$, the energy
per particle reduces to the known expression [30],

$\frac{E}{N}=\frac{3\epsilon_{F}}{5}\left\\{1+(\nu-1)\left[\frac{10}{9\pi}x+\frac{20}{105\pi^{2}}(11-2\ln
2)\,x^{2}\right]\right\\}\ .$ (11)

The energy of the interacting Fermi gas (10) is written in terms of the
occupation of the different spin channels $C_{\lambda}$. This set of values
results in a polarization $P$ of the system, in such a way that when all the
spin states are equally populated the gas is unpolarized $P=0$ and, if only
one of them is populated, $|P|=1$. Keeping the total number of particles $N$
as constant, and for $s=1/2$, there is only one ratio of spin occupations for
a given value of $P$. In contrast, for $s>1/2$ there are more combinations. It
can be shown that the solution which optimizes the energy is the one in which
the increase of particles in one spin state comes for an equal decrease of the
rest with constant $N$ (See App. A). Under these conditions, the
concentrations $C_{\lambda}$ for a given polarization $P$ are

$\displaystyle C_{+}$ $\displaystyle=$ $\displaystyle 1+|P|\,(\nu-1)$ (12)
$\displaystyle C_{\lambda\neq+}$ $\displaystyle=$ $\displaystyle 1-|P|\ ,$
(13)

with subindex $+$ standing for the state with the larger population. As we
have discussed, for $s>1/2$ there are more possible configurations. For
example, Ref. [25] works with a system of $s=3/2$ with population imbalance:
the number of atoms with $s=3/2$, $N_{\pm 3/2}$, is different than the one
with $s=1/2$, $N_{\pm 1/2}$. As we need to make a choice, we choose the one
that minimizes the energy.

It is interesting to check if Eq. (10) converges when the degeneracy
increases. If we set the limit $\nu\rightarrow\infty$ in Eq. (10), we get

$\displaystyle\frac{E}{N}=\frac{2\pi\hbar^{2}}{ma_{0}^{2}}\bigg{[}\frac{3}{20\pi}(6\pi^{2}n)^{2/3}a_{0}^{2}P^{5/3}+na_{0}^{3}(1-P^{2})$
(14)
$\displaystyle+\frac{1}{2\pi^{3}}(6\pi^{2}n)^{4/3}a_{0}^{4}(1-P)P^{4/3}\bigg{]}\
.$

As one can see, Eq. (14) is convergent, as it is finite for any value of $P$.
Particularizing to the unpolarized case $P=0$, one obtains the Hartree-Fock
energy for bosons. And, if we set $P=1$, we obtain the energy of a non-
interacting Fermi gas.

##  III Landau Theory

Upon an increase of the gas parameter, the interacting Fermi gas will
eventually become ferromagnetic. To localize the transition point and the
order of the phase transition we rely on Landau theory. We first explore what
we obtain for the Stoner model, that is, first order approximation in $x$.
Then, we consider the second-order expansion, Eq. (10). In this way, we can
compare the differences that arise due to increasing the order of the
perturbative expansion. At the Hartree-Fock level, first order in $x$, the
Landau expansion of the energy is given by

$f-f_{0}(x)=-\frac{1}{2}A(x-x_{0})P^{2}-\frac{1}{3}B{|P|}^{3}+\frac{1}{4}CP^{4}\
,$ (15)

with $f=5E/(3N\epsilon_{F})$, $x_{0}=\pi/2$, and the rest of constants are
given in App. B. In order to determine the transition point, one imposes two
conditions: i) the energy has to be a minimum, that is, its first derivative
must be zero, and ii) the energy must be smaller than the zero-polarization
energy to have a global minimum. With these two criteria, one finds
$x^{\ast}=x_{0}-2B^{2}/(9AC)$ and $|P^{\ast}|=2B/(3C)$. The constant $B$ is
proportional to $(\nu-2)$, therefore the Stoner model predicts a first-order
phase transition for spin higher than 1/2, while a continuous one for $s=1/2$.

Equipped with the second-order expression of the energy as a function of the
spin-state occupations (10), we can make a Landau expansion of the energy. We
get

$\displaystyle
f-f_{0}(x)=-\frac{1}{2}A(\overline{x}(x)-x_{0})P^{2}-\frac{B(x)}{3}{|P|}^{3}$
(16) $\displaystyle+\frac{C(x)}{4}P^{4}+\frac{L(x)}{4}P^{4}\ln{|P|}\ .$

The coefficients in Eq. (16) depend now on the gas parameter $x$ and have a
more complex expression (See App. B). For the sake of simplicity, from now on,
we will not specify that the coefficients of Eq. (16) depend on the gas
parameter $x$. Applying the same criteria for finding the transition point as
before, we derive the equations for the jump of polarization and transition
density values,

$\ln{|P^{\ast}|}=-\bigg{(}\frac{1}{2}+\frac{C}{L}\bigg{)}+\frac{2B}{3L|P^{\ast}|}$
(17)
$\overline{x}^{\ast}=x_{0}-\frac{B}{3A}|P^{\ast}|-\frac{L}{4A}|P^{\ast}|^{2}\
.$ (18)

We rearrange Eq. (17) to explore the range where a solution exists,

$\frac{2B}{3L|P^{\ast}|}-\ln{|P^{\ast}|}=\frac{1}{2}+\frac{C}{L}\ .$ (19)

The function $2B/(3L|P^{\ast}|)-\ln{|P^{\ast}|}$ is continuously decreasing
between $|P^{\ast}|=0$ and $|P^{\ast}|=1$ as long as $2B/(3L)$ is positive. As
this is the case, in this range, the minimum value will be at $|P^{\ast}|=1$.
Hence, there will exist a solution when the coefficients satisfy the condition

$\frac{1}{2}+\frac{C}{L}-\frac{2B}{3L}\geq 0\ .$ (20)

This condition is fulfilled for any spin value and for densities lower than
the transition one. For spin 1/2, the above equations suffer an important
simplification as the coefficient $B$ is zero,

$\displaystyle|P^{\ast}|=\exp{\bigg{(}-\frac{1}{2}-\frac{C}{L}\bigg{)}}$ (21)
$\displaystyle\overline{x}^{\ast}=x_{0}-\frac{L}{4A}\exp{\bigg{(}-1-\frac{2C}{L}\bigg{)}}\
.$

It is worth noticing the relevance of the coefficient $L$ in Eq. (21), coming
from the term $P^{4}\ln|P|$, for distinguishing the order of the phase
transition: $L$ needs to be different from zero to get a first-order phase
transition [34]. Introducing the values of the constants in Eq. (21), one gets
the transition density at $x=\pi/3(1+0.007)\approx 1.0545$, and a jump in the
polarization equal to $0.545$. Remarkably, the second-order energy changes the
character of the ferromagnetic phase transition becoming now a first-order
one, as it happens for spin $s>1/2$.

##  IV Results

We have applied the above formalism to study fundamental properties of the
Fermi gas, such as the energy or the magnetic susceptibility.

In Fig.1, we show the dependence of the polarization $P$ with the gas
parameter $x$ for a $s=1/2$ Fermi gas. One can observe as the introduction of
second-order corrections to the potential energy modify the character of the
phase transition with respect to the Stoner model. Moreover, the gas parameter
at which the transition occurs is significantly reduced, from $\pi/2$ to
$\sim\pi/3$, approaching the value $x\simeq 1$ where the transition has been
indirectly observed [7].

Figure 1: Itinerant ferromagnetic transition for $s=1/2$. The dashed orange
and solid blue lines stand for the Stoner and second-order approximations,
respectively.

In Fig. 2, we plot the energy of the Fermi gas using the the Stoner one and
the second-order approximation. At very low densities, both models predict the
same energy, which follow the non-polarized behavior. After a certain value of
the gas parameter, the lines become flat, which is the behavior of the fully-
polarized gas. We can see that, when we use the second-order approximation,
the transition happens at a smaller gas parameter.

Figure 2: Energy per particle as a function of the gas parameter for $s=1/2$.
The dashed orange and solid blue lines stand for the Stoner and second-order
approximations, respectively. Figure 3: Energy per particle as a function of
the polarization $P$ for $s=1/2$. The lines correspond to different $x$ values
close to the phase transition $x^{\star}=1.054$. The crosses in each line
indicates the polarization where the energy is minimum.

The first-order phase transition for $s=1/2$ is also observed in Fig. 3.
There, we plot the energy as a function of $P$ and for gas parameter values
close to the phase transition. One can see that the minimum of the energy
jumps from $P=0$ to an intermediate value $P=0.545$ that then progressively
moves to the fully polarized phase (Fig. 1).

Since we are interested in the magnetism of these gases, we proceed to analyze
the magnetic susceptibility, which is inversely proportional to the second
derivative of the energy with respect to the polarization,

$\frac{1}{\chi}=\frac{1}{n}\left(\frac{\partial^{2}(E/N)}{\partial
P^{2}}\right)_{x}\ .$ (22)

If we split the total energy between the kinetic and the potential energy, we
can rewrite $\chi$ as

$\chi=\frac{3}{2}\frac{n}{\epsilon_{F}}\left[\frac{1}{\nu}\sum_{\lambda}C_{\lambda}^{-1/3}(C_{\lambda}^{\prime})^{2}+\frac{3}{2}\frac{1}{\epsilon_{F}}(V/N)^{\prime\prime}\right]^{-1}\
,$ (23)

where the derivatives are with respect to the polarization. In Fig. 4, we show
the magnetic susceptibility around the ferromagnetic transition point for
$s=1/2$. We see again that the transition occurs at a lower value of $x$ [9]
with respect to the Stoner model. One can also notice that $\chi$ changes
behavior, from diverging at the transition point (Stoner) to a large but
finite peak at second order, reflecting the change in the type of phase
transition from a continuous to a first-order one.

Figure 4: Magnetic susceptibility $\chi$ around the transition point for
$s=1/2$. The peak of $\chi$ appears at a lower $x$ value than in the Stoner
model.

From Eq. (23), one can prove that at first order $\chi$ diverges. At this
order, $\chi$ behaves around the transition as

$\chi=\frac{3}{2}\frac{n}{\epsilon_{F}}\left[-\frac{9}{10}A(x-x_{0})\right]^{-1}=\frac{3}{2}\frac{n}{\epsilon_{F}}\frac{5C}{B^{2}}\
.$ (24)

The coefficient $B$ in Eq. (24) is proportional to $(\nu-2)$ and $C$ is finite
for $\nu=2$ (See App. B), hence, the magnetic susceptibility diverges around
the transition for the Stoner model.

We can extend our analysis to larger spin Fermi gases. In this case, the
ferromagnetic transition is first-order as in the Stoner model but the
transition point is, in all cases, observed at smaller $x$ values. This fact
can be seen in Fig. 5 where the energy of the spin 9/2 Fermi gas changes its
behavior at a lower gas parameter than the energy of the spin 5/2 Fermi gas.

Figure 5: Energy per particle as a function of the gas parameter. The solid
blue and dashed orange lines stand for a Fermi gas of $s=5/2$ and $s=9/2$,
respectively. Both lines have been calculated within the second-order
approximation. Figure 6: Energy of the $s=5/2$ (top) and $s=9/2$ (bottom)
Fermi gases as a function of $P$ for $x$ values close to the transition. The
crosses indicate the values of $P$ where the energy is minimum.

In Fig. 6, we show the energies of SU(N) Fermi gases up to second order as a
function of the polarization for $x$ values close to the phase transition. The
top panel is for $s=5/2$ and the bottom one for $s=9/2$ corresponding to Yb
and Sr, respectively. In both cases, the location of the minimum of the energy
as a function of the polarization jumps from $P=0$ to 1 without intermediate
values, in contrast to the case of $s=1/2$. For the sake of completeness, in
Fig. 7 we show the magnetic susceptibility for the two high degenerate Fermi
gases considered: $s=5/2$ and $s=9/2$.

Figure 7: Magnetic susceptibility per particle as a function of the gas
parameter. The slid blue and dashed orange lines stand for a Fermi gas of
$s=5/2$ and $s=9/2$, respectively. Both lines have been computed with the
second-order approximation.

Thanks to having an analytical expression for the energy, we have access to
other important properties as the Tan’s contact. According to Tan relations,
the microscopic behavior of the wave function of the $N$-body system at short
distances ($r\ll n^{-1/3}$) is connected with the behavior of several
macroscopic magnitudes. In particular, the tail of the momentum distribution
for large k values as $k^{-4}$ and the dependence of the energy on the
scattering length. Using the adiabatic sweep theorem [35], the Tan’s contact
can be obtained from the energy of the system by

$C=\frac{8\pi
ma_{0}^{2}}{\nu\hbar^{2}}\frac{N}{V}\frac{\partial(E/N)}{\partial a_{0}}\ .$
(25)

In terms of the gas parameter,

$C=\frac{4\pi
nk_{F}}{\nu}\frac{x^{2}}{\epsilon_{F}}\frac{\partial(E/N)}{\partial x}\ .$
(26)

We have calculated the Tan’s contact using the energy of the Fermi gas at
second order. In Fig. 8, we show $\nu C$ as a function of $x$ for spin
$s=1/2$, $5/2$, and $9/2$. In the paramagnetic phase, $\nu C$ increases
monotonically until it reaches its maximum value at the transition point.
After crossing the ferromagnetic transition, the Tan’s contact becomes zero
because the energy of the fully polarized phase does not depend on the
scattering length $a_{0}$. As one can see in the figure, the $\nu C$ increases
with the value of the spin for a given $x$ value, in agreement with the
increase of interaction energy with spin degeneracy. It is interesting to
notice that, for spin $1/2$, $C$ does not drop abruptly to zero because there
are stable polarizations between 0 and 1. In contrast, for $s>1/2$ the drop is
directly to zero because of the sudden change of the polarization from 0 to 1
at the phase transition.

Figure 8: Tan’s contact in SU(N) Fermi gases. $\nu C$ as a function of the gas
parameter for three different spins: $s=1/2$, $5/2$, and $9/2$. We use the
second-order model for the three cases.

In Fig. 9, we show the evolution of the Tan’s contact as a function of the
spin for the Stoner and second-order models. By increasing the spin, one can
see a tendency to reach a plateau. This plateau, which has been interpreted as
the Bose limit, can be understood by looking at the Tan’s contact behavior
before the transition. For the Stoner model, it is given by

$C=4\pi nk_{F}x^{2}\frac{2}{3\pi}\bigg{(}1-\frac{1}{\nu}\bigg{)}\ ,$ (27)

and in second-order,

$C=4\pi
nk_{F}x^{2}\bigg{[}\frac{2}{3\pi}+\frac{8}{35\pi^{2}}(11-2\ln{2})x\bigg{]}\bigg{(}1-\frac{1}{\nu}\bigg{)}\
.$ (28)

The dependence with $\nu$ for both models is $1-1/\nu$, hence, if we set the
limit $\nu\rightarrow\infty$, we obtain a plateau. In Fig. 9, we also plot
experimental results from Ref. [36] that show a similar behavior to our
theoretical expressions (Eqs. (27) and (28)). Moreover, they lie between both
models, but we cannot state which one fits better due to the experimental
uncertainty. As the experimental points have a different normalization than
our definition, we have to scale them by a factor of $(2\pi)^{2}$.

Figure 9: Tan’s contact in SU(N) Fermi gases. Evolution of $C$ with the spin
of particles at fixed density $x=0.3$ for the Stoner model and the second
order model; the points with error bars are experimental data from Ref. [36]
re-scaled by a factor of $(2\pi)^{2}$. The green and yellow lines follow the
behavior of the Tan’s contact of Eq. (27) and Eq. (28) respectively.

The critical value of the gas parameter depends on the spin of the particles.
In Fig. 10, we report the results obtained for both models (Stoner and second
order). As one can see, the itinerant ferromagnetic transition happens at $x$
values that decrease monotonically with the spin degeneracy. Moreover, the
second order values of $x$ are lower than the ones predicted by the Stoner
model. As for $s>1/2$ the polarization goes from 0 to 1 without intermediate
values ($P^{*}=1$), we can find more easily the laws that the critical values
$x^{*}$ follow. For the Stoner model, the law is given by

Figure 10: Critical gas parameter values of the ferromagnetic transition as a
function of the spin for the Stoner model and the second order model. The
green and yellow lines follow the behavior predicted by Eqs. (29) and (30).

$x^{*}=\frac{9\pi}{10(\nu-1)}\big{(}\nu^{2/3}-1\big{)}\ ,$ (29)

and in second-order,

$x^{*}=\pi\frac{-35+35\sqrt{1+\cfrac{108(11-2\ln{2})}{175(\nu-1)}\big{(}\nu^{2/3}-1\big{)}}}{12(11-2\ln{2})}\
.$ (30)

Notice that for $s=1/2$ Eqs. (29) and (30) do not hold. In this case, the
value of $x^{*}$ is slightly smaller than the value predicted by the laws
(29,30) because the transition is from $P=0$ to $P<1$ (see Fig. 10). The
behavior of the critical gas parameter values with the spin, in Fig. 10, can
be understood by the balance between Fermi kinetic energy and interaction
energy. Increasing the spin degeneracy originates in turn an increase in the
number of interacting pairs of particles because the interatomic potential
acts only between pairs of different $z$-spin component. In other words, the
Fermi gas becomes more interacting and reaches the fully polarized
(ferromagnetic) phase at lower $x$ values.

##  V Conclusions

Summarizing, we present the analytic expression of the energy of a repulsive
SU(N) Fermi gas, in terms of the spin-channel occupations, at second order of
the gas parameter. This analytic derivation allows for an accurate estimation
of the magnetic properties of the Fermi gas for any value of the spin.
Moreover, by using the analytical solution one directly avoids any uncertainty
coming from the numerical integration. This is in fact quite important since
the function to integrate in the second-order term has many singular points
and, if one does not apply a previous analytical treatment to the function,
the numerical method may simply diverge. In order to study the system, we have
chosen the occupational configuration and the polarization protocol that
minimizes the energy: at $P=0$, all the species are equally occupied; and,
when we increase $P$, one species increases and all the rest diminish in the
same manner. However, we point out that our formalism can be applied to any
occupational configuration and any polarization protocol. One just needs to
find the new expressions for the fractions of $\lambda$ particles
$C_{\lambda}$. In fact, in many experiments, the species’ occupation can be
tuned or controlled almost at will [24, 26] and there are studies that have
dealt with imbalanced systems [25]. However, we point out that all these other
configurations are excited states, they have more energy than the
configuration we have chosen. Due to that, if a Fermi system can thermalise,
it will fall to the behavior we have predicted in this work. Using another
occupational configuration and a different polarization protocol leads to
different critical values of the gas parameter $x^{*}$.

Our results show that the ferromagnetic transition turns out to be first-order
for $s=1/2$, in contrast with the continuous transition obtained at the
Hartree-Fock (Stoner) approximation [2]. At second order, and for $s>1/2$, the
phase transition is always of first-order type. Our derivation applies to any
spin of the Fermi particles and this allows for the study of itinerant
ferromagnetism in SU(N) fermions. The critical gas parameter for spin 1/2
decreases significantly with respect to the Stoner model, approaching the
experimental estimation $x\simeq 1$ [7]. Remarkably, the critical $x$
decreases monotonically with the spin of the particles suggesting that the
observation of itinerant ferromagnetism could be favored by working with
highly-degenerate gases as Yb [28] and Sr [29]. In order to verify
experimentally our predictions, one could use a similar method to the one used
in Ref. [7] where the Fermi system can thermalise, but generalized to larger
spin degeneracies.

Beyond the second-order terms analyzed in this work, the energy ceases to be
universal in terms of the gas parameters because other scattering parameters
of the interactions, mainly the $s$-wave effective range and $p$-wave
scattering length get involved [30]. However, these corrections are expected
to be small due to the diluteness of the Fermi gases in experiments that make
scattering be dominated by $s$-wave scattering.

##  VI Acknowledgments

We acknowledge financial support from MCIN/AEI/10.13039/501100011033 (Spain)
grant No. PID2020-113565GB-C21 and from Secretaria d’Universitats i Recerca
del Departament d’Empresa i Coneixement de la Generalitat de Catalunya, co-
funded by the European Union Regional Development Fund within the ERDF
Operational Program of Catalunya (project QuantumCat, ref. 001-P-001644).

##  Appendix A Polarization in SU(N) Fermi gases

We have considered that when the occupation of one of the spin channels grows,
the rest decreases in equal form, keeping the total number of particles as
constant. This results from the fact that the Hamiltonian does not depend on
the angular momentum of each species, hence, all the species are equivalent.
Moreover, this fact can be proved by minimizing the energy,

$E=\frac{3}{5}N\epsilon_{F}\frac{1}{\nu}\sum_{\lambda}C_{\lambda}^{5/3}+V(C_{s},...,C_{-s})\
,$ (31)

with $\sum_{\lambda}C_{\lambda}=\nu$. We define the channel that grows,

$C_{+}=\nu-\sum_{\lambda\neq+}C_{\lambda}\ ,$ (32)

and substitute it in the energy,

$\displaystyle
E=\frac{3}{5}N\epsilon_{F}\frac{1}{\nu}\bigg{[}\big{(}\nu-\sum_{\lambda\neq+}C_{\lambda}\big{)}^{5/3}+\sum_{\lambda\neq+}C_{\lambda}^{5/3}\bigg{]}+V(\nu-\sum_{\lambda\neq+}C_{\lambda},C_{s-1}...,C_{-s})\
.$ (33)

Now, we minimize the energy for all the channels, except for the one that we
have isolated which, due to the constraint (32), depends on the others,

$\frac{\partial E}{\partial C_{\sigma}}=0\quad\forall\sigma\neq+\ .$ (34)

Then,

$\displaystyle\frac{\partial E}{\partial
C_{\sigma}}=N\epsilon_{F}\frac{1}{\nu}\bigg{[}C_{\sigma}^{2/3}-\big{(}\nu-\sum_{\lambda\neq+}C_{\lambda}\big{)}^{2/3}\bigg{]}+\frac{\partial
V}{\partial C_{\sigma}}-\frac{\partial V}{\partial C_{+}}=0\ .$ (35)

Rearranging terms, we obtain

$N\epsilon_{F}\frac{1}{\nu}C_{\sigma}^{2/3}+\frac{\partial V}{\partial
C_{\sigma}}=N\epsilon_{F}\frac{1}{\nu}\big{(}\nu-\sum_{\lambda\neq+}C_{\lambda}\big{)}^{2/3}+\frac{\partial
V}{\partial C_{+}}\ .$ (36)

From Eq. (36), one can see that all the decreasing $C_{\sigma}$ satisfy the
same equation, therefore, they must be equal.

##  Appendix B Parameters of the Landau expansion

In this section we show the expression for all the parameters of the Landau
expansion in terms of the gas degeneracy ($\nu=2s+1$) and the gas parameter
($x=k_{F}a_{0}$). We first show the Landau expansion up to fourth order in $P$
including the logarithmic term that comes from the second order term in
perturbation theory.

$\displaystyle
f-f_{0}=-\frac{A}{2}(\overline{x}-x_{0})P^{2}-\frac{B}{3}{|P|}^{3}+\frac{C}{4}P^{4}+\frac{L}{4}P^{4}\ln{|P|}\
.$ (37)

The parameters that are invariant in any order are

$\displaystyle f=\frac{5E}{3N\epsilon_{F}}\quad;\quad
A=\frac{20}{9\pi}(\nu-1)\quad;\quad x_{0}=\pi/2\ .$ (38)

We split this section into two subsections. The first one contains the
expression of the parameters up to first order in perturbation theory
(Hartree-Fock), and the second one up to second order.

###  B.1 1st order parameters

$\displaystyle
f_{0}=1+\frac{10}{9\pi}(\nu-1)x\quad;\quad\overline{x}=x\quad;\quad
B=\frac{5}{27}(\nu-1)(\nu-2)\quad;\quad
C=\frac{20}{243}(\nu-1)(\nu^{2}-3\nu+3)\quad;\quad L=0$ (39)

###  B.2 2nd order parameters

$\displaystyle
f_{0}=1+\frac{10}{9\pi}(\nu-1)x+\frac{4}{21\pi^{2}}(\nu-1)(11-2\ln{2})x^{2}\quad;\quad\overline{x}=x+\frac{2}{15\pi}(22-7\nu+4(\nu-1)\ln{2})x^{2}$
(40) $\displaystyle
B=\frac{5}{27}(\nu-1)(\nu-2)+\frac{2}{81\pi^{2}}(\nu-1)(\nu-2)(44+\nu+8(\nu-1)\ln{2})x^{2}$
$\displaystyle C=\frac{20}{243}(\nu-1)(\nu^{2}-3\nu+3)+$
$\displaystyle+(\nu-1)\frac{528-336\nu-16\nu^{2}+29\nu^{3}-(96-192\nu+128\nu^{2}-32\nu^{3})\ln{2}}{729\pi^{2}}x^{2}+\frac{20\nu^{3}}{243\pi^{2}}(\nu-1)\ln{\frac{\nu}{6}}x^{2}$
$\displaystyle L=\frac{20\nu^{3}}{243\pi^{2}}(\nu-1)x^{2}$

For spin 1/2, these expressions reduce to

$\displaystyle
f_{0}=1+\frac{10}{9\pi}x+\frac{4(\nu-1)(11-2\ln{2})}{21\pi^{2}}x^{2}\quad;\quad\overline{x}=x+\frac{8(2+\ln{2})}{15\pi}x^{2}\quad;\quad
B=0$ (41) $\displaystyle
C=\frac{20}{243}+\frac{8(3+4\ln{2})}{729\pi^{2}}x^{2}-\frac{160\ln{3}}{243\pi^{2}}x^{2}\quad;\quad
L=\frac{160}{243\pi^{2}}x^{2}\ .$

##  Appendix C Calculation of the second-order term

In the following, we detail the calculation of the second-order term for the
energy,

$\displaystyle\frac{E}{N}=\epsilon_{F}\frac{1}{\nu
k_{F}^{5}}\sum_{\lambda_{1},\lambda_{2}}I_{2}(k_{F,\lambda_{1}},k_{F,\lambda_{2}})a_{0}^{2}(1-\delta_{\lambda_{1},\lambda_{2}});$
(42) $\displaystyle\mbox{where}\quad
I_{2}(k_{F,\lambda_{1}},k_{F,\lambda_{2}})=\frac{3}{16\pi^{5}}\int
d\textbf{l}n_{l}\int d\textbf{k}n_{k}\int
2d\textbf{q}d\textbf{q}^{\prime}\frac{1-(1-n_{q})(1-n_{q^{\prime}})}{q^{2}+q^{\prime
2}-k^{2}-l^{2}}\delta(\textbf{q}+\textbf{q}^{\prime}-\textbf{k}-\textbf{l})\
.$

We need to calculate the multiple integral $I_{2}$. The first step is to
expand the inner part, where the occupation functions appear. One, then,
obtains an expression with $n_{q}+n_{q^{\prime}}-n_{q}n_{q^{\prime}}$. The
terms containing $n_{q}$ and $n_{q^{\prime}}$ correspond to integrals of
spheres (SI) with respect to q or $\textbf{q}^{\prime}$, which can be
integrated. The rest of the integrals with respect to k and l are quite
arduous to integrate, but after some lengthy calculations, one can obtain
them. And finally, the term proportional to $n_{q}n_{q^{\prime}}$,
corresponding to the volume of the intersection of two spheres (SII) is zero
due to symmetry reasons.

We rewrite the expression of $I_{2}$ using the integrated functions SI and
SII,

$I_{2}(k_{F,\lambda_{1}},k_{F,\lambda_{2}})=\frac{3}{16\pi^{5}}\int
d\textbf{l}n_{l}\int d\textbf{k}n_{k}\int
2d\textbf{q}d\textbf{q}^{\prime}\frac{n_{q}+n_{q^{\prime}}-n_{q}n_{q^{\prime}}}{q^{2}+q^{\prime
2}-k^{2}-l^{2}}\delta(\textbf{q}+\textbf{q}^{\prime}-\textbf{k}-\textbf{l})$
$=\frac{3}{16\pi^{5}}\int d\textbf{l}n_{l}\int
d\textbf{k}n_{k}\bigg{[}SI(P,R,k_{F,\lambda_{1}})+SI(P,R,k_{F,\lambda_{2}})-SII(P,R,k_{F,\lambda_{1}},k_{F,\lambda_{2}})\bigg{]}$
(43)

First of all, let’s show that the SII term is zero. In order to follow the
derivation, we need to have in mind that the momenta k and q run over the same
momentum $k_{F,\lambda_{1}}$, and that the momenta l and $\textbf{q}^{\prime}$
do so over $k_{F,\lambda_{2}}$. The procedure is the following. We split the
integral in identical parts. The first part is integrated with respect to l,
the second one with respect to $\textbf{q}^{\prime}$, both running over the
same values. We note that it could have been done integrating with respect to
k and q, instead of l and $\textbf{q}^{\prime}$. Then, we slightly manipulate
the two expressions and we obtain two identical integrals but with opposite
signs, hence, they cancel each other,

$\displaystyle\int d\textbf{l}n_{l}\int d\textbf{k}n_{k}\int
2d\textbf{q}d\textbf{q}^{\prime}\frac{n_{q}n_{q^{\prime}}}{q^{2}+q^{\prime
2}-k^{2}-l^{2}}\delta(\textbf{q}+\textbf{q}^{\prime}-\textbf{k}-\textbf{l})=2\int
d\textbf{l}d\textbf{k}d\textbf{q}d\textbf{q}^{\prime}\frac{n_{l}n_{k}n_{q}n_{q^{\prime}}}{q^{2}+q^{\prime
2}-k^{2}-l^{2}}\delta(\textbf{q}+\textbf{q}^{\prime}-\textbf{k}-\textbf{l})$
(44) $\displaystyle=\int
d\textbf{l}d\textbf{k}d\textbf{q}d\textbf{q}^{\prime}\frac{n_{l}n_{k}n_{q}n_{q^{\prime}}}{q^{2}+q^{\prime
2}-k^{2}-l^{2}}\delta(\textbf{q}+\textbf{q}^{\prime}-\textbf{k}-\textbf{l})+\int
d\textbf{l}d\textbf{k}d\textbf{q}d\textbf{q}^{\prime}\frac{n_{l}n_{k}n_{q}n_{q^{\prime}}}{q^{2}+q^{\prime
2}-k^{2}-l^{2}}\delta(\textbf{q}+\textbf{q}^{\prime}-\textbf{k}-\textbf{l})$
$\displaystyle=\int
d\textbf{k}d\textbf{q}d\textbf{q}^{\prime}\frac{n_{k}n_{q}n_{q^{\prime}}}{-2k^{2}+2\textbf{k}(\textbf{q}+\textbf{q}^{\prime})-2\textbf{q}\textbf{q}^{\prime}}+\int
d\textbf{l}d\textbf{k}d\textbf{q}\frac{n_{l}n_{k}n_{q}}{2q^{2}-2\textbf{q}(\textbf{k}+\textbf{l})+2\textbf{k}\textbf{l}}$
$\displaystyle=-\frac{1}{2}\int
d\textbf{k}d\textbf{q}d\textbf{q}^{\prime}\frac{n_{k}n_{q}n_{q^{\prime}}}{k^{2}-\textbf{k}(\textbf{q}+\textbf{q}^{\prime})+\textbf{q}\textbf{q}^{\prime}}+\frac{1}{2}\int
d\textbf{q}d\textbf{k}d\textbf{l}\frac{n_{q}n_{k}n_{l}}{q^{2}-\textbf{q}(\textbf{k}+\textbf{l})+\textbf{k}\textbf{l}}=0$

With this, $I_{2}$ becomes

$I_{2}(k_{F,\lambda_{1}},k_{F,\lambda_{2}})=\frac{3}{16\pi^{5}}\int
d\textbf{l}n_{l}\int
d\textbf{k}n_{k}\bigg{[}SI(P,R,k_{F,\lambda_{1}})+SI(P,R,k_{F,\lambda_{2}})\bigg{]}$
(45)

We will integrate only one SI; the other one will be the same but
interchanging the Fermi momenta. In the end, we will add both expressions. The
inner part of the integral is

$SI(k,l,k_{F,\lambda})=\int
2d\textbf{q}d\textbf{q}^{\prime}\frac{n_{q}}{q^{2}+q^{\prime
2}-k^{2}-l^{2}}\delta(\textbf{q}+\textbf{q}^{\prime}-\textbf{k}-\textbf{l})=\int
d\textbf{q}\frac{n_{q}}{q^{2}-\textbf{q}\cdot(\textbf{k}+\textbf{l})+\textbf{k}\cdot\textbf{l}}$
$=\int_{0}^{k_{F,\lambda}}q^{2}dq\int_{0}^{\pi}\sin{\theta}d\theta\frac{2\pi}{q^{2}-q|\textbf{k}+\textbf{l}|+\textbf{k}\cdot\textbf{l}}=2\pi\int_{0}^{k_{F,\lambda}}q^{2}dq\frac{1}{q|\textbf{k}+\textbf{l}|}\ln{\bigg{|}\frac{q^{2}+q|\textbf{k}+\textbf{l}|+\textbf{k}\cdot\textbf{l}}{q^{2}-q|\textbf{k}+\textbf{l}|+\textbf{k}\cdot\textbf{l}}\bigg{|}}$
$=\frac{2\pi}{|\textbf{k}+\textbf{l}|}\bigg{\\{}\bigg{(}\frac{k_{F,\lambda}^{2}}{2}-\frac{k^{2}+l^{2}}{4}\bigg{)}\ln{\bigg{|}\frac{k_{F,\lambda}^{2}+k_{F,\lambda}|\textbf{k}+\textbf{l}|+\textbf{k}\cdot\textbf{l}}{k_{F,\lambda}^{2}-k_{F,\lambda}|\textbf{k}+\textbf{l}|+\textbf{k}\cdot\textbf{l}}\bigg{|}}$
$-\frac{|\textbf{k}+\textbf{l}||\textbf{k}-\textbf{l}|}{4}\ln{\bigg{|}\frac{k_{F,\lambda}^{2}+k_{F,\lambda}|\textbf{k}-\textbf{l}|-\textbf{k}\cdot\textbf{l}}{k_{F,\lambda}^{2}-k_{F,\lambda}|\textbf{k}-\textbf{l}|-\textbf{k}\cdot\textbf{l}}\bigg{|}}+k_{F,\lambda}|\textbf{k}+\textbf{l}|\bigg{\\}}$
(46)

We write everything in terms of the modules of k and l and the angle between
them. In the angular part, we change the variable to $x=-\cos{\theta}$,

$\displaystyle|\textbf{k}+\textbf{l}|=\sqrt{k^{2}+l^{2}+2kl\cos{\theta}}=\sqrt{k^{2}+l^{2}-2klx}$
(47)
$\displaystyle|\textbf{k}-\textbf{l}|=\sqrt{k^{2}+l^{2}-2kl\cos{\theta}}=\sqrt{k^{2}+l^{2}+2klx}$
$\displaystyle\textbf{k}\cdot\textbf{l}=kl\cos{\theta}=-klx$
$\displaystyle\int d\textbf{l}n_{l}\int
d\textbf{k}n_{k}=2(2\pi)^{2}\int_{0}^{k_{F,\lambda_{1}}}k^{2}dk\int_{0}^{k_{F,\lambda_{2}}}l^{2}dl\int_{0}^{\pi}\sin{\theta}d\theta=2(2\pi)^{2}\int_{0}^{k_{F,\lambda_{1}}}k^{2}dk\int_{0}^{k_{F,\lambda_{2}}}l^{2}dl\int_{-1}^{1}dx$

We take out the $2\pi$ from SI and integrate over x,

$\displaystyle\frac{1}{2\pi}\int
dxSI(k,l,q)=\frac{2}{3}qx+\frac{1}{kl}\bigg{(}\frac{k^{2}}{4}+\frac{l^{2}}{4}-\frac{q^{2}}{2}\bigg{)}\sqrt{k^{2}+l^{2}-2klx}\ln{\bigg{|}\frac{q^{2}+q\sqrt{k^{2}+l^{2}-2klx}-klx}{q^{2}-q\sqrt{k^{2}+l^{2}-2klx}-klx}\bigg{|}}$
(48)
$\displaystyle-\frac{1}{12kl}(k^{2}+l^{2}+2klx)^{3/2}\ln{\bigg{|}\frac{q^{2}+q\sqrt{k^{2}+l^{2}+2klx}+klx}{q^{2}-q\sqrt{k^{2}+l^{2}+2klx}+klx}\bigg{|}}+\frac{k^{4}+l^{4}+2k^{2}l^{2}-2k^{2}q^{2}-2l^{2}q^{2}+q^{4}}{3kl\sqrt{k^{2}+l^{2}-q^{2}}}\ln{\bigg{|}\frac{q\sqrt{k^{2}+l^{2}-q^{2}}+klx}{q\sqrt{k^{2}+l^{2}-q^{2}}-klx}\bigg{|}}$
$\displaystyle-\frac{q^{3}}{3kl}\bigg{(}\ln{\big{|}q\sqrt{k^{2}+l^{2}-q^{2}}+klx\big{|}}+\ln{\big{|}q\sqrt{k^{2}+l^{2}-q^{2}}-klx\big{|}}\bigg{)}$

The first external integral is

$\displaystyle
SX(k,l,q)=\frac{1}{2\pi}\int_{-1}^{1}dxSI(k,l,q)=\frac{4}{3}qx+\frac{(k-l)(2k^{2}+2l^{2}-kl-3q^{2})}{6kl}\ln{\bigg{|}\frac{q^{2}+q(k-l)-kl}{q^{2}-q(k-l)-kl}\bigg{|}}$
(49)
$\displaystyle-\frac{(k+l)(2k^{2}+2l^{2}+kl-3q^{2})}{6kl}\ln{\bigg{|}\frac{q^{2}+q(k+l)+kl}{q^{2}-q(k+l)+kl}\bigg{|}}+\frac{4}{6kl}(k^{2}+l^{2}-q^{2})^{3/2}\ln{\bigg{|}\frac{q\sqrt{k^{2}+l^{2}-q^{2}}+kl}{q\sqrt{k^{2}+l^{2}-q^{2}}-kl}\bigg{|}}$

The next step is the integration over $k$ and $l$. We will proceed in two ways
(both producing the same result):
1) $\int_{0}^{k}dk\int_{0}^{q}dlk^{2}l^{2}SX(k,l,q)$,
2) $\lim_{k\to q}\int_{0}^{k}dk\int_{0}^{l}dlk^{2}l^{2}SX(k,l,q)$.

1)

$\displaystyle\int
dlk^{2}l^{2}SX(k,l,q)=\frac{k}{120}\bigg{[}4klq(-4k^{2}+7q^{2})+44l^{3}klq+(16k^{5}-40k^{3}q^{2}+30kq^{4}+40k^{3}l^{2}+30kl^{4}-60kl^{2}q^{2})\ln{\bigg{|}\frac{l-q}{l+q}\bigg{|}}$
(50)
$\displaystyle+(16l^{5}+40l^{3}(k-q)(k+q))\ln{\bigg{|}\frac{k-q}{k+q}\bigg{|}}+16(k^{2}+l^{2}-q^{2})^{5/2}\ln{\bigg{|}\frac{q\sqrt{k^{2}+l^{2}-q^{2}}+kl}{q\sqrt{k^{2}+l^{2}-q^{2}}-kl}\bigg{|}}\bigg{]}$

It can be easily checked that the function (50) is zero when $l=0$. Also, we
will need to make use of the following limit,

$\lim_{l\to
q}\ln{\bigg{|}\frac{q\sqrt{k^{2}+l^{2}-q^{2}}+kl}{q\sqrt{k^{2}+l^{2}-q^{2}}-kl}\bigg{|}}=-\lim_{l\to
q}\ln{\bigg{|}\frac{l-q}{l+q}\bigg{|}}-\ln{\bigg{|}\frac{k-q}{k+q}\bigg{|}}-2\ln{\bigg{|}\frac{k+q}{k}\bigg{|}}$
(51)

The definite integral becomes

$\displaystyle\int_{0}^{q}dlk^{2}l^{2}SX(k,l,q)=\frac{k}{120}\bigg{[}-16k^{3}q^{2}+72kq^{4}+(16q^{5}+40q^{3}(k-q)(k+q)-16k^{5})\ln{\bigg{|}\frac{k-q}{k+q}\bigg{|}}-32k^{5}\ln{\bigg{|}\frac{k+q}{k}\bigg{|}}\bigg{]}$
(52)

Finally, the integral over $k$ is

$\int_{0}^{k}dk\int_{0}^{q}dlk^{2}l^{2}SX(k,l,q)=\frac{1}{420}\bigg{[}-8k^{5}q^{2}+66k^{3}q^{4}+30kq^{6}+(-8k^{7}+35k^{4}q^{3}-42k^{2}q^{5}+15q^{7})\ln{\bigg{|}\frac{k-q}{k+q}\bigg{|}}-16k^{7}\ln{\bigg{|}\frac{k+q}{k}\bigg{|}}\bigg{]}$
(53)

2)

$\displaystyle\int_{0}^{k}dk\int_{0}^{l}dlk^{2}l^{2}SX(k,l,q)=\frac{1}{840}\bigg{[}-16k^{5}lq-16kl^{5}q+88k^{3}l^{3}q+44k^{3}lq^{3}+44kl^{3}q^{3}+32klq^{5}$
(54)
$\displaystyle+(16l^{7}-56l^{5}q^{2}+70l^{3}q^{4}+70k^{4}l^{3}+56k^{2}l^{5}-140k^{2}l^{3}q^{2})\ln{\bigg{|}\frac{k-q}{k+q}\bigg{|}}$
$\displaystyle+(16k^{7}-56k^{5}q^{2}+70k^{3}q^{4}+70k^{3}l^{4}+56k^{5}l^{2}-140k^{3}l^{2}q^{2})\ln{\bigg{|}\frac{l-q}{l+q}\bigg{|}}$
$\displaystyle+16(k^{2}+l^{2}-q^{2})^{7/2}\ln{\bigg{|}\frac{q\sqrt{k^{2}+l^{2}-q^{2}}+kl}{q\sqrt{k^{2}+l^{2}-q^{2}}-kl}\bigg{|}}$

As the expression above is symmetrical with respect to $k$ and $l$, it does
not matter which variable we choose to perform the limit to $(k,l)\to q$ since
both results are formally equivalent. As we want to check that this second
method gives the same expression as the one found in method 1, we will set the
limit $l$ going to $q$,

$\displaystyle\lim_{l\to q}\int_{0}^{k}dk\int_{0}^{l}dlk^{2}l^{2}SX(k,l,q)=$
(55)
$\displaystyle\frac{1}{840}\bigg{[}-16k^{5}q^{2}+132k^{3}q^{4}+60kq^{6}+(-16k^{7}+70k^{4}q^{3}-84k^{2}q^{5}+30q^{7})\ln{\bigg{|}\frac{k-q}{k+q}\bigg{|}}-32k^{7}\ln{\bigg{|}\frac{k+q}{k}\bigg{|}}\bigg{]}$

We recover indeed the same expression.

Coming back to the integral we had at the beginning,

$I_{2}(k_{F,\lambda_{1}},k_{F,\lambda_{2}})=\frac{3}{16\pi^{5}}\int
d\textbf{l}n_{l}\int
d\textbf{k}n_{k}\bigg{[}SI(P,R,k_{F,\lambda_{1}})+SI(P,R,k_{F,\lambda_{2}})\bigg{]}$
(56)

We substitute the results we have obtained for the integrals (Eqs. (53) or
(55)) and we recover the factor $2(2\pi)^{3}$ coming from the angular
integrals.

$\displaystyle I_{2}=\frac{3}{16\pi^{5}}2(2\pi)^{3}\bigg{(}\lim_{k\to
q}\int_{0}^{k}dk\int_{0}^{l}dlk^{2}l^{2}SX(k,l,q)+\lim_{l\to
q^{\prime}}\int_{0}^{k}dk\int_{0}^{l}dlk^{2}l^{2}SX(k,l,q^{\prime})\bigg{)}$
(57)
$\displaystyle=\frac{3}{16\pi^{5}}\frac{2(2\pi)^{3}}{420}\bigg{[}-8k^{5}l^{2}-8l^{5}k^{2}+66k^{3}l^{4}+66l^{3}k^{4}+30kl^{6}+30lk^{6}$
$\displaystyle+(7k^{7}+7l^{7}+35k^{4}l^{3}+35l^{4}k^{3}-42k^{2}l^{5}-42l^{2}k^{5})\ln{\bigg{|}\frac{k-l}{k+l}\bigg{|}}-16k^{7}\ln{\bigg{|}\frac{k+l}{k}\bigg{|}}-16l^{7}\ln{\bigg{|}\frac{k+l}{l}\bigg{|}}\bigg{]}$
$\displaystyle=\frac{1}{140\pi^{2}}\bigg{[}2kl(k+l)(15k^{4}-19k^{3}l+52k^{2}l^{2}-19kl^{3}+15l^{4})$
$\displaystyle+7(k+l)(k-l)^{4}(k^{2}+3kl+l^{2})\ln{\bigg{|}\frac{k-l}{k+l}\bigg{|}}-16\bigg{(}k^{7}\ln{\bigg{|}\frac{k+l}{k}\bigg{|}}+l^{7}\ln{\bigg{|}\frac{k+l}{l}\bigg{|}}\bigg{)}\bigg{]}$

After putting everything together and rearranging terms, $I_{2}$ is written in
terms of k and l as

$\displaystyle
I_{2}=\frac{4}{35\pi^{2}}\bigg{[}\frac{1}{8}kl(k+l)(15k^{4}-19k^{3}l+52k^{2}l^{2}-19kl^{3}+15l^{4})$
(58)
$\displaystyle+\frac{7}{16}(k+l)(k-l)^{4}(k^{2}+3kl+l^{2})\ln{\bigg{|}\frac{k-l}{k+l}\bigg{|}}-\bigg{(}k^{7}\ln{\bigg{|}\frac{k+l}{k}\bigg{|}}+l^{7}\ln{\bigg{|}\frac{k+l}{l}\bigg{|}}\bigg{)}\bigg{]}$

Replacing $k$ and $l$ by the Fermi momenta ($k_{F,\lambda_{1}}$ and
$k_{F,\lambda_{2}}$),

$\displaystyle
I_{2}=\frac{4}{35\pi^{2}}\bigg{[}\frac{1}{8}k_{F,\lambda_{1}}k_{F,\lambda_{2}}(k_{F,\lambda_{1}}+k_{F,\lambda_{2}})(15k_{F,\lambda_{1}}^{4}-19k_{F,\lambda_{1}}^{3}k_{F,\lambda_{2}}+52k_{F,\lambda_{1}}^{2}k_{F,\lambda_{2}}^{2}-19k_{F,\lambda_{1}}k_{F,\lambda_{2}}^{3}+15k_{F,\lambda_{2}}^{4})$
(59)
$\displaystyle+\frac{7}{16}(k_{F,\lambda_{1}}+k_{F,\lambda_{2}})(k_{F,\lambda_{1}}-k_{F,\lambda_{2}})^{4}(k_{F,\lambda_{1}}^{2}+3k_{F,\lambda_{1}}k_{F,\lambda_{2}}+k_{F,\lambda_{2}}^{2})\ln{\bigg{|}\frac{k_{F,\lambda_{1}}-k_{F,\lambda_{2}}}{k_{F,\lambda_{1}}+k_{F,\lambda_{2}}}\bigg{|}}$
$\displaystyle-\bigg{(}k_{F,\lambda_{1}}^{7}\ln{\bigg{|}\frac{k_{F,\lambda_{1}}+k_{F,\lambda_{2}}}{k_{F,\lambda_{1}}}\bigg{|}}+k_{F,\lambda_{2}}^{7}\ln{\bigg{|}\frac{k_{F,\lambda_{1}}+k_{F,\lambda_{2}}}{k_{F,\lambda_{2}}}\bigg{|}}\bigg{)}\bigg{]}$

Now, we replace the Fermi momenta by $k_{F}C_{\lambda}^{1/3}$,

$\displaystyle
I_{2}=\frac{4k_{F}^{7}}{35\pi^{2}}\bigg{[}\frac{1}{8}C_{\lambda_{1}}^{1/3}C_{\lambda_{2}}^{1/3}\big{(}C_{\lambda_{1}}^{1/3}+C_{\lambda_{2}}^{1/3}\big{)}\big{(}15C_{\lambda_{1}}^{4/3}-19C_{\lambda_{1}}C_{\lambda_{2}}^{1/3}+52C_{\lambda_{1}}^{2/3}C_{\lambda_{2}}^{2/3}-19C_{\lambda_{1}}^{1/3}C_{\lambda_{2}}+15C_{\lambda_{2}}^{4/3}\big{)}$
(60)
$\displaystyle+\frac{7}{16}\big{(}C_{\lambda_{1}}^{1/3}+C_{\lambda_{2}}^{1/3}\big{)}\big{(}C_{\lambda_{1}}^{1/3}-C_{\lambda_{2}}^{1/3}\big{)}^{4}\big{(}C_{\lambda_{1}}^{2/3}+3C_{\lambda_{1}}^{1/3}C_{\lambda_{2}}^{1/3}+C_{\lambda_{2}}^{2/3}\big{)}\ln{\bigg{|}\frac{C_{\lambda_{1}}^{1/3}-C_{\lambda_{2}}^{1/3}}{C_{\lambda_{1}}^{1/3}+C_{\lambda_{2}}^{1/3}}\bigg{|}}$
$\displaystyle-\bigg{(}C_{\lambda_{1}}^{7/3}\ln{\bigg{|}\frac{C_{\lambda_{1}}^{1/3}+C_{\lambda_{2}}^{1/3}}{C_{\lambda_{1}}^{1/3}}\bigg{|}}+C_{\lambda_{2}}^{7/3}\ln{\bigg{|}\frac{C_{\lambda_{1}}^{1/3}+C_{\lambda_{2}}^{1/3}}{C_{\lambda_{2}}^{1/3}}\bigg{|}}\bigg{)}\bigg{]}$

Finally, in terms of $y=(C_{\lambda_{1}}/C_{\lambda_{2}})^{1/3}$,

$\displaystyle
I_{2}=\frac{4k_{F}^{7}}{35\pi^{2}}C_{\lambda_{1}}C_{\lambda_{2}}\frac{C_{\lambda_{1}}^{1/3}+C_{\lambda_{2}}^{1/3}}{2}\bigg{[}\frac{1}{4}\big{(}15y^{2}-19y+52-19y^{-1}+15y^{-2}\big{)}+\frac{7}{8}y^{-2}\big{(}y-1\big{)}^{4}\big{(}y+3+y^{-1}\big{)}\ln{\bigg{|}\frac{1-y}{1+y}\bigg{|}}$
(61)
$\displaystyle-\frac{2y^{4}}{1+y}\ln{\bigg{|}1+y^{-1}\bigg{|}}-\frac{2y^{-4}}{1+y^{-1}}\ln{\bigg{|}1+y\bigg{|}}\bigg{]}$

## References

  * Giuliani and Vignale [2005] G. F. Giuliani and G. Vignale,  _Quantum Theory of the Electron Liquid_ (Cambridge University Press, The Edinburgh Building, Cambridge, 2005).
  * Stoner [1933] E. Stoner, Atomic moments in ferromagnetic metals and alloys with non-ferromagnetic elements, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 15, 1018 ((1933)).
  * Pfleiderer  _et al._ [2001] C. Pfleiderer, S. Julian, and G. Lonzarich, Non-fermi-liquid nature of the normal state of itinerant-electron ferromagnets, Nature 414, 427 (2001).
  * Leduc [1990] M. Leduc, Spin polarized helium-3, a playground in many domains of physics, Journal de Physique Colloques 51, C6 (1990).
  * Jo  _et al._ [2009] G.-B. Jo, Y.-R. Lee, J.-H. Choi, C. A. Christensen, T. H. Kim, J. H. Thywissen, D. E. Pritchard, and W. Ketterle, Itinerant ferromagnetism in a fermi gas of ultracold atoms, Science 325, 1521 (2009).
  * Sanner  _et al._ [2012] C. Sanner, E. J. Su, W. Huang, A. Keshet, J. Gillen, and W. Ketterle, Correlations and pair formation in a repulsively interacting fermi gas, Phys. Rev. Lett. 108, 240404 (2012).
  * Valtolina  _et al._ [2017] G. Valtolina, F. Scazza, A. Amico, A. Burchianti, A. Recati, T. Enss, and M. Inguscio, Exploring the ferromagnetic behaviour of a repulsive fermi gas through spin dynamics, Nat. Phys. 13, 704 (2017).
  * Conduit  _et al._ [2009] G. J. Conduit, A. G. Green, and B. D. Simons, Inhomogeneous phase formation on the border of itinerant ferromagnetism, Phys. Rev. Lett. 103, 207201 (2009).
  * Pilati  _et al._ [2010] S. Pilati, G. Bertaina, S. Giorgini, and M. Troyer, Itinerant ferromagnetism of a repulsive atomic fermi gas: A quantum monte carlo study, Phys. Rev. Lett. 105, 030405 (2010).
  * Chang  _et al._ [2011] S.-Y. Chang, M. Randeria, and N. Trivedi, Ferromagnetism in the upper branch of the feshbach resonance and the hard-sphere fermi gas, Proc. Natl. Acad. Sci. USA 108, 51 (2011).
  * Conduit and Simons [2009] G. J. Conduit and B. D. Simons, Itinerant ferromagnetism in an atomic fermi gas: Influence of population imbalance, Phys. Rev. A 79, 053606 (2009).
  * Massignan  _et al._ [2013] P. Massignan, Z. Yu, and G. M. Bruun, Itinerant ferromagnetism in a polarized two-component fermi gas, Phys. Rev. Lett. 110, 230401 (2013).
  * Cui and Zhai [2010] X. Cui and H. Zhai, Stability of a fully magnetized ferromagnetic state in repulsively interacting ultracold fermi gases, Phys. Rev. A 81, 041602 (2010).
  * Arias de Saavedra  _et al._ [2012] F. Arias de Saavedra, F. Mazzanti, J. Boronat, and A. Polls, Ferromagnetic transition of a two-component fermi gas of hard spheres, Phys. Rev. A 85, 033615 (2012).
  * Pilati  _et al._ [2021] S. Pilati, G. Orso, and G. Bertaina, Quantum monte carlo simulations of two-dimensional repulsive fermi gases with population imbalance, Phys. Rev. A 103, 063314 (2021).
  * Huang and Yang [1957] K. Huang and C. N. Yang, Quantum-mechanical many-body problem with hard-sphere interaction, Phys. Rev. 105, 767 ((1957)).
  * Lee and Yang [1957] T. D. Lee and C. N. Yang, Many-body problem in quantum mechanics and quantum statistical mechanics, Phys. Rev. 105, 1119 ((1957)).
  * Galitskii [1958] V. M. Galitskii, The energy spectrum of a non-ideal fermi gas, Journal of Experimental and Theoretical Physics 7, 104 ((1958)).
  * Abrikosov and Khalatnikov [1958] A. A. Abrikosov and I. M. Khalatnikov, Concerning a model for a non-ideal fermi gas, Journal of Experimental and Theoretical Physics 6, 888 (1958).
  * Cazalilla  _et al._ [2009] M. A. Cazalilla, A. F. Ho, and M. Ueda, Ultracold gases of ytterbium: ferromagnetism and mott states in an su(6) fermi system, New. J. Phys. 11, 103033 (2009).
  * Cazalilla and Rey [2014] M. A. Cazalilla and A. M. Rey, Ultracold fermi gases with emergent su(n) symmetry, Rep. Prog. Phys. 77, 124401 (2014).
  * Scazza  _et al._ [2014] F. Scazza, C. Hofrichter, M. Höfer, P. C. D. Groot, I. Bloch, and S. Fölling, Observation of two-orbital spin-exchange interactions with ultracold su(n)-symmetric fermions, Nat. Phys. 10, 779 (2014).
  * Xu  _et al._ [2018] S. Xu, J. T. Barreiro, Y. Wang, and C. Wu, Interaction effects with varying n in su(n) symmetric fermion lattice systems, Phys. Rev. Lett. 121, 167205 (2018).
  * He  _et al._ [2020] C. He, Z. Ren, B. Song, E. Zhao, J. Lee, Y.-C. Zhang, S. Zhang, and G.-B. Jo, Collective excitations in two-dimensional su(n) fermi gases with tunable spin, Phys. Rev. Research 2, 012028(R) (2020).
  * Huang  _et al._ [2020] C.-H. Huang, Y. Takasu, Y. Takahashi, and M. A. Cazalilla, Suppression and control of prethermalization in multicomponent fermi gases following a quantum quench, Phys. Rev. A 101, 053620 (2020).
  * Sonderhouse  _et al._ [2020] L. Sonderhouse, C. Sanner, R. B. Hutson, A. Goban, T. Bilitewski, L. Yan, W. R. Milner, A. M. Rey, and J. Ye, Thermodynamics of a deeply degenerate su(n)-symmetric fermi gas, Nat. Phys. 16, 1216 (2020).
  * Ibarra-García-Padilla  _et al._ [2021] E. Ibarra-García-Padilla, S. Dasgupta, H.-T. Wei, S. Taie, Y. Takahashi, R. T. Scalettar, and K. R. A. Hazzard, Universal thermodynamics of an su(n) fermi-hubbard model, Phys. Rev. A 104, 043316 (2021).
  * Pagano  _et al._ [2014] G. Pagano, M. Mancini, G. Cappellini, P. Lombardi, F. Schäfer, H. Hu, X. J. Liu, J. Catani, C. Sias, M. Inguscio, and L. Fallani, A one-dimensional liquid of fermions with tunable spin, Nat. Phys. 10, 198 (2014).
  * Goban  _et al._ [2018] A. Goban, R. B. Hutson, G. E. Marti, S. L. Campbell, M. A. Perlin, P. S. Julienne, J. P. D’Incao, A. M. Rey, and J. Ye, Emergence of multi-body interactions in a fermionic lattice clock, Nature (London) 563, 369 (2018).
  * Bishop [1973] R. F. Bishop, Ground-state energy of a dilute fermi gas, Annals of Physics 77, 106 (1973).
  * Baker-Jr. [1965] G. A. Baker-Jr., Study of the perturbation series for the ground state of a many-fermion system, Phys. Rev. 140B, 9 ((1965)).
  * Chankowski and Wojtkiewicz [2021] P. Chankowski and J. Wojtkiewicz, On the ground-state energy of a mixture of two different oppositely polarized fermionic gases, arXiv:2111.14927 (2021).
  * Kanno [1970] S. Kanno, Criterion for the ferromagnetism of hard sphere fermi liquid. ii, Prog. Theor. Phys. 44, 813 (1970).
  * He and Huang [2012] L. He and X.-G. Huang, Nonperturbative effects on the ferromagnetic transition in repulsive fermi gases, Phys. Rev. A 85, 043624 (2012).
  * Tan [2008] S. Tan, Large momentum part of a strongly correlated fermi gas, Annals of Physics 323, 2971 (2008).
  * Song  _et al._ [2020] B. Song, Y. Yan, C. He, Z. Ren, Q. Zhou, and G.-B. Jo, Evidence for bosonization in a three-dimensional gas of $\mathrm{SU}(n)$ fermions, Phys. Rev. X 10, 041053 (2020).