Structure of the equivalent Newtonian systems in MOND N-body simulations Density profiles and the core-cusp problem

Aims. We investigate the core-cusp problem of the Λ cold dark matter ( Λ CDM) scenario in the context of Modified Newtonian Dynamics (MOND) paradigm exploiting the concept of equivalent Newtonian system (ENS) Methods. By means of particle-mesh N − body simulations in MOND we explore processes of galaxy formation via cold dissipation-less collapse or merging of smaller substructures. From the end states of our simulations, we recover the associated ENS and study the properties of their dark matter halos. We compare the simulation results with simple analytical estimates with a family of γ − models. Results. We find that the dark matter density of ENSs of most spherical cold collapses have a markedly cored structure, in particular for the lowest values of the initial virial ratios. End states of some simulations with clumpy initial conditions have more complex profiles and some of their ENSs exhibit a moderate cusp, with logarithmic density slope always shallower than 1. Conclusions. These results seem to point towards the fact that the absence in most observed galaxies of a central DM cusp, at variance with what one would expect from theoretical and numerical arguments in Λ CDM, would be totally consistent in a MONDian


Introduction
In the Λ cold dark matter scenario (hereafter ΛCDM), theoretical arguments and collisionless N−body simulations (Navarro et al. 1997) predict that galaxies are embedded in dark matter (DM) halos characterized by a ρ(r) ∝ r −1 central cusp.Observational results seem to suggest, from the analysis of the central velocity dispersion profiles of dwarf galaxies, that the DM distribution has a cored density distribution 1 (see Moore 1994;Di Cintio et al. 2014).
Several solutions to this apparent contradiction -often referred to as "the core-cusp problem"-such as self-interacting DM (e.g.see Lovell et al. 2012;Nguyen et al. 2021;Eckert et al. 2022), DM annihilation (e.g.see Vasiliev 2007), baryon feedback (e.g.see Governato et al. 2010;Cole et al. 2011;Pontzen & Governato 2012;Del Popolo & Pace 2016) or simply a misinterpretation of the observational data (McGaugh et al. 2003), have been proposed so far.However, notwithstanding the great 1 Technically speaking, in multi component equilibrium self gravitating systems there exist analytical constraints on the magnitude of a component's density given the logarithmic slope of the other, e.g.see Dubinski & Carlberg 1991;Ciotti & Pellegrini 1992;Ciotti 1996Ciotti , 1999.Moreover for a broad range of spherical density profiles, the central density slope constraints the value of the central anisotropy profile (An & Evans 2006) amount of theoretical and observational work, a clear answer is still far from being obtained.Moreover, given the large interest in alternative theories of gravity, such as among the others f (R) gravities (Buchdahl 1970;Sotiriou & Faraoni 2010); Modified Gravity (MoG, Moffat 2006;Moffat & Rahvar 2013); Retarded Gravity (Raju 2012;Yahalom 2022); Emergent Gravity (Verlinde 2011(Verlinde , 2017)); Refracted Gravity (Cesare et al. 2020;Sanna et al. 2023); Fractional Gravity (Giusti 2020;Benetti et al. 2023), proposed to avoid introducing the DM as a collisionless fluid of exotic particles, it is natural to ask what becomes of the core-cusp problem in those proposals.
In this work we investigate this matter in the modified Newtonian dynamics (hereafter MOND, Milgrom 1983) paradigm.We recall that in the Bekenstein & Milgrom (1984) Lagrangian formulation of MOND (sometimes referred to as AQUAL) the classical Poisson equation for a density-potential pair (ρ; Φ) ∆Φ = 4πGρ (1) is substituted by the non-linear field equation (2) In the equation above where a 0 ≈ 10 −8 cm s −2 is a scale acceleration and µ(x) is the MOND interpolating (monotonic) function Fig. 1.Ratio of the stellar to dark density in the ENS (top) and DM and stellar density profiles (bottom left and bottom right) in units of 3M/4πr 3 c for κ = 100 and γ = 0, 0.5, 1, 1.5, 2 and 2.5.
known only by its asymptotic limits (3) so that for ||∇Φ|| ≫ a 0 Eq. ( 2) one recovers the Newtonian regime, while for ||∇Φ|| ≪ a 0 one obtains the so-called deep-MOND (hereafter dMOND) regime and Eq. ( 2) simplifies to Note that, the non-linear operator in Equation ( 4) is the special case of the p−Laplace operator (see e.g.Stein 1970) for p = 3, while Eq.(1) would correspond to the p = 2 case.In this respect, Equation (2) somewhat "interpolates" between the two regimes via the µ function.Note also that, in both cases, any given baryonic mass density ρ can be taken out from Equation (1) obtaining the relation between the MOND and Newtonian force fields g M and g N , and where S ≡ ∇ × h(ρ) is a density-dependent solenoidal field.It can be proved that the latter is identically null for systems in spherical, cylindrical or planar symmetry, while it is generally non-zero for arbitrary configurations of mass.On which extent the stellar system at hand with mass M is dominated by MOND effects, is usually quantified by the dimensionless parameter where r c is the scale of the baryon distribution.That is, for κ ≫ 1 the system is mainly in Newtonian regime, vice versa for κ ≤ 1 MOND effects become strong at all scales.
For any given stationary model in MOND one can always build the equivalent Newtonian system (hereafter ENS), defined as a system with the same baryonic mass density ρ * plus a DM halo with density ρ DM such that their total potential Φ satisfying Eq. ( 1) is the same as the MOND potential entering Eq. ( 2) for the sole density ρ * .We note that, in principle, the positivity of the DM density of the ENS is not always assured (see Milgrom 1986), in particular for flattened systems (see Ciotti et al. 2006Ciotti et al. , 2012;;Ko 2016).We recall that Milgrom (2010) introduced a Quasi-linear formulation of MOND (hereafter QuMOND) where the modified field equation has the same form as Eq. ( 2), with ν(||g N /a 0 ||) in lieu of ν(||g M /a 0 ||).The QuMOND interpolating function ν(y) can be recovered from µ(x) appearing in Eq. ( 2) as It is easy to show that, from a given baryonic density distribution, one obtains the MONDian potential Φ = Φ N +Φ pDM by first solving a classical Poisson equation for the Newtonian potential Φ N , that trough an algebraic passage involving ν becomes the source for the potential Φ pDM of the so-called "phantom Dark Matter" via a second application of the Poisson equation.Notably, in this alternative bi-potential MOND formulation, the DM is de facto interpreted as the effect of the second potential.As in AQUAL, in QuMOND one can retrieve a dMOND regime, that for spherical systems easily reads g M = a 0 /g N g N (see Milgrom 2021).
If on one hand several workers have investigated the differences between static equilibrium models in Newtonian and MOND gravities, or the interpretation of observations in both theories, on the other much less is known about the formation and evolution of stellar systems.For obvious reasons, in observed systems one has access only to de-projected properties for both stellar and dark components in the Newtonian framework.Numerical experiments, though with their intrinsic limitations, yield information on the full phase-space of the simulated models; in particular the 3D density profiles.In this paper, we explore the structure of the ENS of MOND N−body simulations of galaxy formation, in order to shed some light on the possibility that the core-cusp problem is a MOND artifact in this paradigm of gravity.We stress the fact that the MOND core-cusp problem discussed here, is different from that introduced recently by Eriksen et al. (2021) that deals with the modified gravity versus modified inertia hypothesis (see Milgrom 2022).
The rest of this paper is structured as follows.In Sect. 2 we revise the definition of ENS and discuss their properties.In Sect. 3 we introduce the numerical models and the analysis of the simulations.In Sect. 4 we discuss the properties of the simulations' end states.Finally, in Sect. 5 we summarize and discusses the implications and the relations to previous work.

Equivalent Newtonian systems
As anticipated above, the ENS of a MOND model is the Newtonian system with the same stellar (baryonic) mass distribution ρ * with an additional dark component ρ DM such that the total potential (and thus the force field) is the same of the parent MOND system (see Sanders & Begeman 1994;Angus et al. 2006).For the case of an isolated spherical system one has since the solenoidal term S vanishes.We stress the fact that, Equation (5) in QuMOND can be rewritten exactly as If Equation ( 9) is applied to a spherically symmetric system one has ν(y) = x/y, and the total density of its ENS (baryonic plus phantom DM, see e.g.Hodson et al. 2020;Oria et al. 2021) becomes Let us consider the family of spherical γ−models (Dehnen 1993;Tremaine et al. 1994), with density profile given by where M is the total baryonic mass, 0 ≤ γ < 3 is the logarithmic density slope and r c the scale radius.
If the density profile (11) is substituted in Eq. ( 10) one obtains where with κ defined in Equation ( 6).We note that, for small radii r, Equation ( 13) tends to zero if γ < 1, while it diverges for γ > 1.
In practice, at least for the γ < 1 case, even in central regions the model falls in the MOND regime.Strong MOND corrections in the centre are therefore associated with a dominant DM component ρ DM in the ENS.

Massive galaxies
Let us consider a typical 10 12 M ⊙ massive elliptical galaxy with a scale radius of 3 kpc, modelled with a γ-model.In this case κ ≈ 10 2 .Due to discreteness effects of the underlying stellar system Equation ( 11) can be considered reliable until the radius that contains a fraction of roughly 10 −3 of the total mass M (in this case 10 9 M ⊙ , i.e. the typical mass of its central supermassive black hole).The Lagrangian radius enclosing such mass fraction is The region in MOND regime has a far smaller radius, that for γ ≤ 1 is obtained by y(r 10 −3 ) 10 5 r 2 10 −3 /r 2 c , varying between 2 × 10 2 and 2 × 10 3 .And thus, even in the framework of (Qu)MOND the phantom DM halo does not really dominate in the central region for a cored stellar density profile.In Figure 1 we plot for γ = 0, 0.5, 1, 1.5, 2 and 2.5 the ratio of stellar to phantom DM and their respective radial density profiles for κ = 10 2 .We note that, remarkably, models with a strong cusp (i.e.γ > 1) have phantom DM halos in their ENS characterized by a decreasing density inside the scale radius.Vice versa, cored models are associated with ENS having halos with a weak cusp and several slope changes.

Diffuse galaxies
Let us now consider a diffuse galaxy such that κ ∼ 1; so that its central region can fall in the MOND regime, even for radii bigger than r 10 −3 .Typically, this occurs again if γ < 1.We find that lim r→0 y(r) = 0 in the central region, hence ν(y) ∼ y −1/2 .If substituted in ( 12), this yields that is, the phantom DM component dominates also at small radii.In particular, the latter has a central profile given by The equation above is characterized by a weak cusp with a logarithmic density slope α = − 1+γ 2 > −1.For example, for γ = 0, the DM component in the ENS would have a cusp ∝ r −1/2 .We note that this trend is valid for any spherically symmetric stellar distribution with a central core, and not only for the γ = 0 Dehnen model.We note also it always implies a central weak cusp with logarithmic density slope α = −1/2 for the phantom dark matter.This has the interesting astrophysical implication that a galaxy with a cored stellar density profile could be indeed interpreted in the DM scenario as having a cored halo, due to the fact that week cusps can be often mistaken for cores.
For the cases with γ > 1 and κ = 1, for which the gravitational field diverges in the centre, even though the stellar density is diffuse, the ENS is DM-dominated only in the external region.This can be easily checked by substituting the asymptotic behaviour ν(y , and finding from ( 12) that implying a vanishingly small central DM density.This is summarized in Figure 2 where we plot the same quantities as in the previous Fig.but κ = 1.As expected, in the upper plot showing ρ * /ρ DM , the γ = 0 and 0.5 (red and orange lines) are everywhere below 1, i.e. the system is dominated by the phantom DM distribution at all radii.Vice versa, for the γ ≥ 1 cases, the phantom DM of the ENS dominates only in the external regions.We recall that Sánchez Almeida (2022) showed that galaxies with central regions in MOND regime imply ENS characterized by a decreasing baryon density and a cored DM.

Numerical code and initial conditions
The N−body simulations discussed here have been performed with a modified version of the publicly available nmody particlemesh MOND code (Nipoti et al. 2007a, see also Londrillo & Nipoti 2011 for additional technical details).The latter uses a non-linear Poisson solver to compute Φ from Eq. ( 2) on a N r × N ϑ × N φ spherical grid in polar coordinates, using an iterative relaxation procedure starting from a guess solution (here given by Eq. 5 neglecting S ), as for the linear Poisson methods (see Londrillo & Messina 1990;Londrillo et al. 1991).As a rule, in the simulations discussed here we used a 128 × 32 × 64 grid.
In this work we adopt the following form for the interpolation function Alternative choices can also be implement always leading to qualitatively similar end states.
The equation of motion are integrated using a standard 4 th order leapfrog scheme (see e. We performed two sets of numerical simulations with initial conditions defined as follows. In the first, the particles positions were sampled from Equation (11) while in the second, following Hansen et al. (2006), we first distribute according to a Poissonian distribution inside a larger γ model the centres of N C clumps 2 also described by Eq. ( 11) with different choices of r c , and γ and later populate them with particles.
In both cases, the initial particle velocities are extracted from a position-independent isotropic Maxwell-Boltzmann distribution and normalized to obtain the wanted value of the initial virial ratio 2K/|W|, where K is the total kinetic energy and W the virial function, defined for a (finite mass) continuum system of density ρ and potential Φ as We recall that in isolated dMOND systems of finite mass W = −2 GM 3 a 0 /3 is constant (see Nipoti et al. 2007a).Curiously, even in systems of particles interacting with additive 1/r forces with logarithmic potential the virial function is constant (see Di Cintio et al. 2013Cintio et al. , 2017)).The simulations of this work span a range of N between 10 4 and 10 6 .All simulations were extended up to t = 300t Dyn , where t Dyn ≡ 2r 3 h /GM tot and r h is the radius containing half of the total mass of the system M tot , so that virial oscillations and phasemixing are likely to be complete.
Following Ciotti et al. (2007) in some cases we enforce the spherical symmetry during the collapse by propagating particles only using the radial part of the evaluated force field, so that the system behaves effectively as a spherical shell model introduced in Newtonian gravity by Hénon (1964) and used in MOND among the others by Sanders (2008); Malekjani et al.   2 Clumpy initial conditions were also explored in the context of Newtonian simulations by Nipoti (2015) and Ludlow & Angulo (2017) when investigating the relation of the initial density fluctuation power spectrum with the Sérsic index m (see below) conjectured by Cen (2014).( 2009) and by Di Cintio & Ciotti (2011) for systems interacting with 1/r α forces.

Analysis of the end products
For all simulations presented here we first extract the intrinsic properties of the end products from their phase-space positions.We first evaluate the triaxiality of the final particle distribution (see e.g.Nipoti et al. 2006a;Di Cintio et al. 2013 and references therein) by defining the tensor for the particles with positions r i within the Lagrangian radius r 70 containing the 70% of the stellar mass of the system and evaluating with a standard iterative procedure its three eigenvalues I 1 ≥ I 2 ≥ I 3 .By applying a rotation R to all particles of the system so that the three associated eigenvectors are now oriented along the coordinate axes we then get that the three semiaxes a ≥ b ≥ c from I 1 = Aa 2 , I 2 = Ab 2 and I 3 = Ac 2 , where A is a numerical constant depending on the density profile.Finally, we define the axial ratios b/a = √ I 2 /I 1 and c/a = √ I 3 /I 1 , and the ellipticities in the principal planes Following Nipoti et al. (2007a) and Di Cintio et al. ( 2013) we compare the surface density profiles of the end products with the Sersic (1968) law where Σ e is the projected mass density at effective radius R e , the radius of the circle containing half of the projected mass, and the dimensionless parameters b, m are related by b ≃ 2m − 1/3 + 4/405m as found by Ciotti & Bertin (1999).
Once the projected density in the 3 principal planes is circularized over elliptical shells, we determine the corresponding pair (R e , Σ e ) by particle counts (i.e.we are assuming a constant mass to light ratio for each particle), and fit Eq. ( 21) for the three projections.We find that, in general, all 3 sets of (Σ e , R e , m) are rather similar (differing only for less than the 5%), we therefore chose randomly only one.
In addition, for all simulations we also evaluate the so-called anisotropy index (see Binney & Tremaine 2008) defined by where K r and K t = K θ + K ϕ are the radial and tangential components of the kinetic energy tensor, respectively and read In the expressions above, σ 2 r and σ 2 t are the radial and tangential phase-space averaged square velocity components and are obtained for the end products of the simulations by particle counts over radial shells.
For each simulation we recover the (spherical) DM density of the ENS from Equation ( 8) where the Newtonian force field g N has been evaluated and averaged on the radial coordinate.In practice, we are assuming a "sphericized" system.Finally, for the density distribution ρ DM so obtained we evaluate the logarithmic Table 1.Summary of the simulation properties: After the name of each simulation (Col.1), we report the number of particles (Col.2), the gravity law (MOND or dMOND, Col. 3), the initial density profile (Col.4), the initial virial ratio (Col.5) the axial ratios (Cols.6 and 7) , the Sérsic (Col.8), the final anisotropy index (Col.9) and the virial velocity dispersion (Col.10).

Name
Gravity γ 0 =1 clumpy 10 -1 10 0 10 1 r/r 50 10 -1 10 0 10 1 r/r 50 Fig.3. Final (t = 300t Dyn ) density profiles from MOND simulations (top panels) and the DM halo of the ENS (bottom panels) for cored γ 0 = 0 (left), moderately cuspy γ = 1 (centre) and clumpy (right) initial conditions.The increasing initial values of the virial ratio in the models with spherical initial conditions with γ = 0 and 1 is mapped with increasingly lighter tones of blue and green respectively.All clumpy initial conditions start with 2K 0 /|W| 0 = 0.1 density slope α.We find that the profiles of ρ DM are generally well fitted by the empirical law where r α is a scale radius and ρ α is the associated scale density.Equation ( 24) above recovers the 1/r 2 trend of the density of the ENS as predicted by the logarithmic behaviour of the far field MOND potential.The properties of the simulations and their initial conditions are summarized in Tab. 1 below.

Spherical collapses
One of the main motivations of the present work is to establish whether the end products of MOND dissipationless collapses could, in principle, reproduce the structural properties of elliptical galaxies together with their inferred dark halos.Single component Newtonian collapses with spherical initial conditions, are known to produce flatter end states for increasing values of their initial virial ratio (see Nipoti et al. 2006a,b;Di Cintio et al. 2013 and references therein) at fixed initial density profile.We find that, this (partially) holds true for MOND spherical collapses, as shown in Fig. 3 (top left and top mid panels), where we plot the baryon density distribution at 300t Dyn for γ = 0 and 1 and increasing values of the virial ratio with increasingly lighter tones of blue and green in the range 10 −3 ≤ 2K 0 /|W 0 | ≤ 0.5.Using Equation ( 8) for the angle-averaged final density profile on a spherical grid, we evaluated the density distribution of the DM component of the parent Newtonian model (see bottom panels, same figure).We find that, in qualitative agreement with the structural properties of the ENS (see Figs. 1 and 2 in Sect.2), cuspy end systems can be associated with cored or weakly cuspy phantom halos.In general, the end products of spherical collapses have always inner regions that are baryon dominated when building their ENS, even if the initial conditions are such that κ = 1 (in particular for the γ = 0 cases).
Consistently with Nipoti et al. (2007a), we observe that, independently on the specific value of the initial virial ratio, initial conditions characterized by a moderate density cusp (i.e.0.5 ≤ γ ≤ 2) tend to yield end products that are in general oblate (i.e.0.5 ≲ c/a ≲ b/a), as for Newtonian single component collapses.We typically observe major ellipticities up to ∼ 0.63 (corresponding to the gamma1v0b case, see Tab. 1).Remarkably, MOND collapses with cored initial conditions (i.e.γ = 0) evolve into rather prolate end states for 2K 0 /|W 0 | ≳ 0.1, and markedly triaxial end states for lower values of the initial virial ratio.For both cored and moderately cuspy initial conditions, the inner slope α of the DM halo of the ENS, obtained by fitting with Eq. ( 24) increases for increasing values of the baryon initial virial ratio in the MOND simulation, as shown in Fig. 4 (top panel).The best fit Sérsic index m, measuring the concentration of the projected stellar density profile is always in the range 2 ≤ m ≤ 4.5 for both choices of the initial density profile (mid panel, same figure), while the major ellipticity ϵ = 1 − c/a is typically larger when the initial condition has a lower virial ratio, being smaller for larger values of the initial γ at fixed 2K 0 /|W 0 | (bottom panel, same figure).Remarkably no system is found being more flattened than an E7 galaxy.However, as also found by Nipoti et al. (2007a), dMOND collapses may produce even flatter end states as in the case of the gamma1v0dmd run, for which c/a ∼ 0.24 so that ϵ = 0.76.For fixed initial virial ratio, the end states attain larger values of the central virial velocity dispersion σ vir for increasing values of the initial density slope, while the anisotropy index ξ decreases (cfr.1).At fixed initial density profile, the final values of σ vir have little variation with 2K 0 /|W 0 |, while ξ is usually lower for the relaxed states of hotter initial conditions.In order to clarify whether the properties of the halo in the ENS are or not an artifact of the angle averaging procedure, we also performed a set of simulation in enforced spherical symmetry by propagating particles only using the radial component of the force field.By doing so, the system remains spherically symmetric (as no radial orbit instability is possible) and S = 0 de facto hold true at all times, so that one could apply Eq. ( 8) exactly.In Figure 5 we show the same quantities as in Fig. 4 as function of the initial values of γ for systems starting with a virial ratio of 10 −4 with (empty symbols) and without (filled symbols) enforced spherical symmetry.The trend as well as the values of α 3D and effective 1D simulations are comparable, and the same could be noted also for the Sérsic index m that attains considerably lower values (associated with a more concentrated density profile) for larger values of the initial logarithmic density slope.In all cases (cfr.1), as expected, 1D collapses relax to final stated with rather large values of the orbital anisotropy ξ.
Figure 6 shows the final angle averaged density profiles for γ 0 = 0, 1 and 1.5 in 3D and 1D simulations (solid lines) as well as the density profiles of the ENS halos (dashed lines).Notably, if on one hand the large r behaviour of the baryon density profiles ρ * (where the systems are mostly dominated by radial orbits) does not change significantly, on the other, the inner slope of ρ * is always higher for the end products of the 1D simulations and typically settles around 2.5.With the sole exception of the cored initial conditions (γ 0 = 0), the DM halo of the ENS of the end products is denser (in units of the baryon component density ρ * ,50 evaluated at the half mass radius r 50 ) for the 1D simulations, being in both cases considerably shallower than the parent baryon density.

Clumpy collapses
The numerical studies in MOND carried out so far, have typically explored spherical initial conditions (see Nipoti et al. 2007a;Ciotti et al. 2007;Sanders 2008;Malekjani et al. 2009;Nipoti et al. 2011), disks (Brada & Milgrom 1999;Tiret & Combes 2007, 2008a;Nipoti et al. 2007c;Ghafourian & Roshan 2017;Wittenburg et al. 2020) or galaxy merging (Nipoti et al. 2007b;Tiret & Combes 2008b) and references therein.Here, in addition to the usual spherical collapses we also explored clumpy initial conditions.When starting with such initial states, MOND simulations tend (as expected) to yield markedly triaxial end states with broader ranges of both c/a and b/a.In general, for fixed values of the initial virial ratio, the systems tend to relax at later times with respect to their initially spherical counterparts (the oscillations of 2K/W damp out at about 50t Dyn in spherical collapses, see Nipoti et al. 2007a, while in clumpy systems this happens on average at around 140t Dyn ) for analogous choices of the virial ratio.We report here only the runs corresponding to 2K 0 /|W 0 | = 0.1, see Tab. asymptotic standard error of about, on average, 3% as for the spherical collapses, while the scatter in the m Sérsic parameter is slightly smaller for MOND clumpy systems (see top panel in Fig. 7).For comparison, we also run the same clumpy initial conditions in dMOND finding a larger scatter in m.
As a general trend, the DM halo of the circularized ENS of clumpy collapses are significantly more cored3 than what is typically obtained in spherical collapses.In several cases the inner density slopes are negative, down to ∼ −0.99, corresponding to a DM density profile that decreases in the central regions (middle panels in Fig. 7).Interestingly, no initially clumpy system is found to evolve into a state flatter than an E7 galaxy (thin dashed line in bottom panels of Fig. 7) in MOND simulations.However, some dMOND collapses result in considerably flatter end states (and often prolate) with major ellipticity reaching 0.87 for the clumpy1dmd.
We observe that, final states with larger values of the anisotropy index ξ (i.e. more and more dominated by lowangular momentum orbits) are always associated to larger ellipticities ϵ and Sérsic indexes.A similar, though somewhat weaker, correlation is also found between α and ϵ, that could be read in the DM scenario as steeper inner DM profiles producing flatter stellar distributions.
4.3.The MOND mass-to-light ratio -ellipticity relation Deur (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.Jeans anisotropic modelling, gravitational lensing, X-ray spectra and the dynamics of satellite star clusters) and the ellipticity ϵ, recovered the linear relation where the M/L is normalized such that M/L(ϵ app = 0.3) ≡ 8M ⊙ /L ⊙ ≡ 4M/M * (ϵ app = 0.3), and the intrinsic ellipticity ϵ is extrapolated from its observed 2D projected value ϵ app assuming that all systems are oblate with a Gaussian distribution of projection angles θ so that The Equation above in the context of ΛCDM implies that a larger contribution of the DM mass M DM to the total mass M corresponds to a larger departure from the spherical symmetry (here quantified by larger major ellipticity) for the stellar component.Winters et al. (2023) argue that, if true, such a correlation would be contrasting the standard ΛCDM scenario of galaxy formation, where more massive (and rather spherical) DM halos are embed less flattened stellar systems.We note that some peculiar elliptical galaxies (though excluded by the original sample of Winters et al. 2023) such as the ultrafaint dwarfs (Simon 2019) appear to go against the trend given by Eq. ( 25), having usually ϵ ≲ 0.1 with M/L in some cases up to 10 3 .Using the simulations discussed in the previous sections, we have investigated the relation (25) in the context of MOND, evaluating the effective DM mass M DM in the ENSs of both clumpy and spherical collapses.To do so, after recovering the ρ DM from the angle averaged ENS, we integrate it radially up to the radius containing all simulation particles.
In Figure 8 we show the total to stellar mass (here we have assumed units such that M * /L = 1) ratio M/M * versus major ellipticity ϵ for collapses with both spherical and clumpy initial states, here indicated by circles and diamonds respectively, as well as the observational relation given in Eq. ( 25).We found that the end products of initially clumpy systems fall in (almost) We stress the fact, that none of the simulations discussed above produces final states that could be interpreted as ultrafaint dwarfs (except, possibly, some dMOND collapses), that in the standard cosmological scenario are supposed to be DM dominated at all radii (i.e. even in the central region where our simulations, when interpreted in the context of DM have baryon dominated cores).

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

FedericoFig. 4 .
Fig. 4. Inner density slope of the ENS halo (top panel), best fit Sérsic index (middle panel) and minor ellipticity ϵ = 1 − c/a as function of the initial virial ratio for initial conditions with Dehnen profiles with γ = 0 (circles) and 1 (triangles)

Fig. 5 .
Fig. 5. Inner density slope of the ENS halo (top panel), best fit Sérsic index (middle panel) and minor ellipticity ϵ = 1 − c/a as function of the logarithmic density slope γ of the initial condition for full 3D (filled symbols) and 1D simulations (empty symbols).

Fig. 8 .
Fig.8.Mass ratio against ellipticity relation for the final states of spherical (red circles) and clumpy (green diamonds) initial conditions.The purple dashed line marks the Deur (2014) relation with its uncertainty (blue shaded area), while the orange dotted line marks the linear fit for the models with spherical initial conditions.