# Galaxy rotation curves and universal scaling relations: comparison between
phenomenological and fermionic dark matter profiles

A. Krut  ICRANet, Piazza della Repubblica 10, I-65122 Pescara, Italy  C. R.
Argüelles  ICRANet, Piazza della Repubblica 10, I-65122 Pescara, Italy
Instituto de Astrofísica de La Plata, UNLP-CONICET, Paseo del Bosque s/n
B1900FWA La Plata, Argentina  P.-H. Chavanis  Laboratoire de Physique
Théorique, Université de Toulouse, CNRS, UPS, France  J. A. Rueda  ICRANet,
Piazza della Repubblica 10, I-65122 Pescara, Italy  Dip. di Fisica, Sapienza
Università di Roma, Piazzale Aldo Moro 5, I-00185 Rome, Italy  ICRANet-
Ferrara, Dip. di Fisica e Scienze della Terra, Università degli Studi di
Ferrara, Via Saragat 1, 44122 Ferrara, Italy  INAF, Istituto de Astrofisica e
Planetologia Spaziali, Via Fosso del Cavaliere 100, 00133 Rome, Italy  R.
Ruffini  ICRANet, Piazza della Repubblica 10, I-65122 Pescara, Italy  Dip. di
Fisica, Sapienza Università di Roma, Piazzale Aldo Moro 5, I-00185 Rome, Italy
ICRANet-Ferrara, Dip. di Fisica e Scienze della Terra, Università degli Studi
di Ferrara, Via Saragat 1, 44122 Ferrara, Italy  INAF, Istituto de Astrofisica
e Planetologia Spaziali, Via Fosso del Cavaliere 100, 00133 Rome, Italy

###### Abstract

Galaxies show different halo scaling relations such as the Radial Acceleration
Relation, the Mass Discrepancy Acceleration Relation (MDAR) or the dark matter
Surface Density Relation (SDR). At difference with traditional studies using
phenomenological $\Lambda$CDM halos, we analyze the above relations assuming
that dark matter (DM) halos are formed through a Maximum Entropy Principle
(MEP) in which the fermionic (quantum) nature of the DM particles is dully
accounted for. For the first time a competitive DM model based on first
physical principles, such as (quantum) statistical-mechanics and
thermodynamics, is tested against a large data-set of galactic observables. In
particular, we compare the fermionic DM model with empirical DM profiles: the
NFW model, a generalized NFW model accounting for baryonic feedback, the
Einasto model and the Burkert model. For this task, we use a large sample of
$120$ galaxies taken from the Spitzer Photometry and Accurate Rotation Curves
(SPARC) data-set, from which we infer the DM content to compare with the
models. We find that the Radial Acceleration Relation and MDAR are well
explained by all the models with comparable accuracy, while the fits to the
individual rotation curves, in contrast, show that cored DM halos are
statistically preferred with respect to the cuspy NFW profile. However, very
different physical principles justify the flat inner halo slope in the most
favored DM profiles: while generalized NFW or Einasto models rely on complex
baryonic feedback processes, the MEP scenario involves a quasi-thermodynamic
equilibrium of the DM particles.

dark matter - galaxies: structure - galaxies: fundamental parameters - galaxy:
kinematics and dynamics

tablenum SIXtablenum

##  1 Introduction

How the total gravitating mass distributes with respect to the luminous mass
on galaxy scales is an open question which has regained much attention in the
last decade thanks to the vast data-sets covering broader radial extents
across different Hubble types (de Blok et al., 2008; Cappellari et al., 2011;
Lelli et al., 2016). Several universal relations exist between different pairs
of structural galaxy parameters, which refer either  (i) to the outer regions
of galaxies such as the baryonic Tully-Fisher relation (BTFR) (McGaugh et al.,
2000), the DM Surface Density Relation (SDR) (Donato et al., 2009), the Radial
Acceleration Relation (McGaugh et al., 2016), and the Mass Discrepancy
Acceleration Relation (MDAR) (McGaugh, 2004), which are indeed all closely
related (McGaugh, 2004; Salucci, 2016, 2018);  (ii) to their central regions
such as the $M$-$\sigma$ relation between the bulge’s dispersion velocity and
the central object mass (Ferrarese & Merritt, 2000); or  (iii) to a
combination of both regimes such as the Ferrarese relation (Ferrarese, 2002;
Kormendy & Bender, 2011; Bogdán & Goulding, 2015) between the total halo mass
and its supermassive central object mass.

Actual attempts for a unified understanding of many of the above scaling
relations are typically given in terms of phenomenological halos obtained from
N-body simulations within $\Lambda$CDM (see e.g. Navarro et al., 2017; Ludlow
et al., 2017; Sales et al., 2017; Salucci, 2018).

However, when DM halos are formed through a MEP for collisionless systems of
self-gravitating fermions (Chavanis, 1998; Chavanis et al., 2015a; Chavanis,
2020; Argüelles et al., 2021), it leaves place to novel theoretical
predictions in the phenomenology of real galaxies (Argüelles et al., 2018,
2019; Becerra-Vergara et al., 2020, 2021; Argüelles et al., 2021, 2022b) such
as:

  1. (1)

DM fermions with a finite-temperature can be in a diluted (Boltzmannian-like)
regime or become semi-degenerate. The corresponding DM halos can be,
respectively, King-like or may develop a dense and degenerate compact core at
the center of such a halo (Chavanis et al., 2015a; Argüelles et al., 2021). A
fully relativistic model in which this more general core–halo profiles arise,
is usually referred in the literature as the Ruffini-Argüelles-Rueda (RAR)
model (Ruffini et al., 2015; Argüelles et al., 2018, 2019; Becerra-Vergara et
al., 2020, 2021; Argüelles et al., 2021, 2022a, 2022b). In either case, these
fermionic DM profiles are cored (i.e. develop an extended plateau on halo
scales similar to the Burkert DM profile as shown in Fig. 1), thereby not
suffering from the core-cusp problem associated with the standard $\Lambda$CDM
cosmology (Bullock & Boylan-Kolchin, 2017). This cored feature seems to be a
general conclusion reached for any DM profile which has reached a (quasi)
thermodynamic equilibrium in cosmology (Sánchez Almeida & Trujillo, 2021).

  2. (2)

The relevance of the fermionic solutions with a degenerate DM core surrounded
by a diluted halo (named from now on as core–halo profiles) imply different
consequences: (a) the core might become so densely packed that above a
threshold — the critical mass — the quantum pressure can not support it any
longer against its own weight, leading to the gravitational core-collapse into
a supermassive black hole (SMBH; Alberti & Chavanis, 2020; Argüelles et al.,
2021). For DM particle masses of $\mathcal{O}(10)$ keV, this result provides,
for large enough galaxies, a novel SMBH formation mechanism in the early
Universe as proposed in Argüelles et al. (2021); and (b) for core masses below
its critical value, it exist a set of free parameters in the fermionic model,
such that the quantum cores correlate with their outer halos (Argüelles et
al., 2019) explaining the relation between the total mass $M_{\rm tot}$ of
large enough galaxies with the (presumably) embedded BH mass $M_{\rm BH}$,
i.e. the Ferrarese relation (see point (iii) above).

  3. (3)

For small enough galaxies, as in the case of typical dwarfs, the degenerate
core cannot collapse towards a BH (i.e. the total mass of the galaxy is below
the critical mass of collapse) and thus it remains in a core–halo state where
the central nucleus still mimics the effects of a singularity — with core
masses in the range of Intermediate mass BHs — while the outer halo explains
the rotation curves (RCs) (Argüelles et al., 2019, 2021, 2022a).

  4. (4)

In the more extreme case where the fermions are fully-degenerate (i.e. when
they are treated under the $T\to 0$ approximation), the corresponding halos
are polytropic and may be only applicable to dwarfs (Domcke & Urbano, 2015).

Moreover, it was recently demonstrated that fermionic halos obtained via this
MEP mechanism can arise in a cosmological framework, and remain
thermodynamically and dynamically stable during the life of the Universe
(Argüelles et al., 2021). Indeed, it was there shown the self-consistency of
the approach in the sense that the nature and mass of the DM particles
involved in the linear matter power spectrum — calculated in Argüelles et al.
(2021) within a CLASS code for $\mathcal{O}(10)$ keV fermions — are the very
same building blocks at the basis of the virialized DM configurations with its
inherent effects in the DM profiles. Even if the MEP has been applied in
Argüelles et al. (2021) within a warm DM (WDM) cosmology — and is used here
for finite temperature fermions (see 3.1) — this is not univocal in general.
Other MEP such as the one used in Sánchez Almeida & Trujillo (2021) has been
also applied for classical (self-gravitating) particles within cold DM (CDM)
cosmologies.

Thus, the main purpose of this work is to analyze most of the galaxy
relations, as discussed above in (i)-(iii), together with a large set of
observed RCs provided by the SPARC data. In a forthcoming work we will
reconsider this data-set together with observations coming from the central
regions of galaxies (e.g. supermassive BHs), to include the $M_{\rm BH}-M_{\rm
tot}$ relation in the study.

We assume that DM halos are formed through a MEP in which the fermionic
(quantum) nature of the DM particles is dully accounted for. This takes
special interest since it is the first time a predictive model of this kind,
i.e. based on first physical principles such as (quantum) statistical-
mechanics and thermodynamics, is tested against a large set of galaxy
observables while leading to a good agreement with observations. To see how
well this fermionic model can reproduce the given observables, it will be
compared with most of the commonly used DM models used in the literature.
Further, we check which set of observed data can (or cannot) discriminate the
goodness of the competing models.

For this task, we will first focus on the Radial Acceleration Relation — a
non-linear correlation between the radial acceleration caused by the total
matter and the one generated by its baryonic component only (see section 4.1)
— and on the directly related MDAR.

The main motivation to start studying these acceleration relations is due to
the intense debate they have generated in the past few years about their
underlying physical origin. For instance, McGaugh et al. (2016) argues that
the Radial Acceleration Relation can be explained by a fundamental
acceleration constant $\mathfrak{a}_{0}$ in combination with a modification of
Newtonian gravity and without the necessity of any DM. One potential
explanation for such a fundamental constant comes from Modified Newtonian
Dynamics (MOND) (Kroupa, 2015; Milgrom, 2016; McGaugh et al., 2016; Li et al.,
2018), which has been used to interpret it as evidence against the
$\Lambda$CDM paradigm and in favor to the MOND theory. However, more recent
studies dedicated to analyze this universal relation within the (Bayesian)
posterior distributions on the acceleration scales of individual galaxies
(across a large sample), have provided evidence against the existence of such
a fundamental constant and in favour of $\mathfrak{a}_{0}$ to be an emergent
magnitude (see e.g. Marra et al., 2020 and references therein). On the other
hand, it has been extensively shown that the Radial Acceleration Relation is
consistent with the $\Lambda$CDM paradigm, as found either in hydrodynamical
N-body simulations (Di Cintio & Lelli, 2016; Navarro et al., 2017; Ludlow et
al., 2017; Dutton et al., 2019), or from other more phenomenological
(independent) study based on Universal Rotation Curves (Salucci, 2018).

Thus, in the first part of this work we will use the MEP approach for fermions
to evaluate how competitive it is to reproduce the above relations with
respect to other phenomenological DM halos. The galactic observables are taken
from a filtered SPARC sample of $2369$ data points (corresponding to a total
of $120$ galaxies) and then apply a non-linear least square statistical
analysis in order to check the goodness of fit of:  (a) the above fermionic DM
halos — either in the core–halo regime (Argüelles et al., 2019, 2021) or in
the purely King-like one (Argüelles et al., 2021);  (b) a baryonic feedback
motivated halo model within $\Lambda$CDM according to Di Cintio et al. (2014),
named here as DC14;  (c) the classical NFW model based on early numerical
simulations (Navarro et al., 1997);  (d) the Einasto profile (Einasto & Haud,
1989; Merritt et al., 2006), and  (e) the Burkert model (Burkert, 1995).

Outlining the structure of this paper, in section 2 we give a brief overview
about the SPARC data set, data selection and fitting procedure. In section 3,
we describe the competing DM halo models considered above. In section 4, we
present our results on the Radial Acceleration Relation by performing a
goodness of the fit for each DM model based on a filtered SPARC sample. In
section 5, we focus on fermionic halos and compare their morphology with the
other DM models. Finally, in section 6, we give a brief summary and draw the
conclusions.

We refer to the appendix for further details. In Appendix A, we focus on the
above mentioned fermionic model and analyze the relation between halo mass and
halo radius which is qualitatively consistent with DM-SDR. In Appendix B, we
provide model parameter distributions of the competing DM models.

##  2 Methodology

The data used in this work is obtained from the Spitzer Photometry and
Accurate Rotation Curves (SPARC) data-set. It contains independent
observations of the total velocity ($V_{\rm tot}$) and luminous mass
distributions which allows to infer the bulge ($V_{\rm b}$), disk ($V_{\rm
d}$), and gas ($V_{\rm g}$) velocity contributions. We extract the tangential
velocities of the DM component from the data (for each galaxy of the sample)
by subtracting the inferred baryonic components (as provided in the SPARC
data-set) from the total velocity $V_{\rm tot}$, and statistically compare
with the corresponding velocities predicted by each DM model as detailed in
next.

###  2.1 Data selection

The SPARC data-set includes $3.6\text{\,}\mathrm{\SIUnitSymbolMicro m}$ near-
infrared and $21\text{\,}\mathrm{cm}$ observations. The former traces the
stellar mass distribution (bulge and disk) while the latter traces the atomic
gas distribution and provides velocity fields from which the RCs are derived.
See Lelli et al. (2016) for a complete description of the SPARC data.

The data is distributed in separated files such as Table1.mrt (i.e. Hubble
type, inclination etc.) and Table2.mrt (i.e. RC data) and can be found at
http://astroweb.cwru.edu/SPARC/.

We extract the observed circular velocity $V_{\mathrm{tot}}$ and the baryonic
contribution $V_{\mathrm{bar}}$, composed of a bulge ($V_{\mathrm{b}}$), disk
($V_{\mathrm{d}}$) and gas component ($V_{\mathrm{g}}$). The bulge and disk
components are inferred from surface brightness observations for a given mass-
to-light ratio. The baryonic component is then given by

$V_{\mathrm{bar}}^{2}=\Upsilon_{\mathrm{b}}^{\phantom{2}}V_{\mathrm{b}}^{2}+\Upsilon_{\mathrm{d}}^{\phantom{2}}V_{\mathrm{d}}^{2}+V_{\mathrm{g}}^{2}.$
(1)

For convenience, the given velocities $V_{\mathrm{b}}$ and $V_{\mathrm{d}}$ in
the SPARC data are normalized for a mass-to-light ratio of
$1\,M_{\odot}/L_{\odot}$.

In this work we follow the same data selection criteria as done in McGaugh et
al. (2016). We choose averaged mass-to-light ratios
$\Upsilon_{\mathrm{b}}=0.7$ for all bulges and $\Upsilon_{\mathrm{d}}=0.5$ for
all disks as convenient average representatives. We exclude all galaxies with
a bad quality flag ($Q=3$) and face-on galaxies with an inclination
$i<$30\text{\,}\mathrm{\SIUnitSymbolDegree}$$. The latter is to minimize the
$\sin(i)$ corrections to the observed velocities. For all measurements we
require a minimum precision of $10\text{\,}\mathrm{\char 37\relax}$ in
velocity.

Additionally we reject all points where the baryonic velocity is greater than
$95\text{\,}\mathrm{\char 37\relax}$ of the observed velocity. This condition
is required to avoid negative velocities for the inferred DM components. It
affects mainly data points in the inner region which is dominated by baryonic
matter and strongly depends on the chosen mass-to-light factors. Therefore,
those inner points are less reliable. Finally, we exclude all remaining
galaxies with less than 6 data points to be statistically significant.

We obtain $120$ galaxies (out of $174$) with $2396$ points (of $3355$) in
total. Galaxies not fulfilling the quality criteria have such poor data that
they do not allow to gain any insights. In the worst case (e.g. too few
points, points on a nearly straight line, etc.) it is not possible to fit the
rotation curves.

###  2.2 Data fitting

We infer the circular velocity of the DM component directly from observations
via

$V_{\mathrm{DM}}^{2}=V_{\mathrm{tot}}^{2}-V_{\mathrm{bar}}^{2},$ (2)

where the uncertainty in $V_{\rm DM}$ is calculated from the uncertainties in
the other velocity components within linear error propagation theory. Note
that only the uncertainties in the total RCs ($\Delta V_{\rm tot}$) are
provided within the SPARC data-set. The uncertainty is therefore given by

$\Delta V_{\rm DM}=\left|\frac{\partial V_{\rm DM}}{\partial V_{\rm
tot}}\right|\Delta V_{\rm tot}.$ (3)

The inferred DM contribution of each galaxy then will be fitted by the
competing DM models which are described in section 3.

With this information, we use the Levenberg–Marquardt (LM) algorithm of least-
square error minimization for each dark matter component calculated by

$\chi^{2}(\mbox{\boldmath$p$})=\sum\limits_{i=1}^{N}\left[\vphantom{}\frac{V_{i}-v(r_{i},\mbox{\boldmath$p$})}{\Delta
V_{i}}\right]^{2}.$ (4)

Here, $N$ is the number of data points for a given galaxy, $V_{i}$ is the set
of inferred DM circular velocities from the data at each corresponding
measured radius $r_{i}$, while $v(r_{i},\mbox{\boldmath$p$})$ accounts for the
circular velocity at $r_{i}$ for each model parameter vector $p$ (described
below), and $\Delta V_{i}$ is the uncertainty in $V_{i}$ as given in Eq. 3.

##  3 Dark matter models

Most of the DM halo models in the literature are phenomenological, i.e.
motivated by the phenomenology of rotation curves either from observations or
from numerical N-body simulations. In contrast, we consider a fermionic DM
model based on first physical principles including statistical-mechanics and
thermodynamics.

###  3.1 Fermionic DM halos from MEP

It has been proposed by several authors (see, e.g. Chavanis, 2020 for an
exhaustive list of references) that DM halos could be made of fermions (e.g.
sterile neutrinos) in gravitational interaction. It is usually assumed that
the fermions are in a statistical equilibrium state described by the Fermi-
Dirac distribution function. However, the notion of statistical equilibrium
for systems with long-range interactions is subtle. If the fermions are non-
interacting, apart from gravitational forces, the relaxation time towards
statistical equilibrium due to gravitational encounters scales as $(N/\ln
N)t_{D}$ (see e.g. Binney & Tremaine, 2008) and exceeds the age of the
Universe by many orders of magnitude.

For example, assuming a fermion mass $mc^{2}\sim$50\text{\,}\mathrm{keV}$$, a
DM halo of mass $M\sim${10}^{11}\text{\,}\mathrm{M_{\odot}}$$ and radius
$R\sim$30\text{\,}\mathrm{kpc}$$ contains $N\sim${10}^{72}$$ fermions for a
dynamical time $t_{D}\sim 1/\sqrt{R^{3}/GM}\sim 100\,{\rm Myrs}$.

Therefore, on the Hubble time, the gas of fermions is essentially
collisionless, being described by the Vlasov-Poisson equations. Yet, it can
achieve a form of statistical equilibrium on a coarse-grained scale through a
process of violent relaxation. This concept was introduced by Lynden-Bell
(1967) in the case of collisionless stellar systems and has been exported to
DM by Kull et al. (1996) and Chavanis et al. (2015a).

Assuming ergodicity (efficient mixing), Lynden-Bell (1967) used a MEP and
looked for the most probable equilibrium state consistent with the constraints
of the collisionless dynamics. The maximization of the Lynden-Bell entropy $S$
under suitable constraints leads to a coarse-grained distribution function
$\bar{f}({\mbox{\boldmath$r$}},{\mbox{\boldmath$v$}})$ similar to the Fermi-
Dirac distribution function. Therefore, the process of violent relaxation may
provide a justification of the Fermi-Dirac distribution function for DM halos
without the need of efficient gravitational encounters.

However, when coupled to gravity, this distribution function has an infinite
mass (i.e., there is no maximum entropy state), implying that either violent
relaxation is incomplete or that tidal effects have to be taken into account
(if the system is not isolated). The problem therefore becomes an out-of-
equilibrium problem and it is necessary to develop a kinetic theory of
collisionless relaxation (see e.g. Chavanis, 2021 for a review).

One approach is to use a Maximum Entropy Production Principle (MEPP) and argue
that the most probable evolution of the system on the coarse-grained scale is
the one that maximizes the rate of Lynden-Bell entropy $\dot{S}$ under the
constraints of the collisionless dynamics (Chavanis et al., 1996). This leads
to a generalized Fokker-Planck equation having the form of a fermionic Kramers
equation

$\frac{\partial\bar{f}}{\partial
t}+\mbox{\boldmath$v$}\cdot\frac{\partial\bar{f}}{\partial\mbox{\boldmath$r$}}-\nabla\Phi\cdot\frac{\partial\bar{f}}{\partial\mbox{\boldmath$v$}}=\frac{\partial\mbox{\boldmath$J$}}{\partial\mbox{\boldmath$v$}},$
(5)

where
$\mbox{\boldmath$J$}=D[\partial\bar{f}/\partial\mbox{\boldmath$v$}+(mc^{2}/kT)\bar{f}(1-f/\eta_{0})\mbox{\boldmath$v$}]$
is a diffusion current pushing the system towards statistical equilibrium, $D$
the diffusion coefficient, $T\equiv T(t)$ is the temperature evolving in time
so as to conserve the total energy (Chavanis, 1998), $k$ is the Boltzmann
constant, $c$ is the speed of light, and $m$ is the DM fermion mass. However,
this approach is heuristic and does not determine the expression of the
diffusion coefficient.

An alternative, more systematic, approach is to develop a quasilinear theory
of “gentle” collisionless relaxation (Kadomtsev & Pogutse, 1970; Severne &
Luwel, 1980; Chavanis, 1998, 2004) leading to a fermionic Landau equation of
the form

$\frac{\partial\bar{f}}{\partial
t}+\mbox{\boldmath$v$}\cdot\frac{\partial\bar{f}}{\partial\mbox{\boldmath$r$}}-\nabla\Phi\cdot\frac{\partial\bar{f}}{\partial\mbox{\boldmath$v$}}=\frac{8\pi
G^{2}m^{8}\epsilon_{r}^{3}\epsilon_{v}^{3}\ln\Lambda}{h^{6}}\frac{\partial}{\partial
v_{i}}\int{\rm d}\mbox{\boldmath$v$}^{\prime}\\\
\times\frac{u^{2}\delta_{ij}-u_{i}u_{j}}{u^{3}}\left\\{\bar{f}^{\prime}\left(1-\bar{f}^{\prime}\right)\frac{\partial\bar{f}}{\partial
v_{j}}-\bar{f}\left(1-\bar{f}\right)\frac{\partial\bar{f}^{\prime}}{\partial
v^{\prime}_{j}}\right\\},$ (6)

where
$\bar{f}^{\prime}\equiv\bar{f}({\mbox{\boldmath$r$}},{\mbox{\boldmath$v$}}^{\prime},t)$,
$\ln\Lambda=\ln(R/\epsilon_{r})$ is the Coulomb logarithm, $R$ is the typical
size of the system,
$\mbox{\boldmath$u$}=\mbox{\boldmath$v$}^{\prime}-\mbox{\boldmath$v$}$ is the
relative velocity between the “macro-particles” of mass $m_{\rm eff}\sim
2m^{4}\epsilon_{r}^{3}\epsilon_{v}^{3}/h^{3}\gg m$ (see also below), and
$\epsilon_{r}$, $\epsilon_{v}$ are the correlation lengths in position and
velocity respectively. One can make a connection between the above two kinetic
equations by using a form of thermal bath approximation, i.e., by replacing
$\bar{f}^{\prime}$ in Eq. 6 by its equilibrium (Fermi-Dirac) expression. This
substitution transforms an integro-differential (Landau) equation into a
differential (Kramers) equation. In this manner one can compute the diffusion
coefficient explicitly (Chavanis, 1998).

The timescale of violent relaxation is a few $10-100$ dynamical times
($t_{D}$), which is shorter than the Hubble time $t_{H}$. This is confirmed by
the kinetic theory of violent relaxation that predicts a collisionless
relaxation time $t_{R}^{\rm non-coll.}\sim(M/m_{\rm eff})t_{D}$ which is much
shorter than the collisional relaxation time $t_{R}^{\rm coll.}\sim(M/m)t_{D}$
because $m_{\rm eff}\gg m$ (see formula in the above paragraph).

Indeed, the relaxation of the coarse-grained DF
$\bar{f}({\mbox{\boldmath$r$}},{\mbox{\boldmath$v$}},t)$ towards the Lynden-
Bell distribution (of Fermi-Dirac type, see Eq. 7 below) on a few dynamical
times can be interpreted in terms of “collisions” between “macro-particles” or
“clumps” (i.e. correlated regions) with a large effective mass $m_{\rm eff}$
(Kadomtsev & Pogutse, 1970). These macro-particles considerably accelerate the
relaxation of the system (as compared to ordinary gravitational encounters
between particles of mass $m$) by increasing the diffusion coefficient $D$ in
Eq. 5.

Processes of incomplete relaxation could be taken into account by generalizing
the kinetic approach so that the diffusion coefficient rapidly falls off to
zero in space and time, thereby leading to a sort of kinetic blocking.
Alternatively, if the system is submitted to tidal interactions from
neighboring systems one can look for a stationary solution of Eq. 5 which
accounts for the depletion of the distribution function above an escape
energy.

For classical systems evolving through two-body gravitational encounters like
globular clusters, this procedure leads to the King model (King, 1962). For
fermionic DM halos, one obtains the fermionic King model (Ruffini & Stella,
1983; Chavanis, 1998)

$\bar{f}(r,\epsilon\leq\epsilon_{c})=\frac{1-e^{[\epsilon-\epsilon_{c}(r)]/kT(r)}}{e^{[\epsilon-\mu(r)]/kT(r)}+1},\qquad\bar{f}(r,\epsilon>\epsilon_{c})=0\,,$
(7)

which has been written in the case of general relativistic fermionic systems
for the sake of generality (Argüelles et al., 2018, 2022a). Here,
$\epsilon=\sqrt{p^{2}c^{2}+m^{2}c^{4}}-mc^{2}$ is the particle kinetic energy,
$\mu(r)$ is the chemical potential (with the particle rest-energy subtracted
off), $\epsilon_{c}(r)$ is the escape energy (with the particle rest-energy
subtracted off), and $T(r)$ is the effective temperature. The corresponding
set of three dimensionless parameters (for fixed $m$) are defined by the
temperature, degeneracy and cutoff parameters, $\beta(r)=kT(r)/(mc^{2})$,
$\theta(r)=\mu(r)/[kT(r)]$ and $W(r)=\epsilon_{c}(r)/[kT(r)]$, respectively (a
subscript $0$ is used when the parameters are evaluated at the center of the
configuration).

This distribution function takes into account the Pauli exclusion principle as
well as tidal effects, and can lead to a relevant model of fermionic DM halos
usually referred as the RAR model, which has been successfully contrasted
against galaxy observables (Argüelles et al., 2018, 2019; Becerra-Vergara et
al., 2020, 2021; Argüelles et al., 2022b).

The full family of density $\rho(r)$ and pressure $P(r)$ profiles within this
model can be directly obtained as the corresponding integrals of $\bar{f}(p)$
over momentum space (bounded from above by $\epsilon\leq\epsilon_{c}(r)$) as
detailed in Argüelles et al. (2018). This leads to a four-parametric fermionic
equation of state depending on ($\beta_{0},\theta_{0},W_{0},m$) according to
the parameters in Eq. 7. Once with the fermionic distribution function at
equilibrium as obtained from the MEP explained above, we make use of the fact
that a relaxed system of fermions under self-gravity does admit a perfect
fluid approximation (Ruffini & Bonazzola, 1969). Thus, we use the stress-
energy tensor of a perfect fluid in a spherically symmetric metric,
$g_{\mu\nu}=\operatorname{diag}(\mathrm{e}^{2\nu(r)},-\mathrm{e}^{2\lambda(r)},-r^{2},-r^{2}\sin^{2}\vartheta)$
with $\nu(r)$, $\lambda(r)$ being the temporal and spatial metric functions
and $\vartheta$ the azimutal angle. Such configuration leads to hydrostatic
equilibrium equations of self-gravitating fermions. The local $T(r)$, $\mu(r)$
and $W(r)$ fulfill the Tolman, Klein and particle’s energy conservation
relations, respectively (see Argüelles et al., 2018 for details). Further,
$\nu_{0}$ is here constrained by the Schwarzschild condition $g_{00}g_{11}=-1$
at the surface where the halo pressure (and density) falls to zero.

An illustration of a core-halo ($\theta_{0}>10$) and a corresponding halo-only
($\theta_{0}\ll-1$) solution are shown in Fig. 1 (we refer to section 5.2 and
to the cited works above to get a better understanding of the rich morphology
of the fermionic DM model).

Figure 1: Illustrations of DM models: density profile (top) and rotation curve
(bottom). Shown are typical configuration parameters as obtained from this
work and detailed in section B. The core-halo solution of the fermionic DM
model is illustrated by the configuration parameters $\beta_{0}=${10}^{-7}$$,
$\theta_{0}=30$ and $W_{0}=60$, in combination with a corresponding halo-only
solution. For comparison, other DM models (NEW, DC14, Burkert, Einasto) are
added. All profiles are normalized with respect to the halo radius $r_{h}$,
defined at the halo velocity maximum, with $\rho_{h}=r(r_{h})$ and
$v_{h}=v(r_{h})$.

Mass distributions of that fermionic model are well characterized by the
cutoff difference $W(r_{p})-W(r_{s})$. Since $W(r)$ is defined to be zero at
the surface, i.e. $W(r_{s})=0$ with $r_{s}$ being the surface radius where the
density drops to zero, we need to focus only on the plateau cutoff
$W_{p}=W(r_{p})$ with $r_{p}$ being the plateau radius defined at the first
minimum in the rotation curve.

Other important quantities of the core-halo family of fermionic DM mass
solutions are the core mass $M_{c}=M(r_{c})$ with $r_{c}$ being the core
radius defined at the first maximum in the rotation curve, the halo mass
$M_{h}=M(r_{h})$ with $r_{h}$ being the halo radius defined at the second
maximum in the rotation curve, and the total mass $M_{s}=M(r_{s})$ given at
the surface radius $r_{s}$.

The density profiles in the fermionic model can develop a rich morphological
behaviour: while the halo region is King-like (i.e. from polytropic-like with
$W_{p}\ll 1$ to power law-like with $W_{p}\gtrsim 10$, see 5.2), the inner
region can either develop a dense core at the center of such a halo (i.e. for
large central degeneracy $\theta_{0}>10$), or not (i.e. $\theta_{0}\ll-1$ in
the dilute regime). Both kind of family profiles are thermodynamically and
dynamically stable as well as long lived in a cosmological framework, as
recently demonstrated in Argüelles et al. (2021) for typical galaxies with
total masses of the order
$\sim$5\text{\times}{10}^{10}\text{\,}\mathrm{M_{\odot}}$$.

Remarkably, for core \- halo RAR solutions with fermion masses of
$mc^{2}\approx$50\text{\,}\mathrm{keV}$$, the degenerate and compact DM cores
may work as an alternative to the BH paradigm at the center of non-active
galaxies (Argüelles et al., 2018, 2019; Becerra-Vergara et al., 2020, 2021;
Argüelles et al., 2022b). Furthermore, their eventual gravitational core-
collapse in larger galaxies may offer a novel supermassive BH formation
mechanism from DM (Argüelles et al., 2021).

###  3.2 Other DM halo models

From a phenomenological viewpoint, a DM halo density profile is described by
three characteristics: the inner halo, the outer halo and the transition in
between. Such density profiles are usually described by the
($\alpha,\beta,\gamma$)-model (Hernquist, 1990)

$\frac{\rho(r)}{\rho_{\mathrm{N}}}=\left[\vphantom{}\frac{r}{R_{\mathrm{N}}}\right]^{-\gamma}\left[\vphantom{}1+\left[\vphantom{}\frac{r}{R_{\mathrm{N}}}\right]^{\alpha}\right]^{-\frac{\beta-\gamma}{\alpha}},$
(8)

where $\alpha$ describes the transition, $\beta$ the outer slope and $\gamma$
the inner slope. Following Newtonian dynamics — that is fully sufficient on
halo scales — then the velocity is given by

$\frac{v^{2}(r)}{\sigma_{\mathrm{N}}^{2}}=\frac{R_{\rm
N}}{r}\frac{M(r)}{M_{\rm N}}$ (9)

and the enclosed mass by

$\frac{M(r)}{M_{\rm N}}=\int\limits_{0}^{r}\left[\vphantom{}\frac{r}{R_{\rm
N}}\right]^{2}\frac{\rho(r)}{\rho_{\mathrm{N}}}\frac{\mathrm{d}r}{R_{\rm N}},$
(10)

For the following DM halo models we will use $R_{\mathrm{N}}$,
$\rho_{\mathrm{N}}$, $\sigma_{\mathrm{N}}^{2}=GM_{\mathrm{N}}/R_{\mathrm{N}}$
and $M_{\mathrm{N}}=4\pi\rho_{\mathrm{N}}R_{\mathrm{N}}^{3}$ as scaling
factors for length, density, velocity and mass, respectively.

In Fig. 1 we illustrate the typical morphology of common DM halo models used
here for the SPARC galaxies. For a better comparison the plots are normalized
with respect to the halo located at the velocity maximum on halo scales.

Based on the general ($\alpha,\beta,\gamma$)-model, Di Cintio et al. (2014)
modelled CDM halos including baryonic feedback mechanisms in galaxy formation.
They found that the three parameters ($\alpha,\beta,\gamma$) are related
through

$\displaystyle\alpha$
$\displaystyle=2.94-\log_{10}\left[\vphantom{}(10^{X+2.33})^{-1.08}+(10^{X+2.33})^{2.29}\right],$
(11) $\displaystyle\beta$ $\displaystyle=4.23+1.34X+0.26X^{2},$ (12)
$\displaystyle\gamma$
$\displaystyle=-0.06+\log_{10}\left[\vphantom{}(10^{X+2.56})^{-0.68}+10^{X+2.56}\right],$
(13)

where $X=\log_{10}(M_{*}/M_{\mathrm{halo}})$ describes the stellar-to-dark
matter ratio. For the circular velocity Eq. 9 the enclosed mass Eq. 10 is
given by a hypergeometric function

$\frac{M(r)}{M_{\mathrm{N}}}=\frac{1}{3-\gamma}\left[\vphantom{}\frac{r}{R_{\mathrm{N}}}\right]^{3-\gamma}{}_{2}F_{1}(p_{1},p_{2};\,q_{1};\,-[r/R_{\mathrm{N}}]^{\alpha}),$
(14)

with $p_{1}=(3-\gamma)/\alpha$, $p_{2}=(\beta-\gamma)/\alpha$ and
$q_{1}=1+(3-\gamma)/\alpha$. In the following, we refer this model as DC14.

Alternatively, for $\alpha=1$, $\beta=3$ and $\gamma=1$ the
($\alpha,\beta,\gamma$)-model reduces to the NFW model (Navarro et al., 1996,
1997) as obtained from early DM-only N-body simulations. This DM model
develops cuspy halos of the following type

$\frac{\rho(r)}{\rho_{\mathrm{N}}}=\left[\vphantom{}\frac{r}{R_{\mathrm{N}}}\right]^{-1}\left[\vphantom{}1+\frac{r}{R_{\mathrm{N}}}\right]^{-2},$
(15)

with the circular velocity

$\frac{v^{2}(r)}{\sigma_{\mathrm{N}}^{2}}=\frac{\ln(1+r/R_{\mathrm{N}})}{r/R_{\mathrm{N}}}-\frac{1}{1+r/R_{\mathrm{N}}}.$
(16)

In contrast to NFW, Burkert (1995) proposed a DM density profile with a cored
halo of the following type

$\frac{\rho(r)}{\rho_{\mathrm{N}}}=\left[\vphantom{}1+\frac{r}{R_{\mathrm{N}}}\right]^{-1}\left[\vphantom{}1+\left(\frac{r}{R_{\mathrm{N}}}\right)^{2}\right]^{-1}.$
(17)

For the circular velocity Eq. 9 the enclosed mass Eq. 10 is given by

$\frac{M(r)}{M_{\rm N}}=\frac{1}{4}\ln(1+[r/R_{\rm
N}]^{2})+\frac{1}{2}\ln(1+r/R_{\rm N})\\\ -\frac{1}{2}\arctan(r/R_{\rm N}).$
(18)

With $M_{0}\approx M(R_{N})$ being the mass scale originally interpreted as
the core mass of the halo (e.g. Salucci & Burkert, 2000) we obtain the
relation $M_{N}=8M_{0}$. Further, the density scale $\rho_{N}$ describes the
central density $\rho_{0}$ and the length scale $R_{\rm N}$ can be identified
with the Burkert radius $r_{\rm B}$ fulfilling the condition $\rho(r_{\rm
B})=\rho_{0}/4$.

Another interesting and successful candidate is the Einasto model (Einasto &
Haud, 1989), a purely empirical fitting function with no commonly recognized
physical basis (Merritt et al., 2006). The DM halo density profiles of that
model are of the following type, given in a normalized form,

$\frac{\rho(r)}{\rho_{\mathrm{N}}}=\mathrm{e}^{-[r/R_{\rm N}]^{\kappa}}.$ (19)

The exponent $\kappa$ describes the shape of the density profile. The circular
velocity and the enclosed mass are given by Eqs. 9 and 10. This model develops
mass distributions with a finite mass $M_{\rm tot}/M_{\rm
N}=\Gamma(3/\kappa)/\kappa$ for $r\to\infty$ (see also Retana-Montenegro et
al., 2012). The typical $\kappa$ values obtained in this work for the SPARC
data-set, as well as the comparison with the same values coming from N-body
simulations (either with or without baryonic effects), are given in
subsections 3.3, 3.4 and 4.2, and Fig. 13.

###  3.3 Fitting priors and Monte-Carlo approach

For the fermionic DM model we fix the particle mass $m$ and therefore reduce
the number of free parameters by one, e.g.
$\mbox{\boldmath$p$}=(\beta_{0},\theta_{0},W_{0})$. A particle mass of
$mc^{2}=$50\text{\,}\mathrm{keV}$$ is well motivated by the promising results
obtained in Argüelles et al. (2018); Becerra-Vergara et al. (2020, 2021);
Argüelles et al. (2022b), where the fermionic core-halo DM profile was able to
explain both the S-stars orbits around SgrA*, and the Milky Way rotation
curve. In Argüelles et al. (2019, 2021), in particular, and for the same
particle mass, the fermionic core-halo profiles were successfully applied to
other galaxy types from dwarf to larger galaxy types, providing a possible
explanation for the nature of the intermediate mass BHs, as well as a possible
mechanism for SMBH formation in active galaxies. Moreover, as demonstrated in
Argüelles et al. (2018) there exists a particle mass range between $\sim$ 50
and $\sim$ 350 keV where the compacity of the central core can increase (all
the way to its critical value of collapse) while maintaining the same DM halo-
shape. Therefore, regarding the SPARC RC fitting, as well as all the scaling
relations on halos-scales, our conclusions are not biased by the choice of the
particle mass in the above range.

For the fermionic model we consider solutions which are either in the dilute
regime ($\theta_{0}\ll-1$, i.e. are King-like), or which have developed a
degenerate core (i.e. $\theta_{0}>10$) at the center of such a halo. Fermionic
solutions within only these two families have been shown to be
thermodynamically and dynamically stable when applied to galaxies (Argüelles
et al., 2021).

The NFW and the Burkert models are described by two free scaling parameters,
e.g. $\mbox{\boldmath$p$}=(R_{\rm N},\rho_{\rm N})$. The DC14 model with e.g.
$\mbox{\boldmath$p$}=(X,R_{\rm N},\rho_{\rm N})$ and the Einasto model with
e.g. $\mbox{\boldmath$p$}=(\kappa,R_{\rm N},\rho_{\rm N})$ are described by
three free parameters. Compared to NFW and Burkert, both (Einasto and DC14)
have an additional parameter which affects the sharpness of the transition
from the inner to the outer halo.

In order to find the best-fits we use the LM algorithm (see section 2.2) with
well chosen initial values (i.e. priors) reflecting astrophysical (realistic)
scenarios. Because the LM algorithm finds only local minima, we choose the
parameter sets randomly within a range and follow a Monte-Carlo approach. For
the fermionic core-halo solutions, we choose $\beta_{0}\in[10^{-8},10^{-5}]$,
$\theta_{0}\in[25,45]$ and $W_{0}\in[40,200]$ which correspond to a
conservatively wide range of parameters according to Argüelles et al. (2019).
For the fermionic diluted solutions, we cover the same range for $\beta_{0}$
and $W_{0}$, but with $\theta_{0}\equiv\theta_{p}\in[-40,-20]$. For the other
DM models, the initial scaling factors are chosen from $R_{\rm
N}\in[10^{1},10^{4}]$\mathrm{pc}$$ and $\rho_{\rm
N}\in[10^{-4},10^{-1}]$\mathrm{M_{\odot}}\mathrm{/}\mathrm{pc}^{3}$$.
According to Katz et al. (2017) the additional parameter of the DC14 model is
chosen from $X\in[-3.75,-1.3]$. For the Einasto model, the exponent is chosen
from $\kappa\in[0.1,10]$. This relatively large window has been chosen in
order (i) to account for the broad diversity of rotation curves covered by the
SPARC data, and (ii) not to be limited by any fixed value (e.g. as typically
obtained by CDM-only simulations) since this is an independent analysis to
that of N-body simulations, and may also account for other effects such as
baryonic feedback (see also the discussion in next subsection).

###  3.4 Baryonic feedback

The DM in galaxies evolves together with baryons, thus it is expected some
degree of baryonic feedback on scales where baryons dominate. This also seems
to be the case within the SPARC galaxies here considered, where almost half of
the points in upper-left panel of Fig. 2 roughly fulfill
$a_{\mathrm{tot}}<2a_{\mathrm{bar}}$. One of the main baryonic effects onto
the total gravitational potential is thought to happen due to a sustained
process of stellar bursts, which drives baryonic material from inner-halo
regions while ending in a reduction of DM densities at those halo scales (see
e.g. Governato et al., 2012; Read et al., 2016). However, the quantification
of such baryonic effects have always been calibrated and applied to cuspy CDM
halos. Only recently it was shown that baryonic feedback is DM-model
dependent.

That is, in WDM cosmologies such effects are typically diminished with respect
to CDM (Bozek et al., 2019). The main reasons behind this attenuation are that
in WDM cosmologies DM halos form later, are less centrally dense on inner-halo
scales, and therefore contain galaxies that are less massive with less baryon
content than the CDM counterparts. Since the halos of the Burkert, Einasto and
fermionic DM model correspond to WDM cosmologies, it is expected that baryonic
effects are milder in those cases. However, a thorough quantification of this
feedback has only been worked out for Einasto profiles (Bozek et al., 2019)
and still remains to be solved in the case of fermionic halos before taking
any conclusion, though out of the scope of this work.

In any case, and regardless of a potential feedback of the baryons onto the
fermionic DM profile, it is important to emphasize a key result of this work:
the flatness in the inner-halo slope of MEP profiles is due to a (quasi)
thermodynamic equilibrium reached by the DM particles, thus involving a
different physical principle than the one explained above for the baryons.

Finally, two DM model considered here — DC14 (Di Cintio et al., 2014) and
Einasto (Bozek et al., 2019) — (effectively) account for baryonic feedback.
That is, any baryonic feedback expected to arise for SPARC galaxies should be
reflected in the best-fit parameters of the DM profiles with its consequent
universal shape reflecting such effects. In section 4.2 we compare our
statistical results regarding baryonic effects in the DM model free-parameters
with the ones reported in the literature.

Figure 2: Radial acceleration relation (top) and mass discrepancy acceleration
relation (bottom) for SPARC data and competing DM halo models. Each plot is
divided in 50x50 equal bins. The baryonic centripetal acceleration
$a_{\mathrm{bar}}$ is inferred from luminosity observables while the total
acceleration $a_{\mathrm{tot}}$ is inferred independently from velocity
fields. For DM halo models the total acceleration is composed of the predicted
dark and inferred baryonic components, i.e.
$a_{\mathrm{tot}}=a_{\mathrm{bar}}+a_{\mathrm{DM}}$. The corresponding solid
curves are the best fits characterized by a specific $\mathfrak{a}_{0}$. The
histogram plots (upper row) show a Gaussian distribution of
$\log_{10}(a_{\mathrm{tot}}/\mathfrak{a}_{0})$. The grayscale legend shows the
number of points per bin (of a total 2396 for the 120 SPARC-galaxies used).
Figure 3: continue

##  4 Results

After an insightful analysis of the parameter distribution of each DM model
which best fits the SPARC RCs (given in appendix B), we compare between them
following two complementary approaches.

First, we consider the entire galaxy sample and extract the radial
acceleration information for the total and baryonic components at each
galactocentric radii, and put them all together as in Fig. 2 (SPARC-window).
In this way, we reduce any characteristics of individual galaxies into an
overall (global) picture, i.e. the Radial Acceleration Relation (McGaugh et
al., 2016) and Mass Discrepancy Acceleration Relation (McGaugh, 2004, 2014),
respectively.

In the second approach, we consider each individual galaxy and perform a
goodness of DM-model analysis, showing how well the inferred DM RCs are fitted
by these models.

###  4.1 Acceleration relations

The rotation curve of each component (e.g. bulge, disk, gas) traces its
centripetal acceleration $a=v^{2}/r$, giving access to independent
acceleration measurements. The Radial Acceleration Relation compares the
radial acceleration due to the total mass ($a_{\mathrm{tot}}$) with that due
to the baryonic mass ($a_{\mathrm{bar}}$). It is empirically described by
McGaugh et al. (2016)

$a_{\mathrm{tot}}=\frac{a_{\mathrm{bar}}}{1-\mathrm{e}^{-\sqrt{a_{\mathrm{bar}}/{\mathfrak{a}_{0}}}}},$
(20)

where $\mathfrak{a}_{0}$ is the only adjustable parameter. In the low
acceleration regime ($a_{\mathrm{bar}}\ll\mathfrak{a}_{0}$), where DM
dominates, it clearly shows a deviation from a linear correlation (see top
panels of Fig. 2). While in the high acceleration regime
($a_{\mathrm{bar}}\gg\mathfrak{a}_{0}$), dominated by baryonic matter, the
linear relation is recovered.

This relation is not limited to disk galaxies but also holds for other galaxy
types (e.g. ellipticals, lenticulars, dwarfs spheroidals and even low-surface-
brightness galaxies), what makes it a true universal law among morphology
classification (McGaugh et al., 2016; Lelli et al., 2017; Di Paolo et al.,
2019).

For the SPARC galaxies the different accelerations $a_{\mathrm{bar}}=V_{\rm
bar}^{2}/r$ and $a_{\mathrm{tot}}=V_{\rm tot}^{2}/r$ are inferred from the
SPARC data. See section 2.1 for details how the circular velocities and radii
are obtained. For the competing DM models $a_{\mathrm{bar}}=V_{\rm bar}^{2}/r$
is equally inferred from the SPARC data while $a_{\mathrm{tot}}=[V_{\rm
bar}^{2}+v_{\rm DM}(r)^{2}]/r$ is inferred from the total mass distribution,
composed of the observed baryonic component (taken from SPARC data) and the
best-fitted circular velocities of the dark matter component $v_{\rm DM}(r)$
for each DM model. See sections 2 and 3.3 for details how best-fits are
obtained.

For the SPARC data as well as for each DM model, we applied then a least-
square fitting to obtain $\mathfrak{a}_{0}$. The result obtained here from the
SPARC data only (i.e. without assuming any specific DM model) is fully
consistent (within errors) with
$\mathfrak{a}_{0}\approx$1.2\text{\times}{10}^{-10}\text{\,}\mathrm{m}\mathrm{/}\mathrm{s}^{2}$$
as obtained originally in McGaugh et al. (2016), validating our procedure.
Each DM model is then characterized by a specific best-fit $\mathfrak{a}_{0}$.
The corresponding curves are plotted as solid lines in Fig. 2. In all cases,
we obtain values close to the one obtained in McGaugh et al. (2016). This
allows to conclude that all competing DM models are able to reproduce the
Radial Acceleration Relation. Moreover, they reproduce it equally good without
a clear statistically preferred model.

A closely related relation is the MDAR relation between the baryonic and total
mass components, defined by $D=M_{\rm tot}/M_{\rm bar}$ with $M_{\rm bar}$
being the total baryonic mass and $M_{\rm tot}$ the total galaxy mass
accounting for baryonic and dark matter. For
$a=\frac{\mathrm{d}\Phi}{\mathrm{d}r}$ and $\Phi(r)$ being the gravitational
potential of a spherically symmetry mass distribution the mass discrepancy can
be equivalently written as $D=a_{\mathrm{tot}}/a_{\mathrm{bar}}$. The results
are illustrated in the bottom panels of Fig. 2.

According to some authors, the above two acceleration relations do not imply
the need of any new physics and may be explained within the $\Lambda$CDM
framework (Di Cintio & Lelli, 2016; Katz et al., 2017; Salucci, 2016, 2018).
For smaller disk and LSB galaxies (extending the original SPARC sample), Di
Paolo et al. (2019) found that Eq. 20 is a limiting case of a more complex
relation with the need of adding one extra galaxy parameter. In addition,
based on modern cosmological simulations, Keller & Wadsley (2017) predict even
a redshift dependency of the acceleration parameter $\mathfrak{a}_{0}$,
emphasizing that the correlation is universal only regarding the morphological
classification.

Additionally to the acceleration relations mentioned above (which account for
the entire accelerations distributions among all galaxies and at different
radii), the next strategy is to focus on the DM components of each galaxy and
gather best-fits of the inferred rotation curves allowing for another
(related) quantitative comparison of the DM models.

###  4.2 Goodness of model

The SPARC galaxies show different characteristics in their rotation curve such
as a nearly flat curve through the entire galaxy data; a rising trend in the
inner halo followed by a single maximum; or multiple extrema in the form of
oscillations. See Fig. 8 for three typical examples within the SPARC data-set.
Some galaxies show just a rising trend implying that the rotation curves are
incomplete, likely due to the faintness and/or lack of data for outermost halo
stars. Of interest is therefore a quantitative description about the goodness
of a DM halo model fitting the entire galaxy sample (120 galaxies).

The goodness of a fit for a single galaxy is well described by the $\chi^{2}$
value, see Eq. 4. When competing models with different number of parameters
are compared it is appropriate to consider the reduced $\chi^{2}$ defined as
$\chi_{r}^{2}=\chi^{2}/d$ with the degree of freedom $d=N-p$, $N$ being the
number of observables (for a single galaxy) and $p$ the number of parameters
(of the considered model).

The question now arises how to compare the competing models for a population
of galaxies. In order to find the goodness of a model which is robust against
outliers, we ask how many fitted galaxies have a (reduced) $\chi^{2}$ lower
than a given one. It turns out that the population curve resembling a
cumulative distribution function (CDF) follows nearly a log-normal
distribution. We use the mean value, labelled as $\hat{\chi}^{2}_{r}$, as the
criteria to described the goodness of a model for fitting a population of
galaxies within the SPARC data-set. A parameter analysis of the best-fit
solutions for each DM halo model is detailed in appendix B.

Figure 4: Goodness of model analysis for the entire sample (120 galaxies). We
count the population of fitted galaxies having a reduced $\chi^{2}$ smaller
than a given one. The normalized population (step-like) follows nearly a log-
normal distribution (solid) characterized by the mean value
$\hat{\chi}^{2}_{r}$. The shaded regions span the $95\text{\,}\mathrm{\char
37\relax}$ confidence interval of the best-fitted $\hat{\chi}^{2}_{r}$. Figure
5: Goodness of model analysis for galaxies showing one clear maximum in the
outer halo RC, where DM typically dominates. This condition is fulfilled by 44
galaxies (see e.g. UGC05986 in Fig. 8). We count the population of fitted
galaxies having a reduced $\chi^{2}$ smaller than a given one. The normalized
population (step-like) follows nearly a log-normal distribution (solid)
characterized by the mean value $\hat{\chi}^{2}_{r}$. The shaded regions span
the $95\text{\,}\mathrm{\char 37\relax}$ confidence interval of the best-
fitted $\hat{\chi}^{2}_{r}$.

The goodness analysis for the entire SPARC sample (120 galaxies) is shown in
Fig. 4 and can be summarized as: The NFW model ($p=2$) is statistically
disfavoured with respect to the other DM halo models; the Einasto model
($p=3$) and the DC14 model ($p=3$), on the other hand, are statistically
favoured with similarly good results; the Burkert model ($p=2$) and the
fermionic model ($p=3$) produce comparable results but are somewhat
statistically less favoured.

We recall that the DM models are suited for DM halo rotation curves with only
one maximum. Therefore, we restrict further the SPARC sample to galaxies (44)
showing one clear maximum in their rotation curve (see also section 5.3 for
details about the diversity of SPARC galaxies). This conditions is
equivalently reflected in the fermionic scenario where data supports for a
significant escape of DM particles, that is for energy-cutoff values at
plateau of $W_{p}<10$ (see appendix B for further details). For that sub-
sample the picture changes considerably, see Fig. 5, and can be summarized as:
The NFW model is statistically even more disfavoured while all other models
(including fermionic DM) become more comparable — with a little tendency for
Einasto and against Burkert.

We remind that the DC14 model is based on the analysis of hydro-dynamically
simulated galaxies including complex baryonic feedback processes (Di Cintio et
al., 2014). The good performance of DC14 may indicate the importance of
baryonic feedback in galaxy formation. Indeed, our results regarding DC14 are
in line with the literature as well since, as can be seen from Fig. 12 (bottom
panel), the bulk of our ($\alpha,\beta,\gamma$) parameters for the SPARC data-
set lies within the windows (0,2.6); (2.3,4); (0,2) in rough agreement with Di
Cintio et al. (2014). Interestingly, we obtain a mean for the stellar to DM
mass ratio $X=\log_{10}(M_{*}/M_{\mathrm{halo}})$ of $-2.4$, which is within
the $X$ values ($-2.5,-2.3$) where DM cusps are most effectively flattened due
to baryonic effects as reported in Di Cintio et al. (2014). Besides, our
results are in good agreement with those reported in Behroozi et al. (2013)
from abundance matching, since they obtain X values between (-2.7, -1.7) for
halo masses in the range between $10^{10}M_{\odot}$ and $10^{12}M_{\odot}$,
with $X=-2.4$ for $M_{h}=10^{11}M_{\odot}$, the latter coinciding with the
mean $X$ value in this work.

In the case of Einasto profiles, the comparison with the literature is more
difficult since the effects of baryonic feedback are not explicitly reported
through its free parameters. For example in Bozek et al. (2019) they obtain,
within (zoom-in) hydro-dynamical simulations and for
$\sim${10}^{10}\text{\,}\mathrm{M_{\odot}}$$ halos, a reduction on inner-halo
densities (i.e. at the convergence radius $r=$0.26\text{\,}\mathrm{kpc}$$) of
up to $45\text{\,}\mathrm{\char 37\relax}$ in WDM cosmologies with respect to
the analogous WDM-only simulations. A reduction which is about $15\%$ more
pronounced than comparing WDM to CDM only simulations (Bozek et al., 2019).
Even if all the resulting DM halos are there fitted with the Einasto profile,
these density reductions are only expressed through its concentration
parameter (which is a function of the scale radius and the virial radius).
They found indeed that Einasto profiles in WDM cosmologies have smaller
concentration parameters than the CDM counterparts. This concentration
parameter correlates inversely proportional to the Einasto shape parameter
(denoted here as $\kappa$) as originally shown through Figs. 13 and 14 in
Dutton & Macciò (2014). Thus, considering that $\kappa\approx 0.2$ correspond
to CDM-only cosmologies (Dutton & Macciò, 2014; Fitts et al., 2017), then
larger $\kappa$ values are expected for Einasto profiles having flatter inner-
halo slopes (Dutton & Macciò, 2014). Indeed, our results are qualitatively in
line with those reported in Bozek et al. (2019) since we find a mean value of
$\kappa\approx 0.4$ (see Fig. 13), thus implying smaller concentration
parameters than for CDM profiles. However more work from hydro-dynamical
simulations using Einasto fitting profiles is needed in order to make a proper
quantitative comparison with our statistical results.

We would like also to point out that many galaxies in the SPARC data set are
missing significant information in the outer halo (e.g. due to faint stars) or
show a complex behavior (oscillatory-pattern) in their rotation curves (see
e.g. right box in Fig. 8). In any case, it does not allow to univocally
determine the cutoff parameter (i.e. $W_{p}$) for the fermionic DM model since
any sufficiently large $W_{p}$ would not change the $\chi^{2}$ value (see e.g.
bottom left panel of Fig. 9). Nevertheless, there are some other individual
galaxies where the escaping particles effects are clearly preferred.
Interestingly, all of those galaxies are of magellanic type: NGC0055 (Sm),
UGC05986 (Sm), UGC05750 (Sdm), UGC05005 (Im), F565-V2, (Im), UGC06399 (Sm),
UGC10310 (Sm), UGC07559 (Im), UGC07690 (Im), UGC05918 (Im) and UGC05414 (Im).

Moreover, many galaxies, which are poorly fitted by any of the considered
models, show short range oscillations in their rotation curves with more than
one maximum. None of the models can provide a clear explanation of that
phenomena, found usually in non-magellanic galaxy types: e.g. NGC2403 (Scd),
UGC02953 (Sab), NGC6015 (Scd), UGC09133 (Sab), UGC06787 (Sab), UGC11914 (Sab),
NGC1003 (Scd), NGC0247 (Sd), UGC08699 (Sab) and UGC03205 (Sab).

On phenomenological grounds, in the fermionic DM model it is possible to vary
the width of the maximum bump in the RC through the cutoff parameter in the
strong or moderate cutoff regime ($W_{p}\lesssim 10$). Whether with weak
($W_{p}\gtrsim 10$) or even without cutoff-effects, the RC solutions of the
model show long range oscillations, similar to the isothermal model. In any
case, these RC oscillations have a too long wavelength and therefore do not
offer a convenient explanation. On the other hand, in the case of strong
cutoff, we obtain a narrow maximum bump necessary for many RCs, especially for
galaxies of magellanic type (see above for examples), which usually do not
show those oscillations, but also for some non-magellanic galaxy types, e.g.
NGC5585 (Sd), NGC7793 (Sd), UGC06614 (Sa), ESO079-G014 (Sbc), F571-8 (Sc),
NGC0891 (Sb), UGC06614 (Sa), UGC09037 (Scd), NGC4217 (Sb), UGC04278 (Sd).

NFW and Burkert models cannot explain variations of the inner and outer
rotation curve because the parameters ($\beta$ and $\gamma$) responsible for
such a behaviour (see Eq. 8) are fixed. Additionally a transition from the
inner to the outer halo is generally characterized by $\alpha$ (or $\kappa$).
In contrast to NFW and Burkert, the DC14 and Einasto models have a free
parameter which affects the inner/outer rotation curve steepness and the
sharpness of the halo transition, simultaneously. Such a flexibility is
reflected in generally better $\chi^{2}$ values. Nevertheless, the goodness
for oscillating RCs remains rather poor.

##  5 Fermionic halos

This section is mainly devoted to fermionic halos where we include the
analysis of the DM surface density relation (SDR), originally based on the
Burkert model. Further, we highlight the rich morphology of fermionic DM halos
and show that particular solutions can be associated with different DM models.
As typical examples, we present detailed rotation curve fits and their
corresponding analysis for three selected galaxies indicating the limitations
of the observational data and/or the case where data supports for one clear
maximum in the rotation curve of each DM model.

###  5.1 DM surface density relation

The constant surface density (Donato et al., 2009)

$\Sigma_{0D}=\rho_{\rm 0D}r_{0}\approx
140_{-50}^{+80}$\mathrm{M_{\odot}}\mathrm{/}\mathrm{pc}^{2}$.$ (21)

is valid for about 14 orders of magnitude in absolute magnitude ($M_{B}$)
where $\rho_{0}$ is the central DM halo density and $r_{0}$ the one-halo-
scale-length — both of the Burkert model. At $r_{0}$ the density falls to one-
forth of the central density, i.e. $\rho(r_{0})=\rho_{\rm 0D}/4$.

Note that the center in the Burkert model corresponds to the plateau in the
fermionic DM model, i.e. $\rho_{\rm 0D}\approx\rho_{\rm p}$ where the plateau
density $\rho_{\rm p}$ is defined at the first minimum in the RC. Following
the definition of the Burkert radius $r_{0}$, we identify the one-halo-scale-
length $r_{B}$ of the fermionic model such that $\rho(r_{B})=\rho_{p}/4$. We
thus calculate the product $\rho_{p}r_{B}$ for each galaxy.

The absolute magnitude $M_{B}$ was taken from the Carnegie-Irvine Galaxy
Survey (Ho et al., 2011), providing nine overlapping galaxies with the SPARC
sample. These candidates are well in agreement with the DM surface density
observations, see Fig. 6. The shown candidates include isothermal-like (blue
outlined points) and non-isothermal (green circles) solutions. For comparison,
the results are amended by the MW solution following the fermionic RAR model
(Argüelles et al., 2018).

Figure 6: DM surface density ($\Sigma_{0D}$) predictions in the fermionic
model using the SPARC data-set. The blue region indicates the delimited area
by the $3\sigma$ error bars of all the data points in Donato et al. (2009).
The violet diamond represents the MW as analyzed in Argüelles et al. (2018).
The absolute magnitude was taken from the Carnegie-Irvine Galaxy Survey (Ho et
al., 2011), providing nine overlapping galaxies (blue and green circles). The
full sample (gray bars) and sub-sample (green bars) including only for non-
isothermal solutions ($W_{p}<10$), are shown as histograms, both following
approximately a Gaussian distribution.

Although absolute magnitude information is incomplete, all of the predicted DM
surface densities are within the range of the $3\sigma$ area as well. This is
visualised by a histogram for the full sample (dark grey bars, 120 galaxies)
with comparison to the sub-sample (green bars, 44 galaxies) including non-
isothermal solutions (i.e. $W_{p}<10$). Considering the sub-sample only, we
obtain a mean surface density of about
$148.5\text{\,}\mathrm{M_{\odot}}\mathrm{/}\mathrm{pc}^{2}$, fully inside the
1-$\sigma$ uncertainty in Eq. 21.

The SDR given by Eq. 21 is qualitatively consistent with the scaling relation
$M_{h}\sim r_{h}^{2}$ as given by Eq. A1. Due to the different halo profiles
of the fermionic model, ranging from polytropes of $n=5/2$ to isothermal (see
section 5.2), there is a non-linear relation between the halo radius $r_{h}$
and the one-halo-scale-length $r_{B}$. Nevertheless, considering that the halo
is nearly homogeneous up to approximately the halo radius $r_{h}$ (i.e.
$M_{h}\sim\rho_{p}r_{h}^{3}$), and that $r_{h}\sim r_{B}$ we obtain
$\rho_{p}r_{B}\approx\mathrm{const}$. (see also Argüelles et al., 2019).

###  5.2 Morphology of fermionic halos

We are interested in providing approximate analytic (and semi-analytic)
expressions for the halo morphology associated with the corresponding
fermionic solutions.

On halo scales, i.e. for $r\sim r_{h}$ with $r_{h}$ being the halo radius, the
fermionic DM model resembles the King model (King, 1966). These DM profiles
are characterized by a cored inner halo (i.e. a flat inner slope) followed by
a transition towards a finite mass. When applied to galactic halos, such a
transition can show different behaviours leading to a rich morphology of
density profiles, from polytropic-like to power law-like, mainly depending on
the value of $W_{p}$. Consequently, this cutoff (or particle escape) parameter
controls the sharpness of the inner-outer halo transition as well.

Very similar characteristics are also produced by the Einasto model, although
with a wider spectrum for the sharpness of the transition described by the
parameter $\kappa$.

In the fermionic model, the sharpness of the transition can be quantified as
follows: for $W_{p}\ll 1$, the fermionic density profiles are polytropes of
index $n=5/2$ as clearly shown in Fig. 7 for bluish solutions (see Chavanis et
al., 2015a for a derivation of the polytropic $n=5/2$ equation and its link
with the fermionic solutions, further justifying our results). In contrast,
for negligible escape of particles (e.g. $W_{p}\gtrsim 10$), the fermionic
halos become isothermal-like, with $\rho(r)\propto r^{-2}$ down to the virial
radius as evidenced through the reddish solutions in the same figure. Finally,
for $W_{p}$ values in between the above limiting cases, the fermionic DM halo
profiles start to develop a power-law trend almost matching the Burkert
profile for $W_{p}\approx 7$ (see Fig. 7).

These results agree with those obtained in Sec. VII of Chavanis et al. (2015b)
for the classical King model, which correspond to our fermionic model in the
dilute regime (see also paragraph below). There it is shown that the Burkert
profile gives a good fit of such King distribution for a value of $W_{p}$
close to the point of marginal thermodynamical stability, the latter given at
$W_{p}^{(c)}=7.45$. In Chavanis et al. (2015b) this marginal solution is
labelled through the equivalent parameter $k$.

Statistically, we identified $W_{p}\approx 10$ as a discriminator between two
groups and further detailed in Appendix B. Even if the precise value might be
biased by the data, there is an interesting underlying physical explanation
for it: fermionic DM profiles as obtained from a MEP are thermodynamically and
dynamically stable for $W_{p}\leq 7.45$ if they are in the dilute regime
($\theta_{0}\ll-1$ where $W_{0}\equiv W_{p}$), while in the core-halo regime
($\theta_{0}>10$) the same stability conditions hold for a bounded
$W_{p}\in[W_{p}^{\rm min},W_{p}^{\rm max}]$. For typical SPARC galaxies with a
halo mass of $\sim$5\text{\times}{10}^{10}\text{\,}\mathrm{M_{\odot}}$$, we
find $W_{p}^{\rm max}\sim${10}^{-2}$$.

Figure 7: Morphology of best-fitted fermionic halos (solid). For comparison
the Burkert profile (dashed) and a polytrope of $n=5/2$ (dot-dashed) are
added. Shown are normalized halo density profiles where the density is
normalized by the inner halo density $\rho_{p}$ and the radius by the Burkert
radius $r_{B}$ fulfilling $\rho(r_{B})=\rho_{p}/4$.

The point of stability-change for the classical King model was first obtained
in Chavanis et al. (2015b), and further re-derived here for fermionic DM
profiles but this time for realistic average halo-mass galaxies, as obtained
from the SPARC data-set, following the thermodynamic analysis of Argüelles et
al. (2021). That is, typical galaxies with halo masses of about
$\sim$5\text{\times}{10}^{10}\text{\,}\mathrm{M_{\odot}}$$ and with
appreciable escape of particles (i.e. $W_{p}\ll 1$) are thermodynamically and
dynamically stable, suggesting a deep link between thermodynamics of self-
gravitating systems and galaxy formation. A case-by-case stability analysis in
relation to observed galaxies (e.g. SPARC data-set) is out of the scope of the
present paper and will be the subject of a future work.

Figure 8: Total rotation curves and their composition, shown for UGC05986
(left), DDO161 (center) and NGC6015 (right). The thick lines are the best-fits
of three competing DM models: Fermionic (blue), DC14 (yellow) and NFW (red).
For clarity we have excluded Einasto.

###  5.3 Diversity of SPARC rotation curves

We show in this section a detailed $\chi^{2}$ analysis of the RC fits for
three selected galaxies, each representing some characteristics of given
observational data. We divide the SPARC galaxies in three groups by the
inferred DM component as explained next. This analysis is based on the
fermionic model where such a grouping seems to be appropriate to select
galaxies with valuable predictions about the inner halo.

The first group, represented by UGC05986, shows only a single maximum in its
DM RC, i.e. a rising trend in the inner halo followed by a clear turning
point, as can be seen by the data points in the left plot of Fig. 8. In the
same plot for UGC05986 we show the best-fits of the competing DM models as
thick curves, i.e. the fermionic (blue), DC14 (yellow) and NFW (red).
Regarding the fermionic model, this kind of RC is better fitted by the
solutions with a significant escape of particles ($W_{p}\lesssim 10$), as can
be explicitly seen through the $\chi^{2}$ valleys in top panels of Fig. 9.
However, due to the lack of information in the inner halo structures there is
some uncertainty in the strength of particle escape. The uncertainty is
physically better reflected in the core mass $M_{c}$ which covers about two
orders of magnitude (see middle panel of first row in Fig. 9).

Figure 9: $\chi^{2}$ profiles of the fermionic DM model for three benchmark
galaxies: UGC05986 (top panels), DDO161 (middle panels) and NGC6015 (bottom
panels).

This result goes totally in line with an analogous phenomenological analysis
(Argüelles et al., 2019), developed for typical dwarf, spiral and elliptical
galaxies within the RAR model. According to that analysis (done for
$mc^{2}\approx$50\text{\,}\mathrm{keV}$$), the maximal core mass of larger
galaxies is limited by the critical configuration where the quantum core
becomes unstable and collapses to a BH of mass
$M_{c}^{cr}\approx$2\text{\times}{10}^{8}\text{\,}\mathrm{M_{\odot}}$$.

Among the cases, which are disfavored, are the ones with very large total DM
masses $M_{s}$ corresponding to isothermal-like halos and implying negligible
escape of particles ($W_{p}\gtrsim 10$). These solutions provide a minimal
core mass $M_{c}$ with a huge uncertainty in the total mass.

The second group, represented by DDO161, shows a rising part in the RC towards
a maximum without a clear turning point compared to the first group, see
central plots in Fig. 8. Fitting those galaxies for different $W_{p}$ values
does not favor solutions with or without escaping particles effects. The
variation in the $\chi^{2}$ value remains rather small, see middle panels of
Fig. 9.

Finally, the third group, represented by NGC6015, shows some oscillations in
the RC, mainly in the outer halo, see right plots in Fig. 8. There are various
and speculative reasons for the oscillation, e.g., ongoing merging process,
deviation from equilibrium, etc. In any case, those galaxies are clearly
better fitted by extended isothermal-like halos ($W_{p}\gtrsim 10$) — although
being far from good — see bottom panels of Fig. 9. Such solutions provide a
wide halo maximum followed by a flat RC. In contrast, the non-isothermal
solutions with a cutoff provide only a narrow maximum in the halo, followed by
a Keplerian decreasing tail.

It is worth to recall that different DM models such as the fermionic model,
NFW, DC14 and others are not appropriate to fit the oscillations,
characterized through multiple maxima in the RC. All solutions with a wide
halo are suitable to fit the oscillations well on average, although the best-
fits remain rather poor, leaving almost no insight into the physical
properties of DM on halo scales for those galaxies.

##  6 Summary and conclusion

For the case of disk galaxies, as provided by the SPARC data-set (see section
2.1), we have studied the galactic rotation curves and different galaxy
scaling relations — such as the Radial Acceleration Relation, MDAR and DM
surface density relation (SDR) — from an alternative perspective in which the
halos are formed through a MEP.

Within this paradigm we considered the DM halo as a self-gravitating system of
neutral fermions at finite temperature while the baryonic mass components were
provided from the SPARC data-set.

For comparison, we have taken into account empirical DM fitting models
motivated from DM-only simulations like NFW (within CDM) and Burkert (within
WDM); the DC14 (or generalized NFW) model which contains different physics
such as the influence of baryonic feedback in the morphology of CDM halos; and
the Einasto model as recently studied in Bozek et al. (2019) accounting for
baryonic effects (through hydrodynamical zoom-in simulations) in either
cosmology. Finally, their best-fits to the acceleration relations and SPARC
RCs were compared with the fermionic model (see sections 4.1 and 4.2
respectively).

For all competing DM models, we fitted the DM contribution to the RC as
inferred from the given total rotation curve and the baryonic component (see
section 2.2).

An alternative fitting approach is minimizing the least square errors of the
total rotation curve (e.g. $V_{i}=V_{i,\rm tot}$) where the predictions are a
composition of a theoretical DM halo model and the baryonic component inferred
from observation. On theoretical ground, such an approach is not identical to
the approach explained in section 2 since the propagation of uncertainty
produces a somewhat different weighting. However, we have repeated the same
analysis on both approaches and obtained consistent results despite few minor
numerical variations, and without changing any qualitative conclusion obtained
in this work.

The main results of this work can be summarized as follows, according to three
different issues.

###  6.1 Acceleration relations

The Radial Acceleration Relation as well as MDAR analyzed here are based on an
averaging of many spiral galaxies and hold for different Hubble types. Our
analysis shows that all competing DM models are able to reproduce those
relations, although without a clear favourite because all are similarly good
(see Fig. 2). For instance, for all DM models, we obtain nearly identical
values for $\mathfrak{a}_{0}$ as required in Eq. 20. This result is in line
with a recent analysis done in Kaplinghat et al. (2020).

###  6.2 Individual rotation curve fittings

A deeper understanding of the Radial Acceleration Relation and MDAR is backed
by a goodness of model analysis for $120$ filtered and individual galaxies of
the SPARC data set covering different Hubble types. The DM contribution to the
RCs reflects some diversity in galaxies which, in general, are better fitted
by cored DM halo models instead of cuspy (e.g. NFW, see Fig. 4). This is not
only in agreement with a similar analysis done by Li et al. (2020), but
totally in line with the results of Kaplinghat et al. (2020). That is, it has
been shown here that the DM halo models suitable to explain the acceleration
relations do not necessary explain well the SPARC RCs.

Comparing the fitting goodness of the superior DC14 with the inferior (and
statistically disfavoured) NFW implies that baryonic feedback mechanism is
important in galaxy formation. On the other hand, we found that the fermionic
model, compared to DC14 or Einasto, is equally good in fitting galaxies which
require a significant escape of particles (i.e. $W_{p}<10$, see section 4.2
and Fig. 5). Those galaxies are characterized by a flat inner halo, justified
by very different physical principles in the most favored DM profiles: DC14
and Einasto models rely on complex baryonic feedback processes, while the MEP
scenario involves a quasi-thermodynamic equilibrium of the DM particles. This
may imply that for those galaxies baryonic feedback is less relevant, thus
hinting on the importance of a quasi-thermodynamic equilibrium that may be
reached in those DM halos.

###  6.3 Fermionic halos

We found that for SPARC galaxies the constancy of SDR, originally based on the
Burkert model (Donato et al., 2009), is also achievable within the fermionic
DM model. From the Carnegie-Irvine Galaxy Survey (Ho et al., 2011) we have
further extracted the absolute magnitude for nine overlapping SPARC galaxies.
The surface density predictions of that sub-sample are well in agreement with
observations.

Additionally we demonstrated that particular solutions from the rich
morphology of fermionic DM on halo scales can be associated with different
empirical DM models, depending on the strength of particle evaporation
(described by $W_{p}$). Fermionic DM halos are polytropic-like (with $n=5/2$)
for strong evaporation ($W_{p}\ll 1)$, transition into profiles similar to
Burkert for values of $W_{p}$ close to the point of marginal thermodynamical
stability (given at $W_{p}^{(c)}=7.45$) and finally develop isothermal tails
for negligible evaporation ($W_{p}\gg 10$).

Interestingly, the mean $\kappa$ value obtained here for the Einasto model is
roughly $0.4$, implying a pronounced inner-halo density drop (i.e leading to
an almost flat inner-slope, see Fig. 1), such that the halos look more like
the fermionic profiles.

Indeed, following the work done in Argüelles et al. (2021) for such fermionic
profiles, it can be found that typical galaxies belonging to this sub-group
(with $W_{p}\ll 1$) are thermodynamically and dynamically stable, with an
outer-halo morphology of polytropic nature (see Fig. 7). This may be
evidencing a fundamental and deep link between thermodynamics of self-
gravitating fermions and galaxy formation and morphology.

This work was founded by the Consejo Nacional de Investigaciones
Cient\́mathfrak{i}ficas y Técnicas (CONICET), grant number 11220200102876CO.
We thank B. Famey for useful discussions. We thank the anonymous referee who
helped to improve the presentation of this work.

##  Appendix A Parameter correlations of fermionic halos

We analyzed different pairs of structural galaxy parameters, obtained from
core-halo best-fit solutions of the fermionic model for a particle mass of
$50\text{\,}\mathrm{keV}$.

Of interest here are values at the halo such as the halo radius $r_{h}$ and
halo mass $M_{h}=M(r_{h})$. The halo radius $r_{h}$ is defined at the second
maximum in the rotation curve $v(r)$.

For mass and radius we obtain values located mainly in the intervals
$r_{h}\in[10^{3},10^{5}]$\mathrm{pc}$$ and
$M_{h}\in[10^{8},10^{12}]$\mathrm{M_{\odot}}$$. As shown in Fig. 10 the halo
radius and mass follow a clear relation described by

$\ln\left[\vphantom{}\frac{M_{h}}{M_{\odot}}\right]\approx
2\ln\left[\vphantom{}\frac{r_{h}}{\mathrm{pc}}\right]+5.8\pm 1.7$ (A1)

Interestingly, this relation is independent of particle escape, i.e it holds
for solutions with an isothermal-like halo developing a flat tail as well as
for polytropic non-isothermal halos implying a large escape of particles.

Figure 10: Halo parameter correlation for the best-fits of the fermionic
model. The dashed lines indicate the $95\text{\,}\mathrm{\char 37\relax}$ CI.

##  Appendix B Parameter distributions

Figure 11: Distribution of the best-fit parameters for the fermionic model
with a core-halo (top) and fully diluted (bottom). The blue bars represent the
sub-sample of non-isothermal solutions (44 galaxies) with $W_{p}\lesssim 10$
(top) and $W_{0}\lesssim 10$ (bottom), respectively, while the gray bars
include the full sample (120 galaxies). Note that the values at plateau of the
core-halo solutions (i.e. $\beta_{p}$, $\theta_{p}$, $W_{p}$) can be
identified with the values at the center of the halo-only solutions (i.e.
$\beta_{0}$, $\theta_{0}$, $W_{0}$). In the case of halo-only solutions
(bottom) the dashed line represents the stability change where fermionic DM
halos become unstable for $W_{p}\equiv W_{0}\gtrsim 7.45$. (Chavanis et al.,
2015b)

Figure 12: Distribution of the best-fit parameters for the DC14 model. Above
the histogram a boxplot is shown with a median value at $X\approx-2.4$. The
mean value is very close to the median value. The corresponding values for
$\alpha$, $\beta$ and $\gamma$ are calculated from the best-fitted $X$ value
following Eqs. 11, 12 and 13. Figure 13: Distribution of the best-fit
parameter for the Einasto model. A boxplot (above the histogram) is shown with
a median value of $\kappa\approx 0.42$ and mean value $\kappa\approx 0.46$.

We analysed the distribution of configuration parameters which affect the
shape of the DM RC for each DM model. On this basis, we are interested only in
fermionic, DC14 and Einasto models, which are the only ones where the halo
slopes (usually characterized by $\gamma$) depend on one configuration
parameter.

Starting with the core-halo solutions, see upper plots in Fig. 11, we find
that the majority of central temperature values falls in the range
$\beta_{0}\in[${10}^{-8}$,${10}^{-6}$]$ corresponding to solutions in which
the central DM cores are far from reaching the core-collapse towards a SMBH
(Argüelles et al., 2019, 2021). The distribution of the central degeneracy
values looks like a Gaussian with the mean value at about 36, far in the
degenerate regime $\theta_{0}\gtrsim 10$. The majority is in the range
$\theta_{0}\in[30,40]$. A similar distribution pattern is given for $W_{0}$.
However, for core-halo solutions (i.e. $\theta_{0}\gtrsim 10$) it is better to
look at the plateau cutoff $W_{p}$ which can be identified with $W_{0}$ of a
corresponding halo-only solution.

The plateau of a core-halo solution acts as a connection between the halo and
the embedded core. Therefore, every core-halo solution of the fermionic model
has a corresponding halo-only solution describing the diluted halo without the
embedded, degenerate core. The corresponding values for such halo-only
solutions are given at the plateau which resemble the inner halo, see lower
plots in Fig. 11. The plateau temperature shows a very similar distribution
due to the tiny temperature changes when in the low temperature regime
($\beta_{0}\ll${10}^{-4}$$) where pressure effects are negligible leading to
$\beta_{p}\approx\beta_{0}$. The plateau degeneracy distribution looks also
similar to the central degeneracy but being mirrored and shifted to the
negative (diluted) regime. We find the relation
$\theta_{p}\approx-0.7\theta_{0}-1.2$.

Of great interest is the plateau cutoff $W_{p}$ which is a proxy for the
central cutoff $W_{0}$ parameter and provides better insights about the halo
shape. The plateau cutoff describes the particle escape intensity on halo
scales and is defined as $W_{p}=W(r_{p})$ with the plateau radius $r_{p}$
located at the first minimum in the DM RC. The lower $W_{p}$ the more
truncated is the halo due to evaporation. See also section 5.2 for a
discussion about the halo morphology and comparison with other DM models.

Within the fermionic model we identify two groups in the $W_{p}$ distribution
of core-halo solutions and $W_{0}$ distribution of halo-only solution,
respectively: (1) an isothermal (non-truncated) group with 76 galaxies and (2)
a non-isothermal (truncated) group with 44 galaxies. Both groups are divided
at about $W_{p}\approx 10$ (core-halo) and $W_{0}\approx 10$ (halo-only),
respectively.

For the first group the outer halo seems to be isothermal (i.e. characterized
by a flat tail and $\rho(r)\propto r^{-2}$ for large enough $r$). However, for
those galaxies the exact value of $W_{0}$ cannot be determined due to
insufficient and/or too limited information in the outer halo data, and thus
they are not shown in Fig. 11. In contrast, the second group seems to have a
broader data coverage and/or admitting for a cleaner description of the outer
RC, allowing for a better constrain of the outer halo of the fermionic DM
profiles (i.e. with finite $W_{p}\lesssim 10$ and consequent non-isothermal
halo tails). See Fig. 7 for a comparison between typical solution of both
groups.

For halo-only solutions (i.e. $\theta_{0}\lesssim-5)$ there is an interesting
accumulation of galaxies around the particular value of $W_{0}\approx 7.45$,
see dashed line in Fig. 11. This particular value reflects the point of
stability change where the diluted solutions become unstable, i.e. for
$W_{0}\gtrsim 7.45$ (Chavanis et al., 2015b), and may indicate physical
insights into galaxy formation. Interestingly, for this value of $W_{0}$ the
density profile of the classical King model resembles the Burkert profile (see
section 5.2).

In the case of the DC14 model the majority of best-fitted $X$ values is in the
range $[-5,0]$ with a peak at the median of $X\approx-2.4$ as indicated by the
boxplot in Fig. 12. The crosses may be outliers. When we calculate the
parameters $\alpha$, $\beta$ and $\gamma$ with Eqs. 11, 12 and 13, there seems
to be a grouping in the distribution of $\alpha$ with a separation at
$\alpha=0$, see Fig. 12. However, Eq. 8 is not defined for $\alpha=0$.
Phenomenologically, $\alpha$ describes the transition from the inner to the
outer halo. The larger $\left|\alpha\right|$ the more extended is the
transition, characterized by a long wide maximum. $\beta$ describes the slope
in the outer halo while $\gamma$ describes the slope in the inner halo. For
$\gamma=0$ the DM profiles become cored.

The Einasto model has only a single configuration parameter $\kappa$
describing the transition from the cored inner halo to the outer halo. The
larger $\kappa$ the less extended is the RC maximum. As shown in Fig. 13 the
majority has $\kappa<1$ that corresponds to a rather extended RC maximum. This
distribution represents well the majority of SPARC galaxies showing an
extended outer halo trend without a clear maximum in the outer RC.

## References

  * Alberti & Chavanis (2020) Alberti, G., & Chavanis, P.-H. 2020, European Physical Journal B, 93, 208, doi: 10.1140/epjb/e2020-100557-6
  * Argüelles et al. (2022a) Argüelles, C. R., Becerra-Vergara, E. A., Krut, A., et al. 2022a, International Journal of Modern Physics D, 31, 2230002, doi: 10.1142/S0218271822300026
  * Argüelles et al. (2021) Argüelles, C. R., D\́mathfrak{i}az, M. I., Krut, A., & Yunis, R. 2021, MNRAS, 502, 4227, doi: 10.1093/mnras/staa3986
  * Argüelles et al. (2018) Argüelles, C. R., Krut, A., Rueda, J. A., & Ruffini, R. 2018, Physics of the Dark Universe, 21, 82, doi: 10.1016/j.dark.2018.07.002
  * Argüelles et al. (2019) —. 2019, Physics of the Dark Universe, 24, 100278, doi: 10.1016/j.dark.2019.100278
  * Argüelles et al. (2022b) Argüelles, C. R., Mestre, M. F., Becerra-Vergara, E. A., et al. 2022b, MNRAS, 511, L35, doi: 10.1093/mnrasl/slab126
  * Becerra-Vergara et al. (2020) Becerra-Vergara, E. A., Argüelles, C. R., Krut, A., Rueda, J. A., & Ruffini, R. 2020, A&A, 641, A34, doi: 10.1051/0004-6361/201935990
  * Becerra-Vergara et al. (2021) —. 2021, MNRAS, 505, L64, doi: 10.1093/mnrasl/slab051
  * Behroozi et al. (2013) Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57, doi: 10.1088/0004-637X/770/1/57
  * Binney & Tremaine (2008) Binney, J., & Tremaine, S. 2008, Galactic Dynamics, 2nd edn. (Princeton University Press), 920
  * Bogdán & Goulding (2015) Bogdán, Á., & Goulding, A. D. 2015, ApJ, 800, 124, doi: 10.1088/0004-637X/800/2/124
  * Bozek et al. (2019) Bozek, B., Fitts, A., Boylan-Kolchin, M., et al. 2019, MNRAS, 483, 4086, doi: 10.1093/mnras/sty3300
  * Bullock & Boylan-Kolchin (2017) Bullock, J. S., & Boylan-Kolchin, M. 2017, ARA&A, 55, 343, doi: 10.1146/annurev-astro-091916-055313
  * Burkert (1995) Burkert, A. 1995, ApJ, 447, L25, doi: 10.1086/309560
  * Cappellari et al. (2011) Cappellari, M., Emsellem, E., Krajnović, D., et al. 2011, MNRAS, 413, 813, doi: 10.1111/j.1365-2966.2010.18174.x
  * Chavanis (1998) Chavanis, P.-H. 1998, MNRAS, 300, 981, doi: 10.1046/j.1365-8711.1998.01867.x
  * Chavanis (2004) —. 2004, Physica A Statistical Mechanics and its Applications, 332, 89, doi: 10.1016/j.physa.2003.09.061
  * Chavanis (2020) —. 2020, European Physical Journal Plus, 135, 290, doi: 10.1140/epjp/s13360-020-00268-0
  * Chavanis (2021) —. 2021, arXiv e-prints, arXiv:2112.13664.  https://arxiv.org/abs/2112.13664
  * Chavanis et al. (2015a) Chavanis, P.-H., Lemou, M., & Méhats, F. 2015a, Phys. Rev. D, 92, 123527, doi: 10.1103/PhysRevD.92.123527
  * Chavanis et al. (2015b) —. 2015b, Phys. Rev. D, 91, 063531, doi: 10.1103/PhysRevD.91.063531
  * Chavanis et al. (1996) Chavanis, P. H., Sommeria, J., & Robert, R. 1996, ApJ, 471, 385, doi: 10.1086/177977
  * de Blok et al. (2008) de Blok, W. J. G., Walter, F., Brinks, E., et al. 2008, AJ, 136, 2648, doi: 10.1088/0004-6256/136/6/2648
  * Di Cintio et al. (2014) Di Cintio, A., Brook, C. B., Dutton, A. A., et al. 2014, MNRAS, 441, 2986, doi: 10.1093/mnras/stu729
  * Di Cintio & Lelli (2016) Di Cintio, A., & Lelli, F. 2016, MNRAS, 456, L127, doi: 10.1093/mnrasl/slv185
  * Di Paolo et al. (2019) Di Paolo, C., Salucci, P., & Fontaine, J. P. 2019, ApJ, 873, 106, doi: 10.3847/1538-4357/aaffd6
  * Domcke & Urbano (2015) Domcke, V., & Urbano, A. 2015, J. Cosmology Astropart. Phys, 2015, 002, doi: 10.1088/1475-7516/2015/01/002
  * Donato et al. (2009) Donato, F., Gentile, G., Salucci, P., et al. 2009, MNRAS, 397, 1169, doi: 10.1111/j.1365-2966.2009.15004.x
  * Dutton & Macciò (2014) Dutton, A. A., & Macciò, A. V. 2014, MNRAS, 441, 3359, doi: 10.1093/mnras/stu742
  * Dutton et al. (2019) Dutton, A. A., Macciò, A. V., Obreja, A., & Buck, T. 2019, MNRAS, 485, 1886, doi: 10.1093/mnras/stz531
  * Einasto & Haud (1989) Einasto, J., & Haud, U. 1989, A&A, 223, 89
  * Ferrarese (2002) Ferrarese, L. 2002, ApJ, 578, 90, doi: 10.1086/342308
  * Ferrarese & Merritt (2000) Ferrarese, L., & Merritt, D. 2000, ApJ, 539, L9, doi: 10.1086/312838
  * Fitts et al. (2017) Fitts, A., Boylan-Kolchin, M., Elbert, O. D., et al. 2017, MNRAS, 471, 3547, doi: 10.1093/mnras/stx1757
  * Governato et al. (2012) Governato, F., Zolotov, A., Pontzen, A., et al. 2012, MNRAS, 422, 1231, doi: 10.1111/j.1365-2966.2012.20696.x
  * Hernquist (1990) Hernquist, L. 1990, ApJ, 356, 359, doi: 10.1086/168845
  * Ho et al. (2011) Ho, L. C., Li, Z.-Y., Barth, A. J., Seigar, M. S., & Peng, C. Y. 2011, ApJS, 197, 21, doi: 10.1088/0067-0049/197/2/21
  * Kadomtsev & Pogutse (1970) Kadomtsev, B. B., & Pogutse, O. P. 1970, Phys. Rev. Lett., 25, 1155, doi: 10.1103/PhysRevLett.25.1155
  * Kaplinghat et al. (2020) Kaplinghat, M., Ren, T., & Yu, H.-B. 2020, J. Cosmology Astropart. Phys, 2020, 027, doi: 10.1088/1475-7516/2020/06/027
  * Katz et al. (2017) Katz, H., Lelli, F., McGaugh, S. S., et al. 2017, MNRAS, 466, 1648, doi: 10.1093/mnras/stw3101
  * Keller & Wadsley (2017) Keller, B. W., & Wadsley, J. W. 2017, ApJ, 835, L17, doi: 10.3847/2041-8213/835/1/L17
  * King (1962) King, I. 1962, AJ, 67, 471, doi: 10.1086/108756
  * King (1966) King, I. R. 1966, AJ, 71, 64, doi: 10.1086/109857
  * Kormendy & Bender (2011) Kormendy, J., & Bender, R. 2011, Nature, 469, 377, doi: 10.1038/nature09695
  * Kroupa (2015) Kroupa, P. 2015, Canadian Journal of Physics, 93, 169, doi: 10.1139/cjp-2014-0179
  * Kull et al. (1996) Kull, A., Treumann, R. A., & Boehringer, H. 1996, ApJ, 466, L1, doi: 10.1086/310159
  * Lelli et al. (2016) Lelli, F., McGaugh, S. S., & Schombert, J. M. 2016, AJ, 152, 157, doi: 10.3847/0004-6256/152/6/157
  * Lelli et al. (2017) Lelli, F., McGaugh, S. S., Schombert, J. M., & Pawlowski, M. S. 2017, ApJ, 836, 152, doi: 10.3847/1538-4357/836/2/152
  * Li et al. (2018) Li, P., Lelli, F., McGaugh, S., & Schombert, J. 2018, A&A, 615, A3, doi: 10.1051/0004-6361/201732547
  * Li et al. (2020) —. 2020, ApJS, 247, 31, doi: 10.3847/1538-4365/ab700e
  * Ludlow et al. (2017) Ludlow, A. D., Ben\́mathfrak{i}tez-Llambay, A., Schaller, M., et al. 2017, Phys. Rev. Lett., 118, 161103, doi: 10.1103/PhysRevLett.118.161103
  * Lynden-Bell (1967) Lynden-Bell, D. 1967, MNRAS, 136, 101, doi: 10.1093/mnras/136.1.101
  * Marra et al. (2020) Marra, V., Rodrigues, D. C., & de Almeida, Á. O. F. 2020, MNRAS, 494, 2875, doi: 10.1093/mnras/staa890
  * McGaugh (2014) McGaugh, S. 2014, Galaxies, 2, 601, doi: 10.3390/galaxies2040601
  * McGaugh (2004) McGaugh, S. S. 2004, ApJ, 609, 652, doi: 10.1086/421338
  * McGaugh et al. (2016) McGaugh, S. S., Lelli, F., & Schombert, J. M. 2016, Phys. Rev. Lett., 117, 201101, doi: 10.1103/PhysRevLett.117.201101
  * McGaugh et al. (2000) McGaugh, S. S., Schombert, J. M., Bothun, G. D., & de Blok, W. J. G. 2000, ApJ, 533, L99, doi: 10.1086/312628
  * Merritt et al. (2006) Merritt, D., Graham, A. W., Moore, B., Diemand, J., & Terzić, B. 2006, AJ, 132, 2685, doi: 10.1086/508988
  * Milgrom (2016) Milgrom, M. 2016, arXiv e-prints, arXiv:1609.06642.  https://arxiv.org/abs/1609.06642
  * Navarro et al. (2017) Navarro, J. F., Ben\́mathfrak{i}tez-Llambay, A., Fattahi, A., et al. 2017, MNRAS, 471, 1841, doi: 10.1093/mnras/stx1705
  * Navarro et al. (1996) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563, doi: 10.1086/177173
  * Navarro et al. (1997) —. 1997, ApJ, 490, 493, doi: 10.1086/304888
  * Read et al. (2016) Read, J. I., Agertz, O., & Collins, M. L. M. 2016, MNRAS, 459, 2573, doi: 10.1093/mnras/stw713
  * Retana-Montenegro et al. (2012) Retana-Montenegro, E., van Hese, E., Gentile, G., Baes, M., & Frutos-Alfaro, F. 2012, A&A, 540, A70, doi: 10.1051/0004-6361/201118543
  * Ruffini et al. (2015) Ruffini, R., Argüelles, C. R., & Rueda, J. A. 2015, MNRAS, 451, 622, doi: 10.1093/mnras/stv1016
  * Ruffini & Bonazzola (1969) Ruffini, R., & Bonazzola, S. 1969, Physical Review, 187, 1767, doi: 10.1103/PhysRev.187.1767
  * Ruffini & Stella (1983) Ruffini, R., & Stella, L. 1983, A&A, 119, 35
  * Sales et al. (2017) Sales, L. V., Navarro, J. F., Oman, K., et al. 2017, MNRAS, 464, 2419, doi: 10.1093/mnras/stw2461
  * Salucci (2016) Salucci, P. 2016, arXiv e-prints, arXiv:1612.08857.  https://arxiv.org/abs/1612.08857
  * Salucci (2018) —. 2018, Foundations of Physics, 48, 1517, doi: 10.1007/s10701-018-0209-5
  * Salucci & Burkert (2000) Salucci, P., & Burkert, A. 2000, ApJ, 537, L9, doi: 10.1086/312747
  * Sánchez Almeida & Trujillo (2021) Sánchez Almeida, J., & Trujillo, I. 2021, MNRAS, 504, 2832, doi: 10.1093/mnras/stab1103
  * Severne & Luwel (1980) Severne, G., & Luwel, M. 1980, Ap&SS, 72, 293, doi: 10.1007/BF00639139