Galaxy Rotation Curves and Universal Scaling Relations: Comparison between Phenomenological and Fermionic Dark Matter Profiles

The SPARC data set includes 3.6 μm near-infrared and 21 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_tot and the baryonic contribution V_bar, composed of a bulge (V_b), disk (V_d), and gas component (V_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

$$\begin{eqnarray}&&{V}_{\mathrm{bar}}^{2}={{\rm{\Upsilon }}}_{{\rm{b}}}^{}{V}_{{\rm{b}}}^{2}+{{\rm{\Upsilon }}}_{{\rm{d}}}^{}{V}_{{\rm{d}}}^{2}+{V}_{{\rm{g}}}^{2}.\end{eqnarray} \tag{ 1 }$$

For convenience, the given velocities V_b and V_d in the SPARC data are normalized for a mass-to-light ratio of 1 M_⊙/L_⊙.

In this work, we follow the same data selection criteria as done in McGaugh et al. (2016). We choose averaged mass-to-light ratios ϒ_b = 0.7 for all bulges, and ϒ_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°. The latter is to minimize the $\sin (i)$ corrections to the observed velocities. For all measurements, we require a minimum precision of 10% in velocity.

Additionally we reject all points where the baryonic velocity is greater than 95% 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. The 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 RCs.

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

$$\begin{eqnarray}&&{V}_{\mathrm{DM}}^{2}={V}_{\mathrm{tot}}^{2}-{V}_{\mathrm{bar}}^{2},\end{eqnarray} \tag{ 2 }$$

where the uncertainty in V_DM is calculated from the uncertainties in the other velocity components within the linear error propagation theory. Note that only the uncertainties in the total RCs (ΔV_tot) are provided within the SPARC data set. The uncertainty is therefore given by

$$\begin{eqnarray}&&{\rm{\Delta }}{V}_{\mathrm{DM}}=\left|\displaystyle \frac{\partial {V}_{\mathrm{DM}}}{\partial {V}_{\mathrm{tot}}}\right|{\rm{\Delta }}{V}_{\mathrm{tot}}.\end{eqnarray} \tag{ 3 }$$

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

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

$$\begin{eqnarray}&&{\chi }^{2}({\boldsymbol{p}})=\displaystyle \sum _{i=1}^{N}{\left[\displaystyle \frac{{V}_{i}-v({r}_{i},{\boldsymbol{p}})}{{\rm{\Delta }}{V}_{i}}\right]}^{2}.\end{eqnarray} \tag{ 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 , p ) accounts for the circular velocity at r_i for each model parameter vector p (described below), and ΔV_i is the uncertainty in V_i as given in Equation (3).

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 noninteracting, apart from gravitational forces, the relaxation time toward statistical equilibrium due to gravitational encounters scales as $(N/\mathrm{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² ∼ 50 keV, a DM halo of mass M ∼ 10¹¹ M_⊙ and radius R ∼ 30 kpc contain N ∼ 10⁷² fermions for a dynamical time ${t}_{D}\sim 1/\sqrt{{R}^{3}/{GM}}\sim 100\,\mathrm{Myr}$.

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), Chavanis et al. (2015a).

Assuming ergodicity (efficient mixing), Lynden-Bell (1967) used an 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}({\boldsymbol{r}},{\boldsymbol{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 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 2022 for a review).

One approach is to use a maximum entropy production principle 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

$$\begin{eqnarray}&&\displaystyle \frac{\partial \bar{f}}{\partial t}+{\boldsymbol{v}}\cdot \displaystyle \frac{\partial \bar{f}}{\partial {\boldsymbol{r}}}-{\rm{\nabla }}{\rm{\Phi }}\cdot \displaystyle \frac{\partial \bar{f}}{\partial {\boldsymbol{v}}}=\displaystyle \frac{\partial {\boldsymbol{J}}}{\partial {\boldsymbol{v}}},\end{eqnarray} \tag{ 5 }$$

where ${\boldsymbol{J}}=D[\partial \bar{f}/\partial {\boldsymbol{v}}+({{mc}}^{2}/{kT})\bar{f}(1-f/{\eta }_{0}){\boldsymbol{v}}]$ is a diffusion current pushing the system toward statistical equilibrium, D is the diffusion coefficient, T ≡ 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

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial \bar{f}}{\partial t}+{\boldsymbol{v}}\cdot \displaystyle \frac{\partial \bar{f}}{\partial {\boldsymbol{r}}}-{\rm{\nabla }}{\rm{\Phi }}\cdot \displaystyle \frac{\partial \bar{f}}{\partial {\boldsymbol{v}}}=\displaystyle \frac{8\pi {G}^{2}{m}^{8}{\epsilon }_{r}^{3}{\epsilon }_{v}^{3}\mathrm{ln}{\rm{\Lambda }}}{{h}^{6}}\displaystyle \frac{\partial }{\partial {v}_{i}}\\ \,\times \displaystyle \int d{\boldsymbol{v}}^{\prime} \displaystyle \frac{{u}^{2}{\delta }_{{ij}}-{u}_{i}{u}_{j}}{{u}^{3}}\left\{\bar{f}^{\prime} \left(1-\bar{f}^{\prime} \right)\displaystyle \frac{\partial \bar{f}}{\partial {v}_{j}}-\bar{f}\left(1-\bar{f}\right)\displaystyle \frac{\partial \bar{f}^{\prime} }{\partial v{{\prime} }_{j}}\right\},\end{array}\end{eqnarray} \tag{ 6 }$$

where $\bar{f}^{\prime} \equiv \bar{f}({\boldsymbol{r}},{\boldsymbol{v}}^{\prime} ,t)$, $\mathrm{ln}{\rm{\Lambda }}=\mathrm{ln}(R/{\epsilon }_{r})$ is the Coulomb logarithm, R is the typical size of the system, ${\boldsymbol{u}}={\boldsymbol{v}}^{\prime} -{\boldsymbol{v}}$ is the relative velocity between the "macro-particles" of mass ${m}_{\mathrm{eff}}\sim 2{m}^{4}{\epsilon }_{r}^{3}{\epsilon }_{v}^{3}/{h}^{3}\gg m$ (see also below), and _r , _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 Equation (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}^{\mathrm{non}-\mathrm{coll}.}\,\sim (M/{m}_{\mathrm{eff}}){t}_{D}$, which is much shorter than the collisional relaxation time ${t}_{R}^{\mathrm{coll}.}\sim (M/m){t}_{D}$ because m_eff ≫ m (see formula in the above paragraph).

Indeed, the relaxation of the coarse-grained distribution function $\bar{f}({\boldsymbol{r}},{\boldsymbol{v}},t)$ toward the Lynden-Bell distribution (of Fermi–Dirac type; see Equation (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_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 Equation (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 Equation (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)

$$\begin{eqnarray}&&\bar{f}(r,\epsilon \leqslant {\epsilon }_{c})=\displaystyle \frac{1-{e}^{[\epsilon -{\epsilon }_{c}(r)]/{kT}(r)}}{{e}^{[\epsilon -\mu (r)]/{kT}(r)}+1},\qquad \bar{f}(r,\epsilon \gt {\epsilon }_{c})=0,\end{eqnarray} \tag{ 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, μ(r) is the chemical potential (with the particle rest-energy subtracted off), _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, β(r) = kT(r)/(mc²), θ(r) = μ(r)/[kT(r)], and W(r) = _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, 2022b; Becerra-Vergara et al. 2020, 2021).

The full family of density ρ(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 ≤ _c (r)) as detailed in Argüelles et al. (2018). This leads to a four-parametric fermionic equation of state depending on (β₀, θ₀, W₀, m) according to the parameters in Equation (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 }=\mathrm{diag}({e}^{2\nu (r)},-{e}^{2\lambda (r)},-{r}^{2},-{r}^{2}{\sin }^{2}\vartheta )$ with ν(r), λ(r) being the temporal and spatial metric functions, and ϑ is the azimutal angle. Such configuration leads to hydrostatic equilibrium equations of self-gravitating fermions. The local T(r), μ(r), and W(r) fulfill the Tolman, Klein, and particle's energy conservation relations, respectively (see Argüelles et al. 2018 for details). Further, ν₀ is here constrained by the Schwarzschild condition g₀₀ g₁₁ = −1 at the surface where the halo pressure (and density) falls to zero.

An illustration of a core-halo (θ₀ > 10) and a corresponding halo-only (θ₀ ≪ −1) solution are shown in Figure 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).

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 RC.

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 RC, the halo mass M_h = M(r_h ) with r_h being the halo radius defined at the second maximum in the RC, 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 behavior: while the halo region is King-like (i.e., from polytropic-like with W_p ≪ 1 to power law-like with W_p ≳ 10; see Section 5.2), the inner region can either develop a dense core at the center of such a halo (i.e., for large central degeneracy θ₀ > 10), or not (i.e., θ₀ ≪ −1 in the dilute regime). Both kinds 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 ∼5 × 10¹⁰ M_⊙.

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

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 (α, β, γ)-model (Hernquist 1990)

$$\begin{eqnarray}&&\displaystyle \frac{\rho (r)}{{\rho }_{{\rm{N}}}}={\left[\displaystyle \frac{r}{{R}_{{\rm{N}}}}\right]}^{-\gamma }{\left[1+{\left[\displaystyle \frac{r}{{R}_{{\rm{N}}}}\right]}^{\alpha }\right]}^{-\displaystyle \frac{\beta -\gamma }{\alpha }},\end{eqnarray} \tag{ 8 }$$

where α describes the transition, β is the outer slope, and γ is the inner-slope. Following Newtonian dynamics—that is fully sufficient on halo scales—then the velocity is given by

$$\begin{eqnarray}&&\displaystyle \frac{{v}^{2}(r)}{{\sigma }_{{\rm{N}}}^{2}}=\displaystyle \frac{{R}_{{\rm{N}}}}{r}\displaystyle \frac{M(r)}{{M}_{{\rm{N}}}},\end{eqnarray} \tag{ 9 }$$

and the enclosed mass is given by

$$\begin{eqnarray}&&\displaystyle \frac{M(r)}{{M}_{{\rm{N}}}}=\underset{0}{\overset{r}{\int }}{\left[\displaystyle \frac{r}{{R}_{{\rm{N}}}}\right]}^{2}\displaystyle \frac{\rho (r)}{{\rho }_{{\rm{N}}}}\displaystyle \frac{{dr}}{{R}_{{\rm{N}}}}.\end{eqnarray} \tag{ 10 }$$

For the following DM halo models, we will use R_N, ρ_N, ${\sigma }_{{\rm{N}}}^{2}={{GM}}_{{\rm{N}}}/{R}_{{\rm{N}}}$, and ${M}_{{\rm{N}}}=4\pi {\rho }_{{\rm{N}}}{R}_{{\rm{N}}}^{3}$ as scaling factors for length, density, velocity, and mass, respectively.

In Figure 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 (α, β, γ)-model, Cintio et al. (2014) modeled CDM halos including baryonic feedback mechanisms in galaxy formation. They found that the three parameters (α, β, γ) are related through

$$\begin{eqnarray}&&\alpha =2.94-{\mathrm{log}}_{10}\left[{({10}^{X+2.33})}^{-1.08}+{({10}^{X+2.33})}^{2.29}\right],\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&\beta =4.23+1.34X+0.26{X}^{2},\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&\gamma =-0.06+{\mathrm{log}}_{10}\left[{({10}^{X+2.56})}^{-0.68}+{10}^{X+2.56}\right],\end{eqnarray} \tag{ 13 }$$

where $X={\mathrm{log}}_{10}({M}_{* }/{M}_{\mathrm{halo}})$ describes the stellar-to-DM ratio. For the circular velocity (Equation (9)), the enclosed mass (Equation (10)) is given by a hypergeometric function

$$\begin{eqnarray}&&\displaystyle \frac{M(r)}{{M}_{{\rm{N}}}}=\displaystyle \frac{1}{3-\gamma }{\left[\displaystyle \frac{r}{{R}_{{\rm{N}}}}\right]}^{3-\gamma }{}_{2}{F}_{1}({p}_{1},{p}_{2};{q}_{1};-{[r/{R}_{{\rm{N}}}]}^{\alpha }),\end{eqnarray} \tag{ 14 }$$

with p₁ = (3 −γ)/α, p₂ = (β − γ)/α, and q₁ = 1 + (3 − γ)/α. In the following, we refer this model as DC14.

Alternatively, for α = 1, β = 3, and γ = 1, the (α, β, γ)-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:

$$\begin{eqnarray}&&\displaystyle \frac{\rho (r)}{{\rho }_{{\rm{N}}}}={\left[\displaystyle \frac{r}{{R}_{{\rm{N}}}}\right]}^{-1}{\left[1+\displaystyle \frac{r}{{R}_{{\rm{N}}}}\right]}^{-2},\end{eqnarray} \tag{ 15 }$$

with the circular velocity

$$\begin{eqnarray}&&\displaystyle \frac{{v}^{2}(r)}{{\sigma }_{{\rm{N}}}^{2}}=\displaystyle \frac{\mathrm{ln}(1+r/{R}_{{\rm{N}}})}{r/{R}_{{\rm{N}}}}-\displaystyle \frac{1}{1+r/{R}_{{\rm{N}}}}.\end{eqnarray} \tag{ 16 }$$

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

$$\begin{eqnarray}&&\displaystyle \frac{\rho (r)}{{\rho }_{{\rm{N}}}}={\left[1+\displaystyle \frac{r}{{R}_{{\rm{N}}}}\right]}^{-1}{\left[1+{\left(\displaystyle \frac{r}{{R}_{{\rm{N}}}}\right)}^{2}\right]}^{-1}.\end{eqnarray} \tag{ 17 }$$

For the circular velocity (Equation (9)), the enclosed mass (Equation (10)) is given by

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{M(r)}{{M}_{{\rm{N}}}} & = & \displaystyle \frac{1}{4}\mathrm{ln}(1+{[r/{R}_{{\rm{N}}}]}^{2})+\displaystyle \frac{1}{2}\mathrm{ln}(1+r/{R}_{{\rm{N}}})\\ & & -\,\displaystyle \frac{1}{2}\arctan (r/{R}_{{\rm{N}}}).\end{array}\end{eqnarray} \tag{ 18 }$$

With M₀ ≈ 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₀. Further, the density scale ρ_N describes the central density ρ₀, and the length scale R_N can be identified with the Burkert radius r_B fulfilling the condition ρ(r_B) = ρ₀/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:

$$\begin{eqnarray}&&\displaystyle \frac{\rho (r)}{{\rho }_{{\rm{N}}}}={e}^{-{[r/{R}_{{\rm{N}}}]}^{\kappa }}.\end{eqnarray} \tag{ 19 }$$

The exponent κ describes the shape of the density profile. The circular velocity and the enclosed mass are given by Equations (9) and (10). This model develops mass distributions with a finite mass M_tot/M_N = Γ(3/κ)/κ for r → ∞ (see also Retana-Montenegro et al. 2012). The typical κ 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 Sections 3.3, 3.4, and 4.2, and Figure 12.

For the fermionic DM model, we fix the particle mass m and therefore reduce the number of free parameters by one, e.g., p = (β₀, θ₀, W₀). A particle mass of mc² = 50 keV is well motivated by the promising results obtained in Argüelles et al. (2018, 2022b), Becerra-Vergara et al. (2020, 2021), where the fermionic core-halo DM profile was able to explain both the S-stars orbits around SgrA*, and the Milky Way RC. 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 ∼50 and ∼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 that are either in the dilute regime (θ₀ ≪ −1, i.e., are King-like) or have developed a degenerate core (i.e., θ₀ > 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., p = (R_N, ρ_N). The DC14 model with, e.g., p = (X, R_N, ρ_N), and the Einasto model with, e.g., p = (κ, R_N, ρ_N) are described by three free parameters. Compared to NFW and Burkert, both (Einasto and DC14) have an additional parameter that 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 β₀ ∈ [10⁻⁸, 10⁻⁵], θ₀ ∈ [25, 45], and W₀ ∈ [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 β₀ and W₀, but with θ₀ ≡ θ_p ∈ [ −40, −20]. For the other DM models, the initial scaling factors are chosen from R_N ∈ [10¹, 10⁴] pc and ρ_N ∈ [10⁻⁴, 10⁻¹]M_⊙ pc⁻³. According to Katz et al. (2017), the additional parameter of the DC14 model is chosen from X ∈ [ −3.75, −1.3]. For the Einasto model, the exponent is chosen from κ ∈ [0.1, 10]. This relatively large window has been chosen in order (i) to account for the broad diversity of RCs covered by the SPARC data, and (ii) not to be limited by any fixed value (e.g., as typically obtained by CDM-only simulations) because 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).

The DM in galaxies evolves together with baryons, thus is expected in 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 Figure 2 roughly fulfill a_tot < 2a_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 has always been calibrated and applied to cuspy CDM halos. Only recently it was shown that baryonic feedback is DM-model dependent.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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. Because 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, although it is 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 (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.

The RC of each component (e.g., bulge, disk, gas) traces its centripetal acceleration a = v²/r, giving access to independent acceleration measurements. The radial acceleration relation compares the radial acceleration due to the total mass (a_tot) with that due to the baryonic mass (a_bar). It is empirically described by McGaugh et al. (2016):

$$\begin{eqnarray}&&{a}_{\mathrm{tot}}=\displaystyle \frac{{a}_{\mathrm{bar}}}{1-{e}^{-\sqrt{{a}_{\mathrm{bar}}/{{\mathfrak{a}}}_{0}}}},\end{eqnarray} \tag{ 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 Figure 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; Paolo et al. 2019).

For the SPARC galaxies, the different accelerations ${a}_{\mathrm{bar}}={V}_{\mathrm{bar}}^{2}/r$, and ${a}_{\mathrm{tot}}={V}_{\mathrm{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}_{\mathrm{bar}}^{2}/r$ is equally inferred from the SPARC data, while ${a}_{\mathrm{tot}}=[{V}_{\mathrm{bar}}^{2}+{v}_{\mathrm{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 DM component v_DM(r) for each DM model. See Sections 2 and 3.3 for details on 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\times {10}^{-10}\,{\rm{m}}\,{{\rm{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 Figure 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_tot/M_bar with M_bar being the total baryonic mass, and M_tot is the total galaxy mass accounting for baryonic and DM. For $a=\displaystyle \frac{d{\rm{\Phi }}}{{dr}}$ and Φ(r) being the gravitational potential of a spherically symmetry mass distribution, the mass discrepancy can be equivalently written as D = a_tot/a_bar. The results are illustrated in the bottom panels of Figure 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 ΛCDM framework (Di Cintio & Lelli 2016; Salucci 2016, 2018; Katz et al. 2017). For smaller disk and lower sideband galaxies (extending the original SPARC sample), Paolo et al. (2019) found that Equation (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.

Additional 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 RCs allowing for another (related) quantitative comparison of the DM models.

The SPARC galaxies show different characteristics in their RC 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 Figure 7 for three typical examples within the SPARC data set. Some galaxies show just a rising trend implying that the RCs 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 χ² value; see Equation (4). When competing models with different numbers of parameters are compared, it is appropriate to consider the reduced χ² 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 is the number of parameters (of the considered model).

The question now arises of how to compare the competing models for a population of galaxies. In order to find the goodness of a model that is robust against outliers, we ask how many fitted galaxies have a (reduced) χ² lower than a given one. It turns out that the population curve resembling a cumulative distribution function follows nearly a log-normal distribution. We use the mean value, labeled as ${\hat{\chi }}_{r}^{2}$, 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.

The goodness analysis for the entire SPARC sample (120 galaxies) is shown in Figure 3 and can be summarized as follows: the NFW model (p = 2) is statistically disfavored 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 favored with similarly good results; the Burkert model (p = 2) and the fermionic model (p = 3) produce comparable results but are somewhat statistically less favored.

Zoom In Zoom Out Reset image size

We recall that the DM models are suited for DM halo RCs with only one maximum. Therefore, we restrict further the SPARC sample to galaxies (44) showing one clear maximum in their RC (see also Section 5.3 for details about the diversity of SPARC galaxies). This condition is equivalently reflected in the fermionic scenario where data supports a significant escape of DM particles, which is for energy-cutoff values at plateau of W_p < 10 (see Appendix B for further details). For that subsample, the picture changes considerably—see Figure 4— and can be summarized as follows: the NFW model is statistically even more disfavored while all other models (including fermionic DM) become more comparable—with a little tendency for Einasto and against Burkert.

Zoom In Zoom Out Reset image size

We remind that the DC14 model is based on the analysis of hydrodynamically simulated galaxies including complex baryonic feedback processes (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 because, as can be seen from Figure 11 (bottom panel), the bulk of our (α, β, γ) parameters for the SPARC data set lies within the windows (0, 2.6); (2.3, 4); (0, 2) in rough agreement with Cintio et al. (2014). Interestingly, we obtain a mean for the stellar-to-DM mass ratio $X={\mathrm{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 Cintio et al. (2014). Besides, our results are in good agreement with those reported in Behroozi et al. (2013) from abundance matching, because they obtain X values between (−2.7, −1.7) for halo masses in the range between 10¹⁰ M_⊙ and 10¹² M_⊙, with X = −2.4, for M_h = 10¹¹ M_⊙, 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 because 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) hydrodynamical simulations and for ∼10¹⁰ M_⊙ halos, a reduction on inner-halo densities (i.e., at the convergence radius r = 0.26 kpc) of up to 45% in WDM cosmologies with respect to the analogous WDM-only simulations. The reduction 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 κ) as originally shown through Figures 13 and 14 in Dutton & Macciò (2014). Thus, considering that κ ≈ 0.2 corresponds to CDM-only cosmologies (Dutton & Macciò 2014; Fitts et al. 2017), then larger κ 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) because we find a mean value of κ ≈ 0.4 (see Figure 12), thus implying smaller concentration parameters than for CDM profiles. However more work from hydrodynamical 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 RCs (see, e.g., right box in Figure 7). In any case, it does not allow to univocally determine the cutoff parameter (i.e., W_p ) for the fermionic DM model because any sufficiently large W_p would not change the χ² value (see, e.g., bottom left panel of Figure 8). Nevertheless, there are some other individual galaxies where the escaping particles effects are clearly preferred. Interestingly, all of those galaxies are of Magellanic type: NGC 0055 (Sm), UGC 05986 (Sm), UGC 05750 (Sdm), UGC 05005 (Im), F565-V2, (Im), UGC 06399 (Sm), UGC 10310 (Sm), UGC 07559 (Im), UGC 07690 (Im), UGC 05918 (Im), and UGC 05414 (Im).

Moreover, many galaxies, which are poorly fitted by any of the considered models, show short range oscillations in their RCs 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., NGC 2403 (Scd), UGC 02953 (Sab), NGC 6015 (Scd), UGC 09133 (Sab), UGC 06787 (Sab), UGC 11914 (Sab), NGC 1003 (Scd), NGC 0247 (Sd), UGC 08699 (Sab), and UGC 03205 (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 ≲ 10). Whether with weak (W_p ≳ 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., NGC 5585 (Sd), NGC 7793 (Sd), UGC 06614 (Sa), ESO079-G014 (Sbc), F571-8 (Sc), NGC 0891 (Sb), UGC 06614 (Sa), UGC 09037 (Scd), NGC 4217 (Sb), UGC 04278 (Sd).

NFW and Burkert models cannot explain variations of the inner and outer RC because the parameters (β and γ) responsible for such a behavior (see Equation (8)) are fixed. Additionally a transition from the inner to the outer halo is generally characterized by α (or κ). In contrast to NFW and Burkert, the DC14 and Einasto models have a free parameter, which affects the inner–outer RC steepness and the sharpness of the halo transition, simultaneously. Such a flexibility is reflected in generally better χ² values. Nevertheless, the goodness for oscillating RCs remains rather poor.

The constant surface density (Donato et al. 2009)

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{0{\rm{D}}}={\rho }_{0{\rm{D}}}{r}_{0}\approx {140}_{-50}^{+80}{M}_{\odot }\,{\mathrm{pc}}^{-2}\end{eqnarray} \tag{ 21 }$$

is valid for about 14 orders of magnitude in absolute magnitude (M_B ), where ρ₀ is the central DM halo density, and r₀ is the one-halo-scale-length—both of the Burkert model. At r₀ the density falls to one-forth of the central density, i.e., ρ(r₀) = ρ_0D/4.

Note that the center in the Burkert model corresponds to the plateau in the fermionic DM model, i.e., ρ_0D ≈ ρ_p where the plateau density ρ_p is defined at the first minimum in the RC. Following the definition of the Burkert radius r₀, we identify the one-halo-scale-length r_B of the fermionic model such that ρ(r_B ) = ρ_p /4. We thus calculate the product ρ_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 Figure 5. The shown candidates include isothermal-like (blue outlined points) and nonisothermal (green circles) solutions. For comparison, the results are amended by the MW solution following the fermionic RAR model (Argüelles et al. 2018).

Zoom In Zoom Out Reset image size

Although absolute magnitude information is incomplete, all of the predicted DM surface densities are within the range of the 3σ area as well. This is visualized by a histogram for the full sample (dark gray bars, 120 galaxies) with comparison to the subsample (green bars, 44 galaxies) including nonisothermal solutions (i.e., W_p < 10). Considering the subsample only, we obtain a mean surface density of about 148.5 M_⊙ pc⁻², fully inside the 1σ uncertainty in Equation (21).

The SDR given by Equation (21) is qualitatively consistent with the scaling relation ${M}_{h}\sim {r}_{h}^{2}$ as given by Equation (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 nonlinear 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 ∼ r_B , we obtain ${\rho }_{p}{r}_{B}\approx \mathrm{const}$ (see also Argüelles et al. 2019).

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

On halo scales, i.e., for r ∼ 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 toward a finite mass. When applied to galactic halos, such a transition can show different behaviors 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 κ.

In the fermionic model, the sharpness of the transition can be quantified as follows: for W_p ≪ 1, the fermionic density profiles are polytropes of index n = 5/2 as clearly shown in Figure 6 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 ≳ 10), the fermionic halos become isothermal-like, with ρ(r) ∝ r⁻² 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 ≈ 7 (see Figure 6).

Zoom In Zoom Out Reset image size

These results agree with those obtained in Section 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 labeled through the equivalent parameter k.

Statistically, we identified W_p ≈ 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 an MEP are thermodynamically and dynamically stable for W_p ≤ 7.45 if they are in the dilute regime (θ₀ ≪ −1 where W₀ ≡ W_p ), while in the core-halo regime (θ₀ > 10) the same stability conditions hold for a bounded ${W}_{p}\in [{W}_{p}^{\min },{W}_{p}^{\max }]$. For typical SPARC galaxies with a halo mass of ∼5 × 10¹⁰ M_⊙, we find ${W}_{p}^{\max }\sim {10}^{-2}$.

The point of stability-change for the classical King model was first obtained in Chavanis et al. (2015b), and further rederived 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 ∼5 × 10¹⁰ M_⊙ and with appreciable escape of particles (i.e., W_p ≪ 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.

We show in this section a detailed χ² 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 UGC 05986, 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 Figure 7. In the same plot for UGC 05986, 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 ≲ 10), as can be explicitly seen through the χ² valleys in top panels of Figure 8. 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 2 orders of magnitude (see middle panel of first row in Figure 8).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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² ≈ 50 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}^{\mathrm{cr}}\approx 2\times {10}^{8}{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 a negligible escape of particles (W_p ≳ 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 toward a maximum without a clear turning point compared to the first group; see central plots in Figure 7. Fitting those galaxies for different W_p values does not favor solutions with or without escaping particles effects. The variation in the χ² value remains rather small; see middle panels of Figure 8.

Finally, the third group, represented by NGC 6015, shows some oscillations in the RC, mainly in the outer halo; see right plots in Figure 7. 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 ≳ 10)—although being far from good—see bottom panels of Figure 8. Such solutions provide a wide halo maximum followed by a flat RC. In contrast, the nonisothermal solutions with a cutoff provide only a narrow maximum in the halo, followed by a Keplerian decreasing tail.

It is worth recalling 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.

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 favorite because all are similarly good (see Figure 2). For instance, for all DM models, we obtain nearly identical values for ${{\mathfrak{a}}}_{0}$ as required in Equation (20). This result is in line with a recent analysis done in Kaplinghat et al. (2020).

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 Figure 3). This is not only in agreement with a similar analysis done by Li et al. (2020) but also 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 disfavored) 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 Figure 4). 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.

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 subsample 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 ≪ 1), transition into the 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 ≫ 10).

Interestingly, the mean κ 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 Figure 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 subgroup (with W_p ≪ 1) are thermodynamically and dynamically stable, with an outer-halo morphology of polytropic nature (see Figure 6). This may be evidencing a fundamental and deep link between thermodynamics of self-gravitating fermions and galaxy formation and morphology.