# Can cuspy dark matter dominated halos hold cored stellar mass distributions?

Jorge Sánchez Almeida Instituto de Astrofísica de Canarias, La Laguna,
Tenerife, E-38200, Spain  Departamento de Astrofísica, Universidad de La
Laguna  Angel R. Plastino 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  Ignacio Trujillo
Instituto de Astrofísica de Canarias, La Laguna, Tenerife, E-38200, Spain
Departamento de Astrofísica, Universidad de La Laguna

###### Abstract

According to the current concordance cosmological model, the dark matter (DM)
particles are collision-less and produce self-gravitating structures with a
central cusp which, generally, is not observed. The observed density tends to
a central plateau or core, explained within the cosmological model through the
gravitational feedback of baryons on DM. This mechanism becomes inefficient
when decreasing the galaxy stellar mass so that in the low-mass regime
($M_{\star}\ll 10^{6}\,{\rm M_{\odot}}$) the energy provided by the baryons is
insufficient to modify cusps into cores. Thus, if cores exist in these
galaxies they have to reflect departures from the collision-less nature of DM.
Measuring the DM mass distribution in these faint galaxies is extremely
challenging, however, their stellar mass distribution can be characterized
through deep photometry. Here we provide a way of using only the stellar mass
distribution to constrain the underlying DM distribution. The so-called
Eddington inversion method allows us to discard pairs of stellar distributions
and DM potentials requiring (unphysical) negative distribution functions in
the phase space. In particular, cored stellar density profiles are
incompatible with the Navarro, Frenk, and White (NFW) potential expected from
collision-less DM if the velocity distribution is isotropic and the system
spherically symmetric. Through a case-by-case analysis, we are able to relax
these assumptions to consider anisotropic velocity distributions and systems
which do not have exact cores. In general, stellar distributions with radially
biased orbits are difficult to reconcile with NFW-like potentials, and cores
in the baryon distribution tend to require cores in the DM distribution.

Cold dark matter (265), Dwarf galaxies (416), Galaxy mass distribution (606),
Navarro-Frenk-White profile (1091), Galaxy dark matter halos (1880), Low
surface brightness galaxies (940), Theoretical techniques (2093)

††software: numpy (Harris et al., 2020), scipy (Virtanen et al., 2020)

##  1 Introduction

The current concordance cosmological model assumes the dark matter particles
to be cold and collision-less (e.g., White & Rees, 1978; Blumenthal et al.,
1984; Davis et al., 1985; Smoot et al., 1992; Peebles, 2021; Bechtol et al.,
2022). Thus, the cold dark matter (CDM) particles interact with themselves and
with the baryons through gravitational forces only. Given the initial
conditions set by the cosmological model, the CDM particles evolve under their
own gravity to collapse into halos with cusps (e.g., Cen, 2014; Brown et al.,
2020), i.e., where the density is represented by the iconic NFW profile
(Navarro, Frenk, and White 1997) that grows boundlessly when approaching the
center of the gravitational potential. This prediction contrasts with the fact
that the observed dark matter (DM) haloes often show cores, i.e., their
density tend to be constant as one approaches the center (e.g., Weinberg et
al., 2015; Del Popolo & Le Delliou, 2017; Bullock & Boylan-Kolchin, 2017).
This apparent contradiction is solved within the current CDM paradigm because
the baryon dynamics modifies the global gravitational potential also affecting
the DM distribution and transforming the cusps into cores (Davis et al., 1992;
Governato et al., 2010; Di Cintio et al., 2014b). This mechanism of feedback
of baryons onto DM becomes inefficient when decreasing the galaxy mass,
because the halo-to-stellar mass ratio increases with decreasing mass (e.g.,
Behroozi et al., 2013), reaching a point where the energy provided by star
formation is simply not enough to modify the cusp of the CDM haloes (e.g.,
Peñarrubia et al., 2012; Oñorbe et al., 2015). The larger stellar mass unable
to modify the inner slope of the DM profile is somewhat model dependent (e.g.,
Read et al., 2016), but it roughly corresponds to stellar masses
$M_{\star}<10^{6}\,{\rm M}_{\odot}$ or halo masses $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). Thus, if galaxies
with $M_{h}\ll 10^{10}\,{\rm M_{\odot}}$ show DM cores, they are not due to
baryon feedback but have to reflect the nature of DM: whether it is fuzzy,
self-interacting, warm, or else (e.g., Dodelson & Widrow, 1994; Hu et al.,
2000; Spergel & Steinhardt, 2000; Bechtol et al., 2022).

At these low masses, discerning observationally whether the DM halos have
cores is extremely challenging, if not impossible. DM measurements require
high spectral resolution spectroscopy to infer dynamical masses (whether
optical, infrared, or radio wavelengths are used). The light is spread into
small wavelength bins and so getting high signal-to-noise ratios is expensive
observationally. On the contrary, stellar mass determinations depend on broad-
band photometry, which is orders of magnitude faster than spectroscopy. Thus,
measuring the baryon mass distribution in these low-mass objects is doable
(e.g., Trujillo et al., 2021) and, interestingly, low-mass galaxies tend to
show cores in their stellar mass distribution (e.g., Moskowitz & Walker, 2020;
Carlsten et al., 2021). Since low-mass galaxies are often extremely DM
dominated systems, one could naively think that the cores observed in stars
just reflect the underlying DM mass distribution. If this conjecture turned
out to be correct, it would provide a unique channel to study DM in low-mass
galaxies, in a regime particularly informative to reveal the nature of DM
(e.g., Weinberg et al., 2015; Bullock & Boylan-Kolchin, 2017; Del Popolo & Le
Delliou, 2017). Thus, the question arises as to whether the cores in the
stellar mass distribution of DM dominated systems trace or not cores in the DM
distribution.

The purpose of this work is bringing up the question in the title to try to
give an answer in fairly broad terms. Thus, we show to be unlikely (although
not impossible) that DM dominated systems with a central cusp have a stellar
profile with a central core. Therefore, our work provides a gateway to
investigate the inner shape of the DM distribution in ultra-low mass galaxies
using only their starlight.

We address the question using the so-called Eddington inversion method
(Eddington, 1916; Binney & Tremaine, 2008; Lacroix et al., 2018). Simply put,
it provides the distribution function (DF) in the phase space $f$
corresponding to a stellar mass density distribution $\rho$ immersed in a
gravitational potential $\Phi$. Given two arbitrary $\rho$ and $\Phi$, there
is no guarantee that $f>0$ everywhere, which is the absolutely minimum
requirement for $\rho$ and $\Phi$ to be physically consistent. In this paper,
we study the $f$ resulting from different combinations of $\rho$ (tracing the
stars) and $\Phi$ (dictated only by the DM in ultra-low mass galaxies). We
will show that unless the potential $\Phi$ is created by a mass distribution
with a core, cored $\rho$s often give nonphysical $f<0$. The computations in
the paper neglect the contribution of the baryons to the overall potential,
which we regard as a reasonable working hypothesis for the galaxies of
interest. Thus, the gas in the ultra-low mass galaxies is not treated
explicitly in the paper, but should play only a minor role in the analysis
since it interacts with the stars only through its contribution to the
gravitational potential. Therefore, as soon as the gas mass is much smaller
than the total mass of the system, its presence can be neglected.

The paper is organized as follows: Sect. 2 puts forward the Eddington
inversion method together with the main equations used in our analysis. The
more lengthly derivations are separated in Appendixes A to E. Unphysical pairs
$\rho$ – $\Phi$ yielding $f<0$ somewhere are analyzed in Sect. 3. Among which
one finds the prototypical cored $\rho$ immersed a NFW potential with
isotropic velocities. Examples and particular cases are worked out in Sect. 4
to conclude that most often the cores in baryons trace DM cores in DM
dominated self-gravitating systems. These results and their practical
application to real galaxies are analyzed in Sect. 5, including the effect of
relaxing assumptions like spherical symmetry. Table 1 lists consistent and
inconsistent combinations of $\rho$ and $\Phi$ resulting from our analysis. In
what follows, we use the terms baryons, stars, or particles indistinctly to
refer to the component of the gravitationally bound system that provides the
density $\rho$. Moreover, in the context of this paper, the term low-mass
galaxy is used to describe galaxies where the potential is approximately set
by the DM because the gravity produced by the baryons can be neglected.

##  2 The Eddington inversion method in our context

This section provides a summary of the Eddington inversion method, and so, of
the expressions used in Sects. 3 and 4 to study whether cored baryon density
distributions happen to be inconsistent with the gravitational potential
created by CDM alone. We closely follow the approach and terminology by Binney
& Tremaine (2008, Sect. 4.3), but there are several alternative references on
the subject (e.g., Ciotti & Pellegrini, 1992; Ciotti, 1996; Lacroix et al.,
2018). The main assumptions made when using the Eddington inversion method are
(Binney & Tremaine, 2008, Sect. 4): (1) the gravitational potential is smooth,
(2) the trace particles (e.g., stars) have lifetimes larger than the crossing
time, (3) the trace particles are collision-less, (4) the system is
spherically symmetric, and (5) the system is described by a steady-state DF in
the phase space. We take these assumptions as working hypotheses, which may
not be fulfilled by particular objects but which may be good enough to
describe large populations. For example, after a major merger the steady-state
may require a few Gyr to be recovered (e.g., Lotz et al., 2008), however, most
galaxies only have a few such events during their lifetimes concentrated early
on, therefore, many galaxies should be in a quasi-steady state today. Based on
these premises, we first consider particle systems with an isotropic velocity
distribution. Sect. 2.1 explains how to use the Eddington inversion method to
recover the phase space DF from the three first spatial derivatives of the
baryon density and of the gravitational potential. The general expressions are
particularized to specific mass distributions and gravitational potentials in
Appendixes A.1 and A.2. Section 2.2 relaxes the assumption on the velocity
isotropy, working out the expression of the DF for the Osipkov-Merritt
velocity anisotropy model. Other anisotropic velocity models are considered
too. Even if contrived from a physical stand point, any gravitational
potential is consistent with any density if the particles are arranged in
perfectly circular orbits. The mixing model in Sect. 2.3 describes the linear
superposition of such a DF with circular orbits plus another DF with an
isotropic velocity distribution. Finally, Sect. 2.4 treats the case of
constant velocity anisotropy.

These physical systems and the corresponding DFs were chosen for simplicity,
because they provide clear-cut constraints on the potential with relatively
simple arguments. There are extensions of the Eddington inversion method for
other more general DFs that in principle could be used for similar diagnostics
(e.g., Lynden-Bell, 1962; Dejonghe, 1987; Cuddeford, 1991; Strigari et al.,
2017), but their study remains to be carried out, a task that requires
specific follow-up work (Sect. 5).

###  2.1 Systems with isotropic velocity distribution

For spherically symmetric systems of particles with isotropic velocity
distribution, the phase-space DF $f(\epsilon)$ depends only on the particle
energy $\epsilon$. Then, the space density $\rho(r)$ turns out to be (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)

Here $\epsilon=\Psi-\frac{1}{2}v^{2}$ is the relative energy (per unit mass)
of a particle, and $\Psi(r)=\Phi_{0}-\Phi(r)$ is the relative potential
energy, where $\Phi(r)$ is the gravitational potential energy and $\Phi_{0}$
is the gravitational potential energy evaluated 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$,

$\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 leads to Eddington’s celebrated equation (e.g.,
Binney & Tremaine 2008, Eq. [4.46]) for the phase-space DF $f(\epsilon)$ in
terms of the spatial density $\rho(r)$,

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

Integrating by parts twice,

$f(\epsilon)=\frac{1}{\sqrt{2}\pi^{2}}\,\left[\frac{1}{2\sqrt{\epsilon}}\,\left(\frac{d\rho}{d\Psi}\right)_{\Psi=0}\,+\,\sqrt{\epsilon}\,\left(\frac{d^{2}\rho}{d\Psi^{2}}\right)_{\Psi=0}\,+\,\int_{0}^{\epsilon}\,\frac{d^{3}\rho}{d\Psi^{3}}\,\sqrt{\epsilon-\Psi}\,d\Psi\right].$
(4)

The derivatives at the boundary, $\left(d\rho/d\Psi\right)_{\Psi=0}$ and
$\left(d^{2}\rho/d\Psi^{2}\right)_{\Psi=0}$, are in practice zero (see
Appendix B), therefore,

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

To evaluate numerically the integral appearing in Eq. (5), it is convenient to
change the integration variable from $\Psi$ to $r$, because only $\rho(r)$ and
$\Psi(r)$ are known explicitly. To use $r$ as integration variable, we need to
express the derivatives of $\rho$ with respect to $\Psi$ in terms of the
derivatives of $\rho$ and $\Psi$ with respect to $r$, i.e.,

$\frac{d\rho}{d\Psi}\,=\,\frac{d\rho/dr}{d\Psi/dr},$ (6)
$\frac{d^{2}\rho}{d\Psi^{2}}\,=\,\left(\frac{d\Psi}{dr}\right)^{-3}\,\left[\left(\frac{d^{2}\rho}{dr^{2}}\right)\left(\frac{d\Psi}{dr}\right)-\left(\frac{d\rho}{dr}\right)\left(\frac{d^{2}\Psi}{dr^{2}}\right)\right],$
(7)
$\frac{d^{3}\rho}{d\Psi^{3}}\,=\,\left(\frac{d^{3}\rho}{dr^{3}}\right)\left(\frac{d\Psi}{dr}\right)^{-3}-3\left(\frac{d^{2}\rho}{dr^{2}}\right)\left(\frac{d^{2}\Psi}{dr^{2}}\right)\left(\frac{d\Psi}{dr}\right)^{-4}-\left(\frac{d\rho}{dr}\right)\left(\frac{d^{3}\Psi}{dr^{3}}\right)\left(\frac{d\Psi}{dr}\right)^{-4}$
(8)
$+3\left(\frac{d\rho}{dr}\right)\left(\frac{d^{2}\Psi}{dr^{2}}\right)^{2}\left(\frac{d\Psi}{dr}\right)^{-5}.$

We now change the integration variable in the integral appearing in Eq. (5),

$\int_{0}^{\epsilon}\,\frac{d^{3}\rho}{d\Psi^{3}}\,\sqrt{\epsilon-\Psi}\,d\Psi\,=\,\int_{r_{m}}^{R}\,\frac{d\Psi}{dr}\frac{d^{3}\rho}{d\Psi^{3}}\,\sqrt{\epsilon-\Psi}\,dr\,=\,-\,\int_{R}^{r_{m}}\,\frac{d\Psi}{dr}\frac{d^{3}\rho}{d\Psi^{3}}\,\sqrt{\epsilon-\Psi}\,dr,$
(9)

where $R$ is the value of $r$ such that $\Psi(R)=\epsilon$, and $r_{m}$ is the
maximum value of $r$, corresponding to the outer edge of the system. When the
system has infinite spatial extent $r_{m}\to\infty$. Replacing the expression
(8) for $d^{3}\rho/d\Psi^{3}$ into the integral (9),

$\displaystyle\int_{0}^{\epsilon}\,\frac{d^{3}\rho}{d\Psi^{3}}\,\sqrt{\epsilon-\Psi}\,d\Psi$
$\displaystyle=$
$\displaystyle\int_{R}^{r_{m}}\,\left[-\left(\frac{d^{3}\rho}{dr^{3}}\right)\left(\frac{d\Psi}{dr}\right)^{-2}+3\left(\frac{d^{2}\rho}{dr^{2}}\right)\left(\frac{d^{2}\Psi}{dr^{2}}\right)\left(\frac{d\Psi}{dr}\right)^{-3}\right.$
(10) $\displaystyle+$
$\displaystyle\left.\left(\frac{d\rho}{dr}\right)\left(\frac{d^{3}\Psi}{dr^{3}}\right)\left(\frac{d\Psi}{dr}\right)^{-3}-3\left(\frac{d\rho}{dr}\right)\left(\frac{d^{2}\Psi}{dr^{2}}\right)^{2}\left(\frac{d\Psi}{dr}\right)^{-4}\right]\,\sqrt{\epsilon-\Psi}\,dr.$
(11)

In short, according to Eqs. (5) and (11), the DF $f(\epsilon)$ corresponding
to a density $\rho(r)$ in a potential $\Phi(r)$ can be deduced from the first
three derivatives of $\rho(r)$ and $\Psi(r)$ ($=\Phi_{0}-\Phi$). Appendix A
works them out in various practical cases involving polytropic $\rho$ and NFW
densities and potentials. Examples of $\rho$ and $f(\epsilon)$ will be shown
in Sect. 4, Figs. 1 – 3, 5, and 6.

###  2.2 Systems with anisotropic velocity distribution: the Osipkov-Merritt
model

The systems described in Sect. 2.1 have DFs depending on the particle energy
only, which holds when the dispersion of velocities in the three independent
spatial directions is the same. In terms of the so-called anisotropy
parameter, these systems have $\beta(r)=0$, with

$\beta(r)\,=\,1\,-\frac{\sigma_{\theta}^{2}+\sigma_{\phi}^{2}}{2\sigma_{r}^{2}},$
(12)

where $\sigma_{r}$ is the radial velocity dispersion, and $\sigma_{\theta}$
and $\sigma_{\phi}$ are the tangential velocity dispersions in spherical
coordinates. The velocity isotropy requirement can be relaxed assuming $f$ to
depend not only on $\epsilon$ but also on the modulus of the angular momentum
$L$. This is done in the Osipkov-Merritt model, which assumes a radial
dependence of the anisotropy given by,

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

where the anisotropy radious $r_{b}$ sets the spatial scale of the changes in
anisotropy. For $r\ll r_{b}$ the velocity distribution is isotropic, while for
$r\gg r_{b}$ it is fully anisotropic with $\beta\to 1$ and the orbits becoming
mostly radial ($\sigma_{\theta}^{2}+\sigma_{\phi}^{2}\ll 2\sigma_{r}^{2}$).
The assumption on $\beta(r)$ in the Osipkov-Merritt model (Eq. [13]) may look
artificial driven by analytical simplicity, but it is not quite so. This type
of radial dependence of the anisotropy parameter is obtained in cosmological
numerical simulations of galaxy formation in the low mass end of the mass
spectrum (e.g., El-Badry et al., 2017; Orkney et al., 2023). In these
simulations, the stars tend to have isotropic orbits in the center of the
potential that turn into radial orbits in the outskirts. In the same numerical
simulations, the DM haloes are more isotropic all over (discussed further in
Sect. 5).

Following Binney & Tremaine (2008), the phase-space DF of the Osipkov-Merritt
model depends on the particle position and velocity through the quantity,

$Q\,=\,\epsilon-\frac{L^{2}}{2r_{b}^{2}}.$ (14)

It is convenient to define

$\rho_{\rm OM}(r)\,=\,\left(1+\frac{r^{2}}{r_{b}^{2}}\right)\,\rho(r).$ (15)

The connection between the mass density and the phase space density can be
expressed, in terms of $\rho_{\rm OM}$, in a way similar to the one
corresponding to isotropic systems. Indeed, one has,

$\rho_{\rm OM}(r)=4\pi\sqrt{2}\,\int_{0}^{\Psi(r)}\,f_{\rm
OM}(Q)\sqrt{\Psi(r)-Q}\,dQ,$ (16) $\frac{d\rho_{\rm
OM}}{d\Psi}\,=\,2\sqrt{2}\pi\,\int_{0}^{\Psi}\,\frac{f_{\rm
OM}(Q)}{\sqrt{\Psi-Q}}\,dQ,$ (17)

and

$f_{\rm
OM}(Q)=\frac{1}{2\sqrt{2}\pi^{2}}\,\frac{d}{dQ}\,\int_{0}^{Q}\,\frac{d\rho_{\rm
OM}}{d\Psi}\,\frac{d\Psi}{\sqrt{Q-\Psi}}.$ (18)

Therefore, expressions (5) – (11) also hold in this case replacing $\rho$ with
$\rho_{\rm OM}$ and $\epsilon$ with $Q$.

###  2.3 Systems with anisotropic velocity distribution: the mixing model

A particle system having only circular orbits has always radial velocity
equals zero, and so, $\sigma_{r}=0$, which leads to $\beta=-\infty$
everywhere. A system with such an extreme velocity anisotropy can reproduce
any pair potential – density with a DF, denoted here as $f_{c}$, guaranteed to
be positive everywhere (see Appendix C). A fairly general system with
anisotropic velocity distribution can be constructed as a linear superposition
of a system with circular orbits $f_{c}$ and a system with isotropic velocity
distribution $f_{i}$ (Binney & Tremaine, 2008), so that

$f({\bf r},{\bf
v})=\mu\,f_{i}\left[\Psi(r)-v^{2}/2\right]+(1-\mu)\,f_{c}\left[{\bf
r},v_{r},v_{\theta},v_{\phi}\right],$ (19)

with ${\bf r}$ and ${\bf v}$ the position and velocity in the 6D phase space
and $\mu$ parameterizing the mixing fraction ($0\leq\mu\leq 1$). The symbols
$v_{r},v_{\theta},v_{\phi}$ represent the three coordinates of the velocity
vector in a reference system where $v_{r}$ is the component in the radial
direction set by ${\bf r}$. Equation (19) explicitly shows that $f_{i}$
depends on the velocity through $v^{2}=v_{r}^{2}+v_{\theta}^{2}+v_{\phi}^{2}$,
a property used in Sect. 3 to discuss the feasibility of DFs from the mixing
model. In this case, the anisotropy parameter at a fixed radius is

$\beta=-\frac{1-\mu}{\mu}\,\frac{\sigma_{\theta}^{2}+\sigma_{\phi}^{2}|_{c}}{\sigma_{\theta}^{2}+\sigma_{\phi}^{2}|_{i}}\leq
0,$ (20)

where $\sigma_{\theta}$ and $\sigma_{\phi}$ stand for the velocity dispersion
in the two tangential coordinates and $|_{i}$ and $|_{c}$ point out the
isotropic and the circular velocity DF, respectively. The resulting orbits are
between circularly biased and isotropic, but never radially biased.

###  2.4 Systems with anisotropic velocity distribution: constant velocity
anisotropy

An extension of the above Eddington formalism deals with anisotropic velocity
distributions of constant $\beta$. One starts from the DF,

$f(\epsilon,L)=L^{-2\beta}f_{\epsilon}(\epsilon),$ (21)

which represents a leading order approximation for a wide class of DFs having
the anisotropy parameter $\beta$ constant (An & Evans, 2006; Binney &
Tremaine, 2008). In these systems the DF depends not only on the relative
energy $\epsilon$ but also on the modulus of the angular momentum $L$. Under
this assumption, the mass volume density can be written as (Binney & Tremaine
2008, Eq. [4.66]),

$r^{2\beta}\rho(r)=\kappa_{\beta}\,\int_{0}^{\Psi}\frac{f_{\epsilon}(\epsilon)}{(\Psi-\epsilon)^{\beta-1/2}}\,d\epsilon,$
(22)

where $\kappa_{\beta}$ is a positive numerical value independent of the radius
$r$. Note that this equation is formally quite similar to Eq. (1) provided
$\beta<1/2$, and so will be used in Sect. 3 to point out the inconsistency of
a large number of densities and potentials in way that parallels the isotropic
case. This DF is also closely connected the so-called cusp slope-central
anisotropy theorem by An & Evans (2006), which links the inner slope of a
density profile with the velocity anisotropy. It is examined in our context in
Appendix D.

##  3 Positivity of the phase-space distribution function

Positivity is the basic requirement for any physically sensible phase-space
distribution. Given a relative potential function $\Psi$ and a mass density
profile $\rho$, it is not guaranteed that the phase-space distribution yielded
by the Eddington inversion method is positive everywhere in the phase space. A
negative distribution function implies that the assumptions made when applying
Eddington’s method are physically inconsistent: there is no phase-space DF
that can reproduce the mass density $\rho$ under the assumed potential $\Psi$.
We use this idea here and in Sect. 4 to analyze the consistency of several
combinations of $\Psi$ and $\rho$ that may be of practical importance.

Requiring $f$ to be non-negative constrains the properties of the centers of
low-mass galaxies in fairly general terms. Equation (2) leads to a sufficient
condition for the physical incompatibility between $\rho(r)$ and $\Psi(r)$
(e.g., Lacroix et al., 2018). If, for a given $\Psi$ (and, consequently, a
given $r$), $d\rho/d\Psi$ vanishes, then, it follows from Eq. (2) that the
phase-space density $f(\epsilon)$ yielded by the Eddington method must reach
negative values. (For the integral [2] to be zero with $f(\epsilon)\not=0$,
$f(\epsilon)<0$ somewhere within the interval $0\leq\epsilon\leq\Psi$.) Thus,
if $d\rho/d\Psi=0$ somewhere, then no isotropic distribution is compatible
with the given $\rho(r)$ and $\Psi(r)$. Taking into account the relation,

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

it follows that a cored mass density, defined as having

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

is inconsistent with a NFW background potential, which has

$\lim_{r\to 0}\frac{d\Psi}{dr}=-\frac{V_{c}}{2\,r_{s}^{2}}\neq 0;$ (25)

see Eq. (A8), with the constants $V_{c}$ and $r_{s}$ defined in Appendix A.1.

The condition for $f>0$ derived above has a twist when the requirement of
having a core (Eq. [24]) is somewhat relaxed. Consider a power law baryon
density profile, $\rho\propto r^{-\alpha}$, with $\alpha=0$ for a cored
profile. Consider also a power law for the density profile generating the
potential, $\rho_{p}\propto r^{-\alpha_{p}}$, with $\alpha_{p}=1$ for a NFW
profile. The relative potential $\Psi$ follows from $\rho_{p}$ so that for
$\alpha\not=0$ one finds111The relation is given explicitly in Sect. 4.2, Eqs.
(34) and (36).,

$\frac{d\rho/dr}{d\Psi/dr}\simeq A\,\alpha\,\,r^{-(2+\alpha-\alpha_{p})},$
(26)

with $A>0$ provided $0<\alpha_{p}<3$. (The case $\alpha=0$ is controlled by a
second term to be added to the RHS - right hand side - of Eq. [26], and is
treated in Appendix E.) Thus the inconsistency between baryons and potential
when $r\rightarrow 0$ disappears when $\alpha>0$ since
$(d\rho/dr)\big{/}(d\Psi/dr)\rightarrow\infty$ when $r\rightarrow 0$. Note,
however, that $f$ may still be negative somewhere else with $r\not=0$ even
when $\alpha>0$, as we will show to be often the case (Sect. 4.2). Note also
that all profiles with $\alpha<0$ (i.e., density decreasing toward the center)
are discarded for any potential since the derivative Eq. (26) is either zero
or negative when $r\rightarrow 0$.

The above results hold for systems where the velocity anisotropy is zero,
however, they can be extended to others more general anisotropic systems. The
Osipkov-Merritt model, which has an anisotropy parameter given by Eq. (13),
follows a relation for the DF (Eq. [17]) formally identical to Eq. (2).
Provided the density profile has a core (i.e., provided it follows Eq. [24]),

$\lim_{r\to 0}\frac{d\rho_{OM}}{dr}=\lim_{r\to
0}\,\left(2r\rho/r_{b}^{2}+\left[1+r^{2}/r_{b}^{2}\right]\frac{d\rho}{dr}\right)=0,$
(27)

which implies that cored density profiles are incompatible with a NFW
potential even when the velocity is anisotropic following an Osipkov-Merritt
model. As we stress in Sect. 2.2, this model for the radial variation of the
anisotropy parameter is not as contrived as one may think since it is roughly
followed by the low-mass model galaxies resulting from cosmological numerical
simulations of galaxy formation.

The constraints posed above happen to be a consequence of a more-general cusp
slope-central anisotropy theorem by An & Evans (2006). These authors showed
that for systems with constant velocity anisotropy $\beta$ (i.e., those
described in Sect. 2.4), the need for the DF to be positive provides a
constraint on the inner slope of the density profile $\alpha$ (i.e.,
$\rho\propto r^{-\alpha}{\rm~{}when~{}}r\rightarrow 0$),

$\alpha\geq 2\beta.$ (28)

This holds independently of the gravitational potential $\Psi$. As we show in
Appendix D, when this is combined with the constraint in Eq. (25) set by
having a NFW background potential, it leads to

$\alpha>2\beta.$ (29)

The inequality (29) has a number of implications: (1) cores ($\alpha=0$) are
inconsistent with isotropic velocities ($\beta=0$), as we have shown already,
(2) cores are inconsistent with radially biased orbits (i.e., only $\beta<0$
is allowed for $\alpha=0$), (3) radially biased orbits ($\beta>0$) require
cuspy baryon density profiles ($\alpha>0$), and (3) circular orbits do not
impose any restriction on the inner slope $\alpha$ since their
$\beta=-\infty$. Actually, it is already known that strongly tangentially
biased orbits can reconcile a cored stellar density profile with a cuspy CDM-
like background potential (Breddels & Helmi, 2013, a Schuster-Plummer density
profile in a Hernquist potential).

The constraint on galaxies having radially biased orbits ($\beta>0$) is
particularly important from a practical point of view since these orbits seem
to be the natural outcome of the formation of dwarf galaxies in $\Lambda$CDM
cosmological numerical simulations: see, e.g., El-Badry et al. (2017, Fig. 2)
and Orkney et al. (2023, Fig. 5). Moreover, even if the uncertainties are
large, values of $\beta\gtrsim 0$ are also observed among the DM dominated
satellites of the MW (e.g., Łokas, 2009; Massari et al., 2018, 2020; Read et
al., 2019; Leung et al., 2021; Kowalczyk & Łokas, 2022)

There is also a fairly general family of tangentially biased DFs that can be
discarded right away. It is described by the mixing model (Sect. 2.3) and
covers the whole range of tangentially biased anisotropies from
$\beta=-\infty$ to 0. The mixing model in Sect. 2.3 combines circular orbit
DFs ($f_{c}$ with $\beta=-\infty$) and isotropic velocity DFs ($f_{i}$ with
$\beta=0$) to produce tangentially biased DFs with $\beta<0$ (Eq. [20]). One
may naively think that the always positive $f_{c}$ may compensate $f_{i}<0$ to
yield a positive physically sensible DF $f=\mu\,f_{i}+(1-\mu)\,f_{c}$ (Eq.
[19]). However, all linear combinations can be discarded for any $\mu\not=0$
if $f_{i}<0$ somewhere. The argument goes as follows: assume that $f_{i}<0$ at
${\bf r}={\bf r_{1}}$ and ${\bf v}={\bf v_{1}}=(v_{r1},v_{\theta 1},v_{\phi
1})$ (see the dependencies of the DF on position ${\bf r}$ and velocity ${\bf
v}$ in Eq. [19]). Then $f_{i}$ is also $<0$ at ${\bf r_{1}}$ and ${\bf
v_{2}}=(v_{r2},0,0)$ provided $v^{2}_{r2}=v_{r1}^{2}+v_{\theta 1}^{2}+v_{\phi
1}^{2}$, since $f_{i}$ depends on ${\bf v}$ only through its modulus $v$.
However, $f_{c}({\bf r_{1}},v_{r2},0,0)=0$ because, by definition, $f_{c}$
only represents circular orbits that must have $v_{r}=0$. Thus, Eq. (19) shows
that $f({\bf r_{1}},{\bf v_{2}})<0$ and thus unphysical, with the only
possible workaround of $\mu=0$, and so, of all orbits being circular.

On the basis of the above arguments, ultra-low mass galaxies for which the
stellar mass distribution is well fitted with cored density profiles are
dynamically incompatible with a NFW profile for the dark-mass component, at
least if one assumes that the phase-space DF of the stellar component depends
only on the stellar energy, or that is described by a Osipkov-Merritt model,
or is anisotropic with radially biased orbits, or anisotropic with
tangentially biased orbits following the mixing model. These arguments cannot
rule out a NFW background potential if other types of anisotropic phase-space
distribution for the stars are assumed. For instance, a cored stellar profile
may be compatible with a star distribution having a constant tangentially
biased anisotropy (that is, having a constant and negative $\beta$). To make
the range of compatibilities more clear, Table 1 lists pairs of densities and
potentials together with whether they are consistent or inconsistent.

We note an important property of the consistency of a pair $\rho$ – $\Psi$
based on whether $f(\epsilon)\geq 0$ $\forall\epsilon$. If a particular pair
is consistent or inconsistent, then any global factor affecting the density
profile will not modify this character since $f(\epsilon)$ scales linearly
with a multiplicative factor in $\rho$ (see Eqs. [5] and [11]). Thus, any of
the inconsistencies brought out here hold true independently of the (typically
unknown) mass ratio between the stars and the DM halo creating the potential.

##  4 Numerical results

This section illustrates with specific examples the general results put
forward in Sect. 3, analyzes the behavior outside the core of the system, and
deals with profiles with shapes more complex than the ones considered in Sect.
3. We check whether the DF resulting from particular pairs becomes negative at
some point, which would discard the combination. $\beta=0$ is assumed, so the
DF $f(\epsilon)$ follows from Eqs. (5) and (11). Equation (5) is integrated
numerically for every $\epsilon$ applying a Simpson’s rule. The radial
derivatives of $\rho$ and $\Psi$ in Eq. (11) are computed analytically
whenever possible using the equations in Appendix A (Sect. 4.1). Otherwise we
compute them numerically (Sect. 4.2).

Figure 1:  Cored baryon distribution in a NFW gravitational potential. Left
panel: a polytrope of order $m=5$ gives $\rho(r)$, with the central density
$\rho(0)=10^{5}\,{\rm M_{\odot}\,kpc^{-3}}$ and a core radius $r_{0}=1.4\,{\rm
kpc}$, chosen so that $M_{\star}\sim 10^{6}\,{\rm M_{\odot}}$ (the red solid
line). The NFW density profile that defines the overall gravitational
potential through the Poisson equation (the gray line) has
$\rho_{s}/\rho(0)=10$ and $r_{s}/r_{0}=4$, which provides $M_{\star}/M_{\rm
p}\simeq 1.4\times 10^{-4}$. The bullet symbol points out the total density at
the virial radius, assuming a concentration of 30 ($r_{vir}=30\,r_{s}$). Right
panel: baryon DF $f$ needed to match baryon density and potential according to
Eddington’s inversion method which reaches negative values (the dashed line)
implying that this particular combination is unphysical.

###  4.1 Cored density in a NFW potential

The computation of $f(\epsilon)$ is straightforward when $\rho$ is a Schuster-
Plummer profile,

$\rho(r)=\frac{\rho(0)}{\left[1+\left(r/r_{0}\right)^{2}\right]^{5/2}},$ (30)

and $\Psi$ is described by a NFW potential (Appendix A), a combination used
here as reference of cored density profile immersed in a CDM-only potential.
The Schuster-Plummer profile is the polytrope of order $m=5$, and was chosen
as reference because it provides a fair representation of the stellar mass
distribution in real dwarf galaxies (e.g., Sánchez Almeida et al., 2021). As
the rest of polytropes, this density profile has a core, therefore, it is not
consistent with the potential derived from the cuspy NFW profile (Sect. 3). An
example is shown in Fig. 1. The parameters that define this polytrope and the
potential have been tuned to represent a realistic galaxy with stellar mass
$M_{\star}\simeq 10^{6}\,{\rm M_{\odot}}$, core radius $r_{0}=1.4~{}{\rm
kpc}$, and total mass around $10^{4}$ times the stellar mass (e.g., Behroozi
et al., 2013; Kormendy & Freeman, 2016). The stars are immersed in the NFW
potential generated by the matter distribution represented as the gray solid
line in the left panel of Fig. 1. (Note that this density fully defines $\Psi$
through Poisson’s equation and independently of the velocity distribution of
the particles creating the potential.) This component completely dominates the
mass and the potential of the system: the mass within the gray solid line,
$M_{p}$, has $M_{p}/M_{\star}\simeq 10^{4}$. As expected, $f<0$ for some
$\epsilon$ signaling that this combination of baryons and potential is
unphysical.

For self gravitating systems, Poisson equation guarantees that
$f(\epsilon)\geq 0$ $\forall\epsilon$. For the stellar density profile shown
in Fig. 1, $f(\epsilon)$ is analytic (Eq. [A19]). We use this fact to check
the numerical integration scheme used to derive $f(\epsilon)$.

###  4.2 Double power law density and potential

A cored density in a NFW potential are inconsistent, as we have showed. Here
we expand the range of shapes to figure out how much the conditions for a core
and a NFW potential can be relaxed and still getting inconsistent results. In
our study, we use a family of density profiles commonly used in the literature
(e.g., Hernquist, 1990; Merritt et al., 2006; Di Cintio et al., 2014a),

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

with $x=r/r_{s}$, that encompasses both the NFW profile
($a=1,b=3,{\rm~{}and~{}}c=1$) and the Schuster-Plummer profile
($a=2,b=5,{\rm~{}and~{}}c=0$) shown in Fig. 1. Actually, for $a=2$, $b=m$, and
$c=0$, $\rho_{abc}$ approximately accounts for the inner region of a polytrope
of index $m$ (e.g., Sánchez Almeida, 2022) which is important in this context
since polytropes describe density profiles of self-gravitating N-body systems
when reaching thermodynamical equilibrium (see, Plastino & Plastino, 1993;
Sánchez Almeida et al., 2020; Sánchez Almeida & Trujillo, 2021). The constants
$r_{s}$ and $\rho_{s}$ in Eq. (31) provide the global scaling for radius and
density, respectively. The parameter $c$ gives the inner logarithmic slope,

$\lim_{r\to 0}\frac{d\log\rho_{abc}}{d\log r}=-c.$ (32)

The three-parameter function $\rho_{abc}$ can be folded into a single
parameter family using $a=2-c$ and $b=5-2c$,

$\rho_{c}(r)=\frac{\rho_{s}}{x^{c}(1+x^{2-c})^{(5-3c)/(2-c)}},$ (33)

which seamlessly scans from Schuster-Plummer to NFW when $c$ goes from 0 to 1
(see Fig. 2).

Figure 2:  Doble power law density profile (Eq. [33]) that goes seamlessly
from a Schuster-Plummer profile to a NFW profile when the inner slope ($-c$)
goes from 0 to -1.

In order to compute the DF of the baryons, one needs the derivatives of the
density profile and the potential (Eqs. [5] and [11]). Using the Poisson
equation for a spherically symmetric system (e.g., An & Zhao, 2013), the
potential and the required derivatives can be obtained in terms of the inner
mass,

$M_{p}(<r)=4\pi\,\int_{0}^{r}\,t^{2}\,\rho_{p}(t)\,dt,$ (34)

so that

$\Psi=\frac{G\,M_{p}(<r)}{r}+4\pi G\,\int_{r}^{\infty}\,t\,\rho_{p}(t)\,dt,$
(35) $\frac{d\Psi}{dr}=-\frac{G\,M_{p}(<r)}{r^{2}},$ (36)
$\frac{d^{2}\Psi}{dr^{2}}=\frac{2GM_{p}(<r)}{r^{3}}-4\pi G\rho_{p},$ (37)

and

$\frac{d^{3}\Psi}{dr^{3}}=-\frac{6GM_{p}(<r)}{r^{4}}+8\pi
G\frac{\rho_{p}}{r}-4\pi G\frac{d\rho_{p}}{dr},$ (38)

which for $\rho_{p}=\rho_{abc}$ can be computed analytically only for certain
values of $a,b$ and $c$ (An & Zhao, 2013).

Employing, Eq. (31) and Eqs. (34) – (38), one can integrate numerically Eq.
(11) to obtain $f(\epsilon)$ via Eq. (5). Using this approach, we have scanned
through a large number of pairs baryon densities and potentials both
characterized by double exponential density profiles but with different
parameters. Unless otherwise stated explicitly, we employ the simplified
version of the density given in Eq. (33). The main results of our numerical
exercise will be discussed next and are also summarized in Table 1. To prevent
confusion during the description, the parameters corresponding to the density
profile that creates the potential are labeled with the subscript $p$ whereas
those of the baryon density do not have any subscript.

  1. 1.

If the baryons have a core (i.e., if $c=0$) then the density generating the
potential must also have a core to be consistent (i.e., $c_{p}=0$). This is
shown by the numerical simulations (a counter-example with $c=0$, $c_{p}=0.05$
producing $f<0$ is shown in Fig. 3), but it also follows analytically from the
study carried out in Appendix E and discussed in the item 6 below. The
graphical summary in Fig. 4 shows that when $c=0$, $c_{p}$ must be zero for
$f$ to be $>0$ everywhere.

Figure 3: Similar to Fig. 1 but this time only with a hint of cusp in the
density that generates the potential: $c=0$ and $c_{p}=0.05$. Note that $f<0$
at large $\epsilon$, implying that density and potential are inconsistent with
each other.

Figure 4:  Summary of the allowed (blue bullet symbols) and forbidden (red
plus symbols) regions of the parameter space when both the baryon density
profile and the density defining the underlying potential follow a profile
described by Eq. (33) and illustrated in Fig. 2. The symbols $c$ and $c_{p}$
represent the inner slope of the density profile and the density generating
the potential, respectively. The ordinate axis is split into two (logarithmic
scale on top and linear scale at the bottom ) to show the whole range of
values of interest. The lines $x=0$ and $y=0$ are shown in light grey. Left
panel: $r_{s}/r_{sp}=1/4$, with the baryons more centrally concentrated that
the potential. Right panel: $r_{s}/r_{sp}=2$, with the baryons more spread out
than the potential. Baryons do not contribute to the overall potential and
their velocity distribution is assumed to be isotropic ($\beta=0$) at all
radii. In both cases, the sampling in $c_{p}$ is $\Delta c_{p}=0.02$ whereas
the sampling in $c$ is $\Delta c=10^{-4}$ in the linear scale and $\Delta\log
c=0.05$ in the logarithmic scale.

  2. 2.

If the baryon core is not perfect ($c\gtrsim 0$; denoted as soft-core in Table
1), then a NFW profile may or may not be compatible with it. Figure 5 shows
examples of incompatible (top panels) and compatible (bottom panels). We have
scanned a range of values for $c$ and $c_{p}$ ($-0.001<c<1$ and
$-0.1<c_{p}<1.1$) with the compatibility summarized in Fig. 4. Roughly
speaking, NFW potentials ($c_{p}=1$) are inconsistent with densities having
$c\lesssim 0.1$.

Figure 5:  Top panels: similar to Fig. 1 but with a hint of cusp in the
density: $c=0.005$ and $c_{p}=1$. The DF $f<0$ but not at the largest
$\epsilon$. This combination of baryon soft-core and NFW potential is still
inconsistent. However, the inconsistency goes away as soon as the inner slope
of the stellar density profile increases, as shown in the bottom panels, where
$c=0.1$ and $c_{p}=1$. The global picture of compatibility – incompatibility
is summarized in Fig. 4. The shape of the profiles defining the baryon
distribution and the potential is given by Eq. (33).

  3. 3.

Any density profile with $c=0$ and $a>2$ is physically irrealizable (see Fig.
6), independently of the potential. The inconsistency remains even with the
potential created by the self-gravity of the density, and means that no
$\beta=0$ DF is able to reproduce $a>2$ profiles. The behavior, summarized in
Fig. 7, is predicted analytically in Appendix E and discussed further in item
6.

Figure 6:  Similar to Fig. 5 with $c=c_{p}=0$ and $a=a_{p}=2.1$. Even if the
two density profiles have the same shape, $f<0$ somewhere (the dashed line).
This behavior for $a>2$ is predicted analytically in Appendix E and discussed
further in Sect 4.2, item 6.  Figure 7:  Similar to Fig. 4 but using
$a=2.5-c$, so that Eq. (39) (or Eq. [E4]) can be tested. Note that $c=c_{p}=0$
is unphysical even though $\rho$ and the density profile producing the
potential are identical.

  4. 4.

A density profile significantly broader than the potential also yield
inconsistent distribution functions. According to Fig. 8,
$r_{s}/r_{sp}\lesssim 2$ for the DF to be non-negative, a constraint that 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. As we mention in Sect. 3, the density
contrast between the density and the density producing the potential
($\rho_{s}/\rho_{sp}$) is irrelevant since it cannot change the sign of $f$.

Figure 8:  Diagnostic plot similar to Fig. 4, except that this time we
represent simulations where both $c$ and $r_{s}/r_{sp}$ are varied. Note how
$r_{s}$ cannot be larger than $\sim 2\,r_{sp}$ for $f$ to remain positive, a
constraint that may be used in real galaxies to set a lower limit to the size
of the DM halo. In these simulations the $abc$ profile shapes of density and
potential are identical, explicitly, $a=a_{p}$, $b=b_{p}$, and $c=c_{p}$, the
three of them varying as for Eq. (33). The sampling in $c$ is $\Delta
c=10^{-4}$ in the linear scale and $\Delta\log c=0.1$ in the logarithmic
scale. The relative radii are shown in a logarithmic scale with a sampling of
$\Delta\log(r_{s}/r_{sp})=0.1$.

  5. 5.

Note that the region where $c>c_{p}$, i.e., where the baryons are more cuspy
than the halo, presents no inconsistency in the summary plots of the Figs. 4
and 7. We bring this fact up because some numerical simulations of ultra-low
mass galaxies seem to show compact stellar concentrations having
$c>c_{p}\simeq 1$ (e.g., Orkney et al., 2021, 2023). These structures are
physically feasible within the logical framework of our work which, among
others, assumes negligible stellar mass ($M_{\star}\ll M_{p}$) and spherical
symmetry.

  6. 6.

For the $abc$ densities that we are considering, the derivative used to
diagnose the positivity of the DF in Sect. 3 (Eq. [23]) turns out to be (Eq.
[E3], Appendix E),

$\frac{d\rho/dr}{d\Psi/dr}\simeq\frac{D}{r^{2+c-c_{p}}}\,\Big{[}c+\frac{b-c}{r_{s}^{a}}\,r^{a}\Big{]},$
(39)

with $D>0$ for $c_{p}<3$. If $c\not=0$, the first term in the RHS of Eq. (39)
dominates the behavior of the ratio when $r\rightarrow 0$. This term is
identical to Eq. (26) with $\alpha=c$ and $\alpha_{p}=c_{p}$, and to allow for
the derivative to differ from zero (and so for the DF to be positive), it only
demands $c>c_{p}-2$. This is a very loose constraint and actually, most of the
forbidden (red) region in Figs. 4 and 7 actually meets this requirement. These
two results, i.e., having $d\rho/d\Psi\not=0$ and but $f<0$ somewhere, are
consistent because the first one is more demanding than the second since it
requires the (weighted) integral of $f$ to be positive (Eq. [2]), which can be
met even when $f<0$ for some values of $\epsilon$ (see the example in the top
panel of Fig. 5). When $c=0$, the second term in the RHS of Eq. (39) rules,
and then the potential and the density profiles would be inconsistent when
$2-c_{p}-a>0$ since the ratio of derivatives goes to zero. For $a=2$, as
expected for polytropes, $c_{p}>0$ is ruled out and the potential must have a
core to be consistent with the core in the density profile. This condition for
$c=0$ is truly restrictive, and is strictly followed by the simulations in
Fig. 4.

##  5 Discussion and Conclusions

Table 1: Summary of the compatibility between baryon density profile ($\rho$)
and potential

Baryons & Potential, Velocity | Consistency | Comments | Section
---|---|---|---
(1) | (2) | (3) | (4)
Core † & NFW ‡, isotropic |  ✗ | Eqs. (24) and (25). $\beta=0\,^{*}$. Fig. 1 | Sect. 3
Power law § & Power law, isotropic |  ✍ |  $\alpha>0$ § ✓ $\alpha<0$ ✗. Eq. (26). $\beta=0$ | Sects. 3, 4.2
Core & Soft-core #, isotropic |  ✗ |  $\beta=0$. Fig. 3. Fig. 4 | Sect. 4.2, App. E
Core & Core, isotropic |  ✍ |  $\beta=0$. $a\leq 2$ ✓$a>2$ ✗. Fig. 6 | Sects. 4.1, 4.2, App. E
Soft-core & NFW, isotropic |  ✍ |  $\beta=0$. Figs. 4, 5. $c\gtrsim 0.1$ ✓$c\lesssim 0.1$ ✗. | Sects. 3, 4.2
Soft-core & Soft-core, isotropic |  ✍ |  $\beta=0$. Figs. 4, 5 | Sects. 3, 4.2
|  |  $r_{s}\gtrsim 2\,r_{sp}$ ✗, $c>c_{p}$ ✓ | Sect. 4.2
Core & NFW, O-M model |  ✗ |  $\beta(\not=0)$ in Eq. (13) | Sect. 3
Core & NFW, radially biased |  ✗ | Constant $\beta$. $\beta>0$ | Sect. 3, App. D
Core & Any, radially biased |  ✗ | Constant $\beta$. $\beta>0$ | Sect. 3, App. D
Power-law & Any, anisotropic |  ✍ | Constant $\beta$. $\alpha>2\beta$ | Sect. 3, App. D
Core & NFW, circular |  ✓ | $\beta=-\infty$ | App. C
Any & Any, circular |  ✓ | $\beta=-\infty$ | App. C
Any & Any, tangentially biased |  ✍ |  $\beta<0$. Eq. (19). ✗$f_{i}<0$ | Sects. 2.3, 3

Note. —

† Core $\equiv\,d\log\rho/d\log r\to 0$ when $r\to 0$.

‡ Navarro, Frenk, and White potential (Eq. [A6]) produced by a NFW profile
(Eq. [A5]).

∗ Velocity anisotropy parameter $\beta$ defined in Eq. (12).

§ $\rho\propto r^{-\alpha}$.

# Soft-cores defined in Eqs. (31) and (33), and illustrated in Fig. 2. Power
laws § are a particular type of those.

(1) Description of the baryon density, the gravitational potential, and the
velocity distribution.

(2) The symbols ✓, ✗, and ✍ stand for compatible, incompatible, and may or may
not, respectively.

(3) Additional comments and keywords.

(4) Section of the text where the combination described in (1) is discussed.

According to the current concordance cosmological model, the DM particles are
collision-less and, evolving under their own gravity, produce self-gravitating
structures that approximately follow the iconic NFW profile with a cusp in its
center (CDM haloes). These cusps ($\rho\propto r^{-1}$) are generally not
observed in galaxies. The total density often presents a central plateau or
core ($\rho\sim{\rm constant}$), which is believed to be produced by the
coupling with baryons through gravity. Star-formation driven outbursts modify
the overall gravitational potential, affecting the CDM distribution too. This
mechanism of baryon feedback becomes inefficient when decreasing the galaxy
stellar mass, reaching a point where the energy provided by baryons is simply
not enough to modify the cusp of the CDM haloes (see, Sect. 1 for references
and details). Despite all uncertainties and model dependencies, this threshold
mass roughly corresponds to isolated galaxies with stellar masses
$<10^{6}\,{\rm M}_{\odot}$ or halo masses $<10^{10}\,{\rm M_{\odot}}$. Thus,
if these ultra-low mass galaxies show cores, they are not due to baryon
feedback processes but have to reflect the nature of DM: whether it is fuzzy,
self-interacting, warm, or any of the other possibilities put forward in the
literature.

Direct measurements of the DM mass distribution in these faint galaxies are
difficult since they require high spectral resolution spectroscopy, which is
observationally extremely challenging. However, there may be a shortcut if the
starlight somehow follows the DM since, even in low-mass low-luminosity
galaxies, deep photometry is doable (e.g., Trujillo et al., 2021). One may
naively think that stars must trace DM in these systems whose potential is
fully dominated by DM. Nevertheless, stars are so weakly coupled with the DM
that can potentially maintain a mass distribution differing from the DM
distribution for longer than the age of the Universe (e.g., Binney & Tremaine,
2008). Thus, in order to use the observable stellar mass distribution as a
proxy for the elusive DM distribution, one has to show that somehow starlight
traces DM in this DM dominated systems. More specifically, we know that low-
mass galaxies often show cores in their stellar mass distribution (Sect. 1).
The question arises as whether this cored baryon distribution is or not
consistent with the DM distribution expected from CDM particles (aka NFW
profile). We address the question using the so-called Eddington inversion
method. Under mildly restrictive assumptions (gravity from baryons negligible,
stationary-state, smooth potential, and spherical symmetry; see, Sect. 2), the
method provides the DF in the phase space $f$ corresponding to a mass density
distribution immersed in a gravitational potential. Given two arbitrary
density and potential, there is no guarantee that $f>0$ everywhere, which is
required for them to be physically consistent.

In this paper, we have studied different combinations of baryon density and
gravitational potential that may help us to discern whether DM profiles in
ultra-low mass galaxies have or not a core. We focus on the consistency of the
various gravitational potentials with baryon density profiles showing a core
(Eq. [24]) or soft-core (Eq. [32], with $c\gtrsim 0$). The main conclusions of
our analysis are summarized in Table 1 and can be expanded as follows:

  * -

Stellar cores in a NFW potential are incompatible provided the velocity
distribution is isotropic ($\beta=0$).

  * -

Stellar cores and potentials stemming from a density with a quasi-core
($c_{p}>0$) are incompatible too. This result holds for isotropic velocities
($\beta=0$).

  * -

As expected for physical consistency, stellar cores and potentials resulting
from cored density profiles 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 provided they follow the Osipkov-Merritt model. Even if
artificial, it approximately describes the global trend expected in ultra-low
mass galaxies, with $\beta\sim 0$ in the center and then increasing outwards
($\beta>0$).

  * -

Stellar cores and NFW potentials are incompatible in systems with radially
biased orbits (${\rm constant~{}}\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 baryion
density and potential, including a cored stellar density in a NFW potential.
These configuration is very artificial, though. Unlikely to happen in real
dwarf galaxies where orbits are expected to be radially biased (see the
discussion below).

  * -

The linear superposition of two DFs is also a DF. Thus, one may think that the
addition of a positive DF for circular orbits may compensate the negative DF
for isotropic orbits to yield a positive physically sensible DF. However, this
is not the case. Independently of the relative weight, the mixing of an
unphysical DF for isotropic velocities ($\beta=0$) with a physically
realizable DF for circular orbits ($\beta=-\infty$) always yields unphysical
DFs (Sect. 2.3).

  * -

We denote as soft-cores those profiles where inner slope is not exactly zero
but close to it ($c\gtrsim 0$). Soft-cores are inconsistent with NFW profiles
when $c\lesssim 0.1$ while they are consistent when $c\gtrsim 0.1$. When the
density profile that characterizes the potential also has a soft core (i.e.,
when $0\leq c_{p}\leq 1$), then the situation is more complicated as shown in,
e.g., Fig. 4. This statements hold for isotropic velocity distributions.

  * -

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

  * -

Positive inner slope in the stellar distribution, where the density grows
outwards, is discarded in every way.

  * -

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 equivalent to 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.

  * -

Pairs of density and potential where the inner slope of the density is larger
than that of the potential ($c\geq c_{p}$) are not inconsistent. This result
refers to isotropic velocities too.

  * -

The above conclusions do not depend on a scaling factor on stellar density
profile, therefore, they do not depend on the (unknown) ratio between the
stellar mass and the total mass of the system.

  * -

The functions used to represent the density and the potential are flexible
enough to describe the central region in any polytrope of arbitrary index $m$
(Eq. [31], with $a=2$, $b=m$, and $c=0$). Polytropes are important in the
context of self-gravitating systems since they describe the density expected
in N-body systems reaching thermodynamical equilibrium (see, Plastino &
Plastino, 1993; Sánchez Almeida et al., 2020). In other words, they portray
the DM density distribution expected if the DM were not collision-less (e.g.,
Sánchez Almeida & Trujillo, 2021).

How useful the above constraints are very much depends on the anisotropy of
the velocity field $\beta$ (Eq. [12]). In general, radially biased ($\beta>0$)
and isotropic ($\beta=0$) orbits are more difficult to reconcile with a cuspy
gravitational potential than tangentially biased orbits ($\beta<0$). The
question arises as what is the anisotropy to be expected in real galaxies.
This issue can be addressed from two complementary directions, namely, what is
the anisotropy observed in the smallest galaxies, and what is the anisotropy
recovered for the smallest galaxies formed in cosmological numerical
simulations. Even if the uncertainties are large because the estimates rely on
measuring velocities of individual stars, the DM dominated satellites of the
Milky Way (MW) tend to have $\beta\gtrsim 0$ (e.g., Łokas, 2009; Massari et
al., 2018, 2020; Read et al., 2019; Leung et al., 2021; Kowalczyk & Łokas,
2022). Note that these objects are not isolated galaxies and their internal
baryon structure may be strongly mediated by the presence of the MW and its
circum-galactic medium through tidal forces, ram-pressure, and starvation
(e.g., Combes, 2004; Benson, 2010; Sánchez Almeida et al., 2017). However, the
observed trend is consistent with numerical simulations. Radial anisotropies
seems to be the natural outcome of the formation of dwarf galaxies in
$\Lambda$CDM cosmological numerical simulations: see, e.g., El-Badry et al.
(2017, Fig. 2) and Orkney et al. (2023, Fig. 5). Moreover, $\beta$ tends to
zero when approaching the center of the gravitational potential, where the
stellar cores may be present and have to be observed. Thus, $\beta\gtrsim 0$
at the centers seems to be a sensible conjecture when interpreting stellar
mass distributions in real galaxies.

One of the seemingly more restrictive assumption leading to the constraints in
Table 1 is the spherical symmetry of the density and potential. As it happens
with the isotropy of the velocity field, the question of whether this is a
good assumption for real ultra-low mass galaxies arises. Actually, the two
issues are closely connected since, in real galaxies, both are set by the
history of star-formation driven by cosmological gas accretion and mergers
(e.g., Dekel et al., 2009; Sánchez Almeida et al., 2014). In general, the
smallest simulated galaxies tend to be rounded, although not perfectly
spherical, with the DM component closer to sphericity (e.g., Bullock, 2002;
Orkney et al., 2023). On the other hand, the observed dwarf isolated galaxies
are triaxial, but with three axes of similar lengths (e.g, Roychowdhury et
al., 2013; Sánchez-Janssen et al., 2016; Putko et al., 2019). In addition to
whether real ultra-low mass galaxies are or not well fitted by spherically
symmetric models, independent theoretical arguments point out that this
assumption is not so critical since the incompatibilities may still hold when
dropped. The extensions of the Eddington inversion method for axi-symmetric
systems (Lynden-Bell, 1962; Binney & Tremaine, 2008) lead to expressions for
the DF similar to Eqs. (2), (17), and (D3). They are expected to lead to
restrictions similar to those worked out in this paper. We are presently
exploring them with promising results. There are also extensions or variants
of the Eddington inversion approach, suitable for other more general
spherically symmetric DFs, that in principle could be used for diagnostics and
would be worth considering (e.g., Dejonghe, 1987; Cuddeford, 1991; Strigari et
al., 2017), but their analysis remains to be carried out.

The constraints in Table 1 result from treating particular cases, each one
with its own peculiarities. The analysis of other cases (e.g., the study of
axi-symmetric systems mentioned above) will enlarge a list which at present
contains only a fraction of the constraints yet to be discovered. In this
sense, our work is only a pathfinder that shows how the traditional Eddington
method can be used to study DM haloes in ultra-low mass galaxies. Given the
observed stellar distribution, the method seriously limits the properties of
the DM halo where it resides. Moreover, its interest probably exceeds the
original scope that motivated the present study, and may be of application to
other astrophysical systems where the stars represent only a minor fraction of
the total mass, for example, the intra-cluster light as tracer of the DM
galaxy cluster potential (e.g., Montes & Trujillo, 2019, 2022).

In short, the question in the title of the paper, Can CDM matter halos hold
cored stellar mass distributions?, has no simple yes or no answer. Instead, we
find it to be unlikely, although not imposible, than cored stellar mass
distributions can be hosted in NFW DM haloes, provided the system is
spherically symmetric. Thus, our work supports the interest of determining
surface brightness profiles of ultra-low mass galaxies to constrain the nature
of DM. This work can be used as a guide to interpret observations so that the
closer the observed galaxies to the hypotheses (spherically symmetry,
stationarity, velocity isotropy, etc.) the more useful the constraints in
Table 1.

Our ultimate goal is applying the mathematical tools developed in this paper
to observed dwarf galaxies with masses low enough to constrain the nature of
DM (Sect. 1). This challenging task still requires several intermediate steps
to be completed. In our roadmap, we would like to test the machinery with the
few local group galaxies for which independent information on the DM halo and
on the stellar distribution is available (e.g., Battaglia & Nipoti, 2022), to
see whether the constraints imposed by the Eddington inversion method and by
the kinematical measurements are consistent. We also need to know what is the
signal-to-noise ratio and the number of targets required to make firm claims.
Having a hundred targets with surface brightness profiles reaching down to 30
mag arcsec-2 seems to be doable (e.g., Carlsten et al., 2021) but, does it
suffice? Finally, we have to carefully select the actual data set of faint
isolated dwarf galaxies. The two requirements are in tension since
intrinsically faint galaxies are nearby and so tend to be satellites, but both
are needed. One obvious possibility is waiting for better data (e.g., Ivezić
et al., 2019; Abbott et al., 2021; Trujillo et al., 2021). Alternatively, one
can also think of studying the ultra-faint dwarfs of the local group (e.g.,
Moskowitz & Walker, 2020) cherry-picking those where the tidal forces and
other environmental effects may be minimal (e.g., with large pericentic
passage) and which truly proceed from low mass progenitors (e.g., Grebel et
al., 2003). Tidal forces change the internal structure of satellites and
reduce their stellar mass content, thus blurring any clear-cut interpretation
of the observed DM distribution in terms of the DM nature, a caveat to keep in
mind if this pathway is chosen. All these works are currently ongoing or
planed.

Thanks are due to Claudio Dalla-Vechhia for insightful discussions during the
early stages of the work, and to Giuseppina Battaglia, Arianna Di cintio, and
Ruben Sánchez-Janssen for references. Thanks are due to Matthew Orkney and
Justin Read for discussions and clarifications on the velocity anisotropy and
mass profile of the galaxies in their simulations. JSA acknowledges financial
support from the Spanish Ministry of Science and Innovation (MICINN), project
PID2019-107408GB-C43 (ESTALLIDOS). 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). JSA also wants
to explicitly thank Angel Luis Platino and the Facultad de Ciencias Económicas
de La Universidad Nacional de La Plata for their hospitality during this
visit. ARP acknowledges support to visit the IAC from the Fundación Jesús
Serra and the IAC under their Visiting Researcher Programme 2020–2022. IT
acknowledges support from the Project PCI2021-122072-2B, financed by
MICIN/AEI/10.13039/501100011033, and the European Union NextGenerationEU/RTRP
and 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
PID2019-107427GB-C32 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.

##  Appendix A Analytic derivatives of the density and the potential

According to Eqs. (5) and (11), the DF $f(\epsilon)$ corresponding to a
density $\rho(r)$ in a potential $\Psi(r)$ can be deduced from the first three
derivatives of $\rho(r)$ and $\Psi(r)$. This appendix works them out for
various practical cases that involve polytropes and NFW potentials. They all
are used in the main text.

###  A.1 Distribution function for a Schuster-Plummer stellar mass density in
a NFW potential

The Schuster-Plummer density (Eq. [30]) is defined as,

$D(r)=\rho(0)\left[1+\frac{r^{2}}{r_{0}^{2}}\right]^{-\frac{5}{2}},$ (A1)

so that,

$\frac{dD}{dr}=D_{1}(r)=-\frac{5\rho(0)}{r_{0}^{2}}\,r\,\left[1+\frac{r^{2}}{r_{0}^{2}}\right]^{-\frac{7}{2}},$
(A2)
$\frac{d^{2}D}{dr^{2}}=D_{2}(r)=-\frac{5\rho(0)}{r_{0}^{2}}\left[1+\frac{r^{2}}{r_{0}^{2}}\right]^{-\frac{7}{2}}+\frac{35\rho(0)}{r_{0}^{4}}\,r^{2}\,\left[1+\frac{r^{2}}{r_{0}^{2}}\right]^{-\frac{9}{2}},$
(A3)

and

$\frac{d^{3}D}{dr^{3}}=D_{3}(r)=\frac{35\rho(0)}{r_{0}^{4}}\,r\,\left[1+\frac{r^{2}}{r_{0}^{2}}\right]^{-\frac{9}{2}}+\frac{70\rho(0)}{r_{0}^{4}}\,r\,\left[1+\frac{r^{2}}{r_{0}^{2}}\right]^{-\frac{9}{2}}-\frac{315\rho(0)}{r_{0}^{6}}\,r^{3}\,\left[1+\frac{r^{2}}{r_{0}^{2}}\right]^{-\frac{11}{2}}.$
(A4)

On the other hand, the NFW density profile is defined as,

$\rho_{\rm NFW}(r)=\frac{\rho_{s}}{(r/r_{s})(1+r/r_{s})^{2}},$ (A5)

with $r_{s}$ and $\rho_{s}$ two constants. It creates a potential given by
(e.g., Binney & Tremaine, 2008),

$\Phi_{\rm NFW}(r)=-\frac{V_{c}}{r}\ln\left(1+\frac{r}{r_{s}}\right),$ (A6)

with $V_{c}=4\pi G\rho_{s}r_{s}^{3}$. Then the relative potential $\Psi(r)$,
denoted for the NFW profile as $V(r)$, turns out to be,

$V(r)=\Phi_{\rm NFW}(\infty)-\Phi_{\rm
NFW}(r)=\frac{V_{c}}{r}\ln\left(1+\frac{r}{r_{s}}\right),$ (A7)

with its derivatives given by,

$\frac{dV}{dr}=V_{1}(r)=\frac{V_{c}}{r^{2}}\left[\frac{r}{r+r_{s}}-\ln\left(1+\frac{r}{r_{s}}\right)\right],$
(A8)
$\frac{d^{2}V}{dr^{2}}=V_{2}(r)=-\frac{2V_{c}}{r^{3}}\left[\frac{r}{r+r_{s}}-\ln\left(1+\frac{r}{r_{s}}\right)\right]-\frac{V_{c}}{r(r+r_{s})^{2}},$
(A9)

and

$\frac{d^{3}V}{dr^{3}}=V_{3}(r)=\frac{6V_{c}}{r^{4}}\left[\frac{r}{r+r_{s}}-\ln\left(1+\frac{r}{r_{s}}\right)\right]+\frac{2V_{c}}{r^{2}(r+r_{s})^{2}}+\frac{V_{c}(3r+r_{s})}{r^{2}(r+r_{s})^{3}}.$
(A10)

Using Eqs. (5) and (11), the DF corresponding to a density given by Eq. (A1)
and a potential set by Eq. (A7) turns out to be,

$f(\epsilon)=\frac{1}{\pi^{2}\sqrt{2}}\int_{R}^{\infty}\,dr\,\sqrt{\epsilon-V(r)}\,\,\left[\frac{-D_{3}}{V_{1}^{2}}+\frac{3D_{2}V_{2}}{V_{1}^{3}}+\frac{D_{1}V_{3}}{V_{1}^{3}}-\frac{3D_{1}V_{2}^{2}}{V_{1}^{4}}\right],$
(A11)

with the limit $R$ implicitly defined as $\epsilon=V(R)$.

###  A.2 Distribution function for a Schuster-Plummer stellar mass density in
a Schuster-Plummer potential

The gravitational potential corresponding to the mass density in Eq. (A1) is
(e.g., Binney & Tremaine, 2008),

$\Phi_{\rm SP}(r)=-W_{c}\left(1+\frac{r^{2}}{r_{0}^{2}}\right)^{-1/2},$ (A12)

with

$W_{c}=\left(\frac{4\pi}{3}\right)Gr_{0}^{2},$ (A13)

so that the corresponding relative potential becomes,

$W(r)=\Phi_{\rm SP}(\infty)-\Phi_{\rm
SP}(r)=W_{c}\left(1+\frac{r^{2}}{r_{0}^{2}}\right)^{-1/2},$ (A14)

with its derivatives given by,

$\frac{dW}{dr}=W_{1}(r)=-\frac{W_{c}}{r_{0}^{2}}r\left(1+\frac{r^{2}}{r_{0}^{2}}\right)^{-3/2},$
(A15)
$\frac{d^{2}W}{dr^{2}}=W_{2}(r)=-\frac{W_{c}}{r_{0}^{2}}\left(1+\frac{r^{2}}{r_{0}^{2}}\right)^{-3/2}+\frac{3W_{c}}{r_{0}^{4}}r^{2}\left(1+\frac{r^{2}}{r_{0}^{2}}\right)^{-5/2},$
(A16)

and

$\frac{d^{3}W}{dr^{3}}=W_{3}(r)=\frac{9W_{c}}{r_{0}^{4}}r\left(1+\frac{r^{2}}{r_{0}^{2}}\right)^{-5/2}-\frac{15W_{c}}{r_{0}^{6}}r^{3}\left(1+\frac{r^{2}}{r_{0}^{2}}\right)^{-7/2}.$
(A17)

Using Eqs. (5) and (11), the DF corresponding to a density given by Eq. (A1)
and a potential set by Eq. (A14) turns out to be,

$f(\epsilon)=\frac{1}{\pi^{2}\sqrt{2}}\int_{R}^{\infty}\,dr\,\sqrt{\epsilon-W(r)}\,\,\left[\frac{-D_{3}}{W_{1}^{2}}+\frac{3D_{2}W_{2}}{W_{1}^{3}}+\frac{D_{1}W_{3}}{W_{1}^{3}}-\frac{3D_{1}W_{2}^{2}}{W_{1}^{4}}\right],$
(A18)

with the radius $R$ implicitly defined as $\epsilon=W(R)$. In the case of a
self-gravitating system, so that the Schuster-Plummer potential is the one
created by the Schuster-Plummer mass density, then Eq. (A18) can be integrated
analytically to yield,

$f(\epsilon)=\frac{\rho(0)}{W_{c}^{5}}\frac{120}{(2\pi)^{3/2}\Gamma(9/2)}\,\epsilon^{7/2},$
(A19)

an expression used to check our numerical evaluations of $f(\epsilon)$.

##  Appendix B The terms at $\Psi=0$ in the Eddington inversion method

The relative potential $\Psi$ generated by a spherically symmetric system of
finite total mass behaves as $\Psi\propto r^{-1}$ for $r\to\infty$ (e.g., Eq.
[35]). Consider objects where $\rho\propto r^{-b}$ for $r\to\infty$ (e.g., Eq.
[31]). Combining the asymptotic behaviors of $\Psi$ and $\rho$, one finds that
$d\rho/d\Psi\propto r^{1-b}$ and $d^{2}\rho/d\Psi^{2}\propto r^{2-b}$.
Consequently,
$\left(d\rho/d\Psi\right)_{\Psi=0}=\left(d^{2}\rho/d\Psi^{2}\right)_{\Psi=0}=0$
provided $b>2$. The above argument is not strictly valid if $\Psi$ stands for
the NFW potential, because it does not correspond to a mass distribution with
finite total mass, and it behaves as $\Psi\propto r^{-1}\,\ln r$ for
$r\to\infty$. However, if $b>2$, using Eq. (7) one can show that still
$\left(d\rho/d\Psi\right)_{\Psi=0}=\left(d^{2}\rho/d\Psi^{2}\right)_{\Psi=0}=0$
when $r\to\infty$, in spite of the logarithmic factor appearing in the NFW
potential.

##  Appendix C Distribution function of an spherically-symmetric system with
circular orbits

Any density $\rho(r)$ can be reproduced with a system of spherically-symmetric
circular orbits (Binney & Tremaine, 2008, Sect. 4.3.2). By definition, their
radial velocity is zero, $v_{r}=0$, and their tangential velocity is equal to
the circular velocity, $v_{t}=v_{c}(r)$, with

$v_{c}^{2}(r)=\frac{G\,M_{p}(<r)}{r}.$ (C1)

The symbol $M_{p}(<r)$ stands for the inner mass creating the potential
$\Phi$. A general DF with the required properties is,

$f_{c}(r,v_{t},v_{r})=F(r)\,\delta(v_{t}-v_{c})\,\delta(v_{r}),$ (C2)

where $\delta$ represents a Dirac-delta function and $F$ is a function to be
set by the density. Since, $\rho$ is recovered from the integral of $f_{c}$
over all velocities,

$\rho(r)=2\pi\int\\!\\!\\!\\!\int f_{c}\,v_{t}\,dv_{t}\,dv_{r},$ (C3)

then

$F(r)=\frac{\rho(r)}{2\pi\,v_{c}(r)},$ (C4)

which, together with Eq. (C1), uniquely defines $F$ for any combination of
$\rho$ and $\Phi$. Note that even if the DF in Eq. (C2) is not explicitly
written in terms of $\epsilon$ and $L$, it is straightforward to verify that
it is a stationary DF.

##  Appendix D The theorem by An & Evans in our context

Section 2.4 puts forward the DF,

$f(\epsilon,L)=L^{-2\beta}f_{\epsilon}(\epsilon),$ (D1)

which represents a leading order approximation for a wide class of DFs having
the anisotropy parameter $\beta$ constant (An & Evans, 2006; Binney &
Tremaine, 2008). In this case the DF depends not only on the energy $\epsilon$
but also on the modulus of the angular momentum $L$. Under this assumption,
the mass volume density can be written as (Binney & Tremaine 2008, Eq.
[4.66]),

$r^{2\beta}\rho(r)=\kappa_{\beta}\,\int_{0}^{\Psi}\frac{f_{\epsilon}(\epsilon)}{(\Psi-\epsilon)^{\beta-1/2}}\,d\epsilon,$
(D2)

where $\kappa_{\beta}$ is a positive numerical value independent of the radius
$r$. As we argued in Sect. 3, for the integral in the RHS of Eq. (D2) to be
zero, $f_{\epsilon}<0$ somewhere which through Eq. (D1) makes $f$ unphysical.
Assuming $\rho(r)\propto r^{-\alpha}$ when $r\rightarrow 0$, then the left-
hand-side of Eq. (D2) differs from zero if $2\beta-\alpha\leq 0$, which is the
theorem proved by An & Evans. Here we go a step further and provided
$\beta<1/2$ (Eq. [D2] diverges when taking derivatives and $\beta>1/2$), one
obtains (Binney & Tremaine 2008 Eq. [4.67]),

$\frac{d\big{[}r^{2\beta}\rho(r)\big{]}\big{/}dr}{d\Psi/dr}=\kappa_{\beta}\,\left(\frac{1}{2}-\beta\right)\int_{0}^{\Psi}\frac{f_{\epsilon}(\epsilon)}{(\Psi-\epsilon)^{\beta+1/2}}\,d\epsilon.$
(D3)

In the case of $\Psi$ given by a NFW potential, its radial derivative at $r=0$
differs from zero and is negative (Eq. [25]) therefore, to avoid the RHS of
Eq. (D3) to be less or equal to zero (and so to avoid an unphysical
$f_{\epsilon}<0$), $2\beta-\alpha\not=0$ and $2\beta-\alpha-1\leq 0$. The 2nd
condition is automatically met because $2\beta-\alpha$ is already $\leq 0$
according to An & Evans. Together with this inequality, the first condition
implies that for $f>0$ then

$\alpha>2\beta.$ (D4)

Our derivation assumes $\beta<1/2$, however, it is not difficult to show that
the inequality still holds in the limit case when $\beta=1/2$ and Eq. (D3) is
not valid.

There are several obvious consequences of the inequality in Eq. (D4): (1)
cores ($\alpha=0$) are inconsistent with isotropic velocities ($\beta=0$), (2)
cores are inconsistent with radially biased velocities (i.e., only $\beta<0$
is allowed), (3) radially biased orbits ($0<\beta<1/2$) require cuspy baryon
density profiles ($\alpha>0$), and (3) circular orbits do not pose any problem
since $\beta=-\infty$.

##  Appendix E Value of $(d\rho/dr)\big{/}(d\Psi/dr)$ when $r\rightarrow 0$
and $d\rho/dr\rightarrow 0$

Starting out from the definition of $\rho_{abc}$ in Eq. (31), one finds for
$r\rightarrow 0$,

$\frac{d\rho_{abc}}{dr}\simeq-\frac{\rho_{s}}{r_{s}}\,\frac{c+(b-c)\,x^{a}}{x^{1+c}},$
(E1)

with $x=r/r_{s}$. Similarly, Eq. (36) provides the 1st order approximation,

$\frac{d\Psi}{dr}\simeq-B\,r^{1-c_{p}},$ (E2)

where $c_{p}$ is the value of $c$ of the density profile assumed to generate
the potential $\Psi$ and $B$ is a positive constant. Putting together the two
previous equations, one finds,

$\frac{d\rho/dr}{d\Psi/dr}\simeq\frac{D}{r^{2+c-c_{p}}}\,\Big{[}c+\frac{b-c}{r_{s}^{a}}\,r^{a}\Big{]},$
(E3)

with $D>0$ for $c_{p}<3$. In the range of interest for galaxies, the
parameters $c$ and $c_{p}$ are from $\sim 0$ to $\sim 1$, with $c\leq c_{p}$,
whereas $a$ is between $\sim 1$ and $\sim 2$. If $c\not=0$, the first term in
the RHS of Eq. (E3) dominates the behavior of the ratio when $r\rightarrow 0$.
This term is identical to Eq. (26) with $\alpha=c$ and $\alpha_{p}=c_{p}$, and
its behavior is discussed in detail in Sect. 3. The case when $c=0$ is
particularly interesting since it represents a cored density profile. Only the
2nd term in the RHS of Eq. (E3) is not zero and it turns out to be,

$\frac{d\rho/dr}{d\Psi/dr}\simeq\frac{D\,b}{r_{s}^{a}}\frac{1}{r^{2-c_{p}-a}}.$
(E4)

Thus, the potential and the density profiles would be inconsistent when
$2-c_{p}-a<0$ since the ratio of derivatives goes to zero when $r\rightarrow
0$. It implies that when $a=2$ (e.g., Schuster-Plummer profile; Appendix A.2),
all $c_{p}>0$ are physically irrealizable (see Fig. 4 for $c=0$). It also
implies that when $c_{p}=0$, and so the potential has the same core as the
density, the pair density and potential are unphysical for $a>2$. This
somewhat surprising behavior has been checked numerically (Figs. 6 and 7).

## References

  * Abbott et al. (2021) Abbott, T. M. C., Adamów, M., Aguena, M., et al. 2021, ApJS, 255, 20, doi: 10.3847/1538-4365/ac00b3
  * An & Zhao (2013) An, J., & Zhao, H. 2013, MNRAS, 428, 2805, doi: 10.1093/mnras/sts175
  * An & Evans (2006) An, J. H., & Evans, N. W. 2006, ApJ, 642, 752, doi: 10.1086/501040
  * Battaglia & Nipoti (2022) Battaglia, G., & Nipoti, C. 2022, Nature Astronomy, 6, 659, doi: 10.1038/s41550-022-01638-7
  * Bechtol et al. (2022) Bechtol, K., Birrer, S., Cyr-Racine, F.-Y., et al. 2022, arXiv e-prints, arXiv:2203.07354.  https://arxiv.org/abs/2203.07354
  * Behroozi et al. (2013) Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57, doi: 10.1088/0004-637X/770/1/57
  * Benson (2010) Benson, A. J. 2010, Phys. Rep., 495, 33, doi: 10.1016/j.physrep.2010.06.001
  * Binney & Tremaine (2008) Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition
  * Blumenthal et al. (1984) Blumenthal, G. R., Faber, S. M., Primack, J. R., & Rees, M. J. 1984, Nature, 311, 517, doi: 10.1038/311517a0
  * Breddels & Helmi (2013) Breddels, M. A., & Helmi, A. 2013, A&A, 558, L3, doi: 10.1051/0004-6361/201322449
  * Brown et al. (2020) Brown, S. T., McCarthy, I. G., Diemer, B., et al. 2020, MNRAS, 495, 4994, doi: 10.1093/mnras/staa1491
  * Bullock (2002) Bullock, J. S. 2002, in The Shapes of Galaxies and their Dark Halos, ed. P. Natarajan, 109–113, doi: 10.1142/9789812778017_0018
  * Bullock & Boylan-Kolchin (2017) Bullock, J. S., & Boylan-Kolchin, M. 2017, ARA&A, 55, 343, doi: 10.1146/annurev-astro-091916-055313
  * Carlsten et al. (2021) Carlsten, S. G., Greene, J. E., Greco, J. P., Beaton, R. L., & Kado-Fong, E. 2021, ApJ, 922, 267, doi: 10.3847/1538-4357/ac2581
  * Cen (2014) Cen, R. 2014, ApJ, 790, L24, doi: 10.1088/2041-8205/790/2/L24
  * Chan et al. (2015) Chan, T. K., Kereš, D., Oñorbe, J., et al. 2015, MNRAS, 454, 2981, doi: 10.1093/mnras/stv2165
  * Ciotti (1996) Ciotti, L. 1996, ApJ, 471, 68, doi: 10.1086/177954
  * Ciotti & Pellegrini (1992) Ciotti, L., & Pellegrini, S. 1992, MNRAS, 255, 561, doi: 10.1093/mnras/255.4.561
  * Combes (2004) Combes, F. 2004, in Recycling Intergalactic and Interstellar Matter, ed. P.-A. Duc, J. Braine, & E. Brinks, Vol. 217, 440, doi: 10.48550/arXiv.astro-ph/0308293
  * Cuddeford (1991) Cuddeford, P. 1991, MNRAS, 253, 414, doi: 10.1093/mnras/253.3.414
  * Davis et al. (1985) Davis, M., Efstathiou, G., Frenk, C. S., & White, S. D. M. 1985, ApJ, 292, 371, doi: 10.1086/163168
  * Davis et al. (1992) —. 1992, Nature, 356, 489, doi: 10.1038/356489a0
  * Dejonghe (1987) Dejonghe, H. 1987, MNRAS, 224, 13, doi: 10.1093/mnras/224.1.13
  * Dekel et al. (2009) Dekel, A., Birnboim, Y., Engel, G., et al. 2009, Nature, 457, 451, doi: 10.1038/nature07648
  * Del Popolo & Le Delliou (2017) Del Popolo, A., & Le Delliou, M. 2017, Galaxies, 5, 17, doi: 10.3390/galaxies5010017
  * Di Cintio et al. (2014a) Di Cintio, A., Brook, C. B., Dutton, A. A., et al. 2014a, MNRAS, 441, 2986, doi: 10.1093/mnras/stu729
  * Di Cintio et al. (2014b) Di Cintio, A., Brook, C. B., Macciò, A. V., et al. 2014b, MNRAS, 437, 415, doi: 10.1093/mnras/stt1891
  * Dodelson & Widrow (1994) Dodelson, S., & Widrow, L. M. 1994, Phys. Rev. Lett., 72, 17, doi: 10.1103/PhysRevLett.72.17
  * Eddington (1916) Eddington, A. S. 1916, MNRAS, 76, 572, doi: 10.1093/mnras/76.7.572
  * El-Badry et al. (2017) El-Badry, K., Wetzel, A. R., Geha, M., et al. 2017, ApJ, 835, 193, doi: 10.3847/1538-4357/835/2/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, doi: 10.1093/mnras/stac3799
  * Governato et al. (2010) Governato, F., Brook, C., Mayer, L., et al. 2010, Nature, 463, 203, doi: 10.1038/nature08640
  * Grebel et al. (2003) Grebel, E. K., Gallagher, John S., I., & Harbeck, D. 2003, AJ, 125, 1926, doi: 10.1086/368363
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Hayashi et al. (2020) Hayashi, K., Chiba, M., & Ishiyama, T. 2020, ApJ, 904, 45, doi: 10.3847/1538-4357/abbe0a
  * Hernquist (1990) Hernquist, L. 1990, ApJ, 356, 359, doi: 10.1086/168845
  * Hu et al. (2000) Hu, W., Barkana, R., & Gruzinov, A. 2000, Phys. Rev. Lett., 85, 1158, doi: 10.1103/PhysRevLett.85.1158
  * Ivezić et al. (2019) Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111, doi: 10.3847/1538-4357/ab042c
  * Jackson et al. (2021) Jackson, R. A., Martin, G., Kaviraj, S., et al. 2021, MNRAS, 502, 4262, doi: 10.1093/mnras/stab077
  * Kormendy & Freeman (2016) Kormendy, J., & Freeman, K. C. 2016, ApJ, 817, 84, doi: 10.3847/0004-637X/817/2/84
  * Kowalczyk & Łokas (2022) Kowalczyk, K., & Łokas, E. L. 2022, A&A, 659, A119, doi: 10.1051/0004-6361/202142212
  * Lacroix et al. (2018) Lacroix, T., Stref, M., & Lavalle, J. 2018, J. Cosmology Astropart. Phys, 2018, 040, doi: 10.1088/1475-7516/2018/09/040
  * Leung et al. (2021) Leung, G. Y. C., Leaman, R., Battaglia, G., et al. 2021, MNRAS, 500, 410, doi: 10.1093/mnras/staa3107
  * Łokas (2009) Łokas, E. L. 2009, MNRAS, 394, L102, doi: 10.1111/j.1745-3933.2009.00620.x
  * Lotz et al. (2008) Lotz, J. M., Jonsson, P., Cox, T. J., & Primack, J. R. 2008, MNRAS, 391, 1137, doi: 10.1111/j.1365-2966.2008.14004.x
  * Lynden-Bell (1962) Lynden-Bell, D. 1962, MNRAS, 123, 447, doi: 10.1093/mnras/123.5.447
  * Massari et al. (2018) Massari, D., Breddels, M. A., Helmi, A., et al. 2018, Nature Astronomy, 2, 156, doi: 10.1038/s41550-017-0322-y
  * Massari et al. (2020) Massari, D., Helmi, A., Mucciarelli, A., et al. 2020, A&A, 633, A36, doi: 10.1051/0004-6361/201935613
  * Merritt et al. (2006) Merritt, D., Graham, A. W., Moore, B., Diemand, J., & Terzić, B. 2006, AJ, 132, 2685, doi: 10.1086/508988
  * Montes & Trujillo (2019) Montes, M., & Trujillo, I. 2019, MNRAS, 482, 2838, doi: 10.1093/mnras/sty2858
  * Montes & Trujillo (2022) —. 2022, ApJ, 940, L51, doi: 10.3847/2041-8213/ac98c5
  * Moskowitz & Walker (2020) Moskowitz, A. G., & Walker, M. G. 2020, ApJ, 892, 27, doi: 10.3847/1538-4357/ab7459
  * Navarro et al. (1997) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493, doi: 10.1086/304888
  * Oñorbe et al. (2015) Oñorbe, J., Boylan-Kolchin, M., Bullock, J. S., et al. 2015, MNRAS, 454, 2092, doi: 10.1093/mnras/stv2072
  * Orkney et al. (2023) Orkney, M. D. A., Taylor, E., Read, J. I., et al. 2023, arXiv e-prints, arXiv:2302.12818.  https://arxiv.org/abs/2302.12818
  * Orkney et al. (2021) Orkney, M. D. A., Read, J. I., Rey, M. P., et al. 2021, MNRAS, 504, 3509, doi: 10.1093/mnras/stab1066
  * Peñarrubia et al. (2012) Peñarrubia, J., Pontzen, A., Walker, M. G., & Koposov, S. E. 2012, ApJ, 759, L42, doi: 10.1088/2041-8205/759/2/L42
  * Peebles (2021) Peebles, P. J. E. 2021, arXiv e-prints, arXiv:2106.02672.  https://arxiv.org/abs/2106.02672
  * Plastino & Plastino (1993) Plastino, A. R., & Plastino, A. 1993, Physics Letters A, 174, 384, doi: 10.1016/0375-9601(93)90195-6
  * Putko et al. (2019) Putko, J., Sánchez Almeida, J., Muñoz-Tuñón, C., et al. 2019, ApJ, 883, 10, doi: 10.3847/1538-4357/ab365a
  * Read et al. (2016) Read, J. I., Agertz, O., & Collins, M. L. M. 2016, MNRAS, 459, 2573, doi: 10.1093/mnras/stw713
  * Read et al. (2019) Read, J. I., Walker, M. G., & Steger, P. 2019, MNRAS, 484, 1401, doi: 10.1093/mnras/sty3404
  * Roychowdhury et al. (2013) Roychowdhury, S., Chengalur, J. N., Karachentsev, I. D., & Kaisina, E. I. 2013, MNRAS, 436, L104, doi: 10.1093/mnrasl/slt123
  * Sánchez Almeida (2022) Sánchez Almeida, J. 2022, Universe, 8, 214, doi: 10.3390/universe8040214
  * Sánchez Almeida et al. (2014) Sánchez Almeida, J., Elmegreen, B. G., Muñoz-Tuñón, C., & Elmegreen, D. M. 2014, A&A Rev., 22, 71, doi: 10.1007/s00159-014-0071-1
  * Sánchez Almeida et al. (2017) Sánchez Almeida, J., Filho, M. E., Dalla Vecchia, C., & Skillman, E. D. 2017, ApJ, 835, 159, doi: 10.3847/1538-4357/835/2/159
  * Sánchez Almeida & Trujillo (2021) Sánchez Almeida, J., & Trujillo, I. 2021, MNRAS, 504, 2832, doi: 10.1093/mnras/stab1103
  * Sánchez Almeida et al. (2020) Sánchez Almeida, J., Trujillo, I., & Plastino, A. R. 2020, A&A, 642, L14, doi: 10.1051/0004-6361/202039190
  * Sánchez Almeida et al. (2021) —. 2021, ApJ, 921, 125, doi: 10.3847/1538-4357/ac1ba8
  * Sánchez-Janssen et al. (2016) Sánchez-Janssen, R., Ferrarese, L., MacArthur, L. A., et al. 2016, ApJ, 820, 69, doi: 10.3847/0004-637X/820/1/69
  * Smoot et al. (1992) Smoot, G. F., Bennett, C. L., Kogut, A., et al. 1992, ApJ, 396, L1, doi: 10.1086/186504
  * Spergel & Steinhardt (2000) Spergel, D. N., & Steinhardt, P. J. 2000, Phys. Rev. Lett., 84, 3760, doi: 10.1103/PhysRevLett.84.3760
  * Strigari et al. (2017) Strigari, L. E., Frenk, C. S., & White, S. D. M. 2017, ApJ, 838, 123, doi: 10.3847/1538-4357/aa5c8e
  * Trujillo et al. (2021) Trujillo, I., D’Onofrio, M., Zaritsky, D., et al. 2021, A&A, 654, A40, doi: 10.1051/0004-6361/202141603
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Weinberg et al. (2015) Weinberg, D. H., Bullock, J. S., Governato, F., Kuzio de Naray, R., & Peter, A. H. G. 2015, Proceedings of the National Academy of Science, 112, 12249, doi: 10.1073/pnas.1308716112
  * White & Rees (1978) White, S. D. M., & Rees, M. J. 1978, MNRAS, 183, 341, doi: 10.1093/mnras/183.3.341