Dissipationless collapse and the dynamical mass-ellipticity relation of elliptical galaxies in Newtonian gravity and MOND

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 +/- 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\'ersic 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.


Introduction
In the ΛCDM scenario, stellar systems such as galaxies and clusters are thought to be embedded in dark matter (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, Navarro-Frenk-White;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 disc galaxies, the interplay between the DM density distribution ρ DM and that of the stellar component ρ * has attracted a lot of interest, in particular with respect to the widely studied cusp-core problem (Moore 1994).That is, 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 → 0) distributions predicted by DM N-body simulations, such as the NFW profile (e.g.see Oh et al. 2015).In addition, it can also be proved analytically that flat cored models cannot 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 for the core-cusp problem have been proposed, ranging 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, implying that the flattening of the central DM density distribution is due the intrinsic nature of the DM (e.g. the so-called self-interacting DM; Baldi & Salucci 2012;Kahlhoefer et al. 2019;Sánchez Almeida & Trujillo 2021;Burger et al. 2022).Some authors also suggested 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, alternatives to the DM hypothesis, such as the modified Newtonian dynamics (MOND; Milgrom 1983), explain the core-cusp problem as an artefact of the gravity model (see e.g.Eriksen et al. 2021;Sánchez Almeida 2022;Re & Di Cintio 2023).
Although 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, much less is known about the relationship between the DM profile and the intrinsic shape of the stellar component, namely its ellipticity .Deur (2014Deur ( , 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.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 disc dynamics, Bertola et al. 1991Bertola et al. , 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), observed that the two quantities are related by the linear relation M/L = (14.1 ± 5.4) . (1) In the equation above, the mass-to-light ratio is normalised such that M/L( 2D = 0.3) ≡ 8 M /L ≡ 4M/M * ( 2D = 0.3), and the intrinsic ellipticity is inferred from its apparent 2D-projected value 2D in the assumption of oblateness and a Gaussian distribution of projection angles θ from Equation ( 1) implies that (at least for the galaxies examined) a larger contribution of the DM mass M DM to the total mass M corresponds to a larger departure from 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 formation1 , requiring that the baryons aggregate and form structure in dark matter halos that had collapsed at earlier times.However, notably, 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 characterised by a rather spherical shape ( ≤ 0.1) while being DM dominated with M/L up to 102 or more.A factor governing the mass ratioellipticity relation could be related to the departure from spherical symmetry of the halos themselves and its relation to their total DM content.Indeed, 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 et al. 2015;Evslin 2016;Gonzalez et al. 2022) point out that more massive halos have significantly higher departure from 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) to the end products of cold gravitational collapses are consistent with Eq. ( 1) if the initial conditions are sufficiently clumpy.The origin of the triaxiality of (single-component) gravitational systems, forming via direct collapse from nearly spherical initial conditions, is usually ascribed to the process of radial-orbit instability 2 (ROI; see Polyachenko & Shukhman 1981;Palmer & Papaloizou 1987;Bertin et al. 1994;Maréchal & Perez 2010), where an initially minor departure from spherical symmetry grows by accreting particles on lowangular-momentum orbits (Allen et al. 1990).This explanation is supported by 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 models with a larger degree of radial anisotropy3 , which is quantified (see e.g.Binney & Tremaine 2008, see also Fridman & Poliachenko 1984) by the parameter where K r and K t = K θ + K φ are the radial and tangential components of the kinetic energy tensor, respectively, given by where σ 2 r and σ 2 t are the radial and tangential phase-space averaged square velocity components.Typically, single-component models have been found to be (numerically) stable for ξ 1.7 (e.g.see Nipoti et al. 2002), while some authors put a larger threshold for stability at ξ 2.3 (see Meza & Zamorano 1997, see also Maréchal & Perez 2011).
In stellar systems with a DM halo, the latter has often been indicated as having a stabilising effect against ROI, however numerical experiments seem to weaken this hypothesis (Stiavelli & Sparke 1991;Nipoti et al. 2011).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).Moreover, a locally fluctuating halo potential has, in principle, different effects on the ignition of the ROI at fixed ξ for different baryon profiles, as conjectured by Di Cintio & Casetti (2021).
Simple N-body simulations of dissipationless collapse, merger, or instability growth in anisotropic equilibria have been rather successful in reproducing the scaling laws of elliptical galaxies and the fundamental plane itself (see e.g. the work by Nipoti et al. 2002Nipoti et al. , 2003Nipoti et al. , 2006)).It is therefore natural to ask whether highly idealised collisionless numerical simulations could reproduce Eq. (1).In this work, I use collisionless N-body simulations in Newtonian and MOND gravities to investigate the implications of dissipationless collapse in the emergence of the M/L-relation.The paper is structured as follows.In Sect.2, I discuss the initial conditions, the numerical codes, and the analysis of the simulation outputs.In Sect.3, the results of the simulations are presented and interpreted.In Sect.4, I draw conclusions.

Properties of initial conditions
The numerical simulations of dissipationless collapse presented in this paper (in both Newtonian and MOND paradigms) are characterised by two types of initial conditions, namely spherically symmetric initial states or clumpy baryon distributions.In all Newtonian simulations, when present, the DM halos are always initially spherical.Such halos have been implemented either with particles or with a smooth non-evolving density distribution, exerting a fixed potential.All simulated systems are assumed to be isolated so that any effect of the environment in their shape can be excluded, and in accordance with the explicit exclusion of group galaxies in the samples analysed by Deur (2014) and Winters et al. (2023).
The positions for the particles of the given component (baryons or DM) in spherical systems were sampled from the socalled γ-models (Dehnen 1993, see also Tremaine et al. 1994), which are defined by the density profile The density ρ i above is therefore defined by the total mass M i , scale radius r c,i , and logarithmic density slope γ.In this work, I discuss mainly the γ = 0 and γ = 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 the Einasto 1965 profile, which can model a finite cusp or core; see also Retana-Montenegro et al. 2012), keeps the number of model parameters limited to the density slope γ , the mass ratio M/M * , and the scale radii, which are fixed, while preserving both the "cosmological" central density ρ 0 (r) ∝ 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 positionindependent isotropic Maxwell-Boltzmann distribution and were then normalised in order to get the wanted value of the virial ratio in the range 0 ≤ 2K/|W| ≤ 0.2.The simulations shown here are limited to 2K/|W| = 10 −3 .When considering a virialised live halo with (initial) density ρ DM , the velocities of the DM particles are obtained by sampling with the usual rejection method, that is, the ergodic phase-space distribution function f (E) evaluated numerically with the standard Eddington (1916) inversion as where, E = v 2 /2 − Φ(r) is the specific energy for unit mass and the total potential Φ = Φ * + Φ DM , where the stellar and DM potentials Φ * and Φ DM are given for a γ-model by Clumpy initial conditions have been implemented following Hansen et al. (2006).First a root γ = 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 γ are generated with centres having Poissonian displacements from the root Hernquist model.When present, the DM halo is initialised as above for the spherical collapse.Again, stellar particle velocities are selected from an isotropic Maxwellian and are later normalised to obtain the desired value of 2K/|W|.MOND systems, as they lack a DM component, are characterised by the dimensionless κ parameter, defined by κ ≡ GM/a 0 r 2 c , where a 0 ≈ 10 −8 cm s −2 is the MOND scale acceleration (Milgrom 1983), which means that, for κ 1, one recovers the Newtonian regime, while for κ 1 the system is in the 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, Eq. ( 1) is hereafter expressed in terms of M/M * .All properties of the initial conditions presented in this work are summarised in Table 1.

Numerical codes
The Newtonian N-body simulations discussed here were 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 = N * + N DM (i.e. the total number of simulation particles) between 10 4 and 10 6 (see Cols. 2  and 3 in Table 1 below).In the lowest-resolution cases (here N = 104 ), only stars are simulated with particles and the DM halo (when present) is always modelled as an external fixed potential (see 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 ≤ soft ≤ 0.05 in units of the initial r c, * .
MOND systems were 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 grid 4 , starting from a seed guess solution (in the same fashion as 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): The MOND interpolation function µ(x) is assumed hereafter to be yielding the usual asymptotic limits which for ||∇Φ|| a 0 Eq. ( 8) essentially recover the Poisson equation of the Newtonian regime, while for ||∇Φ|| a 0 the system is in the so-called deep-MOND (often labelled as dMOND) regime5 with Eq. ( 8) simplifying to All simulations were extended up to t = 300t Dyn , where, as above, we define the half mass dynamical time, 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, it is ensured that the collective oscillations are sufficiently damped out and the virial ratio of the bound matter 2K/|W| 1.

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 references therein).Once the end products are translated to the centre-of-mass frame, the particles undergo three random rotations, and the rotated system is then projected in the plane perpendicular to the three putative lines of sight.The 2D ellipticities are obtained as 2D = 1− , where in each of the three projections, = 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 Σ(R) were fitted with the Sérsic (1968) law where Σ e is the projected mass density at effective radius R e , which is the radius of the circle containing half of the projected mass.As the two dimensionless parameters b and m are related by b 2m − 1/3 + 4/405m (see Ciotti & Bertin 1999), Eq. ( 13) can be fitted with the single parameter m.I reiterate 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 to shallow central density profiles with steeper outer slopes.The intrinsic triaxiality of the end products is recovered by evaluating the second-order tensor6 for the particles inside the spheres of radius r 50 , r 70 , and r 90 containing 50%, 70%, and 90% of the stellar mass of the system, respectively.
The matrix I i j is diagonalised 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 ≈ 10 5 , the process requires about ten iterations.Once the three eigenvalues I 1 ≥ I 2 ≥ I 3 are recovered, a rotation is applied to the system such that the three eigenvectors are oriented along the coordinate axes.For 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 0.5 are defined as oblate; while models with b/a > 0.5 and c/a < 0.5 are found to be triaxial.
For the end products of the MOND simulations, the (bona fide) ENS is recovered with the same procedure as that used by Re & Di Cintio (2023), that is, by evaluating the angle-averaged (spherical) density profile on a radial grid from which the Newtonian and MOND force fields g N and g M are recovered.The DM density of the ENS is then obtained as I note that Eq. ( 15) above is only valid for isolated spherical systems, as in any other geometry, substituting the source term in Eq. ( 8) with the Poisson equation ∆Φ N = 4πGρ and integrating out the divergence term on both sides yields where S ≡ ∇ × h(ρ) is a density-dependent solenoidal field.
The DM content M DM of the ENS is finally obtained by integrating Eq. ( 15) from 0 up to the radius of the farmost particle so that the total dynamical mass of the model is again

Projected properties
For three random projections (indicated with differently sized symbols), Fig. 1 (left panel) shows 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 2D (or vice versa), though the linear relation and its associated uncertainty -given in Eq. ( 1) and marked in the figure by the dashed line and the (orange) shaded area -is not produced.Collapses in initially virialized live halos seem to produce systems with a smaller departure from the spherical symmetry at low or rather large (up to ∼10 2 ) M/M * , with a maximum value of 2D for M/M * ∼ 3. Remarkably, for the random projections shown here, at fixed mass ratio, 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 2D ∼ 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 3D recovered by inverting Eq. ( 2) for the intermediate value of 2D whilst assuming oblateness for a Gaussian distribution of the assumed projection angle θ (points).The mean value of 3D for 30 independent realisations of θ 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).Although the point distributions present an increasing trend of M/M * with 3D , the deprojected ellipticities still fail to reproduce the linear relation observed by Deur (2014) and Winters et al. (2023).However, it is worth noting that some MOND systems and Newtonian systems with a frozen halo are intrinsically prolate (see discussion below), which contradicts the oblateness assumption in applying Eq. ( 2).
The final states attained after starting from different sets of Newtonian initial conditions (i.e.clumpy, γ = 0 flat cored A254, page 5 of 11 initial baryon profiles, not shown here) do not show 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 of 1.1 m 7, which is 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 of up to ∼9 and 12.3 for initially cuspy γ = 1 or flat cored γ = 0 baryon profiles, respectively.For M/M * 40, the γ = 1 initial condition in a frozen halo produces final values of m ∼ 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 ∼ 2. Notably, MOND clumpy initial conditions could produce extremely centrally shallow final projected profiles with m down to ≈0.89 for M/M * ≈ 50 (corresponding to κ ∼ 10).Sonnenfeld et al. (2019) found for their sample of spheroids that m ∝ M 0.46 * and M DM ∝ M 1.7 * , which roughly correspond to m ∝ (M/M * − 1) 0.66 , indicated in Fig. 2 by the dotted-dashed line.With the notable exception of the spherical collapses in frozen halos, whose end products exhibit an increasing trend of m with M/M * , the other sets of simulations fall outside the trend of the Sérsic index with M/M * when the latter exceeds ≈30.Notably, M/M * ≈ 30 is also the maximum value of the mass ratio explored in the samples used by Sonnenfeld et al. (2019).

Intrinsic properties
As in the N-body simulations discussed here the total dynamicalto-stellar mass ratio is essentially the control parameter, it is convenient to show the intrinsic minimum ellipticity as a function of M/M * .In Fig. 3, is plotted against M/M * for Newtonian simulations with spherical (γ = 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 γ = 1 model.
Initially spherical systems collapsing in a frozen halo with central ρ DM ∝ 1/r (indicated in figure by empty and filled squares) relax to flatter end states for increasing M/M * , eventually tending to ≈ 0.8.For mass ratios M/M * 5, such end states are significantly more flattened than an E7 galaxy (indicated by the dashed line at = 0.7) at the Lagrangian radii enclosing both 70% (empty symbols) and 90% (filled symbols) of M * .When the halo is live (i.e.modelled using particles, shown by the circles in Fig. 3) the picture is rather different and the ellipticity shows a non-monotonic trend with the mass ratio M/M * .Curiously, starts to decrease at around M/M * ≈ 5, the value for which (M/M * ) settles to a somewhat constant value for the frozen halo simulations.Spherical collapses in cored halos (i.e. with γ = 0, not shown here) show essentially the same behaviour as their counterparts with a cusp; although for the cases where the halo is frozen, the limit value of at increasing M/M * is at around 0.7.
Clumpy initial conditions in spherical halos yield end products that show a qualitative trend of with M/M * , similar to their spherical counterparts when the DM halo is frozen (though with a larger scatter of at bigger mass ratios).Vice versa, a live halo almost always induces mildly flattened end states, around ≈ 0.3 for M/M * 6, while a somewhat non-monotonic trend is evident at low mass ratios.for the end products of collapses with frozen and live DM halos, as well as dynamical to baryonic matter in MOND simulations for the simulations with parameters given in rows (7, 9, 18, 19, and 22) of Table 1.The cyan dot-dashed line marks the relation emerging from the observational data of Sonnenfeld et al. (2019).
The final DM halo ellipticities DM have been evaluated for all live halo simulations.In general, for N * = N DM = 5 × 10 4 one has 0.93 DM 0.98, implying that the collapse of the baryon distribution does not significantly alter the sphericity of the virialised halo, while its central regions become significantly more "cored" at later times, when M/M * 6.This could, in principle, be a numerical collisionality artefact induced by the simulation resolution (i.e.number of particles used for DM and stars, and particle mass ratio µ = m DM /m * ).For example, Fig. 4 shows the final DM (left panel) and stellar (right panel) density profiles for different realisations of a model with initial γ = 1 halo and star density profiles with M/M * = 5 and various choices of particle resolution ranging from µ = 5 to 0.167 (i.e. the individual mass of DM particles spans from values larger than that of the 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, * ∼ 10 −2 .In particular, the halo density ρ DM becomes increasingly cored for higher DM resolutions 7 (i.e.larger N DM and smaller µ at fixed M/M * ).Surprisingly, in that case, the associated ρ * is strikingly similar to that obtained for the frozen halo simulation (downward-and upwardpointing triangles in the right panel).We see similar behaviour for lower values of M/M * and clumpy initial conditions for the initial stellar distribution.
Figure 5 shows as a function of M/M * in MOND simulations (corresponding to rows 22 and 23 in Table 1) where the effective total dynamical mass M is obtained as discussed in Sect.Frozen halo (r 70% ) Frozen halo (r 90% ) Live halo (r 70% ) Live halo (r 90% ) Fig. 3. Intrinsic ellipticity as a 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 evaluated inside r 70% and r 90% , respectively.The simulation parameters are summarised in rows (14-17) of Table 1.Baryons (t=0) µ=0 (Frozen halo) Fig. 4. Density profiles at t = 300t Dyn for halo (left panel) and baryons (right panel) for = 5 (squares), 2.5 (circles), 1 (diamonds), 0.167 (triangles), and 0 (downward triangles).The dashed and solid lines mark the initial profiles (γ = 1 in both cases) for DM and baryons, respectively.The simulation parameters are summarised in row (20) of Table 1.
1 and ≈14, 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).However, the latter simulations were limited to κ = 1 and 100 cases only, but with different values of the initial virial ratio 2K/|W|.Clumpy systems (as for their Newtonian counterparts) have a rather flat trend of with a few outliers for values of the mass ratio smaller than M/M * ≈ 5, corresponding to simulations with an initial value of κ ≈ 500.Curiously, for spherical 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 summarised 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 also observed in Nipoti et al. (2007Nipoti et al. ( , 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 Fig. 6), 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 = 0.7 for an E7 galaxy.Newtonian collapses in live halos (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 end states Fig. 7. Anisotropy parameter inside r 70% as a function of for the end products of the simulations.As in Fig. 6, the colour code maps the (logarithmic) mass ratio M/M * listed in Col. 4 of Table 1.
for effective M/M * 6 and markedly prolate end states for 6 M/M * 35.Larger values of the effective mass ratio (associated to values of κ < 10 in the initial condition) generally produce strongly triaxial systems, without any evident differences between cored, cuspy, and clumpy initial mass distributions.
Independently of the specific gravity model used (Newtonian or MOND), as a general trend, the simulations performed in this work yield increasingly radially anisotropic states for larger values of , up to ξ ∼ 26.5 for Newtonian models collapsing in frozen halos.Typically, increasing values  1) and MOND models (dashed lines).
of M/M * lead to larger final ξ.In Fig. 7 (cf.also Fig. 7 in Re & Di Cintio 2023), the anisotropy parameter ξ is shown as function of for various choices of M/M * , as indicated by the colour map.MOND models stand out as being able to attain rather large values of ξ for intermediate values of (around 0.45) and small values of their effective dynamical-to-stellar mass ratio, M/M * .As shown by the time evolution of ξ (top panel) and the axial ratios (bottom panels) of Fig. 8, for γ = 1 systems, MOND models (red and purple dashed lines) reaching comparable values of the anisotropy index to that of Newtonian systems with a frozen halo8 (line 7 in Table 1) become considerably more flattened over a shorter period of dynamical time.We note that, for a given stellar mass M * , MOND systems have dynamical timescales that are significantly different from those of comparable Newtonian models with a DM halo (cf.Eq. ( 12)), with MONDian collective timescales being typically longer (Nipoti et al. 2007).I also note that MOND systems attaining similar values of ξ to the Newtonian models with M/M * = 2 and 1 have slightly larger dynamical-to-stellar mass ratios (≈4 and 1.7, respectively).It should also be borne in mind that, when converted into physical units, the timescales of Fig. 8 may differ significantly from one another.For example, using the halo-to-stellar-mass relation of Sonnenfeld et al. (2019) (see also Tinker et al. 2017), with a scale ratio of 100 kpc, the dynamical time for Newtonian (baryons plus DM) systems would range from ≈1.3 Gyr down to ≈0.06 Gyr for mass ratios in the range of 2 M/M * 30, spanning an interval of luminous masses of between 10 10 and 10 12 M .
As expected, in both paradigms of gravity explored here, initially cold spherical models (2K/|W| 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 higher the mass ratio M/M * , the higher the initial (and final) value of ξ (top panel in Fig. 8).This leads to the suggestion that a spherical or almost spherical preformed DM halo induces strongly radially unstable initial states for the stellar or baryonic component, at least within an interval of M/M * .For fixed mass ratios and a fixed halo density profile, the end products of cold clumpy A254, page 9 of 11 initial conditions 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 ξ than those that are spherical to begin with, with comparable "phantom" DM evaluated from Eq. ( 15), reaching similar final values of .

Summary and discussions
With an aim to shed light on the origin of the ellipticitydynamical mass relation, I performed N-body simulations of dissipationless collapse in both Newtonian gravity with dark matter and MOND.In general, the products of most numerical experiments, in agreement with previous studies with analogous initial conditions, have Sérsic index values, radial anisotropy, and density profiles in ranges comparable to those of the observed elliptical galaxies.
The Newtonian and MOND simulations discussed here both show that, at least for values of the mass ratio M/M * (intrinsic or effective in the case of MOND) of between 5 and 6, disipationless collapse alone produces a certain increasing trend with the ellipticity of the spheroid.However, the linear proportionality (cf.Eq. ( 1)) discussed by Deur (2014) and later by Winters et al. (2023) could not be recovered in either of the two paradigms of gravity considered here, that is, with intrinsic ellipticities or their deprojected values estimated with Eq. ( 2).
In particular, in Newtonian simulations with a frozen DM halo, the values of at the Lagrangian radii enclosing 90% of the stellar mass "saturate" at around 0.75 for M/M * larger than 6.Of course, this is most likely a numerical effect due to the rather unphysical fixed halo assumption.For the Newtonian simulations where a live DM halo was considered, an unequivocal non-monotonic relation between 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 ≈20 is observed for spherical initial conditions.It is therefore tempting to hypothesise that at even larger mass ratios the end products of Newtonian collapses should 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).However, we must be aware that for a given mass ratio, the end products of live halo simulations could be influenced by resolution.Indeed, with increasing resolution (i.e.decreasing values of µ) when regarding spherical collapses, the final models tend to depart more from spherical symmetry.MOND models interpreted in the context of a DM scenario (i.e. when the halo of the associated ENS is accounted for) also present an increasing trend of with M/M * below M/M * ≈ 10; this is in partial contrast with the simulations of Re & Di Cintio (2023), which are limited to a single value of κ = GM * /a 2 0 r 2 s , yielding a monotonic and quasilinear 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, which is at variance with the Newtonian simulations of Nipoti (2015), where clumpy initial conditions were used with the power-spectrum index n as a control parameter.Clumpy initial conditions in Newtonian simulations with frozen halos can also attain low values of m; however, in such cases this is likely a numerical artefact.Independently of the specific nature of the initial baryon distribution, live DM halos flatten down to ∼ 0.06 (retaining an oblate 3D structure) and have their central cusp significantly lowered for systems starting with initial conditions characterized by a clumpy stellar component, which qualitatively confirms what has been noted for the two-component Newtonian simulations of Cole et al. (2011); see also Pascale et al. (2023).The increasing trend of m with M/M * deducible from Sonnenfeld et al. (2019) is observed only below M/M * ∼ 30.However, systems with larger dynamical-to-luminous matter ratios were not in the samples analysed by these latter authors.
In summary, the results listed above point towards 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 could be a consequence of the dissipationless collapse.The fact that an unequivocal linear relation is not recovered can be ascribed to the uncertainty on the measurements or the deprojection procedure for , and/or to the strong assumptions made on the mass-to-light ratio of the visible matter.Moreover, the non-monotonicity in the trend between and M/M * in some simulation setups is not at all in contrast with the observed ellipticals.The qualitative interpretation that emerges is that, the deeper the (total) potential well during the collapses, the stronger the radial orbit instability and therefore the stronger the departure of the final state from spherical symmetry.For extremely large M/M * , the baryonic matter behaves essentially as tracer particles within the halo potential.If the latter has a rather spherical shape, the stellar component barely departs from it, as no ROI takes place.
Most importantly, one has to bear in mind that Eq. ( 1) is defined for an ensemble that explicitly excludes certain systems, such as ultrafaint dwarfs, which by definition fall outside an increasing trend of M/M * with .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 also be supported in a modified dynamics scenario.As a by-product of this work, it appears that in MOND, sufficiently clumped initial conditions can yield flattened and centrally shallow spheroids even in the absence of dissipation, as opposed to Newtonian collapses with DM requiring some form of dissipation.
A natural follow-up to this study, in the context of Newtonian gravity, would focus necessarily on the effects of an initially nonvirialised halo as well as clumpy DM distributions, which are not considered here for reasons of simplicity.For what concerns MOND (and possibly other alternative theories without DM), the interplay between the clumpyness of the initial (baryonic) density and the initial acceleration regime (i.e.how below or above a 0 ) should also be explored.

Fig. 1 .
Fig.1.Mass ratio M/M * as a function of the 2D ellipticity 2D measured on three random projections (indicated with increasingly large symbol size) for the end products of models with initially cold (K 0 = 0) Hernquist profiles (γ = 1) with frozen (squares) and live (circles) Hernquist DM halos and in MOND (diamonds) (see rows 7, 9, and 22 in Table1), and as a function of the deprojected ellipticity for different choices of the inclination angle θ (dots) and 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.ratios are defined by b/a = √ I 2 /I 1 and c/a = √ I 3 /I 1 , meaning that the ellipticities in the principal planes are 1 = 1 − √ I 2 /I 1 and 2 = 1 − √ I 3 /I 1 .Models with c/a ∼ b/a 0.5 are defined as prolate; models for which c/a ∼ b/a 0.5 are defined as oblate; while models with b/a > 0.5 and c/a < 0.5 are found to be triaxial.For the end products of the MOND simulations, the (bona fide) ENS is recovered with the same procedure as that used byRe & Di Cintio (2023), that is, by evaluating the angle-averaged (spherical) density profile on a radial grid from which the Newtonian and MOND force fields g N and g M are recovered.The DM density of the ENS is then obtained as

Fig. 2 .
Fig.2.Sérsic index m as a 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 for the simulations with parameters given in rows(7, 9, 18, 19, and 22)  of Table1.The cyan dot-dashed line marks the relation emerging from the observational data ofSonnenfeld et al. (2019).
2.3.Spherical collapses starting from cuspy initial conditions with γ = 1 and various values of the parameter κ always produce relaxed end states with 0.4 0.7.Between M/M * =

Fig. 8 .
Fig.8.Evolution of the anisotropy parameter ξ (top right panel) and axial ratios c/a (bottom left) and c/a (bottom right) for different Newtonian models in frozen halos (solid lines, row 7 in Table1) and MOND models (dashed lines).

Table 1 .
Summary of the simulation properties.