11institutetext: Instituto de Astrofísica de Canarias, c/ Vía Láctea s/n,
E-38205, La Laguna, Tenerife, Spain 22institutetext: Departamento de
Astrofísica, Universidad de La Laguna, E-38203, La Laguna, Tenerife, Spain
33institutetext: CeBio y Departamento de Ciencias Básicas, Universidad
Nacional del Noroeste de la Prov. de Buenos Aires, UNNOBA, CONICET, Roque
Saenz Peña 456, Junin, Argentina

# Application of the Eddington inversion method to constrain the dark matter
halo of galaxies using only observed surface brightness profiles

Jorge Sánchez Almeida 1122 Angel R. Plastino 33 Ignacio Trujillo 1122

(Received ; accepted March, 14, 1592)

###### Abstract

Context. The halos of low-mass galaxies may allow us to constrain the nature
of dark matter (DM), but the kinematic measurements to diagnose the required
properties are technically extremely challenging. However, the photometry of
these systems is doable.

Aims. Using only stellar photometry, constrain key properties of the DM haloes
in low-mass galaxies.

Methods. Unphysical pairs of DM gravitational potentials and starlight
distributions can be identified if the pair requires a distribution function
$f$ that is negative somewhere in the phase space. We use the classical
Eddington inversion method (EIM) to compute $f$ for a battery of DM
gravitational potentials and $\sim 100$ observed low-mass galaxies with
$M_{\star}$ between $10^{6}$ and $10^{8}\,{\rm M_{\odot}}$. The battery
includes NFW potentials (expected from cold DM) and potentials stemming from
cored mass distributions (expected in many alternatives to cold DM). The
method requires fitting the observed profiles with analytic functions, for
which we use polytropes (with zero inner slope, a.k.a. core) and profiles with
variable inner and outer slopes.

Results. In general, the polytropes fit well the observed starlight profiles.
If they were the correct fits (which could be the case) then all galaxies are
inconsistent with NFW-like potentials. Alternatively, when the inner slope is
allowed to vary for fitting, between 40 % and 70 % of the galaxies are
consistent with cores in the stellar mass distribution and thus inconsistent
with NFW-like potentials.

Conclusions. Even though the stellar mass of the observed galaxies is still
not low enough to constrain the nature of DM, this work shows the practical
feasibility of the EIM technique to infer DM properties only from photometry.

###### Key Words.:

dark matter — Galaxies: dwarf — Galaxies: fundamental parameters — Galaxies:
halos — Galaxies: photometry — Galaxies: structure

##  1 Introduction

The classical Eddington Inversion Method (EIM; Eddington, 1916; An & Evans,
2006; Binney & Tremaine, 2008; Lacroix et al., 2018; Ciotti, 2021) can be used
to constrain the shape of the dark matter (DM) halo of a galaxy using only the
observed stellar mass distribution. It provides the distribution function (DF)
in the 6D phase space to be followed by a given mass density profile immerse
in a given gravitational potential. The minimal requirement for any feasible
DF is to be non-negative everywhere in the phase space, and this seemingly
simple constrains is enough to discard as unphysical many combinations of
density and potential leading to a DF $<0$. Recently, Sánchez Almeida et al.
(2023) use the EIM in the regime of low mass galaxies, where the contribution
of the baryons to the overall potential can be neglected and the potential is
set only by the DM. Among others, they showed that cored stellar mass
distributions (i.e., with constant density when the distance to the center
goes to zero) are inconsistent with a Navarro, Frenk, and White potential,
expected in DM dominated systems with collision-less cold DM (the so-called
NFW profiles after Navarro et al., 1997). Strictly speaking, this result holds
for spherically symmetric particle systems with an isotropic velocity
distribution, however, extensions that relax these assumptions and yield
similar constraints exist (a detailed account is given in Sect. 2.1).

The technique has great potential since, using only photometry, it may be able
to constrain the shape of the DM halos in ultra-low mass galaxies. In this
mass regime the baryon feedback is unable to turn cuspy DM halos into cored
halos. Thus, if the DM haloes of these galaxies happen to have a core, it will
indicate a departure of the DM from being collision-less, reflecting the much
sought-after presently-unknown true nature of DM (whether it is fuzzy, self-
interacting, warm, or other alternatives; e.g., Dodelson & Widrow, 1994; Hu et
al., 2000; Spergel & Steinhardt, 2000; Bechtol et al., 2022; Carr et al.,
2023). The limit defining ultra-low mass galaxies, i.e., the largest stellar
mass unable to modify the DM profile, is model dependent (e.g., Read et al.,
2016), however, it roughly corresponds to a stellar mass of
$M_{\star}<10^{6}\,{\rm M}_{\odot}$ or a DM halo mass of $M_{h}<10^{10}\,{\rm
M_{\odot}}$ (e.g., Di Cintio et al., 2014b; Chan et al., 2015; Hayashi et al.,
2020; Jackson et al., 2021; Expósito-Márquez et al., 2023). Traditionally,
halo shapes are deduced from kinematical measurements, which require high
spectral resolution spectroscopy and which, keeping in mind the need for large
statistics to derive firm conclusions on the nature of DM, are virtually
imposible in this ultra-low mass regime. However, the broad-band photometry
needed to infer their stellar mass profiles starts to be doable (e.g.,
Trujillo et al., 2021) and will become routinely simple in the near future
with instruments like the Rubin Observatory (Ivezić et al., 2019). Here is
where the EIM comes into play, constraining the properties of the halo mass
distribution from the stellar light alone.

The present work represents the first application to real galaxies the EIM-
based technique devised by Sánchez Almeida et al. (2023). We put forward
diagnostic diagrams and then discuss the kind of constraints to be expected in
practice. The target galaxies were collected by Carlsten et al. (2021, 2022)
who select them as dwarf satellites of Milky Way-like galaxies. This sample is
not ideal for ascertaining the nature of DM since (1) their masses exceed the
$10^{6}\,{\rm M}_{\odot}$ limit and (2) the fact that they are satellites
complicates the interpretation of the DM halo shape in terms of the DM nature
alone. The evolution of satellites may be influenced by the central galaxy
through tidal forces and gas pressure gradients (e.g., Lan et al., 2016;
Sánchez Almeida et al., 2017). However, this galaxy set is a good testbed
since it comprises a large number of objects ($\gtrsim$ 100) with density
profiles comparable in lack of symmetry and noise with those to be expected
soon for the appropriate galaxies.

The paper is structured as follows: Sect. 2 summarizes the technique to be
employed, with Sect. 2.1 giving an overview of the EIM, including an extension
that goes beyond the spherical symmetry assumption (Appendix B). The actual
algorithm is described in Sect. 2.2 and tested in Appendix C. The data are
presented in Sect. 3, with the new diagnostic diagrams spelled out in Sect. 4.
These diagrams depend on the functions used for fitting the observed profiles.
For this task, we firstly use projected polytropes in Sect. 5.1. Polytropes
have cores, an assumption relaxed in Sect. 5.2 (free inner slope fits) and
Sect. 5.3 (free inner and outer slope fits). The conclusions arising from the
results are analyzed in Sect. 6 including the impact of assumptions like
spherical symmetry and isotropic velocity distributions made when we applied
the EIM.

##  2 The Eddington inversion method (EIM)

For the sake of comprehensiveness, and to set up the scene, Sect. 2.1 gives a
brief account of the EIM and the implications for low mass galaxies found by
Sánchez Almeida et al. (2023). Then the practical application of the method
carried out in the present paper is explained in Sect. 2.2.

###  2.1 Contextual briefing of the EIM

For spherically symmetric systems of particles with isotropic velocity
distribution, the phase-space DF $f(\epsilon)$ depends only on the particle
energy $\epsilon$. The space density $\rho(r)$ is an integral over all
velocities $v$ at a fixed distance $r$. If $\Psi(r)=\Phi_{0}-\Phi(r)$ is the
relative potential energy, then the relative energy (per unit mass) is
$\epsilon=\Psi-\frac{1}{2}v^{2}$, so that the density becomes (e.g., Binney &
Tremaine, 2008, Sect. 4.3),

$\rho(r)=4\pi\sqrt{2}\,\int_{0}^{\Psi(r)}\,f(\epsilon)\sqrt{\Psi(r)-\epsilon}\,d\epsilon,$
(1)

with $\Phi(r)$ the gravitational potential and $\Phi_{0}$ its value at the
edge of the system. For realistic systems, the relative potential $\Psi$ is a
monotonically decreasing function of the distance from the center $r$.
Consequently, $\rho$ can be regarded as a function of $\Psi$. Differentiating
$\rho$ with respect to $\Psi$ yields,

$\frac{d\rho}{d\Psi}\,=\,2\pi\sqrt{2}\,\int_{0}^{\Psi}\,\frac{f(\epsilon)}{\sqrt{\Psi-\epsilon}}\,d\epsilon.$
(2)

Inverting this Abel integral and after some mathematical manipulations, one
finds the equation for $f(\epsilon)$ in terms of the spatial density,

$f(\epsilon)=\frac{1}{\sqrt{2}\pi^{2}}\,\int_{0}^{\epsilon}\,\frac{d^{3}\rho}{d\Psi^{3}}\,\sqrt{\epsilon-\Psi}\,d\Psi,$
(3)

which provides the renowned EIM equation. It gives the DF $f(\epsilon)$
consistent with $\rho$ and $\Phi$.

However, given an arbitrary pair of $\rho$ and $\Phi$ (or $\rho$ and $\Psi$),
it is not guaranteed that the DF yielded by Eq. (3) is positive everywhere in
the phase space. If $f<0$ somewhere it means that this particular combination
of $\rho$ and $\Phi$ cannot exist in practice. Sánchez Almeida et al. (2023)
analyze a sizable set of pairs $\rho$ – $\Phi$ showing a number of
inconsistencies that may be important for real galaxies. The simplest one, yet
very powerful, directly follows from Eq. (2). If $d\rho/d\Psi=0$ then the
integral in the right hand side of Eq. (2) has to be zero so that
$f(\epsilon)<0$ somewhere within the interval $0\leq\epsilon\leq\Psi$. This is
precisely what happens with any cored density profile immersed in a NFW
potential. Since the only variable is $r$,

$\frac{d\rho}{d\Psi}\,=\,\frac{d\rho/dr}{d\Psi/dr}.$ (4)

A cored mass density is defined to have

$\lim_{r\to 0}\frac{d\rho}{dr}=0,$ (5)

which is inconsistent with a NFW background potential, which has

$\lim_{r\to 0}\frac{d\Psi}{dr}\neq 0.$ (6)

The above conclusion holds for spherically symmetric systems with isotropic
velocity distribution, however, this restriction can be relaxed to accomodate
more general systems. Surpassing the simplest assumption requires treating
particular cases individually, each one with its own peculiarities, therefore,
only a handful of them have been analyzed so far, probably representing the
tip of the iceberg. In order to discuss these individual cases, it is
necessary to define the velocity field anisotropy parameter $\beta$,

$\beta\,=\,1\,-\frac{\sigma_{vt}^{2}}{2\sigma_{vr}^{2}},$ (7)

where $\sigma_{vr}$ and $\sigma_{vt}$ are the radial and tangential velocity
dispersion, respectively. The isotropic systems leading to the EIM Eq. (3)
have $\sigma_{vt}^{2}=2\sigma_{vr}^{2}$ and so $\beta=0$. Systems with
radially biased orbits have $\beta>0$ whereas systems with tangentially biased
orbits have $\beta<0$. Systems with circular orbits are a limiting case of the
latter having $\sigma_{vr}=0$ and therefore having $\beta=-\infty$. The case
where

$\beta(r)=\frac{r^{2}}{r^{2}+r_{b}^{2}},$ (8)

is called Osipkov-Merritt model and introduces a new characteristic length
scale $r_{b}$. It is tractable analytically and approximately describes the
global trend expected in low mass galaxies, with $\beta\sim 0$ in the center
and then increasing outwards ($\beta>0$).

The constraints on the potential imposed by the stellar mass distribution are
best described in terms of the inner logarithmic slope of the stellar density
profile $c$,

$\lim_{r\to 0}\frac{d\log\rho}{d\log r}=-c,$ (9)

and the inner slope of the total mass density profile $\rho_{p}$ defining the
gravitational potential,

$\lim_{r\to 0}\frac{d\log\rho_{p}}{d\log r}=-c_{p}.$ (10)

The univocal association between $\Psi$ and $\rho_{p}$ is granted through the
Poisson’s equation which, for spherically symmetric systems, yields (e.g., An
& Zhao, 2013),

$\Psi=4\pi
G\,\left[\frac{1}{r}\,\int_{0}^{r}\,t^{2}\,\rho_{p}(t)\,dt+\int_{r}^{\infty}\,t\,\rho_{p}(t)\,dt\right].$
(11)

In terms of $c$ and $c_{p}$, a cored stellar density profile has $c=0$ whereas
a NFW potential has $c_{p}=1$. From now on, the term soft-core is used to
denote those profiles where the inner slope is not exactly zero but close to
it ($c$ or $c_{p}\gtrsim 0$). Equipped with these definitions, the main known
constraints imposed by the EIM are (Sánchez Almeida et al., 2023):

  * •

Stellar cores ($c=0$) and NFW potentials are incompatible, provided the
velocity distribution is isotropic ($\beta=0$). This was already proven above.

  * •

Stellar cores are also incompatible with potentials stemming from a total
density with a quasi-core ($c_{p}>0$) for $\beta=0$.

  * •

As expected for physical consistency, stellar cores ($c=0$) and potentials
resulting from cored density profiles ($c_{p}=0$) are consistent in isotropic
($\beta=0$) and radially biased systems ($\beta>0$).

  * •

Stellar cores and NFW potentials are also incompatible in systems with
anisotropic velocities following the Osipkov-Merritt model.

  * •

Stellar cores and NFW potentials are incompatible in systems with radially
biased orbits (${\rm constant\leavevmode\nobreak\ }\beta>0$). Actually, a
stellar core is incompatible with any potential without a core in systems with
constant radially biased orbits ($\beta>0$).

  * •

Circular orbits ($\beta=-\infty$) can accommodate any combination of baryon
density and potential, including a cored stellar density in a NFW potential.
However, this configuration is very artificial from a physical standpoint.

  * •

A convex linear combination of two DFs that are both consistent with a given
potential constitutes a new DF which is also consistent with the same
potential. Thus, one may think that adding a positive DF for circular orbits
may compensate a negative DF for isotropic orbits to yield a positive
physically sensible DF. However, this is not the case. Independently of the
relative weight, the superposition of an unphysical DF for isotropic
velocities ($\beta=0$) with a physically realizable DF for circular orbits
($\beta=-\infty$) always yields unphysical DFs.

  * •

Soft-cores are inconsistent with NFW potentials when $c\lesssim 0.1$ while
they are consistent when $c\gtrsim 0.1$.

  * •

The inner slope of a soft stellar core and the radial anisotropy are related
so that $c>2\beta$ for the system to be consistent. In other words, systems
with radially biased orbits are inconsistent with soft stellar cores.

  * •

Positive inner slope in the stellar distribution ($c<0$), where the density
grows outwards, is discarded in every way for isotropic velocity.

  * •

For stellar densities and potentials with the same shape (whether cored or
not), the stellar density distribution cannot be broader than twice the width
of the density characterizing the potential. This result refers to isotropic
velocities and may be used in real galaxies to set a lower limit to the size
of the DM halo from the size of the observed starlight.

  * •

The above results do not depend on a global factor scaling the stellar density
profile, therefore, they do not depend on the (unknown) ratio between the
stellar mass and the total mass of the system. Moreover, it implies that light
profiles or star number counts could replace mass density in the EIM.

To the above list of known constraints, here we add another one relaxing the l
symmetry assumption. In Appendix B, we prove that the incompatibility between
NFW potentials and stellar cores also holds for axi-symmetric systems of
arbitrary axial ratio. It is particularly relevant because this fact shows how
the incompatibility goes beyond the spherical symmetry assumption, being more
fundamental than, and not attached to, this particular hypothesis.

###  2.2 Procedure followed in this paper

As summarized in Sect. 2.1, the EIM gives the DF $f$ needed for the stellar
density profile $\rho$ to reside in a gravitational potential $\Phi$. If two
arbitrary $\rho$ and $\Phi$ yield $f<0$ somewhere in the phase space, this
particular combination is physically inconsistent and can be discarded. Here
we apply the technique to constrain the properties of the haloes in a number
of observed galaxies. The procedure works as follows:

  1. 1.

Measure $\Sigma_{o}$ in a real galaxy, i.e., the projection along the line-of-
sight (LOS) of the density $\rho_{o}$.

  2. 2.

Fit $\Sigma_{o}$ with synthetic $\Sigma$ to estimate $\rho_{o}$, i.e., the
$\rho$ whose $\Sigma$ best-fits $\Sigma_{o}$. For spherically symmetric
systems, the relation between $\rho$ and $\Sigma$ is univocal (e.g., Binney &
Tremaine, 2008),

$\Sigma(R)=2\,\int_{R}^{\infty}\frac{r\,\rho(r)}{\sqrt{R^{2}-r^{2}}}\,dr,$
(12)

with $R$ the projected distance to center of the system.

  3. 3.

Using the EIM (Eq. [3]), compute the $f_{oi}$ corresponding to the $\rho$ best
fitting $\rho_{o}$ and a battery of $\Phi_{i}$, with $i=1,2,\dots$

  4. 4.

Discards as unphysical, all the $\Phi_{i}$ that yield $f_{oi}<0$ somewhere in
the phase space.

The above algorithm can be applied to individual galaxies as well as to sets
of them ($o=1,2,\dots$). For example, one can compute the fraction of times
that a fixed $\Phi_{i}$ is unphysical considering all the galaxies in a
particular set. This is the approach followed in our work, which requires the
computation of $f_{oi}$ for all galaxies in the set
($\forall\leavevmode\nobreak\ o$) and for all potentials
($\forall\leavevmode\nobreak\ i$). Then we use as diagnostic for the
reliability of the potential $\Phi_{i}$ the fraction of times that it is
consistent, i.e., the fraction of times where $f_{oi}\geq 0$ everywhere.
Practical examples of the procedure are given in Sect. 4.

This scheme has been tested for consistency using the set of Globular Clusters
(GCs) collected by Miocchi et al. (2013). In general, GCs are known to have
little or no DM (e.g., Sollima et al., 2012; Ibata et al., 2013), therefore,
GC can be regarded as self-gravitating systems so that the gravitational
potential is set directly by the observed stellar distribution. Moreover, the
stellar surface density of GCs usually show cores. Thus, GCs are excellent
testbeds for our EIM-based approach since it provides constraints on the
properties of the potential which can be compared directly with the true
values obtained from the surface density profile. The tests are detailed in
Appendix C, and they lead to the conclusion that our scheme gives consistent
results with the gravitational potential of the GCs inferred from their
stellar distribution.

Figure 1:  Surface brightness profiles in the $g$ band for the early-type
galaxies in Carlsten et al. (2021, 2022). The left panel shows stellar masses
between $10^{6}$ and $10^{7}\,{\rm M_{\odot}}$ whereas the right panel
corresponds to masses between $10^{7}$ and $10^{8}\,{\rm M_{\odot}}$. Note how
the profiles tend to present a central plateau or core. The color code is
maintained in Figs. 6, 10, and 14.

Figure 2: Examples of projected polytrope fits to some of the profiles in Fig.
1. The first three rows contain god fits. The bottom row shows examples of bad
fits, usually profiles with bumps resulting from deviations from axial
symmetry of the original galaxies. The blue dots represent the observed points
whereas the orange line shows the fit, which is solid or dashed depending on
whether it is within or outside the range of fitted radii. The secondary sub-
panel in each panel shows the residuals of the fit. The insets provide
information on the galaxy and the fit.

##  3 Observed profiles and fits

This section describes the practical realization of steps number 1 and 2 in
the procedure described in Sect. 2.2.

The observed surface density mass distributions analyzed in the paper come
from Carlsten et al. (2021, 2022). They provide the surface brightness profile
in the $g$ band ($\mu_{g}$) for a set of dwarf satellites of Milky Way-like
galaxies. They use photometry from various existing surveys including the DESI
Legacy Imaging Surveys (Dey et al., 2019) and the Beijing–Arizona Sky Survey
(Zou et al., 2017, 2018); see Carlsten et al. (2021) for details. The galaxies
are low mass dwarfs with stellar masses between $10^{6}$ and $10^{8}\,{\rm
M_{\odot}}$ as inferred from integrated photometry (Into & Portinari, 2013).
The profiles are shown in Fig. 1, with the faintest signals in the outskirts
reaching down to 30 mag arcsec-2. Even though the sky subtraction was
carefully carried out using ad-hoc custom-made procedures (Carlsten et al.,
2019), one cannot discard over or under corrections in the faint end of the
profiles ($\mu_{g}\gtrsim 29$ mag arcsec-2).

In principle, the EIM requires knowing the stellar mass density to be applied,
but we use $\mu_{g}$ directly because the shape of the DF is independent of a
scaling factor in $\Sigma_{o}$ (Sect. 2.2) and

$\log\Sigma_{o}=-\mu_{g}/2.5+\log\aleph=\log(\aleph F),$ (13)

with $\aleph$ the mass-to-light ratio. Obviously, for $\Sigma_{o}$ to scale as
the surface brightness $F$, $\aleph$ has to be constant along the profile,
which is a fair assumption for early-type galaxies, and the reason why we only
use the early-type galaxies in Carlsten et al. (2021) for our analysis111The
early-type galaxies in the sample have colors consistent with stellar
populations of the order of $6$ Gyr (Carlsten et al., 2021). Assuming a 10 %
age gradient between centers and outskirts (i.e., 0.6 Gyr), one expects a
gradient in $g-r$ around 0.03 mag (e.g., Bruzual & Charlot, 2003) which
translates into a change in $\log\aleph$ of 0.05 (e.g., Into & Portinari,
2013). Such change is negligible small compared with the observed range of
$-\mu_{g}/2.5$ ($\gtrsim 2$; see Fig. 2).. Carlsten et al. (2021) classified
early and late-type galaxies by eyeball inspection of their morphology so that
smooth featureless objects are early types while objects with star-forming
regions, blue clumps, or dust lanes are late types. The original set contains
111 early-type galaxies, but it was additionally purged to leave 93 as
described in Sect. 4.

The volume density profiles required to compute the phase-space distribution
function through the EIM (Eq. [3]) are found by fitting various surface
density profiles to the observed $-\mu_{g}/2.5$. The fits are carried out
using a least-squares minimization routine based on the Levenberg–Marquardt
algorithm, reproducing the procedure employed in Sánchez Almeida et al.
(2021). We start off from a family of volume density profiles $\rho_{\xi}(r)$,
with $\xi$ representing the free parameters that characterize the profile.
Even when $\rho_{\xi}(r)$ is analytical, its 2D projection is not, therefore,
Eq. (12) has to be evaluated numerically to carry out the fits. To speed up
the procedure, the surface densities required for fitting result from
interpolating on a precomputed grid $\Sigma_{\xi}(R)$, related to
$\rho_{\xi}(r)$ through Eq. (12). We employ a linear interpolation in as many
dimensions as needed to account for the number of free parameters defining
$\rho_{\xi}(r)$. The fits do not include the effect of the point-spread-
function (PSF) on the shape of the synthetic profiles since the observed
profiles are well resolved, even at the innermost radius (100 pc, Fig. 1). At
the average distance of the observed galaxies ($\sim 7$ Mpc), 100 pc
corresponds to 3″, which is significantly larger than the typical width of the
PSF ($\lesssim 1\arcsec$; Carlsten et al., 2021).

Figure 2 shows examples of this kind of fit to profiles in Fig. 1, where the
density is fitted using polytropes. In this case the profiles depend on three
parameters: a global density scaling, a spatial scaling $r_{s}$, and the
polytropic index $m$ (see, Sánchez Almeida et al., 2021). Polytropes are
relevant in this context because they represent self-gravitating systems of
maximum Tsallis entropy and so they describe the thermodynamical equilibrium
of systems experiencing long-range forces (Plastino & Plastino, 1993; Sánchez
Almeida, 2022b). Polytropes are know to describe the distribution of DM
observed in dwarf galaxies (Sánchez Almeida et al., 2020), the starlight
observed in galaxies of various masses (Sánchez Almeida et al., 2021), as well
as the DM distribution in numerical simulations of self-gravitating systems
reaching thermodynamical equilibrium (Sánchez Almeida & Trujillo, 2021).

Figure 2 shows both good (top three rows) and bad (bottom row) fits to some of
the observed profiles. The quality control parameter, called “Relative RMS” in
Fig. 2 and throughout the text, is defined as,

${\rm Relative\leavevmode\nobreak\ RMS}=\frac{\rm RMS}{\max(\log F)-\min(\log
F)},$ (14)

with RMS the root mean square of the residuals of the fit and $\log F$ the
logarithm of the observed flux (Eq. [13]). It represents the error in the fit
relative to the amplitude of the observed signal. It is a good proxy for the
quality of the fits so that ${\rm Relative\leavevmode\nobreak\ RMS}\simeq
0.02$ separates the good fits ($<0.02$) form the fair fits ($>0.02$) as judged
rather arbitrarily from the eye-ball inspection of the individual fits (see
Fig. 3).

Figure 3: Residuals of the fits to projected polytropes of the profiles in
Fig. 1. The color code refers to the maximum Relative RMS defining each set of
profiles. We define as good fits those represented in blue and orange.

All polytropes have cores in the sense given in Eq. (5), therefore, to broaden
the range of possibilities, we also try other more flexible density profiles,
namely, the $\rho_{abc}$ profiles defined as

$\rho_{abc}=\frac{\rho_{s}}{x^{c}(1+x^{a})^{(b-c)/a}},$ (15)

where $x=r/r_{s}$, and $\rho_{s}$ and $r_{s}$ are scaling constant in density
and radius, 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,{\rm\leavevmode\nobreak\ and\leavevmode\nobreak\ }c=1$) and
the $m=5$ polytrope (a.k.a. Schuster-Plummer profile, with
$a=2,b=5,{\rm\leavevmode\nobreak\ and\leavevmode\nobreak\ }c=0$). Moreover,
for $a=2$, $b=m$, and $c=0$, i.e.,

$\rho_{2m0}=\frac{\rho_{s}}{\left[1+(r/r_{s})^{2}\right]^{m/2}},$ (16)

it approximately accounts for the inner region of a polytrope of arbitrary
index $m$ (e.g., Sánchez Almeida, 2022b). The $\rho_{abc}$ profile is more
flexible than the polytropes because the inner and outer slopes can be fitted
independently: the parameter $b$ encodes the outer logarithmic slope of the
profile

$\lim_{r\to\infty}\frac{d\log\rho_{abc}}{d\log r}=-b,$ (17)

whereas $c$ encodes the inner logarithmic slope (as in Eq. [9]). Examples of
these alternative fits will be given later on in connection with the results
(Sect. 5).

Figure 4: Constraints on the underlying potential inferred from the polytropic
fit of a representative galaxy. The corresponding fit is shown in Fig. 2
(first column second row). Each potential is characterized by two parameters
as shown in Eq. (18): the inner slope, $c_{p}$, and the global spatial scale,
$r_{sp}$. Red dots correspond to $f<0$ and so mark unphysical potentials
whereas blue dots correspond to $f\geq 0$ and so trace potentials consistent
with the density profile. The parameter $r_{s}$ gives the spatial scale of the
stellar light distribution. The horizontal lines define NFW potentials
($c_{p}=1$; the solid gray line) and cored potentials ($c_{p}=0$; the dashed
gray line). Figure 5: Constraints on the underlying potential imposed by the
set of polytropic fits like those in Fig. 2. Every potential is characterized
by $c_{p}$ and $r_{sp}$ and the color code shows the frequency of having
$f\geq 0$ considering all good and fair fits (93 galaxies). The whole region
of $c_{p}>0$ is discarded whereas $c_{p}\leq 0$ is permitted provided
$r_{sp}\gtrsim r_{s}/2$. The diagram does not change much when only good fits
are considered (Relative RMS $<0.02$) or only low mass galaxies are selected
($M_{\star}\leq 10^{7}\,{\rm M_{\odot}}$). The parameter $r_{s}$ gives the
spatial scale of the stellar light distribution. Horizontal lines defining NFW
potentials ($c_{p}=1$; the solid gray line) and cored potentials ($c_{p}=0$;
the dashed gray line) are included.

##  4 Diagnostic diagrams

This section describes the practical realization of steps 3 and 4 in the
procedure sketched out in Sect. 2.2.

To constrain the properties of the gravitational potential, we compute for
each observed profile the DFs in the phase-space corresponding to a battery of
potentials. To make the interpretation intuitive, the potentials are described
in terms the mass density profile that creates the potential $\rho_{p}(r)$
(Eq. [11]). In general the relation between $\rho_{p}$ and $\Psi$ is not
analytical and has to be computed numerically. For $\rho_{p}$ we use a
particular type of $\rho_{abc}$ function (Eq. [15]),

$\rho_{p}(r)=\frac{\rho_{sp}}{x^{c_{p}}(1+x^{2-{c_{p}}})^{(5-3{c_{p}})/(2-{c_{p}})}},$
(18)

that scans seamlessly from a Schuster-Plummer profile to a NFW profile when
$c_{p}$ goes from 0 to 1. Thus, the potentials in the battery are defined in
terms of the two free parameters that define $\rho_{p}$, i.e., $c_{p}$ and
$r_{sp}$. (The scaling factor $\rho_{sp}$ is irrelevant for the EIM, as
pointed out in Sect. 2.1.) The key question of the analysis is knowing whether
the DF resulting from the observed profile and each potential is or not
negative. The answer to the question can be summarized in the kind of
diagnostic diagram shown in Fig. 4. Red dots correspond to $f<0$ and so mark
unphysical potentials whereas blue dots corresponds to $f\geq 0$ and so trace
gravitational potentials consistent with the density profile.

Figure 4 shows a typical example corresponding to galaxy # 15, with good
polytropic fit in Fig. 2. Here an throughout the paper $f$ is computed using
the methods and tools developed by Sánchez Almeida et al. (2023). These tools
only work for $\rho_{abc}$ profiles, therefore, for this calculation
polytropes were approximated using Eq. (16). Note that all $c_{p}>0$ are
discarded in Fig. 4. All potentials have to derive from a cored density
profile to be consistent with the fitted polytropes. This comes with no
surprise since all polytropes have a core (e.g., Sánchez Almeida et al., 2020)
and we proved that a core in stars requires a core in the potential for
spherically symmetric systems and isotropic velocity distributions (see
details in Sect. 2.1). In addition, some of the models with $r_{s}\geq r_{sp}$
are also incompatible. We also knew from Sánchez Almeida et al. (2023) of the
inconsistency of potentials with characteristic radii ($r_{sp}$) significantly
smaller than the stellar core radii ($r_{s}$). Perhaps more surprising is the
fact that no incompatibility is found for some potentials with $c_{p}<0$,
i.e., where the total density goes to cero when $r\to 0$. These potentials are
included in the battery for comprehensiveness although they are expected to be
Rayleigh–Taylor unstable (e.g., Sánchez Almeida, 2022a) and so of difficult
practical realization. Once again, Sánchez Almeida et al. (2023) showed that
in this case ($c=0$, $a=2$, and $c_{p}<0$) the derivative $d\rho/d\Psi>0$ when
$r\to 0$ and so there is no need for $f$ to be negative (see their Appendix
E).

From the 111 galaxies in Fig. 1, we discard those having bad fits (like those
of the bottom row in Fig. 2), and then are left with 93 good and fair fits.
For every one of these profiles, we have a diagnostic diagram like the example
shown in Fig. 4. We construct a second kind of diagnostic diagram to show the
constraints on the underlaying potential imposed for the whole set of observed
profiles. For every potential, characterized by $c_{p}$ and $r_{s}/r_{sp}$, we
count the frequency of $f\geq 0$ considering all the observed galaxies. The
result is shown in Fig. 5. As we knew, the whole region of $c_{p}>0$ is
discarded by all (100 %) of the fits. The region with $c_{p}\leq 0$ is
permitted provided $r_{sp}\gtrsim r_{s}/2$. The results do not change much
when only the good fits are used (Relative RMS $<0.02$) or when only the low
mass galaxies are considered ($M_{\star}<10^{7}\,{\rm M_{\odot}}$).

Formal error bars can be assigned to this kind of diagnostic diagram as
follows. It is based on counting the number of galaxies that, for the same
potential, do or do not have a particular property (i.e., whether or not
$f\geq 0$ everywhere). Thus, the actual number of counts should follow a
binomial distribution. Using this idea, we work out in Appendix A the error
bars to be expected. They are largest when the percentage is 50 %, and tend to
zero at 0 % and 100 %. The actual error depends on the number of galaxies
employed in the diagram. For $\sim 50$ galaxies (typical number of good
galaxies or low-mass galaxies in our sample), it is smaller than 6 % (Fig.
15).

##  5 Results

Armed with the fits described in Sect. 3 and the diagnostic diagrams of Sect.
4, here we describe the constraints on the potentials set by the observed
surface brightness profiles. The degrees of freedom in the fits increase from
polytropes (Sect. 5.1) to variable inner slope (Sect. 5.2) to allow both the
inner and outer slopes to vary simultaneously (Sect. 5.3).

Figure 6: Quantitative summary with the results of the polytropic fit to the
profiles in Fig. 1. Low ($\leq 10^{7}\,{\rm M_{\odot}}$) and high
($>10^{7}\,{\rm M_{\odot}}$) stellar mass galaxies are represented with a
different color as indicated in the bottom panel. Top panel: dependence of the
inferred polytropic index $m$ on the quality of the fits, as parameterized by
the Relative RMS. There seems to be a trend for the best quality fits to have
$m\sim 5$ (marked with the gray solid horizontal line). Middle panel: scatter
plot with the size of the core versus the goodness of the fit. Bottom panel:
scatter plot with the stellar mass versus the goodness of the fit. More
massive galaxies lead to slightly better fits.  Figure 7: Quantitative summary
with the results of the polytropic fit to the profiles in Fig. 1. Medium
(Relative RMS $>0.02$) and high ($\leq 0.02$) quality fits are represented
with different colors as indicated in the bottom panel. Top panel: scatter
plot of the polytropic index $m$ as a function the stellar mass of the galaxy.
Middle panel: variation of the core radius with stellar mass. More massive
galaxies tend to have larger cores. Bottom panel: scatter plot with the
variation of the absolute RMS with stellar mass. Good fits (blue dots) have
absolute RMS $\lesssim 0.03$ dex.

Figure 8: Fits similar to Fig. 2 but allowing the inner slope of the stellar
mass profile to vary (green dashed lines), and the inner and outer slopes to
vary independently (orange lines). The density used for fitting was assumed to
be $\rho_{abc}$ in Eq. (15) with different bonds between $a$, $b$, and $c$ for
the green and the orange lines. See the main text and the caption of Fig. 2
for further details. Figure 9: Scatter plot of the Relative RMS resulting from
variable inner slope fits and from polytropic fits. Different colors designate
different quality of the polytropic fits (see the inset). The points
corresponding to low mass galaxies are encircled in red. The absolute number
of good polytropic fits is larger as is also the relative fraction of good
over bad fits.

###  5.1 Polytropic fits, with zero inner slope

Examples of polytropic fits are shown in Fig. 2. The quality of these fits is
quantified using the Relative RMS (Eq. [14]) and overall, the quality is
satisfactory (see the residuals in Fig. 3). As we mentioned in Sect. 4, from
the 111 original profiles in Fig. 1, 18 are discarded because the profile is
not smooth and the fits are bad (like those of the bottom row of Fig. 2). From
the remaining 93, 41 are good (44 %, with Relative RMS$<0.02$) and the rest
fair (64 %).

The results of the fits are summarized in Figs. 6 and 7. The top panel of Fig.
6 shows the dependence of the inferred polytropic index $m$ as a function the
quality of the fits parameterized by the Relative RMS. There seems to be a
trend for the best quality polytropic fits to have $m\sim 5$, a behavior also
found in other galaxy samples (Sánchez Almeida et al., 2021). This trend does
not depend on the stellar mass, divided in low ($\leq$ $10^{7}M_{\odot}$) and
high ($>$ $10^{7}M_{\odot}$) mass bins for the plot. The middle panel of Fig.
6 contains the size of the core versus the goodness of the fit. No trend is
found. The bottom panel of Fig. 6 shows the stellar mass versus the goodness
of the fit, and there is a hint of trend so that fits are slightly better for
the more massive galaxies. This effect may be an observational bias since more
massive objects tend to be brighter and thus present better observations that
leave smaller residuals. The top panel of Fig. 7 shows the polytropic index as
a function the stellar mass of the galaxy. Overall, there is no clear trend of
$m$ with $\log M_{\star}$, which may be due to the fact that the lowest mass
galaxies happen to have the worst fits where $m$ is more uncertain. However,
the good fits seem to indicate a trend for $m$ to grow with increasing
$M_{\star}$ (the blue symbols in Fig. 7, top panel). The middle panel of Fig.
7 displays the variation of the core radius with stellar mass. More massive
galaxies tend to have larger cores. This trend does not depend on the quality
of the fits, represented with two colors depending on whether the Relative RMS
is smaller or larger than 0.02. The bottom panel of Fig. 7 shows the variation
of the absolute RMS with stellar mass. Note that the good fits have an
absolute RMS better than 0.03 dex, equivalent to 7 % in flux.

With the polytropic fits described above, we constructed the diagnostic
diagrams shown in Figs. 4 and 5. They were discussed as examples to present
the diagnostic diagrams in Sect. 4, and there is little else to add except to
emphasize that all kinds of cuspy potentials are discarded ($c_{p}>0$) and the
characteristic radius of the potential cannot be significantly smaller than
the radius of the stellar distribution ($r_{sp}\gtrsim r_{s}/2$). The first
result was to be expected since given the hypotheses made in our application
of the EIM, a core in stars necessarily needs of a core in the potential, and
all polytropes are cored profiles independently of their index.

Figure 10: Quantitative summary with the results of the variable inner slope
fit. Except for the inner slope $c$, the rest is similar to Fig. 6 describing
polytropic fits. Note how $c$ is often negative, a trend enhanced as the
quality of the fits worsens. Figure 11: Quantitative summary with the results
of the variable inner slope fit. Except for the inner slope $c$ replacing the
polytropic index, the figure is similar to Fig. 7 describing the polytropic
fits.

Figure 12: Constraints on the underlying gravitational potential imposed by
the variable inner slope fits. Top panel: diagnostic diagram, equivalent to
Fig. 5 for the polytropic fits, which includes all galaxies together. Note
that around half of the galaxies have unphysical $f<0$ irrespectively of the
inner slope of the potential ($c_{p}$), a percentage that increases as
$r_{s}/r_{sp}$ increases. This is mostly due to the inner slope $c$ being
negative for a number of galaxies. Middle panel: similar to the top panel but
including only good fits (Relative RMS $<0.02$). Note how the fraction of
unphysical potentials with $f<0$ has been significantly reduced (the tone
becomes bluish rather than greenish). Bottom panel: similar to the top panel
but including only low mass galaxies ($M_{\star}<$ $10^{7}\,{\rm M_{\odot}}$).
As expected, the fraction of potentials with $f<0$ increases (the tone becomes
reddish) due to the increase of fits with $c<$ 0\. The color bar shows % of
galaxies with $f\geq 0$.

###  5.2 Variable inner slope fits

Going a step further, we no longer force a core in the stellar density profile
and the fits allows for a inner slope different from zero. The shape of the
stellar density profiles is taken from Eq. (18) which scans from cored to NFW
profiles when the inner slope $c$ goes from 0 to 1. Examples of these fits are
shown as green dashed lines in Fig. 8, which contains exactly the same
galaxies as Fig. 2. In general, the fits using projected polytropes seem to do
a slightly better job than these other fits, even though the number of free
parameters is the same in both cases: $m$ has been replaced with $c$. This may
be due to the fact that polytropes, rather than $\rho_{abc}$ profiles, are to
be expected when self-gravitating systems reach thermodynamic equilibrium
(Plastino & Plastino, 1993). Figure 9 shows a scatter plot between the
Relative RMS obtained with the two types of fit. Even if the differences are
not overwhelming, the polytropic fits tend to be better than for these other
fits. The fraction of good fits (Relative RMS $<$ 0.02) is larger for
polytropes than for the variable inner slope fits (44 % versus 32 %).
Considering only good and fair fits, polytropic fits have lower Relative RMS
in 59 % of the cases. Finally, the number of times where the polytropic fits
are significanly better than the alternative (defined as ratio of Relative RMS
$>$ 1.5) is much larger than the oposite (19 galaxies versus 3 galaxies).

The properties of the variable inner slope fits are summarized in Figs. 10 and
11. The top panel of Fig. 10 shows the dependence of the inner slope $c$ as a
function the Relative RMS. Note that often the best fitting $c$ is negative.
Most of these negative slopes may be produced by systematic errors in the fit,
a conjecture supported by two facts: (1) the number of galaxies with $c<0$
increases with increasing Relative RMS (Fig. 10, top panel) and (2) polytropic
fits, where $c=0$ is imposed, often provide better fits, as discussed in the
previous paragraph. The middle panel of Fig. 10 contains the size of the core
versus the goodness of the fit. No obvious trend is found. The bottom panel of
Fig. 10 shows the stellar mass versus the goodness of the fit indicating a
hint of trend where more massive galaxies show better fits. The top panel of
Fig. 11 shows a scatter plot of the inner slope $c$ as a function the stellar
mass. There is no clear trend. The middle panel of Fig. 11 displays the
variation of the core radius with stellar mass. More massive galaxies tend to
have larger cores, as it happened with polytropic fits (Fig. 7, middle panel).
The bottom panel of Fig. 11 shows the variation of the absolute RMS with
stellar mass.

The diagnostic diagrams to constrain the properties of the gravitational
potential consistent with these fits are shown in Fig. 12. Contrary to what
happens with the polytropic fits, where most of the potentials are either
allowed or forbidden (Fig. 5), most of the potentials are sometimes allowed
and sometimes forbidden. The results for all galaxies are shown in the top
panel of Fig. 12. Around 50 % of the potentials are not allowed, with a formal
error in this estimate of the order of 5 – 10 % (Appendix A). This happens
irrespectively of the inner slope of the potential $c_{p}$, with the fraction
increasing with increasing $r_{s}/r_{sp}$. As we discussed above (Fig. 10, top
panel), a significant part of the fitted profiles have $c\leq 0$, and these
profiles are inconsistent with any potential. The middle panel of Fig. 12
includes only good fits (Relative RMS $<0.02$). The fraction of unphysical
potentials is significantly reduced as expected following the reduction of
objects with $c<0$ (Fig. 10, top panel). The bottom panel of Fig. 12 shows
only low mass galaxies ($M_{\star}<$ $10^{7}\,{\rm M_{\odot}}$). This other
time the fraction of unphysical potentials increases due to the increase of
the percentage of fits with $c<$ 0 among low-mass objects. Thus, we suspect
the diagnostic plots to be significantly biased by the large uncertainty in
$c$, which either allows for many potentials to be consistent ($c$
significantly larger than 0) or inconsistent ($c<0$). The shape chosen for
fitting does not constrain well enough the central slope of the observed
galaxy. Or equivalently, the errors (stochastic and systematic) are still not
small enough to reach a solid conclusion. We note that the scenario does not
change significantly even when the Relative RMS is chosen as small as $<0.01$
(Fig. 3). In the next section, we try another type of fit to check whether it
improves the situation. The fitting functions used here tie the inner and
outer slopes of the density profile, and this artificial bound may be causing
significant systematic errors.

Figure 13: Top panel: Scatter plot Relative RMS from the variable inner and
outer slope fits versus the variable inner slope fits. It is clear that the
variable inner and outer slope fits generally has lower RMS, as expected since
the underlying shape is the same and the former has one more degree of
freedom. Bottom panel: Scatter plot Relative RMS from the variable inner and
outer slope fits versus polytropic fits. The colors and symbols are the same
as in Fig. 9.  Figure 14: Summary of the results with free inner ($c$) and
outer ($b$) slopes presented as a function of the Relative RMS of the fit. To
be compared with Fig. 6, for the polytropes, and Fig. 10, for fits when only
the inner slope is free. The color code separates low and high stellar mass
galaxies, as indicated in the inset. Top panel: Inner slope. See how often
$c<0$, particularly for poor fits. Middle panel: Outer slope. For reference,
$b=5$ in the case of polytropic fits with $m=5$, which is the limit $m$ for
good polytropic fits (Fig. 6). Bottom panel: Stellar mass. Lower mass galaxies
tend to have poorer fits.

###  5.3 Variable inner and outer slopes

Even though the variable inner slope profiles have the same number of free
parameters as the polytropic fits, they tend to provide slightly worst fits
(Fig. 9). The problem could be due to the fact that the inner slope in these
fits, $c$, also sets the outer slope, $b=5-2c$ (c.f. Eqs. [15] and [18]). This
tie minimizes the number of free parameters in the fit, but there is no clear
physical reason to expect it. To overcome the potential problem created by the
artificial link between $b$ and $c$, we carried out additional fits where both
the inner and outer slopes in Eq. (15) are free parameters. The third
parameter of the $abc$ profiles, $a$, which determines the transition between
the inner and outer parts, are set to 2, which is the value corresponding to
the polytropes222We also tried $a=1$ but the fits are worse than those for
$a=2$ reported here. (Eq. [16]). Examples of fits with free inner and outer
slopes are shown as orange lines in Fig. 8. The figure also includes the fits
with variable inner slope only (green dashed lines). Even though the Relative
RMS is lower when both $c$ and $b$ are free parameters (top panel, Fig. 13),
the fitted profiles look very similar, and also similar to the observed
profiles (Fig. 8). Figure 14 shows how $c$ and $b$ depend on the Relative RMS
of the fit.

The constraints on the gravitational potential imposed by these fits are
similar to those when only the inner slope is free (Fig. 12). There is no
major qualitative difference. A significant part of the fits still render
$c\leq 0$ and therefore unphysical $f<0$ DFs. The galaxies can be divided into
those with $c\gtrsim 0.1$, which are consistent with a NFW potential
($c_{p}=1$) and those with $c\lesssim 0.1$ which are inconsistent (see Sect.
2.1). However, part of the scatter in $c$ is probably artificial as deduced
from the fact that the scatter increases with increasing Relative RMS (Fig.
14, top panel). The uncertainties in the fits are partly due to errors in the
photometry but, most likely, the are systematic arising from deviations from
the assumptions implicit in the fits, including spherical symmetry,
$\rho_{abc}$ profile, and $a=2$. We have tried to estimate the impact of these
systematic errors on the inner slope $c$ under the working hypothesis that all
negative slopes $c<0$ are due to noise. We select all those profiles with good
fits (Relative RMS $<0.02$) and $c<0$. If these values are random excursions
produced by noise of profiles with $c\simeq 0$, then the standard deviation of
the noise would be $\sigma_{c}\simeq 0.3$. (The same exercise with the fits in
Sect. 5.2, where only the inner slope was fitted, renders $\sigma_{c}\geq
0.25$.) Thus, we find that 38 % of the good fits have $|c|\leq\sigma_{c}$ and
therefore are consistent with $c=0$ within errors. (The same exercise with the
inner slope free fits renders 60 %.) In other words, around 40 % of the fits
are consistent with cores in the stellar mass distribution (i.e., with $c=0$)
and so likely inconsistent with NFW-like potentials. This figure is also
similar to the limit steaming from the diagnostic plot (e.g., Fig. 12, middle
panel).

The Relative RMS of the free inner and outer slope fits is in general smaller
than for polytropes (Fig. 13, bottom panel). However, the comparison of the
two types of fits is not direct since two-slopes fits have one more degree of
freedom than polytropes, and since the two shapes are different, an increase
of the number of free parameters my artificially reduce the RMS because of
overfitting. This being said, we find that the number of profiles with
Relative RMS $<0.02$ is very similar in both cases (41 for polytropic fits
versus 39 for $\rho_{2bc}$) so is the number of profiles with significantly
better fit of one relative to the other (10 versus 12). The similitude between
the numbers characterizing both types of fit supports the above conclusion
that the observed profiles are often consistent with $c=0$. The core ($c=0$)
is imposed for polytropes but inferred for fits with free inner and outer
slopes.

##  6 Conclusions

Sánchez Almeida et al. (2023) proposed the use of the classical Eddington
Inversion Method (EIM) to constrain the properties of the gravitational
potential in galaxies using only the observed distribution of starlight. As we
explain in Sect. 1, the procedure may be able to decipher the inner shape of
the DM halos in the ultra-low mass regime where (1) the potential is set by
the DM, and (2) the baryon feedback processes are unable to transform cuspy DM
halos (typical of collisionless DM) into cored halos. Thus, in this mass
regime, the finding of cores in the DM distribution would imply that the DM is
not collisionless. Even though dynamical measurements of the DM distribution
in these galaxies are unfeasible, the EIM may allow them, thus opening a new
possibility to explore the nature of DM.

The algorithm employed in the paper is detailed in Sect. 2.2 and tested using
GCs in Appendix C. Given an observed stellar mass distribution, we check
whether or not it is consistent with a battery of gravitational potentials.
The inconsistency manifest itself when, for the stellar mass to reside in a
particular potential, the distribution function in the phase space $f$ has to
adopt negative values in some regions of the phase space. The EIM allows us to
compute $f$ under the assumption that the system is spherically symmetric with
isotropic velocities. Here we present the first practical application of the
technique using surface brightness profiles of around 100 low mass galaxies
($10^{6}\,{\rm M_{\odot}}\leq M_{\star}\leq 10^{8}\,{\rm M_{\odot}}$) from
Carlsten et al. (2021, 2022), who selected them as dwarf satellites of Milky
Way-like galaxies (see Sect. 3). Their masses are still too large to constrain
the nature of DM but the set provides an excellent testbed for the technique
in terms of number of targets and noise level. The technique depends on the
type of function used for fitting the observed profiles. We try three types:
polytropes, profiles with free inner slope (Eq. [18]), and profiles with free
inner and outer slopes (Eq. [15], with $a=2$). The application of the
technique allows us to reach the following conclusions:

  * •

The observed profiles are generally well fitted by projected polytropes, with
41 galaxies ($44$ % of the total) having pretty low residuals with a Relative
RMS (Eq. [14]) $<0.02$. The polytropic index tends to 5 for the best fits and
the derived core radius scales with galaxy stellar mass (Fig. 6). The good
fits have a RMS $\sim 0.03$ dex, equivalent to 7 % in flux (Fig. 7).

  * •

As expected (Sect. 2.1), NFW-like potentials are discarded if the polytropic
shape is good enough to reproduce the inner parts of the galaxies (Fig. 5).

  * •

More massive (more luminous) galaxies tend to have better fits, supporting
that noise (random and systematic) is affecting the results. The trend is
present for the three types of fits (Figs. 6, 10, and 14, bottom panels).

  * •

When a parameterized profile with variable inner slope is used for fitting,
then the quality of the fits worsens with respect to the polytropes, even
though both profiles have the same number of free parameters (Fig. 9). In this
case, the inner slope $c$ is often negative (Fig. 10), and the diagnostic
diagrams basically discard all potentials if $c\lesssim 0.1$ (Fig. 12).

  * •

When that inner and outer slopes are allowed to vary simultaneously (Sect.
5.3), the quality of the fits improves (Fig. 13, left panel), but the fact
that $c$ is often negative remains (Fig. 14, top panel). Some of theses
negative slopes may be produced by systematic errors in the fit, a conjecture
supported by the fact that the number of galaxies with $c<0$ increases with
increasing Relative RMS (Fig. 14, top panel). Assuming that the negative
excursions of $c$ are due to errors in the fit, we estimate the error in $c$
to be $\sigma_{c}\simeq 0.3$, which implies that around 40 % of the good fits
are consistent with $c=0$ (Sect. 5.3).

  * •

If polytropes were the correct fits (which is likely the case; see Sect. 5.3)
then all galaxies are inconsistent with NFW-like potentials (Fig. 5). In the
case of allowing the inner slope to vary (Sects. 5.2 and 5.3), between 40 %
and 70 % of the galaxies are consistent with cores in the stellar mass
distribution ($c=0$) thus inconsistent with NFW-like potentials (Fig. 12).

  * •

The EIM tools that we developed and use assume isotropic velocities. However,
the inner slope of a stellar core and the radial anisotropy $\beta$ (Eq. [7])
are related so that $c>2\beta$ for the system to be physically realizable. In
other words, large radially biased orbits ($\beta>0$) are strongly
inconsistent with soft stellar cores even if $c$ is significantly larger than
zero. In low mass galaxies, $\beta>0$ is to be expected (e.g., El-Badry et
al., 2017; Orkney et al., 2021) and so the limits on $c$ in the previous item
are likely to be conservative.

The above results have been derived under the assumption that the observed
stellar systems are spherically symmetric. However, dwarfs tend to be oblate
or even triaxials (e.g., Sánchez-Janssen et al., 2016; Putko et al., 2019).
The construction of the 1D profiles in Fig. 1 includes symmetrizing the
observed profiles, but it is unclear to what extent this artifice is enough to
grant the application of an approach which, strictly speaking, is valid only
for spherically symmetric systems. Fortunately, this assumption does not seem
to be critical since it does not compromise the main observational constraint
leading to the results, namely, the existence of stellar cores. Firstly, all
good observed profiles converge to the same shape, a polytrope of index $m\sim
5$ (Sect. 5.1). Thus, since this shape seems to be independent of the original
axial ratio, it is also the one expected in purely spherical stellar system,
for which our analysis applies. Secondly, one can show (Appendix B) how the
incompatibility between NFW potentials and stellar cores also remains for axi-
symmetric systems (e.g., oblate), therefore, it is more fundamental than, and
not attached to, the spherical symmetry assumption.

In short, the EIM-based technique seems to work in the sense of providing
constraints on the gravitational potential using only good but typical
observed surface brightness profiles. Our tests also indicate that good data
are needed to break the degeneracies in the determination of the inner slope
of the profiles. However, since the key parameter is the inner slope, reaching
the innermost regions is more important than decreasing noise. Thus, observing
with enough spatial resolution is critical. Having large samples of objects is
also important, so that those deviating from the EIM assumptions can be
discarded without compromising the large statistics needed to infer general
conclusions on the nature of DM.

Even if this is not mentioned explicitly, the gravitational potential to be
used in the EIM also includes the contribution from baryons. All the above
conclusions still hold provided the total density of the system (i.e., stars,
gas, and DM included) is well approximated by the density profile used to
construct the battery of potentials. Equation (15) was used in this paper, but
this is not a serious restriction of the technique since other shapes can be
incorporated if needed.

###### Acknowledgements.

Thanks are due to Scott Carlsten who generously provided the surface density
profiles analyzed in the paper and to Barbara Lanzoni and Francesco Ferraro
for the globular cluster profiles used in Appendix C. Thanks are also due to
an anonymous referee who suggested clarifications leading to the extension to
non-spherically symmetric systems in Appendix B and the sanity check in
Appendix C. JSA acknowledges financial support from the Spanish Ministry of
Science and Innovation, project PID2022-136598NB-C31 (ESTALLIDOS), 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. We
acknowledge the use of the Python packages numpy (Harris et al., 2020), scipy
(Virtanen et al., 2020), and matplotlib (Hunter, 2007).

## References

  * An & Zhao (2013) An, J. & Zhao, H. 2013, MNRAS, 428, 2805
  * An & Evans (2006) An, J. H. & Evans, N. W. 2006, ApJ, 642, 752
  * Bechtol et al. (2022) Bechtol, K., Birrer, S., Cyr-Racine, F.-Y., et al. 2022, arXiv e-prints, arXiv:2203.07354
  * Binney & Tremaine (2008) Binney, J. & Tremaine, S. 2008, Galactic Dynamics: Second Edition
  * Bruzual & Charlot (2003) Bruzual, G. & Charlot, S. 2003, MNRAS, 344, 1000
  * Carlsten et al. (2019) Carlsten, S. G., Beaton, R. L., Greco, J. P., & Greene, J. E. 2019, ApJ, 879, 13
  * Carlsten et al. (2022) Carlsten, S. G., Greene, J. E., Beaton, R. L., Danieli, S., & Greco, J. P. 2022, ApJ, 933, 47
  * Carlsten et al. (2021) Carlsten, S. G., Greene, J. E., Greco, J. P., Beaton, R. L., & Kado-Fong, E. 2021, ApJ, 922, 267
  * Carr et al. (2023) Carr, B., Clesse, S., Garcia-Bellido, J., Hawkins, M., & Kuhnel, F. 2023, arXiv e-prints, arXiv:2306.03903
  * Chan et al. (2015) Chan, T. K., Kereš, D., Oñorbe, J., et al. 2015, MNRAS, 454, 2981
  * Ciotti (2021) Ciotti, L. 2021, Introduction to Stellar Dynamics
  * Dey et al. (2019) Dey, A., Schlegel, D. J., Lang, D., et al. 2019, AJ, 157, 168
  * Di Cintio et al. (2014a) Di Cintio, A., Brook, C. B., Dutton, A. A., et al. 2014a, MNRAS, 441, 2986
  * Di Cintio et al. (2014b) Di Cintio, A., Brook, C. B., Macciò, A. V., et al. 2014b, MNRAS, 437, 415
  * Dodelson & Widrow (1994) Dodelson, S. & Widrow, L. M. 1994, Phys. Rev. Lett., 72, 17
  * Eddington (1916) Eddington, A. S. 1916, MNRAS, 76, 572
  * El-Badry et al. (2017) El-Badry, K., Wetzel, A. R., Geha, M., et al. 2017, ApJ, 835, 193
  * Expósito-Márquez et al. (2023) Expósito-Márquez, J., Brook, C. B., Huertas-Company, M., et al. 2023, MNRAS, 519, 4384
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357
  * Hayashi et al. (2020) Hayashi, K., Chiba, M., & Ishiyama, T. 2020, ApJ, 904, 45
  * Hernquist (1990) Hernquist, L. 1990, ApJ, 356, 359
  * Hu et al. (2000) Hu, W., Barkana, R., & Gruzinov, A. 2000, Phys. Rev. Lett., 85, 1158
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90
  * Ibata et al. (2013) Ibata, R., Nipoti, C., Sollima, A., et al. 2013, MNRAS, 428, 3648
  * Into & Portinari (2013) Into, T. & Portinari, L. 2013, MNRAS, 430, 2715
  * Ivezić et al. (2019) Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111
  * Jackson et al. (2021) Jackson, R. A., Martin, G., Kaviraj, S., et al. 2021, MNRAS, 502, 4262
  * Lacroix et al. (2018) Lacroix, T., Stref, M., & Lavalle, J. 2018, J. Cosmology Astropart. Phys., 2018, 040
  * Lan et al. (2016) Lan, T.-W., Ménard, B., & Mo, H. 2016, MNRAS, 459, 3998
  * Lynden-Bell (1962) Lynden-Bell, D. 1962, MNRAS, 123, 447
  * Martin (1971) Martin, B. R. 1971, Statistics for physicists
  * Merritt et al. (2006) Merritt, D., Graham, A. W., Moore, B., Diemand, J., & Terzić, B. 2006, AJ, 132, 2685
  * Miocchi et al. (2013) Miocchi, P., Lanzoni, B., Ferraro, F. R., et al. 2013, ApJ, 774, 151
  * Navarro et al. (1997) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493
  * Orkney et al. (2021) Orkney, M. D. A., Read, J. I., Rey, M. P., et al. 2021, MNRAS, 504, 3509
  * Plastino & Plastino (1993) Plastino, A. R. & Plastino, A. 1993, Physics Letters A, 174, 384
  * Putko et al. (2019) Putko, J., Sánchez Almeida, J., Muñoz-Tuñón, C., et al. 2019, ApJ, 883, 10
  * Read et al. (2016) Read, J. I., Agertz, O., & Collins, M. L. M. 2016, MNRAS, 459, 2573
  * Sánchez Almeida (2022a) Sánchez Almeida, J. 2022a, ApJ, 940, 46
  * Sánchez Almeida (2022b) Sánchez Almeida, J. 2022b, Universe, 8, 214
  * Sánchez Almeida et al. (2017) Sánchez Almeida, J., Filho, M. E., Dalla Vecchia, C., & Skillman, E. D. 2017, ApJ, 835, 159
  * Sánchez Almeida et al. (2023) Sánchez Almeida, J., Plastino, A. R., & Trujillo, I. 2023, ApJ, 954, 153
  * Sánchez Almeida & Trujillo (2021) Sánchez Almeida, J. & Trujillo, I. 2021, MNRAS, 504, 2832
  * Sánchez Almeida et al. (2020) Sánchez Almeida, J., Trujillo, I., & Plastino, A. R. 2020, A&A, 642, L14
  * Sánchez Almeida et al. (2021) Sánchez Almeida, J., Trujillo, I., & Plastino, A. R. 2021, ApJ, 921, 125
  * Sánchez-Janssen et al. (2016) Sánchez-Janssen, R., Ferrarese, L., MacArthur, L. A., et al. 2016, ApJ, 820, 69
  * Sollima et al. (2012) Sollima, A., Bellazzini, M., & Lee, J. W. 2012, ApJ, 755, 156
  * Spergel & Steinhardt (2000) Spergel, D. N. & Steinhardt, P. J. 2000, Phys. Rev. Lett., 84, 3760
  * Trujillo et al. (2021) Trujillo, I., D’Onofrio, M., Zaritsky, D., et al. 2021, A&A, 654, A40
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261
  * Zou et al. (2018) Zou, H., Zhang, T., Zhou, Z., et al. 2018, ApJS, 237, 37
  * Zou et al. (2017) Zou, H., Zhou, X., Fan, X., et al. 2017, PASP, 129, 064101

##  Appendix A Error bars for the diagnostic plots

The diagnostic plot introduced in Sect. 4 and shown in Figs. 5 and 12 is based
on counting the number of galaxies that, for the same potential, do or do not
have a particular property (i.e., whether or not $f\geq 0$ everywhere). Thus,
the number of counts $n$ should follow a binomial distribution characterized
by the probability $q$ of having the property (e.g., Martin, 1971). In this
case, the expected value of counts is $qN_{g}$, with $N_{g}$ the number of
galaxies. The error in the number of counts, computed as the square root of
the variance of the distribution, is $\sqrt{q(1-q)N_{g}}$. Thus, the
percentage of galaxies represented in the diagnostic plot, $g=100\times
n/N_{g}$, has an error $\sigma_{g}$ of

$\sigma_{g}=100\times\sqrt{\frac{q(1-q)}{N_{g}}}.$ (19)

The error in the percentage versus the measured percentage is shown in Fig. 15
for three representative values of $N_{g}$ (25, 50, and 100). Errors are
largest when the percentage is 50 % , and tend to zero at 0 and 100 %. Their
values for the characteristic $N_{g}\sim 50$ are typically $\lesssim$ 6 %.

Figure 15: Errors expected in the diagnostic diagrams showing percentage of
galaxies with $f\geq 0$ everywhere. It depends on the number of galaxies used
for diagnostics, as indicated in the inset.

##  Appendix B Incompatibility between stellar cores and NFW potentials beyond
spherically symmetric systems

This appendix proves that the incompatibility between NFW potentials and
stellar cores goes beyond the assumption of spherical symmetry, being more
fundamental than and not attached to this particular hypothesis.

Lynden-Bell (1962) extended the EIM to axi-symmetric systems. The axis of
symmetry of the stars is assumed to be the same as the potential. In this
case, the even component of the distribution function, $f_{+}$, depends on two
integrals of motion, the relative energy $\epsilon$ and the angular momentum
along the axis of symmetry. The density depends on the cylindrical coordinates
$R$ (radius) and $z$ (distance along the axis of symmetry). Using the
potential $\Psi$ to parameterize the dependence on $z$ at a fixed $R,$ one can
write down an equation (see Binney & Tremaine, 2008, Sect. 4.4.1),

$\frac{\partial\rho(R,\Psi)}{\partial\Psi}=2\pi\sqrt{2}\int_{0}^{\Psi}\frac{f_{+}\left[\epsilon,2(\Psi-\epsilon)R^{2}\right]}{\sqrt{\Psi-\epsilon}}d\epsilon,$
(20)

which is formally identical to Eq. (2). This equation, particularized to
$R=0$, can be used to compute the derivative of the density along $z$,
rendering,

$\frac{\partial\rho(0,\Psi)/\partial z}{\partial\Psi/\partial
z}=2\pi\sqrt{2}\int_{0}^{\Psi}\frac{f_{+}\left[\epsilon,0\right]}{\sqrt{\Psi-\epsilon}}d\epsilon.$
(21)

As it happens with $f$, $f_{+}$ has to be non-negative everywhere in the phase
space (Binney & Tremaine, 2008, Eq. 4.39b), therefore, the arguments given in
Sect. 2.1 apply here as well. If there is a stellar core along the $z$
direction, i.e.,

$\partial\rho(0,\Psi)/\partial z=0{\rm\leavevmode\nobreak\
when\leavevmode\nobreak\ }z\to 0,$ (22)

then for $f_{+}\geq 0$ everywhere,

$\partial\Psi/\partial z=0{\rm\leavevmode\nobreak\ when\leavevmode\nobreak\
}z\to 0,$ (23)

which is incompatible with $\Psi$ being a NFW potential and therefore stellar
cores and NFW potential are physically incompatible also in axi-symmetric
systems. As we stress above, this fact proves that the incompatibility between
NFW potentials and stellar cores does not arise from the spherical symmetry
assumption.

##  Appendix C Sanity check based on Globular Clusters

Globular clusters (GC) are known to be self-gravitating stellar systems, with
little or no DM (e.g., Sollima et al., 2012; Ibata et al., 2013), therefore,
we know their gravitational potential directly from the observed stellar
distribution. Moreover, GCs usually show inner cores produced by the
thermalization of the system through two-body gravitational interactions
between stars (e.g., Binney & Tremaine, 2008). Thus, GCs provide excellent
testbeds for the EIM-based approach followed in the paper. This appendix
presents one such test using the GCs observed with the Hubble Space Telescope
by Miocchi et al. (2013), with their mass surface density profiles derived
from resolved stellar counts. The tools introduced in Sect. 2.2 are applied to
them. Specifically, the fits allow for free inner and outer slopes, as
described in Sect. 5.3. From the original 26 GCs, we keep 17 having good fits
to the innermost radii with cores. An example is shown in the top panel of
Fig. 16, which is equivalent to any panels in Figs. 2 or 8. The corresponding
diagnostic diagram (as in Fig. 4) is shown in the middle panel of Fig 16,
where the large star symbol shows the inner slope and characteristic radius of
the GC (upper panel in Fig. 16). If the GC is self-gravitating, these two
parameters from the stars define the gravitational potential. One can see how
the parameters from the fit appear in the allowed region of the diagnostic
diagram (indicated by the blue bullets). Moreover, the diagnostic plot
completely discards NGC 288 residing in a NFW potential (i.e., discards
$c_{p}=1$ for $r_{s}=r_{sp}$). Repeating the exercise with the other GCs
renders the lower panel in Fig. 16, which is equivalent to Figs. 4 and 5 in
the main text. Again, the star symbols representing the true potentials fall
well within the allowed region in the diagnostic plot, and have harder times
to be consistent with NFW potentials. Thus, the tools employed in the paper
give consistent results when applied to GCs, thought to be self-gravitating
structures and therefore having known gravitational potentials.

Figure 16:  Sanity check employing globular clusters (from Miocchi et al.,
2013), thought to be self-gravitating systems. Top panel: fit to one of them,
NGC 288. Middle panel: diagnostic diagram for NGC 288 inferred from the above
fit. Because it is self-gravitating, the parameters of the stellar
distribution (indicated by the star symbol) are those of the true
gravitational potential ($r_{sp}$ and $c_{p}$), and they appear in the allowed
region of the diagnostic diagram. Bottom panel: diagnostic diagrams combining
all GCs together. As in the middle panel, the star symbols indicated the true
potentials, and they all fall well within the allowed region in the diagnostic
plot.