11institutetext: Dipartimento di Fisica ”Giuseppe Occhialini”, Universitá di
Milano Bicocca, Piazza della Scienza 3 20126, Milano, Italy 22institutetext:
INFN-Sezione di Milano Via Celoria 15 20133, Milano, Italy
22email<EMAIL_ADDRESS>33institutetext: CNR-ISC, via Madonna del Piano
17 50022 Sesto Fiorentino, Italy 44institutetext: INAF-Osservatorio
Astronomico di Arcetri, Largo Enrico Fermi 5 50125 Firenze Italy
55institutetext: INFN-Sezione di Firenze, via Sansone 1 50022 Sesto
Fiorentino, Italy
55email<EMAIL_ADDRESS>
# Structure of the equivalent Newtonian systems in MOND N-body simulations

Density profiles and the core-cusp problem

Federico Re 1122 Pierfrancesco Di Cintio 334455

(Received ??; accepted ??)

###### Abstract

Aims. We investigate the core-cusp problem of the $\Lambda$ cold dark matter
($\Lambda$CDM) scenario in the context of Modified Newtonian Dynamics (MOND)
paradigm exploiting the concept of equivalent Newtonian system (ENS)

Methods. By means of particle-mesh $N-$body simulations in MOND we explore
processes of galaxy formation via cold dissipationless collapse or merging of
smaller substructures. From the end states of our simulations, we recover the
associated ENS and study the properties of their dark matter halos. We compare
the simulation results with simple analytical estimates with a family of
$\gamma-$models.

Results. We find that the dark matter density of ENSs of most spherical cold
collapses have a markedly cored structure, in particular for the lowest values
of the initial virial ratios. End states of some simulations with clumpy
initial conditions have more complex profiles and some of their ENSs exhibit a
moderate cusp, with logarithmic density slope always shallower than 1.

Conclusions. These results seem to point towards the fact that the absence in
most observed galaxies of a central DM cusp, at variance with what one would
expect from theoretical and numerical arguments in $\Lambda$CDM, would be
totally consistent in a MONDian description.

###### Key Words.:

Galaxies: kinematics and dynamics - Galaxies: formation - Gravitation -
Methods: numerical - Methods: analytical

##  1 Introduction

In the $\Lambda$ cold dark matter scenario (hereafter $\Lambda$CDM),
theoretical arguments and collisionless $N-$body simulations (Navarro et al.
1997) predict that galaxies are embedded in dark matter (DM) halos
characterized by a $\rho(r)\propto r^{-1}$ central cusp. Observational results
seem to suggest, from the analysis of the central velocity dispersion profiles
of dwarf galaxies, that the DM distribution has a cored density
distribution111Technically speaking, in multi component equilibrium self
gravitating systems there exist analytical constraints on the magnitude of a
component’s density given the logarithmic slope of the other, e.g. see
Dubinski & Carlberg 1991; Ciotti & Pellegrini 1992; Ciotti 1996, 1999.
Moreover for a broad range of spherical density profiles, the central density
slope constraints the value of the central anisotropy profile (An & Evans
2006) (see Moore 1994; Di Cintio et al. 2014).
Several solutions to this apparent contradiction -often referred to as “the
core-cusp problem”- such as self-interacting DM (e.g. see Lovell et al. 2012;
Nguyen et al. 2021; Eckert et al. 2022), DM annihilation (e.g. see Vasiliev
2007), baryon feedback (e.g. see Governato et al. 2010; Cole et al. 2011;
Pontzen & Governato 2012; Del Popolo & Pace 2016) or simply a
misinterpretation of the observational data (McGaugh et al. 2003), have been
proposed so far. However, notwithstanding the great amount of theoretical and
observational work, a clear answer is still far from being obtained. Moreover,
given the large interest in alternative theories of gravity, such as among the
others $f(R)$ gravities (Buchdahl 1970; Sotiriou & Faraoni 2010); Modified
Gravity (MoG, Moffat 2006; Moffat & Rahvar 2013); Retarded Gravity (Raju 2012;
Yahalom 2022); Emergent Gravity (Verlinde 2011, 2017); Refracted Gravity
(Cesare et al. 2020; Sanna et al. 2023); Fractional Gravity (Giusti 2020;
Benetti et al. 2023), proposed to avoid introducing the DM as a collisionless
fluid of exotic particles, it is natural to ask what becomes of the core-cusp
problem in those proposals.
In this work we investigate this matter in the modified Newtonian dynamics
(hereafter MOND, Milgrom 1983) paradigm. We recall that in the Bekenstein &
Milgrom (1984) Lagrangian formulation of MOND (sometimes referred to as AQUAL)
the classical Poisson equation for a density-potential pair $(\rho;\Phi)$

$\Delta\Phi=4\pi G\rho$ (1)

is substituted by the non-linear field equation

$\nabla\cdot\left[\mu\left(\frac{||\nabla\Phi||}{a_{0}}\right)\nabla\Phi\right]=4\pi
G\rho.$ (2)

In the equation above where $a_{0}\approx 10^{-8}$cm s-2 is a scale
acceleration and $\mu(x)$ is the MOND interpolating (monotonic) function known
only by its asymptotic limits

$\mu(x)\sim\cases{\displaystyle}1,\quad x\gg 1,\\\ \displaystyle x,\quad x\ll
1;{}$ (3)

Figure 1: Ratio of the stellar to dark density in the ENS (top) and DM and
stellar density profiles (bottom left and bottom right) in units of $3M/4\pi
r_{c}^{3}$ for $\kappa=100$ and $\gamma=0$, $0.5$, $1$, $1.5$, $2$ and $2.5$.

so that for $||\nabla\Phi||\gg a_{0}$ Eq. (2) one recovers the Newtonian
regime, while for $||\nabla\Phi||\ll a_{0}$ one obtains the so-called deep-
MOND (hereafter dMOND) regime and Eq. (2) simplifies to

$\nabla\cdot\left[||\nabla\Phi||\nabla\Phi\right]=4\pi G\rho a_{0}.$ (4)

Note that, the non-linear operator in Equation (4) is the special case of the
$p-$Laplace operator (see e.g. Stein 1970) for $p=3$, while Eq. (1) would
correspond to the $p=2$ case. In this respect, Equation (2) somewhat
”interpolates” between the two regimes via the $\mu$ function. Note also that,
in both cases, any given baryonic mass density $\rho$ can be taken out from
Equation (1) obtaining the relation

$\mu\left(\frac{||\mathbf{g}_{M}||}{a_{0}}\right)\mathbf{g}_{M}=\mathbf{g}_{N}+\mathbf{S}$
(5)

between the MOND and Newtonian force fields $\mathbf{g}_{M}$ and
$\mathbf{g}_{N}$, and where $\mathbf{S}\equiv\nabla\times\mathbf{h}(\rho)$ is
a density-dependent solenoidal field. It can be proved that the latter is
identically null for systems in spherical, cylindrical or planar symmetry,
while it is generally non-zero for arbitrary configurations of mass. On which
extent the stellar system at hand with mass $M$ is dominated by MOND effects,
is usually quantified by the dimensionless parameter

$\kappa\equiv\frac{GM}{r_{c}^{2}a_{0}},$ (6)

where $r_{c}$ is the scale of the baryon distribution. That is, for $\kappa\gg
1$ the system is mainly in Newtonian regime, vice versa for $\kappa\leq 1$
MOND effects become strong at all scales.
For any given stationary model in MOND one can always build the equivalent
Newtonian system (hereafter ENS), defined as a system with the same baryonic
mass density $\rho_{*}$ plus a DM halo with density $\rho_{DM}$ such that
their total potential $\Phi$ satisfying Eq. (1) is the same as the MOND
potential entering Eq. (2) for the sole density $\rho_{*}$. We note that, in
principle, the positivity of the DM density of the ENS is not always assured
(see Milgrom 1986), in particular for flattened systems (see Ciotti et al.
2006, 2012; Ko 2016). We recall that Milgrom (2010) introduced a Quasi-linear
formulation of MOND (hereafter QuMOND) where the modified field equation has
the same form as Eq. (2), with $\nu(||\mathbf{g}_{N}/a_{0}||)$ in lieu of
$\nu(||\mathbf{g}_{M}/a_{0}||)$. The QuMOND interpolating function $\nu(y)$
can be recovered from $\mu(x)$ appearing in Eq. (2) as

$\nu=\frac{1}{\mu}.$ (7)

It is easy to show that, from a given baryonic density distribution, one
obtains the MONDian potential $\Phi=\Phi_{N}+\Phi_{pDM}$ by first solving a
classical Poisson equation for the Newtonian potential $\Phi_{N}$, that trough
an algebraic passage involving $\nu$ becomes the source for the potential
$\Phi_{pDM}$ of the so-called ”phantom Dark Matter” via a second application
of the Poisson equation. Notably, in this alternative bi-potential MOND
formulation, the DM is de facto interpreted as the effect of the second
potential. As in AQUAL, in QuMOND one can retrieve a dMOND regime, that for
spherical systems easily reads $g_{M}=\sqrt{a_{0}/g_{N}}g_{N}$ (see Milgrom
2021).
If on one hand several workers have investigated the differences between
static equilibrium models in Newtonian and MOND gravities, or the
interpretation of observations in both theories, on the other much less is
known about the formation and evolution of stellar systems. For obvious
reasons, in observed systems one has access only to de-projected properties
for both stellar and dark components in the Newtonian framework. Numerical
experiments, though with their intrinsic limitations, yield information on the
full phase-space of the simulated models; in particular the 3D density
profiles. In this paper, we explore the structure of the ENS of MOND $N-$body
simulations of galaxy formation, in order to shed some light on the
possibility that the core-cusp problem is a MOND artifact in this paradigm of
gravity. We stress the fact that the MOND core-cusp problem discussed here, is
different from that introduced recently by Eriksen et al. (2021) that deals
with the modified gravity versus modified inertia hypothesis (see Milgrom
2022).
The rest of this paper is structured as follows. In Sect. 2 we revise the
definition of ENS and discuss their properties. In Sect. 3 we introduce the
numerical models and the analysis of the simulations. In Sect. 4 we discuss
the properties of the simulations’ end states. Finally, in Sect. 5 we
summarize and discusses the implications and the relations to previous work.

Figure 2: Same as in Fig. 1 but for $\kappa=1$

##  2 Equivalent Newtonian systems

As anticipated above, the ENS of a MOND model is the Newtonian system with the
same stellar (baryonic) mass distribution $\rho_{*}$ with an additional dark
component $\rho_{DM}$ such that the total potential (and thus the force field)
is the same of the parent MOND system (see Sanders & Begeman 1994; Angus et
al. 2006). For the case of an isolated spherical system one has

$\rho_{DM}=(4\pi G)^{-1}\nabla\cdot(\mathbf{g}_{M}-\mathbf{g}_{N}),$ (8)

since the solenoidal term $\mathbf{S}$ vanishes. We stress the fact that,
Equation (5) in QuMOND can be rewritten exactly as

$\mathbf{g}_{M}=\nu\left(\frac{||\mathbf{g}_{N}||}{a_{0}}\right)\mathbf{g}_{N}.$
(9)

If Equation (9) is applied to a spherically symmetric system one has
$\nu(y)=x/y$, and the total density of its ENS (baryonic plus phantom DM, see
e.g. Hodson et al. 2020; Oria et al. 2021) becomes

$\rho_{*}(r)+\rho_{DM}(r)=\frac{d(y\nu)}{dy}\rho_{*}(r)-\frac{yd\nu}{dy}\frac{2}{r^{3}}\int_{0}^{r}\rho_{*}(r)r^{2}dr.$
(10)

Let us consider the family of spherical $\gamma-$models (Dehnen 1993; Tremaine
et al. 1994), with density profile given by

$\rho_{*}(r)=\frac{3-\gamma}{4\pi}\frac{Mr_{c}}{r^{\gamma}(r+r_{c})^{4-\gamma}},$
(11)

where $M$ is the total baryonic mass, $0\leq\gamma<3$ is the logarithmic
density slope and $r_{c}$ the scale radius.
If the density profile (11) is substituted in Eq. (10) one obtains

$\rho_{*}+\rho_{DM}=\rho_{*}\left[\frac{r_{c}}{r+r_{c}}\frac{d(y\nu)}{dy}-\frac{2}{3-\gamma}\frac{yd\nu}{dy}\right],$
(12)

where

$y=\frac{||\mathbf{g}_{N}||}{a_{0}}=\kappa\left(\frac{r}{r_{c}}\right)^{1-\gamma}\left(1+\frac{r}{r_{c}}\right)^{\gamma-3},$
(13)

with $\kappa$ defined in Equation (6). We note that, for small radii $r$,
Equation (13) tends to zero if $\gamma<1$, while it diverges for $\gamma>1$.
In practice, at least for the $\gamma<1$ case, even in central regions the
model falls in the MOND regime. Strong MOND corrections in the centre are
therefore associated with a dominant DM component $\rho_{DM}$ in the ENS.

###  2.1 Massive galaxies

Let us consider a typical $10^{12}M_{\odot}$ massive elliptical galaxy with a
scale radius of 3 kpc, modelled with a $\gamma$-model. In this case
$\kappa\approx 10^{2}$. Due to discreteness effects of the underlying stellar
system Equation (11) can be considered reliable until the radius that contains
a fraction of roughly $10^{-3}$ of the total mass $M$ (in this case
$10^{9}M_{\odot}$, i.e. the typical mass of its central supermassive black
hole). The Lagrangian radius enclosing such mass fraction is

$r_{10^{-3}}=\frac{r_{c}}{10^{\frac{3}{3-\gamma}}-1}.$ (14)

The region in MOND regime has a far smaller radius, that for $\gamma\leq 1$ is
obtained by $y(r_{10^{-3}})\cong 10^{5}r_{10^{-3}}^{2}/r_{c}^{2}$, varying
between $\cong 2\times 10^{2}$ and $\cong 2\times 10^{3}$. And thus, even in
the framework of (Qu)MOND the phantom DM halo does not really dominate in the
central region for a cored stellar density profile. In Figure 1 we plot for
$\gamma=0$, 0.5, 1, 1.5, 2 and 2.5 the ratio of stellar to phantom DM and
their respective radial density profiles for $\kappa=10^{2}$. We note that,
remarkably, models with a strong cusp (i.e. $\gamma>1$) have phantom DM halos
in their ENS characterized by a decreasing density inside the scale radius.
Vice versa, cored models are associated with ENS having halos with a weak cusp
and several slope changes.

###  2.2 Diffuse galaxies

Let us now consider a diffuse galaxy such that $\kappa\sim 1$; so that its
central region can fall in the MOND regime, even for radii bigger than
$r_{10^{-3}}$. Typically, this occurs again if $\gamma<1$. We find that
$\lim_{r\rightarrow 0}y(r)=0$ in the central region, hence $\nu(y)\sim
y^{-1/2}$. If substituted in (12), this yields

$\frac{\rho_{*}+\rho_{DM}}{\rho_{*}}\sim\frac{5-\gamma}{6-2\gamma}\sqrt{\frac{a_{0}r_{c}^{2}}{GM}}\left(\frac{r_{c}}{r}\right)^{\frac{1-\gamma}{2}};$
(15)

that is, the phantom DM component dominates also at small radii. In
particular, the latter has a central profile given by

$\rho_{DM}\sim\frac{(5-\gamma)a_{0}}{8\pi
Gr_{c}}\sqrt{\frac{a_{0}r_{c}^{2}}{GM}}\left(\frac{r_{c}}{r}\right)^{\frac{1+\gamma}{2}}.$
(16)

The equation above is characterized by a weak cusp with a logarithmic density
slope $\alpha=-\frac{1+\gamma}{2}>-1$. For example, for $\gamma=0$, the DM
component in the ENS would have a cusp $\propto r^{-1/2}$. We note that this
trend is valid for any spherically symmetric stellar distribution with a
central core, and not only for the $\gamma=0$ Dehnen model. We note also it
always implies a central weak cusp with logarithmic density slope
$\alpha=-1/2$ for the phantom dark matter. This has the interesting
astrophysical implication that a galaxy with a cored stellar density profile
could be indeed interpreted in the DM scenario as having a cored halo, due to
the fact that week cusps can be often mistaken for cores.
For the cases with $\gamma>1$ and $\kappa=1$, for which the gravitational
field diverges in the centre, even though the stellar density is diffuse, the
ENS is DM-dominated only in the external region. This can be easily checked by
substituting the asymptotic behaviour $\nu(y)\sim
1+\frac{1}{y}-\frac{1}{y^{2}}+o(y^{-2})$, and finding from (12) that

$\frac{\rho_{DM}}{\rho_{*}}\sim\frac{2}{3-\gamma}\frac{a_{0}r_{c}^{2}}{GM}\left(\frac{r}{r_{c}}\right)^{\gamma-1}-\frac{r}{r_{c}}+O(r^{\gamma}),$
(17)

implying a vanishingly small central DM density. This is summarized in Figure
2 where we plot the same quantities as in the previous Fig. but $\kappa=1$. As
expected, in the upper plot showing $\rho_{*}/\rho_{DM}$, the $\gamma=0$ and
0.5 (red and orange lines) are everywhere below 1, i.e. the system is
dominated by the phantom DM distribution at all radii. Vice versa, for the
$\gamma\geq 1$ cases, the phantom DM of the ENS dominates only in the external
regions. We recall that Sánchez Almeida (2022) showed that galaxies with
central regions in MOND regime imply ENS characterized by a decreasing baryon
density and a cored DM.

##  3 Numerical code and models

###  3.1 Numerical code and initial conditions

The $N-$body simulations discussed here have been performed with a modified
version of the publicly available nmody particle-mesh MOND code (Nipoti et al.
2007a, see also Londrillo & Nipoti 2011 for additional technical details). The
latter uses a non-linear Poisson solver to compute $\Phi$ from Eq. (2) on a
$N_{r}\times N_{\vartheta}\times N_{\varphi}$ spherical grid in polar
coordinates, using an iterative relaxation procedure starting from a guess
solution (here given by Eq. 5 neglecting $S$), as for the linear Poisson
methods (see Londrillo & Messina 1990; Londrillo et al. 1991). As a rule, in
the simulations discussed here we used a $128\times 32\times 64$ grid. In this
work we adopt the following form for the interpolation function

$\mu(x)=\frac{x}{\sqrt{1+x^{2}}}.$ (18)

Alternative choices can also be implement always leading to qualitatively
similar end states.
The equation of motion are integrated using a standard 4th order leapfrog
scheme (see e.g. Dehnen & Read 2011) with an adaptive timestep $\Delta t$
conditioned by the stability threshold $\Delta t=C/\sqrt{{\rm
max}|\nabla\cdot\mathbf{g}|}$, where the Courant-Friedrichs-Lewy condition $C$
was taken in the range $0.01\leq C\leq 0.1$.
We performed two sets of numerical simulations with initial conditions defined
as follows.
In the first, the particles positions were sampled from Equation (11) while in
the second, following Hansen et al. (2006), we first distribute according to a
Poissonian distribution inside a larger $\gamma$ model the centres of $N_{C}$
clumps222Clumpy initial conditions were also explored in the context of
Newtonian simulations by Nipoti (2015) and Ludlow & Angulo (2017) when
investigating the relation of the initial density fluctuation power spectrum
with the Sérsic index $m$ (see below) conjectured by Cen (2014). also
described by Eq. (11) with different choices of $r_{c}$, and $\gamma$ and
later populate them with particles.
In both cases, the initial particle velocities are extracted from a position-
independent isotropic Maxwell-Boltzmann distribution and normalized to obtain
the wanted value of the initial virial ratio $2K/|W|$, where $K$ is the total
kinetic energy and $W$ the virial function, defined for a (finite mass)
continuum system of density $\rho$ and potential $\Phi$ as

$W=-\int\rho(\mathbf{r})\langle\mathbf{r},\nabla\Phi\rangle d^{3}\mathbf{r}.$
(19)

We recall that in isolated dMOND systems of finite mass
$W=-2\sqrt{GM^{3}a_{0}}/3$ is constant (see Nipoti et al. 2007a). Curiously,
even in systems of particles interacting with additive $1/r$ forces with
logarithmic potential the virial function is constant (see Di Cintio et al.
2013, 2017).
The simulations of this work span a range of $N$ between $10^{4}$ and
$10^{6}$. All simulations were extended up to $t=300t_{\rm Dyn}$, where
$t_{\rm Dyn}\equiv\sqrt{2r_{h}^{3}/GM_{\rm tot}}$ and $r_{h}$ is the radius
containing half of the total mass of the system $M_{\rm tot}$, so that virial
oscillations and phase-mixing are likely to be complete.
Following Ciotti et al. (2007) in some cases we enforce the spherical symmetry
during the collapse by propagating particles only using the radial part of the
evaluated force field, so that the system behaves effectively as a spherical
shell model introduced in Newtonian gravity by Hénon (1964) and used in MOND
among the others by Sanders (2008); Malekjani et al. (2009) and by Di Cintio &
Ciotti (2011) for systems interacting with $1/r^{\alpha}$ forces.

Figure 3: Final ($t=300t_{\rm Dyn}$) density profiles from MOND simulations (top panels) and the DM halo of the ENS (bottom panels) for cored $\gamma_{0}=0$ (left), moderately cuspy $\gamma=1$ (centre) and clumpy (right) initial conditions. The increasing initial values of the virial ratio in the models with spherical initial conditions with $\gamma=0$ and 1 is mapped with increasingly lighter tones of blue and green respectively. All clumpy initial conditions start with $2K_{0}/|W|_{0}=0.1$ Table 1: Summary of the simulation properties: After the name of each simulation (Col. 1), we report the number of particles (Col. 2), the gravity law (MOND or dMOND, Col. 3), the initial density profile (Col. 4), the initial virial ratio (Col. 5) the axial ratios (Cols. 6 and 7) , the Sérsic (Col. 8), the final anisotropy index (Col. 9) and the virial velocity dispersion (Col. 10). Name | Gravity | Initial profile | $N$ | $2K_{0}/|W_{0}|$ | $c/a$ | $b/a$ | $m$ | $\xi$ | $\sigma_{\rm vir}$ | $\alpha$
---|---|---|---|---|---|---|---|---|---|---
gamma0v0 | MOND | $\gamma=0$ | $3\times 10^{4}$ | $10^{-4}$ | $0.36$ | $0.56$ | $2.45$ | $3.04$ | $0.83$ | $0.25$
gamma05v0 | MOND | $\gamma=0.5$ | $3\times 10^{4}$ | $10^{-4}$ | $0.42$ | $0.62$ | $2.33$ | $2.84$ | $0.84$ | $-0.15$
gamma1v0 | MOND | $\gamma=1$ | $3\times 10^{4}$ | $10^{-4}$ | $0.54$ | $0.98$ | $2.07$ | $2.73$ | $0.89$ | $-0.11$
gamma15v0 | MOND | $\gamma=1.5$ | $3\times 10^{4}$ | $10^{-4}$ | $0.53$ | $0.94$ | $0.89$ | $2.91$ | $1.03$ | $-0.23$
gamma2v0 | MOND | $\gamma=2$ | $3\times 10^{4}$ | $10^{-4}$ | $0.57$ | $0.92$ | $0.72$ | $2.26$ | $1.05$ | $-0.35$
gamma0ve1m3 | MOND | $\gamma=0$ | $3\times 10^{4}$ | $10^{-3}$ | $0.36$ | $0.56$ | $3.34$ | $3.13$ | $0.83$ | $0.21$
gamma0ve3m3 | MOND | $\gamma=0$ | $3\times 10^{4}$ | $3\times 10^{-3}$ | $0.33$ | $0.50$ | $4.28$ | $3.15$ | $0.85$ | $0.50$
gamma0ve1m2 | MOND | $\gamma=0$ | $3\times 10^{4}$ | $10^{-2}$ | $0.34$ | $0.57$ | $2.53$ | $3.10$ | $0.84$ | $0.40$
gamma0ve3m2 | MOND | $\gamma=0$ | $3\times 10^{4}$ | $3\times 10^{-2}$ | $0.32$ | $0.52$ | $3.06$ | $3.33$ | $0.86$ | $0.52$
gamma0ve1m1 | MOND | $\gamma=0$ | $3\times 10^{4}$ | $0.1$ | $0.32$ | $0.46$ | $2.23$ | $2.98$ | $0.84$ | $0.70$
gamma0ve2m1 | MOND | $\gamma=0$ | $3\times 10^{4}$ | $0.2$ | $0.33$ | $0.38$ | $3.07$ | $3.03$ | $0.83$ | $0.75$
gamma0ve3m1 | MOND | $\gamma=0$ | $3\times 10^{4}$ | $0.3$ | $0.36$ | $0.37$ | $3.29$ | $3.02$ | $0.83$ | $0.85$
gamma0ve4m1 | MOND | $\gamma=0$ | $3\times 10^{4}$ | $0.4$ | $0.37$ | $0.38$ | $2.91$ | $3.01$ | $0.82$ | $0.89$
gamma0ve5m1 | MOND | $\gamma=0$ | $3\times 10^{4}$ | $0.5$ | $0.40$ | $0.42$ | $2.50$ | $2.92$ | $0.83$ | $1.10$
gamma1ve1m3 | MOND | $\gamma=1$ | $3\times 10^{4}$ | $10^{-3}$ | $0.49$ | $0.87$ | $1.90$ | $2.86$ | $0.90$ | $-0.21$
gamma1ve3m3 | MOND | $\gamma=1$ | $3\times 10^{4}$ | $3\times 10^{-3}$ | $0.50$ | $0.91$ | $3.75$ | $2.96$ | $0.91$ | $-0.14$
gamma1ve1m2 | MOND | $\gamma=1$ | $3\times 10^{4}$ | $10^{-2}$ | $0.48$ | $0.67$ | $2.48$ | $3.18$ | $0.91$ | $-0.12$
gamma1ve3m2 | MOND | $\gamma=1$ | $3\times 10^{4}$ | $3\times 10^{-2}$ | $0.47$ | $0.82$ | $4.10$ | $3.20$ | $0.93$ | $-0.10$
gamma1ve1m1 | MOND | $\gamma=1$ | $3\times 10^{4}$ | $0.1$ | $0.42$ | $0.62$ | $2.66$ | $3.32$ | $0.90$ | $0.00$
gamma1ve2m1 | MOND | $\gamma=1$ | $3\times 10^{4}$ | $0.2$ | $0.45$ | $0.46$ | $3.10$ | $2.99$ | $0.87$ | $0.15$
gamma1ve3m1 | MOND | $\gamma=1$ | $3\times 10^{4}$ | $0.3$ | $0.51$ | $0.52$ | $2.25$ | $3.33$ | $0.88$ | $0.21$
gamma1ve4m1 | MOND | $\gamma=1$ | $3\times 10^{4}$ | $0.4$ | $0.54$ | $0.54$ | $2.87$ | $3.46$ | $0.88$ | $0.25$
gamma1ve5m1 | MOND | $\gamma=1$ | $3\times 10^{4}$ | $0.5$ | $0.95$ | $0.96$ | $3.04$ | $3.71$ | $0.87$ | $0.45$
gamma1v1em1 | MOND | $\gamma=1$ | $5\times 10^{4}$ | $0.1$ | $0.51$ | $0.95$ | $1.24$ | $2.43$ | $0.88$ | $-0.15$
gamma1v0b | MOND | $\gamma=1$ | $2.1\times 10^{5}$ | $0$ | $0.37$ | $0.69$ | $2.58$ | $2.90$ | $0.89$ | $-0.15$
gamma1v0dmd | dMOND | $\gamma=1$ | $2.1\times 10^{5}$ | $0$ | $0.24$ | $0.41$ | $2.45$ | $3.49$ | $0.82$ | $0.50$
clumpy1 | MOND | clumpy | $8.7\times 10^{4}$ | $0.1$ | $0.46$ | $0.85$ | $1.91$ | $2.72$ | $0.87$ | $0.70$
clumpy2 | MOND | clumpy | $8.7\times 10^{4}$ | $0.1$ | $0.46$ | $0.67$ | $2.86$ | $2.48$ | $0.84$ | $0.50$
clumpy3 | MOND | clumpy | $8.7\times 10^{4}$ | $0.1$ | $0.60$ | $0.97$ | $3.32$ | $1.44$ | $1.11$ | $-0.40$
clumpy4 | MOND | clumpy | $8.7\times 10^{4}$ | $0.1$ | $0.57$ | $0.67$ | $1.98$ | $2.07$ | $0.88$ | $0.05$
clumpy5 | MOND | clumpy | $8.7\times 10^{4}$ | $0.1$ | $0.66$ | $0.94$ | $3.47$ | $1.62$ | $1.18$ | $-0.75$
clumpy6 | MOND | clumpy | $8.7\times 10^{4}$ | $0.1$ | $0.41$ | $0.70$ | $1.75$ | $1.98$ | $0.85$ | $-0.80$
clumpy7 | MOND | clumpy | $8.7\times 10^{4}$ | $0.1$ | $0.30$ | $0.54$ | $1.67$ | $2.22$ | $0.83$ | $1.10$
clumpy8 | MOND | clumpy | $8.7\times 10^{4}$ | $0.1$ | $0.33$ | $0.59$ | $3.39$ | $2.02$ | $0.84$ | $0.95$
clumpy9 | MOND | clumpy | $8.7\times 10^{4}$ | $0.1$ | $0.80$ | $0.98$ | $1.64$ | $1.41$ | $1.08$ | $-0.40$
clumpy10 | MOND | clumpy | $8.7\times 10^{4}$ | $0.1$ | $0.56$ | $0.71$ | $0.98$ | $1.89$ | $0.85$ | $-0.99$
clumpy1dmd | dMOND | clumpy | $8.7\times 10^{4}$ | $0.1$ | $0.13$ | $0.28$ | $3.46$ | $7.45$ | $0.77$ | $1.01$
clumpy2dmd | dMOND | clumpy | $8.7\times 10^{4}$ | $0.1$ | $0.36$ | $0.61$ | $3.16$ | $2.69$ | $0.83$ | $1.00$
clumpy3dmd | dMOND | clumpy | $8.7\times 10^{4}$ | $0.1$ | $0.35$ | $0.38$ | $1.40$ | $1.94$ | $0.82$ | $-1.99$
clumpy4dmd | dMOND | clumpy | $8.7\times 10^{4}$ | $0.1$ | $0.30$ | $0.37$ | $1.74$ | $2.34$ | $0.81$ | $0.70$
clumpy5dmd | dMOND | clumpy | $8.7\times 10^{4}$ | $0.1$ | $0.55$ | $0.56$ | $0.82$ | $2.08$ | $0.82$ | $-0.50$
clumpy6dmd | dMOND | clumpy | $8.7\times 10^{4}$ | $0.1$ | $0.33$ | $0.62$ | $1.55$ | $2.00$ | $0.83$ | $0.90$
clumpy7dmd | dMOND | clumpy | $8.7\times 10^{4}$ | $0.1$ | $0.26$ | $0.57$ | $1.44$ | $2.24$ | $0.82$ | $0.91$
clumpy8dmd | dMOND | clumpy | $8.7\times 10^{4}$ | $0.1$ | $0.29$ | $0.59$ | $1.89$ | $2.00$ | $0.81$ | $0.99$
clumpy9dmd | dMOND | clumpy | $8.7\times 10^{4}$ | $0.1$ | $0.35$ | $0.39$ | $1.56$ | $2.10$ | $0.81$ | $0.80$
clumpy10dmd | dMOND | clumpy | $8.7\times 10^{4}$ | $0.1$ | $0.42$ | $0.60$ | $1.76$ | $1.70$ | $0.75$ | $0.61$
gamma0v5em51D | MOND (1D) | $\gamma=0$ | $3\times 10^{4}$ | $10^{-4}$ | $0.96$ | $0.98$ | $2.26$ | $20.5$ | $0.71$ | $-0.01$
gamma05v5em51D | MOND (1D) | $\gamma=0.5$ | $3\times 10^{4}$ | $10^{-4}$ | $0.98$ | $0.99$ | $2.27$ | $33.1$ | $0.98$ | $-0.45$
gamma1v5em51D | MOND (1D) | $\gamma=1$ | $3\times 10^{4}$ | $10^{-4}$ | $0.97$ | $0.97$ | $1.62$ | $88.0$ | $1.63$ | $-0.40$
gamma15v5em51D | MOND (1D) | $\gamma=1.5$ | $3\times 10^{4}$ | $10^{-4}$ | $0.97$ | $0.99$ | $0.71$ | $44.8$ | $1.48$ | $-0.21$
gamma2v5em51D | MOND (1D) | $\gamma=2$ | $3\times 10^{4}$ | $10^{-4}$ | $0.96$ | $0.97$ | $0.51$ | $451$ | $4.72$ | $-0.35$
gamma1v3em31D | MOND (1D) | $\gamma=1$ | $3\times 10^{4}$ | $3\times 10^{-3}$ | $0.97$ | $1.00$ | $2.10$ | $4.47$ | $0.66$ | $0.01$
gamma1v1em21D | MOND (1D) | $\gamma=1$ | $3\times 10^{4}$ | $10^{-2}$ | $0.97$ | $0.99$ | $3.78$ | $15.8$ | $0.60$ | $-0.10$
gamma1v3em21D | MOND (1D) | $\gamma=1$ | $3\times 10^{4}$ | $3\times 10^{-2}$ | $0.96$ | $0.98$ | $3.11$ | $5.19$ | $0.63$ | $-0.51$
gamma1v1em11D | MOND (1D) | $\gamma=1$ | $3\times 10^{4}$ | $0.1$ | $0.97$ | $0.99$ | $2.42$ | $4.02$ | $0.64$ | $-0.62$
gamma1v2em11D | MOND (1D) | $\gamma=1$ | $3\times 10^{4}$ | $0.2$ | $0.98$ | $0.99$ | $3.61$ | $3.19$ | $0.62$ | $-0.45$
gamma1v3em11D | MOND (1D) | $\gamma=1$ | $3\times 10^{4}$ | $0.3$ | $0.99$ | $0.99$ | $4.09$ | $2.57$ | $0.63$ | $-0.35$

###  3.2 Analysis of the end products

For all simulations presented here we first extract the intrinsic properties
of the end products from their phase-space positions. We first evaluate the
triaxiality of the final particle distribution (see e.g. Nipoti et al. 2006a;
Di Cintio et al. 2013 and references therein) by defining the tensor

$I_{ij}\equiv m\sum_{k=1}^{N}r_{i}^{(k)}r_{j}^{(k)}$ (20)

for the particles with positions $\mathbf{r}_{i}$ within the Lagrangian radius
$r_{70}$ containing the 70% of the stellar mass of the system and evaluating
with a standard iterative procedure its three eigenvalues $I_{1}\geq I_{2}\geq
I_{3}$. By applying a rotation $\mathcal{R}$ to all particles of the system so
that the three associated eigenvectors are now oriented along the coordinate
axes we then get that the three semiaxes $a\geq b\geq c$ from $I_{1}=Aa^{2}$,
$I_{2}=Ab^{2}$ and $I_{3}=Ac^{2}$, where $A$ is a numerical constant depending
on the density profile. Finally, we define the axial ratios
$b/a=\sqrt{I_{2}/I_{1}}$ and $c/a=\sqrt{I_{3}/I_{1}}$, and the ellipticities
in the principal planes $\epsilon_{1}=1-\sqrt{I_{2}/I_{1}}$ and
$\epsilon_{2}=1-\sqrt{I_{3}/I_{1}}$.
Following Nipoti et al. (2007a) and Di Cintio et al. (2013) we compare the
surface density profiles of the end products with the Sersic (1968) law

$\Sigma(R)=\Sigma_{e}e^{-b\left[\left(\frac{R}{R_{e}}\right)^{1/m}-1\right]},$
(21)

where $\Sigma_{e}$ is the projected mass density at effective radius $R_{e}$,
the radius of the circle containing half of the projected mass, and the
dimensionless parameters $b,m$ are related by $b\simeq 2m-1/3+4/405m$ as found
by Ciotti & Bertin (1999).
Once the projected density in the 3 principal planes is circularized over
elliptical shells, we determine the corresponding pair $(R_{e},\Sigma_{e})$ by
particle counts (i.e. we are assuming a constant mass to light ratio for each
particle), and fit Eq. (21) for the three projections. We find that, in
general, all 3 sets of ($\Sigma_{e}$, $R_{e}$, $m$) are rather similar
(differing only for less than the 5%), we therefore chose randomly only one.
In addition, for all simulations we also evaluate the so-called anisotropy
index (see Binney & Tremaine 2008) defined by

$\xi=\frac{2K_{r}}{K_{t}},$ (22)

where $K_{r}$ and $K_{t}=K_{\theta}+K_{\phi}$ are the radial and tangential
components of the kinetic energy tensor, respectively and read

$K_{r}=2\pi\int\rho(r)\sigma^{2}_{r}(r)r^{2}{\rm d}r,\quad
K_{t}=2\pi\int\rho(r)\sigma^{2}_{t}(r)r^{2}{\rm d}r.$ (23)

In the expressions above, $\sigma^{2}_{r}$ and $\sigma^{2}_{t}$ are the radial
and tangential phase-space averaged square velocity components and are
obtained for the end products of the simulations by particle counts over
radial shells.
For each simulation we recover the (spherical) DM density of the ENS from
Equation (8) where the Newtonian force field $g_{N}$ has been evaluated and
averaged on the radial coordinate. In practice, we are assuming a
”sphericized” system.

Figure 4: Inner density slope of the ENS halo (top panel), best fit Sérsic
index (middle panel) and minor ellipticity $\epsilon=1-c/a$ as function of the
initial virial ratio for initial conditions with Dehnen profiles with
$\gamma=0$ (circles) and 1 (triangles)

Finally, for the density distribution $\rho_{DM}$ so obtained we evaluate the
logarithmic density slope $\alpha$. We find that the profiles of $\rho_{DM}$
are generally well fitted by the empirical law

$\rho(r)=\frac{\rho_{\alpha}r_{\alpha}^{2}}{r^{\alpha}(r^{2}+r_{\alpha}^{2})^{\frac{2-\alpha}{2}}},$
(24)

where $r_{\alpha}$ is a scale radius and $\rho_{\alpha}$ is the associated
scale density. Equation (24) above recovers the $1/r^{2}$ trend of the density
of the ENS as predicted by the logarithmic behaviour of the far field MOND
potential. The properties of the simulations and their initial conditions are
summarized in Tab. 1 below.

##  4 $N-$body simulations

###  4.1 Spherical collapses

One of the main motivations of the present work is to establish whether the
end products of MOND dissipationless collapses could, in principle, reproduce
the structural properties of elliptical galaxies together with their inferred
dark halos. Single component Newtonian collapses with spherical initial
conditions, are known to produce flatter end states for increasing values of
their initial virial ratio (see Nipoti et al. 2006a, b; Di Cintio et al. 2013
and references therein) at fixed initial density profile.
We find that, this (partially) holds true for MOND spherical collapses, as
shown in Fig. 3 (top left and top mid panels), where we plot the baryon
density distribution at $300t_{\rm Dyn}$ for $\gamma=0$ and $1$ and increasing
values of the virial ratio with increasingly lighter tones of blue and green
in the range $10^{-3}\leq 2K_{0}/|W_{0}|\leq 0.5$. Using Equation (8) for the
angle-averaged final density profile on a spherical grid, we evaluated the
density distribution of the DM component of the parent Newtonian model (see
bottom panels, same figure). We find that, in qualitative agreement with the
structural properties of the ENS (see Figs. 1 and 2 in Sect. 2), cuspy end
systems can be associated with cored or weakly cuspy phantom halos. In
general, the end products of spherical collapses have always inner regions
that are baryon dominated when building their ENS, even if the initial
conditions are such that $\kappa=1$ (in particular for the $\gamma=0$ cases).
Consistently with Nipoti et al. (2007a), we observe that, independently on the
specific value of the initial virial ratio, initial conditions characterized
by a moderate density cusp (i.e. $0.5\leq\gamma\leq 2$) tend to yield end
products that are in general oblate (i.e. $0.5\lesssim c/a\lesssim b/a$), as
for Newtonian single component collapses. We typically observe major
ellipticities up to $\sim 0.63$ (corresponding to the gamma1v0b case, see Tab.
1). Remarkably, MOND collapses with cored initial conditions (i.e. $\gamma=0$)
evolve into rather prolate end states for $2K_{0}/|W_{0}|\gtrsim 0.1$, and
markedly triaxial end states for lower values of the initial virial ratio. For
both cored and moderately cuspy initial conditions, the inner slope $\alpha$
of the DM halo of the ENS, obtained by fitting with Eq. (24) increases for
increasing values of the baryon initial virial ratio in the MOND simulation,
as shown in Fig. 4 (top panel). The best fit Sérsic index $m$, measuring the
concentration of the projected stellar density profile is always in the range
$2\leq m\leq 4.5$ for both choices of the initial density profile (mid panel,
same figure), while the major ellipticity $\epsilon=1-c/a$ is typically larger
when the initial condition has a lower virial ratio, being smaller for larger
values of the initial $\gamma$ at fixed $2K_{0}/|W_{0}|$ (bottom panel, same
figure). Remarkably no system is found being more flattened than an E7 galaxy.
However, as also found by Nipoti et al. (2007a), dMOND collapses may produce
even flatter end states as in the case of the gamma1v0dmd run, for which
$c/a\sim 0.24$ so that $\epsilon=0.76$. For fixed initial virial ratio, the
end states attain larger values of the central virial velocity dispersion
$\sigma_{\rm vir}$ for increasing values of the initial density slope, while
the anisotropy index $\xi$ decreases (cfr. 1). At fixed initial density
profile, the final values of $\sigma_{\rm vir}$ have little variation with
$2K_{0}/|W_{0}|$, while $\xi$ is usually lower for the relaxed states of
hotter initial conditions.
In order to clarify whether the properties of the halo in the ENS are or not
an artifact of the angle averaging procedure, we also performed a set of
simulation in enforced spherical symmetry by propagating particles only using
the radial component of the force field. By doing so, the system remains
spherically symmetric (as no radial orbit instability is possible) and $S=0$
de

Figure 5: Inner density slope of the ENS halo (top panel), best fit Sérsic
index (middle panel) and minor ellipticity $\epsilon=1-c/a$ as function of the
logarithmic density slope $\gamma$ of the initial condition for full 3D
(filled symbols) and 1D simulations (empty symbols).

facto hold true at all times, so that one could apply Eq. (8) exactly. In
Figure 5 we show the same quantities as in Fig. 4 as function of the initial
values of $\gamma$ for systems starting with a virial ratio of $10^{-4}$ with
(empty symbols) and without (filled symbols) enforced spherical symmetry. The
trend as well as the values of $\alpha$ 3D and effective 1D simulations are
comparable, and the same could be noted also for the Sérsic index $m$ that
attains considerably lower values (associated with a more concentrated density
profile) for larger values of the initial logarithmic density slope. In all
cases (cfr. 1), as expected, 1D collapses relax to final stated with rather
large values of the orbital anisotropy $\xi$.
Figure 6 shows the final angle averaged density profiles for $\gamma_{0}=0$, 1
and 1.5 in 3D and 1D simulations (solid lines) as well as the density profiles
of the ENS halos (dashed lines). Notably, if on one hand the large $r$
behaviour of the baryon density profiles $\rho_{*}$ (where the systems are
mostly dominated by radial orbits) does not change significantly, on the
other, the inner slope of $\rho_{*}$ is always higher for the end products of
the 1D simulations and typically settles around $2.5$. With the sole exception
of the cored initial conditions ($\gamma_{0}=0$), the DM halo of the ENS of
the end products is denser (in units of the baryon component density
$\rho_{*,50}$ evaluated at the half mass radius $r_{50}$) for the 1D
simulations, being in both cases considerably shallower than the parent baryon
density.

Figure 6: Final baryon density profiles (coloured solid lines) and ENS halos
(coloured dashed lines) for $\gamma=0,$ 1 and 1.5. The black lines refer to
the 1D the cases with the same initial conditions. Figure 7: Matrix plot of
the Sérsic index, slope of the DM profile in the ENS, major ellipticity and
anisotropy index for simulations with clumpy (diamonds) and spherical
(circles) initial conditions. Empty symbol mark the dMOND runs. Figure 8: Mass
ratio against ellipticity relation for the final states of spherical (red
circles) and clumpy (green diamonds) initial conditions. The purple dashed
line marks the Deur (2014) relation with its uncertainty (blue shaded area),
while the orange dotted line marks the linear fit for the models with
spherical initial conditions.

###  4.2 Clumpy collapses

The numerical studies in MOND carried out so far, have typically explored
spherical initial conditions (see Nipoti et al. 2007a; Ciotti et al. 2007;
Sanders 2008; Malekjani et al. 2009; Nipoti et al. 2011), disks (Brada &
Milgrom 1999; Tiret & Combes 2007, 2008a; Nipoti et al. 2007c; Ghafourian &
Roshan 2017; Wittenburg et al. 2020) or galaxy merging (Nipoti et al. 2007b;
Tiret & Combes 2008b) and references therein. Here, in addition to the usual
spherical collapses we also explored clumpy initial conditions. When starting
with such initial states, MOND simulations tend (as expected) to yield
markedly triaxial end states with broader ranges of both $c/a$ and $b/a$. In
general, for fixed values of the initial virial ratio, the systems tend to
relax at later times with respect to their initially spherical counterparts
(the oscillations of $2K/W$ damp out at about $50t_{\rm Dyn}$ in spherical
collapses, see Nipoti et al. 2007a, while in clumpy systems this happens on
average at around $140t_{\rm Dyn}$) for analogous choices of the virial ratio.
We report here only the runs corresponding to $2K_{0}/|W_{0}|=0.1$, see Tab.
1.
The final three dimensional (angle averaged) density profiles (see top right
panel in Fig. 3) are strikingly more complex than those obtained from
spherical initial conditions and bare individually more slope changes. The
projected 2D density profiles are fitted by the Sérsic law with roughly the
same (percentage) asymptotic standard error of about, on average, 3% as for
the spherical collapses, while the scatter in the $m$ Sérsic parameter is
slightly smaller for MOND clumpy systems (see top panel in Fig. 7). For
comparison, we also run the same clumpy initial conditions in dMOND finding a
larger scatter in $m$.
As a general trend, the DM halo of the circularized ENS of clumpy collapses
are significantly more cored333Notably, in Newtonian simulations of clumps in
fall in DM halo Cole et al. (2011) found that the central DM cusp is
considerably weakened by the collapsing clumpy satellites. than what is
typically obtained in spherical collapses. In several cases the inner density
slopes are negative, down to $\sim-0.99$, corresponding to a DM density
profile that decreases in the central regions (middle panels in Fig. 7).
Interestingly, no initially clumpy system is found to evolve into a state
flatter than an E7 galaxy (thin dashed line in bottom panels of Fig. 7) in
MOND simulations. However, some dMOND collapses result in considerably flatter
end states (and often prolate) with major ellipticity reaching 0.87 for the
clumpy1dmd.
We observe that, final states with larger values of the anisotropy index $\xi$
(i.e. more and more dominated by low-angular momentum orbits) are always
associated to larger ellipticities $\epsilon$ and Sérsic indexes. A similar,
though somewhat weaker, correlation is also found between $\alpha$ and
$\epsilon$, that could be read in the DM scenario as steeper inner DM profiles
producing flatter stellar distributions.

###  4.3 The MOND mass-to-light ratio - ellipticity relation

Deur (2014, 2020) and more recently Winters et al. (2023), using a broad
sample of elliptical galaxies from independent surveys, and different methods
to evaluate the mass to light ratio $M/L$ (i.e. Jeans anisotropic modelling,
gravitational lensing, X-ray spectra and the dynamics of satellite star
clusters) and the ellipticity $\epsilon$, recovered the linear relation

$M/L=(14.1\pm 5.4)\epsilon,$ (25)

where the $M/L$ is normalized such that $M/L(\epsilon_{\rm app}=0.3)\equiv
8M_{\odot}/L_{\odot}\equiv 4M/M_{*}(\epsilon_{\rm app}=0.3)$, and the
intrinsic ellipticity $\epsilon$ is extrapolated from its observed 2D
projected value $\epsilon_{\rm app}$ assuming that all systems are oblate with
a Gaussian distribution of projection angles $\theta$ so that

$\epsilon_{\rm app}=1-\sqrt{(1-\epsilon)^{2}\sin^{2}\theta+\cos^{2}\theta}.$
(26)

The Equation above in the context of $\Lambda$CDM implies that a larger
contribution of the DM mass $M_{\rm DM}$ to the total mass $M$ corresponds to
a larger departure from the spherical symmetry (here quantified by larger
major ellipticity) for the stellar component. Winters et al. (2023) argue
that, if true, such a correlation would be contrasting the standard
$\Lambda$CDM scenario of galaxy formation, where more massive (and rather
spherical) DM halos are embed less flattened stellar systems. We note that
some peculiar elliptical galaxies (though excluded by the original sample of
Winters et al. 2023) such as the ultrafaint dwarfs (Simon 2019) appear to go
against the trend given by Eq. (25), having usually $\epsilon\lesssim 0.1$
with $M/L$ in some cases up to $10^{3}$.
Using the simulations discussed in the previous sections, we have investigated
the relation (25) in the context of MOND, evaluating the effective DM mass
$M_{\rm DM}$ in the ENSs of both clumpy and spherical collapses. To do so,
after recovering the $\rho_{DM}$ from the angle averaged ENS, we integrate it
radially up to the radius containing all simulation particles.
In Figure 8 we show the total to stellar mass (here we have assumed units such
that $M_{*}/L=1$) ratio $M/M_{*}$ versus major ellipticity $\epsilon$ for
collapses with both spherical and clumpy initial states, here indicated by
circles and diamonds respectively, as well as the observational relation given
in Eq. (25). We found that the end products of initially clumpy systems fall
in (almost) all cases within Winters et al. (2023)’s relation and its error
range (indicated in figure by the shaded area), while for the spherical
collapses fall on a rather steeper relation. We performed a linear fit (marked
in figure by the orange dotted line) obtaining

$M/M_{*}=(23.24\pm 0.59)\epsilon.$ (27)

We stress the fact, that none of the simulations discussed above produces
final states that could be interpreted as ultrafaint dwarfs (except, possibly,
some dMOND collapses), that in the standard cosmological scenario are supposed
to be DM dominated at all radii (i.e. even in the central region where our
simulations, when interpreted in the context of DM have baryon dominated
cores).

##  5 Discussion and conclusions

In this work we have investigated the structure of the dark matter density
profiles of the (angular averaged) equivalent Newtonian systems of the end
states of MOND dissipationless collapse simulations. We studied a broader
range of initial conditions than those discussed by Nipoti et al. (2007a),
including non spherical ones.
The main results of this work can be summarized as follows: Simple analytical
estimates in spherical symmetry suggest that the presence of a core or even
centrally decreasing DM distribution in ENS of MOND models with cuspy stellar
profiles. Vice versa, cored stellar profiles are associated with ENS DM
central density profiles $\rho_{DM}\propto 1/r^{\alpha}$ with $\alpha\lesssim
1$. Our MOND $N-$body simulations and the angle averaged ENS of their end
states nicely confirm this. This established, we can conclude that he flat-
cored halos invoked by some observational studies, can be reasonably
considered in agreement with our numerical finding, as the dynamical effect of
a weak cusp, independently on the specific value of the central logarithmic
density slope of the baryons, can be easily mistaken for that of a cored dark
mass distribution in the DM paradigm.
In general, we observe that as for the simulations in Newtonian gravity, in
MOND the stronger is the collapse (i.e. lower initial virial ratio and/or
larger initial density slope), the steeper is the final density profile, and
thus the dark halo of the ENS has a markedly cored, or sometimes even
depleted, inner density. Obviously, the end product of simplified MOND
$N-$body simulations with enforced spherical symmetry have ENS with markedly
flat cores, for a broad spectrum of initial values of density slope and virial
ratio, with baryon density always dominated by a rather strong cusp at inner
radii. Moreover, we also find that, if interpreted in the context of DM, the
relaxed end states with smaller values of the ellipticity (i.e. less
flattened) should have cuspier DM halos. In general, independently on the
specific form of the initial density profile, colder initial conditions are
always associated to flatter end states.
As a by-product of this simulation study on ENSs, we have also recovered a
numerical confirmation of the claimed Deur (2014) observational linear
correlation between $M/L$ (or $M/M_{*}$) and $\epsilon$, though with seemingly
different slope, when evaluating the dark matter content of ENSs in units of
the baryon mass (the latter being a pre-defined simulation parameter).
Our findings lead us to speculate that in the context of MOND the core cusp
problem could be a “MOND artifact” in the same sense as rings and DM shells
discussed by Milgrom & Sanders (2008). Moreover, we stress the fact that in
the DM halos reconstructed from observational data using the line-of-sight
velocity dispersion of a given tracer stellar population, the effect of the
velocity anisotropy profiles $\beta(r)=1-\sigma_{t}^{2}(r)/2\sigma_{r}^{2}(r)$
(and the intrinsic departure from the spherical symmetry) is neglected, as
noted by Evans et al. (2009) for the case of dwarf spheroids. In fact, since
the central stellar $\beta$ profile imposes a constraint on the slope of the
DM component in the form of the inequality $\beta\leq\alpha/2$ (see An & Evans
2006; Ciotti & Morganti 2009, 2010) the entity of the central density cusp or
core inferred for observed galaxies is likely to bare a rather large
uncertainty. In the context of (single component) MOND models, the relation
between anisotropy and central density cusps has not been explored in detail,
neither analytically nor in simulations. Simple numerical experiments (Di
Cintio et al. 2013) with inverse power-law radial forces seem to suggest that
the density slope anisotropy inequality is a rather general property of the
relaxed states of collapses with long-range interactions.
A natural follow up of this work will be a systematic study of the interplay
of the $\beta$ profiles in MOND systems and the DM density profiles of the
parent ENSs.

###### Acknowledgements.

We would like to express gratitude to Carlo Nipoti for the assistance with the
simulation in nmody and Michal Bílek for the discussions at an early stage of
this work. One of us (PFDC) wishes to acknowledge funding by “Fondazione Cassa
di Risparmio di Firenze” under the project HIPERCRHEL for the use of high
performance computing resources at the university of Firenze.

## References

  * An & Evans (2006) An, J. H. & Evans, N. W. 2006, ApJ, 642, 752
  * Angus et al. (2006) Angus, G. W., Famaey, B., & Zhao, H. S. 2006, MNRAS, 371, 138
  * Bekenstein & Milgrom (1984) Bekenstein, J. & Milgrom, M. 1984, ApJ, 286, 7
  * Benetti et al. (2023) Benetti, F., Lapi, A., Gandolfi, G., Salucci, P., & Danese, L. 2023, ApJ, 949, 65
  * Binney & Tremaine (2008) Binney, J. & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton University Press)
  * Brada & Milgrom (1999) Brada, R. & Milgrom, M. 1999, ApJ, 519, 590
  * Buchdahl (1970) Buchdahl, H. A. 1970, MNRAS, 150, 1
  * Cen (2014) Cen, R. 2014, ApJ, 790, L24
  * Cesare et al. (2020) Cesare, V., Diaferio, A., Matsakos, T., & Angus, G. 2020, A&A, 637, A70
  * Ciotti (1996) Ciotti, L. 1996, ApJ, 471, 68
  * Ciotti (1999) Ciotti, L. 1999, ApJ, 520, 574
  * Ciotti & Bertin (1999) Ciotti, L. & Bertin, G. 1999, A&A, 352, 447
  * Ciotti et al. (2006) Ciotti, L., Londrillo, P., & Nipoti, C. 2006, ApJ, 640, 741
  * Ciotti & Morganti (2009) Ciotti, L. & Morganti, L. 2009, MNRAS, 393, 179
  * Ciotti & Morganti (2010) Ciotti, L. & Morganti, L. 2010, MNRAS, 408, 1070
  * Ciotti et al. (2007) Ciotti, L., Nipoti, C., & Londrillo, P. 2007, in Collective Phenomena in Macroscopic Systems, 177–186
  * Ciotti & Pellegrini (1992) Ciotti, L. & Pellegrini, S. 1992, MNRAS, 255, 561
  * Ciotti et al. (2012) Ciotti, L., Zhao, H., & de Zeeuw, P. T. 2012, MNRAS, 422, 2058
  * Cole et al. (2011) Cole, D. R., Dehnen, W., & Wilkinson, M. I. 2011, Monthly Notices of the Royal Astronomical Society, 416, 1118
  * Dehnen (1993) Dehnen, W. 1993, MNRAS, 265, 250
  * Dehnen & Read (2011) Dehnen, W. & Read, J. I. 2011, European Physical Journal Plus, 126, 55
  * Del Popolo & Pace (2016) Del Popolo, A. & Pace, F. 2016, Astrophys. Space Sci., 361, 162, [Erratum: Astrophys.Space Sci. 361, 225 (2016)]
  * Deur (2014) Deur, A. 2014, MNRAS, 438, 1535
  * Deur (2020) Deur, A. 2020, arXiv e-prints, arXiv:2010.06692
  * Di Cintio et al. (2014) Di Cintio, A., Brook, C. B., Macciò, A. V., et al. 2014, MNRAS, 437, 415
  * Di Cintio & Ciotti (2011) Di Cintio, P. & Ciotti, L. 2011, International Journal of Bifurcation and Chaos, 21, 2279
  * Di Cintio et al. (2013) Di Cintio, P., Ciotti, L., & Nipoti, C. 2013, MNRAS, 431, 3177
  * Di Cintio et al. (2017) Di Cintio, P., Ciotti, L., & Nipoti, C. 2017, MNRAS, 468, 2222
  * Dubinski & Carlberg (1991) Dubinski, J. & Carlberg, R. G. 1991, ApJ, 378, 496
  * Eckert et al. (2022) Eckert, D., Ettori, S., Robertson, A., et al. 2022, A&A, 666, A41
  * Eriksen et al. (2021) Eriksen, M. H., Frandsen, M. T., & From, M. H. 2021, A&A, 656, A123
  * Evans et al. (2009) Evans, N. W., An, J., & Walker, M. G. 2009, MNRAS, 393, L50
  * Ghafourian & Roshan (2017) Ghafourian, N. & Roshan, M. 2017, MNRAS, 468, 4450
  * Giusti (2020) Giusti, A. 2020, Phys. Rev. D, 101, 124029
  * Governato et al. (2010) Governato, F., Brook, C., Mayer, L., et al. 2010, Nature, 463, 203
  * Hansen et al. (2006) Hansen, S., Moore, B., Zemp, M., & Stadel, J. 2006, JCAP, 0601, 014
  * Hénon (1964) Hénon, M. 1964, Annales d’Astrophysique, 27, 83
  * Hodson et al. (2020) Hodson, A. O., Diaferio, A., & Ostorero, L. 2020, A&A, 640, A26
  * Ko (2016) Ko, C.-M. 2016, ApJ, 821, 111
  * Londrillo & Messina (1990) Londrillo, P. & Messina, A. 1990, MNRAS, 242, 595
  * Londrillo et al. (1991) Londrillo, P., Messina, A., & Stiavelli, M. 1991, MNRAS, 250, 54
  * Londrillo & Nipoti (2011) Londrillo, P. & Nipoti, C. 2011, N-MODY: A Code for Collisionless N-body Simulations in Modified Newtonian Dynamics, Astrophysics Source Code Library, record ascl:1102.001
  * Lovell et al. (2012) Lovell, M. R., Eke, V., Frenk, C. S., et al. 2012, MNRAS, 420, 2318
  * Ludlow & Angulo (2017) Ludlow, A. D. & Angulo, R. E. 2017, MNRAS, 465, L84
  * Malekjani et al. (2009) Malekjani, M., Rahvar, S., & Haghi, H. 2009, ApJ, 694, 1220
  * McGaugh et al. (2003) McGaugh, S. S., Barker, M. K., & de Blok, W. J. G. 2003, ApJ, 584, 566
  * Milgrom (1983) Milgrom, M. 1983, ApJ, 270, 365
  * Milgrom (1986) Milgrom, M. 1986, ApJ, 306, 9
  * Milgrom (2010) Milgrom, M. 2010, MNRAS, 403, 886
  * Milgrom (2021) Milgrom, M. 2021, Phys. Rev. D, 103, 044043
  * Milgrom (2022) Milgrom, M. 2022, Phys. Rev. D, 106, 064060
  * Milgrom & Sanders (2008) Milgrom, M. & Sanders, R. H. 2008, ApJ, 678, 131
  * Moffat (2006) Moffat, J. W. 2006, J. Cosmology Astropart. Phys., 2006, 004
  * Moffat & Rahvar (2013) Moffat, J. W. & Rahvar, S. 2013, MNRAS, 436, 1439
  * Moore (1994) Moore, B. 1994, Nature, 370, 629
  * Navarro et al. (1997) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, Astrophys. J., 490, 493
  * Nguyen et al. (2021) Nguyen, Q. L., Mathews, G. J., Phillips, L. A., et al. 2021, Modern Physics Letters A, 36, 2130001
  * Nipoti (2015) Nipoti, C. 2015, ApJ, 805, L16
  * Nipoti et al. (2011) Nipoti, C., Ciotti, L., & Londrillo, P. 2011, MNRAS, 414, 3298
  * Nipoti et al. (2006a) Nipoti, C., Londrillo, P., & Ciotti, L. 2006a, MNRAS, 370, 681
  * Nipoti et al. (2006b) Nipoti, C., Londrillo, P., & Ciotti, L. 2006b, N-body simulations of dissipationless galaxy formation, Science and Supercomputing at CINECA, 2005, p.122-126
  * Nipoti et al. (2007a) Nipoti, C., Londrillo, P., & Ciotti, L. 2007a, ApJ, 660, 256
  * Nipoti et al. (2007b) Nipoti, C., Londrillo, P., & Ciotti, L. 2007b, MNRAS, 381, L104
  * Nipoti et al. (2007c) Nipoti, C., Londrillo, P., Zhao, H., & Ciotti, L. 2007c, MNRAS, 379, 597
  * Oria et al. (2021) Oria, P. A., Famaey, B., Thomas, G. F., et al. 2021, ApJ, 923, 68
  * Pontzen & Governato (2012) Pontzen, A. & Governato, F. 2012, MNRAS, 421, 3464
  * Raju (2012) Raju, C. K. 2012, in American Institute of Physics Conference Series, Vol. 1483, The Sixth International School on Field Theory and Gravitation-2012, ed. J. Alves Rodrigues, Waldyr, R. Kerner, G. O. Pires, & C. Pinheiro, 260–276
  * Sánchez Almeida (2022) Sánchez Almeida, J. 2022, ApJ, 940, 46
  * Sanders (2008) Sanders, R. H. 2008, MNRAS, 386, 1588
  * Sanders & Begeman (1994) Sanders, R. H. & Begeman, K. G. 1994, MNRAS, 266, 360
  * Sanna et al. (2023) Sanna, A. P., Matsakos, T., & Diaferio, A. 2023, A&A, 674, A209
  * Sersic (1968) Sersic, J. L. 1968, Atlas de Galaxias Australes
  * Simon (2019) Simon, J. D. 2019, ARA&A, 57, 375
  * Sotiriou & Faraoni (2010) Sotiriou, T. P. & Faraoni, V. 2010, Reviews of Modern Physics, 82, 451
  * Stein (1970) Stein, E. M. 1970, Singular Integrals and Differentiability Properties of Functions (PMS-30) (Princeton University Press)
  * Tiret & Combes (2007) Tiret, O. & Combes, F. 2007, A&A, 464, 517
  * Tiret & Combes (2008a) Tiret, O. & Combes, F. 2008a, A&A, 483, 719
  * Tiret & Combes (2008b) Tiret, O. & Combes, F. 2008b, in Astronomical Society of the Pacific Conference Series, Vol. 396, Formation and Evolution of Galaxy Disks, ed. J. G. Funes & E. M. Corsini, 259
  * Tremaine et al. (1994) Tremaine, S., Richstone, D. O., Byun, Y.-I., et al. 1994, AJ, 107, 634
  * Vasiliev (2007) Vasiliev, E. 2007, Phys. Rev. D, 76, 103532
  * Verlinde (2011) Verlinde, E. 2011, Journal of High Energy Physics, 2011, 29
  * Verlinde (2017) Verlinde, E. P. 2017, SciPost Physics, 2, 016
  * Winters et al. (2023) Winters, D. M., Deur, A., & Zheng, X. 2023, MNRAS, 518, 2845
  * Wittenburg et al. (2020) Wittenburg, N., Kroupa, P., & Famaey, B. 2020, ApJ, 890, 173
  * Yahalom (2022) Yahalom, A. 2022, International Journal of Modern Physics D, 31, 2242018