© The Authors 2024

1. Introduction

The gravitational interaction between a galaxy and a satellite (a dwarf galaxy or a globular cluster) modifies both systems: The host galaxy strips stars from the satellite at a rate that depends on their density profiles and on the orbit, and the density profile of the host shows reorganized matter (due to dynamical friction) in the vicinity of the satellite orbit. The ensemble of stars that are tidally stripped from the satellite constitute the so-called stellar stream (also known as tidal stream). Such streams have been detected in the Milky Way (MW), in Andromeda, and in the Local Volume (Martínez-Delgado et al. 2010).

Stellar streams currently constitute one of the main MW observables that are related to the dynamics, together with other baryonic observables such as the Galaxy rotation curve, the radial surface density profile of the disk, and the vertical density profile of the disk at the solar radius. Tidal streams probe the acceleration field produced by the Galaxy (Johnston et al. 1999; Zhao et al. 1999; Law et al. 2009; Lux et al. 2013; Johnston & Carlberg 2016; Ibata et al. 2016, 2017; Thomas et al. 2017; Reino et al. 2021) as well as its formation history, according to Helmi et al. (1999), Helmi (2020), Ramos et al. (2022), and Cunningham et al. (2023), among others. A stellar stream could be the closest realization of a galactic orbit that can be observed in nature. Nevertheless, the larger the progenitor, the greater the discrepancy between its orbit and the stream phase-space configuration. The taxonomy of streams is very rich (Amorisco 2015) because different gravitational configurations occur, from almost 1D streams whose progenitor is a small globular cluster, to wide shells produced by the partially radial sinking of a large progenitor into the MW center. In addition, an accreted stream whose progenitor is a globular cluster orbiting a MW satellite galaxy is theoretically possible, giving rise to a stellar cocoon around the stream track (Carlberg 2018; Malhan et al. 2019, 2021; Gialluca et al. 2021; Qian et al. 2022). Moreover, perturbations of the streams due to dark matter subhalos are theoretically possible as well. They form off-track features such as the detected spur (Price-Whelan & Bonaca 2018) in the GD-1 stream (Grillmair & Dionatos 2006).

Stream data together with other measurements from baryonic tracers can substantiate claims about an unknown aspect of the gravitational field: Whether it needs to be modeled with a dark matter (DM) halo component or with modified Newtonian dynamics (MoND) theory. For example, within the DM paradigm, spatial densities might have different profiles (e.g., spherical, axisymmetric, or triaxial) depending on the number of conserved components of the angular momentum, which also influences the stream properties (Vera-Ciro & Helmi 2013; Price-Whelan et al. 2016; Mestre et al. 2020). In particular, Malhan & Ibata (2019) constrained the MW dark matter halo shape using Gaia DR2 data of the GD-1 stream assuming an axisymmetric generalization of the Navarro-Frenk-White profile (NFW, Navarro et al. 1996) and obtain a flattening halo parameter of $q=0.{82}_{-0.13}^{+0.25}$, which is compatible with spherical symmetry. Moreover, Palau & Miralda-Escudé (2023) measured the oblateness in an axisymmetric generalization of the NFW profile using three stellar streams: NGC 3201, M68, and Palomar 5. They obtain consistency with a spherical halo.

In addition to the above traditional halo models, which arise from ΛCDM cosmologies or from various ad hoc symmetry considerations, some DM profiles take the quantum nature of the DM candidate into account (see, e.g., Schive et al. 2014 for bosonic profiles that are composed of axion-like particles, and Ruffini et al. 2015; Chavanis et al. 2015; Argüelles et al. 2021 for fermionic profiles that are typically composed of sterile neutrinos). A relevant aspect of this type of profiles is the source of quantum pressure acting in the innermost regions of the halos: While in the boson case, the profiles develop a highly dense soliton, in the fermion case, the profiles develop a degenerate compact core surrounded by a more diluted halo that is self-bounded in radius. In this work, we focus on the latter. A dense core - diluted halo fermionic DM profile like this was obtained from the fully relativistic Ruffini-Argüelles-Rueda (RAR) model, which was successfully applied to different galaxy types (Argüelles et al. 2019; Krut et al. 2023), including the MW from Sagittarius A* until the entire halo (Argüelles et al. 2018; Becerra-Vergara et al. 2020, 2021; Argüelles et al. 2022). This type of model built in terms of a Fermi-Dirac-like phase-space distribution function (including central degeneracy and cutoff in the particle energy) is also known as the relativistic fermionic King (RFK) model (Chavanis 2022b).

We model the 5D track of the GD-1 stellar stream inside a MW with a fermionic dark matter core-halo distribution. At the same time, we aim to explain the dynamics of the S-cluster stars through the high-density fermion core and without assuming a central BH. Finally, we compare the virial mass of the Galaxy as predicted by the fermionic model with that obtained from recent Gaia DR3 rotation-curve data. In Sect. 2 we explain the methods; in Sect. 3 we present the best-fit results, and in Sect. 4 we conclude.

2. Methods

In this section, we explain the observables and the methods we used in this research. The exact pipeline applied in order to obtain the results and plots of this paper can be found at the corresponding GitHub repository¹.

2.1. Observables and assumed measurements

The main observables that we used were computed with the polynomial fits found by Ibata et al. (2020) for the GD-1 stream using astrometric (Gaia DR2) and high-precision spectroscopic datasets, together with the analysis of the STREAMFINDER algorithm. The polynomials are the following:

$$\begin{array}{cc}\hfill {\varphi }_2& =0.008367{\varphi }_1^3-0.05332{\varphi }_1^2-0.07739{\varphi }_1-0.02007,\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill {\stackrel{\sim }{\mu }}_{\alpha }& =3.794{\varphi }_1^3+9.467{\varphi }_1^2+1.615{\varphi }_1-7.844,\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill {\mu }_{\delta }& =-1.225{\varphi }_1^3+8.313{\varphi }_1^2+18.68{\varphi }_1-3.95,\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill v_h& =90.68{\varphi }_1^3+204.5{\varphi }_1^2-254.2{\varphi }_1-261.5,\hfill \end{array}$$(4)

with ϕ₁ and ϕ₂ in radians, ${\stackrel{\sim }{\mu }}_{\alpha }={\mu }_{\alpha }cos\delta$ and μ_δ in mas yr⁻¹, and v_h in km s⁻¹. These quantities correspond to the longitude and latitude in the GD-1 celestial frame of reference (Koposov et al. 2010), the proper motion in right ascension and declination, and the heliocentric radial velocity, respectively. The domain of the polynomials was limited to −90° < ϕ₁ < 10°. To obtain our observable data set, we sampled the domain with 100 equidistant points (${\varphi }_1^{(i)},i=1,\dots ,100$) and evaluated the polynomials at these points, thus obtaining the sets ${\varphi }_1^{(i)}$, ${\varphi }_2^{(i)}$, ${\stackrel{\sim }{\mu }}_{\alpha }^{(i)}$, ${\mu }_{\delta }^{(i)}$, and $v_h^{(i)}$. We did not include among the observables in Eqs. (1–4) the photometric heliocentric distance D, which is justified by the posterior analysis performed in Sect. 3.1. In one of the experiments, we considered an observable of a different nature, the core mass (M_core), which is defined as the mass enclosed at the radius when the circular velocity reaches its first maximum. The constraint for the mass of the core of the distribution was assumed to be M_core = 3.5 × 10⁶ M_⊙, in agreement with Becerra-Vergara et al. (2020, 2021), and Argüelles et al. (2022). The core radius does not include the entire mass inside the innermost S-star pericenter we considered (S2) because the first maximum of the circular velocity corresponds to a shorter distance at which the core density-region is still falling. In Sect. 3 we obtain a DM mass inside the S2 star pericenter of M(r_{peri − S2}) = 4.03 × 10⁶ M_⊙, which agrees very well with the Schwarzschild BH mass constraints of 4.1 × 10⁶ M_⊙ and 3.97 × 10⁶ M_⊙ obtained from the same S2 star in GRAVITY Collaboration (2018) and Do et al. (2019), respectively. The parameters assumed in this paper are the galactocentric distance of the Sun R_⊙ = 8.122 kpc (GRAVITY Collaboration 2018) and the peculiar solar velocity² v_⊙p = (11.1 km s⁻¹,12.24 km s⁻¹,7.25 km s⁻¹) (Schönrich et al. 2010).

2.2. Fermionic dark matter halo model

Our fermionic DM model is a spherical and isotropic distribution of fermions at finite temperature in hydrostatic equilibrium, subject to the laws of General Relativity (GR), that is, the T.O.V. equation complemented with the Tolman and Klein thermodynamics equilibrium conditions and the particle energy conservation along a geodesic, as defined in Argüelles et al. (2018), while using a notation from Chavanis (2020) as detailed below.

We started with a spherically symmetric metric, which is defined as

$$\begin{array}{c}\hfill \mathrm{d}{s}^2=g_{00}(r)\mathrm{d}{t}^2+g_{11}(r)\mathrm{d}{r}^2-r^2\mathrm{d}{\vartheta }^2-r^{2}sin\vartheta \mathrm{d}{\phi }^2,\end{array}$$(5)

with g₀₀(r) = e^ν(r)c² and g₁₁(r) = −e^λ(r), where c is the speed of light, t stands for the time, (r, ϑ, φ) are spherical coordinates, and ν and λ are metric exponents whose radial dependence is computed below.

Our first differential equation is that of the mass versus radius for a spherical system of density ρ,

$$\begin{array}{c}\hfill \frac{\mathrm{d}M}{\mathrm{d}r}=4\pi r^2\rho (r),\end{array}$$(6)

from which we obtained the enclosed mass M(r) at a given radius r by simple integration,

$$\begin{array}{c}\hfill M(r)=4\pi {\displaystyle {\int }_0^r}{r}^{\prime 2}\rho (r^{\prime })\mathrm{d}{r}^{\prime }.\end{array}$$(7)

From the Einstein equations, the relation between M(r) and the metric exponent λ can be found,

$$\begin{array}{c}\hfill {\mathrm{e}}^{-\lambda (r)}=1-\frac{2G}{c^2}\frac{M(r)}{r},\end{array}$$(8)

where G is the gravitational constant.

Additionally, the following version of the T.O.V. equation for the metric exponent ν can be deduced:

$$\begin{array}{c}\hfill \frac{\mathrm{d}\nu }{\mathrm{d}r}=\frac{1+\frac{{\displaystyle 4\pi r^{3}P(r)}}{{\displaystyle M(r)c^2}}}{r(\frac{{\displaystyle r{c}^2}}{{\displaystyle 2GM(r)}}-1)},\end{array}$$(9)

where P is the pressure. The two quantities ρ and P are defined as the following integrals over momentum space,

$$\begin{array}{cc}\hfill \rho (r)& ={\displaystyle \int \frac{E(p)}{c^2}}f(r,p)\mathrm{d}\mathit{p},\hfill \end{array}$$(10)

$$\begin{array}{cc}\hfill P(r)& =\frac{1}{3}{\displaystyle \int p}\frac{\mathrm{d}E(p)}{\mathrm{d}p}f(r,p)\mathrm{d}\mathit{p}=\frac{1}{3}{\displaystyle \int \frac{p^2{c}^2}{E(p)}}f(r,p)\mathrm{d}\mathit{p},\hfill \end{array}$$(11)

where p is the spatial momentum vector, p is its norm, $E(p)=\sqrt{p^2{c}^2+m^2{c}^4}$ is the total energy of a particle of mass m, and f is the phase-space distribution of the system, given by a Fermi-Dirac distribution with an energy cutoff. This number density f can be obtained from a maximum entropy principle computed from a kinetic theory that includes self-gravity and violent relaxation, as shown in Chavanis 2004 (for a review see also Chavanis 2022a), and was recently applied to a vast sample of disk galaxies in Krut et al. (2023). It can be expressed as

$$\begin{array}{c}\hfill f(r,p)=\frac{g}{h^3}\frac{1-{\mathrm{e}}^{[E(p)-E_{\mathrm{c}}(r)]/kT(r)}}{1+{\mathrm{e}}^{[E(p)-\mu (r)]/kT(r)}}\mathrm{if}E(p)\le E_{\mathrm{c}}(r),\\ \hfill \end{array}$$(12)

and f(r, p) = 0 otherwise, where k is the Boltzmann constant, h is the Planck constant, g = 2s + 1 is the spin multiplicity of quantum states, with s = 1/2, and the following local quantities were used: the chemical potential μ(r), the cutoff energy E_c(r), and the effective temperature T(r). The coefficient g/h³ is the maximum accessible value of the distribution function fixed by the Pauli exclusion principle.

The above equations were complemented with two thermodynamic equilibrium conditions given by Tolman (1930) and Klein (1949), together with the condition of energy conservation along the geodesic given in Merafina & Ruffini (1989),

$$\begin{array}{c}\hfill \frac{1}{T}\frac{\mathrm{d}T}{\mathrm{d}r}=\frac{1}{\mu }\frac{\mathrm{d}\mu }{\mathrm{d}r}=\frac{1}{E_{\mathrm{c}}}\frac{\mathrm{d}{E}_{\mathrm{c}}}{\mathrm{d}r}=-\frac{1}{2}\frac{\mathrm{d}\nu }{\mathrm{d}r}.\end{array}$$(13)

Thus, we built a system of five differential equations given by Eqs. (6), (9), and (13), with initial conditions at the center of the distribution M(0) = 0, ν(0) = 0, T(0) = T₀, μ(0) = μ₀, and E_c(0) = E_c0. The differential equations do not depend on ν, but on its radial derivative, so that the system can be solved starting with an arbitrary initial value ν(0) = 0, adding a finite value later, namely ν₀, in such a way that the solution satisfies the condition of continuity with the Schwarzschild metric at the border of the fermionic distribution. In Appendix A we solve the system of equations numerically. Following Argüelles et al. (2018), we use adimensional versions of the initial conditions throughout,

$$\begin{array}{cc}\hfill {\theta }_0& =\frac{{\mu }_0-m{c}^2}{k{T}_0},\hfill \\ \hfill W_0& =\frac{E_{c0}-m{c}^2}{k{T}_0},\hfill \\ \hfill {\beta }_0& =\frac{k{T}_0}{m{c}^2},\hfill \end{array}$$(14)

subtracting the rest-energy in the first two cases.

Although our fermionic system is univocally determined by the four free parameters m, θ₀, W₀, and β₀, in some convenient situations, we use ω₀ = W₀ − θ₀ instead of W₀.

2.3. Milky Way and stream models

We modeled our Galaxy by combining the fermionic dark halo described above, whose parameters are determined in this work, with a fixed baryonic component identical to that in Model I of Pouliasis et al. (2017). We call this full Galaxy model Fermionic-MW.

In addition to this model and for qualitative comparative purposes, we used the Galactic model fit by Malhan & Ibata (2019), which is the MWPotential2014 model with an axisymmetric NFW profile. This model was evaluated by means of the Galpy package (Bovy 2015). A robust statistical comparison between the RAR and other models will be performed in a future paper, but here, we intend to verify that our GD-1 fit agrees with the latest fit performed in the literature. This fitted model has a circular velocity at the position of the Sun v_c(R_⊙) = 244 ± 4 km s⁻¹ and a z-flattening of the DM density distribution $q_{\rho }=0.{82}_{-0.13}^{+0.25}$. We call this Galaxy model NFW-MW.

GD-1 is a dynamically cold stream, and its stars are highly correlated, although its progenitor has not been detected with certainty (de Boer et al. 2018; Price-Whelan & Bonaca 2018; Malhan et al. 2018). Its almost one-dimensional distribution in phase space could be well approximated with the orbit of a theoretical progenitor, however, as previously reported by Malhan & Ibata (2019), Price-Whelan & Bonaca (2018) and Koposov et al. (2010).

The initial conditions of the orbit were given in the spherical equatorial coordinates of the ICRS frame: right ascension α, declination δ, D, ${\stackrel{\sim }{\mu }}_{\alpha }$, μ_δ, and v_h. The code uses the Astropy ecosystem (Astropy Collaboration 2022, 2018, 2013) in order to transform these initial conditions to galactocentric coordinates assuming a galactocentric reference frame with the Sun at the position x_⊙ = (−R_⊙,0,0) and the solar velocity given by the sum of the circular velocity at the position of the Sun and the peculiar solar velocity: v_⊙ = v_⊙p + (0,v_c(R_⊙),0). The circular velocity depends on the model and position and is given by

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill v_{\mathrm{c}}^2(R_{\odot })& =R_{\odot }{|\mathrm{\nabla }\mathrm{\Phi }(\mathit{x})|}_{\mathit{x}={\mathit{x}}_{\odot }}\hfill \\ \hfill & =R_{\odot }{|\mathrm{\nabla }{\mathrm{\Phi }}_{\mathrm{B}}(\mathit{x})|}_{\mathit{x}={\mathit{x}}_{\odot }}+G\frac{M_{\mathrm{DM}}(R_{\odot })}{R_{\odot }},\hfill \end{array}\end{array}$$(15)

where Φ is the total potential, Φ_B is the potential generated by the three baryonic components, and M_DM is the enclosed fermionic DM mass. In the last term, we used the spherical symmetry of the DM distribution in order to relate the acceleration with the enclosed mass. For the NFW-MW model, the circular velocity was computed using the gradient of the total potential, that is, the first line of Eq. (15), performed with the Galpy code.

We integrated the orbit forward and backward in time during a time interval of Δt = 0.2 Gyr, starting in both cases from a given initial condition for the progenitor. In the next sections we explain how these initial conditions were chosen for some simulations and fitted for others. The integrator we used was an eighth-order Runge-Kutta (DOP853 called from the SciPy solve_ivp function) with relative and absolute tolerance parameters given by rtol = 5 × 10⁻¹⁴ and atol = 0.5 × 10⁻¹⁴, respectively.

The resulting orbit was successively transformed into the ICRS and then into the GD-1 coordinate frames. For the latter, we used the GD1Koposov10 class defined in the Gala package (Price-Whelan 2017; Price-Whelan et al. 2020), which uses the transformation matrix defined by Koposov et al. (2010). After these transformations, we obtained the orbit expressed in the observable variables ϕ₁, ϕ₂, ${\stackrel{\sim }{\mu }}_{\alpha }$, μ_δ, and v_h. Finally, we computed these variables at the points ${\varphi }_1^{(i)}$, by interpolation. These values, together with the observed GD-1 data defined in Sect. 2.1, were used to evaluate the following stream function:

$$\begin{array}{cc}\hfill {\chi }_{\mathrm{stream}}^2& ={\chi }_{{\varphi }_2}^2+{\chi }_{{\stackrel{\sim }{\mu }}_{\alpha }}^2+{\chi }_{{\mu }_{\delta }}^2+{\chi }_{v_h}^2\hfill \end{array}$$(16)

$$\begin{array}{cc}\hfill {\chi }_{\eta }^2& =\frac{1}{{\sigma }_{\eta }^2}{\displaystyle \sum _{i=1}^{100}}{\left({\eta }^{(i)}-\eta ({\varphi }_1^{(i)})\right)}^2,\hfill \end{array}$$(17)

where $\eta \in \{{\varphi }_2,v_h,{\stackrel{\sim }{\mu }}_{\alpha },{\mu }_{\delta }\}$, and σ_η are the corresponding dispersions of the stream data points estimated by inspection of Figs. 1, 3, and 4 of Ibata et al. (2020), ${\sigma }_{{\varphi }_2}=0\stackrel{°}{.}5$, σ_vh = 10 km s⁻¹, and ${\sigma }_{{\stackrel{\sim }{\mu }}_{\alpha }}={\sigma }_{{\mu }_{\delta }}=2$ mas yr⁻¹. Thus, ${\chi }_{\text{stream}}^2$ measures the departure of the model from the observed stream.

For some fits, we also considered the departure of the model from a dark mass constraint in the core of the distribution,

$$\begin{array}{c}\hfill {\chi }_{\mathrm{core}}^2=\frac{{(m_{\mathrm{c}}-M_{\mathrm{core}})}^2}{{\sigma }_m^2},\end{array}$$(18)

where the value of M_core was defined in Sect. 2.1, m_c is the core mass of the model (i.e., the variable), and σ_m was fixed at 0.01 M_core. For these fits, we then used the following compound function:

$$\begin{array}{c}\hfill {\chi }_{\mathrm{full}}^2={\chi }_{\mathrm{stream}}^2+{\chi }_{\mathrm{core}}^2.\end{array}$$(19)

We note that in order to compute the estimated core mass for each model, m_c in Eq. (18), we calculated the first local maximum of the circular velocity in GR. This is compactly expressed as

$$\begin{array}{c}\hfill V_{\mathrm{circ},\mathrm{DM}}(r)=c{\left(\frac{r}{2}\frac{g_{00}^{\prime }(r)}{g_{00}(r)}\right)}^{1/2}=c{\left(\frac{r}{2}\frac{\mathrm{d}\nu (r)}{\mathrm{d}r}\right)}^{1/2},\end{array}$$(20)

where the velocity has components Vⁱ = dxⁱ/dt (i = 1, 2, 3) measured in a local frame that is fixed in space at a distance r from the Galaxy center. From the T.O.V. Eq. (9), it is possible to show that far from the core, this relativistic formula for the velocity tends to the classical law, that is, $\sqrt{GM(r)/r}$.

2.4. Optimization algorithms

Our goal was to fit our Fermionic-MW model by minimizing the full χ² function given by Eq. (19) (when fitting the NFW-MW model, we only used the function given by Eq. (16)). To this end, we used two optimization algorithms that belong to the family of black box algorithms, which perform very well when the function to be minimized presents many relative minimums or when the landscape is complex, as in our case. One algorithm is an implementation of a differential evolution algorithm, which is a metaheuristic algorithm that finds the solution of an optimization problem by iteratively updating generations of candidate solutions until a certain tolerance is met. Generally, a few best candidates from each generation survive in order to create the descendants, that is, the next generation, by making stochastic combinations from them. We used the SciPy implementation of this algorithm, called optimize.differential_evolution algorithm, with the metaparemeters given by strategy=“best2bin”, maxiter = 200, popsize = 200, tol = 5 × 10⁻⁸ and atol = 0, unless otherwise stated. This method can be run in parallel with shared memory.

The other algorithm is an implementation of the mesh adaptive direct search (MADS) algorithm called NOMAD (Audet et al. 2021). Audet & Dennis (2006) explained the method in detail. The JULIA (Bezanson et al. 2017) wrapper of this algorithm, https://bbopt.github.io/NOMAD.jl/stable/, was used through the package https://docs.sciml.ai/Optimization/stable/. We used default values of all the metaparameters, except for maxiters = 700.

3. Results

3.1. Fitting the Fermionic-MW model

We fit both the Fermionic-MW model parameters and the initial conditions (IC) of the orbit of the progenitor through four steps that consisted of (i) obtaining an order zero value of the orbit IC using the NFW-MW model, (ii) fitting the Fermionic-MW parameters with fixed IC values, (iii) polishing the IC values with fixed Fermionic-MW potential, and (iv) polishing the Fermionic-MW parametes with fixed IC values.

As a first step, we searched for a provisional but good enough set of IC that was able to reproduce the orbit in our Fermionic-MW model. To this end, we fit the initial conditions in the fixed NFW-MW model by using the ${\chi }_{\text{stream}}^2$ function as defined in Eq. (16) and the differential evolution algorithm specified in Sect. 2.4. Since the optimization algorithm needs bounds for the variables, we used boxes centered near the midpoint of the observable data, that is, η⁽⁵¹⁾, which corresponds to $\alpha =149\stackrel{°}{.}24$, $\delta =36\stackrel{°}{.}61$, ${\stackrel{\sim }{\mu }}_{\alpha }=-5.70$ mas yr⁻¹, μ_δ = −12.48 mas yr⁻¹, and v_h = −18.81 km s⁻¹ plus D = 7.69 kpc³, and with sides of length equal to the absolute values of the variables at their centers (except for α, where we used a side of length α/5).

The differential evolution algorithm converged to the solution $\alpha =149\stackrel{°}{.}25$, $\delta =36\stackrel{°}{.}59$, D = 8.01 kpc, ${\stackrel{\sim }{\mu }}_{\alpha }=-5.56$ mas yr⁻¹, μ_δ = −12.40 mas yr⁻¹, and v_h = −19.15 km s⁻¹, with a value of ${\chi }_{\text{stream}}^2$ = 11.58. The orbit that corresponds to these IC is displayed in the observable space in Fig. 1 (top: ϕ₂; middle: ${\stackrel{\sim }{\mu }}_{\alpha }$, μ_δ; bottom: v_h) with a dotted (green) line. The solid (black) line shows the corresponding observable data η, while the shaded (gray) area demarcates the corresponding 1σ_η regions.

Fig. 1. Stream fits in observable space: Sky position (top: ϕ2), proper motions (middle: ${\stackrel{\sim }{\mu }}_{\alpha }$, μδ), and heliocentric velocity (bottom: vh).

Previous works (Argüelles et al. 2018, 2019; Krut et al. 2023) have shown that the parameters of the fermionic DM model can be split into two sets for the family of fermionic DM profiles with highly degenerate cores (i.e., θ₀ ≳ 15 Argüelles et al. 2019). On one hand, m and β₀ control the core of the distribution in the sense that for given values of M_core and m, it is possible to find a consistent value of β₀ with little influence on the other two parameters θ₀ and W₀. This implies a partial degeneracy between m and β₀. On the other hand, for the same type of fermionic core-halo solutions with a positive central degeneracy as mentioned above, the main behavior of the distribution in the halo is determined by ω₀, as reported in Argüelles et al. (2019) for different galaxy types.

As a starting point for the fitting of the Fermionic-MW parameters (second step), we took the values obtained by Becerra-Vergara et al. (2020): (m, β₀, θ₀, W₀) = (56 keV c⁻², 1.1977  ×  10⁻⁵, 37.765, 66.3407), which allowed their MW model (DM plus baryons) to satisfy the geodesic motions of both S2 and G2 at Sagittarius A*, and the rotation curve from Sofue (2013). From these values, fixing m = 56 keV c⁻² and taking the already fitted IC of the stream progenitor for the NFW-MW model, we performed a differential evolution minimization of the ${\chi }_{\text{full}}^2$ function in the window (θ₀, ω₀, β₀)∈[35, 40]×[25, 30]×[10⁻⁵, 1.5 × 10⁻⁵]. Using metaparameter values maxiter  =  popsize  =  300, the algorithm converged to (θ₀, ω₀, β₀) = (36.094, 27.368, 1.252 × 10⁻⁵), giving ${\chi }_{\text{stream}}^2$ = 16.190 and ${\chi }_{\text{core}}^2=7.676\times {10}^{-10}$.

The third step consisted of polishing the IC of the orbit by using the differential evolution algorithm with fixed Fermionic-MW parameters, with metaparameter values maxiter  =  popsize  =  400, which gave a result very similar to that of the FMW-MW case: $\alpha =149\stackrel{°}{.}39$, $\delta =36\stackrel{°}{.}87$, D = 8.02 kpc, ${\stackrel{\sim }{\mu }}_{\alpha }=-5.55$ mas yr⁻¹, μ_δ = −12.33 mas yr⁻¹ and v_h = −20.84 km s⁻¹, with an improved value of ${\chi }_{\text{stream}}^2$ = 3.59.

The last step consisted of polishing the fermionic parameters using the second optimization algorithm described in Sect. 2.4, that is, NOMAD. We divided a macroscopic orthohedron in parameter space, (θ₀, ω₀, β₀)∈[35.8, 36.3]×[27.0, 27.6]×[1.2 × 10⁻⁵, 1.3 × 10⁻⁵], into 17³ = 4913 smaller orthohedrons. In each subregion, we performed an independent optimization in a parallel distributed scheme and searched for the parameters that minimized ${\chi }_{\text{full}}^2$. Then, we selected the global minimum by comparing the results of each distributed process and obtained the following final fit parameters of the model: (θ₀, ω₀, β₀) = (36.0704, 27.3501, 1.2527 × 10⁻⁵), giving ${\chi }_{\text{full}}^2$ = 13.53. The corresponding orbit is displayed in the observable space in Fig. 1 with a dashed (amber) line. The Fermionic-MW and NFW-MW models fit the GD-1 stream very well. We did not perform any statistically rigorous comparison of the two models to determine which was more consistent with the data. In a future paper, we will compute the posterior distribution of the fit parameters and will then be able to give error bounds.

For completeness, in Fig. 2 we plot the heliocentric distance (D) using the same line types as in Fig. 1. The solid (black) line corresponds to the fifth-order polynomial fit in Ibata et al. (2020) to the photometric distance measured there,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill D& =-4.302{\varphi }_1^5-11.54{\varphi }_1^4-7.161{\varphi }_1^3+5.985{\varphi }_1^2\hfill \\ \hfill & +8.595{\varphi }_1+10.36.\hfill \end{array}\end{array}$$(21)

Fig. 2. Photometric distance (D). This was not included as an observable in the stream-fitting procedure.

The two theoretical D(ϕ₁) curves agree with each other, but they both differ considerably with respect to the observed curve (i.e., the polynomial). Because the polynomial presents a suspicious constant value for ϕ₁ ∈ [ − 70, −30], we plotted the galactocentric distance versus ϕ₁ (not shown), and we found that this latter curve presents an unphysical wobbling behavior. We therefore did not include the photometric distance as a fitting observable.

It is instructive to see how the value of ${\chi }_{\text{stream}}^2$ is modified when the parameters θ₀ and ω₀ (i.e., those that control the halo) are varied. To this end, we fixed (m, β₀) = (56 keV c⁻², 1.2527 × 10⁻⁵) and varied (θ₀, ω₀) in a grid spanning the rectangle [34, 38]×[25.5, 28.2]. Fig. 3 shows a contour plot of this grid, colored by the value of their corresponding ${\chi }_{\text{stream}}^2$; the black point corresponds to our fit solution. The lowest values of the function are located in a thin valley, which shows that the solution of the fitting problem is locally degenerate along a straight line in the (θ₀, ω₀) plane. To show the behavior of ${\chi }_{\text{stream}}^2$ around the solution, we first fit a line to the points that satisfied ${\chi }_{\text{stream}}^2$ < 50. We obtained h(x) = 0.7939 + 0.7362x. We then subtracted this line from ω₀, and finally remade a plot of ${\chi }_{\text{stream}}^2$ in the window (θ₀, ω₀ − h(θ₀)) ∈ [35, 37]×[−0.011, 0.011]. The result is shown in Fig. 4. It is noticeable that an absolute minimum indeed exists (nondegenerate problem) and that the variance along ω₀ − h(θ₀) is two orders of magnitude smaller than the variance along θ₀.

Fig. 3. Values of ${\chi }_{\text{stream}}^2$ for (m, β0) = (56 keV c−2, 1.254 × 10−5) in the window (θ0, ω0)∈[34, 38]×[25.5, 28.2]. The black point corresponds to our solution, and the gray region corresponds to ${\chi }_{\text{stream}}^2$ > 104. The lowest values of the function are located along a thin valley.

Fig. 4. Values of ${\chi }_{\text{stream}}^2$ for (m, β0) = (56 keV c−2, 1.254 × 10−5) in the window (θ0, ω0 − h(θ0)) ∈ [35, 37]×[−0.011, 0.011]. The black point corresponds to our solution, and the gray region corresponds to ${\chi }_{\text{stream}}^2$ > 500.

3.2. Rotation curves, accelerations, and virial quantities

We computed the resulting rotation curves for our two models and compared them with three observed rotation curves (Fig. 5).

Fig. 5. Rotation curves of the Fermionic-MW model (this work, in dashed amber) and the NFW-MW model (Malhan & Ibata 2019, in dotted green), which fit the GD-1 stream, are compared a posteriori with different observed rotation curves (Eilers et al. 2019 with purple triangles, Sofue 2020 with light blue rhombi, and Jiao et al. 2023 with black circles). Only the Fermionic MW model can account for the GD-1 stream data and the sharp drop of the recent Gaia DR3 rotation curve.

In order to build a unified rotation curve, Sofue (2020) computed the running average of many rotation curves resulting from different dynamical tracers according to the galactocentric distance. In the central parts of the Galaxy, the tracers used were the molecular gas and the infrared stellar motion, while in the outer parts (beyond r ∼ 30 kpc), the tracers were the radial motions of satellite galaxies and globular clusters. Rotation curves resulting from Galactic disk objects were also used. On the other hand, Eilers et al. (2019) used a selection of red giant branch stars as tracers of the disk dynamics.

It is worth noticing that Sofue (2020) assumed (R_⊙, v_c(R_⊙)) = (8 kpc, 238 km s⁻¹), while Eilers et al. (2019) assumed R_⊙ = 8.122 kpc and a solar galactocentric velocity v_⊙ = (11.1, 245.8, 7.8) km s⁻¹, with which they estimated to be v_c(R_⊙) = 229 ± 0.2 km s⁻¹.

It is interesting to note that the two theoretical models, Fermionic-MW and NFW-MW, both give v_c(R_⊙)≈244 km s⁻¹ at the solar radius, which agree excellently with the estimate found by Bovy (2020), v_c(R_⊙) = 244 ± 8 km s⁻¹ for R_⊙ = 8.275 kpc, or v_c(R_⊙) = 242 ± 8 km s⁻¹ for our adopted value of R_⊙. Although this velocity is higher than the standard value (220–230 km s⁻¹), it should be mentioned that according to Table 1 in Sofue (2020), our computed velocity is not an outlier (see also Section 6.2 of Honma et al. 2012). Our solution has an average slope s = −4.15 ± 0.015 km s⁻¹ kpc⁻¹, fit for 14.5 kpc < r < 26.5 kpc, which is comparable to the corresponding value of s = −3.93 ± 0.15 km s⁻¹ kpc⁻¹ measured from the rotation curve of Jiao et al. (2023), and agrees better than the corresponding slope of the NFW-MW model (see also Fig. 5 for comparison).

We computed the galactocentric distance of the GD-1 observables projected onto the plane z = 0 and found that it lies inside the interval 11.5 kpc  ≲  r  ≲  16.4 kpc, which is displayed as a vertical shaded (gray) band in the Galaxy RC of Fig. 5. The location of the GD-1 stream orbit corresponds to z ∈ [2.6, 9.7] kpc and thus explores the nonsphericality of the full MW models (due to the axisymmetry of baryons and NFW). It is noticeable that both models approximately agree in their circular velocity in the GD-1 region. The galactocentric distance (not projected) of the stream corresponds to the interval 13.9 kpc ≲ r ≲ 16.6 kpc (subject to errors in the photometric distance, as commented in Section 3.1). The theoretical orbit of the stream in the Fermionic-MW has a pericenter of 14.3 kpc and an apocenter of 24.5 kpc, and it is currently in its pericentric passage.

We computed the cylindrical components of the acceleration, a_r and a_z, for the two MW models along the stream and obtained a maximum difference of |Δa_r| ≲ 0.08 km s⁻¹ Myr⁻¹ and |Δa_{r| ≲ 0.15 km s⁻¹ Myr⁻¹}, respectively. This agreement supports the idea that cold tidal streams are excellent accelerometers (Ibata et al. 2016; Nibauer et al. 2022; Craig et al. 2023).

With respect to virial quantities, the core-halo dark matter solution has a finite virial radius r_DM, vir = 27.4 kpc and a virial mass M_DM, vir = 1.4 × 10¹¹ M_⊙. The total baryon mass of our model is M_b = 0.9 × 10¹¹ M_⊙, and the total virial mass therefore amounts to M_vir = 2.3 × 10¹¹ M_⊙. The value of the MW total mass at 50 kpc reported in Table 3 of Gibbons et al. (2014) is 2.9 × 10¹¹ M_⊙ with (σ, 2σ) = (0.4, 0.9)×10¹¹ M_⊙, which means that our solution lies in the 2σ region. The mass of the fermionic solution is constant for radii larger than r_DM, vir, while the mass of the model studied by Gibbons et al. (2014) continues to increase with radius according to their Table 3, although their estimates at large radii have relatively high error bounds, for instance, 2σ = 3 × 10¹¹ M_⊙ for M(200 kpc) = 5.6 × 10¹¹ M_⊙.

The most recent MW mass estimates were obtained by Jiao et al. (2023) and Ou et al. (2023) from Gaia DR3 data. These authors reported data that are compatible with even lower values of the MW virial mass of ≈2 × 10¹¹ M_⊙ (however, see Koop et al. 2024), in agreement with our fermionic model predictions. These MW mass estimates correspond to a sharp Keplerian decline in the MW rotation that ends at ≈26.5 kpc. This again agrees remarkably well with the virial radius predicted by our Fermionic-MW model of ≈27 kpc.

The fitted fermionic DM model has a density in the solar neighborhood of ρ_DM, ⊙ = 1.46 × 10⁷ M_⊙ kpc⁻³ = 0.55 GeV cm⁻³c⁻², which falls within the 2σ region of the estimate made by Salucci et al. 2010 (0.43 ± 0.21 GeV cm⁻³c⁻²), but is higher than the one obtained by Eilers et al. 2019 (0.30 ± 0.03 GeV cm⁻³c⁻²) or by Ou et al. 2024 (0.447 ± 0.004 GeV cm⁻³c⁻²).

3.3. Example of the S-cluster fit: Paradigmatic case of the S2 orbit

In this section, we answer the relevant question of how well the fermionic model that fits the stream according to the procedure of Sect. 3.1 can fit the iconic S2 star orbit with a focus in Sagittarius A*. Even if a good fit is expected since the core mass of the DM distribution M_core = 3.5 × 10⁶ M_⊙ was fit together with the stream constraint, it is important to note that the free parameters of the Fermionic-MW model are not precisely the same as those obtained in Becerra-Vergara et al. (2020). That is, while in Becerra-Vergara et al. (2020) the free DM model parameters were obtained to explain the S2 star geodesic together with the MW RC as given in Sofue (2013), here, the DM halo region was instead fit in order to reproduce the GD-1 stream, with somewhat different (β₀,θ₀,W₀) values.

We therefore performed a least-squares fit following Becerra-Vergara et al. (2020) for the case of the S2 orbit as an example. As commented above, in this case, we used the Fermionic-MW DM model, which fits the GD-1 stream best, that is, (θ₀, ω₀, β₀) = (36.0704, 27.3501, 1.2527 × 10⁻⁵). We obtained excellent results. In Fig. 6, we show the projected S2 orbit in the plane of the sky, while in Fig. 7, we show the time evolution of the redshift function z (which is directly related to the heliocentric radial velocity according to Equation (C.17a) in Becerra-Vergara et al. 2020), α and δ, for the best-fitting values of the osculating orbital parameters. These values are given in Table 1, along with the model-predicted value of the periapsis precession per orbital period, Δϕ, and the orbital period, P.

Fig. 6. Modeled and observed projected orbit in the sky for a Fermionic-MW DM model, shown as a solid (amber) line, and a BH of Mbh = 4.075 × 106 M⊙, shown as a dotted (green) line.

Fig. 7. Redshift (z), right ascension (α), and declination (δ) as a function of time epoch for the same models as displayed in Fig. 6.

Our fitting procedure was applied in the gravitational field of two different scenarios: a Fermionic-MW DM model for m = 56 keV c⁻² and M_core = 3.5 × 10⁶ M_⊙, and a Schwarzschild BH model with a central mass of M_BH = 4.075 × 10⁶ M_⊙. The resulting values for the ${\chi }_{\text{S2}}^2$ minimization presented in Table 2 agree perfectly with those obtained in Becerra-Vergara et al. (2020)⁴. We used the latest publicly accessible data from Do et al. (2019).

For the exemplified case of the S2 star, the models differ in the predicted value of the relativistic precession of the S2 periapsis Δϕ. This interesting relativistic effect in the case of a regular (i.e., nonsingular) DM core was compared with the publicly available astrometric data of S2 and was compared with the BH case in Argüelles et al. (2022). These authors showed that higher particle masses (i.e., leading to more compact DM cores, as detailed in the section below) imply that less extended DM mass fills the S2 orbit. Thus, precession growths from retrograde to prograde as it tends to the unique value predicted by the BH model. Argüelles et al. (2022) showed that already for particle masses of m = 60 keV c⁻², which is slightly above the value considered here, m = 56 keV c⁻², the periapsis precession is very close the one predicted by the Schwarzschild BH. Conversely, for particle masses m ≲ 56 keV c⁻², the DM core is too extended in radius and produces high values of retrograde S2 periapsis precession and poorer orbit fits (Becerra-Vergara et al. 2020; Argüelles et al. 2022).

3.4. Varying the fermion mass to reach more compact cores

As already mentioned, we found a fermionic solution that agrees with both GD-1 data and the geodesic motion of the best-studied S-cluster star around Sagittarius A*, the S2 star. However, this S-star constraint is also satisfied by any fermionic DM profile with a more compact core than the solution corresponding to m = 56 keV c⁻², as long as the core mass is m_c ≈  M_core. It is therefore interesting to determine how much compactness can be reached while keeping both the S2 star and GD-1 constraints, in the light of the new observations of the Event Horizon Telescope, Event Horizon Telescope Collaboration (2022), where a shadow angular diameter of 48.7 ± 7.0 μas has been measured. This diameter corresponds to a shadow radius of ∼2.46 Schwarzschild radii assuming a BH mass of M_bh = 4.075 × 10⁶ M_⊙. In order to extend the fermionic solutions to other values of the fermion mass (m), we used the second optimization algorithm described in Sect. 2.4 for m = 100, 200, 300, and 360 keV c⁻². For each fermion mass, we divided a given macroscopic orthohedron⁵ in parameter space in 20³ = 8000 smaller orthohedrons. In each subregion, we performed an independent optimization in a parallel distributed scheme to search for the parameters that minimize ${\chi }_{\text{full}}^2$ with the NOMAD algorithm. Afterward, we selected the global minimum by comparing the results of each distributed process. The result is that for all the fermion masses, it is possible to find values of the other parameters in such a way that both the GD-1 stream and the core mass constraints are respected with the same precision as in the initial (m = 56 keV c⁻²) case. As shown in Fig. 8, all the solutions have the same density profile in the halo region, while their difference is limited to the compactness of the core. In Table 2 we show the values of ${\chi }_{\text{stream}}^2$ for all the fermion masses we studied. We also give the value of ${\chi }_{\text{S2}}^2$ for the Fermionic-MW and the BH models. All the fermionic models we analyzed are statistically indistinguishable⁶, and more data (e.g., a central shadow feature or fainter S-stars closer to SgrA* than S2) are needed in order to further constrain the particle mass range. Both projects are currently being developed within our group. The values of the core radii of these solutions are approximately 1097, 232, 35, 10, and 5 Schwarzschild radii for m = 56, 100, 200, 300, and 360 keV c⁻², respectively. The last analyzed value of m = 360 keV c⁻² corresponds to a DM core very close to the last stable solution according to the stability criterion of Argüelles et al. (2021), leading to a gravitational core-collapse into a BH of about 4 × 10⁶ M_⊙.

Fig. 8. Fermionic DM density profiles with different core compacities (i.e., different m) fitting the GD-1 stream and the DM core mass, which agree with the S2 star data orbiting Sagittarius A*. The Schwarzschild radius is computed assuming a BH mass Mbh = 4.075 × 106 M⊙.

4. Conclusions

We have fitted the GD-1 stream located at about 14 − 15 kpc from the Galaxy center and the S2 star orbit located at miliparsec scales in a MW potential consisting of a fermionic core-halo DM distribution (Argüelles et al. 2018; Becerra-Vergara et al. 2020, 2021; Argüelles et al. 2022), plus a fixed baryonic distribution (Pouliasis et al. 2017). Remarkably, the total resulting MW mass and the virial radius of the Galaxy predicted by the fermionic DM model agree very well with the virial mass of ≈2 × 10¹¹ M_⊙ and the galactocentric radial range of ≈20 − 26 kpc in which the MW rotation curve sharply drops, as recently meassured by Gaia DR3. This relatively low-mass MW is also consistent with the independent estimate given by Gibbons et al. (2014), although it is considerably lower than typical values given in the literature, for instance, Watkins et al. (2010) and references therein.

We obtained the free parameters of the fermionic model by fixing the fermion mass and by simultaneously fitting two astrophysical constraints: the stream observables, and a DM core mass of 3.5 × 10⁶ M_⊙, the latter taken from previous fits of the S-stars cluster at the center of the Galaxy without a central BH (Becerra-Vergara et al. 2020, 2021). We could thus reproduce the polynomials fit by Ibata et al. (2020), which correspond to the observed sky position, proper motion, and radial velocity of the stream, with a high accuracy and for fermion masses ranging from 56 to 360 keV c⁻².

In order to compare our results with other GD-1 fits in the literature, we also fit the stream orbit with the axisymmetric generalization of the NFW distribution from Malhan & Ibata (2019). We obtained agreement in the GD-1 phase-space track of both the Fermionic-MW and the NFW-MW models. Additionally, we obtained agreement of the two models in the rotation curves at projected (z = 0) galactocentric distances, r, corresponding to the stream observables, that is, 11.5 ≤ r ≤ 16.4 kpc. The average slope of the rotation curve between 14.5 and 26.5 kpc obtained from the Fermionic-MW model was s = −4.18 ± 0.02 km s⁻¹ kpc⁻¹, which agrees much better with the recent observations of Jiao et al. (2023) than the NFW-MW model. We showed that the two MW models perfectly agree in their acceleration vectors as a function of the position along the stream.

We found a circular velocity at the position of the Sun of v_c(R_⊙) = 244 km s⁻¹, which is in line with the value independently obtained by Malhan & Ibata (2019).

The fermionic DM solution has a finite radius of r_DM, vir = 27.4 kpc and a virial mass of M_DM, vir = 1.4 × 10¹¹ M_⊙, implying a total (DM plus baryons) virial mass of the Galaxy of M_vir = 2.3 × 10¹¹ M_⊙, which is at 2σ from the value reported in Table 3 of Gibbons et al. (2014) for a radius of 50 kpc. The value obtained for the DM density in the solar neighborhood is ρ_DM, ⊙ = 1.46 × 10⁷ M_⊙ kpc⁻³ = 0.55 GeV cm⁻³c⁻², which falls inside the 2σ region of a previous estimate by Salucci et al. (2010).

Finally, we showed that it is possible to find a 1D family of solutions parameterized by the fermion mass with the same halo that fits the GD-1 stream, but with a different compactness of the central core. This always reproduces the S2 star orbit (see Figs. 6 and 7 for the case of m = 56 keV c⁻²). For the limiting case studied here (m = 360 keV c⁻²), we obtained a core radius of r_c ≈ 5 Schwarzschild radii. A precise relativistic ray-tracing study about simulated ring-like images of the central cores of fermionic distributions is in progress. We try to place strict constraints on the minimum compactness needed for an agreement with the EHT observations.

In summary, the findings of this work not only support the idea that cold tidal streams are excellent probes of the acceleration field of the MW, but they show that the (spherical) fermionic model is able to agree with a set of independent observables covering three totally different Galaxy distance scales: ∼10⁻⁶ kpc (S-cluster), ∼14 kpc (GD-1), and ∼30 kpc (Gaia DR3 RC mass estimates).

Acknowledgments

MFM would like to dedicate this work to the memory of a great friend since childhood, Gustavo Fabián Larrion (1980-2021), Q.E.P.D. We thank the anonymous referee for the helpful comments that improved the quality of the paper. MFM thanks Khyati Malhan for his comments about the NFW-MW model. MFM thanks Leandro Martínez, Ian Weaver, Joaquín Pelle and the JULIA community for their great support with the language and workflow. We thank Jorge A. Rueda for his comments about the manuscript. We thank Juan Ignacio Rodriguez for his great support with hardware and software issues regarding the IALP server and personal computers. We thank Federico Bareilles and the informatic support team of the FCAGLP. This work used computational resources from CCAD – Universidad Nacional de Córdoba (https://ccad.unc.edu.ar/), which are part of SNCAD – MinCyT, República Argentina. This work also used computational resources from the HPC center DIRAC, funded by Instituto de Física de Buenos Aires (UBA-CONICET) and part of SNCAD-MinCyT initiative, Argentina. DDC and MFM acknowledge support from CONICET (PIP2169) and from the Universidad Nacional de La Plata (G178). C.R.A. acknowledges support from CONICET, the ANPCyT (grant PICT-2018-03743), and ICRANet. VC thanks financial support from CONICET, Argentina. The figures presented in this work where made with the AlgebraOfGraphics.jl (https://aog.makie.org/dev/), Makie.jl (Danisch & Krumbiegel 2021) and Matplotlib (Hunter 2007) packages. Some of our optimization results were initially guided by the use of the LN_NELDERMEAD (Nelder & Mead 1965; Box 1965; Shere 1974) algorithm, from the NLopt.jl (Johnson 2007) package. In order to run JULIA in a parallel SLURM environment we made use of the Distributed.jl and https://github.com/kleinhenz/SlurmClusterManager.jl packages.

References

Appendix A: Numerical solution of the system of differential equations

In this section, we explain how we numerically solved the physical equations defined in Sec. 2.2. We start by defining some constants,

$$\begin{array}{cc}\hfill {\rho }_{\bullet }& =g{\pi }^{3/2}{m}^4{c}^3{h}^{-3},\hfill \end{array}$$(A.1)

$$\begin{array}{cc}\hfill r_{\bullet }& =c/{(8\pi G{\rho }_{\bullet })}^{1/2},\hfill \end{array}$$(A.2)

and introduce the following changes of variables:

$$\begin{array}{cc}\hfill \zeta (r)& =ln(r/r_{\bullet }),\hfill \end{array}$$(A.3)

$$\begin{array}{cc}\hfill z(\zeta )& =ln\psi (r(\zeta )),\hfill \end{array}$$(A.4)

$$\begin{array}{cc}\hfill \stackrel{\sim }{\nu }(\zeta )& =\nu (r(\zeta )),\hfill \end{array}$$(A.5)

$$\begin{array}{cc}\hfill \beta (\zeta )& =\frac{kT(r(\zeta ))}{m{c}^2},\hfill \end{array}$$(A.6)

$$\begin{array}{cc}\hfill \alpha (\zeta )& =\frac{\mu (r(\zeta ))}{m{c}^2},\hfill \end{array}$$(A.7)

$$\begin{array}{cc}\hfill {ϵ}_c(\zeta )& =\frac{E_c(r(\zeta ))}{m{c}^2},\hfill \end{array}$$(A.8)

$$\begin{array}{cc}\hfill ϵ(p)& =\frac{E(p)}{m{c}^2},\hfill \end{array}$$(A.9)

where

$$\begin{array}{c}\hfill \psi (r)=1-{\text{e}}^{-\lambda (r)}=\frac{2G}{c^2}\frac{M(r)}{r},\end{array}$$(A.10)

into equations (6), (9) and (13), obtaining

$$\begin{array}{cc}\hfill \frac{\mathrm{d}z}{\mathrm{d}\zeta }& =-1+{\mathrm{e}}^{2\zeta -z}\frac{\stackrel{\sim }{\rho }(\zeta )}{{\rho }_{\bullet }},\hfill \end{array}$$(A.11)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}\stackrel{\sim }{\nu }}{\mathrm{d}\zeta }& =({\mathrm{e}}^z+{\mathrm{e}}^{2\zeta }\frac{\stackrel{\sim }{P}(\zeta )}{{\rho }_{\bullet }{c}^2}){(1-{\mathrm{e}}^z)}^{-1},\hfill \end{array}$$(A.12)

$$\begin{array}{cc}\hfill \frac{1}{\beta }\frac{\mathrm{d}\beta }{\mathrm{d}\zeta }& =\frac{1}{\alpha }\frac{\mathrm{d}\alpha }{\mathrm{d}\zeta }=\frac{1}{{ϵ}_{\mathrm{c}}}\frac{\mathrm{d}{ϵ}_{\mathrm{c}}}{\mathrm{d}\zeta }=-\frac{1}{2}\frac{\mathrm{d}\stackrel{\sim }{\nu }}{\mathrm{d}\zeta }.\hfill \end{array}$$(A.13)

The thermodynamical quantities, density and pressure, are given by

$$\begin{array}{cc}\hfill \stackrel{\sim }{\rho }(\zeta )& =\frac{4{\rho }_{\bullet }}{\sqrt{\pi }}{\displaystyle {\int }_1^{\infty }}{ϵ}^2{[{ϵ}^2-1]}^{1/2}\stackrel{\sim }{f}(\zeta ,ϵ)\mathrm{d}ϵ,\hfill \end{array}$$(A.14)

$$\begin{array}{cc}\hfill \stackrel{\sim }{P}(\zeta )& =\frac{4c^2{\rho }_{\bullet }}{3\sqrt{\pi }}{\displaystyle {\int }_1^{\infty }}{[{ϵ}^2-1]}^{3/2}\stackrel{\sim }{f}(\zeta ,ϵ)\mathrm{d}ϵ,\hfill \end{array}$$(A.15)

where the fermionic King distribution as a function of ϵ = E/mc² in the new variables is given by

$$\begin{array}{c}\hfill \stackrel{\sim }{f}(\zeta ,ϵ)\equiv \frac{h^3}{g}f(r(\zeta ),p(ϵ))=\frac{1-{\mathrm{e}}^{[ϵ-{ϵ}_{\mathrm{c}}(\zeta )]/\beta (\zeta )}}{1+{\mathrm{e}}^{[ϵ-\alpha (\zeta )]/\beta (\zeta )}}\mathrm{if}ϵ\le {ϵ}_{\mathrm{c}}(\zeta ),\end{array}$$(A.16)

and $\stackrel{\sim }{f}(\zeta ,ϵ)=0$ otherwise.

Equations (A.13) can be analytically integrated to obtain

$$\begin{array}{cc}\hfill \beta (\zeta )& ={\beta }_0{\mathrm{e}}^{-\stackrel{\sim }{\nu }(\zeta )/2},\hfill \\ \hfill \alpha (\zeta )& ={\alpha }_0{\mathrm{e}}^{-\stackrel{\sim }{\nu }(\zeta )/2},\hfill \\ \hfill {ϵ}_c(\zeta )& ={ϵ}_{c0}{\mathrm{e}}^{-\stackrel{\sim }{\nu }(\zeta )/2},\hfill \end{array}$$(A.17)

thus transforming the original system of five differential equations, that is, (A.11-A.13), into a system of just two differential equations that are to be solved numerically subject to the constraints (A.17).

It is not possible to integrate these equations from r = 0 because ζ(r) diverges at the origin. Therefore, the following approximations for the initial conditions at a value r_min ≳ 0 were used:

$$\begin{array}{cc}\hfill \nu (r_{\min })& =\frac{1}{3}\frac{{\rho }_0}{{\rho }_{\bullet }}{\left[\frac{r_{\min }}{r_{\bullet }}\right]}^2\equiv \tau ,\hfill \end{array}$$(A.18)

$$\begin{array}{cc}\hfill \psi (r_{\min })& =\frac{1}{3}\frac{{\rho }_0}{{\rho }_{\bullet }}{\left[\frac{r_{\min }}{r_{\bullet }}\right]}^2=\tau ,\hfill \end{array}$$(A.19)

which implies

$$\begin{array}{c}\hfill \frac{r_{\min }}{r_{\bullet }}=\sqrt{3\tau \frac{{\rho }_{\bullet }}{{\rho }_0}},\end{array}$$(A.20)

where τ ≡ 2 × 10⁻¹⁵ and

$$\begin{array}{c}\hfill {\rho }_0\equiv \rho (0)=\frac{4{\rho }_{\bullet }}{\sqrt{\pi }}{\displaystyle {\int }_1^{\infty }}{ϵ}^2{[{ϵ}^2-1]}^{1/2}{\stackrel{\sim }{f}}_0(ϵ)\mathrm{d}ϵ,\end{array}$$(A.21)

where

$$\begin{array}{c}\hfill {\stackrel{\sim }{f}}_0(ϵ)=\frac{1-{\mathrm{e}}^{[ϵ-{ϵ}_{\mathrm{c}0}]/{\beta }_0}}{1+{\mathrm{e}}^{[ϵ-{\alpha }_0]/{\beta }_0}}\mathrm{if}ϵ\le {ϵ}_{\mathrm{c}0},\end{array}$$(A.22)

and ${\stackrel{\sim }{f}}_0(ϵ)=0$ otherwise.

In this way, the initial conditions of our numerical system are given by ζ_min = ζ(r_min), ${\stackrel{\sim }{\nu }}_{\min }=\tau$, and z_min = lnτ, and the system parameters to be varied are m, β₀, α₀, and ϵ_c0. In agreement with Eqs. (14), we used as parameters the following normalized quantities: m, β₀, θ₀ = (α₀ − 1)/β₀, and W₀ = (ϵ_c0 − 1)/β₀, or ω₀ = W₀ − θ₀ instead of W₀ in some cases.

Equations (A.11) and (A.12) were solved with the LSODA algorithm using a PYTHON (Van Rossum & Drake 1995) script⁷ that makes use of the NumPy (Harris et al. 2020) and SciPy (Virtanen et al. 2020) libraries. We used relative and absolute tolerance parameters given by rtol = 5 × 10⁻¹⁴ and atol = 0, respectively.

After obtaining the numerical solution, since the right sides of Eqns. (A.11) and (A.12) do not depend on the metric, but only on its radial derivative, we can add a constant ${\stackrel{\sim }{\nu }}_0$ to the solution in order to satisfy the condition of continuity with the Schwarzschild metric at the border of the fermion distribution, obtaining

$$\begin{array}{c}\hfill {\stackrel{\sim }{\nu }}_0=2ln\left(\frac{{\beta }_{\mathrm{b}}}{{\beta }_0}\sqrt{1-{\psi }_{\mathrm{b}}}\right),\end{array}$$(A.23)

where ψ_b and β_b are quantities evaluated at the border, that is, when ρ(r_b)→0.⁸

All Tables

All Figures

	Fig. 1. Stream fits in observable space: Sky position (top: ϕ2), proper motions (middle: ${\stackrel{\sim }{\mu }}_{\alpha }$, μδ), and heliocentric velocity (bottom: vh).
In the text

	Fig. 2. Photometric distance (D). This was not included as an observable in the stream-fitting procedure.
In the text

	Fig. 3. Values of ${\chi }_{\text{stream}}^2$ for (m, β0) = (56 keV c−2, 1.254 × 10−5) in the window (θ0, ω0)∈[34, 38]×[25.5, 28.2]. The black point corresponds to our solution, and the gray region corresponds to ${\chi }_{\text{stream}}^2$ > 104. The lowest values of the function are located along a thin valley.
In the text

	Fig. 4. Values of ${\chi }_{\text{stream}}^2$ for (m, β0) = (56 keV c−2, 1.254 × 10−5) in the window (θ0, ω0 − h(θ0)) ∈ [35, 37]×[−0.011, 0.011]. The black point corresponds to our solution, and the gray region corresponds to ${\chi }_{\text{stream}}^2$ > 500.
In the text

Fig. 5. Rotation curves of the Fermionic-MW model (this work, in dashed amber) and the NFW-MW model (Malhan & Ibata 2019, in dotted green), which fit the GD-1 stream, are compared a posteriori with different observed rotation curves (Eilers et al. 2019 with purple triangles, Sofue 2020 with light blue rhombi, and Jiao et al. 2023 with black circles). Only the Fermionic MW model can account for the GD-1 stream data and the sharp drop of the recent Gaia DR3 rotation curve.

In the text

	Fig. 6. Modeled and observed projected orbit in the sky for a Fermionic-MW DM model, shown as a solid (amber) line, and a BH of Mbh = 4.075 × 106 M⊙, shown as a dotted (green) line.
In the text

	Fig. 7. Redshift (z), right ascension (α), and declination (δ) as a function of time epoch for the same models as displayed in Fig. 6.
In the text

	Fig. 8. Fermionic DM density profiles with different core compacities (i.e., different m) fitting the GD-1 stream and the DM core mass, which agree with the S2 star data orbiting Sagittarius A*. The Schwarzschild radius is computed assuming a BH mass Mbh = 4.075 × 106 M⊙.
In the text