Fermionic Dark Matter: Physics, Astrophysics, and Cosmology

The nature of dark matter (DM) is one of the most relevant questions in modern astrophysics. We present a brief overview of recent results that inquire into a possible fermionic quantum nature of the DM particles, focusing mainly on the interconnection between the microphysics of the neutral fermions and the macrophysical structure of galactic halos, including their formation both in the linear and non-linear cosmological regimes. We discuss the general relativistic Ruffini-Arg\"uelles-Rueda (RAR) model of fermionic DM in galaxies, its applications to the Milky Way, the possibility that the Galactic center harbors a DM core instead of a supermassive black hole (SMBH), the S-cluster stellar orbits with an in-depth analysis of the S2's orbit including precession, the application of the RAR model to other galaxy types (dwarf, elliptic, big elliptic and galaxy clusters), and universal galaxy relations. All the above focusing on the model parameters constraints, most relevant to the fermion mass. We also connect the RAR model fermions with particle physics DM candidates, self-interactions, and galactic observables constraints. The formation and stability of core-halo galactic structures predicted by the RAR model and their relation to warm DM cosmologies are also treated. Finally, we briefly discuss how gravitational lensing, dynamical friction, and the formation of SMBHs can also probe the DM nature.


Introduction
The main evidence for the existence of dark matter (DM) is implied by its gravitational 17 effects in a plethora of astrophysical and cosmological environments, including the cosmic 18 microwave background (CMB), baryon acoustic oscillations (BAO), galactic structures (its 19 formation, evolution, and morphology), gravitational lensing, stellar streams, and many others. 20 However, understanding its nature and precise overall mass distribution in galactic structures 21 within a particle DM paradigm are still open questions. 22 Several attempts have been made to explain this phenomenon through ordinary matter, 23 including active neutrinos [1,2] or macroscopic objects such as MACHOS [3]. However, a 24 microscopic origin of the DM particles regarding a new particle species not included in the 25 Standard Model remains the most likely hypothesis [4,5].
More recent cosmological observations obtained in the last three decades have favored 27 the adoption of the ΛCDM paradigm [6,7]: in the standard scheme, the DM is assumed to 28 be produced in thermal equilibrium via weak interactions with the primordial plasma and 29 modeled as collisionless after its decoupling from the other particle species. In this scenario, 30 the DM decoupling is assumed to occur at a temperature smaller than the DM rest mass, so 31 the distribution corresponds to non-relativistic particles. The traditional DM candidate within 32 this paradigm is the so-called weakly interacting massive particle (WIMP), typically a heavy 33 neutral lepton with masses around 100 GeV/c 2 [5]. Nonetheless, many other DM candidates 34 exist either of bosonic or fermionic nature, such as axion-like particles or sterile neutrinos, 35 respectively, with somewhat different early Universe decoupling regimes relative to WIMPs, 36 though still in agreement with current cosmological observables (see, e.g., [8,9] for extensive 37 reviews in the case of axions and sterile neutrinos respectively). 38 Besides their different early Universe physics, these different DM particles imply different 39 outcomes in the non-linear regime of structure formation, which may allow favoring one 40 candidate over the other. These non-linear processes typically correspond to the gravitational 41 collapse of primordial self-gravitating DM structures such as DM halos, in which the quantum 42 nature of the DM particles (either bosonic or fermionic) can cause distinguishable patterns 43 which data can test. Their effects on the precise shape and stability of the DM density profiles, 44 including the distinctive quantum effects (e.g., quantum pressure) through the central regions 45 of the halos, may open new important avenues of research in the field. Indeed, in [10], it has 46 been shown that these interesting quantum imprints in the DM halos exist for bosonic DM, and 47 in [11] for fermionic DM. 48 A key motivation to include the quantum nature of the particles in the study of DM 49 halos is the limited inner spatial resolution obtained within cosmological N-body simulations 50 for such formed structures. It encompasses huge uncertainties in (1) the central DM mass 51 distributions; (2) the relation/effects with the supermassive black hole (SMBH) at the center of 52 large galaxies; and (3) the relationship with the baryonic matter on inner-intermediate halo 53 scales, among others. Indeed, using classical (rather than quantum) massive pseudo-particles 54 as the matter building blocks in N-body simulations within the ΛCDM cosmology does not 55 allow testing quantum pressure effects in DM halos. 56 Increasing attention has been given in the last decade to the study of DM halos in terms of 57 quantum particles, given they may alleviate/resolve many of the drawbacks still present in the 58 traditional CDM paradigm on small-scales (i.e., typically below ∼ 10 kpc scales). The bulk of 59 these models is comprised of the following three categories: 60 (i) Ultralight bosons with masses 1 m b c 2 ∼ 1-10 × 10 −22 eV, known as ultralight DM, Fuzzy 61 DM, or even scalar field DM [12][13][14][15][16][17][18][19][20][21][22]. 62 (ii) The case of fully-degenerate fermions (i.e., in the zero temperature approximation 63 under the Thomas-Fermi approach) with masses m df c 2 ∼ few ×10 2 eV [23 -26]. Or the case 64 of self-gravitating fermions but distributed in the opposite limit, that is, in the dilute regime 65 (i.e., in Boltzmannian-like fashion) which, however, do not imply an explicit particle mass 66 dependence when contrasted with halo observables (see, e.g., [27]). 67 (iii) The more general case of self-gravitating fermions in a semi-degenerate regime (i.e., 68 at finite temperature), which can include both regimes in the same system: that is, to be highly 69 degenerate in the center and more diluted in the outer region (see [26,[28][29][30][31] for a list of generic 70 works). Recently, the phenomenology of this theory to the study of DM in real galaxies (using 71 specific boundary conditions from observations), was developed in full general relativity either 72 including for escape of particles [32][33][34][35][36][37] or not ( [38]), and leading to particle masses in the range 73 m f c 2 ∼ few 10-100 keV. The latter model is usually referred to in the literature as the Ruffini- 74 Argüelles-Rueda (RAR) model (it has been sometimes called the relativistic fermionic-King 75 model). 76 A relevant aspect of the above models is the DM particle mass dependence on the density 77 profiles, which differs from phenomenological profiles used to fit results from classical N-body 78 numerical simulations. Moreover, these type of self-gravitating systems of DM opens the 79 possibility of having access to the very nature, mass, and explicit dependence on their phase-80 space distributions at the onset of DM halo formation in real galaxies (see, e.g., [20] and refs. 81 therein, for bosons and, e.g. [11,26,33] and refs. therein, for fermions). 82 In this review, we will focus on fermionic candidates as described in (iii), describe the 83 huge progress in the last decade in their theory and phenomenology, and comment on future 84 perspectives. This work is organized as follows. In Section 2, we recall the main features of 85 the RAR model, the theoretical framework including the fermion equation of state (EOS), the 86 general relativistic equations of equilibrium, and the properties of the solutions such as the DM 87 density profile and galactic rotation curves. Section 3 is devoted to applying the RAR model 88 to the description of the MW, so the constraint that the Galactic data of the rotation curves 89 impose on the fermion mass. Special emphasis is also given to the high-quality data of the 90 motion of the S-cluster stars in the proximity of Sgr A*. We outline in Section 4 the application 91 of the RAR model to other galaxy types, specifically to dwarf spheroidal, spiral, and elliptical 92 galaxies, as well as galaxy clusters. The extension of the model to various galaxies allows us to 93 verify its agreement with universal observational correlations. Section 5 discusses this topic. 94 Section 6 discusses the possibility of the sterile neutrinos being the massive fermions of the 95 RAR model, including the constraints on possible fermion self-interactions. Section 7 discusses 96 the formation and stability of the fermionic DM from a cosmological evolution viewpoint, e.g., 97 the non-linear structure and the relation with warm DM cosmologies. Section 8 outlines some 98 additional astrophysical systems and observations that can help to constrain DM models, e.g., 99 the gravitational lensing, the DM dynamical friction on the motion of binaries, and the link 100 between DM halos and SMBHs at the centers of active galaxies. Finally, Section 9 concludes. 101 2. The RAR model: theoretical framework 102 The RAR model was proposed to evaluate the possible manifestation of the dense quantum 103 core -classical halo distribution in the astrophysics of real galaxies, i.e., as a viable possibility to 104 establish a link between the dark central cores to DM halos within a unified approach in terms 105 of DM fermions [38]. The RAR model equilibrium equations consist of the Einstein equations 106 in spherical symmetry for a perfect fluid energy-momentum tensor. The Fermi-Dirac statistics 107 gives the pressure and density, while the closure relations are determined by the Klein and 108 Tolman thermodynamic equilibrium conditions [38]. The solution to this system of equations 109 leads to continuous and novel profiles for galactic dark matter halos, whose global morphology 110 depends on the fermionic particle mass. Such morphology has a universal behavior of the type 111 dense core -dilute halo that extends from the center to the galactic halo, which allows providing 112 solutions to various tensions faced by standard cosmological paradigms on galactic scales. The 113 outermost part of such distributions makes it possible to explain the galactic rotation curves 114 (in a similar way as traditional dark matter profiles do), while their central morphology is 115 characterized by high concentrations of semi-degenerate fermions (due to the Pauli exclusion 116 principle) with important astrophysical consequences for galactic nuclei (see [39][40][41] for its 117 applications). Similar core-halo profiles with applications to fermionic DM were obtained in 118 [30] and more recently in [31] from a statistical approach within Newtonian gravity.

119
The above corresponds to the original version of the RAR model, with a unique family 120 of density profile solutions that behaves as ρ(r) ∝ r −2 at large radial distances from the 121 center. This treatment was extended in [32] by introducing a cutoff in momentum space in the 122 distribution function (DF) (i.e., accounting for particle-escape effects) that allows defining the 123 galaxy border. The extended RAR model conceives the DM in galaxies as a general relativistic 124 self-gravitating system of massive fermions (spin 1/2) in hydrostatic and thermodynamic 125 equilibrium. Following [32], we solve the Tolman-Oppenheimer-Volkoff (TOV) equations 126 using an equation of state (EOS) that takes into account (i) relativistic effects of the fermionic 127 constituents, (ii) finite temperature effects, and (iii) particle escape effects at large momentum 128 (p) through a cutoff in the Fermi-Dirac distribution f c , where = c 2 p 2 + m 2 c 4 − mc 2 is the particle kinetic energy, c is the cutoff particle energy, µ 130 is the chemical potential from which the particle rest-energy is subtracted, T is the temperature, 131 k B is the Boltzmann constant, c is the speed of light, and m is the fermion mass.

132
This extended version of the RAR model differs from the original one presented in [38] 133 only in condition (iii). The inclusion of the cutoff parameter allows a more realistic description 134 of galaxies since it allows us (1) to model their finite size and (2) to apply relaxation mechanisms 135 that realize in nature the Fermi-Dirac distribution in Eq. (1), where the entropy can reach a 136 maximum, unlike the original RAR model with no particle escape condition, i.e., c → ∞ (see 137 [31] for a discussion). In fact, it is important to highlight that quantum phase-space distribution 138 (Eq. 1) can be obtained from a maximum entropy production principle [29]. It has been shown 139 that there is a stationary solution of a generalized Fokker-Planck equation for fermions that 140 includes the physics of violent relaxation and evaporation, appropriate within the non-linear 141 stages of galactic DM halo structure formation [29]. The full set of dimensionless parameters of 142 the model are where β, θ, and W are the temperature, degeneracy, and cutoff parameters, respectively. We do 144 not consider anti-fermions because the temperature fulfills T mc 2 /k B .

145
The stress-energy tensor is that of fermionic gas modeled as a perfect fluid,

146
T αν = diag c 2 ρ(r), P(r), P(r), P(r) , 147 whose density and pressure are associated with the distribution function f c , where the integration is carried out over the momentum space bounded by ≤ c .

149
The system is considered to be spherically symmetric, so we adopt the metric 150 ds 2 = e ν(r) dt 2 − e λ(r) dr 2 − r 2 dθ 2 + r 2 sin 2 θdφ 2 , where (r,θ,φ) are the spherical coordinates, the metric functions e ν(r) and e −λ(r) = 1 − 2M(r)/r 151 are functions of r, and M(r) is the mass function. The Tolman [42] We set the constants on the right-hand side of Eqs. (7a)-(7c) by evaluating the equations 155 at the boundary radius, say r b . For instance, the escape energy condition (7c) becomes being . This condition approaches the Newtonian escape 157 velocity condition, v 2 e = −2Φ, in the weak-field, non-relativistic limit, c → ∞ and e ν/2 ≈ 158 1 + Φ/c 2 , where Φ is the Newtonian gravitational potential, and setting Φ(r b ) = 0. 159 Therefore, the Einstein equations together with the equilibrium conditions (7a)-(8) form 160 the following system of coupled non-linear ordinary integrodifferential equations where the subscript '0' stands for variable evaluated at r = 0. We have introduced di-162 mensionless quantities:r = r/χ,M = GM/(c 2 χ),ρ = Gχ 2 ρ/c 2 ,P = Gχ 2 P/c 4 , with 163 χ = 2π 3/2 (m Pl /m)h/(mc), being m Pl = √h c/G the Planck mass. Equations (9a) and (9b) 164 are the relevant Einstein equations, Eq. (9c) is a convenient combination of the Klein and 165 Tolman relations for the gradient of θ = µ/(k B T), and Eq. (9e) is a direct combination of the 166 Klein and cutoff energy conditions. Note that in the limit W → ∞ (no particle escape: c → ∞), 167 the above equation system reduces to the equations obtained in the original RAR model [38]. 168 Argüelles et al. [32] have shown that including the cutoff parameter, besides its theoretical 169 relevance, allows handling more stringent outer halo constraints and leads to central cores 170 with higher compactness compared to the ones produced by the unbounded solutions of the 171 original RAR model [38].

173
Therefore, in the RAR model, the distribution of DM in the galaxy is calculated self-174 consistently by solving the Einstein equations, subjected to thermodynamic equilibrium condi-175 tions. The system of equations must be subjected to boundary conditions to satisfy galactic 176 observables. In this section, we summarize the most relevant example, the case of the Milky 177 Way (MW).

The Milky Way rotation curves 179
Thanks to the vast amount of rotation curve data (from the inner bulge to outer halo) [45], 180 our Galaxy is the ideal scenario to test the RAR model. The observational data used in [45] 181 varies from a few pc to a few 100 kpc, covering various orders of magnitude of the radial extent 182 with different baryonic and dark mass structures, and can go down to the ∼ 10 −4 pc scale 183 when including the S-cluster stars [46]. To fit the observed rotation curve, and according to 184 [45], the following matter components of the MW must be assumed: A central region (r ∼ 10 −4 -2 pc) of young stars and molecular gas whose dynamic is 186 dictated by a dark and compact object centered at Sgr A*.

187
ii). An intermediate spheroidal bulge structure (r ∼ 2-10 3 pc) composed mostly of older 188 stars, with inner and main mass distributions explained by the exponential spheroid 189 model.

190
iii). An extended flat disk (r ∼ 10 3 -10 4 pc) including star-forming regions, dust, and gas, 191 whose surface mass density is described by an exponential law.

194
Considering the previous mass distributions for the Galaxy, the extended RAR model was 195 successfully applied to explain the MW rotation curve as shown in Figs. 1 and 2, implying a 196 more general dense core -diluted halo behavior for the DM distribution: A DM core with radius r c (defined at the first maximum of the twice-peaked rotation 198 curve), whose value is shown to be inversely proportional to the particle mass m, in which 199 the density is nearly uniform. This central core is supported against gravity by the fermion 200 degeneracy pressure, and general relativistic effects are appreciable.

205
• Finally, the DM density reaches a Boltzmann regime supported by thermal pressure with 206 negligible general relativistic effects. It shows a behavior ρ ∝ r −n with n > 2 that is due 207 to the phase-space distribution cutoff 2 . This leads to a DM halo bounded in radius (i.e., 208 ρ ≈ 0 occurs when the particle escape energy approaches zero).

209
The different regimes in the density profiles are also revealed in the DM rotation curve 210 showing (see right panel of Fig. 1): A linearly increasing circular velocity v ∝ r reaching a first maximum at the quantum core 212 radius r c .

213
• a Keplerian power law, v ∝ r −1/2 , with decreasing behavior representing the transition 214 from quantum degeneracy to the dilute regime. After a minimum, highlighting the 215 plateau, the circular velocity follows a linear trend until reaching the second maximum, 216 which is adopted as the one-halo scale length in the fermionic DM model.

217
• A decreasing behavior consistent with the power-law density tail ρ ∝ r −n due to the cutoff 218 constraint.  The Galaxy data analysis allowed us to rule out the fermion mass range mc 2 < 10 keV 221 because the corresponding rotation curve starts to exceed the total velocity observed in the 222 baryon-dominated region. On the other hand, by focusing only on the quantum core, it was 223 possible to derive constraints that further limit the allowed fermion mass to mc 2 ≈ 48-345 224 keV. The dynamics of the S-cluster give the mass a lower bound. The S-stars analysis made 225 through a simplified circular velocity analysis in general relativity showed that the fermion 226 mass should be mc 2 ≥ 48 keV. Namely, the quantum core radius of the solutions for mc 2 < 48 227 keV are always greater than the radius of the S2 star pericenter, i.e., r c > r p (S2) = 6 × 10 −4 228 pc, which rules out fermions lighter than 48 keV. The mass upper bound of mc 2 = 345 keV 229 corresponds to the last stable configuration before reaching the critical mass for gravitational 230 collapse, M cr c ∝ m 3 Pl /m 2 [32,48,49].

231
As was explicitly shown in [32,33,50], this new type of dense core -diluted halo density 232 profile suggests that the DM might explain the mass of the dark compact object in Sgr A* as 233 well as the halo mass. It applies not only to the MW but also to other galactic structures, from 234 dwarfs and ellipticals to galaxy clusters [33]. Specifically, an MW analysis [32] has shown 235 that this DM profile can indeed explain the dynamics of the closest S-cluster stars (including 236 S2) around Sgr A*and to the halo rotation curve without changing the baryonic bulge-disk 237 components (see Figs. 1 and 2). Therefore, for a fermion mass between 48-345 keV, the RAR 238 solutions explain the Galactic DM halo and, at the same time, provide a natural alternative to 239 the central BH scenario.  Fig. 1. The graphic allows us to appreciate better the difference between the diverse DM models in the radial window where the Rotation Curve is most relevant. Reprinted from [32], Copyright (2018), with permission from Elsevier.

241
The extensive and continuous monitoring of the closest stars to the Galactic center has 242 produced, over decades, a large amount of high-quality data on their positions and velocities. 243 The explanation of these data, especially the S2 star motion, reveals a compact source, Sagittar-244 ius A* (Sgr A*), whose mass must be about 4 × 10 6 M . This result has been the protagonist 245 of the Nobel Prize in Physics 2020 to Reinhard Genzel and Andrea Ghez for the discovery of a 246 supermassive compact object at the centre of our galaxy. Traditionally, the nature of Sgr A* has been 247 attributed to a supermassive black hole (SMBH), even though direct proof of its existence is ab-248 sent. Further, recent data on the motion of the G2 cloud show that its post-peripassage velocity 249 is lower than expected from a Keplerian orbit around the hypothesized SMBH. An attempt 250 to overcome this difficulty has used a friction force, produced (arguably) by an accretion flow 251 whose presence is also observationally unconfirmed. 252 We have advanced in Argüelles et al. [32] and in Becerra-Vergara et al. [34,35], Argüelles 253 et al. [51] an alternative scenario that identifies the nature of the supermassive compact object 254 in the MW center, with a highly concentrated core of DM made of fermions (referred from 255 now on as darkinos). The existence of a high-density core of DM at the center of galaxies had 256 been demonstrated in Ruffini et al. [38], where it was shown that core-halo profiles are obtained 257 from the RAR fermionic DM model. The DM galactic structure is calculated in the RAR model 258 treating the darkinos as a self-gravitating system at finite temperatures, in thermodynamic 259 equilibrium, and in general relativity. It has already been shown that this model, for darkinos of 260 48-345 keV, successfully explains the observed halo rotation curves of the MW [32] and other 261 galaxy types [33]. 262 Therefore, since 2020 we move forward by performing first in Becerra-Vergara et al. 263 [34] and then in Becerra-Vergara et al. [35], Argüelles et al. [36,51], observational tests of the 264 theoretically predicted dense quantum core at the Galactic center within the DM-RAR model. 265 Namely, to test whether the quantum core of darkinos could work as an alternative to the SMBH 266 scenario for SgrA*. For this task, the explanation of the multiyear accurate astrometric data of 267 the S2 star around Sgr A*, including the relativistic redshift that has recently been verified, is 268 particularly important. Another relevant object is G2, whose most recent observational data, as 269 we have recalled, challenge the scenario of an SMBH.   For both objects, it was found that the RAR model fits the data better than the BH scenario. 278 The reduced chi-squares of the time-dependent orbit and z data are χ 2 S2,RAR ≈ 3.1 and 279 χ 2 S2,BH ≈ 3.3 for S2 and χ 2 G2,RAR ≈ 20 and χ 2 G2,BH ≈ 41 for G2. The fit of the z data 280 shows that, while for S2 the fits are comparable, i.e.χ 2 z,RAR ≈ 1.28 andχ 2 z,BH ≈ 1.04, for G2 281 only the RAR model fits the data:χ 2 z,RAR ≈ 1.0 andχ 2 z,BH ≈ 26. Therefore, the sole DM core, 282 for 56 keV fermions, explains the orbits of S2 and G2. No drag force nor other external agents 283 are needed, i.e., their motion is purely geodesic.

284
The above robust analysis (and detailed in [34]) has tested the extended RAR model [32] 285 with the precise astrometric data of the S2 star (and also of G2), showing an excellent agreement 286 for a particle mass of mc 2 = 56 keV. Thus, it is confirmed that the elliptic orbits of these two 287 objects moving in the surroundings of the DM-core made of 56 keV fermions are compact 288 enough so that the core radius r c is always smaller than the pericenter of the stars and have 289 their foci coinciding with the Galactic center as obtained by observations (see, e.g., [52]). The 290 above results invalidate recently claimed drawbacks of the RAR model raised in [53] since 291 none of the working hypotheses therein apply to the RAR model. More recently, in Becerra-Vergara et al. [35], we have extended the above results to all the 294 best observationally resolved S-cluster stars, namely to the up-to-date astrometry data of the 17 295 S-stars orbiting Sgr A*, achieving to explain the dynamics of the S-stars with similar (and some 296 cases better) accuracy compared to a central BH model (see Fig. 5).
297 Table 1 in Becerra-Vergara et al. [35] summarizes the best-fit model parameters and reduced-298 χ 2 for the position and the line-of-sight radial velocity of the 17 S-cluster stars for the central 299 BH and the RAR model. The average reduced-χ 2 of the RAR model was 1.5741 and the 300 corresponding value of the central BH model was 1.6273.

301
Therefore, a core of fermionic DM at the Galactic center explains the orbits of the S-stars 302 with similar accuracy compared to a central BH model. The same core -halo distribution of 56 303 Figure 5. Best-fit orbits for the 17 best-resolved S-star orbiting Sgr A*. It shows the projected orbit in the sky, X vs. Y, where X is right ascension and Y is declination. The black dashed curves correspond to the BH model and the colored curves to the RAR model for mc 2 ≈ 56 keV fermions. Reproduced from [35] with the authors' permission.  keV fermions also explains the MW rotation curves [32,34]. Data of the motion of objects near 304 Sgr A*, if accurate enough, could put additional constraints on the fermion mass. The recently 305 detected hot-spots apparently in a circular motion at 7-23 GM BH /c 2 radius [54,55], however, fail 306 in this task because the highly-model dependence of the objects real orbit inference, including 307 the lack of information on the object's nature. In addition, the quality of their astrometry data 308 is low relative to the S-stars data. We hope the data quality of these spots or similar objects can 309 increase so that they can put relevant constraints on Sgr A* models. We keep our eyes open to 310 the data of newly observed S-stars (S62, S4711-S4714). Their motion models based on a central 311 BH predict they have pericenter distances ∼ 400 GM BH /c 2 [56,57]. If confirmed, these new 312 S-stars might constrain the DM core size, so the fermion mass. Furthermore, [36] focused on the periapsis precession of the S2 orbit. We have there 315 quantified, for the first time within the RAR-DM model, the effects on the S2-star periapsis 316 (precession) shift due to an extended DM mass filling the S2 orbit, in contrast with the vacuum 317 solution of the traditional Schwarzschild BH. The main result is as follows. While the central BH 318 scenario predicts a unique prograde precession, in the DM scenario, it can be either retrograde 319 or prograde, depending on the amount of DM mass enclosed within the S2 orbit, which in turn 320 is a function of the fermion mass (see Fig. 6 and Table 6).

321
Therefore, for larger and larger particle masses, the behavior of the RAR model tends to be 322 that of the BH. As Fig. 6 shows, the precession tends nearly asymptotically to 12 arcmin (≈ 0.2 323 deg), the predicted precession of a central Schwarzschild BH (see, e.g., [58]). The RAR model 324 produces the same prograde precession of a central BH for 345 keV fermions, corresponding to 325 the unstable DM core for gravitational collapse into a BH. For 56.4 keV fermions, the prograde 326 and retrograde contributions balance each other leading to a zero net precession. For lower 327 masses, the net precession is retrograde. It has been shown in [36] that currently, available 328 data constrains the amount of retrograde precession, imposing a lower limit to the fermion 329 mass of ≈ 57 keV, hence an upper limit to the amount of mass enclosed in the orbit of ≈ 0.1% 330 of the core's mass (see Table 1). Indeed, the latter limit agrees with the one obtained by [58]. 331 For fermion masses above 57 keV, the prograde precession of the RAR model and that of a 332 central BH agree within the experimental uncertainties. The reason has been explained in [36]: 333 Table 1. Comparison of the BH and RAR DM models that best fit all the publicly available data of the S2 orbit. The 2nd column shows the central object mass, M CO . For the Schwarzschild BH model, M CO = M BH , while for the RAR model, M CO = M c , with M c the DM core mass. The 3rd column shows the radius of the central object, r c . For the Schwarzschild BH model, r c is given by the event horizon radius, R Sch = 2GM BH /c 2 . The 4th column shows the DM mass enclosed within the S2 orbit, ∆M DM /M CO . The best-fitting pericenter and apocenter radii of the S2 orbit are given in the 5th and 6th columns. The values of the average reduced-χ 2 of the best fits, defined as in 34, are given in the 7th column. The last two columns show the model predictions of the periapsis precession of the real orbit, ∆φ, and of the sky-projected orbit, ∆φ sky . Reproduced from [36] with the authors' permission.  59], while the best place to 334 analyze orbital precession is around the apocenter.

335
The above is shown in Fig. 7 that plots the relativistic precession of S2 projected orbit in a 336 right ascension -declination plane. It can be seen that while the positions in the plane of the sky 337 nearly coincide about the last pericenter passage in the three models, they can be differentiated 338 close to the next apocenter. Specifically, the upper right panel evidences the difference at the 339 apocenter between the prograde case (as for the BH and RAR model with mc 2 = 58 keV) and 340 the retrograde case (i.e., RAR model with mc 2 = 56 keV).

341
The bad news is that, as evidenced by Fig. 7, all the current and publicly available data 342 of S2 can not discriminate between the two models. The good news is that the upcoming S2 343 astrometry data close to the next apocentre passage could potentially establish if a classical BH 344 or a quantum DM system governs Sgr A*.

345
A further interesting consequence of this scenario is that a core made of darkinos becomes 346 unstable against gravitational collapse into a BH for a threshold mass of ∼ 10 8 M . Collapsing 347 DM cores can provide the BH seeds for forming SMBHs in active galaxies (such as M87) without 348 the need for prior star formation or other BH seed mechanisms involving super-Eddington 349 accretion rates as demonstrated in [11] from thermodynamic arguments. This topic is of major 350 interest, and further consequences and ramifications are currently being studied, including:    . Relativistic precession of S2 in the projected orbit on the plane of the sky, predicted in the BH and RAR DM models. While it is prograde for the BH and RAR (mc 2 = 58 keV) (in dashed black and green respectively), it is retrograde for the RAR DM model (mc 2 = 56 keV) (in dashed red). The solid (theoretical) curves and gray (data) points correspond to the first period (≈ 1994-2010), while the dashed (theoretical) curves and cyan (data) points to the second period (≈ 2010-2026). Right panels: zoom of the region around apocenter (top panel) and pericenter (bottom panel). Reproduced from [36] with the authors' permission. Table 2. Typical halo radius and mass for dsPh, spirals, elliptical galaxies, and galaxy cluster, adopted in [32] for testing the RAR model. typical dSph typical spiral typical elliptical typical galaxy cluster 3 × 10 7 1 × 10 12 5 × 10 12 3 × 10 14 types, such as dwarfs, spirals, ellipticals, and galaxy clusters, as detailed in [33]. Thus, for the 366 relevant case of a fermion mass of mc 2 ≈ 50 keV (motivated by the MW phenomenology as 367 shown in 3), we make full coverage of the remaining free RAR model parameters (β 0 ,θ 0 ,W 0 ) and 368 present the complete family of density profiles which satisfy realistic halo boundary conditions 369 as inferred from observations (see Fig. 8).

370
For galaxy types located far away from us, the observational inferences of their DM 371 distributions are limited to a narrow window of galaxy radii, typically from a few up to several 372 half-light radii. Generally, we have no access to observations for possible detection of a central 373 dark compact object (as for SgrA* in our Galaxy) nor to constrain the system's total (or virial) 374 mass. Thus, we adopt here (see [33] for further details) a similar methodology as applied 375 to the MW (as shown in Fig. 3 and detailed in [32]), but only limited to halo-scales where 376 observational data is available (i.e., we do not set any boundary mass conditions in the central 377 core). In particular, we select as the boundary conditions a characteristic halo radius r h with the 378 corresponding halo M h ≡ M(r h ). The halo radius is defined as the location of the maximum 379 in the halo rotation curve, which is defined as the one-halo scale length in the RAR model. 380 We list below the parameters (r h , M h ) adopted for the different DM halos as constrained from 381 observations in typical dwarf spheroidal (dSph), spiral, elliptical galaxies, and bright galaxy 382 clusters (BCG). We only exemplified the first three cases; see [33] for the BCGs. The halo parameters in Table 2 are an average of the eight best-resolved dSphs of the MW 385 studied in [60] by solving the Jeans equations and adopting a cored-Hernquist DM profile, 386 which is similar to the mid-outer region of the RAR profiles [33]. The halo parameters shown in Table 2 are obtained from a group of nearby disk galaxies 389 taken from the THINGS data sample [61], where it was possible to obtain a DM model-390 independent evidence for a maximum in the rotation curves. This was obtained by accounting 391 for baryonic (stars and gas) components in addition to the (full) observed rotation curve from 392 HI tracers (see [33], and refs. therein for details). The halo parameters in Table 2 are obtained from (i) a sample of elliptical galaxies from 395 [62], studied via weak lensing, from which in [63] it was provided the corresponding halo mass 396 models for the tangential shear of the distorted images; and (ii) the largest and closest elliptical 397 M87 as studied in [64], accounting for combined halo mass tracers like stars, globular clusters, 398 and X-ray sources (see [33], and refs. therein for further details). The halo parameters in Table 2 are obtained from a sample of 7 BCGs from [65]. There, 401 the luminous and dark components were disentangled to obtain DM distributions which were 402 reproduced by a generalized NFW (gNFW) model [66], developing a maximal velocity at the 403 one-halo length scale r max(bcg) . In [65], the data among all these 7 cases (i.e., weak lensing and 404 The continuous-magenta curves, occurring only for spiral and elliptical galaxies, indicate the critical solutions which develop compact, critical cores (before collapsing to a BH) of M crit c = 2.2 × 10 8 M . The dashed-magenta curves for dwarfs are limited (instead) by the astrophysical necessity of a maximum in the halo rotation curve. The bounding black solutions correspond to the ones having the minimum core mass (or minimum ρ 0 ), which in turn imply larger cutoff parameters (implying ρ ∝ r −2 when W 0 → ∞. Reprinted from [33], Copyright (2019), with permission from Elsevier. stellar kinematics) support such a maximum within a radial extent from 10 kpc up to 3 Mpc; 405 whose corresponding averages in r h and M h are given in Table 2.

407
The main conclusions of having applied the RAR model for mc 2 ≈ 50 keV to different 408 galaxy types are summarized as follows (see [33] for a more detailed explanation): Typical dwarf galaxies can harbor dense and compact DM cores with masses from 410 M c ∼ 10 3 M up to M c ∼ 10 6 M (see Fig. 8), offering a natural explanation for the 411 so-called intermediate-mass BHs (IMBH). Since the total mass of the typical dSphs here 412 analyzed is below the critical mass of core-collapse (i.e., M tot(d) ∼ 10 7 M < M cr c ≈ 413 2 × 10 8 M ), the core can never become critical and thus will never collapse to a BH. 414 Therefore, the RAR model predicts (for a particle mass of mc 2 ≈ 50 keV) that dSph 415 galaxies can never develop a BH at their center, a result that may explain why these 416 galaxies never become active.

(II)
Typical spiral and elliptical galaxies can harbor denser and more compact DM cores 418 (w.r.t dSphs) with masses from M c ∼ 10 5 M up to M cr c ≈ 2 × 10 8 M (see Fig. 8). 419 Thus, they offer a natural alternative to the supermassive BHs hypothesis (see Fig. 3 420 for the MW). Since the total mass in spirals and ellipticals is much larger than M cr c , the 421 core mass can become critical and eventually collapse towards an SMBH of ∼ 10 8 M 422 which may then grow even larger by accretion.

(III)
Typical bright cluster of galaxies (BCGs) can harbor dense and compact DM cores with 424 masses from M c ∼ 10 6 M up to M cr c ≈ 2 × 10 8 M (see Fig. 8). The implications of 425 this prediction for BCGs are still unclear, mainly given the limited spatial resolution 426 achieved by actual observational capabilities below the central kpc. More work is 427 needed (for example, using strong lensing observations) to evaluate if galaxy clusters 428 show an enhancement in DM density similar to the one predicted by the RAR model. 429 (IV) By combining the range of DM core masses, inner halo densities, and total halo masses 430 as predicted by the RAR model across all of these systems, it is possible to test whether 431 or not this model can answer different Universal scaling relations. A point which is 432 studied in the next section (and further detailed in [33] and [37]).

434
Galaxies follow different scaling relations (or Universal scaling relations) such as the 435 DM surface density relation (SDR) [63], the radial acceleration relation [67], the mass dis-436 crepancy acceleration relation (MDAR) [68] and the Ferrarese relation [69,70] between the 437 total halo mass and its supermassive central object mass. Thus, this section's main goal is to 438 illustrate/summarize the ability of the RAR model to agree with all these relations, covering 439 a really broad window of radial scales across very different galaxy types. The analyses are 440 exemplified for a fermion of rest mass-energy of ≈ 50 keV.

441
For given halo parameters (r h , M h ) as inferred from observations in the smallest up to 442 the largest galaxy types (see Table 2), the RAR model predicts a window of halo values on 443 other radial-scales: (a) the plateau density (or inner-halo constant density) ρ pl ; (b) the total halo 444 mass M tot ; and (c) the fermion-core mass M c . We will then use this predictive power of the 445 RAR model (shown Fig. 8) to test if it can answer for the above Universal scaling relations as 446 reported in the literature. We will focus first on the Ferrarese relation between the total halo 447 mass and its supermassive central object mass and in the SDR.  Fig. 8). The dashed horizontal line represents the Universal relation from the best fit of the data as found by [63]. The dark-gray region indicates the delimited area by the 3 − σ error bars of all the data points. Reprinted from [33], Copyright (2019), with permission from Elsevier.

449
This relation establishes that, for large enough galaxies 3 , the more massive the halo, the 450 larger (in mass) is the supermassive compact object lying at its center. This relation is shown 451 in dashed in the left panel Fig. 9, as taken from [69], together with other more updated 452 versions of such relation such as the one found in [70] (shown in dotted line in the same figure). 453 Interestingly enough, the whole family of RAR solutions, with its corresponding window of 454 predicted core and total DM masses M c , M tot (see left panel of Fig. 9), from typical dwarf to 455 ellipticals galaxies (see green area), do overlap with the Ferrarese relation. That is, the RAR 456 model contains the observational strip (in light blue), while at the same time, it extends out such 457 a Ferrarese strip, indicating a yet-unseen or otherwise unphysical cases (see [33] for further 458 details).

459
Even if no observational data exist yet at the down-left corner of Fig. 9 (left panel), special 460 attention has to be given to the RAR model predictions for dwarf galaxies: recent observations 461 towards the center of some ultra-compact dwarf galaxies of total mass few ∼ 10 7 M (e.g., 462 [72][73][74]), indicate the existence of putative massive BH of few ∼ 10 6 M . Interestingly, the RAR 463 model naturally allows for a slight difference (less than an order of magnitude) between M c 464 and M tot . However, more work is needed on this particular dwarf galaxy to make a more 465 definite statement. The DSR relation establishes that the central surface DM density in galaxies is roughly 468 constant, spanning more than 14 orders in absolute magnitude (M B ): ρ 0D r 0 ≈ 140M pc −2 469 (with ρ 0D the inner-halo -or sometimes called central-DM halo density measured at the Burkert 470 halo radius r 0 ) [63].

471
Since the Burkert central density corresponds to the plateau density of the RAR density 472 profiles (see right panel of Fig. 11), i.e., ρ 0D ≡ ρ pl , and r 0 ≈ 2/3 r h (as shown in [33] connecting 473 the Burkert and RAR one-halo scale-lengths), it is possible to check the DSR by calculating the 474 product 2/3 ρ pl r h along the entire family of astrophysical RAR solutions of Section 4 (including 475 the Milky Way). This is shown in the right panel of Fig. 9.

476
It can be seen that the RAR model predictions agree with the observed relation, the latter 477 displayed in Fig. 9 (right panel) in the dark-grey region delimited (or enveloped) within the 478 3 − σ error bars along all the data points considered in Donato et al. [63]. That figure shows 479 that each typical galaxy type's predicted RAR surface density (vertical solid lines) is within 480 the expected 3 − σ data region. We further notice that typical bright clusters are beyond 481 the observed window as reported by Donato et al. [63], who considered only up to elliptical 482 structures. 483 484 We now confront the RAR model predictions with two closely related Universal relations, 485 though this time not based on typical galaxies with given (averaged) parameters (r h , M h ) 486 inferred from observables (as done before), but from a sample of 120 rotationally supported 487 galaxies taken from the SPARC data set [75].

488
The radial acceleration relation is a non-linear correlation between the radial acceleration 489 exerted on the total matter distribution and the one caused by the baryonic matter component 490 only (see Eq. 10 below). Thus it offers an important restriction to any model that explains the 491 DM halo in galaxies.

492
The different mass components in a galaxy, like a bulge, disk, gas, etc., trace their own 493 contribution to the centripetal or radial acceleration a = v 2 /r. As originally shown in [67], it 494  Figure 10. Radial acceleration relation (top) and mass discrepancy acceleration relation (bottom) for SPARC data and competing DM halo models. Each plot is divided into 50 × 50 equal bins. The baryonic centripetal acceleration a bar is inferred from luminosity observables while the total acceleration a 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 tot = a DM + a bar . The corresponding solid curves are the best fits characterized by a specific a 0 . The histogram plots (upper row) show a Gaussian distribution of log 10 (a 0 /a 0 ). The grayscale legend shows the number of points per bin (2396 for the 120 SPARC-galaxies used). Reproduced from [37].
exists a relatively tight relation between the radial acceleration owing to the total mass, say a tot , 495 and that produced by the baryonic mass component, say a bar , which the empirical formula can 496 well describe with a 0 the only adjustable parameter. Notice that for a bar a 0 , the DM dominates, and the 498 relation deviates from linear, while for a bar a 0 , the baryonic component dominates and 499 the linear relation is recovered (see Fig. 10). The relation has been shown to apply applies 500 to different galaxy types, e.g., disk, elliptical, lenticular, dwarf spheroidal, and low-surface-501 brightness galaxies, so it looks like a real universal law [67,76,77]. 502 We also include in the analysis a close relation: the mass discrepancy acceleration relation 503 (MDAR) relation between baryonic and total mass components. It is defined by D = M tot /M bar 504 with M bar the total baryonic mass and M tot the total mass content of the galaxy including the 505 DM component. If we define the radial acceleration through its definition in terms of the 506 gravitational potential of a spherically symmetric mass distribution, then the MDAR can be 507 written as D = a tot /a bar . In Fig. 10, we show the results of the best fits for the radial acceleration relation and MDAR 509 for both the RAR model (central panels) and the DC14 DM model (right panels) accounting for 510 baryonic-feedback [78] (see also the full Figure 2 in [37] for best-fits to both relations by other 511 DM typical halo models used in the literature). Since both models achieve very close fitted 512 values of a 0 ≈ 1.2 × 10 −10 m s −2 as the one obtained from the SPARC data only (i.e., done in 513 DM model-independent fashion, see left panels), then it is concluded that this Universal relation 514 does not help to prefer one model over the other statistically (see [37] for details). Instead, as 515 demonstrated in [37], the individual fitting of each rotation curve (covering the full SPARC 516 sample studied) allows to favor of the cored profiles over the cuspy Navarro-Frenk-White 517 profile (not shown here).

518
Finally, we have also compared the DM halo models by their density profiles and rotation 519 curves for typical SPARC galaxies in the left panel of Fig. 11. Most models show similarities in 520 the halo tails, but the RAR model is the only one showing a dense-fermion core at the galaxy's 521 center, which acts as an alternative to the SMBH scenario. It is also relevant to mention that 522 while many of the models (such as Einasto or DC14) need the baryonic feedback to produce the 523 statistically favored cored density profiles, the RAR model naturally achieves a cored behavior 524 (i.e., the inner-halo plateau) due to the quasi-thermodynamic equilibrium of the particles at 525 formation (see [37] and references therein for further discussions).

526
Although we do not expect significant qualitative or quantitative differences in the galaxy 527 structure parameters for fermion masses around the above-explored value of 50 keV, studying 528 the implications of different fermion masses for the universal relations will be interesting. We 529 plan to have new results on this topic soon. 6. Fermionic DM and particle physics 531 6.1. Are the sterile neutrinos the fermions of the RAR model?

532
In this section, we will discuss a possible connection to particle physics (i.e., beyond the 533 standard model of particles), analyzing the possibility that the DM particles are self-interacting 534 DM. Namely, DM particles interact among themselves via some unknown fundamental in-535 teraction besides gravity. It is a very active field of research within the DM community since 536 self-interacting DM (SIDM) has been proposed as a possible solution to several challenges 537 which is facing the standard cosmological model on galactic scales (see [41,79] for different 538 reviews on the subject).

539
In 2016, our group presented an extension of the RAR model to include fermion self-540 interactions, referred to as the Argüelles, Mavromatos, Rueda and Ruffini (AMRR) model [80]. 541 In the AMRR model, it has been advanced that the darkinos might be the right-handed sterile 542 neutrinos introduced in the minimum standard model extension (νMSM) paradigm [81]. The 543 AMRR model adopts right-handed neutrinos self-interacting via dark-sector massive (axial) 544 vector mediators.

545
In 2020, [82] explored the radiative decay channel of such sterile neutrinos into X-rays 546 due to the Higgs portal interactions of the νMSM. This work shows that such generalized 547 RAR profiles, including fermion self-interactions, agree with the overall MW rotation curve. In 548 addition, the window of self-interacting DM cross-sections that satisfy the known bullet cluster 549 constraints has been identified.

550
To further constrain the AMRR-νMSM model, an indirect detection analysis has been 551 performed using X-ray observations from the Galactic center by the Nustar mission [82]. 552 Figure 12 summarizes all the observational constraints.

553
It has also advanced a new generation mechanism based on vector-meson decay, able to 554 produce these sterile neutrinos in the early Universe.

555
Summarizing, by considering a DM profile that self-consistently accounts for the particle-556 physics model, the analysis of NuSTAR X-ray data shows how sterile-neutrino self-interactions 557 affect the νMSM parameter-space constraints. The decay of the massive vector field that 558 mediates the self-interactions affects standard production mechanisms in the early Universe. 559 This mechanism might broaden the allowed parameter space compared to the standard νMSM 560 scenario. The formation and stability of collisionless self-gravitating systems are long-standing 564 problems that date back to the work of D. Lynden-Bell in 1967 on violent relaxation and 565 extend to the virialization of DM halos. Such a relaxation process predicts that a Fermi-Dirac 566 phase-space distribution can describe spherical equilibrium states when the extremization of 567 a coarse-grained entropy is reached. In the case of DM fermions, the most general solution 568 develops a degenerate compact core surrounded by a diluted halo. As we have recently shown 569 [32][33][34], the core-halo profiles obtained within the fermionic DM-RAR model explain the galaxy 570 rotation curves, and the DM core can mimic the effects of a central BH. A yet open problem 571 is whether these astrophysical core-halo configurations can form in nature and if they remain 572 stable within cosmological timescales. These issues have been recently assessed in [11]. 573 Specifically, we performed a thermodynamic stability analysis in the microcanonical 574 ensemble for solutions with given particle numbers at halo virialization in a cosmological 575 framework. For the first time, we demonstrate that the above core-halo DM profiles are stable 576 (i.e., maxima of entropy) and extremely long-lived. We find a critical point at the onset of 577 instability of the core-halo solutions, where the fermion-core collapses towards a supermassive 578 black hole. For particle masses in the keV range, the core-collapse can only occur for M vir  Sterile neutrino parameter space limits obtained for Galactic Center observations using the AMRR profiles (continuous red line) when assuming DM production due to self-interactions through a massive vector field mediator. The light-red shaded region above the continuous red line corresponds to the AMRR limits given by X-ray bounds (i.e. indirect detection analysis), while the vertical shaded region below 48 keV indicates the smallest DM mass compatible with S-cluster stars' rotation curve data. The upper shaded region corresponds to production mechanism bounds: non-resonant production under no lepton asymmetry. Other dotted lines refer to several X-ray bounds for different DM halo profiles, including 0-bounce photon analysis. Reprinted from [82], Copyright (2020), with permission from Elsevier. 10 9 M starting at z vir ≈ 10 in the given cosmological framework. This is a key result since 580 the fermionic DM system can provide a novel mechanism for SMBH formation in the early 581 Universe, offering a possible solution to the yet-open problem of how SMBHs grow so massive 582 so fast. We recall the result of Fig. 13, which evidences the existence of a last stable configuration 583 before the core collapse into an SMBH. We refer the reader to [11] for details on this relevant 584 result.

585
This thermodynamic approach allows a detailed description of the relaxed halos from the 586 very center to the periphery, which N-body simulations do not allow due to finite inner-halo 587 resolution. In addition, it includes richer physical ingredients such as (i) general relativity -588 necessary for a proper gravitational DM core-collapse to an SMBH; (ii) the quantum nature of 589 the particles -allowing for an explicit fermion mass dependence in the profiles; (iii) the Pauli 590 principle self-consistently included in the phase-space distribution function -giving place to 591 novel core-halo profiles at (violent) relaxation.

592
Such treatment allows linking the behavior and evolution of the DM particles from the 593 early Universe to the late stages of non-linear structure formation. We obtain the virial halo 594 mass, M vir , with associated redshift z vir . The fermionic halos are assumed to be formed by 595 fulfilling a maximum entropy production principle at virialization. It allows obtaining a most 596 likely distribution function of Fermi-Dirac type, as first shown in Chavanis [29] (generalizing 597 Lynden-Bell results), here applied to explain DM halos. Finally, the stability and typical lifetime 598 of such equilibrium states and their possible astrophysical applications are studied with a 599 thermodynamic approach.

600
For the first time, we have calculated the caloric curves for self-gravitating, tidally-601 truncated matter distributions of O(10) keV fermions at finite temperatures within general 602 relativity. We applied this framework to realistic DM halos (i.e., sizes and masses). With the 603 precise shape of the caloric curve, we establish the families of stable as well as astrophysical 604 DM profiles (see Figs. 13 and 14). They are either King-like or develop a core-halo morphology 605 that fits the rotation curve in galaxies [32,33]. In the first case, the fermions are in the dilute 606 regime and correspond to a global entropy maximum. In the second case, the degeneracy 607 pressure (i.e., Pauli principle) holds the quantum core against gravity and corresponds to a 608 local entropy maximum. Those metastable states are extremely long-lived, and, as such, they 609 are the more likely to arise in nature. Thus, these results prove that DM halos with a core-halo 610 morphology are a very plausible outcome within nonlinear stages of structure formation. The traditional ΛCDM paradigm of cosmology is in remarkable agreement with large-613 scale cosmological observations and galaxy properties. However, there are increasing tensions 614 of the ΛCDM with observations on smaller scales, such as the so-called missing DM sub-halo 615 problem and the core-cusp discrepancy.

616
High-resolution cosmological simulations of average-sized halos in ΛCDM predict an 617 overproduction of small-scale structures significantly larger than the observed number of small 618 satellite galaxies in the Local Group. Moreover, N-body simulations of CDM predict a cuspy 619 density profile for virialized halos, while observations show dSphs having flattened smooth 620 density profiles in their central regions.

621
A possible alternative to alleviate or try to resolve such tensions is to consider warm dark 622 matter (WDM) particles, meaning that they are semi-relativistic during the earliest stages of 623 structure formation with non-negligible free-streaming particle length.

624
WDM models feature an intermediate velocity dispersion between HDM and CDM that 625 results in a suppression of structures at small scales due to free-streaming. If this free streaming 626 scale today is smaller than the size of galaxy clusters, it can solve the missing satellites problem. 627 However, thermally produced WDM suffers from the so-called catch-22 problem when studied 628   Fig. 13 with corresponding fixed total halo mass. Only profiles (1) (resembling a King distribution) and the core-halo one (3) are stable, while profile (2) is thermodynamically unstable. Interestingly, solutions like (3) were successfully applied to explain the DM halo in the MW in [32]. They are stable, extremely long-lived, and fulfill the observed surface DM density relation and the expected value of the DM dispersion velocity. Reproduced from [11] with the authors' permission.
within N-body simulations [83]. Such WDM-only simulations either show unrealistic core sizes 629 for particle masses above the keV range or acquire the right halo sizes for sub-keV masses in 630 direct conflict with phase-space and Lyman-α constraints.

631
Another compelling alternative to collisionless CDM, apart from WDM, is to consider 632 interactions in CDM. This consideration relaxes the assumption that CDM interacts only 633 gravitationally after early decoupling and includes interactions between DM and SM particles, 634 additional hidden particles, or among DM particles. These later models are denominated as 635 self-interacting DM models (SIDM).

636
To shed light on this matter, we have recently provided a general framework for self-637 interacting WDM in cosmological perturbations by deriving from first principles a Boltzmann 638 hierarchy that retains certain independence from an interaction Lagrangian [84]. Elastic inter-639 actions among the massive particles were considered to obtain a more general hierarchy than 640 those usually obtained for non-relativistic (cold DM) or ultra-relativistic (neutrinos) approxi-641 mations. The more general momentum-dependent kernel integrals in the Boltzmann collision 642 terms are explicitly calculated for different field-mediator models, including a scalar or massive 643 vector field.

644
In particular, if the self-interactions maintain the DM fluid in kinetic equilibrium until 645 the fluid becomes non-relativistic, the background distribution function at that moment will 646 switch into a non-relativistic form. This constitutes the scenario known as non-relativistic 647 self-decoupling (a.k.a. late kinetic decoupling). The consequences of this scenario are poorly 648 explored in the literature, and only some preliminary results have been recently obtained 649 within simplified DM fluid approximations [85]. However, more recently [86], this late kinetic 650 decoupling physics was fully explored with self-interactions treated from first principles using 651 interaction Lagrangian (i.e., superseding the fluid approximation), following the formalism 652 developed in [84]. There, it was found that if one imposes continuity of the limiting expressions 653 for the energy density, the non-relativistic distribution function can be found in an analytic 654 expression (see [86] for details). Figure 15 illustrates the effects of self-interactions in the matter power spectrum for the 656 case of a massive scalar field-mediator while including the late kinetic decoupling case. We 657 have used an extended version of the CLASS code incorporating our results for SI-WDM 658 with particle masses in the ∼ keV range. There, we see some of the particular features of the 659 models. With the inclusion of self-interactions, for models with non-relativistic self-decoupling 660 (i.e., late kinetic decoupling), the resulting power spectra may differ significantly from their 661 relativistic counterparts. Indeed, we find that in this regime, the models are "colder" (i.e., as 662 if they correspond to a higher particle mass) and show, even at smaller k values, a distinctive 663 oscillatory pattern (see, e.g., dot-dashed curves in Fig. 15. This increases the small structures for 664 these models, implying that the few keV (traditional) WDM models excluded from phase-space 665 arguments now agree with observations in this new SI-WDM scenario.

666
Most tensions inherent to νMSM WDM models arise from structure formation, i.e., the 667 MW satellite counts and Lyman-α observations: the preferred parameter ranges may under-668 produce small structures and almost rule out the available parameter space. So, including 669 self-interactions can significantly relax the existing bounds on this family of models. The pre-670 diction of these models for the number of MW satellites and the observations of the Lyman-α 671 forest was presented in [88] (see Fig. 17).

672
While in the above paragraphs, we have discussed the effects of self-interacting O(1 − 673 10) keV light fermionic candidates on the linear-structure formation (i.e., its suppression effects 674 on the linear matter power spectrum and its consequences on Lyman-α forest and satellite 675 counts [88]), we will discuss in next its effects on DM halo profiles. For this, we will follow the 676 specific self-interacting model where the darkinos are right-handed neutrinos self-interacting 677 via massive (axial) vector-boson mediators (also considered as a possible case in [88]). The 678 for a vector field SI-WDM model, using a modification to CLASS. We assume the relaxation time approximation and consider two values of the DM particle mass: 1 and 10 keV. Also plotted are the power spectra of CDM and WDM models with DM mass of 1 and 10 keV. All WDM and SI-WDM models consider a nonresonant production scenario (Dodelson-Widrow mechanism, [87]) with T ∼ (4/11) 1/3 T γ . Importantly, the effect of large enough self-interactions increases the amount of small structure formed for these models (see dot-dashed curves), implying that the few keV (traditional) WDM models which were excluded in the past, can now be back to an agreement with observations in this new SI-WDM scenarios. Reproduced from [86], with permission from SNCSC.
main consequences of this model for DM halo structures were studied in [82], for a particle 679 mass mc 2 ≈ 50 keV. In particular, it mainly investigated the consequences of the Milky Way 680 DM halo and the bullet cluster. The massive boson mediator adds a pressure term in the 681 equilibrium equations (see, e.g., [40,82]) which, for normalized interaction constants up to 682 C V ∼ 10 12 (in Fermi constant units), causes no appreciable effects in the Milky Way rotation 683 curve. However, values larger thanC V ∼ 10 13 are ruled out since the additional pressure 684 term is enough to push forward the halo, spoiling the fit to the data, see Fig. 16. Interestingly, 685 values of interaction strengths between our light fermionic candidates ofC V ∼ 10 8 agree with 686 the bullet cluster measurements. Indeed, on cluster scales, it was demonstrated in [82] (see 687 section 3.2 therein) that such interaction constants resolve the tensions between predictions of 688 ΛCDM-based numerical simulations and observations since the corresponding self-interacting 689 DM (SIDM) cross section 4 (σ SIDM ) lies in the expected range [89] 690 0.1 cm 2 g −1 ≤ σ SIDM m ≤ 0.47 cm 2 g −1 . For each point (θ, m) in the parameter space, we consider a self-interacting model under a vector field mediator, with its interaction constant given by σ/m ∼ 0.144C 2 v /m 3 = 0.1 cm 2 g −1 , the upper limit given by Bullet Cluster constraints (see [88] for details). For comparison, we plot the Lyman-alpha bounds for the non-interacting case for a comparable analysis, plus other bounds to the νMSM parameter space for informative purposes, namely X-Ray indirect detection bounds (in blue) and sterile neutrino production bounds (in grey). We also plot the sterile neutrino model compatible with a tentative 3.5 keV DM signal, the subject of debate in recent years, as a purple triangle. The complete list of references of all such bounds can be found in [88]. Reprinted from [88]. Copyright IOP Publishing. Reproduced with permission. Over the last decade, gravitational lensing has become a very effective tool to probe the 694 DM distribution of systems from galaxy to cluster scales. The availability of high-quality HST 695 imaging and integral field spectroscopy with the VLT has allowed high-precision strong lensing 696 models to be developed in recent years, with tight constraints on the inner mass distribution of 697 galaxy clusters (e.g., [90]). By comparing high-resolution mass maps of galaxy cluster cores 698 obtained by such lens models with cosmological hydrodynamical simulations in the LCDM 699 framework, a tension has emerged in the mass profile of the sub-halo population associated 700 with cluster galaxies [91]. Namely, the observed sub-halos appear more compact than those 701 in simulations, with circular velocities inferred from the lens models, which are higher, for a 702 given sub-halo halo mass, than those derived from simulations. It is still unclear whether this 703 tension is due to some limitations of the numerical simulations in the modeling of baryonic 704 physics (back-reaction effects on DM) or rather a more fundamental issue related to the CDM 705 paradigm.

706
In [92], the lensing effects of the core-halo DM distribution were computed. The NFW and 707 nonsingular isothermal sphere (NSIS) DM models and the central Schwarzschild BH model 708 were compared at distances near the Galaxy center. The DM density profiles lead to small 709 deviations of light (of 0.1 arcsec) in the halo part (∼ 8 kpc). However, the RAR profile produces 710 strong lensing effects when the dense core becomes highly compact. For instance, the density 711 distribution of the RAR model for a fermion mass of ∼ 100 keV generates strong lensing 712 features at distances ∼ 0.1 mpc (see [92] for additional details). The DM quantum core has no 713 photon sphere inside or outside but generates multiple images and Einstein rings. Thus, it will 714 be interesting to compare and contrast lensing images produced by the RAR model with the 715 ones produced by phenomenological profiles.

716
It is also interesting to construct the shadow-like images produced by the fermion DM 717 cores of the RAR model and analyze them in light of the EHT data from Sgr A*. Indeed, 718 horizonless objects can be compact, lack a hard surface, and cast a shadow surrounded by a 719 ring-like feature of lensed photons [93]. This has been shown for boson stars for Sgr A* in 720 [94] and further studied by [95] within Numerical Relativity simulations. Using full general 721 relativistic ray tracing techniques [96], we have recently computed relativistic images of the 722 fermionic DM cores of the RAR model under the assumption that photons are emitted from a 723 surrounding accretions disc (self-consistently computed in the given metric) and have shown 724 that there exists a particle mass range in which the shadow feature acquires the typical sizes as 725 resolved for the EHT in Sgr A* (Pelle, Argüelles, et al., to be submitted). The measured orbital period decay of relativistic compact-star binaries (e.g., the famous 728 Hulse-Taylor binary pulsar) has been explained with very-high precision by the gravitational-729 wave emission predicted by general relativity of an inspiraling binary, assuming the binary 730 is surrounded by empty space and within the point-like approximation of the two bodies. 731 However, the presence of DM around the binary might alter the orbital dynamics because of a 732 traditionally neglected phenomenon: DM dynamical friction (DMDF). The binary components 733 interact with their own gravitational wakes produced by the surrounded DM, leading to an 734 orbital evolution that can be very different from the orbital evolution solely driven by the 735 emission of gravitational waves.

736
In [97,98], the effect of the DMDF on the motion of NS-NS, NS-WD, and WD-WD was 737 evaluated. Quantitatively, the crucial parameters are the orbital period and the value of the 738 DM density at the binary location in the galaxy. A comparison among the DMDF produced by 739 the NFW, the NSIS, and the RAR model was presented.
For NS-NS, NS-WD, and WD-WD with measured orbital decay rates, the energy loss by 741 gravitational waves dominates over the DMDF effect. However, there are astrophysically viable 742 conditions for which the two effects become comparable. For example, 1.3-0.2 M NS-WD, 743 1.3-1.3 M NS-NS, and 0.25-0.50 M WD-WD, located at 0.1 kpc, the two effects compete with 744 each other for a critical orbital period in the range 20-30 days for the NFW model. For the RAR 745 model, it occurs at an orbital period ∼ 100 days (see [98] for details). Closer to the Galactic 746 center, the DMDF effect keeps increasing, so the above values for the critical orbital period 747 become shorter. For certain system parameters, the DMDF can lead to an orbital widening 748 rather than a shrinking.

749
The DMDF depends on the density profile, velocity distribution function, and velocity 750 dispersion profile. Therefore, measuring with high accuracy the orbital decay rate of compact-751 star binaries at different Galactic locations and as close as possible to the Galactic center might 752 reveal a powerful tool to test DM models.

753
Indeed, during the peer-reviewing process of this review, Chan and Lee [99] published 754 an analysis of the DMDF contribution to the orbital decay rate of the X-ray binaries formed 755 by stellar-mass BHs and stellar companions. Specifically, they analyzed the case of A0620-00 756 and XTE J1118+480, finding that the DMDF can explain their fast orbital decay. The analysis in 757 [99] uses a phenomenological DM profile, so it will be interesting to apply the DMDF with the 758 RAR model, as done in [98], to check whether or not X-ray binaries can further constrain the 759 fermion mass. The explanation of the formation, growth, and nature of SMBHs observed at galaxy centers 762 are amongst the most relevant problems in astrophysics and cosmology. Relevant questions 763 yet waiting for an answer are: how can they increase their size in a relatively short time to 764 explain the farthest quasars [100,101]?; where do the BH seeds come from and how large they 765 must be to form SMBHs of ∼ 10 8 -10 9 M at high cosmological redshift [102]?; and what is the 766 connection between the host galaxy and central SMBH masses [ Recently, we have proposed in Arguelles et al. (submitted) an SMBH formation channel 782 conceptually different from cases (I) and (II). The new scenario is based on the gravitational 783 collapse and subsequent growth of dense fermionic DM cores. The cores originate at the center 784 of the halos and start to grow as they form. As we have recalled in this work, those dense core-785 diluted halo DM density distributions are predicted by maximum entropy production principle 786 models of halo formation [11,113]. We have recalled how the core made of fermions of ∼ 50 787 keV rest mass-energy becomes unstable against gravitational collapse into a BH for a threshold 788 mass of ∼ 10 8 M . Therefore, the new scenario produces larger BH seeds that comfortably 789 grow to SMBH mass values of 10 9 M in a relatively short time without super-Eddington 790 accretion (Arguelles et al., submitted).

792
Possibly one of the most relevant steps has been to develop a framework that allows 793 obtaining the DM density (and other related physical properties) profiles from first principles. 794 The distribution of DM in the galaxy is given by the solution of the Einstein equations subjected 795 to appropriate boundary conditions to fulfill observational data. The nature of the particle 796 candidate sets the equation of state, so the energy-momentum tensor. It turns out that neutral 797 massive fermion DM particles distribute throughout the galaxy following a dense core-dilute 798 halo density profile. See Section 2 for details.

799
Thanks to the above, the nature of the DM particle, like the rest mass, can be constrained 800 using stellar dynamics data, e.g., the MW rotational curves data and the most accurate data of 801 the S-cluster stellar orbits around Sgr A*. See Section 3 for details.

802
A remarkable and intriguing result is that the orbits of the S-stars are correctly described 803 by the presence of the dense core of fermions without the need for the presence of a massive BH 804 at the MW center. We have shown that even the most accurate available data of the S2 star can 805 not distinguish between the two models since both models describe the data with comparable 806 accuracy. However, we have shown that the data of the periapsis precession of the S2 star in 807 the forthcoming three years (i.e., by 2026) could help to discriminate the two scenarios. See 808 Section 3 for details.

809
In addition to the MW, the RAR model describes the DM component observed in other 810 galaxy types, namely dSphs, ellipticals, and galaxy clusters. This point includes the explanation 811 of existing observational Universal relations and the comparison with other DM models 812 regarding the same observables. See Section 4 for details.

813
All the above constrains the DM particle nature limiting the fermion mass to the ∼ 50-350 814 keV range, constituting a relevant starting point for complementary particle physics analyses. 815 Among the plethora of particle-physics candidates, the beyond-standard model sterile neutrinos 816 (also including self-interactions) remain an interesting candidate for being the fermions of the 817 RAR model. See Section 5 for details.

818
None DM discussion is complete without examining the cosmological implications. In this 819 order of ideas, it has been first shown that the core-halo configurations can indeed form in the 820 Universe in appropriate cosmological timescales thanks to the mechanism of violent relaxation 821 that predicts equilibrium states described by the Fermi-Dirac phase-space distribution. Second, 822 it has been shown that those equilibrium states are long-lasting, with a lifetime of several 823 cosmological timescales, so well over the Universe's lifetime. Third, it has been shown that 824 those small-scale DM core-halo substructures can form from non-linear cosmological density 825 perturbations during the cosmological evolution. See Section 6 for details.

826
Last but not least, we have discussed some additional theoretical and observational 827 scenarios that can help to probe DM models. Specifically, if DM permeates galaxies, stellar 828 objects do not sit in a perfect vacuum. Thus, the distribution of DM can cause dynamical 829 friction, a purely gravitational effect that can alter a purely Keplerian motion of binaries, and 830 which, under certain conditions, can become as large as gravitational wave emission losses. 831 Strong gravitational lensing is potentially sensitive to the density profile of DM at sufficiently 832 small scales where DM model profiles differ. We have also outlined how the SMBHs observed 833 at the center of active galaxies can be formed from BH seeds from the gravitational collapse of 834 dense cores of fermionic DM. See Section 6 for details.

835
In summary, we express that only theoretical models that join microphysics and macro-836 physics, such as the RAR model, can lead to a comprehensive set of predictions ranging from 837 participle physics to galactic stellar dynamics and to cosmology, which can be put to direct 838 observational scrutiny. We hope that future observations and a further refined analysis of DM 839 models, including the quantum nature of the particles, will lead to a breakthrough in revealing 840 the DM nature. In general relativity, the Lagrangian for a free particle in a gravitational field can be 847 expressed as

877
Due to the spherical symmetry, the metric is invariant under rotations of the polar coor-878 dinate. Therefore, we can assume without loss of generality θ = π/2. In this case, our set of 879 EOM is the following where E and L are the conserved energy and the angular momentum of particle per unit mass. 881 From the condition for mass particle geodesics, g µνẋ µẋν = 1, and the equation of motion 882 for t(τ) and φ(τ), we can obtain the effective potential. For massive particles must be fulfilled 883 that 884 g µνẋ µẋν = g 00 (r)ṫ 2 − g 11 (r)ṙ 2 − r 2θ2 − r 2 sin 2 θφ 2 = 1, so for θ = π 2 and replacing Eqs. (A10) and (A12) in Eq. (A13) we obtain 886 g 00 (r) g 11 (r)ṙ 2 = E 2 − g 00 (r) 1 + L r 2 (A14) 887 where the second term of the right-hand side of the previous equation is the well-known 888 effective potential

889
U 2 e f f (r) ≡ g 00 (r) 1 + L 2 r 2 . (A15) Appendix B. Projection of orbit onto the plane of sky 890 When a telescope measures the motion of a star, it does not measure the real dynamics but 891 rather an apparent one, i.e., it measures the orbit and velocity data projected on the plane that 892 lies perpendicular to the line of sight of the star. For this reason, to compare the theoretical orbit 893 with the observational data, we must project the real orbit on the observation plane in the sky as 894 shown in Fig. A1. This plane-of-sky is described in coordinates (X, Y) defined by the observed 895 angular positions (the declination δ and the right ascension α), where X = R (δ − δ SgrA * ) and 896 Y = R (α − α SgrA * ) being R the distance to the Galactic center [59,114,115]. According to the 897 above description, the apparent theoretical orbit can then be obtained from: 898 X = x cos(X, x) + y cos(X, y), (A16) where the coordinates (x, y) are determined in the plane of the orbit from x = r cos φ and 899 y = r sin φ, while the direction cosines can be found by Eulerian rotation of axes, this is: cos(X, y) = − cos Ω sin ω − sin Ω cos ω cos i, cos(Y, x) = sin Ω cos ω + cos Ω sin ω cos i, cos(Y, y) = − sin Ω sin ω + cos Ω cos ω cos i, finally we find that orbit on plane-of-sky is given by: 902 X = r[cos(φ + ω) cos Ω − sin(φ + ω) sin Ω cos i], (A22) Y = r[cos(φ + ω) sin Ω + sin(φ + ω) cos Ω cos i].
(A23) Figure A1. Projection of the real orbit onto the plane of the sky. The axes originate at Sgr A* (the focus of the ellipse). The picture illustrates the orbital parameters: φ is the azimuth angle of the spherical system of coordinates associated with the x, y, z Cartesian coordinates, i.e., for an elliptic motion in the x-y plane, it is the true anomaly, i is the angle of inclination between the real orbit and the observation plane, Ω is the angle of the ascending node, and ω is the argument of pericenter. It is worth noting that the Z-axis of the coordinate system is defined by the vector pointing from the Solar System to the Galactic center. Reproduced from [34] with permission from Astronomy & Astrophysics, Copyright ESO.

903
Similarly to orbit, the star's radial velocity must also be projected onto the plane of the 904 sky. The radial velocity on the observation plane is defined as the velocity in the observer's 905 direction along the line of sight. Therefore, if we adopt the X-Y plane as the observation plane, 906 then the line of sight is in the direction of the Z-axis given by 907 Z = x cos(Z, x) + y cos(Z, y), where 908 cos(Z, x) = sin ω sin i, cos(Z, y) = cos ω sin i, in terms of φ, ω and i; Z is given by The radial velocity is defined as the velocity along the observer's line of sight, and since Z 911 is the direction along the observer's line of sight, the derivative of Z relative to time gives us 912 the apparent radial velocity of the star. This is 913 dZ dt = u Z = [rφ cos(φ + ω) +ṙ sin(φ + ω)] sin(i).
It is worth noting that in the coordinate system (X, Y, Z), the direction of Z-axis is defined 915 by the vector pointing from the Solar System to the Galactic center and that the axes originate 916 at Sgr A*, which is considered the focus of the orbit [46,52,59].