11institutetext: National Council of Research - Institute of Complex Systems
(CNR-ISC), Via Madonna del piano 10, I-50019 Sesto Fiorentino, Italy
22institutetext: National Institute of Nuclear Physics (INFN) - Florence unit,
via G. Sansone 1, I-50019 Sesto Fiorentino, Italy 33institutetext: Department
of Physics and Astronomy, University of Florence, Piazzale E. Fermi 2, I-50125
Firenze, Italy 44institutetext: National Institute of Astrophysics - Arcetri
Astrophysical Observatory (INAF-OAA), Piazzale E. Fermi 5, I-50125 Firenze,
Italy
44email<EMAIL_ADDRESS>
# Dissipationless collapse and the dynamical mass-ellipticity relation of
elliptical galaxies in Newtonian gravity and MOND

Pierfrancesco Di Cintio 11223344

(Received October xx, xxxx; accepted March xx, xxxx)

###### Abstract

Context. Deur (2014) and Winters et al. (2023) proposed an empirical relation
between the dark to total mass ratio and ellipticity in elliptical galaxies
from their observed total dynamical mass-to-light ratio data $M/L=(14.1\pm
5.4)\epsilon$. In other words, the larger is the content of dark matter in the
galaxy, the more the stellar component would be flattened. Such observational
claim, if true, appears to be in stark contrast with the common intuition of
the formation of galaxies inside dark halos with reasonably spherical
symmetry.

Aims. Comparing the processes of dissipationless galaxy formation in different
theories of gravity, and emergence of the galaxy scaling relations therein is
an important frame where, in principle one could discriminate them.

Methods. By means of collisionless $N$-body simulations in modified Newtonian
dynamics (MOND) and Newtonian gravity with and without active dark matter
halos, with both spherical and clumpy initial structure, I study the trends of
intrinsic and projected ellipticities, Sérsic index and anisotropy with the
total dynamical to stellar mass ratio.

Results. It is shown that, the end products of both cold spherical collapses
and mergers of smaller clumps depart more and more from the spherical symmetry
for increasing values of the total dynamical mass to stellar mass, at least in
a range of halo masses. The equivalent Newtonian systems of the end products
of MOND collapses show a similar behaviour. The $M/L$ relation obtained from
the numerical experiments in both gravities is however rather different from
that reported by Deur and coauthors.

###### Key Words.:

galaxies: formation – galaxies: evolution – gravitation – methods: numerical

##  1 Introduction

In the $\Lambda$CDM scenario, stellar systems such as galaxies and clusters
are thought to be embedded in Dark Matter (hereafter DM) halos, accounting for
the missing mass fraction evidenced by velocity dispersion or rotational
velocity measures (e.g. see Bertin 2014). DM halos are usually assumed to be
spherically symmetric structures in collisionless equilibrium with radial
densities well fitted by the Navarro et al. (1997) (hereafter NFW) profile,
originally obtained from cosmological dissipationless collapse simulations
(see Dubinski & Carlberg 1991).
In the context of spheroids such as elliptical galaxies and the bulges of disk
galaxies the interplay between the DM density distribution $\rho_{\rm DM}$ and
that of the of the stellar component $\rho_{*}$ has attracted a lot of
interest, in particular with respect to the well known cusp-core problem for
which the observed rotation curves of dwarf galaxies imply a flat-cored halo
(i.e. vanishing logarithmic density slope at inner radii) in contrast with the
cuspy (i.e. diverging DM density at $r\rightarrow 0$) distributions predicted
by DM $N-$body simulations such as for example the NFW profile (e.g. see Oh et
al. 2015). In addition, it can also be proved analytically that flat cored
models can not admit consistent phase-space distributions when embedded in
cuspy halos, except in a small range of values of the central anisotropy
profile (see Ciotti 1999, see also Ciotti & Pellegrini 1992).
Several explanations to the core-cusp problem have been proposed so far, from
the effect of stellar evolution feedback (Di Cintio et al. 2014; Brook et al.
2014; Genina et al. 2018; Katz et al. 2018; Burger & Zavala 2021) to extra
channels for dissipation, and hence the possibility of flattening their
central density distribution, for the DM component due to its intrinsic nature
(such as for example self-interacting DM, Baldi & Salucci 2012; Kahlhoefer et
al. 2019; Sánchez Almeida & Trujillo 2021; Burger et al. 2022) or even a
misinterpretation of the observational data on dwarf galaxies (McGaugh et al.
2003; Valenzuela et al. 2007). Moreover, in the context of modified theories
of gravity, as alternatives to the DM hypothesis, such as for example the
modified Newtonian dynamics (hereafter MOND, Milgrom 1983) explain the core-
cusp problem as an artifact of the gravity model (see e.g. Eriksen et al.
2021; Sánchez Almeida 2022; Re & Di Cintio 2023).
If on one hand much work has been devoted to the interplay between the slopes
of the dark and visible matter density profiles and their implications for the
anisotropy profile, on the other much less is known on the relationship of the
DM profile and the intrinsic shape of the stellar component. 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$ such as for example Virial theorem (Bacon et al. 1985);
anisotropic Jeans modelling (Cappellari 2008; Pechetti et al. 2017);
gravitational lensing (Bolton et al. 2008; Auger et al. 2010; Shu et al.
2017); gas disk dynamics (Bertola et al. 1991, 1993; Pizzella et al. 1997);
X-ray emission spectra (Jin et al. 2020) and the dynamics of satellite star
clusters or companion galaxies (Alabi et al. 2017; Harris et al. 2019; Chen &
Hwang 2020), and the ellipticity $\epsilon$ observed that the two are related
by the linear relation

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

where the mass to light ratio is normalized such that $M/L(\epsilon_{\rm
2D}=0.3)\equiv 8M_{\odot}/L_{\odot}\equiv 4M/M_{*}(\epsilon_{\rm 2D}=0.3)$,
and the intrinsic ellipticity $\epsilon$ is inferred from its apparent
2D-projected value $\epsilon_{\rm 2D}$ in the assumption of oblateness and a
Gaussian distribution of projection angles $\theta$ from

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

Equation (1) implies that (at least for the galaxies examined) a larger
contribution of the DM mass $M_{\rm DM}$ to the total mass $M$ corresponds to
a larger departure from the spherical symmetry (i.e. larger ellipticity). The
lower DM fraction in rounder galaxies could be, in principle, perceived as
being in stark contrast with the standard scenario of galaxy formation111Note
that, on the contrary, globular clusters are essentially spherical while not
containing a significant fraction of DM, requiring that the baryons aggregates
and forms structure in dark matter halos that had collapsed at earlier times.
Moreover, some peculiar elliptical galaxies (though excluded by the original
sample of Winters et al. 2023) such as the dwarf (Mateo 1998) or ultrafaint
dwarf galaxies (Simon 2019) appear to go against the trend given by Eq. (1) as
they are characterized by a rather spherical shape ($\epsilon\leq 0.1$) while
being DM dominated with $M/L$ up to $10^{2}$ or more. Of course, one of the
reasons behind the mass ratio - ellipticity relation could be related to the
departure from spherical symmetry of the halos themselves and its relation to
their total DM content. In fact, several studies (both observational and
numerical, see e.g. Jing & Suto 2002; Ragone-Figueroa et al. 2010; Vera-Ciro
et al. 2011; Bett 2012; Schneider et al. 2012; Rojas-Niño_2015; Evslin 2016;
Gonzalez et al. 2022) point out that more massive halos have significantly
higher departure from the spherical symmetry. More recently, Re & Di Cintio
(2023) performed $N-$body simulations in MOND showing that the equivalent
Newtonian systems (ENS, i.e. stellar systems with a dark halo such that the
total Newtonian gravitational potential is the same as the parent MOND model)
of the end products of cold gravitational collapses are consistent with Eq.
(1) if the initial conditions are sufficiently clumpy.
In general, the origin of the triaxiality of (single component) gravitational
systems is usually ascribed to the process of radial-orbit instability222Some
authors, (see Joyce et al. 2009; Worrakitpoonpon 2015; Sylos Labini et al.
2015; Benhaiem et al. 2016) also suggest the angular momentum loss via
particle escape during the violent phases of gravitational collapse as a
viable channel for the loss of spherical symmetry. (hereafter ROI, see
Polyachenko & Shukhman 1981; Palmer & Papaloizou 1987; Bertin et al. 1994;
Maréchal & Perez 2010), where an initially minor departure from spherical
symmetry accretes particles on low angular momentum orbits, as shown in
numerical experiments; Merritt & Aguilar (1985); Bellovary et al. (2008);
Barnes et al. (2009); Gajda et al. (2015); Di Cintio et al. (2017); Di Cintio
& Casetti (2020); Weinberg (2023). The ROI is stronger in model with a larger
degree of radial anisotropy, quantified (see e.g. Binney & Tremaine 2008, see
also Fridman & Poliachenko 1984) by the parameter

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

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

$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,$ (4)

where $\sigma^{2}_{r}$ and $\sigma^{2}_{t}$ are the radial and tangential
phase-space averaged square velocity components. Typically, single component
models have been fount to be (numerically) stable for $\xi\lesssim 1.7$ (e.g.
see Nipoti et al. 2002) while some authors put a larger threshold for
stability at $\xi\lesssim 2.3$ (see Meza & Zamorano 1997, see also Maréchal &
Perez 2011).
In stellar systems with a DM halo, the latter has been often indicated as
having a stabilizing effect against ROI, however numerical experiments seem to
weaken this hypothesis (Stiavelli & Sparke 1991; Nipoti et al. 2011), while in
general, even if the total amount of anisotropy is large, an extended
isotropic core does indeed reduce the ROI, as shown by Trenti & Bertin (2006),
and a locally fluctuating halo potential has in principle different effects on
the ignition of the ROI at fixed $\xi$ for different baryon profiles, as
conjectured by Di Cintio & Casetti (2021).
In this work, I investigate by means of collisionless $N-$body simulations in
Newtonian and MOND gravities the implications of the dissipationless collapse
in the emergence of the $M/L$ \- $\epsilon$ relation. The paper is structured
as follows. In Section 2 I discuss the initial conditions, the numerical codes
and the analysis of the simulation outputs. In Section 3 the results of the
simulations are presented and interpreted. Section 4 draws the conclusions.

##  2 Methods and models

###  2.1 Properties of initial conditions

In this work two different sets of numerical experiments have been performed,
cold dissipationless collapses from spherically symmetric initial states (for
baryons and DM, when present) and cold dissipationeless collapses of clumpy
baryon distributions in spherically symmetric DM halos and without DM in MOND.
In both sets of simulations, in the Newtonian models the halos have been
implemented either with particles or with a smooth density distribution
exerting a fixed potential.
The positions for the particles of the given component in spherical systems
were sampled from the so called $\gamma-$models (Dehnen 1993, see also
Tremaine et al. 1994) defined by the density profile

$\rho_{i}(r)=\frac{3-\gamma}{4\pi}\frac{M_{i}r_{c,i}}{r^{\gamma}(r+r_{c,i})^{4-\gamma}},$
(5)

with total mass $M_{i}$, scale radius $r_{c,i}$ and logarithmic density slope
$\gamma$. In this work I discuss mainly the $\gamma=0$ and $\gamma=1$ (i.e.
the Hernquist 1990 model) cases, corresponding to a flat-cored and a
moderately cuspy distribution, respectively. Using the density distribution
(5), instead of the widely adopted NFW (or its alternatives, such as for
example the Einasto 1965 profile, that can model a finite cusp or core, see
also Retana-Montenegro et al. 2012), keeps the number model parameters limited
to the density slope $\gamma$ one the mass ratio $M/M_{*}$ and scale radii are
fixed, while preserving both the ”cosmological” central density
$\rho_{0}(r)\propto 1/r$ and the large radii $1/r^{3}$ fall off.
Similarly to previous studies (e.g. see Londrillo et al. 2003; Nipoti et al.
2006; Di Cintio et al. 2013 and references therein) baryonic particle
velocities were sampled from a position independent isotropic Maxwell-
Boltzmann distribution and then normalized in order to get the wanted value of
the virial ratio in the range $0\leq 2K/|W|\leq 0.2$. When considering a
virialized live halo with (initial) density $\rho_{DM}$, the DM particles’
velocities are obtained sampling with the usual rejection method the ergodic
phase-space distribution function $f(\mathcal{E})$ evaluated numerically with
the standard Eddington (1916) inversion as

$f(\mathcal{E})=\frac{1}{\sqrt{8}\pi^{2}}\int_{\mathcal{E}}^{0}\frac{{\rm
d}^{2}\rho_{DM}}{{\rm d}\Phi^{2}}\frac{{\rm d}\Phi}{\sqrt{\Phi-\mathcal{E}}},$
(6)

where, $\mathcal{E}=v^{2}/2-\Phi(r)$ is the specific energy for unit mass and
the total potential $\Phi=\Phi_{*}+\Phi_{\rm DM}$, where the stellar and DM
potentials $\Phi_{*}$ and $\Phi_{\rm DM}$ are given for a $\gamma$-model by

$\Phi_{i}(r)=\frac{GM_{i}}{r_{c,{i}}(2-\gamma)}\left[\left(\frac{r}{r+r_{c,{i}}}\right)^{2-\gamma}-1\right];\quad
i=*,{\rm DM}.$ (7)

Clumpy initial conditions have been implemented following Hansen et al.
(2006). First a root $\gamma=1$ model is generated and then $N_{C}$ clumps
also described by the density (5) with different choices of the scale radius
$r_{c}$, and density slope $\gamma$ are generated with centres having
Poissonian displacements from the root Hernquist model. When present, the DM
halo is initialized as above for the spherical collapse. Again, stellar
particles velocities are selected from a isotropic Maxwellian and later
normalized to obtain the desired value of $2K/|W|$. MOND systems, as they lack
of a dark component, are characterized by the dimensionless $\kappa$
parameter, defined by $\kappa\equiv GM/a_{0}r_{c}^{2}$, where $a_{0}\approx
10^{-8}$cm s-2 is the MOND scale acceleration (Milgrom 1983), so that for
$\kappa\gg 1$ one recovers the Newtonian regime, while for $\kappa\lesssim 1$
the system is in MONDian regime.
Throughout this work a constant mass to light ratio for the stellar component
is adopted so that, in computer units $M_{*}/L=1$ for all models. As a
consequence, Equation (1) is hereafter expressed in terms of $M/M_{*}$.

###  2.2 Numerical codes

The Newtonian $N-$body simulations discussed here have been performed with the
publicly available fvfps code (fortran version of a fast Poisson solver,
Londrillo et al. 2003), using a parallel version of the classical Barnes & Hut
(1986) tree scheme for the force evaluation combined with the Dehnen (2002)
fast multipole method (see also Dehnen & Reed 2011). The simulations span a
range of $N$ (accounting for both stellar and DM components) between $10^{4}$
and $10^{6}$. In the lowest resolution cases (here $N=10^{4}$) only stars are
simulated with particles and the DM halo (when present) is modelled as an
external fixed potential (cfr. Eq. 7).
The gravitational forces are smoothed below a cut-off distance given by the
so-called softening length (see Dehnen 2001) in the range
$0.02\leq\epsilon_{\rm soft}\leq 0.05$ in units of the initial $r_{c,*}$.
MOND systems have been simulated using the nmody modified Poisson Solver (for
the details see Londrillo & Nipoti 2011, see also Nipoti et al. 2007) applying
an iterative relaxation procedure on a spherical grid333As a rule, here we
always used a $128\times 32\times 64$ grid, starting from a seed guess
solution (in the same fashion of the standard Newtonian Poisson solvers see
Londrillo & Messina 1990; Londrillo et al. 1991) to evaluate the gravitational
potential from the non-linear Poisson equation (Bekenstein & Milgrom 1984)

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

The MOND interpolation function $\mu(x)$ is assumed hereafter to be

$\mu_{M}(x)=\frac{x}{\sqrt{1+x^{2}}},$ (9)

yielding the usual asymptotic limits

$\mu_{M}(x)\sim\cases{\displaystyle}1,\quad x\gg 1,\\\ \displaystyle x,\quad
x\ll 1;{}$ (10)

that for $||\nabla\Phi||\gg a_{0}$ Eq. (8) essentially recovers the Poisson
Equation of the Newtonian regime, while for $||\nabla\Phi||\ll a_{0}$ the
system is in the so-called deep-MOND (often labelled as dMOND) regime444It is
worth noticing that it is not guaranteed that a given system is at all radii
in dMOND regime, nor that when evolving from an initial condition in such
state, it will continue to be so, Nipoti et al. 2007. with Eq. (8) simplifying
to

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

All simulations were extended up to $t=300t_{\rm Dyn}$ where as usual we
define the half mass dynamical time $t_{\rm Dyn}\equiv\sqrt{2r_{h}^{3}/GM}$,
for which $r_{h}$ is the radius containing half of the total system mass
$M=M_{*}+M_{DM}$ in Newtonian models or $M=M_{*}$ in MOND models and single
component Newtonian models. By doing so, one is ensured that the collective
oscillations are sufficiently damped out and the virial ratio of the bound
matter $2K/|W|\simeq 1$. Particles are propagated with a second order leapfrog
scheme with adaptive time step (see Hut et al. 1995) $\Delta t$ that varies
during each run as $\Delta t\equiv\eta/\sqrt{4\pi G\rho_{\rm max}(t)}$, where
$\rho_{\rm max}(t)$ is the time-dependent maximal density and $\eta$ is the
so-called Courant-Friedrichs-Lewy condition, that in the simulations discussed
here was fixed to 0.3.

Figure 1: Mass ratio $M/M_{*}$ as function of the 2D ellipticity
$\epsilon_{2D}$ measured on three random projections (indicated with
increasingly bigger symbol size) for the end products of models with initially
cold ($K_{0}=0$) Hernquist profiles ($\gamma=1$) with frozen (squares) and
live (circles) Hernquist DM halos and in MOND (diamonds), and as function of
the deprojected ellipticity for different choices of the inclination angle
$\theta$ (dots) an their average value (right panel). For comparison, in both
panels the dashed line and the orange shaded area mark the Deur 2014 relation,
while the vertical red dashed lines indicate the limit ellipticity 0.7. Figure
2: Sérsic index $m$ as function of the total to baryonic matter ratio for the
end products of collapses with frozen and live DM halos, as well as dynamical
to baryonic matter in MOND simulations. The cyan dotted-dashed line marks the
relation emerging from the observational data of Sonnenfeld et al. 2019.
Figure 3: Intrinsic ellipticity $\epsilon$ as function of $M/M_{*}$ for
initially cold spherical systems in live (circles) and frozen (squares) halos
(left panel); and for initially cold clumpy systems in live (downward
triangles) and frozen (upward triangles) halos (right panel). The empty and
filled symbols refer to $\epsilon$ evaluated inside $r_{70\%}$ and $r_{90\%}$,
respectively. Figure 4: Density profiles at $t=300t_{\rm Dyn}$ for halo (left
panel) and baryons (right panel) for $\mu=5$ (squares), 2.5 (circles), 1
(diamonds), 0.167 (triangles) and 0 (downward triangles). The dashed and solid
lines mark the initial profiles ($\gamma=1$ in both cases) for DM and baryons,
respectively. Figure 5: Same as Fig. 3 for MOND systems with $\gamma=1$ (left)
and clumpy (right) initial conditions. Figure 6: Minimum vs maximum axial
ratios for the end products of Newtonian collapses in live (circles) and
frozen (squares) halos, and MOND collapses for the different values of the
mass ratio $M/M_{*}$ (colour bar). Figure 7: Anisotropy parameter inside
$r_{70\%}$ as function of $\epsilon$ for the end products of the cold collapse
simulations. The color maps of the points marks the value of $M/M_{*}$.
Symbols have the same meaning as in Fig. 6.

###  2.3 Structural analysis of the end products

For the sets of simulations introduced above, projected and intrinsic
properties were evaluated in the standard way (see e.g. Nipoti et al. 2006; Di
Cintio et al. 2013 and reference therein). Once the end products are
translated to the centre of mass frame, the particles undergo three random
rotations, the rotated system is then projected in the plane perpendicular to
the three putative lines of sight. The 2D ellipticities are obtained as
$\epsilon_{2D}=1-\varrho$, where in each of the three projections
$\varrho=b_{2D}/a_{2D}$ is the ratio between the minimum and maximum projected
semiaxis.
To have a measure of the concentration of the model, the angle averaged two
dimensional density profiles $\Sigma(R)$ were fitted with the Sérsic (1968)
law

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

In the equation above $\Sigma_{e}$ is the projected mass density at effective
radius $R_{e}$, the radius of the circle containing half of the projected
mass. Since the two dimensionless parameters $b$ and $m$ are related by
$b\simeq 2m-1/3+4/405m$ (see Ciotti & Bertin 1999), Eq.(12) needs effectively
a simple one parameter fit. I recall that, for high values of $m$ the density
profile is steep in the central regions and shallow in the outer, while low
values of $m$ correspond shallow central density profiles with steeper outer
slopes.
The intrinsic triaxiality of the end products is recovered evaluating the
second order tensor555Note that $I_{ij}$ is not the inertia tensor, which is
given instead by ${\rm Tr}(I_{ij})\delta_{ij}-I_{ij}$.

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

for the particles inside the spheres of radius $r_{50}$ , $r_{70}$ and
$r_{90}$ containing the 50%; 70% and the 90% of the stellar mass of the
system, respectively.
The matrix $I_{ij}$ is diagonalized iteratively, with the requirement that the
percentage difference of the largest eigenvalue between two consequent
iterations does not exceed $10^{-3}$. On average, for $N\approx 10^{5}$ the
process requires about 10 iterations. Once the three eigenvalues $I_{1}\geq
I_{2}\geq I_{3}$ are recovered, a rotation is applied to the system such that
the three eigenvectors oriented along the coordinate axes. For of a
heterogeneous ellipsoid of semiaxes $a,b$ and $c$, one has $I_{1}=Aa^{2}$,
$I_{2}=Ab^{2}$ and $I_{3}=Ac^{2}$, where $A$ is a constant depending on the
density profile. The axial ratios are defined by $b/a=\sqrt{I_{2}/I_{1}}$ and
$c/a=\sqrt{I_{3}/I_{1}}$, so that the ellipticities in the principal planes
are $\epsilon_{1}=1-\sqrt{I_{2}/I_{1}}$ and
$\epsilon_{2}=1-\sqrt{I_{3}/I_{1}}$. Models with $c/a\sim b/a\lesssim 0.5$ are
defined prolate while models for which $c/a\sim b/a\gtrsim 0.5$ are defined
oblate, while models having $b/a>0.5$ and $c/a<0.5$ result evidently triaxial.
For the end products of the MOND simulations, the (bona fide) ENS is recovered
with the same procedure of Re & Di Cintio (2023) by evaluating the angle
averaged (spherical) density profile on a radial grid from which the Newtonian
and MOND force fields $\mathbf{g}_{N}$ and $\mathbf{g}_{M}$ are recovered. The
DM density of the ENS is then obtained as

$\tilde{\rho}_{DM}=(4\pi G)^{-1}\nabla\cdot(\mathbf{g}_{M}-\mathbf{g}_{N}).$
(14)

Note that, Equation (14) above, is valid only for isolated spherical systems,
as substituting the source term in Eq. (8) with the Poisson Equation
$\Delta\Phi_{N}=4\pi G\rho$ and integrating out the divergence term on both
sides yields in any other geometry

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

where $\mathbf{S}\equiv\nabla\times\mathbf{h}(\rho)$ is a density-dependent
solenoidal field.
The DM content $M_{DM}$ of the ENS is finally obtained by integrating Eq. (14)
from 0 up to the radius of the farmost particle so that the total dynamical
mass of the model is again $M=M_{*}+M_{DM}$.

##  3 Simulations and results

###  3.1 Projected properties

Figure 1 (left panel) shows for three random projections (indicated with
differently sized symbols) the 2D ellipticites of the end products of
spherical collapses in frozen (squares) and live (circles) halos as well as
baryon-only MOND (diamonds) collapses versus the mass ratio on a log scale. In
all cases both baryons and DM (when present) start from Hernquist density
profiles with the same scale radius $r_{s}$ and $N=10^{5}$. The stellar
systems produced in collapses within frozen halos show a somewhat increasing
trend of $M/M_{*}$ with $\epsilon_{2D}$ (or vice versa), though not
reproducing the linear relation and its associated uncertainty, given in Eq.
(1) and marked in figure by the dashed line and the (orange) shaded area.
Collapses in initially virialized live halos, seem to produce systems with
lesser departure from the spherical symmetry at low or rather large (up to
$\sim 10^{2}$) $M/M_{*}$, with maximum value of $\epsilon_{2D}$ for
$M/M_{*}\sim 3$. Remarkably, for the random projections shown here, at fixed
mass ratio $\epsilon_{2D}$ is systematically smaller for the products of the
collapses in live halos. MOND collapses show a somewhat intermediate behaviour
attaining large values of the projected ellipticity (around $\epsilon_{2D}\sim
0.6$) for dynamical to baryon mass ratios of order $10^{2}$.
For comparison, in the right panel of Fig. 1 the mass ratio is plotted against
the estimated intrinsic minimal ellipticity $\epsilon_{3D}$ recovered
inverting Eq. (2) for the intermediate value of $\epsilon_{2D}$ in the
assumption of oblateness for a Gaussian distribution of the assumed projection
angle $\theta$ (points). The mean value of $\epsilon_{3D}$ for 30 independent
realizations of $\theta$ is marked by the symbols, with the same coding as in
the left panel. Again, the dashed curve and the orange shaded area highlight
Eq. (1). If on one hand the point distributions presents an increasing trend
of $M/M_{*}$ with $\epsilon_{3D}$, on the other, the deprojected ellipticities
again fail to reproduce the linear relation observed by Deur (2014) and
Winters et al. (2023). It is worth noting though, that, some MOND systems and
Newtonian systems with a frozen halo are intrinsically prolate (see discussion
below), thus contradicting the oblateness assumption in applying Eq. (2).
Final states attained starting by different sets of Newtonian initial
conditions (i.e. clumpy, $\gamma=0$ flat cored initial baryon profiles, not
shown here) do not present radically different projected ellipticity trends.
Figure 2 presents the Sersić index $m$ as a function of $M/M_{*}$. Most
models, independently on the specific value of the mass ratio have indexes in
the range $1.1\lesssim m\lesssim 7$, compatible with the values for the
observed elliptical galaxies (e.g. see Zahid & Geller 2017; Sonnenfeld et al.
2019; Sonnenfeld 2019 and references therein). Collapses of initially
spherical distributions in frozen halos, produce remarkably large Sersić
indexes at high $M/M_{*}$, with values up to $\sim 9$ and 12.3 for initially
cuspy $\gamma=1$ or flat cored $\gamma=0$ baryon profiles, respectively. For
$M/M_{*}\lesssim 40$ the $\gamma=1$ initial condition in a frozen halo produce
final values of $m\sim 4$, corresponding to a de Vaucouleurs (1948) profile.
The Sérsic indexes of MOND models starting from both spherical (diamonds) and
clumpy (pentagons) initial conditions do not show a clear trend with the ratio
of dynamical to stellar mass, with initially spherical systems yielding a
narrower span of values around $m\sim 2$. Notably, MOND clumpy initial
conditions could produce extremely centrally shallow final projected profiles
with $m$ down to $\approx 0.89$ for $M/M_{*}\approx 50$ (corresponding to
$\kappa\sim 10$). Sonnenfeld et al. (2019) found for their sample of spheroids
that $m\propto M_{*}^{0.46}$ and $M_{DM}\propto M_{*}^{1.7}$, that would
correspond roughly to $m\propto(M/M_{*}-1)^{0.66}$, indicated in Figure by the
dotted-dashed line. Unfortunately, no set of simulations appears to reproduce
such trend of the Sérsic index, except perhaps the spherical collapses in
frozen halos whose end products exhibit and increasing trend of $m$ with
$M/M_{*}$.

###  3.2 Intrinsic properties

As in the $N-$body simulations discussed here the ratio of the total dynamical
mass to stellar mass is essentially the control parameter, it is convenient to
show the intrinsic minimum ellipticity $\epsilon$ as function of $M/M_{*}$. In
Fig. 3 $\epsilon$ is plotted against $M/M_{*}$ for Newtonian simulations with
spherical ($\gamma=1$, left panel) and clumpy (right panel) initial conditions
for the baryons. In all cases the DM halo is, at least initially, simulated
with a $\gamma=1$ model.
Initially spherical systems collapsing in a frozen halo with central
$\rho_{DM}\propto 1/r$ (indicated in figure by empty and filled squares) relax
to flatter end states for increasing $M/M_{*}$, eventually tending to
$\epsilon\approx 0.8$. For mass ratios $M/M_{*}\gtrsim 5$ such end states are
significantly more flattened than a E7 galaxy (indicated by the dashed line at
$\epsilon=0.7$) both at the Lagrangian radii enclosing 70% (empty symbols) and
90% (filled symbols) of $M_{*}$. When the halo is live (i.e., modelled using
particles, circles in figure) the picture is rather different and the
ellipticity shows a non monotonic trend with the mass ratio $M/M_{*}$.
Curiously, $\epsilon$ starts to decrease with at around $M/M_{*}\approx 5$,
value for which $\epsilon(M/M_{*})$ settles to a somewhat constant value for
the frozen halo simulations. Spherical collapses in cored halos (i.e. with
$\gamma=0$, not shown here) have essentially the same behaviour as their
counterparts with a cusp, though for the cases where said halo is frozen, the
limit value of $\epsilon$ at increasing $M/M_{*}$ is at around 0.7.
Clumpy initial conditions in spherical halos yield end products that have a
qualitative trend of $\epsilon$ with $M/M_{*}$ as their spherical counterparts
when the DM halo is frozen (though with a larger scatter of $\epsilon$ at
bigger mass ratios). Vice versa, a live halo almost always induces mildly
flattened end states, around $\epsilon\approx 0.3$ for $M/M_{*}\gtrsim 6$
while a somewhat non-monotonic trend is evident a low mass ratios.

Figure 8: Evolution of the anisotropy parameter $\xi$ (top right panel) and
axial ratios $c/a$ (bottom left) and $c/a$ (bottom right) for different
Newtonian models in frozen halos (solid lines) and MOND models (dashed lines)

The final DM halos’ ellipticities $\epsilon_{DM}$ have been evaluated for all
live halo simulations. In general, for $N_{*}=N_{DM}=5\times 10^{4}$ one has
$0.93\lesssim\epsilon_{DM}\lesssim 0.98$ implying that the collapse of the
baryon distribution does not alter significantly the sphericity of the
virialized halo, while its central regions become significantly more ”cored”
at later times when $M/M_{*}\lesssim 6$. This could be in principle a
numerical collisionality artifact induced by the simulation resolution (i.e.
number of particles used for DM and stars and particle mass ratio
$\mu=m_{DM}/m_{*}$) as shown in Fig. 4 where the final DM (left panel) and
stellar (right panel) density profiles are plotted for different realizations
of a model with initial $\gamma=1$ halo and stars density profiles with
$M/M_{*}=5$ and various choices of particle resolution ranging from $\mu=5$ to
0.167 (i.e. the individual mass of DM particles spans from values larger than
that of particles representing the stellar component to significantly lower).
In both sets of density profiles the curves differ significantly from one
another at radii smaller than $r/r_{c,*}\sim 10^{-2}$. In particular, the halo
density $\rho_{DM}$ becomes increasingly cored for higher DM resolutions666In
principle, the DM resolution-relate issues could be probed using the effective
multi-component models introduced by Nipoti et al. 2021. (i.e. larger $N_{DM}$
and and smaller $\mu$ at fixed $M/M_{*}$). Surprisingly, in that case, the
associated $\rho_{*}$ is strikingly similar to that obtained for the frozen
halo simulation (downward and upward triangles in the right panel). A similar
behaviour has been verified for lower values of $M/M_{*}$ and clumpy initial
conditions for the initial stellar distribution.
Figure 5 shows $\epsilon$ as function of $M/M_{*}$ in MOND simulations where
the effective total dynamical mass $M$ is obtained as discussed in Sect. 2.3.
Spherical collapses starting from cuspy initial conditions with $\gamma=1$ and
various values of the parameter $\kappa$ always produce relaxed end states
with $0.4\lesssim\epsilon\lesssim 0.7$. Between $M/M_{*}=1$ and $\approx 14$,
$\epsilon$ has a markedly increasing trend, though nowhere linear as claimed
in Deur (2014) and Winters et al. (2023) and observed in the MOND $N-$body
simulations by Re & Di Cintio (2023) limited to $\kappa=1$ and 100 cases only,
but with different values of the initial virial ratio $2K/|W|$. Clumpy systems
have a similar behaviour but with the increasing trend of $\epsilon$ breaking
at a lower value of $M/M_{*}\approx 5$, corresponding to simulations with an
initial value of $\kappa\approx 500$. Curiously, for spherical and clumpy
initial conditions, the minimal ellipticity is obtained in correspondence of
the trend inversion.
Newtonian and MOND systems starting from similar initial conditions have
qualitatively different intrinsic triaxiality as summarized in Fig. 6 where
the final ratio of the minimum to maximum semiaxis $c/a$ is plotted against
the ratio of the intermediate to maximum semiaxis $b/a$. In general, as
observed also in Nipoti et al. 2007, 2011, single component Newtonian
collapses mostly produce marginally oblate systems while the end products of
MOND collapses are in general slightly prolate or triaxial. Here we observe
that when a (frozen) DM halo is added (squares in figure), Newtonian end
states often become markedly prolate, in particular for $M/M_{*}>20$ (as
colour coded in Fig. 6), in particular when $c/a$ falls below 0.3 (as
indicated by the red dashed line), corresponding to ellipticities larger than
the threshold value of $\epsilon=0.7$ for an E7 galaxy. Newtonian collapses in
live halos (cfr. circles) mostly produce oblate (at large $M/M_{*}$) or mildly
triaxial systems (for low values of $M/M_{*}$). The MOND simulations performed
here typically produce oblate ens states for effective $M/M_{*}\lesssim 6$ and
markedly prolate end states for $6\lesssim M/M_{*}\lesssim 35$, larger values
of the effective mass ratio (associated to values of $\kappa<10$ in the
initial condition) are generally producing strongly triaxial systems, without
any evident trend with the cored, cuspy or clumpy nature of the initial mass
distribution.
Independently on the specific gravity model at hand (Newtonian or MOND), as a
general trend the simulations performed in this work yield more and more
radially anisotropic states for larger values of $\epsilon$, up to $\xi\sim
26.5$ for Newtonian models collapsing in frozen halos. Typically, increasing
values of $M/M_{*}$ reflect in larger final $\xi$. In Figure 7 (cfr. also Fig.
7 in Re & Di Cintio 2023) the anisotropy parameter $\xi$ is shown as function
of $\epsilon$ for various choices of $M/M_{*}$ indicated by the color map.
MOND models stand out as being able to attain rather large values of $\xi$ for
intermediate values of $\epsilon$ (around $0.45$) and small values of their
effective dynamical to stellar mass ratio $M/M_{*}$. As shown by the time
evolution of $\xi$ (top panel) and the axial ratios (bottom panels) of Fig. 8,
MOND models (red and purple dashed lines) reaching comparable values of the
anisotropy index to that of Newtonian systems become considerably more
flattened in lesser time in units of their dynamical time. As expected, in
both paradigms of gravity, initially cold spherical models ($2K/|W|\ll 1$),
become rapidly (i.e. below $1t_{Dyn}$) radially anisotropic and undergo a
process akin to the ROI for unstable equilibrium systems and once relaxed
appear triaxial or flattened. In the Newtonian cases, the bigger the mass
ratio $M/M_{*}$ the higher the initial value of $\xi$ (cfr. top panel in Fig.
8). This leads to conjecture that a (almost) spherical pre-formed DM halo
induces strongly radially unstable initial states for the stellar/baryonic
component, at least within an interval of $M/M_{*}$. For fixed mass ratios and
halo density profile, the end products of cold clumpy initial conditiosn are
systematically less anisotropic than their initially spherical counterparts.
In MOND collapses, on the contrary, clumpy initial conditions tend to produce
systems with larger values of $\xi$ than those starting spherical, with
comparable ”phantom” DM evaluated from Eq. (14), reaching similar final values
of $\epsilon$.

##  4 Summary and discussions

Aiming at shedding some light on the origin the ellipticity - dynamical mass
relation I have performed $N-$body simulations of dissipationless collapse in
both Newtonian gravity with dark matter and MOND. The Newtonian and MOND
simulations presented in the previous Section point towards the fact that, at
least for values of the mass ratio $M/M_{*}$ (intrinsic or effective in the
case of MOND) between 5 and 6, a certain increasing trend between the latter
and the ellipticity of the spheroid is present. The linear proportionality
(cfr. Eq. 1) discussed by Deur (2014) and later by Winters et al. (2023) could
not be recovered in any of the two paradigms of gravity considered here,
neither for the intrinsic ellipticities nor the deprojected values estimated
with Eq. (2).
In particular, in Newtonian simulations with a frozen DM halo the values of
$\epsilon$ at the Lagrangian radii enclosing the 90% of the stellar mass
”saturate” at around 0.75 for $M/M_{*}$ larger than 6. For the Newtonian
simulations where a live DM halo was considered, an unequivocal non-monotonic
relation between $\epsilon$ and $M/M_{*}$ (and hence $M/L$ if a constant
$L/M_{*}$ is assumed) is observed. In particular, a tendency to produce less
flattened end states when $M/M_{*}$ exceeds $\approx 20$ is observed for
spherical initial conditions. It is therefore tantalizing to infer that at
even larger mass ratios the end products of Newtonian collapses should
essentially be almost spherical, as one would expect for the case of
ultrafaint dwarf galaxies where $M/M_{*}$ may exceed $10^{2}$ (e.g. Zoutendijk
et al. 2021). One should however bear in mind that for a given mass ratio, the
end products of live halo simulations could be influenced by the resolution,
in general, for spherical collapses for increasing resolution (i.e.,
decreasing values of $\mu$) the final models tend to depart more from the
spherical symmetry. MOND models interpreted in the context of a DM scenario
(i.e. when the halo of the associated ENS is accounted) also present an
increasing trend of $\epsilon$ with $M/M_{*}$ below $M/M_{*}\approx 10$, in
partial contrast with the simulations of Re & Di Cintio (2023), limited to a
single value of $\kappa=GM_{*}/a_{0}^{2}r_{s}^{2}$, yielding a monotonic and
quasi-linear trend between the quantities.
Prompted by recent results on the relation between stellar and halo masses and
the Sérsic index of the stellar component (see Sonnenfeld et al. 2019), $m$
was evaluated for the Newtonian and MOND simulation outputs. Remarkably, MOND
collapses (both clumpy and spherical) are found to be able to produce systems
with $m$ order unity (corresponding to an extremely shallow core) without
invoking any dissipative mechanism, at variance with the Newtonian simulations
of Nipoti (2015), where clumpy initial conditions where used with the power-
spectrum index $n$ as control parameter. Clumpy initial conditions in
Newtonian simulations with frozen halos also can also attain low values of
$m$, however in such cases this is likely a numerical artifact. In general,
independently on the specific nature of the initial baryon distribution, live
DM halos flatten down to $\epsilon\sim 0.06$ (retaining a oblate 3D structure)
and have their central cusp significantly lowered for clumpy initial
conditions on the stellar mass, thus qualitatively confirming what noted in
the two component Newtonian simulations of Cole et al. (2011), see also
Pascale et al. (2023).
To summarize, the above listed results point toward the fact that in a certain
span of halo masses in units of the visible matter total mass $M_{*}$ the
increasing trend between $M/M_{*}$ and the flattening $\epsilon$ could be a
consequence of the dissipationless collapse. The fact that the linear relation
previously discussed in the literature is not recovered, can be likely
ascribed to the uncertainty on the measurements of the deprojection procedure
itself. From the point of view of modified Newtonian dynamics (MOND),
evaluating the dark mass content of the ENS of the end states of MONDian
simulations reveals that a relation akin to Eq. (1) could be also supported in
a modified dynamics scenario. Moreover in MOND sufficiently clumped initial
conditions can yield flattened and centrally shallow spheroids even in absence
of dissipation, as opposed to Newtoian collapses with DM requiring some form
of dissipation.

###### Acknowledgements.

I would like to express gratitude to Federico Re, Carlo Nipoti and Stefano
Zibetti for the useful discussions at an early stage of this work. I am also
wishing 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

  * Alabi et al. (2017) Alabi, A. B., Forbes, D. A., Romanowsky, A. J., et al. 2017, MNRAS, 468, 3949
  * Auger et al. (2010) Auger, M. W., Treu, T., Bolton, A. S., et al. 2010, ApJ, 724, 511
  * Bacon et al. (1985) Bacon, R., Monnet, G., & Simien, F. 1985, A&A, 152, 315
  * Baldi & Salucci (2012) Baldi, M. & Salucci, P. 2012, J. Cosmology Astropart. Phys., 2012, 014
  * Barnes et al. (2009) Barnes, E. I., Lanzel, P. A., & Williams, L. L. R. 2009, ApJ, 704, 372
  * Barnes & Hut (1986) Barnes, J. & Hut, P. 1986, Nature, 324, 446
  * Bekenstein & Milgrom (1984) Bekenstein, J. & Milgrom, M. 1984, ApJ, 286, 7
  * Bellovary et al. (2008) Bellovary, J. M., Dalcanton, J. J., Babul, A., et al. 2008, ApJ, 685, 739
  * Benhaiem et al. (2016) Benhaiem, D., Joyce, M., Sylos Labini, F., & Worrakitpoonpon, T. 2016, A&A, 585, A139
  * Bertin (2014) Bertin, G. 2014, Dynamics of Galaxies (Cambridge, UK: Cambridge University Press)
  * Bertin et al. (1994) Bertin, G., Pegoraro, F., Rubini, F., & Vesperini, E. 1994, ApJ, 434, 94
  * Bertola et al. (1991) Bertola, F., Bettoni, D., Danziger, J., et al. 1991, ApJ, 373, 369
  * Bertola et al. (1993) Bertola, F., Pizzella, A., Persic, M., & Salucci, P. 1993, ApJ, 416, L45
  * Bett (2012) Bett, P. 2012, MNRAS, 420, 3303
  * Binney & Tremaine (2008) Binney, J. & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton University Press)
  * Bolton et al. (2008) Bolton, A. S., Burles, S., Koopmans, L. V. E., et al. 2008, ApJ, 682, 964
  * Brook et al. (2014) Brook, C. B., Di Cintio, A., Knebe, A., et al. 2014, ApJ, 784, L14
  * Burger & Zavala (2021) Burger, J. D. & Zavala, J. 2021, ApJ, 921, 126
  * Burger et al. (2022) Burger, J. D., Zavala, J., Sales, L. V., et al. 2022, MNRAS, 513, 3458
  * Cappellari (2008) Cappellari, M. 2008, MNRAS, 390, 71
  * Chen & Hwang (2020) Chen, C.-Y. & Hwang, C.-Y. 2020, ApJ, 903, 38
  * Ciotti (1999) Ciotti, L. 1999, ApJ, 520, 574
  * Ciotti & Bertin (1999) Ciotti, L. & Bertin, G. 1999, A&A, 352, 447
  * Ciotti & Pellegrini (1992) Ciotti, L. & Pellegrini, S. 1992, MNRAS, 255, 561
  * Cole et al. (2011) Cole, D. R., Dehnen, W., & Wilkinson, M. I. 2011, MNRAS, 416, 1118
  * de Vaucouleurs (1948) de Vaucouleurs, G. 1948, Annales d’Astrophysique, 11, 247
  * Dehnen (1993) Dehnen, W. 1993, MNRAS, 265, 250
  * Dehnen (2001) Dehnen, W. 2001, Mon. Not. R. Astron. Soc., 324, 273
  * Dehnen (2002) Dehnen, W. 2002, Journal of Computational Physics, 179, 27
  * Dehnen & Reed (2011) Dehnen, W. & Reed, J. 2011, 2011, Eur. Phys. J. Plus 126, 55 [arXiv:1105.1082]
  * 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 & Casetti (2020) Di Cintio, P. & Casetti, L. 2020, MNRAS, 494, 1027
  * Di Cintio & Casetti (2021) Di Cintio, P. & Casetti, L. 2021, in IAUS proceedings series, Vol. 364, Multi-scale (time and mass) Dynamics of Space Objects, ed. A. Celletti, C. Beaugé, C. Gales, & A. Lemâitre
  * 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
  * Eddington (1916) Eddington, A. S. 1916, MNRAS, 76, 572
  * Einasto (1965) Einasto, J. 1965, Trudy Astrofizicheskogo Instituta Alma-Ata, 5, 87
  * Eriksen et al. (2021) Eriksen, M. H., Frandsen, M. T., & From, M. H. 2021, A&A, 656, A123
  * Evslin (2016) Evslin, J. 2016, ApJ, 826, L23
  * Fridman & Poliachenko (1984) Fridman, A. M. & Poliachenko, V. L. 1984, Physics of gravitating systems. II - Nonlinear collective processes: Nonlinear waves, solitons, collisionless shocks, turbulence. Astrophysical applications (Springer, Berlin)
  * Gajda et al. (2015) Gajda, G., Łokas, E. L., & Wojtak, R. 2015, MNRAS, 447, 97
  * Genina et al. (2018) Genina, A., Benítez-Llambay, A., Frenk, C. S., et al. 2018, MNRAS, 474, 1398
  * Gonzalez et al. (2022) Gonzalez, E. J., Hoffmann, K., Gaztañaga, E., et al. 2022, MNRAS, 517, 4827
  * Hansen et al. (2006) Hansen, S., Moore, B., Zemp, M., & Stadel, J. 2006, JCAP, 0601, 014
  * Harris et al. (2019) Harris, G. L. H., Babyk, I. V., Harris, W. E., & McNamara, B. R. 2019, ApJ, 887, 259
  * Hernquist (1990) Hernquist, L. 1990, ApJ, 356, 359
  * Hut et al. (1995) Hut, P., Makino, J., & McMillan, S. 1995, ApJ, 443, L93
  * Jin et al. (2020) Jin, Y., Zhu, L., Long, R. J., et al. 2020, MNRAS, 491, 1690
  * Jing & Suto (2002) Jing, Y. P. & Suto, Y. 2002, ApJ, 574, 538
  * Joyce et al. (2009) Joyce, M., Marcos, B., & Sylos Labini, F. 2009, Mon. Not. R. Astron. Soc., 397, 775
  * Kahlhoefer et al. (2019) Kahlhoefer, F., Kaplinghat, M., Slatyer, T. R., & Wu, C.-L. 2019, J. Cosmology Astropart. Phys., 2019, 010
  * Katz et al. (2018) Katz, H., Desmond, H., Lelli, F., et al. 2018, MNRAS, 480, 4287
  * 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
  * Londrillo et al. (2003) Londrillo, P., Nipoti, C., & Ciotti, L. 2003, Memorie della Societa Astronomica Italiana Supplementi, 1, 18
  * Maréchal & Perez (2010) Maréchal, L. & Perez, J. 2010, MNRAS, 405, 2785
  * Maréchal & Perez (2011) Maréchal, L. & Perez, J. 2011, Transport Theory and Statistical Physics, 40, 425
  * Mateo (1998) Mateo, M. L. 1998, ARA&A, 36, 435
  * McGaugh et al. (2003) McGaugh, S. S., Barker, M. K., & de Blok, W. J. G. 2003, ApJ, 584, 566
  * Merritt & Aguilar (1985) Merritt, D. & Aguilar, L. A. 1985, Mon. Not. R. Ast. Soc, 217, 787
  * Meza & Zamorano (1997) Meza, A. & Zamorano, N. 1997, ApJ, 490, 136
  * Milgrom (1983) Milgrom, M. 1983, ApJ, 270, 365
  * Navarro et al. (1997) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, Astrophys. J., 490, 493
  * Nipoti (2015) Nipoti, C. 2015, ApJ, 805, L16
  * Nipoti et al. (2021) Nipoti, C., Cherchi, G., Iorio, G., & Calura, F. 2021, MNRAS, 503, 4221
  * Nipoti et al. (2011) Nipoti, C., Ciotti, L., & Londrillo, P. 2011, MNRAS, 414, 3298
  * Nipoti et al. (2002) Nipoti, C., Londrillo, P., & Ciotti, L. 2002, MNRAS, 332, 901
  * Nipoti et al. (2006) Nipoti, C., Londrillo, P., & Ciotti, L. 2006, MNRAS, 370, 681
  * Nipoti et al. (2007) Nipoti, C., Londrillo, P., & Ciotti, L. 2007, ApJ, 660, 256
  * Oh et al. (2015) Oh, S.-H., Hunter, D. A., Brinks, E., et al. 2015, The Astronomical Journal, 149, 180
  * Palmer & Papaloizou (1987) Palmer, P. L. & Papaloizou, J. 1987, MNRAS, 224, 1043
  * Pascale et al. (2023) Pascale, R., Calura, F., Lupi, A., et al. 2023, MNRAS, 526, 1428
  * Pechetti et al. (2017) Pechetti, R., Seth, A., Cappellari, M., et al. 2017, ApJ, 850, 15
  * Pizzella et al. (1997) Pizzella, A., Amico, P., Bertola, F., et al. 1997, A&A, 323, 349
  * Polyachenko & Shukhman (1981) Polyachenko, V. L. & Shukhman, I. G. 1981, Soviet Ast., 25, 533
  * Ragone-Figueroa et al. (2010) Ragone-Figueroa, C., Plionis, M., Merchán, M., Gottlöber, S., & Yepes, G. 2010, MNRAS, 407, 581
  * Re & Di Cintio (2023) Re, F. & Di Cintio, P. 2023, A&A, 678, A110
  * Retana-Montenegro et al. (2012) Retana-Montenegro, E., van Hese, E., Gentile, G., Baes, M., & Frutos-Alfaro, F. 2012, A&A, 540, A70
  * Sánchez Almeida (2022) Sánchez Almeida, J. 2022, ApJ, 940, 46
  * Sánchez Almeida & Trujillo (2021) Sánchez Almeida, J. & Trujillo, I. 2021, MNRAS, 504, 2832
  * Schneider et al. (2012) Schneider, M. D., Frenk, C. S., & Cole, S. 2012, J. Cosmology Astropart. Phys., 2012, 030
  * Sérsic (1968) Sérsic, J. L. 1968, Atlas de Galaxias Australes (Cordoba, Argentina: Observatorio Astronomico)
  * Shu et al. (2017) Shu, Y., Brownstein, J. R., Bolton, A. S., et al. 2017, ApJ, 851, 48
  * Simon (2019) Simon, J. D. 2019, ARA&A, 57, 375
  * Sonnenfeld (2019) Sonnenfeld, A. 2019, in Linking Galaxies from the Epoch of Initial Star Formation to Today, 64
  * Sonnenfeld et al. (2019) Sonnenfeld, A., Wang, W., & Bahcall, N. 2019, A&A, 622, A30
  * Stiavelli & Sparke (1991) Stiavelli, M. & Sparke, L. S. 1991, ApJ, 382, 466
  * Sylos Labini et al. (2015) Sylos Labini, F., Benhaiem, D., & Joyce, M. 2015, MNRAS, 449, 4458
  * Tremaine et al. (1994) Tremaine, S., Richstone, D. O., Byun, Y.-I., et al. 1994, AJ, 107, 634
  * Trenti & Bertin (2006) Trenti, M. & Bertin, G. 2006, ApJ, 637, 717
  * Valenzuela et al. (2007) Valenzuela, O., Rhee, G., Klypin, A., et al. 2007, ApJ, 657, 773
  * Vera-Ciro et al. (2011) Vera-Ciro, C. A., Sales, L. V., Helmi, A., et al. 2011, MNRAS, 416, 1377
  * Weinberg (2023) Weinberg, M. D. 2023, MNRAS, 525, 4962
  * Winters et al. (2023) Winters, D. M., Deur, A., & Zheng, X. 2023, MNRAS, 518, 2845
  * Worrakitpoonpon (2015) Worrakitpoonpon, T. 2015, MNRAS, 446, 1335
  * Zahid & Geller (2017) Zahid, H. J. & Geller, M. J. 2017, ApJ, 841, 32
  * Zoutendijk et al. (2021) Zoutendijk, S. L., Brinchmann, J., Bouché, N. F., et al. 2021, A&A, 651, A80