# The formation of early-type galaxies through monolithic collapse of gas
clouds in Milgromian gravity

Robin Eappen1, Pavel Kroupa1,2, Nils Wittenburg2, Moritz Haslbauer2,3 and
Benoit Famaey4
1Charles University in Prague, Faculty of Mathematics and Physics,
Astronomical Institute, V Holešovickách 2, CZ-180 00 Praha 8, Czech Republic
2Helmholtz-Institut für Strahlen- und Kernphysik (HISKP), Universität Bonn,
Nussallee 14–16, 53115 Bonn, Germany
3Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn,
Germany
4Universit´e de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR
7550, F-67000 Strasbourg, France E-mail<EMAIL_ADDRESS>mail<EMAIL_ADDRESS>
(Accepted 2022 August 4. Received 2022 August 2; in original form 2021
September 22)

###### Abstract

Studies of stellar populations in early-type galaxies (ETGs) show that the
more massive galaxies form earlier and have a shorter star formation history
(SFH). In this study, we investigate the initial conditions of ETG formation.
The study begins with the collapse of non-rotating post-Big-Bang gas clouds in
Milgromian (MOND) gravitation. These produce ETGs with star-forming timescales
(SFT) comparable to those observed in the real Universe. Comparing these
collapse models with observations, we set constraints on the initial size and
density of the post-Big-Bang gas clouds in order to form ETGs. The effective-
radius–mass relation of the model galaxies falls short of the observed
relation. Possible mechanisms for later radius expansion are discussed. Using
hydrodynamic MOND simulations this work thus for the first time shows that the
SFTs observed for ETGs may be a natural occurrence in the MOND paradigm. We
show that different feedback algorithms change the evolution of the galaxies
only to a very minor degree in MOND. The first stars have, however, formed
more rapidly in the real Universe than possible just from the here studied
gravitational collapse mechanism. Dark-matter-based cosmological structure
formation simulations disagree with the observed SFTs at more than 5 sigma
confidence.

###### keywords:

galaxies: formation – galaxies: evolution – galaxies: elliptical – galaxies:
star formation – galaxies: stellar content

††pubyear: 2022††pagerange: The formation of early-type galaxies through
monolithic collapse of gas clouds in Milgromian gravity–D

##  1 Introduction

Present-day early-type galaxies (ETGs) host old stars without much gas and
with a very low star formation activity. Studying the stellar properties of
these galaxies is very important for understanding the cosmic history of the
Universe. Glazebrook et al. (2017) notes that extreme star-formation events in
the early Universe are not rare events and that they play a significant role
in the early mass assembly of these galaxies. They report the discovery of a
quiescent galaxy with stellar mass $\approx 1.7\times 10^{11}M_{\odot}$ at z =
3.17. Martín-Navarro et al. (2018), in their stellar population analysis of
ETGs, report that the bulk of the stars form in a short timescale with star
formation lasting longer in the central regions, and they suggest that an
early monolithic formation is highly likely. The observational evidence is
reviewed in Yan et al. (2021) and Kroupa et al. (2020), who show how the
formation of super-massive black holes is a natural outcome of the formation
of early-type galaxies within the integrated galaxy-wide IMF (IGIMF, Kroupa &
Weidner 2003; Weidner & Kroupa 2006; Jeřábková et al. 2018; Yan et al. 2021)
theory.

The observations thus point towards putative progenitors of ETGs undergoing
very rapid early formation (Cowie et al., 1996; Thomas et al., 2005; Nelan et
al., 2005; Recchi et al., 2009; McDermid et al., 2015; Liu et al., 2016).
Studies based on elemental abundances of the ETGs indeed constrain them to
have formed rapidly on a timescale of $\approx$ 1 Gyr. Also, studies from
stellar population models with observed line indices show that ETGs must have
formed the bulk of their stars in a short timescale. Cowie et al. (1996),
suggested the name "downsizing" to describe that the less massive galaxies
have more extended SFHs compared to the massive ones.

The duration of the formation, i.e. the star-forming timescale (SFT) for ETGs
is expressed in Thomas et al. (2005), Recchi et al. (2009) and McDermid et al.
(2015) in terms of the downsizing time which is a function of the mass of the
galaxy. These downsizing timescales correspond to approximately 0.34
$\mathrm{Gyr}$ for a galaxy with a present-day baryonic mass of
$10^{12}M_{\rm\odot}$, suggesting that it should have formed under a
monolithic cloud collapse scenario rather than hierarchical merging which
requires longer time-spans (Wuyts et al., 2010; Ricciardelli et al., 2010). On
average, more massive galaxies have been found to have formed the bulk of
their stars on a shorter timescale and have completed their star-formation
activity at a higher redshift (Daddi et al., 2005; Wiklind et al., 2008;
Castro-Rodríguez & López-Corredoira, 2012) which points to a strong dependence
of the SFT in a galaxy on its present-day stellar-mass content.

Massive ETGs are thought to be important tracers of the cosmic history of the
stellar mass assembly and of galaxy evolution. A theoretical model that
describes ETGs should explain how and why they have these observed SFTs. The
standard model of cosmology (SMoC) has been successful in forming ETGs (see
Naab & Ostriker 2017 for a detailed review on the theory and numerical
simulations in the SMoC) but has not been able to reproduce ETGs with SFTs
(see Fig. 9) similarly short as the observed ETGs (Thomas et al., 2005; Recchi
et al., 2009; McDermid et al., 2015). With the existence of dark matter (DM)
remaining a hypothesis, an alternate theory of cosmology should be considered
(Kroupa et al., 2012; Kroupa, 2015).

By combining observational constraints on the dynamics in galaxies with
constraints from the Solar System, Milgrom (1983a) corrected the theory of
gravitation at low acceleration, which could be a consequence of the quantum
vacuum (Milgrom, 1999; Pazy, 2013; Verlinde, 2017; Smolin, 2017). In analogy
with Newtonian dynamics, a non-relativistic Milgromian gravity theory (MOND)
can be constructed by setting up a Lagrangian which, upon extremization of the
action, yields a generalised Milgromian Poisson equation. Two such Lagrangians
have been proposed, called AQUAL (Bekenstein & Milgrom, 1984) and QUMOND
(Milgrom, 2010). The latter, which we adopt hereafter, yields the following
generalized Poisson equation,

$\Delta\Phi(\vec{x})=4\pi G\rho_{\rm
b}(\vec{x})+\vec{\nabla}\cdot[\tilde{v}(\left\lvert\vec{\nabla}\phi\right\rvert/a_{0})\vec{\nabla}\phi(\vec{x})],$
(1)

or,

$\Delta\Phi(\vec{x})=4\pi G[\rho_{\rm b}(\vec{x})+\rho_{\rm ph}(\vec{x})],$
(2)

where $\rm\rho_{\rm b}(\it\vv{x})$ is the baryonic density,
$\rm\phi(\it\vv{x})$ is the Newtonian potential which fulfils the standard
Poisson equation, $\rm\Delta\phi(\it\vv{x})$ = 4$\rm\pi$G$\rm\rho_{\rm
b}(\it\vv{x})$ and Milgrom’s constant $a_{\rm 0}$ $\approx\rm 1.2\times
10^{-10}\ m\ s^{-2}\approx 3.8\ pc\ Myr^{-2}$ . The phantom dark matter (PDM)
density, $\rm\rho_{ph}(\it\vv{x})$, is not a real density distribution but a
mathematical function that arises out of the non-linearity of the Poisson
equation. $\Phi(\it\vv{x})$ is the total gravitational potential from which
the accelerations follow, $\vv{a}=-\vv{\nabla}\vv{\Phi}$, and $\tilde{v}(y)$
is a transition function characterizing the theory (see Milgrom 2008, 2010,
2014; Famaey & McGaugh 2012 and Banik & Zhao 2022 for detailed reviews on the
theory). $\tilde{v}(y)$ has the limits,

$\tilde{v}(y)\rightarrow 0\;\rm{for}\;\it y\;\gg\rm
1\;\rm{and}\;\it\tilde{v}(\it y)\;\rightarrow\it y^{\rm-1/2}\;\rm{for}\;\it
y\;\ll\rm 1.$ (3)

The above formulation deals only with linear differential equations and is
shown to emerge as a natural modification of a Palatini-type formulation of
Newtonian gravity. It is a member of a larger class of bi-potential theories
(quasi-linear formulation of MOND, QUMOND Milgrom 2010). Bekenstein & Milgrom
(1984) also noted a correspondence with some theories of quark confinement
using a different form of the function $\tilde{v}$.

MOND predicted galaxy scaling relations obeyed by galaxies (Milgrom, 1983b)
such as the Baryonic Tully Fisher Relation (BTFR) (McGaugh et al., 2000;
McGaugh, 2005, 2012; Liu et al., 2016) and the Radial Acceleration Relation
(RAR) (Sanders, 1990; McGaugh, 2004; Lelli et al., 2017).

Any realistic theory of galaxy formation has to address the question of when
and how the stars form in them. There exists a large body of literature on the
formation of ETGs in the SMoC (Naab & Ostriker, 2017) for a review. Here we
concentrate on the question whether ETGs can form in Milgromian gravitation
from collapsing post-Big-Bang gas clouds and how these models compare to the
observed ETGs. The aim of this work is to use simulations of collapsing and
non-rotating post-Big-Bang gas clouds to study the duration of their SFTs and
the resulting size–mass relation. Galaxy scale fluctuations in a MOND Universe
should grow by the application of MOND to the peculiar accelerations allowing
for large density contrasts and for large galaxies to form as early as
$z=10-25$ from the collapse of almost isolated gas clouds. Nipoti et al.
(2007) and Sanders (2008) have therefore studied disipationless collapse in a
MOND context. Here we model for the first time the collapse of massive non-
rotating gas clouds including the physics of star-formation. The results are
found to be very close to the observations. Previous work has shown the
important bulk properties of disk galaxies to follow naturally from the
collapse of initially rotating post-Big-Bang gas clouds (Wittenburg et al.,
2020).

The layout of this paper is as follows. In Section 2, the numerical hydro-
dynamical code used is discussed briefly and in Section 3, the details of the
simulations are described. The results from the work are presented in Section
4. Section 5 contains the discussion and Section 6 has conclusions along with
potential future work.

##  2 Phantom Of RAMSES (POR)

The Phantom of RAMSES (POR) code (Lüghausen et al., 2015; Nagesh et al., 2021)
is a customized version of RAMSES (Teyssier, 2002) applying the adaptive mesh
refinement (AMR) method to solve the Milgromian Poisson equation (Eq. 1). It
uses a multi-grid and a conjugate gradient solver to solve the generalised
Poisson equation 111Candlish et al. (2015) developed a similar code
independently called RAyMOND..

There are several transition functions used in the literature to accommodate
the non-linearity of the Poisson equation in MOND. POR uses,

$\tilde{v}(y)=-\frac{1}{2}+\left(\frac{1}{4}+\frac{1}{y}\right)^{\frac{1}{2}},$
(4)

where $y=|\vv{\nabla}\phi|/a_{0}$ (see Lüghausen et al. 2015).

The POR working scheme involves first solving the standard Poisson equation to
compute the Newtonian potential, and then the PDM density is calculated using
a discrete scheme (see Lüghausen et al. 2015). Then, the Poisson equation is
solved for a second time with the Newtonian density and PDM density to compute
the total gravitational potential.

As for now, POR has been successfully applied in the simulations of Antennae-
like galaxies (Renaud et al., 2016), simulations of the Sagittarius satellite
galaxy (Thomas et al., 2017), simulations of the Local Group producing the
planes of satellites (Bílek et al., 2018, 2021b; Banik et al., 2022),
simulations of streams from globular clusters (Thomas et al., 2018),
simulations of the formation of exponential disk galaxies (Wittenburg et al.,
2020), the global stability of M33 (Banik et al., 2020), the evolution of
globular-cluster systems of ultra-diffuse galaxies due to dynamical friction
(Bílek et al., 2021a) and polar-ring galaxies have also been successfully
modelled with a pre-POR code (Lüghausen et al., 2013). A detailed user-guide
can be found in Nagesh et al. (2021).

##  3 Models

We use the POR numerical code introduced in Section 2 to set-up post-Big-Bang
gas clouds in an isolated environment in-order to allow a first assessment of
how ETGs would have formed in MOND. Here, we start with initial conditions
comparable to those in Wittenburg et al. (2020), but with no initial rotation.

Stellar particles are created in POR if the gas mass–volume–density exceeds a
user-defined threshold value. The code checks at every time step whether any
cell exceeds the density threshold, and then a stellar particle is created.
Details on the conditions for the star formation threshold value and various
other code parameters are described in Lüghausen et al. (2015), Wittenburg et
al. (2020) and Nagesh et al. (2021).

Table 1: The initial conditions and results of all the models.

Model | $M_{\rm initial}$ | $R_{\rm initial}$ | $\Delta\tau_{\rm m}$ | $t_{\rm start}$ | $t_{\rm peak}$ | $r_{\rm eff}$ | $\it{SFR}_{\rm peak}$
---|---|---|---|---|---|---|---
Name | (109 $M_{\odot}$) | (kpc) | (Gyr) | (Gyr) | (Gyr) | (kpc) | ($M_{\odot}\rm yr^{-1}$)
e1 | 0.6 | 50 | 0.62 | 1.18 | 2.15 | 0.43 | 6.71E+00
e2 | 1 | 50 | 0.53 | 1.02 | 1.86 | 0.47 | 1.34E+01
e3 | 6.4 | 50 | 0.17 | 0.51 | 0.82 | 0.39 | 2.66E+02
e4 | 10 | 50 | 0.13 | 0.44 | 0.71 | 0.6 | 5.20E+02
e5 | 30 | 50 | 0.07 | 0.29 | 0.49 | 0.84 | 2.06E+03
e6 | 50 | 50 | 0.05 | 0.23 | 0.44 | 0.73 | 4.40E+03
e7 | 64 | 50 | 0.05 | 0.21 | 0.41 | 0.67 | 7.22E+03
e8 | 70 | 50 | 0.04 | 0.2 | 0.39 | 0.77 | 7.94E+03
e9 | 100 | 50 | 0.04 | 0.17 | 0.35 | 0.83 | 1.45E+04
e10 | 0.6 | 100 | 0.91 | 2.46 | 3.5 | 0.4 | 4.88E+00
e11 | 1 | 100 | 0.78 | 2.11 | 3.95 | 0.53 | 1.03E+01
e12 | 5 | 100 | 0.37 | 1.19 | 1.73 | 0.47 | 9.67E+01
e13 | 10 | 100 | 0.26 | 0.95 | 1.41 | 0.47 | 2.77E+02
e14 | 20 | 100 | 0.19 | 0.73 | 1.21 | 1.41 | 6.05E+02
e15 | 30 | 100 | 0.16 | 0.63 | 1.07 | 0.56 | 1.12E+03
e16 | 50 | 100 | 0.13 | 0.52 | 0.92 | 0.68 | 2.52E+03
e17 | 64 | 100 | 0.11 | 0.48 | 0.87 | 0.75 | 3.97E+03
e18 | 70 | 100 | 0.11 | 0.46 | 0.85 | 0.77 | 4.40E+03
e19 | 100 | 100 | 0.09 | 0.4 | 0.76 | 0.99 | 7.78E+03
e20 | 0.6 | 200 | 1.62 | 5.12 | 6.4 | 0.4 | 3.15E+00
e21 | 1 | 200 | 1.38 | 4.37 | 7.57 | 0.56 | 6.38E+00
e22 | 10 | 200 | 0.53 | 1.95 | 2.98 | 0.55 | 1.27E+02
e23 | 30 | 200 | 0.34 | 1.31 | 2.18 | 0.68 | 6.25E+02
e24 | 50 | 200 | 0.28 | 1.11 | 1.90 | 1.21 | 1.32E+03
e25 | 70 | 200 | 0.25 | 0.99 | 1.76 | 0.93 | 2.18E+03
e26 | 100 | 200 | 0.22 | 0.87 | 1.59 | 0.93 | 3.62E+03
e27 | 0.6 | 300 | 1.65 | 7.91 | 9.29 | 0.31 | 2.54E+00
e28 | 5 | 300 | 0.79 | 3.81 | 5.09 | 0.46 | 4.51E+01
e29 | 10 | 300 | 0.65 | 2.99 | 4.48 | 1.03 | 9.41E+01
e30 | 30 | 300 | 0.48 | 2.02 | 3.31 | 0.76 | 4.48E+02
e31 | 50 | 300 | 0.42 | 1.68 | 2.9 | 0.73 | 8.96E+02
e32 | 70 | 300 | 0.39 | 1.49 | 2.66 | 0.83 | 1.38E+03
e33 | 100 | 300 | 0.35 | 1.32 | 2.42 | 1.36 | 2.34E+03
e34 | 6.4 | 500 | 1.41 | 5.97 | 8.44 | 0.41 | 3.29E+01
e34c | 6.4 | 500 | 1.43 | 5.96 | 8.41 | 0.39 | 3.25E+01
e35 | 10 | 500 | 1.18 | 5.09 | 7.45 | 0.47 | 6.41E+01
e36 | 30 | 500 | 0.75 | 3.43 | 5.58 | 0.6 | 2.85E+02
e37 | 50 | 500 | 0.62 | 2.85 | 4.88 | 0.84 | 5.95E+02
e38 | 64 | 500 | 0.55 | 2.58 | 4.48 | 0.95 | 8.60E+02
e39 | 70 | 500 | 0.54 | 2.48 | 4.4 | 1.04 | 9.39E+02
e39c | 70 | 500 | 0.53 | 2.47 | 4.44 | 1.09 | 9.48E+02
e40 | 1000 | 500 | 0.24 | 1.03 | 1.9 | 1.7 | 1.48E+04

Note: $M_{\rm initial}$ is the initial mass of the post-Big-Bang gas cloud,
$R_{\rm initial}$ is the initial radius of the post-Big-Bang gas cloud,
$\Delta\tau_{\rm m}$ is the resulting SFT, $t_{\rm peak}$ is the time after
start of the simulation when the first stellar particle is formed, $t_{\rm
peak}$ is the resulting time when the peak of the SFH is observed in the
simulation, $r_{\rm eff}$ is the resulting projected effective-radius and
$\it{SFR}_{\rm peak}$ is the resulting SFR at the peak of the SFH. All the
models here have a maximum resolution of 0.24 kpc except for models e34c and
e39c which have a maximum resolution of 0.06 kpc and 0.12 kpc, respectively.
Models e34c and e39c are computations with complex feedback (see Section 4.3).

In total 42 different model galaxies were calculated for this study to
understand how different initial masses, $M_{\rm initial}$, and initial radii,
$R_{\rm initial}$, of the gas cloud can affect the SFT (Table 1). Only the
simple cooling/heating feedback algorithm is used in the simulation of all the
model galaxies except models e34c and e39c, where a more complex feedback
algorithm is used (see Section 4.3). The computations are run for 10 Gyr. The
starting temperature of the gas is $\mathrm{10^{4}K}$ and the size of the
simulation box is 1000 kpc. The minimum and maximum refinement levels for the
models 222Models e34c and e39c have maximum refinement levels of 14 and 13,
respectively, which sets the maximum spatial resolution to 0.06 kpc and 0.12
kpc. are 7 and 12 respectively which sets a limit to the minimum spatial
resolution of 7.81 kpc and maximum spatial resolution of 0.24 kpc,
respectively, for each simulation (Teyssier, 2002; Lüghausen et al., 2015).

Figure 1: SFH of model e39 (Table 1). The red dots are the SFR at each $\delta
t$ and the dashed grey line is the calculated best-fit (Eq. 5). The x-axis is
the time in the computational domain where the computation and thus the
collapse starts at 0 and the calculation ends at 10 Gyr. The fit is done only
to obtain the FWHM of the SFH and the discrepancy at large $t$ is not
important for the purpose of this study.

In all computations, the star formation rate (SFR) is calculated by separating
all stellar particles in bins of $\delta t$ = 10 Myr according to their age,
summing up the stellar mass in every time bin and dividing this by the length
of $\delta t$. The SFR increases sharply at the time of the collapse of the
sphere until the maximum is reached and then it decreases approximately
exponentially (Fig. 1). The SFH of a galaxy is constructed by the consecutive
SFRs of each $\delta t$ time step. The SFH can be represented by a Lorentz
function (dashed grey line in Fig. 1),

$\it{SFR}/{M_{\odot}\mathrm{yr^{-1}}}=\frac{\it{SFR}_{\mathrm{peak}}/M_{\odot}\mathrm{yr^{-1}}}{\mathrm{\pi}}+\\\
\left(\frac{0.5\cdot\Delta\tau_{\mathrm{m}}/\mathrm{Gyr}}{(t/\mathrm{Gyr}-t_{\mathrm{peak}}/\mathrm{Gyr})^{2}+(0.5\cdot\Delta\tau_{\mathrm{m}}/\mathrm{Gyr})^{2}}\right),$
(5)

where $t_{\rm peak}$ is the time in the simulation where the peak of the SFH
is observed and $\Delta\tau_{\rm m}$ is the SFT for the model. The full-width
at half maximum (FWHM) of the SFH gives the duration for which the bulk of the
star formation takes place, so we use the FWHM to define the SFT (i.e.
$\Delta\tau_{\rm m}$) for the models in this work.

The collapse for models with the same $M_{\rm initial}$ and varying $R_{\rm
initial}$ shows that $\Delta\tau_{\rm m}$ increases with increasing $R_{\rm
initial}$. The sum of all formed stellar particles at the end of each
simulation is the final stellar mass 333We do not list the $M_{\rm final}$
since the final mass of the model galaxy is equal to the initial mass of the
gas cloud., $M_{\rm final}$, of each model, independently of the chosen
feedback algorithm. The models here do not loose mass through outflows.

We also calculate the time taken for the first stellar particle to form,
$t_{\rm start}$, and compare it to the MOND free-fall time,

$t_{\mathrm{ff,mond}}=\left(\frac{\pi}{2}\right)^{0.5}\cdot\frac{R_{\mathrm{initial}}}{(G\cdot
M_{\mathrm{initial}}\cdot a_{\rm 0})^{0.25}},$ (6)

(eq. 27 in Zonoozi et al. 2021) in Fig. 11. The free-fall time can be used as
an approximation for the time the gas needs to collapse into stars. The models
have $t_{\rm start}$ comparable to the MOND free-fall time (Fig. 11).

##  4 Results

We compare the SFTs, ages, and effective-radii of our simulated galaxies with
the observed galaxy properties in this section. The average ages estimated by
Thomas et al. (2005) from absorption line indices of 124 ETGs, $t_{\rm
thomas}$, for both high-density and low-density environments (quantities in
brackets are the values derived for the high-density environment) can be
written,

$t_{\mathrm{thomas}}=1.33(2.67)\cdot(M_{\mathrm{final}}/M_{\mathrm{\odot}})^{0.07(0.05)},$
(7)

(eq. 3 in Thomas et al. 2005) and the relation deduced for the SFTs is,

$\Delta\tau_{\rm
thomas}/\mathrm{Gyr}=3.67-0.37\cdot\mathrm{log_{10}}(M_{\mathrm{final}}/M_{\mathrm{\odot}})$
(8)

(eq. 5 in Thomas et al. 2005), being the same for high- and low-density
environments.

###  4.1 Downsizing SFH

The SFTs are found to follow a similar downsizing behaviour as documented in
Thomas et al. (2005, 2010), where, for a given $R_{\rm initial}$, the galaxies
with larger masses (higher density) form at an earlier epoch compared to
galaxies with lower masses (lower density). Fig. 2 shows that the model
galaxies with high initial cloud density, $\rho_{\rm initial}$, form earlier
and quicker (shorter $\Delta\tau_{\rm m}$) than model galaxies with low
$\rho_{\rm initial}$. We use a power law function to obtain the scaling
relation in Fig. 2,

$\Delta\tau_{\rm m}/\mathrm{Gyr}=3.03\cdot(\rho_{\rm
initial}/M_{\odot}\mathrm{kpc^{-3}})^{-0.30}.$ (9)

Figure 2: SFT–density relation for all models listed in Table 1. The coloured
markers are model galaxies and the dashed grey line is the calculated best-fit
(Eq. 9). Figure 3: SFT–mass relation. The coloured markers are model galaxies
listed in Table 1 with different $R_{\rm initial}$ as indicated by the legend.
The observationally constrained relations are from Thomas et al. (2005)
(dashed green line), Recchi et al. (2009) (solid green line) and Yan et al.
(2021) (which is a combination of eq. 3 of Thomas et al. 2005 and eq. 3 of
McDermid et al. 2015, dotted green line). The red left-pointing triangles are
the models with complex feedback (see Section 4.3).

A downsizing SFH can be well explained using Fig. 3 where the SFT is plotted
against $M_{\rm final}$ for each model galaxy along with the observations from
Thomas et al. (2005), Recchi et al. (2009) and McDermid et al. (2015).

We use a power law function to fit the models and obtain the scaling relation
for the models in Fig. 3,

$\Delta\tau_{\rm m}/\mathrm{Gyr}=A\cdot(M_{\rm final}/M_{\odot})^{B},$ (10)

where $A$ and $B$ are fit parameters for each $R_{\rm initial}$ given in Table
4. From Fig. 3, it can be seen that models with $R_{\rm initial}$ = 500 kpc
lead to similar $\Delta\tau_{\rm m}$ values as the observational results from
Thomas et al. (2005) and models with $R_{\rm initial}$ = 200 kpc and $R_{\rm
initial}$ = 300 kpc are similar to the observational data from Recchi et al.
(2009) but steeper. The models with $R_{\rm initial}$ = 50 kpc and $R_{\rm
initial}$ = 100 kpc do not agree with any of the observed relations. We note
that our models do not agree with the new downsizing relation obtained by Yan
et al. (2021) and this is discussed in Section 5.

Figure 4: A plot showing the ratio of SFTs of the model galaxies and observed
galaxies ($\Delta\tau_{\rm m}/\Delta\it\tau_{\rm{thomas}}$) vs final stellar
mass of the galaxy. The coloured markers are model galaxies as listed in Table
1 with different $R_{\rm initial}$ as indicated by the legend. The dashed
black line is where $\Delta\tau_{\rm m}$ = $\Delta\tau_{\rm{thomas}}$.

Fig. 4 shows that only models with $R_{\rm initial}$ = 500 kpc
($\Delta\tau_{\rm m}/\Delta\it\tau_{\rm{thomas}}$ $\approx$ 1) are comparable
to the results from Thomas et al. (2005). From here, we only consider models
with $R_{\rm initial}$ = 500 kpc for discussions. Using the basic assumption
that the standard age of the Universe is 13.8 Gyr and the age-dating deduced
by Thomas et al. (2005) for ETGs is valid, we translate the SFH of our models
to a time related to the start of the Big-Bang. Since our models have SFTs
comparable to Thomas et al. (2005), we assume that the peaks of the SFHs of
our models should coincide with the average ages computed by Thomas et al.
(2005). This allows us to place our models to a time relative to the Big-Bang
(Eq. 7). We shift the SFHs of our model galaxies to a time corresponding to
the average ages as deduced by Thomas et al. (2005). To achieve this, we first
transform the observed average formation times (Eq. 7) to time since the Big-
Bang,

$t_{\rm{av,th}}/{\rm Gyr}=13.8-t_{\rm{thomas}}/{\rm Gyr},$ (11)

and the time the first stellar particles would have formed according to the
models presented here becomes,

$t_{\rm{first}}/{\rm Gyr}=t_{\rm{av,th}}/{\rm Gyr}-(t_{\rm{peak}}/{\rm
Gyr}-t_{\rm{start}}/{\rm Gyr}),$ (12)

where $t_{\rm{start}}$ is the time when the first stellar particle is formed
in the simulation and $t_{\rm{first}}$ is the time when the first star would
have formed in the ETG relative to the Big-Bang that is consistent with the
Thomas et al. (2005) age-mass relation. Fig. 5 shows the average age and when
the first star forms in the models for both low- and high-density
environments. Fig. 6 is a representation of the SFH of the model galaxies
which follow the downsizing timescales as deduced by Thomas et al. (2005).
That is, a massive ($\approx\mathrm{10^{12}}\ M_{\odot}$) post-Big-Bang gas
cloud would have started to form an ETG around 1.5 Gyr after the Big-Bang in a
high-density environment and 3.6 Gyr after the Big-Bang in a low-density
environment. The observed quasar activity only $\approx$ 0.70 Gyr after the
Big-Bang suggests that there may have been significantly more over-dense
regions (Kroupa et al., 2020). The transparent SFHs in Fig. 6 are on the time
axis where $t$ = 0 is the start of the simulation. Thus, the real Universe
appear to have accelerated the start of the collapse even in the low-density
regions. Rather than the real Universe having accelerated the onset of star
formation, another interpretation of Fig. 6 is possible though. Fig. 6 shows
that the low-mass models agree with the observed galaxies in terms of the time
when the first stars formed in low-density environments, while the massive
models agree with the observed timings in high-density regions. This suggests
that low-density regions form preferentially low-mass galaxies, while high-
density regions form preferentially massive galaxies. This will be testable
with cosmological structure formation simulations in MOND.

Figure 5: $t_{\rm av,th}$ and $t_{\rm first}$ vs. $M_{\rm final}$ for models
with $R_{\rm initial}$ = 500 kpc in the low-density (_left_) and high-density
(_right_) environment (Eq. 11 $\&$ 12). The times are computed since the
nominal Big-Bang. The unfilled upward pointing purple triangles in the second
panel are relative to the start of the computation at $t$ = 0.

Figure 6: A plot showing the SFH of model e40, e39, e38, e37, e36 and e35
(Table 1) in the low-density (_left_) and high-density (_right_) environment.
The black line is the time when the first star is formed in the models. Each
model has been shifted in time according to the Thomas et al. (2005) age-
dating (Eqs. 11 & 12). The nearly transparent SFHs are relative to the start
of the computation at $t$ = 0.

###  4.2 Effective-radii of the model galaxies

Figure 7: The effective-radius–mass relation. The green triangles are the real
ETGs from Dabringhausen & Fellhauer (2016), purple triangles are the model
galaxies e34 to e40 with $R_{\rm initial}$ = 500 kpc listed in Table 1 and the
blue M’s are the merged models listed in Table 2. The black dashed line is the
best fit for the real galaxies (Eq. 14) and the grey dashed line is the best
fit for the model galaxies (Eq. 15). The solid red line is the MOND radius
(Eq. 13) and the red left-pointing triangles are the models with complex
feedback (see Section 4.3).

The effective-radius is defined as the projected radius within which half of
the final stellar mass of the model galaxy is found, assuming that the mass
density follows the luminosity density, this can be compared to the projected
half-light radius from observations. The effective-radius is calculated by
projecting the model galaxy into the XY plane (Eappen et al. 2022c will
further discuss the effective-radius in other projections as well). The shape
of the model galaxy is found to be disk-like similar to Wittenburg et al.
(2020). The collapse occurs along the $Z$-direction such that the mid-plane of
the formed galaxy is the XY plane. The effective-radius, $r_{\rm eff}$, is
plotted against $M_{\rm final}$ of the model galaxies and compared to the real
galaxies from Dabringhausen & Fellhauer (2016) in Fig. 7.

The MOND radius, $r_{\rm M}$, at which the Newtonian radial acceleration
equals $a_{\rm 0}$ (Milgrom, 2014) is,

$r_{\rm{M}}=\left(\frac{\mathrm{G}\cdot M_{\rm final}}{a_{0}}\right)^{0.5},\\\
$ (13)

and is plotted as the solid red line in Fig. 7. Sanders (2008) notes that the
effective radius of the dissipationless models in his work are comparable to
$r_{\rm M}$ that is, galaxy mass ($\approx 10^{11}\ M_{\odot}$) objects
naturally collapse to a radius of about 10 kpc. The hydro-dynamical
dissipative simulations of the models in this work form model galaxies similar
in shape, size and kinematical properties to the compact relic galaxies like
NGC 1277 (Trujillo et al., 2014) and this will be discussed in a companion
paper (Eappen et al. 2022b in preparation).

The best-fit relation for the Dabringhausen & Fellhauer (2016) galaxies in
Fig. 7 is given by,

$r_{\rm{D}}/{\rm pc}=0.46\cdot(M_{\rm\star}/M_{\odot})^{0.35},$ (14)

where $r_{\rm D}$ is the observed effective-radius calculated by Dabringhausen
& Fellhauer (2016) and $M_{\star}$ is the total stellar mass of the galaxy.
The best-fit for the model galaxies in Fig. 7 is calculated using a power law
function,

$r_{\rm{eff}}/{\rm pc}=0.12\cdot(M_{\rm final}/M_{\odot})^{0.36}.$ (15)

It can be seen in Fig. 7 that the effective-radius – mass relation (Eq. 15) of
the model galaxies has a power-law index comparable to the effective-radius –
mass relation (Eq. 14) of the observed ETGs from Dabringhausen & Fellhauer
(2016) but that it falls below the observed relation by a factor of 4.

It is possible that the monolithic collapse which forms individual compact
ETGs occur in the inner regions of local density enhancements that form
neighbouring gas clouds which collapse monolithically, forming binary or even
more ETGs (similar to collapse of stars in a very dense configuration e.g.
Joseph & Wright 1985, Wright et al. 1988 and Stolte et al. 2014). If these are
on radial orbits they would merge, leaving a merger remnant which is likely to
have a larger effective-radius (For further empirical evidence that early-type
galaxies formed first at their innermost regions see section 2.2 in Kroupa et
al. 2020). Trujillo et al. (2007) study an effective size evolution of compact
galaxies where they find that the massive ETGs were formed in the early
Universe, and have subsequently evolved passively until today. They find an
effective size evolutionary mechanism would be able to evolve their compact
galaxies to the observed local relation with just two major mergers. To
investigate this we study the merger scenario where two model galaxies listed
in Table 1 are placed 100 kpc apart in a box of 1 Mpc. The models only consist
of stellar particles. The spin angular momentum of each model galaxies is
retained as they are placed in the new box and they merge under MOND gravity
to form a new merged model galaxy in 1.5 Gyr. Table 2 lists the initial
conditions and results of the merged models. Five different merger studies
were done placing the models at different axes and at different orientation.
M1, M2, M3 and M4 are merged models formed from the merger of the same model
galaxies and M5 is formed by the merger of model galaxies of different masses
(e39 and e36) that merge at different angles with respect to the XY plane
(where parent models e36 and e39 are placed tilted at an angle of 45 degrees
with respect to the XY plane). The effective-radius–mass relation of these
merged models is plotted in Fig. 7. It can be seen that the merged models lie
around the same dashed black line as the real ETGs in Fig. 7. A detailed study
of the morphology and other structural properties of these merged galaxies
will be discussed in a follow up paper (Eappen et al. 2022c).

Another possible mechanism for the model galaxies to reach large radii is
through stellar evolution mass loss which is directly linked to the galaxy-
wide IMF assumed. It has been shown that the stellar remnant population is
comparable to the mass of living stars in the IGIMF theory but using a
canonical invariant IMF would decrease the remnant mass to be much lower than
the stellar mass (Yan et al., 2019). Yan et al. (2021) finds that the ratio of
dynamical mass and stellar mass, $M_{\rm dyn}/M_{\star}$, is higher for more
massive galaxies because they have a more top-heavy galaxy-wide IMF and more
stellar remnants. Given the high SFRs (> $10^{2}M_{\rm\odot}/\rm yr$) the
model galaxies reach (Table 1), it is expected that the IGIMF is top-heavy
(Yan et al., 2017; Jeřábková et al., 2018; Yan et al., 2021) leading to
significant mass loss from the evolving stars. This mass lost would be in the
form of gas and if it is heated sufficiently enough it would expand and lead
to the expansion of the galaxy (Suzuki et al., 2022).

Table 2: Initial conditions of the merger models

Model | Parent | $r_{\rm eff}$ | $M_{\rm final}$
---|---|---|---
Name | models | (kpc) | ($M_{\rm\odot}$)
M1 | e35 + e35 | 2.6 | 0.17E11
M2 | e36 + e36 | 2.9 | 0.38E11
M3 | e37 + e37 | 3.2 | 0.9E11
M4 | e39 + e39 | 3.0 | 1.25E11
M5 | Re36 + Re39 | 2.4 | 0.75E11

Note: Column 1: name, Column 2: The models from Table 1 that merge, Column 3
($r_{\rm eff}$): the effective-radius of the merged model and Column 4
($M_{\rm final}$): the mass of the merged model. Models Re36 and Re39 are
models e36 and e39 but rotated (see text for details).

###  4.3 Effects of complex feedback mechanism

Figure 8: A plot comparing the SFH of models with simple and complex feedback.
Models e34c and e39c are models with complex feedback.

The above results (except for models e34c and e39c) have been obtained using
the simple feedback algorithm (simple cooling/heating and no additional more
complex baryonic physics, such as supernovae). Models e34c and e39c were
computed with additional feedback process to investigate how complex baryonic
physics affects the results obtained.

The simple star-formation algorithm only considers the heating/cooling modules
in PoR. Cooling and heating of the gas are computed by using tables that are
included in the code, which describe the cooling and heating rates due to
several physical processes. PoR uses Courty & Alimi (2004) look-up tables for
the different cooling/heating processes. The complex feedback algorithm
includes more initial inputs such as supernovae, radiative transfer, sink
particles and a higher refinement level which increases the resolution of the
simulation. The supernova (and in fact all feedback scheme) scheme is
implemented as in RAMSES where the kinetic part of the supernova energy is
injected as a spherical blast wave with the size of galactic superbubbles of
radius 150 pc. The radiative transfer option allows to compute the radiative
transfer between stellar particles with the addition of photon fluxes which
uses a first-order Godunov solver (Rosdahl et al., 2013). A sink particle
scheme is used to stabilize the simulation by artificially stopping the
collapsing gas cloud if a certain density is exceeded and the gas is condensed
into a point mass (sink particle, see Teyssier 2002; Lüghausen et al. 2015;
Wittenburg et al. 2020; Nagesh et al. 2021 for a detailed description on how
RAMSES and PoR handles feedback).

Table 3: Effects of feedback

| Simple Feedback | Complex Feedback |
---|---|---|---
| e34 | e34c | PC
$r_{\rm eff}$ | 0.41 | 0.39 | \- 4.87%
$\Delta\tau_{\rm m}$ | 1.41 | 1.43 | \+ 1.41%
| e39 | e39c | PC
$r_{\rm eff}$ | 1.04 | 1.09 | \+ 4.80%
$\Delta\tau_{\rm m}$ | 0.54 | 0.53 | \- 1.85%

Note: PC shows the percentage change (Eq. 16) in the results obtained from
models with same $R_{\rm initial}$ and $M_{\rm initial}$ but with simple (e34
$\&$ e39) and complex (e34c $\&$ e39c) feedback.

The model e39c has the same $R_{\rm initial}$ and $M_{\rm initial}$ as model
e39 but is computed with complex baryonic feedback of supernovae, sink
particles, and radiative transfer with a higher resolution (maximum refinement
level set to 13 instead of 12). Similarly model e34c has the same initial
conditions as model e34 but is computed with a higher resolution (maximum
refinement level set to 14 instead of 12). Fig. 8 shows the SFH of the models
with simple and complex feedback. Table 3 shows the effect of feedback on the
models which have the same $R_{\rm initial}$ and $M_{\rm initial}$. The
percentage change (PC) is calculated using Eq. 16 in Appendix D and the
results obtained (Table 3) show that the percentage change of the results from
the complex feedback is within 5$\%$ of the results obtained from simple
feedback. The complex feedback is significantly more CPU intensive but is
found to affect the results negligibly (see Fig. 3, 7 and 8 and Table 3), as
already shown in the disk-galaxy formation simulations by Wittenburg et al.
(2020). This demonstrates again that differences in the feedback processes
play a minor role in the evolution of galaxies and that the physics of galaxy
evolution are pre-dominantly defined by the Milgromian law of gravitation
(Kroupa, 2015).

##  5 Discussion

The monolithic post-Big-Bang collapse of non-rotating gas clouds in MOND
results in a SFT–mass relation (Eq. 10) which is comparable to the observed
relation for ETGs (Thomas et al., 2005, 2010). Collapsing post-Big-Bang gas
clouds in Milgromian gravitation thus behave as the observed ETGs in terms of
the collapse timing (i.e more massive clouds collapse faster and form their
stellar particles more quickly as shown in Fig. 2). The grid of models with
different $R_{\rm initial}$ and $M_{\rm initial}$ allows us to constrain
$R_{\rm initial}\approx$ 500 kpc as best fitting the observed SFT–mass
relation. This suggests that typically a region spanning almost a Mpc across
collapsed to form ETGs if sufficiently massive. We note that field elliptical
galaxies are slightly younger (consistent with the results found by Thomas et
al. 2005) but otherwise do not differ significantly from elliptical galaxies
in clusters (Saracco et al., 2017). This suggests that also in the field the
bulk stellar population of elliptical galaxies formed through monolithic
collapse on the downsizing time scale.

How do these results compare to the expectations from the standard
hierarchical structure formation theory as quantified through the Illustris
TNG (Pillepich et al., 2018a; Nelson et al., 2018, 2019) and Eagle (Schaye et
al., 2015; McAlpine et al., 2016) projects? It is noted that Fontanot & Monaco
(2010) write: "The so-called downsizing trend in galaxy formation (e.g. more
massive galaxies forming on a shorter time-scale and at higher redshift with
respect to lower mass counterparts) is not fully recovered." On the other
hand, based on their analysis of the Millenium Simulation, De Lucia et al.
(2006) claim "These findings are consistent with recent observational results
that suggest ‘down-sizing’ or ‘anti-hierarchical’ behaviour for the star
formation history of the elliptical galaxy population, despite the fact that
our model includes all the standard elements of hierarchical galaxy formation
and is implemented on the standard, $\Lambda$CDM cosmogony." Given these
discrepant results based on the same theory and observational data, we
evaluate the formation time scales of early type galaxies that form in the
Illustris TNG and EAGLE projects (see Figs. 9 and 10), noting that both of
these rest on very different simulation codes (adaptive mesh refinement versus
SPH, respectively) and also very different baryonic physics algorithms. The
details are provided in Appendix B and the result is, in summary, that the
formation time scales of early type galaxies formed in both of these
simulation projects do not concur with the observed population.

The MOND radius (Eq. 13) proposed by theory is larger but follows the same
trend as the observed $r_{\rm eff}$. This similarity of the MOND radius and
the effective-radius was first noted by Sanders (2008) where he concludes that
a dissipationless initially expanding gas cloud collapse would naturally form
galaxies with effective-radii that are comparable to their MOND radii. The
models in our study do not consider an expanding gas cloud collapse and fall
short of this expected effective-radius–mass relation (see Fig. 7) but this is
possibly also due to a non-variant galaxy-wide IMF which makes the post-Big-
Bang clouds to collapse too deeply. If, on the other hand, the galaxy-wide IMF
varies such that it becomes top-heavy at a high SFR (eg. Jeřábková et al.
2018; Yan et al. 2021) then the collapse is likely less deep and $r_{\rm eff}$
will end up larger. We do not address the chemical evolution of the models
here as the chemical enrichment has not yet been incorporated into PoR and
would also require a systematically varying galaxy-wide IMF (the current
simulationss effectively assume the canonical IMF for the formation of stellar
particles). The SFHs of the models presented here are however comparable, by
virtue of them matching the formation timescales, with the constraints from
abundance and stellar population observations as analysed by Yan et al.
(2021). It is to be expected that the bulk chemical enrichment of the models
computed here would follow the closed box models of Yan et al. (2021). The
problem of spatially-resolved chemical enrichment of the forming model ETGs
will be studied in the future in connection with the systematically varying
galaxy-wide IMF.

The model galaxies in this work do not agree with the SFT–mass relation found
in Yan et al. (2021) (dotted green curve in Fig. 3). This is possibly due to
the non-variant IMF that is assumed in the simulation as mentioned above. The
other method through which the model galaxies could reach the Yan et al.
(2021) SFT–mass relation is if the gas cloud collapse calculations are
considered in an expanding cosmological volume which will be investigated in
the future.

In Fig. 7 we show that mergers of the model galaxies form models that lie on
the observationally deduced effective-radius–mass relation. These kinds of
mergers could happen in the initially densest parts of the Universe, for e.g.
in central regions of galaxy clusters, where the galaxies would have initially
formed through monolithic collapse of post-Big-Bang gas clouds and then
slightly later merge with another galaxy of comparable mass in a timescale
somewhat larger than the formation timescale of these galaxies to reach the
observed effective-radius–mass relation. The morphology and structural
properties of the model galaxies and the mergers using different projections
will be discussed in a follow up paper (Eappen et al. 2022c in preparation)
where we find that the shape of these merged models is triaxial and very
similar to the observed ETGs. Thus, the formation timescales and other
observed chemical properties of ETGs could have been embedded in them before
they undergo one or two major mergers such that only the structural properties
and the morphology of these ETGs changed through the mergers. Trujillo et al.
(2007) reports a similar result in their study of size evolution of massive
galaxies at different redshifts. They find that at z $\approx$ 1.5 massive
spheroid like objects were a factor of 4 smaller than we see today and suggest
that just two major mergers of equal mass would evolve the size of the galaxy
to what we observe today. Similarly, Pipino & Matteucci (2008) in their study
of ETGs through chemical evolution models find that a series of dry mergers of
galaxies is not the way to recover downsizing and suggest 1–3 major-dry
mergers to be in agreement with the observations. Our results show that it is
possible to evolve the size and mass of the galaxy without changing the
intrinsic properties of the galaxy such as the SFT without a large number of
mergers.

##  6 Conclusion

In this work, we discuss the formation of ETGs in comparison with
observational results. The isolated monolithic collapse of post-Big-Bang
clouds of different initial radii at fixed mass is computed in Milgromian
gravitation and is found to produce model galaxies with SFTs similar to the
observed ETGs.

The main results are:

  * •

This work for the first time shows that the SFTs observed for ETGs is a
natural occurrence in the MOND paradigm if ETGs form from monolithically
collapsing non-rotating post-Big-Bang gas clouds. We were able to reproduce
the SFTs as noted by observations (Fig. 3). The models with $R_{\rm initial}$
= 500 kpc are similar to the relation of Thomas et al. (2005, 2010) and the
models with $R_{\rm initial}$ = 200 kpc fit fairly closely to the relation of
Recchi et al. (2009). This work thereby gives a rough estimate of the initial
conditions of the post-Big-Bang clouds in order to produce ETGs. The results
indicate that the observed Universe started to form ETGs earlier than the here
studied collapse models, although the SFTs of the models match those observed
for gas clouds with initial radii in the range 200 to 500 kpc. Higher
resolution simulations will be needed to address this issue.

  * •

The effective-radius–mass relation of the model galaxies falls short of the
expected value but the models with only one major merger (Fig. 7) are found to
be consistent with observations (Dabringhausen & Fellhauer, 2016).

  * •

We find that different feedback processes affect the results negligibly (as
already shown by Wittenburg et al. 2020 in their formation of disk-galaxies)
suggesting that the physics of galaxy formation is pre-dominantly defined by
the Milgromian law of gravitation (Kroupa, 2015).

The future aspects of this project will analyze the structural properties and
shapes of both the model galaxy and merged galaxies using different
projections to calculate the sizes (Eappen et al. 2022c, in preparation). Then
the next step would be to investigate how a possible inclusion of a
systematically varying galaxy-wide IMF (Yan et al., 2017) would lead to the
expansion of the models also taking into account the cosmological expansion of
space.

## Acknowledgements

We would like to thank Zhiqiang Yan for providing insightful suggestions. The
numerical simulations were performed using the Tiger–Cluster of the
Astronomical Institute of Charles University, Prague. R.E. is supported by the
Grant Agency of Charles University under grant No. 234122. B.F. acknowledges
funding from the Agence Nationale de la Recherche (ANR projects
ANR-18-CE31-0006 and ANR-19-CE31-0017), and from the European Research Council
(ERC) under the European Union’s Horizon 2020 Framework programme (grant
agreement number 834148). We thank the DAAD Eastern European grant at Bonn
University for supporting the research visits.

## Data Availability

All data used here have been generated as described using the publicly
available POR code and the cited literature.

## References

  * Banik & Zhao (2022) Banik I., Zhao H., 2022, Symmetry, 14, 1331
  * Banik et al. (2020) Banik I., Thies I., Famaey B., Candlish G., Kroupa P., Ibata R., 2020, ApJ, 905, 135
  * Banik et al. (2022) Banik I., Thies I., Truelove R., Candlish G., Famaey B., Pawlowski M. S., Ibata R., Kroupa P., 2022, MNRAS, 513, 129
  * Bekenstein & Milgrom (1984) Bekenstein J., Milgrom M., 1984, ApJ, 286, 7
  * Bhattacharyya & Johnson (1977) Bhattacharyya G. K., Johnson R. A., 1977, John Wiley & Sons, New York
  * Bílek et al. (2018) Bílek M., Thies I., Kroupa P., Famaey B., 2018, A&A, 614, A59
  * Bílek et al. (2021a) Bílek M., Zhao H., Famaey B., Müller O., Kroupa P., Ibata R., 2021a, A&A, p. 537
  * Bílek et al. (2021b) Bílek M., Thies I., Kroupa P., Famaey B., 2021b, Galaxies, 9, 100
  * Candlish et al. (2015) Candlish G. N., Smith R., Fellhauer M., 2015, MNRAS, 446, 1060
  * Castro-Rodríguez & López-Corredoira (2012) Castro-Rodríguez N., López-Corredoira M., 2012, A&A, 537, A31
  * Courty & Alimi (2004) Courty S., Alimi J. M., 2004, A&A, 416, 875
  * Cowie et al. (1996) Cowie L. L., Songaila A., Hu E. M., Cohen J. G., 1996, AJ, 112, 839
  * Crain et al. (2015) Crain R. A., et al., 2015, MNRAS, 450, 1937
  * Dabringhausen & Fellhauer (2016) Dabringhausen J., Fellhauer M., 2016, MNRAS, 460, 4492
  * Daddi et al. (2005) Daddi E., et al., 2005, ApJ, 626, 680
  * De Lucia et al. (2006) De Lucia G., Springel V., White S. D. M., Croton D., Kauffmann G., 2006, MNRAS, 366, 499
  * Famaey & McGaugh (2012) Famaey B., McGaugh S. S., 2012, Living Reviews in Relativity, 15, 10
  * Fontanot & Monaco (2010) Fontanot F., Monaco P., 2010, MNRAS, 405, 705
  * Glazebrook et al. (2017) Glazebrook K., et al., 2017, Nature, 544, 71
  * Jeřábková et al. (2018) Jeřábková T., Hasani Zonoozi A., Kroupa P., Beccari G., Yan Z., Vazdekis A., Zhang Z. Y., 2018, A&A, 620, A39
  * Joseph & Wright (1985) Joseph R. D., Wright G. S., 1985, MNRAS, 214, 87
  * Kroupa (2015) Kroupa P., 2015, Canadian Journal of Physics, 93, 169
  * Kroupa & Weidner (2003) Kroupa P., Weidner C., 2003, ApJ, 598, 1076
  * Kroupa et al. (2012) Kroupa P., Pawlowski M., Milgrom M., 2012, International Journal of Modern Physics D, 21, 1230003
  * Kroupa et al. (2020) Kroupa P., Subr L., Jerabkova T., Wang L., 2020, MNRAS, 498, 5652
  * Lelli et al. (2017) Lelli F., McGaugh S. S., Schombert J. M., Pawlowski M. S., 2017, ApJ, 836, 152
  * Liu et al. (2016) Liu Y., et al., 2016, ApJ, 818, 179
  * Lüghausen et al. (2013) Lüghausen F., Famaey B., Kroupa P., Angus G., Combes F., Gentile G., Tiret O., Zhao H., 2013, MNRAS, 432, 2846
  * Lüghausen et al. (2015) Lüghausen F., Famaey B., Kroupa P., 2015, Canadian Journal of Physics, 93, 232
  * Marinacci et al. (2018) Marinacci F., et al., 2018, MNRAS, 480, 5113
  * Martín-Navarro et al. (2018) Martín-Navarro I., Vazdekis A., Falcón-Barroso J., La Barbera F., Yıldırım A., van de Ven G., 2018, MNRAS, 475, 3700
  * McAlpine et al. (2016) McAlpine S., et al., 2016, Astronomy and Computing, 15, 72
  * McDermid et al. (2015) McDermid R. M., et al., 2015, MNRAS, 448, 3484
  * McGaugh (2004) McGaugh S. S., 2004, ApJ, 609, 652
  * McGaugh (2005) McGaugh S. S., 2005, ApJ, 632, 859
  * McGaugh (2012) McGaugh S. S., 2012, AJ, 143, 40
  * McGaugh et al. (2000) McGaugh S. S., Schombert J. M., Bothun G. D., de Blok W. J. G., 2000, ApJ, 533, L99
  * Milgrom (1983a) Milgrom M., 1983a, ApJ, 270, 365
  * Milgrom (1983b) Milgrom M., 1983b, ApJ, 270, 371
  * Milgrom (1999) Milgrom M., 1999, Physics Letters A, 253, 273
  * Milgrom (2008) Milgrom M., 2008, arXiv e-prints, p. arXiv:0801.3133
  * Milgrom (2010) Milgrom M., 2010, MNRAS, 403, 886
  * Milgrom (2014) Milgrom M., 2014, Scholarpedia, 9, 31410
  * Naab & Ostriker (2017) Naab T., Ostriker J. P., 2017, ARA&A, 55, 59
  * Nagesh et al. (2021) Nagesh S. T., Banik I., Thies I., Kroupa P., Famaey B., Wittenburg N., Parziale R., Haslbauer M., 2021, Canadian Journal of Physics, 99, 607
  * Naiman et al. (2018) Naiman J. P., et al., 2018, MNRAS, 477, 1206
  * Nelan et al. (2005) Nelan J. E., Smith R. J., Hudson M. J., Wegner G. A., Lucey J. R., Moore S. A. W., Quinney S. J., Suntzeff N. B., 2005, ApJ, 632, 137
  * Nelson et al. (2018) Nelson D., et al., 2018, MNRAS, 475, 624
  * Nelson et al. (2019) Nelson D., et al., 2019, Computational Astrophysics and Cosmology, 6, 2
  * Nipoti et al. (2007) Nipoti C., Londrillo P., Ciotti L., 2007, ApJ, 660, 256
  * Pazy (2013) Pazy E., 2013, Phys. Rev. D, 87, 084063
  * Pillepich et al. (2018a) Pillepich A., et al., 2018a, MNRAS, 473, 4077
  * Pillepich et al. (2018b) Pillepich A., et al., 2018b, MNRAS, 475, 648
  * Pipino & Matteucci (2008) Pipino A., Matteucci F., 2008, A&A, 486, 763
  * Planck Collaboration et al. (2014) Planck Collaboration et al., 2014, A&A, 571, A1
  * Planck Collaboration et al. (2016) Planck Collaboration et al., 2016, A&A, 594, A13
  * Recchi et al. (2009) Recchi S., Calura F., Kroupa P., 2009, A&A, 499, 711
  * Renaud et al. (2016) Renaud F., Famaey B., Kroupa P., 2016, MNRAS, 463, 3637
  * Ricciardelli et al. (2010) Ricciardelli E., Trujillo I., Buitrago F., Conselice C. J., 2010, MNRAS, 406, 230
  * Rodriguez-Gomez et al. (2015) Rodriguez-Gomez V., et al., 2015, MNRAS, 449, 49
  * Rosdahl et al. (2013) Rosdahl J., Blaizot J., Aubert D., Stranex T., Teyssier R., 2013, MNRAS, 436, 2188
  * Sanders (1990) Sanders R. H., 1990, A&ARv, 2, 1
  * Sanders (2008) Sanders R. H., 2008, MNRAS, 386, 1588
  * Saracco et al. (2017) Saracco P., Gargiulo A., Ciocca F., Marchesini D., 2017, A&A, 597, A122
  * Schaye et al. (2015) Schaye J., et al., 2015, MNRAS, 446, 521
  * Smolin (2017) Smolin L., 2017, Phys. Rev. D, 96, 083523
  * Springel et al. (2018) Springel V., et al., 2018, MNRAS, 475, 676
  * Stolte et al. (2014) Stolte A., et al., 2014, ApJ, 789, 115
  * Suzuki et al. (2022) Suzuki T. L., et al., 2022, arXiv e-prints, p. arXiv:2206.14238
  * Teyssier (2002) Teyssier R., 2002, A&A, 385, 337
  * Thomas et al. (2005) Thomas D., Maraston C., Bender R., Mendes de Oliveira C., 2005, ApJ, 621, 673
  * Thomas et al. (2010) Thomas D., Maraston C., Schawinski K., Sarzi M., Silk J., 2010, MNRAS, 404, 1775
  * Thomas et al. (2017) Thomas G. F., Famaey B., Ibata R., Lüghausen F., Kroupa P., 2017, A&A, 603, A65
  * Thomas et al. (2018) Thomas G. F., Famaey B., Ibata R., Renaud F., Martin N. F., Kroupa P., 2018, A&A, 609, A44
  * Trujillo et al. (2007) Trujillo I., Conselice C. J., Bundy K., Cooper M. C., Eisenhardt P., Ellis R. S., 2007, MNRAS, 382, 109
  * Trujillo et al. (2014) Trujillo I., Ferré-Mateu A., Balcells M., Vazdekis A., Sánchez-Blázquez P., 2014, ApJ, 780, L20
  * Verlinde (2017) Verlinde E., 2017, SciPost Physics, 2, 016
  * Weidner & Kroupa (2006) Weidner C., Kroupa P., 2006, MNRAS, 365, 1333
  * Weinberger et al. (2017) Weinberger R., et al., 2017, MNRAS, 465, 3291
  * Wiklind et al. (2008) Wiklind T., Dickinson M., Ferguson H. C., Giavalisco M., Mobasher B., Grogin N. A., Panagia N., 2008, ApJ, 676, 781
  * Wittenburg et al. (2020) Wittenburg N., Kroupa P., Famaey B., 2020, ApJ, 890, 173
  * Wright et al. (1988) Wright G. S., Joseph R. D., Robertson N. A., James P. A., Meikle W. P. S., 1988, MNRAS, 233, 1
  * Wright et al. (2019) Wright R. J., Lagos C. d. P., Davies L. J. M., Power C., Trayford J. W., Wong O. I., 2019, MNRAS, 487, 3740
  * Wuyts et al. (2010) Wuyts S., Cox T. J., Hayward C. C., Franx M., Hernquist L., Hopkins P. F., Jonsson P., van Dokkum P. G., 2010, ApJ, 722, 1666
  * Yan et al. (2017) Yan Z., Jerabkova T., Kroupa P., 2017, A&A, 607, A126
  * Yan et al. (2019) Yan Z., Jerabkova T., Kroupa P., Vazdekis A., 2019, A&A, 629, A93
  * Yan et al. (2021) Yan Z., Jeřábková T., Kroupa P., 2021, A&A, 655, A19
  * Zonoozi et al. (2021) Zonoozi A. H., Lieberz P., Banik I., Haghi H., Kroupa P., 2021, MNRAS, 506, 5468

##  Appendix A Fit Parameters

The fit parameters for Eq.10 are provided in Table 4.

Table 4: Fit parameters for Eq.10 $R_{\rm initial}$ | $A$ | $B$
---|---|---
(kpc) |  |
50 | 40747.40 | -0.55
100 | 6102.39 | -0.43
200 | 2989.50 | -0.37
300 | 227.87 | -0.25
500 | 6253.16 | -0.37

##  Appendix B SFTs for $\Lambda$CDM simulations

Using the TNG50-1 and TNG100-1 simulation of the Illustris The Next Generation
(TNG) project (Weinberger et al., 2017; Naiman et al., 2018; Springel et al.,
2018; Pillepich et al., 2018a, b; Nelson et al., 2018; Marinacci et al., 2018;
Nelson et al., 2019) and the Ref-L100N1504 simulation (hereafter EAGLE100) of
the EAGLE project (Crain et al., 2015; Schaye et al., 2015; McAlpine et al.,
2016), the star-formation timescales, SFTs, of the formed model galaxies are
investigated within the $\Lambda$CDM framework. These projects are sets of
different self-consistent cosmological galaxy formation simulations. The TNG
project is consistent with the Planck-2015 (Planck Collaboration et al., 2016)
results: $H_{0}=67.74\,\rm{km\,s^{-1}\,Mpc^{-1}}$,
$\Omega_{\mathrm{b},0}=0.0486$, $\Omega_{\mathrm{m},0}=0.3089$,
$\Omega_{\Lambda,0}=0.6911$, $\sigma_{8}=0.8159$, and $n_{\mathrm{s}}=0.9667$.
Here, the TNG50-1 has a box side of 51.7 co-moving Mpc (cMpc) an initial
baryonic mass resolution of $m_{\mathrm{b}}=8.5\times 10^{4}\,\rm{M_{\odot}}$,
and a dark matter particle mass of $m_{\mathrm{dm}}=4.5\times 10^{5}\
\rm{M_{\odot}}$. The TNG100-1 simulation has a box side of $110.7$ co-moving
Mpc (cMpc), $m_{\mathrm{b}}=1.4\times 10^{6}\,\rm{M_{\odot}}$, and
$m_{\mathrm{dm}}=7.5\times 10^{6}\,\rm{M_{\odot}}$. An overview of the
physical and numerical parameters is given in table 1 of Nelson et al. (2019).
The EAGLE project assumes the Planck-2013 (Planck Collaboration et al., 2014)
measurements with the cosmological parameters being
$H_{0}=67.77\,\rm{km\,s^{-1}\,Mpc^{-1}}$, $\Omega_{\mathrm{b},0}=0.04825$,
$\Omega_{\mathrm{m},0}=0.307$, $\Omega_{\Lambda,0}=0.693$,
$\sigma_{8}=0.8288$, and $n_{\mathrm{s}}=0.9611$ (see also table 1 of Schaye
et al., 2015). The EAGLE100-1 run has a box size of $100$ cMpc, an initial
baryonic particle mass of $m_{\mathrm{b}}=1.81\times 10^{6}\rm{M_{\odot}}$,
and a dark matter particle mass of $m_{\mathrm{dm}}=9.70\times
10^{6}\rm{M_{\odot}}$ (table 2 of Schaye et al., 2015).

We select present-day ($z=0$) subhaloes with a stellar mass
$M_{*}>10^{9}\,\rm{M_{\odot}}$ and a star formation rate of
$\it{SFR}<\mathrm{10^{-3}}\,\rm{M_{\odot}\,yr^{-1}}$. In addition, we remove
subhaloes with a Subhaloflag value of 0 in the TNG simulation, ensuring that
only galaxies or satellites of a cosmological origin are considered. This
gives a final sample of $418$ (45 central galaxies), $4384$ (803 central
galaxies) and $2916$ (373 central galaxies) galaxies in the TNG50-1, TNG100-1
and EAGLE100 run, respectively.

The SFT in the cosmological $\Lambda$CDM simulations is calculated by using
three different methods. In the first approach, we extract the formation time
and the masses of the stellar particles at $z=0$ of the above selected
galaxies. The SFT, $\Delta\tau$, for the TNG100-1 and EAGLE100 data is
calculated using the same method as described in Section 3 where the FWHM of
the SFH gives the duration for which the bulk of the star formation takes
place. The star formation rate (SFR) is calculated by separating all stellar
particles in bins of $\delta t$ = 10 Myr according to their age, summing up
the stellar mass in every time bin and dividing this by the length of $\delta
t$. Their SFT–mass relation is shown in the top panels of Fig. 9.

In the second method, we extract again the stellar particles at $z=0$ of the
selected galaxies but using the mass of the star particle at the time of its
formation in order to calculate the SFT (bottom panels of the Fig. 9). We use
the same method as described above to calculate the SFRs and obtain the
SFT–mass relation.

The mean of $\Delta\tau$ for all three simulation runs lies around 5 Gyr and a
trend in $M_{*}$ is absent for both methods (except for EAGLE) but is in
complete disagreement to observed ETGs (Thomas et al., 2005; Recchi et al.,
2009; McDermid et al., 2015). The Wilcoxon Signed-Rank test (Bhattacharyya &
Johnson, 1977) shows the SMoC models to disagree with the observed SFTs at
more than 5 sigma confidence because the latter lie consistently and
significantly below the former.

Figure 9:  SFT–mass relation of model galaxies formed in the TNG100-1
(_left_), TNG50-1 (_centre_) and EAGLE100 (_right_) (SMoC) simulation by using
the method 1 (top panels) and method 2 (bottom panels) compared with
observationally constrained relations for ETGs. The orange coloured dots
labelled as CG are the central galaxies in each simulation. There are a total
of 30 bins on the x- and y-axis and the binsize is 0.08 dex. Figure 10: Same
as Fig. 9 but using method 3 as described in Appendix B.

In the third method, we trace the selected galaxies back in cosmic time by
using the merger tree catalogs (Rodriguez-Gomez et al., 2015) and extract
their SFR from the Subfind Subhalo catalogs at different snapshots. For this,
we only select subhaloes with a dark matter mass of $M_{\mathrm{dm}}>0$ in
order to guarantee that these galaxies are included in the merger trees
catalogs. This gives $417$ and $4296$ subhaloes in the TNG50-1 and TNG100-1
run, respectively. Their SFT–mass relation is shown in Fig. 10.

The third method cannot be applied for the EAGLE run because only 29 snapshots
outputs (see table C.1 of McAlpine et al. (2016) are publicly available and a
fit cannot be made to get the SFT.

Given the test performed here which relies on quantifying the present-day (z =
0) galaxies in the cosmological simulations, the results unambiguously show
that the models do not have the observed downsizing. This is in contradiction
to other $\Lambda$CDM models (see section 2.3 of Wright et al. 2019) where
they use their own quantification of measuring quenching times of model
galaxies. All 3 methods result in a 5 sigma tension (Figs. 9 and 10).

##  Appendix C Collapse time

The MOND free-fall time is the time taken for the gas cloud to collapse under
its own gravitational attraction to form stars. This is calculated using Eq. 6
(see Zonoozi et al. 2021). We compare the time taken for the first stellar
particle to form in the simulation, $t_{\rm start}$, with this theoretical
value (Fig. 11). The $t_{\mathrm{start}}$, and the MOND free fall time are
comparable in our simulation demonstrating consistency between the numerical
and analytical results.

Figure 11: Collapse time – cloud density relation for all models with $R_{\rm
initial}$ = 500 kpc. The purple triangles are the time when the first stellar
particle is formed in the simulation, $t_{\rm start}$ and the red triangles
are the MOND free-fall time (Eq. 6).

##  Appendix D Percentage change

The percentage change (PC) quantifies the change from an initial value to the
new value and expresses the change as an increase or decrease,

$\mathrm{PC}=\frac{(X_{2}-X_{1})}{|X_{1}|}\times 100.$ (16)

Here, $X_{1}$ is the initial value and $X_{2}$ is the new value. If the PC
obtained is a positive value then it is quantified as an increase from the
initial value and if the obtained PC is negative then it is quantified as a
decrease from the initial value. In Table 3, the results obtained for the
simple feedback (SF) models are used as the initial values and the results
obtained for the complex feedback (CF) are used as the new values.