© The Authors 2025

1. Introduction

In the standard cosmological model, dark matter (DM) is made of cold collisionless particles (CDM). They evolve under their own gravity to form halos following the canonical NFW profile (after Navarro, Frenk, and White 1997), where the mass density profile increases with decreasing radius (r) approximately as r⁻¹. These predicted cuspy profiles are seldom observed (e.g., Del Popolo & Le Delliou 2017; Bullock & Boylan-Kolchin 2017; Salucci 2019) since the inferred DM haloes tend to show a constant central density or core. The formation of cores is naturally accommodated within the standard cosmological model since purely baryonic processes move gas around, leading to the transformation of the overall potential and redistribution of the DM particles. In the case of dwarf galaxies, the energy turning DM cusps into cores is provided by star formation (e.g., Governato et al. 2010; Pontzen & Governato 2012); therefore, when the formed stellar mass is too small, the baryon feedback alone cannot transform cuspy DM halos into cored halos and the DM haloes should remain NFW-like. Even though the limiting mass characterizing these halo unevolved galaxies (HUGs) is model dependent (e.g., Read et al. 2016; Koudmani et al. 2024), it approximately corresponds to M_⋆ < 10⁶ M_⊙ (e.g., Peñarrubia et al. 2012; Di Cintio et al. 2014b; Chan et al. 2015; Hayashi et al. 2020; Jackson et al. 2021). Thus, if the DM haloes of galaxies with M_⋆ ≪ 10⁶ M_⊙ show cores, it would indicate the DM not being collisionless but fuzzy, self-interacting, warm, or other alternatives to CDM (Dodelson & Widrow 1994; Hu et al. 2000; Spergel & Steinhardt 2000; Bechtol et al. 2022; Carr et al. 2024).

Traditionally, the DM halo shapes are deduced from spatially resolved kinematical measurements, which require time-consuming high spectral resolution spectroscopy in the optical, infrared, or radio band. Keeping in mind the need for large statistics to reach reliable conclusions, through this approach, it is nearly impossible to measure enough DM halos in the HUG regime to address the DM nature issue. However, the broad-band photometry needed to infer the stellar mass distribution in the HUG regime starts to be achievable (e.g., Trujillo et al. 2021; Carlsten et al. 2021; Richstein et al. 2024; Zaritsky et al. 2024) and will become routinely simple in the near future with instruments such as the Rubin Observatory (Ivezić et al. 2019) or the Euclid satellite (e.g., Laureijs et al. 2011). Fortunately, one can use photometry alone to constrain the DM halo mass distribution using the classical Eddington inversion method (EIM; An & Evans 2006; Ciotti & Morganti 2010; Sánchez Almeida et al. 2023).

The EIM (Eddington 1916; Binney & Tremaine 2008; Lacroix et al. 2018; Ciotti 2021) provides the distribution function (DF) in the phase space needed if an observed mass density profile happens to be immersed in an assumed gravitational potential. If the required DF becomes negative somewhere in the phase space, it proves the observed density to be physically inconsistent with the assumed potential. Such inconsistency between potential and density happens for a particularly interesting combination in the context of deciphering the nature of DM, namely, when a stellar density with a core resides in a NFW potential (see An & Evans 2006; Ciotti & Morganti 2010; Sánchez Almeida et al. 2023). Stellar cores are common in dwarfs (Moskowitz & Walker 2020; Carlsten et al. 2021; Sánchez Almeida et al. 2021; Battaglia & Nipoti 2022; Richstein et al. 2024), and if this fact remains in the critical HUG range, it would evidence the need to go beyond the standard cold DM model. We note that the original inconsistency was worked out for a particularly simplistic combination where the stars form a spherical system with isotropic velocities and residing in a NFW potential. However, this particular case seems to reflect a more general and profound inconsistency since the assumptions can be substantially relaxed and the inconsistency remains: it still holds for (1) quasi stellar cores and quasi NFW potentials, where the central density is not exactly constant and the inner slope of the potential is not −1 (Sánchez Almeida et al. 2023); (2) anisotropic orbits of the type expected in dwarfs (isotropic at the center and radially biased in the outskirts, of the type Osipkov-Merritt or Cuddeford; Ciotti & Morganti 2010; Sánchez Almeida et al. 2023); (3) Einasto potentials, also characteristic of CDM without the mathematical singularity at r = 0 hampering NFW potentials (Sánchez Almeida 2024); and (4) axi-symmetric systems, proving the inconsistency to go beyond the spherical symmetry assumption (Sánchez Almeida et al. 2024a).

A first attempt to constrain the DM halo of real galaxies using EIM was carried out by Sánchez Almeida et al. (2024a). They analyzed noisy and incomplete data of around 100 low-mass (M_⋆ between 10⁶ and 10⁸ M_⊙) satellites of MW-like galaxies taken from Carlsten et al. (2021). Fits to the observed surface brightness profiles were compared with a battery of gravitational potentials including NFW potentials as well as potentials stemming from cored mass distributions (expected in many alternatives to CDM). The method requires fitting the observed profiles with several analytic functions with variable inner and outer slopes. Between 40% and 70% of the galaxies are consistent with pure cores in the stellar mass distribution and thus inconsistent with NFW-like potentials. Unfortunately, the fitted galaxies are still too massive to be in the HUG regime and thus to be conclusive on the nature of the DM issue. Moreover, this first technique has the drawback of not providing the DF that best fit the data but just pointing out incompatibilities.

A second alternative approach was followed by Sánchez Almeida et al. (2024b), using EMI to claim deviations of the real DM from the CDM paradigm. They analyzed six ultra-faint dwarfs in the HUG regime (M_⋆ ∼ 10³–10⁴ M_⊙), all showing a clear stellar core incompatible with a NFW potential and fully compatible with the cored potentials predicted by many alternatives to CDM. Also based on the EIM, they used a new method to directly compute the DF by fitting the observed mass surface density profile as a superposition of density profiles characteristic of the assumed potential. In essence, each potential defines a set of density profiles, and it is the expansion of the observed profile in this database that provides the DF. The new method is outlined in the Letter by Sánchez Almeida et al. (2024b), but we feel compelled to describe the procedure in detail, which is the main purpose of the present paper. The technique is for general application, provided the underlying assumptions are met, and we have chosen the galaxy Nube, recently discovered by Montes et al. (2024), to showcase the operation and difficulties of the new tool. The choice of Nube is not accidental. Although they have similar stellar masses, it is seven times larger and 100 times dimmer that the Small Magellanic Cloud. Nube is an outlier of most global scaling relations and does not seem to be predicted by modern CDM cosmological simulations. Thus, the characterization of the DM halo of Nube is particularly interesting in the context of CDM tests.

The paper is organized as follows: Sect. 2 describes the galaxy Nube and the stellar mass surface density data analyzed in the present work. Part of this description includes a Monte Carlo simulation (Appendix A) indicating how the determination of the center of the galaxy does not influence the inferred stellar mass density profile. This section also includes fitting the observed profile with several analytic forms (Sect. 2.1). Section 3 puts forward the method to derive the DF. The general mathematical formulation (main Sect. 3) is specified to particular potentials in Sect. 3.1 (Schuster-Plummer potential), Sect. 3.2 (NFW potential), and Appendix B (abc potential). The interpretation of our formulation in terms of the classical statistical mechanics of gravitating systems is carried out in Appendix C. The actual implementation of the algorithm, which follows a Bayesian approach, is presented in Sect. 3.3, with a number of sanity checks presented in Appendices E.1, E.2, and E.3. Section 4 describes the application to Nube showing how NFW potentials (cuspy) are highly disfavored compared with Schuster-Plummer potentials (cored). It is split into two subsections, showing the constraints imposed by the analytic function fits (Sect. 4.1) and the actual application of the new algorithm (Sect. 4.2). A short Sect. 5 outlines simple extensions of the algorithm that relax the assumption of isotropic velocities. Finally, the conclusions are summarized in Sect. 6 and general guidelines to improve the method are given.

2. Data and first analysis of Nube

Nube was discovered and characterized by Montes et al. (2024). It has a stellar mass similar to the Small Magellanic Cloud but happens to be unprecedentedly large. Nube is 7 times larger (effective radius of 6.9 kpc) and 100 times fainter (V-band central surface brightness of 26.2 mag arcsec⁻²) than the typical galaxies of its mass (M_⋆ ≃ 3.9 × 10⁸ M_⊙). This makes Nube the most diffuse object of its class and, at present, a dwarf without any clear counterpart in the CDM cosmological simulations aimed at reproducing ultra-diffuse galaxies (Montes et al. 2024, and references therein). Thus, the DM halo of Nube is particularly interesting so that Nube represents an excellent target to test the usefulness of the new tools.

For the sake of comprehensiveness, here we summarize the main properties of the observation and reduction. The stellar mass profile in Fig. 1 comes from a series of deep images in the Sloan u, g, r, i and z filters taken with HiPERCAM (Dhillon et al. 2018, 2021) operated at the 10-m GTC telescope (Gran Telescopio Canarias). The apparent size of Nube is smaller than the field of view of the camera, allowing a reliable background subtraction. Simultaneously, it is large enough to grant the inner core of the galaxy (∼10 arcsec) to be well resolved. After observing for ∼70 min, the g-band image reaches a surface brightness limit ∼31 mag arcsec⁻² (3σ in areas equivalent to 10 × 10 arcsec²). The surface density stellar mass was inferred from the photometry in the g-band, with a mass-to-light ratio inferred from g − r using the calibration in Roediger & Courteau (2015), which assumes a Chabrier (2003) initial mass function. The stellar mass surface density profile, Σ(R), was computed from ring averages at different radial distances, up to 30 arcsec from the center of the galaxy.

Fig. 1. Stellar mass surface density profile of Nube as observed by Montes et al. (2024, the symbols with error bars). The figure includes three different fits to this observation: a projected polytrope (PP; the black line, with its index m given in the inset) and two projected abc profiles (Eq. (1)). The inner slope (−c) is forced to be zero in the PP fit whereas becomes positive (i.e., c < 0) when allowed to vary (see the inset). The gray error bars are those provided by Montes et al. (symmetric in a linear scale and so asymmetric in the logarithmic representation used in the figure) whereas the symmetrized ones1 (the blue bars) are equivalent but symmetric in the logarithmic representation.

The resulting stellar mass surface density profile is shown in Fig. 1. The error bars were calculated as a combination of Poisson noise and errors involved in subtracting the background (for further details, see Montes et al. 2024). The individual points in the profile come from averages in rings that do not overlap; therefore, their errors are independent. Since we will fit Σ(R) in a logarithmic scale, for convenience, the original error bars provided by Montes et al. (2024) are used after symmetrization in logΣ¹. Both the original error bars (in gray color) and the symmetric ones (in blue) are given in Fig. 1.

The derivation of the profile shown in Fig. 1 depends on several assumptions, whose influence on the analysis presented in the following sections was evaluated and determined to be secondary. The selection of the galaxy center was examined and found to have a negligible impact on Σ(R) (Appendix A). The uncertainties in the error bar estimate was addressed in Appendix E.2 where we also study the effect of changing the mass-to-light ratio.

2.1. Fitting Nube with simple profiles

The first step of the analysis was fitting Σ(R) with simple profiles. The result is also included in Fig. 1. Polytropes are expected to describe the mass distribution when self-gravitating systems reach either thermodynamic equilibrium or another long lasting meta-stable state (Plastino & Plastino 1993; Sánchez Almeida et al. 2020), and they reproduce the mass distribution in many practical instances (Sánchez Almeida et al. 2020, 2021). We use the tool described in Sánchez Almeida et al. (2021) to fit Nube with a projected polytrope (PP). The best PP fit (the black line in Fig. 1) does a good job reproducing the observation. It is important to note, however, that polytropes have cores (i.e., dρ/dr → 0 when r → 0) and so they are unable to follow the mild but noticeable drop of the mean Σ(R) towards the innermost radii of Nube (see Fig. 1). To be able to reproduce this drop, we also tried with plane of the sky projections of abc profiles, defined as

$$\begin{array}{c}\hfill {\rho }_{\mathit{abc}}=\frac{{\rho }_s}{x^c{(1+x^a)}^{(b-c)/a}},\end{array}$$(1)

where x = r/r_s, and ρ_s and r_s are scaling constants setting the volume density and the size, respectively. These profiles are commonly used to model the density of baryons or DM (e.g., Hernquist 1990; Merritt et al. 2006; Di Cintio et al. 2014a) and have the advantage of encompassing the iconic NFW profile (a = 1, b = 3,  and c = 1) and the m = 5 polytrope (a.k.a. Schuster-Plummer profile, with a = 2, b = 5,  and c = 0). Note that −c and −b give the logarithmic slope of the profile in the inner and outer radii, respectively. We fit the mass surface density of Nube with c as a free parameters and setting a = 2 − c and b = 5 − 2c, which allows the profile to seamlessly scan from a Schuster-Plummer profile to a NFW profile when c varies from 0 to 1. The result is shown as the orange line in Fig. 1, which has c ≃ −0.15 and improves the root-mean-square (RMS) of the residuals with respect to the PP fit (see the inset in Fig. 1). We also try fits allowing both the inner and outer slopes c and b to vary (the red line in Fig. 1, which assumes a = 2, fixed to minimize the number of free parameters). The fit is even better, also yielding a positive inner slope (c ≃ −1). The fact that the inner slope tends to be positive is makes it difficult to reproduce Σ(R) self-consistently within any potential, as we will discuss in detail in Sect. 4.

3. Method to derive the distribution function f(ϵ)

For a spherically symmetric system of particles with isotropic velocity distribution, the phase-space DF f(ϵ) depends only on the particle energy ϵ. Then, the volume density ρ(r) turns out to be (e.g., Binney & Tremaine 2008, Sect. 4.3)

$$\begin{array}{c}\hfill \rho (r)=4\pi \sqrt{2}{\displaystyle {\int }_0^{\mathrm{\Psi }(r)}}f(ϵ)\sqrt{\mathrm{\Psi }(r)-ϵ}dϵ,\end{array}$$(2)

with $ϵ=\mathrm{\Psi }-\frac{1}{2}{v}^2$ the relative energy per unit mass of each particle, v the particle’s velocity, and Ψ(r) = Φ₀ − Φ(r) the relative potential energy. The symbol Φ(r) stands for the gravitational potential energy and Φ₀ is the gravitational potential energy evaluated at the edge of the system. The previous equation can be rewritten as

$$\begin{array}{c}\hfill \rho (r)={\displaystyle {\int }_0^{{ϵ}_{max}}}f(ϵ)\xi (ϵ,r)dϵ,\end{array}$$(3)

with

$$\begin{array}{c}\hfill \xi (ϵ,r)=4\pi \sqrt{2}\sqrt{{ϵ}_{\max }}\sqrt{\frac{\mathrm{\Psi }(r)}{\mathrm{\Psi }(0)}-\frac{ϵ}{{ϵ}_{max}}}\mathrm{\Pi }(X-r),\end{array}$$(4)

where ϵ_max = Ψ(0), X is the radius implicitly defined as Ψ(X)/Ψ(0) = ϵ/ϵ_max, and Π(x) represents the step function,

$$\begin{array}{c}\hfill \mathrm{\Pi }(x)=\{\begin{array}{cc}0\hfill & \mathrm{if}x<0,\hfill \\ 1\hfill & \mathrm{if}x\ge 0.\hfill \end{array}\end{array}$$(5)

Equation (3) admits a physically revealing interpretation. The function ξ(ϵ, r) parameterizes a family of energy dependent volume densities characteristic of the potential Ψ². Then, the volume density is just the superposition of these other characteristic densities with the DF parameterizing the contribution of each energy (see Eq. (3)). Examples of these characteristic densities for NFW and Schuster-Plummer potentials are given in Figs. 2 and 3, with their derivation worked out in Sects. 3.2 and 3.1, respectively. The general case of an abc potential is treated in Appendix B.

Fig. 2. Characteristic density profiles corresponding to each relative energy α = ϵ/ϵmax in a Schuster-Plummer potential (the solid lines; Eq. (16)) and in a NFW potential (the dashed lines; Eq. (23)). The symbol β stands for the normalized radial coordinate r/rsp. All energies contribute to the innermost regions (β ≪ 1) whereas only the smallest energies contribute to the outer halo (β ≫ 1).

Fig. 3. Similar to Fig. 2 but representing the contribution to the total density of the different energies (α) at each constant radial position (β).

Integrating Eq. (3) over all the volume, the total mass of the system M turns out to be an integral of the masses corresponding to the different ϵ, explicitly,

$$\begin{array}{c}\hfill M={\displaystyle {\int }_0^{{ϵ}_{max}}}f(ϵ)\varpi (ϵ)dϵ,\end{array}$$(6)

with

$$\begin{array}{c}\hfill \varpi (ϵ)=4\pi {\displaystyle {\int }_0^{\infty }}\xi (ϵ,r)r^{2}dr.\end{array}$$(7)

The variable ϖ(ϵ) is shown in Fig. 4 for a Schuster-Plummer potential and a NFW potential. As we will show in Sects. 3.1 and 3.2 (and appears in Fig. 4) ϖ(ϵ) diverges when ϵ → 0 thus posing some general restriction on f(ϵ) when ϵ ≪ ϵ_max. We note that the quantity ϖ(ϵ) coincides with the quantity that in classical statistical mechanics is known as the density of states. For the sake of clarity, the connections of our approach with the classical interpretation are pointed out and discussed in Appendix C.

Fig. 4. Specific mass corresponding to each energy (Eq. (7) with α = ϵ/ϵmax) for Schuster-Plummer and NFW potentials. Top panel: log-linear representation. Bottom panel: log-log representation. Note that the specific mass diverges for low energies as a power-law of exponent −2.5, for the Schuster-Plummer potential, and around −2.81, for the NFW potential (the dashed lines labeled in the inset of the bottom panel).

In principle, f(ϵ) could be retrieved using Eq. (3) by fitting ρ(r) with a linear superposition of ξ(ϵ, r). In practice, however, there is no unique way to discretize Eq. (3) for such purpose. We approach the practical problem expanding f(ϵ) as a polynomial of order n,

$$\begin{array}{c}\hfill f(ϵ)\simeq {\displaystyle \sum _{i=3}^n}\frac{a_i}{{ϵ}_{max}^{3/2}}{(ϵ/{ϵ}_{max})}^i,\end{array}$$(8)

so that

$$\begin{array}{c}\hfill \rho (r)\simeq {\displaystyle \sum _{i=3}^n}{a}_i{F}_i(r),\end{array}$$(9)

with

$$\begin{array}{c}\hfill F_i(r)={\displaystyle {\int }_0^1}{\alpha }^i\frac{\xi (\alpha {ϵ}_{max},r)}{\sqrt{{ϵ}_{max}}}d\alpha .\end{array}$$(10)

It is important to note that the polynomial expansion in Eq. (8) lacks the three first terms (it begins at i = 3). This is a constraint imposed by the need to have a finite total mass since for ϵ → 0, ϖ(ϵ) diverges as ϵ^−γ with 2 < γ < 3 (see Fig. 4 and Sects. 3.1 and 3.2). We also note that the normalization in Eq. (10) was chosen so that F_i(r) does not depend on ϵ_max (see Eq. (4)). The discretization also holds for the projection of the 3D densities in the plane of the sky, that is,

$$\begin{array}{c}\hfill \mathrm{\Sigma }(R)\simeq {\displaystyle \sum _{i=3}^n}{a}_i{S}_i(R),\end{array}$$(11)

$$\begin{array}{c}\hfill S_i(R)={\displaystyle {\int }_0^1}{\alpha }^i\frac{{\xi }_{\mathrm{\Sigma }}(\alpha {ϵ}_{max},R)}{\sqrt{{ϵ}_{max}}}d\alpha ,\end{array}$$(12)

where Σ(R) and ξ_Σ(ϵ, R) are the 2D projection (the Abel transform) of ρ(r) and ξ(ϵ, r), respectively. The symbol R stands for the radial coordinate in the plane of the sky projection, as used in Sect. 2 to describe the observation of Nube.

In the next subsections, the general expressions given above are particularized for the extreme potentials used in the work, namely, a NFW potential (Sect. 3.1) and a Schuster-Plummer potential (Sect. 3.2). The general case of the potential created by an abc density is worked out in Appendix B.

3.1. Case of a Schuster-Plummer potential

The polytrope of index 5 is usually called Schuster-Plummer profile,

$$\begin{array}{c}\hfill {\rho }_{\mathrm{SP}}(r)=\frac{{\rho }_{\mathit{sp}}}{{[1+{(r/r_{\mathit{sp}})}^2]}^{5/2}},\end{array}$$(13)

where ρ_sp and r_sp are the central density and the characteristic radial scale, respectively³. The potential produced by this cored density profile is the following (e.g., Sánchez Almeida et al. 2023, Eq. (A.14)):

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_{\mathrm{SP}}(r)=\frac{{\mathrm{\Psi }}_{\mathrm{SP}}(0)}{{[1+{(r/r_{\mathit{sp}})}^2]}^{1/2}},\end{array}$$(14)

with Ψ_SP(0) = ϵ_max = 4πGρ_spr_sp²/3. Using these values, Eq. (4) renders

$$\begin{array}{c}\hfill {\xi }_{\mathrm{SP}}(ϵ,r)=\sqrt{{ϵ}_{max}}h(ϵ/{ϵ}_{max},r/r_{\mathit{sp}}),\end{array}$$(15)

with

$$\begin{array}{c}\hfill h(\alpha ,\beta )=4\pi \sqrt{2}\sqrt{[{(1+{\beta }^2)}^{-1/2}-\alpha ]}\mathrm{\Pi }({\beta }_X-\beta ),\end{array}$$(16)

and β_X defined as

$$\begin{array}{c}\hfill {\beta }_X=\sqrt{1-{\alpha }^2}/\alpha .\end{array}$$(17)

Figure 2 shows h(α, β) as a function of the radial coordinate β for a number of energies α (the solid lines). All energies contribute to the innermost regions (β ≪ 1) whereas only the smallest energies contribute to the outer halo (β ≫ 1). Figure 3 also shows h(α, β) but this time as a function of the energy for a number of radii (the solid lines).

The mass corresponding to each relative energy (Eq. (7)) happens to be

$$\begin{array}{c}\hfill \varpi (\alpha {ϵ}_{max})=4\pi r_{\mathit{sp}}^3\sqrt{{ϵ}_{max}}{\displaystyle {\int }_0^{\infty }}h(\alpha ,\beta ){\beta }^{2}d\beta .\end{array}$$(18)

This mass diverges at low energies since one can prove that

$$\begin{array}{c}\hfill {\displaystyle \underset{\alpha \to 0}{lim}}\varpi (\alpha {ϵ}_{max})\simeq r_{\mathit{sp}}^3\sqrt{{ϵ}_{max}}\frac{\sqrt{2}{\pi }^3}{{\alpha }^{5/2}}.\end{array}$$(19)

The solid orange line in Fig. 4 shows ϖ(ϵ) computed numerically from ξ_SP using the Simpson’s rule form Scipy (Virtanen et al. 2020). The numerical calculation was tested against the low-energy trend given by Eq. (19), which is also shown in Fig. 4 as the dashed orange line.

3.2. Case of a NFW potential

In the case of a NFW density setting the potential,

$$\begin{array}{c}\hfill {\rho }_{\mathrm{NFW}}(r)=\frac{{\rho }_{\mathit{sp}}}{(r/r_{\mathit{sp}}){(1+r/r_{\mathit{sp}})}^2},\end{array}$$(20)

so that the potential turns out to be (e.g., Sánchez Almeida et al. 2023, Eq. (A.7))

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_{\mathrm{NFW}}(r)={\mathrm{\Psi }}_{\mathrm{NFW}}(0)\frac{ln(1+r/r_{\mathit{sp}})}{r/r_{\mathit{sp}}},\end{array}$$(21)

with Ψ_NFW(0) = ϵ_max = 4πGρ_spr_sp². Using this gravitational potential, Eq. (4) renders

$$\begin{array}{c}\hfill {\xi }_{\mathrm{NFW}}(ϵ,r)=\sqrt{{ϵ}_{max}}g(ϵ/{ϵ}_{max},r/r_{\mathit{sp}}),\end{array}$$(22)

where

$$\begin{array}{c}\hfill g(\alpha ,\beta )=4\pi \sqrt{2}\sqrt{[ln(1+\beta )/\beta -\alpha ]}\mathrm{\Pi }({\beta }_X-\beta ),\end{array}$$(23)

with β_X implicitly defined as

$$\begin{array}{c}\hfill ln(1+{\beta }_X)/{\beta }_X=\alpha .\end{array}$$(24)

The characteristic density g(α, β) is shown as a function of the radial coordinate β in Fig. 2 (the dashed lines) and as a function of the energy α in Fig. 3 (the dashed lines).

The mass corresponding to each relative energy (Eq. (7)) would be

$$\begin{array}{c}\hfill \varpi (\alpha {ϵ}_{max})=4\pi r_{\mathit{sp}}^3\sqrt{{ϵ}_{max}}{\displaystyle {\int }_0^{\infty }}g(\alpha ,\beta ){\beta }^{2}d\beta ,\end{array}$$(25)

which is shown in Fig. 4 (the blue solid line). (The numerical integration scheme is the same as for the Schuster-Plummer potential sketched in Sect. 3.1.) As it happens with the Schuster-Plummer potential, ϖ(ϵ) diverges when ϵ → 0. In this case, the numerical integration gives a mass that approximately scales as ϵ^−2.81 (the blue dashed line in Fig. 4).

3.3. Algorithm to infer f(ϵ) from Σ(R)

Except for the arbitrary scaling parameterized by ϵ_max, Eqs. (8) and (11) provide a method to infer the DF f(ϵ) needed for a galaxy of observed mass surface density Σ(R) to reside in a given gravitational potential. A fitting algorithm using Eq. (11) provides the coefficients a_i determining f(ϵ) through Eq. (8). The characteristic densities in Eq. (12) have to be computed numerically starting from the potential in a chain requiring at least two integrations: the Abel transform that projects the volume densities on the plane of the sky and the integral over all energies expressed by Eq. (12). We compute the Abel transform using the direct method implemented in the PyAbel Python package (Hickstein et al. 2019). Then the 2nd integration is carried out using the Simpson’s rule from Scipy (Virtanen et al. 2020). Several of the monomials for the Schuster-Plummer and NFW potentials are given in Fig. 5. All functions S_i(R) show a central plateau with a power-law drop in the outskirts being more steep as the index of the monomial increases. The numerical method was tested using the analytic solution for $S_{\frac{7}{2}}(R)$ worked out in Appendix D.

Fig. 5. Characteristic functions corresponding to the polynomial expansion of the DF: see Eqs. (8), (11), and (12). The symbol i stands por the exponent of the monomial, and the solid and dashed lines correspond to the Schuster-Plummer and the NFW potentials, respectively. It begins at i = 3 since for i ≲ 3 the total mass of the resulting density profile diverges (see main text).

The free parameters retrieved from fitting Σ(R) are the amplitudes a_i together with the global radial scaling factor setting the width of the potential r_sp, the latter making the fit nonlinear. The fits were carried out using a Bayesian approach, with the log-likelihood defined as −χ²/2 where

$$\begin{array}{c}\hfill {\chi }^2={\displaystyle \sum _j}{\left[\frac{log\mathrm{\Sigma }(R_j)-log{\mathrm{\Sigma }}_m(R_j)}{\mathrm{\Delta }log\mathrm{\Sigma }(R_j)}\right]}^2,\end{array}$$(26)

with Σ(R_j) the observed Σ(R) at the j-th radial position, ΔlogΣ(R_j) its error, and Σ_m(R_j) the corresponding model at R_j. The sum includes all radii. The posterior is explored using the ensemble sampler for Markov Chain Monte Carlo (MCMC) emcee (Foreman-Mackey et al. 2013). The best fit provided by a least squares routine that minimizes χ² was used to initialize the exploration (least_squares from scipy; Virtanen et al. 2020). We carried out a first unconstrained fit, allowing f(ϵ) to vary freely. This is used as reference in all the forthcoming discussion. However, the exploration was initialized forcing the least squares routine to yield physically sensible solutions with f(ϵ)≥0 for all ϵ. Several trial-and-error tests led us to set the final hyper-parameters used for fitting as described below. The order of the polynomial was chosen to n = 10, large enough to provide the flexibility needed to reproduce the inner plateau observed in Nube (Fig. 1). The priors in the Bayesian analysis were chosen to be as uninformative as possible. The radial scaling r_sp was forced to be non-negative and r_sp ≤ 10³ kpc whereas the amplitudes a_i were let to vary unconstrained. Also as a prior, we asked the outermost slope of the fitted logΣ(R) to be less than -2, thus preventing Σ(R) to have infinite mass outside the observed radii. In addition, we force f(ϵ)≥0. The posterior was explored with 32 walkers and 6000 samples.

The results reported in this paper do not depend on the exact values of the used hyper-parameters, as deduced from a number of tests detailed in Appendix E.1. In these sanity checks, the hyper-parameters are modified to assess the effect on the interpretation of Nube. Explicitly, we tried: (1) initializations not forced to have f(ϵ)≥0, (2) changing the number of walkers and samples, (3) changing the order of the polynomial used for f(ϵ) (n in Eq. (8)), (5) constraining the amplitudes a_i relative to the values of the least squares best fit, and (6) other variations referred to uncertainties in the data itself. See Appendix E.1 for a complete account.

4. The gravitational potential of Nube

4.1. Constraints from simple fits to Σ(R)

The observation of Nube in Sect. 2.1 shows a seemingly positive inner slope,

$$\begin{array}{c}\hfill \omega ={\displaystyle \underset{R\to 0}{lim}}\frac{dlog\mathrm{\Sigma }(R)}{dlogR}>0,\end{array}$$(27)

which is well reproduced by profiles that also have positive inner slopes in their 3D mass distribution. In the case of ρ_abc profiles (Eq. (1)) this is achieved with c < 0 (Fig. 1, the red and the orange lines). The so-called cusp slope-central anisotropy theorem by An & Evans (2006) (see also Ciotti & Morganti 2010; Sánchez Almeida et al. 2023) applies to spherically symmetric systems with constant velocity anisotropy β_u and it states that

$$\begin{array}{c}\hfill c\ge 2{\beta }_u,\end{array}$$(28)

with

$$\begin{array}{c}\hfill {\beta }_u=1-\frac{{\sigma }_t^2}{2{\sigma }_r^2},\end{array}$$(29)

where σ_r and σ_t are the radial and tangential velocity dispersions, respectively. The form of the EIM adopted in this work (Eq. (2)) assumes β_u = 0 so that no potential is able to reproduce the required inner positive slope (i.e., c < 0 is inconsistent with β_u = 0, according to Eq. (28)). This is a issue that permeates the analysis presented in subsequent sections. Obviously, the observational error bars are so large that they allow for c = 0 or even slightly positive slope (see Fig. 1, the green line), but the fact that the best least squares fit is unphysical marks the study. Expanding on this argument, the light profile of Nube in different colors do not show the drop (Fig. 4 in Montes et al. 2024, and also Fig. 3 below), which seems to appear when transforming the observed photometry into stellar mass. The impact of the uncertainty in the used mass-to-light ratio is analyzed in Appendix E.2.

4.2. Constraints from the full DFs

Using the procedure described in Sect. 3 and detailed in Sect. 3.3, we fit the Σ(R) of Nube assuming two extreme gravitational potentials; one with a core (a Schuster-Plummer potential) and another with a cusp (a NFW potential). In addition, we also fit Nube assuming the gravitational potential to be generated by a ρ₂₃₀ profile (Eq. (1) with a = 2, b = 3 and c = 0), which has a core but with the outskirts approaching a NFW profile (see Fig. 10).

The result for the cored Schuster-Plummer potential are included in Fig. 6. The red solid line represents the least-squares best fit allowing for any f(ϵ) whereas the other thin lines are the fits derived from the MCMC exploration of the posterior, forcing f(ϵ)≥0 to render physically sensible solutions and initialized with f(ϵ)≥0 least squares solution. We note that the best fit to the observed Σ(R) is very good but requires f(ϵ) < 0 when ϵ → ϵ_max. This fact reflects the issue discussed in Sect. 4.1 that the assumed isotropic velocity distribution is not consistent with positive inner slopes (Eq. (27)). The same exercise with a NFW potential is included in Fig. 7. The best fit is as good as the one for the Schuster-Plummer potential. The arrows in Fig. 8c show the value of the merit function χ² (Eq. (26)) of the best fits, and both are very similar (cf. the blue and orange arrows). However, the physically meaningful fits forcing f(ϵ)≥0 are significantly worst in the case of a NFW potential; compare Figs. 7a and 6a, the orange and blue points in Fig. 8a, and the orange and blue histograms in Fig. 8c. Moreover, the innermost slope ω becomes too negative in the case of the NFW profile to be compatible with zero (Fig. 8b). For the sake of completeness, Fig. 9 also includes the fits resulting from assuming a ρ₂₃₀ potential. The corresponding histograms of ω and χ² are in Fig. 8, the green points and histograms. These fits are similar and of similar quality as the fits provided by the Schuster-Plummer potential (cf. the orange with the green points and histograms in Fig. 8).

Fig. 6. Diagnostic plot for the technique. (a) Fit to the stellar mass surface density observed in Nube (the blue symbols) assumed to reside in a Schuster-Plummer potential. The red thick solid line represents the least-squares best fit with an unconstrained DF f(ϵ), whereas the other thin lines are the fits derived from the MCMC exploration of the posterior forcing f(ϵ) to be ≥0∀ϵ and beginning the exploration from the f(ϵ)≥0 least squares solution. These other fits are color coded according to the innermost slope of the surface density profile (ω in Eq. (27)), which always happens to be negative but close to zero. (b) DFs required by the best fit and by the fits forced to have f(ϵ)≥0. Note how the best fit requires f(ϵ) < 0 and so is unphysical. The color code is the same as in (a).

Fig. 7. Similar to Fig. 6 except that Nube is assumed to reside in a NFW potential. The best fitting function (the red thick solid line) is similar to the best fit assuming a Schuster-Plummer potential (Fig. 6a), but the physically realizable f(ϵ)≥0 fits are clearly worst and have inner slopes differing from zero significantly (see the color code of the thin lines).

Fig. 8. (a) Scatter plot χ2 versus innermost slope ω of the fits for the three potentials resulting from the MCMC navigation of the posterior. The large times symbols represent the best fit obtained with unconstrained f(ϵ), which yield unphysical f(ϵ) < 0. The vertical gray line marks ω = −0.01. (b) Histograms of the distribution of innermost slopes for the points in (a). (c) Histograms of χ2 for the points in (a). The arrows represent the χ2 of the best fit obtained with unconstrained f(ϵ). Panels a, b, and c share the same color code as indicated in the insets.

Fig. 9. Similar to Fig. 6 except that Nube is assumed to reside in a ρ230 potential, which has a core and approaches a NFW profile in the outskirts. The best fitting function (the red thick solid line) is similar to the best fit in a NFW potential (Fig. 7), but the physically sensible f(ϵ)≥0 are clearly better and have a inner slope closer to zero.

A way of quantifying the significant difference between the physically sensible fits based on core (Schuster-Plummer and ρ₂₃₀) and cusp (NFW) potentials comes from the histograms in Fig. 8. In the case of the NFW potential, only 0.5% of the points exploring the posterior have inner slope reasonably close to zero (>  − 0.01; Fig. 8a). This fraction increases to 76% and 66% for the Schuster-Plummer and ρ₂₃₀ potentials, respectively. Similarly, the mean value of the χ² distribution is around 5 for Schuster-Plummer and ρ₂₃₀ and more than 8 for the NFW potential based fits. Finally, the minimum χ² is 1.9 for both the Schuster-Plummer and the ρ₂₃₀ potentials whereas it is 1.5 times larger for the NFW potential. If this minimum χ² for NFW is assumed to follow a χ² probability distribution function with 11 degrees of freedom⁴, reaching a value 1.5 times larger than the true value, assumed to be set by the Schuster-Plummer χ², has a small probability of 11%. In other words, the probability that the best Schuster-Plummer f(ϵ)≥0 fit is better than the best NFW f(ϵ)≥0 fit is around 89%. All these statistical tests combined indicate that Nube is much more likely to reside in a gravitational potential with a core than with a cusp.

A number of sanity checks support that the method presented in Sect. 3.3 works as expected when applied to known profile-potential pairs, thus supporting the above conclusions. They are discussed in Appendix E. The impact of the assumed hyper-parameters of the bayesian fit, including the priors, are analyzed in Appendix E.1. The effect of the uncertainties in the error assigned to Nube, which enter into the definition of the merit function (Eq. (26)), is examined in Appendix E.2. We repeat te analysis considering only photometric errors, constant mass-to-light ratios, and rearranging the true surface density profile to force a monotonic decrease of Σ(R) in the inner region. None of these modifications alter than main conclusion that the resulting fits are quite good for a Schuster-Plummer potential and outrageous for a NFW potential.

One final outcome of the analysis is the radial extension of the potential parameterized by the characteristic radius r_sp. Considering the MCMC exploration of the posterior, the range of values is log(r_sp/1 kpc)≃0.79 ± 0.12 for the Schuster-Plummer potential, 0.47 ± 0.17 for the NFW potential, and 0.54 ± 0.18 for the ρ₂₃₀ potential. We note that their values cannot be compared directly since they are just scaling factors of different functional forms. One can compare the different potentials through the surface density giving rise to them. This comparison is included in Fig. 10 which, together with the profile of Nube, shows the mass surface density profiles corresponding to the three alternative potentials; the Schuster-Plummer potential (the blue line), the NFW potential (the green line), and the ρ₂₃₀ potential (the magenta line). Even though the spatial scaling of the potentials is set by the fit, the vertical scaling (i.e., the depth of the potential or its total mass) is arbitrary, and it was arbitrarily chosen in Fig. 10 to yield a stellar mass 30 times the stellar mass of Nube, which has M_⋆ ≃ 3.9 × 10⁸ M_⊙ as inferred from the best fitting profiles.

Fig. 10. Comparison between the stellar surface density profile of Nube and the best-fitting potentials with f ≥ 0. The plot shows the mass surface density that gives rise to the best-fitting Schuster-Plummer potential (blue line), NFW potential (green line), and ρ230 potential (magenta line). The spatial scaling of the potentials is set by the fit whereas the vertical scaling is arbitrary, and it was arbitrarily chosen to be 30 times the stellar mass of Nube.

We can also compare the observed density profile and the potentials using the core radius. Defined as the radius where the surface density is 1/2 the maximum value, Σ(R_c) = Σ(0)/2, it yields

$$\begin{array}{c}\hfill log[R_{\mathit{cp}}/R_c]=\{\begin{array}{cc}0.06\pm 0.08\hfill & \mathrm{Schuster}-\mathrm{Plummer},\hfill \\ -0.01\pm 0.13\hfill & {\rho }_{230},\hfill \end{array}\end{array}$$(30)

where R_c and R_cp stand for the core radius of the stars and the corresponding potential, respectively. (The NFW potential does not show a core and therefore is not included in Eq. (30).) The error bars comes from the standard deviation of the MCMC sampling of the posterior. We note that the core radius of stars and potential are similar.

5. Simple extensions of the present formalism for anisotropic velocities

Even though they are not used in the present work, the formalism for isotropic velocity distributions detailed in Sect. 3 still holds, mutatis mutandis, for two particularly interesting anisotropic velocity distribution. In the so-called Osipkov-Merrit model, the velocity is assumed to have an anisotropy similar to the one expected in isolated dwarfs, with the orbits nearly isotropic in the center and tending to be radial in the outskirts (e.g., El-Badry et al. 2017; Orkney et al. 2023). In this case one can write down (e.g., Binney & Tremaine 2008) an equation formally identical to Eq. (2) replacing the energy ϵ with Q,

$$\begin{array}{c}\hfill Q=ϵ-\frac{L^2}{2r_b^2},\end{array}$$(31)

and the volume density ρ(r) with

$$\begin{array}{c}\hfill {\rho }_{\mathrm{OM}}(r)=(1+\frac{r^2}{r_b^2})\rho (r),\end{array}$$(32)

where L stands for the norm of the angular momentum vector and r_b is the characteristic radius where the transition between isotropic orbits (center) and radial orbits (outskirts) occurs. Thus, the tool developed in Sect. 3.3 can be applied directly to ρ_OM to retrieve f(Q). It would need to assume a value for r_b but, given the speed of the fitting procedure, one can easily treat it as an hyper-parameter and derive f(Q) given r_b.

Something similar happens with the case of an arbitrary but constant velocity anisotropy β_u (defined in Eq. (29)). In this other case, the DF can be split as

$$\begin{array}{c}\hfill f(ϵ)=L^{-2{\beta }_u}{f}_{ϵ}(ϵ),\end{array}$$(33)

and Eq. (2) can be replaced with the following (e.g., Binney & Tremaine 2008; Sánchez Almeida et al. 2023):

$$\begin{array}{c}\hfill r^{2{\beta }_u}\rho (r)=\kappa {\displaystyle {\int }_0^{\mathrm{\Psi }(r)}}\frac{f_{ϵ}(ϵ)}{{[\mathrm{\Psi }(r)-ϵ]}^{{\beta }_u-1/2}}dϵ,\end{array}$$(34)

where κ is a positive numerical value independent of the radius. Note that this equation is formally quite similar to Eq. (2) provided β_u < 1/2, and so can be used to constraint f(ϵ) assuming a value for β_u < 1/2. Obviously, the equivalent to the characteristic densities (Eq. (4)) would have to be recomputed according to β_u, but preparing a battery of these functions for different β_u is doable.

6. Conclusions

We present a new technique to constrain some properties of the gravitational potential of a galaxy only from photometry, that is, only from the distribution of stellar mass inferred from the observed starlight. Under a number of simplifying assumptions (spherical symmetry, isotropic velocity distribution, identical stars, and stationarity), the classical EIM (Sect. 3) allows us to infer the DF in the phase space needed for the observed stars to reside in an assumed gravitational potential. Thus, gravitational potential and starlight can be shown to be inconsistent if the required DF is negative somewhere. This seemingly simple idea has a far-reaching diagnostic capability in the context of understanding the nature of DM. The gravitational potential expected from CDM is inconsistent with the central plateau or core often observed in the starlight distribution of dwarf galaxies (Sect. 1).

This new technique allowed Sánchez Almeida et al. (2024b) to point out possible deviations from the CDM paradigm. The implementation as a specific tool is detailed in Sect. 3, with several consistency tests provided in Appendix E. Sections 3, 3.1, and 3.2 spell out the mathematical formulation, whereas Sect. 3.3 describes the actual numerical implementation of the general technique using a Bayesian approach. The low-surface brightness dwarf galaxy Nube recently discovered by Montes et al. (2024) has, among other properties, a conspicuous and large inner core (Fig. 1 and Sect. 2) used by Montes et al. (2024) to work out the constraints on fuzzy DM imposed by the existence of such a core. Nube is used in this paper to showcase the application of our tool and to illustrate the kind of physical information it provides.

The actual application to Nube is described in Sect. 4. Provided Nube complies with the assumptions underlying EMI, cuspy NFW potentials are strongly disfavored compared to those with cores (Schuster-Plummer or ρ₂₃₀). As we explain in Sect. 4.1, the mild inner positive slope of Nube (Fig. 1 and Eq. (27)) cannot be reproduced with f(ϵ)≥0. However, when f(ϵ) is forced to be positive, the resulting fits assuming a core potential have smaller χ² and innermost slopes closer to zero than the fits assuming NFW potentials. According to the statistical tests carried out in Sect. 4, the probability that the best Schuster-Plummer f(ϵ)≥0 fit is better than the best NFW f(ϵ)≥0 fit is around 89%. The fact that Nube resides in a potential that is not cuspy is not fully unexpected. Its stellar mass, ∼3.9 × 10⁸ M_⊙, is large enough for the baryon feedback to modify the inner region of the global potential, turning cusps into cores (Sect. 1). The large size of core is surprising, though, a fact difficult to explain by the current cosmological CDM numerical simulations of low surface brightness galaxies (Montes et al. 2024). Another potential possibility to explain why Nube lives in a cored potential would be that the stars significantly contribute to the total mass of the system, so that the overall gravitational potential automatically follows the stellar distribution. However, this explanation is unlikely since the DM content estimated by Montes et al. (2024) is from 20 to 120 times larger than the stellar mass. We also studied the extent to which these conclusions depend on the estimated errors in the observed profile and the assumed mass-to-light ratio (Appendix E.2). Neither of these two issues jeopardizes the conclusion that cuspy profiles are disfavored.

The tool also allowed us to constrain the length scale of the potential (Eq. (30)). In terms of the core radius (i.e., when the density drops to half the central value), we find the cored potentials to be similar to the large stellar core shown by Nube (effective radius of 6.9 kpc). The possibility of constraining the relation between the core size of stars and DM happens to be one of the interesting outcomes of EIM (Sánchez Almeida et al. 2023).

The EIM-based tool described in the paper has considerable room for improvement. The use of other potentials to represent the DM distribution is as simple as computing the required characteristic densities (Eq. (4)). Section 5 outlines simple extensions that relax the need for isotropic velocities, so that the same kind of tool should work for systems with constant but anisotropic velocities and systems with gradients of anisotropy, from isotropic in their center to radially biased in the outskirts (also known as Osipkov-Merrit). We note that quasi-isotropic orbits and Osipkov-Merrit-like velocity anisotropies are indeed preferred by the model dwarf galaxies that formed in realistic cosmological numerical simulations (El-Badry et al. 2017; González-Samaniego et al. 2017; Orkney et al. 2023) and this is also found in dwarf spheroidal galaxies with observed kinematics (Massari et al. 2020; Kowalczyk & Łokas 2022). Moreover, Sánchez Almeida et al. (2024a) show how stellar cores are also inconsistent with NFW potentials in axi-symmetric systems using an extension of the original EIM. This extension represents a solid starting point to develop the tool further, so that we can constrain the gravitational potential dropping the spherical symmetry assumption.

Acknowledgments

Thanks are due to Andrés Asensio for guiding us on the use of the Bayesian tools. Thanks are also due to Ignacio Ferreras for insightful discussions on how to quantify the goodness of fits based on the two competing potentials. JSA acknowledges financial support from the Spanish Ministry of Science and Innovation, project PID2022-136598NB-C31 (ESTALLIDOS8), and from Gobierno de Canarias through EU FEDER funding, project PID2020010050. His visit to La Plata was partly covered by the MICINN through the Spanish State Research Agency, under Severo Ochoa Centers of Excellence Programme 2020-2023 (CEX2019- 000920-S). IT acknowledges support from the ACIISI, Consejería de Economía, Conocimiento y Empleo del Gobierno de Canarias and the European Regional Development Fund (ERDF) under grant with reference PROID2021010044 and from the State Research Agency (AEI-MCINN) of the Spanish Ministry of Science and Innovation under the grant PID2022-140869NB-I00 and IAC project P/302302, financed by the Ministry of Science and Innovation, through the State Budget and by the Canary Islands Department of Economy, Knowledge and Employment, through the Regional Budget of the Autonomous Community. MM acknowledges support from the project RYC2022-036949-I financed by the MICIU/AEI/10.13039/501100011033 and by Fondo Social Europeo Plus (FSE+). We acknowledge the use of the Python packages numpy (Harris et al. 2020), scipy (Virtanen et al. 2020), matplotlib (Hunter 2007), and emcee (Foreman-Mackey et al. 2013).

References

Appendix A: Effects on Nube of errors in determining the galaxy center

This appendix studies whether the inner drop in the density profile of Nube (Fig. 1) could be an artifact produced by selecting the wrong center when creating the surface density profile. To analyze the possibility, we carried out a series of Monte Carlo (MC) simulations like the one shown in Fig. A.1.

Fig. A.1. MC simulation to show the effect of using a wrong center to compute the stellar mass surface density radial profile of a galaxy. Top left panel: the dots represent individual stars. It shows only 103 of them to avoid overcrowding but the actual simulation has 104. The red circle is centered in the true distribution whereas the blue circle is offset by 4 kpc. Top right panel: same as the left panel color coded with the residual left when subtracting a radial profile computed with the 4 kpc offset (the thick orange line in the bottom panel). Bottom panel: surface density profiles (the symbols with error bars) computed using ∼50 rings offset from the true center as indicated in the top left panel. The error bars give the Poisson noise arising when counting stars. The true centered distribution is shown in blue.

Mock stars are randomly produced (the dots in the top left panel of Fig. A.1) following a surface density mimicking Nube (the blue solid line in the bottom panel of Fig. A.1). Nube is represented by a Schuster-Plummer profile (Eq. 1], where a = 2, b = 5, and c = 0) with a core radius of 8 kpc. Then the mock stars were counted in off-centered rings (e.g., the blue ring in Fig. A.1, top left panel) to produce density profiles like the red symbols with error bars in Fig. A.1, bottom panel. (The error bars account for the Poisson error when counting.) Note that the resulting density profile shows a drop in the innermost regions of the simulated profile (Fig. A.1, bottom panel, red symbols) which is not present in the original density profile (the blue solid line). A function like the one used to reproduce Nube (Eq. [1]), which allows for a variable inner slope, was used to fit the mock density profile (the orange thick line in Fig. A.1, bottom panel). The fitted profile has positive inner slope (c < 0, in the 3D profile).

Simulations like the ones shown in Fig. A.1 were repeated many times to produce the summary plot represented in Fig. A.2. Each point is the mean of a 100 different realizations of a mock galaxy with the same parameters (i.e, same profile, same offset, and same number of stars), with the error bars showing the standard deviation among all these realizations. Different colors represent different number of stars in the mock galaxy, which is a proxy for the error bars in the radial profile. The profile shown in the bottom panel of Fig. A.1 corresponds to MC realizations with 10⁴ stars per galaxy. The effective noise with this number of stars is similar to the one shown by the observations of Nube (Fig. 1) but, for the sake of comprehensiveness, we also show values for galaxies simulated with 10³ stars (the orange symbols; too noisy for Nube) and with 10⁵ stars (the green symbols; too good for Nube). Figure A.2 also includes a point representing the observation of Nube (the star symbol). We have been very generous with the error bars assigned to Nube. The formal error bar assigned to its center by the algorithm to compute the surface density profile (Sect. 2) is tiny (∼0.6 pix, which render ∼0.025 kpc considering a plate scale of 0.08 arcsec pix⁻¹ and a scale of 0.5 kpc arcsec⁻¹). Figure A.2 includes error bars corresponding to 1/10 of the core radius, which is very conservative. As far as the inner slope, we use the mean of the two values found when the last point of the observed profile is or not included in the fit (see Fig. 1 and Sect. 2.1), with the error bar being the semi-difference between the two.

Given the MC simulations described above, a number of arguments discard that the negative inner slope of Nube is caused by an error in the center used to compute the density profile. (1) The residuals of the fitting leave a dipole-like pattern (Fig. A.1, top right panel) which is not the residual observed in Nube (Fig. 11, right panel, in Montes et al. 2024). In other words, it is difficult to reconcile the residuals left by GALFIT on Nube with a significant global shift of the fitted function, in agreement with the negligible error that GALFIT provides (∼0.025 kpc). (2) Considering the MC simulation that represents Nube best (blue symbols in Fig. A.2) one needs an artificial offset of some 5 kpc to reproduce the observed slope, which is comparable with the core radius of Nube and, therefore, unrealistically large. (3) The results are robust in the sense that changing the hyper-parameters that define the MC simulation do not alter the conclusions (Appendix E.1). (4) Other bias may create it (Appendix E.2).

Fig. A.2. Summary of the MC to show the effect of using wrong centers to compute the stellar mass surface density radial profile of a galaxy. Inferred inner slope versus artificial offset (in kpc). The true inner slope is zero, as is the slope inferred from the fits when the offset goes to zero. The different colors represent different number of stars used to construct the mock galaxies (see the inset for the actual values). Those with errors closest to Nube are represented with blue symbols. Nube is inconsistent with them.

Appendix B: Characteristic densities for the potential arising from ρ_abc profiles

As we did in Sánchez Almeida et al. (2023), the case of a potential where the inner slope of the corresponding density is not zero (Schuster-Plummer) or minus one (NFW) can be treated in quite general terms using a ρ_abc profile as defined in Eq. (1). Using the Poisson equation for a spherically symmetric system (e.g., An & Zhao 2013), the potential is,

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_{\mathit{abc}}(r)=\frac{G{M}_{\mathit{abc}}(<r)}{r}+4\pi G{\displaystyle {\int }_r^{\infty }}t{\rho }_{\mathit{abc}}(t)dt,\end{array}$$(B.1)

with

$$\begin{array}{c}\hfill M_{\mathit{abc}}(<r)=4\pi {\displaystyle {\int }_0^r}{t}^2{\rho }_{\mathit{abc}}(t)dt,\end{array}$$

so that

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_{\mathit{abc}}(r)={ϵ}_{\mathit{max}}\frac{K(r/r_{\mathit{sp}},a,b,c)}{K(0,a,b,c)},\end{array}$$(B.2)

with

$$\begin{array}{c}\hfill {ϵ}_{\mathit{max}}=4\pi G{\rho }_{\mathit{sp}}{r}_{\mathit{sp}}^{2}K(0,a,b,c),\end{array}$$(B.3)

and

$$\begin{array}{c}\hfill K(x,a,b,c)=\frac{1}{x}{\displaystyle {\int }_0^x}\frac{t^{2-c}}{{(1+t^a)}^{(b-c)/a}}dt+{\displaystyle {\int }_x^{\infty }}\frac{t^{1-c}}{{(1+t^a)}^{(b-c)/a}}dt.\end{array}$$(B.4)

Then the characteristic density ξ_abc follows from Eq. (4). We use this approach to carry out the fits for ρ₂₃₀ potentials described in the main text (e.g., Fig. 9). In this case, a = 2, b = 3, and c = 0, which leads to a density with a core like a Schuster-Plummer profile and an outskirt similar to a NFW profile; see Fig. 10.

Appendix C: Classical interpretation of Eqs. (6) and (7)

The variable ϖ(ϵ), defined in Eq.(7) and appearing in Eq.(6), is interpreted in Sect. 3 as the mass corresponding to each relative energy ϵ. This variable, ϖ(ϵ), basically coincides with the quantity that in classical statistical mechanics is known as the density of states. It is usually denoted as g(E), where, in the context of stellar dynamics, E is the energy per unit mass (see Binney & Tremaine 2008, page 292, Eq. (4.56)). The density of states g(E) represents the phase-space volume per unit energy. That is, g(E)dE is the volume in phase-space corresponding to particles with energies in the range (E, E + dE). The quantity defined by us in Eq. (7) represents the density of states expressed in terms of the relative energy. In fact, our Fig. 4, which depicts ϖ(ϵ) for various potentials relevant for the present work, looks, qualitatively, as a mirror image of Fig. 4.3 of Binney & Tremaine (2008) because we plot the density of states against the relative energy, while Binney & Tremaine plot it against energy. The density of states plays a key role in classical statistical mechanics, and also in galactic dynamics, where it allows to compute the differential energy distribution, N(E) = g(E)f(E), defined in such a way that N(E)dE is the number of stars with energies in the range E + dE. Note that the integrand appearing in our Eq. (6) is basically the differential energy distribution expressed in terms of the relative energy ϵ. In line with the connection between ϖ(ϵ) and the density of states, the quantity defined in Eq.(4) would be related to the volume in phase-space per unit energy and per unit radius r. In fact, 4πr²ξ(ϵ, r)dϵdr is the volume in phase space corresponding to particles with energies in the range (ϵ, ϵ + dϵ) and radii in the range (r, r + dr). For our purpose, in order to implement the procedure for inferring the distribution function f(ϵ), we find it convenient to give the quantities ξ(ϵ, r) and ϖ(ϵ) an alternative interpretation, as we have already explained. Our interpretation, although different form the classical one, is consistent from a formal point of view and more useful in practice for our purpose.

Appendix D: Testing the numerical calculation of the eigendensities (proper densities)

Sánchez Almeida et al. (2023, Eq. A19) showed that the DF corresponding to a self-gravitating Schuster-Plummer density profile is analytic,

$$\begin{array}{c}\hfill f(ϵ)=\frac{120}{{(2\pi )}^{3/2}\mathrm{\Gamma }(9/2)}\frac{{\rho }_{\mathit{sp}}}{{ϵ}_{\mathit{max}}^{3/2}}{(ϵ/{ϵ}_{\mathit{max}})}^{7/2},\end{array}$$(D.1)

with ϵ_max and ρ_sp defined in Sect. 3.1 and Γ(x) the gamma function. Formally, it is like the polynomial expansion we use for f(ϵ) (Eq. [8]) with a single term,

$$\begin{array}{c}\hfill f(ϵ)=\frac{a_{\frac{7}{2}}}{{ϵ}_{\mathit{max}}^{3/2}}{(ϵ/{ϵ}_{\mathit{max}})}^{7/2},\end{array}$$(D.2)

where

$$\begin{array}{c}\hfill a_{\frac{7}{2}}={\rho }_{\mathit{sp}}\frac{120}{{(2\pi )}^{3/2}\mathrm{\Gamma }(9/2)}.\end{array}$$(D.3)

On the other hand, the surface density corresponding to a Schuster-Plummer density profile is also analytic (e.g., Binney & Tremaine 2008; Sánchez Almeida 2022),

$$\begin{array}{c}\hfill \mathrm{\Sigma }(R)={\rho }_{\mathit{sp}}{r}_{\mathit{sp}}\frac{4/3}{{[1+{(R/r_{\mathit{sp}})}^2]}^2},\end{array}$$(D.4)

which must be equal to the surface density provided by Eq. (11),

$$\begin{array}{c}\hfill \mathrm{\Sigma }(R)=a_{\frac{7}{2}}{S}_{\frac{7}{2}}(R),\end{array}$$(D.5)

which implies

$$\begin{array}{c}\hfill S_{\frac{7}{2}}(R)=\frac{{(2\pi )}^{3/2}\mathrm{\Gamma }(9/2)}{90}\frac{r_{\mathit{sp}}}{{[1+{(R/r_{\mathit{sp}})}^2]}^2}.\end{array}$$(D.6)

We have used the expression (D.6) to test the numerical algorithms developed to compute S_i(R). The result is shown in Fig. D.1, which compares the analytic and numerical expressions. The relative difference between them at each radius (|ΔS_i/S_i|; the red dashed line) is always smaller than 1 %, and it is typically smaller than 0.1 % in the core (R/r_sp < 1). In addition, the difference relative to the maximum of the profile is always smaller than 0.1 % (|ΔS_i|/maxS_i; the red dotted line).

Fig. D.1. Tests for the numerical calculations employed to compute the characteristic densities Si(R) in Eq. (12). The case of a self-gravitating Schuster-Plummer profile provides the analytic expression for i = 7/2 given in Eq. (D.6). The analytic and numerical expressions are shown as indicated in the inset (marked as i = 3.5). The differences between them (ΔSi/Si and ΔSi/maxSi) are also included as the red dotted and dashed lines. The cases i = 3 and i = 4 are shown for reference.

Appendix E: Sanity checks to test the algorithm that retrieves DFs

This Appendix collects a number of sanity checks that support the robustness and consistency of the diagnostic method used in the paper (Sect. 3.3).

Fig. E.1. Scatter plot summarizing how changing the hyper-parameters of the fit affect the interpretation of Nube’s surface density profile. The symbols with error bars represent the mean and the standard deviation of χ2 and ω inferred from the different posteriors. Each color corresponds to a different potential, with blue, orange, and green symbols representing NFW, Schuster-Plummer, and ρ230 potentials, respectively. Each type of symbol corresponds to a different set of hyper-parameters as indicated in the inset (see Appendixes E.1 and E.2 for details). The values for the nominal hyper-parameters used in the main text are portrayed as bullet symbols and denoted as “Reference” in the inset. The figure also includes the χ2 and ω from the best fits obtained with unconstrained f(ϵ), which are the same as the thick times symbols shown in Fig. 8. The vertical gray line is the same as that in Fig. 8.

E.1. Varying the hyper-parameters that define the algorithm

In order to test the dependence of the results on the hyper-parameters defining the algorithm, Nube’s profile was refitted changing them. Among others, these hyper-parameters define the priors of the Bayesian fits. Since the main result stemming from the application of the algorithm is the fact that the gravitational potential of Nube arises from a mass distribution with a core rather than a cusp (Sects. 4 and 6), we study whether this result is modified by the use of hyper-parameters different from the nominal ones in Sect. 3.3. The relationship between χ² and the innermost slope ω is used as diagnostics tool. These two quantities are used in the main text to argue that cuspy NFW potentials provide worst fits than the cored Plummer-Schuster or ρ₂₃₀ potentials. Simply put, χ² is larger for cuspy potentials that also provide ω farther away from the observed value (Sect. 4.2 and Fig. 8). Figure E.1 shows χ² versus ω for different hyper-parameters. The symbols with error bars represent the mean and the standard deviation of the corresponding distribution of χ² and ω inferred from the posterior. Each color corresponds to a different potential, with blue, orange, and green symbols representing NFW, Schuster-Plummer, and ρ₂₃₀ potentials, respectively. Each type of symbol corresponds to a different set of hyper-parameters. The nominal values used in the main text are portrayed as bullet symbols. These values are changed one at a time to produce other alternative hyper-parameters. The square symbols represent initializing the sampling of the posterior with the unconstrained best fitting least squares, which has f(ϵ) < 0. These fits also restrict the amplitudes defining the DF (a_i in Eq. [8]) to within 10⁻² and 10² of the best fit values. The times symbols correspond to doubling the number of wakers when sampling the posterior. The diamond symbols correspond to doubling the number of samples when sampling the posterior. The hexagon represents fits where the order of the polynomial used for f(ϵ) (n in Eq. [8]) differs from the nominal value (7 rather than 10). Other orders of the polynomial give similar results. For reference, Fig. E.1 includes the χ² and ω from the best fits obtained with unconstrained f(ϵ) which are the same as the thick times symbols shown in Fig. 8.

The main conclusion arising from inspecting Fig. E.1 is that NFW potential fits (in blue) are always worst than Schuster-Plummer potential fits (in orange) and ρ₂₃₀ fits (in green). Specifically, the NFW fits exhibit larger χ² and ω values deviating farther from those of Nube, represented in Fig. E.1 by the solid time symbols. This systematic preference for cored over cuspy fits is independent of the chosen set of hyper-parameters, reinforcing the robustness of the results against specific details of their selection.

E.2. Impact of the uncertainties in the mass profile of Nube

Fig. E.2. Similar to the fits corresponding to Figs. 6 and 7 but with smaller error bars that consider only photometric errors and errors in the mass-to-light ratio calibration. The fits assuming a Schuster-Plummer potential are shown in the top panels whereas those assuming a NFW potential are in the bottom panels. The result discussed in the main text remains. Compared with the NFW fits, the Schuster-Plummer potential fits with f(ϵ) > 0 (colored lines) have both an inner slope closer to the observed one and a smaller χ2. The color code is the same in the two panels and also the same of the code used in the figures of the main text. For further details, see the caption of Fig. 6.

Here we study two aspects of the uncertainties in the mass profile of Nube (Fig. 1) that may potentially have an impact on the conclusions. The first one has to do with the error bars employed in the calculation of χ² (Eq. [26]) which affect the posterior and may potentially influence the conclusions. The ones used in the main text and shown in Fig. 1 are larger than the scatter of observed points because, together with the Poisson errors associated with the photometry, they include the systematic errors associated to the sky subtraction and the estimate of mass-to-light ratio (Montes et al. 2024). We also study what happens if these other errors are disregarded leaving error bars closer to the scatter of the individual points in the radial profile. The result of the exercise is shown in Fig. E.2. Obviously, the overall value of the χ² increases, but the conclusions reached in the main text remain. Compared with the NFW fits, Schuster-Plummer potential fits with f(ϵ) > 0 (the colored lines in Fig. E.2) have an inner slope closer to the observed one together with a significantly smaller χ². This fact can be better appreciated in the diagnostic plot shown in Fig. E.1, where the corresponding χ² and ω are shown as star symbols labeled with "No Sky Err."

Fig. E.3. Equivalent to Figs. 6 and 7 but using a scaled version of the r-band surface band profile of Nube as mass surface density. This is equivalent to assuming a constant mass-to-light ratio throughout the galaxy. Note that the central drop in mass has gone away (cf. Fig. 1 with the data points in figures a). The fits assuming a Schuster-Plummer potential are shown in the top panels whereas those assuming a NFW potential are in the bottom panels. This other alternative calibration of the mass profile does not modify the conclusions based on the original calibration. The color code is the same in the two panels and also the same of the code used in the figures of the main text. For further details, see the caption of Fig. 6. Note that the EIM technique is insensitive to a global scaling of the mass profile, and this arbitrary factor has been chosen here so that the level of Σ(R) is similar to the original one in Fig. 1.

The second study refers to the impact of the mass-to-light ratio calibration. As it is mentioned in Sect. 4.1, the central drop in mass in Fig. 1 is not present in the light profile of Nube and may be an artifact appearing when transforming the observed photometry into stellar mass. We consider the impact of the used mass-to-light ratio on the conclusion by using a constant value rather than the color-varying ratio employed by Montes et al. (2024) and used in our study. In this test we use,

$$\begin{array}{c}\hfill \mathrm{\Sigma }(R)\propto {10}^{-0.4\mathrm{SB}},\end{array}$$(E.1)

with SB the observed surface brightness. A profile thus compute is shown in Fig. E.3, where we have chosen the r-band photometry because the dependence of the mass-to-light ratio on colors is smaller in the red, yet the r-band exhibits low noise. The chosen errors are somewhat arbitrary trying to account for photometric errors and matching those in Fig. E.3. This mass profile do not show the central drop of Fig. 1. Using this data, we repeat the analysis and the resulting fits and DFs are shown in Fig. E.3. As for the original data set, Schuster-Plummer potential fits are better than the NFW potential fits. They have an inner slope closer to the observed one and their χ² is significantly smaller. Their values are included in Fig. E.1 using the symbol O and labeled as 10^−0.4 SB. In addition, the Schuster-Plummer potential best fit does not require the negative distribution function needed for the original Nube profile: compare the red thick line in Fig. E.3b (top panel) with Fig. 6b.

Continuing with the impact of changing the mass-to-light ratio, we construct another mock profile using the Nube data but rearranging the order of the observed Σ(R) so that the profile monotonically decreases outward in the inner part. Note that such rearrangements leaves a profile consistent with the original data set keeping in mind the large error bars (Fig. 1). The resulting profile is shown in Fig. E.4. This figure is similar to Fig. E.5, and evidences how the Schuster-Plummer potential does an excellent job with f(ϵ)≥0 whereas the profiles forced to have f(ϵ) > 0 in NFW potentials provide much worst fits.

Fig. E.4. Mock profiles using all the data from Nube but rearranging the inner points so that the resulting profile decreases monotonically outward (the symbols labeled as Mock Obs). This re-arrangement is consistent with the error bar of the observation but clears out the problem of the algorithm we employ to explain positive inner slopes. Top panels: case of a Schuster-Plummer potential. Bottom panels: case of a NFW potential. This figure is similar to Figs. E.5, with an identical color code.

E.3. Recovering properties of known core profiles

Fig. E.5. Summary of the analysis carried out with a mock mass surface density profile corresponding to a Schuster-Plummer profile (the symbols labeled as Mock Obs). The figures are similar to Figs. 6, 7, and 9 in the main text. The result assuming a Schuster-Plummer potential is shown in the top panels whereas the result assuming a NFW potential is in the bottom panels. Note how the free-f(ϵ) fits (the thick red lines) yield f(ϵ) > 0 for the Schuster-Plummer potential whereas it requires f(ϵ) < 0 for the NFW potential. The fits forced to have f(ϵ) > 0 (color thin lines) are much worst in the case of the NFW profile. The color code is the same in the two panels and also the same of the code used in the figures of the main text.

We repeat the analysis with a mock mass surface density profile corresponding to a Schuster-Plummer profile with error similar to those of Nube. The result assuming a Schuster-Plummer potential is shown in Fig. E.5, top panels. Obviously, a Schuster-Plummer profile is fully consistent with a Schuster-Plummer potential and, as expected, the inferred f(ϵ) (the red solid line) is always positive and in agreement with its expected form (Sánchez Almeida et al. 2023, given explicitly in Eq. D.1 above). Moreover, the surface density profiles and DFs resulting from the MCMC exploration of the posterior forcing f(ϵ)≥0 are also compatible with this best fit. This test shows what is to be expected in case of assuming a potential fully consistent with the observed Σ(R). On the other hand, the same mock surface density is analyzed assuming a NFW potential (Fig. E.5, bottom panels). The required f(ϵ) is negative, and so unphysical, and the fits forcing f ≥ 0 are way off the best fitting profile.

All Figures

Fig. 1. Stellar mass surface density profile of Nube as observed by Montes et al. (2024, the symbols with error bars). The figure includes three different fits to this observation: a projected polytrope (PP; the black line, with its index m given in the inset) and two projected abc profiles (Eq. (1)). The inner slope (−c) is forced to be zero in the PP fit whereas becomes positive (i.e., c < 0) when allowed to vary (see the inset). The gray error bars are those provided by Montes et al. (symmetric in a linear scale and so asymmetric in the logarithmic representation used in the figure) whereas the symmetrized ones1 (the blue bars) are equivalent but symmetric in the logarithmic representation.

In the text

	Fig. 2. Characteristic density profiles corresponding to each relative energy α = ϵ/ϵmax in a Schuster-Plummer potential (the solid lines; Eq. (16)) and in a NFW potential (the dashed lines; Eq. (23)). The symbol β stands for the normalized radial coordinate r/rsp. All energies contribute to the innermost regions (β ≪ 1) whereas only the smallest energies contribute to the outer halo (β ≫ 1).
In the text

	Fig. 3. Similar to Fig. 2 but representing the contribution to the total density of the different energies (α) at each constant radial position (β).
In the text

Fig. 4. Specific mass corresponding to each energy (Eq. (7) with α = ϵ/ϵmax) for Schuster-Plummer and NFW potentials. Top panel: log-linear representation. Bottom panel: log-log representation. Note that the specific mass diverges for low energies as a power-law of exponent −2.5, for the Schuster-Plummer potential, and around −2.81, for the NFW potential (the dashed lines labeled in the inset of the bottom panel).

In the text

	Fig. 5. Characteristic functions corresponding to the polynomial expansion of the DF: see Eqs. (8), (11), and (12). The symbol i stands por the exponent of the monomial, and the solid and dashed lines correspond to the Schuster-Plummer and the NFW potentials, respectively. It begins at i = 3 since for i ≲ 3 the total mass of the resulting density profile diverges (see main text).
In the text

Fig. 6. Diagnostic plot for the technique. (a) Fit to the stellar mass surface density observed in Nube (the blue symbols) assumed to reside in a Schuster-Plummer potential. The red thick solid line represents the least-squares best fit with an unconstrained DF f(ϵ), whereas the other thin lines are the fits derived from the MCMC exploration of the posterior forcing f(ϵ) to be ≥0∀ϵ and beginning the exploration from the f(ϵ)≥0 least squares solution. These other fits are color coded according to the innermost slope of the surface density profile (ω in Eq. (27)), which always happens to be negative but close to zero. (b) DFs required by the best fit and by the fits forced to have f(ϵ)≥0. Note how the best fit requires f(ϵ) < 0 and so is unphysical. The color code is the same as in (a).

In the text

	Fig. 7. Similar to Fig. 6 except that Nube is assumed to reside in a NFW potential. The best fitting function (the red thick solid line) is similar to the best fit assuming a Schuster-Plummer potential (Fig. 6a), but the physically realizable f(ϵ)≥0 fits are clearly worst and have inner slopes differing from zero significantly (see the color code of the thin lines).
In the text

Fig. 8. (a) Scatter plot χ2 versus innermost slope ω of the fits for the three potentials resulting from the MCMC navigation of the posterior. The large times symbols represent the best fit obtained with unconstrained f(ϵ), which yield unphysical f(ϵ) < 0. The vertical gray line marks ω = −0.01. (b) Histograms of the distribution of innermost slopes for the points in (a). (c) Histograms of χ2 for the points in (a). The arrows represent the χ2 of the best fit obtained with unconstrained f(ϵ). Panels a, b, and c share the same color code as indicated in the insets.

In the text

	Fig. 9. Similar to Fig. 6 except that Nube is assumed to reside in a ρ230 potential, which has a core and approaches a NFW profile in the outskirts. The best fitting function (the red thick solid line) is similar to the best fit in a NFW potential (Fig. 7), but the physically sensible f(ϵ)≥0 are clearly better and have a inner slope closer to zero.
In the text

Fig. 10. Comparison between the stellar surface density profile of Nube and the best-fitting potentials with f ≥ 0. The plot shows the mass surface density that gives rise to the best-fitting Schuster-Plummer potential (blue line), NFW potential (green line), and ρ230 potential (magenta line). The spatial scaling of the potentials is set by the fit whereas the vertical scaling is arbitrary, and it was arbitrarily chosen to be 30 times the stellar mass of Nube.

In the text

Fig. A.1. MC simulation to show the effect of using a wrong center to compute the stellar mass surface density radial profile of a galaxy. Top left panel: the dots represent individual stars. It shows only 103 of them to avoid overcrowding but the actual simulation has 104. The red circle is centered in the true distribution whereas the blue circle is offset by 4 kpc. Top right panel: same as the left panel color coded with the residual left when subtracting a radial profile computed with the 4 kpc offset (the thick orange line in the bottom panel). Bottom panel: surface density profiles (the symbols with error bars) computed using ∼50 rings offset from the true center as indicated in the top left panel. The error bars give the Poisson noise arising when counting stars. The true centered distribution is shown in blue.

In the text

Fig. A.2. Summary of the MC to show the effect of using wrong centers to compute the stellar mass surface density radial profile of a galaxy. Inferred inner slope versus artificial offset (in kpc). The true inner slope is zero, as is the slope inferred from the fits when the offset goes to zero. The different colors represent different number of stars used to construct the mock galaxies (see the inset for the actual values). Those with errors closest to Nube are represented with blue symbols. Nube is inconsistent with them.

In the text

Fig. D.1. Tests for the numerical calculations employed to compute the characteristic densities Si(R) in Eq. (12). The case of a self-gravitating Schuster-Plummer profile provides the analytic expression for i = 7/2 given in Eq. (D.6). The analytic and numerical expressions are shown as indicated in the inset (marked as i = 3.5). The differences between them (ΔSi/Si and ΔSi/maxSi) are also included as the red dotted and dashed lines. The cases i = 3 and i = 4 are shown for reference.

In the text

Fig. E.1. Scatter plot summarizing how changing the hyper-parameters of the fit affect the interpretation of Nube’s surface density profile. The symbols with error bars represent the mean and the standard deviation of χ2 and ω inferred from the different posteriors. Each color corresponds to a different potential, with blue, orange, and green symbols representing NFW, Schuster-Plummer, and ρ230 potentials, respectively. Each type of symbol corresponds to a different set of hyper-parameters as indicated in the inset (see Appendixes E.1 and E.2 for details). The values for the nominal hyper-parameters used in the main text are portrayed as bullet symbols and denoted as “Reference” in the inset. The figure also includes the χ2 and ω from the best fits obtained with unconstrained f(ϵ), which are the same as the thick times symbols shown in Fig. 8. The vertical gray line is the same as that in Fig. 8.

In the text

Fig. E.2. Similar to the fits corresponding to Figs. 6 and 7 but with smaller error bars that consider only photometric errors and errors in the mass-to-light ratio calibration. The fits assuming a Schuster-Plummer potential are shown in the top panels whereas those assuming a NFW potential are in the bottom panels. The result discussed in the main text remains. Compared with the NFW fits, the Schuster-Plummer potential fits with f(ϵ) > 0 (colored lines) have both an inner slope closer to the observed one and a smaller χ2. The color code is the same in the two panels and also the same of the code used in the figures of the main text. For further details, see the caption of Fig. 6.

In the text

Fig. E.3. Equivalent to Figs. 6 and 7 but using a scaled version of the r-band surface band profile of Nube as mass surface density. This is equivalent to assuming a constant mass-to-light ratio throughout the galaxy. Note that the central drop in mass has gone away (cf. Fig. 1 with the data points in figures a). The fits assuming a Schuster-Plummer potential are shown in the top panels whereas those assuming a NFW potential are in the bottom panels. This other alternative calibration of the mass profile does not modify the conclusions based on the original calibration. The color code is the same in the two panels and also the same of the code used in the figures of the main text. For further details, see the caption of Fig. 6. Note that the EIM technique is insensitive to a global scaling of the mass profile, and this arbitrary factor has been chosen here so that the level of Σ(R) is similar to the original one in Fig. 1.

In the text

Fig. E.4. Mock profiles using all the data from Nube but rearranging the inner points so that the resulting profile decreases monotonically outward (the symbols labeled as Mock Obs). This re-arrangement is consistent with the error bar of the observation but clears out the problem of the algorithm we employ to explain positive inner slopes. Top panels: case of a Schuster-Plummer potential. Bottom panels: case of a NFW potential. This figure is similar to Figs. E.5, with an identical color code.

In the text

Fig. E.5. Summary of the analysis carried out with a mock mass surface density profile corresponding to a Schuster-Plummer profile (the symbols labeled as Mock Obs). The figures are similar to Figs. 6, 7, and 9 in the main text. The result assuming a Schuster-Plummer potential is shown in the top panels whereas the result assuming a NFW potential is in the bottom panels. Note how the free-f(ϵ) fits (the thick red lines) yield f(ϵ) > 0 for the Schuster-Plummer potential whereas it requires f(ϵ) < 0 for the NFW potential. The fits forced to have f(ϵ) > 0 (color thin lines) are much worst in the case of the NFW profile. The color code is the same in the two panels and also the same of the code used in the figures of the main text.

In the text