<figure><img src="image_1.png" /></figure>

# The Formation of Exponential Disk Galaxies in MOND

Nils Wittenburg<fnref n="1" />,*, Pavel Kroupa<fnref n="1" />,<fnref n="2" />, and Benoit Famaey<fnref n="3" />

Received 2019 October 18; revised 2020 January 15; accepted 2020 January 18; published 2020 February 26

## Abstract

The formation and evolution of galaxies are highly dependent on the dynamics of stars and gas, which is governed by the underlying law of gravity. To investigate how the formation and evolution of galaxies take place in Milgromian gravity (MOND), we present full hydrodynamical simulations with the Phantom of Ramses code. These are the first-ever galaxy formation simulations done in MOND with detailed hydrodynamics, including star formation, stellar feedback, radiative transfer, and supernovae. These models start from simplified initial conditions, in the form of isolated, rotating gas spheres in the early universe. These collapse and form late-type galaxies obeying several scaling relations, which was not a priori expected. The formed galaxies have a compact bulge and a disk with exponentially decreasing surface mass density profiles and scale lengths consistent with observed galaxies, as well as vertical stellar mass distributions with distinct exponential profiles (thin and thick disk). This work thus shows for the first time that disk galaxies with exponential profiles in both gas and stars are a generic outcome of collapsing gas clouds in MOND. These models have a slight lack of stellar angular momentum because of their somewhat compact stellar bulge, which is connected to the simple initial conditions and the negligible later gas accretion. We also analyze how the addition of more complex baryonic physics changes the resulting main properties of the models and find this to be negligibly so in the Milgromian framework.

*Unified Astronomy Thesaurus concepts: Galaxies (573); Disk galaxies (391); Hydrodynamical simulations (767); Gravitation (661); Non-standard theories of gravity (1118)*

## 1. Introduction

Galaxies, groups of galaxies, and clusters of galaxies are the largest gravitationally bound baryonic structures of our universe. Yet, the formation and evolution of galaxies are still among the most challenging phenomena in astrophysics, partly because it is, even today, not possible to perform high-resolution simulations that contain all dynamical and hydrodynamical effects from galaxy scales down to the scale of a single star.

On the other hand, concerning their general properties, galaxies seem to be simple objects. The vast majority of galaxies with stellar mass above ≈10¹⁰ M☉ are late-type galaxies at all redshifts (e.g., Delgado-Serrano et al. 2010; Tamburri et al. 2014), with pure ellipticals representing less than 5% (e.g., Delgado-Serrano et al. 2010). Disney et al. (2008) analyzed hundreds of galaxies with respect to six of their properties, which describe the galaxy completely. Correlation matrices for all parameters were analyzed with a principal component analysis. Surprisingly, they found that all parameters are strongly correlated and that they are determined by only one of them, which they could not determine with their method; thus, only one principal component is evident from observations. Furthermore, the mass discrepancy acceleration relation/radial acceleration relation (MDAR/RAR; McGaugh 2004; Lelli et al. 2017), which generalizes several other galaxy scaling relations such as the mass asymptotic speed relation (MASR, aka baryonic Tully-Fisher relation, BTFR)<fnref n="4" /> (Sanders 1990; McGaugh et al. 2000), the “dichotomy” between HSB and LSB galaxies (de Blok & McGaugh 1997; Tully & Verheijen 1997), the central density relation (Lelli et al. 2013, 2016b), the “baryon-halo

conspiracy” (van Albada & Sancisi 1986), Renzo’s rule (Sancisi 2004), the Faber-Jackson relation (Faber & Jackson 1976), and the $σ_⋆ - V_{H1}$ relation (Serra et al. 2016), while having little to no intrinsic scatter, strengthens the argument that galaxies follow a tight and seemingly simple law. This is further emphasized by late-type (i.e., disk) galaxies lying on a star-forming main sequence with a small dispersion, over a broad range of redshifts (Speagle et al. 2014).

Especially late-type galaxies show ordered structures, regardless of, e.g., the total mass or morphology, as the majority of them have an exponentially decreasing stellar surface mass density profile. This phenomenon has been known for more than 40 yr (Freeman 1970), but its origin is still not understood. For example, Herpich et al. (2017) show that these profiles might appear because of the dynamics within the disk, but there are also explanations that use the initial conditions of the formation of the galaxy for their reasoning (Fall & Efstathiou 1980; Dalcanton et al. 1997; Dutton 2009).

The order and simplicity of galaxies, among other observations on galactic scales (Kroupa 2015), are difficult to explain in the standard model of cosmology, as the evolution of galaxies is, in that framework, based on many subsequent mergers that happen stochastically. Fossils left over from this merging history orbit the galaxy as satellite galaxies, but the observed very significant phase-space correlation of satellite galaxies around the Milky Way (MW; Kroupa et al. 2005; Pawlowski et al. 2012; Pawlowski & Kroupa 2020), Andromeda (Metz et al. 2007; Ibata et al. 2013), and Cen A (Müller et al. 2018) is inconsistent with the observed satellites being such fossils (Pawlowski 2018). An explanation in standard dark matter models for the mutual alignment of the MW and Andromeda satellite population has not been found, but it appears to arise naturally in Milgromian dynamics (Banik et al. 2018b;

<footnote n="1">¹ University of Bonn, Helmholtz-Institut für Strahlen- und Kernphysik, Nussallee 14-16, D-53115 Bonn, Germany; nwittenbrug@astro.uni-bonn.de</footnote>

<footnote n="2">² Charles University in Prague, Faculty of Mathematics and Physics, Astronomical Institute, V Holešovičkách 2, CZ-180 00 Praha 8, Czech Republic</footnote>

<footnote n="3">³ Université de Strasbourg, Observatoire Astronomique de Strasbourg, CNRS UMR 7550, 11 rue de l'Université, F-67000 Strasbourg, France</footnote>

<footnote n="4">⁴ To prevent confusion about the velocity that has to be used here, we choose to use the name MASR.</footnote>
Bílek et al. 2018). Moreover, given the nondetection of dark matter particles despite a significant effort to find them (Bertone & Tait 2018) and astronomical arguments suggesting that dark matter does not exist (Kroupa et al. 2010; Kroupa 2012, 2015), it is necessary to also begin to investigate whether non-dark-matter-based models lead to objects that resemble observed galaxies. This is, of course, an immense task, which we merely start taking on here, with rather simplified initial configurations allowing us to investigate some very generic trends. As we will show, some interesting regularities will nevertheless appear.

On galaxy scales Milgromian dynamics (MOND; Milgrom 1983, 2014; Famaey & McGaugh 2012) accounts for the observations regarding the shape of the rotation curves of galaxies and galaxy scaling relations such as the MASR (BTFR; Sanders 1990; McGaugh et al. 2000; McGaugh 2012). As the baryonic content is the only source for the gravitational potential, MOND is also capable of predicting the rotation curve of a galaxy solely from its baryonic content, and furthermore, due to its predictive nature, it led to the discovery of new scaling relations, e.g., the MDAR (Sanders 1990; McGaugh 2004, 2005, 2012, 2014), which is mathematically the same as the RAR (see also Equation (2) of Milgrom 1983 as a first prediction of the MDAR/RAR).

For completeness we note that Di Cintio & Lelli (2016), Ludlow et al. (2017), Keller & Wadsley (2017), and Navarro et al. (2017) argue that these relations can be reproduced within the standard Lambda cold dark matter (ΛCDM) models, although they were not predicted and despite certain questions remaining open about their scatter (see, e.g., Di Cintio & Lelli 2016; Desmond 2017a, 2017b; Lelli et al. 2017).

Although the MOND framework is analytically successful on galaxy scales, it has hardly been tested in simulations. Part of the reason is its nonlinearity and therefore the lack of simulation codes that are able to calculate the MOND potential. Especially, self-consistent simulations with star formation and a full hydrodynamical treatment of gas were impossible to perform until the recent development of the Phantom of Ramses (POR; Lüghausen et al. 2015) and RAyMOND (Candlish et al. 2015) codes, which are customized versions of the publicly available simulation code Ramses (Teyssier 2002).

Previous applications of POR to state-of-the art problems are available in the very high resolution simulations of Antennae-like galaxies (Renaud et al. 2016) and (without hydro) in the simulation of the Sagittarius satellite galaxy and its stream by Thomas et al. (2018), as well as the computations of streams from globular clusters (Thomas et al. 2018) and the computation of the Andromeda–MW encounter by Bílek et al. (2018). Prior to this, Milgromian simulations with gas (Tiret & Combes 2008; Combes 2014) were carried out by using sticky particles and a MOND Poisson solver (Brada & Milgrom 1999; Tiret & Combes 2007). These publications were concerned with the evolution of preexisting disk galaxies and/or their satellites, but none were concerned with the formation of galaxy disks.

This work presents the first fully self-consistent hydrodynamical simulations of galaxy formation done in the Milgromian framework. It is, however, worth insisting here that the MOND paradigm, as it stands today, is essentially mute on cosmology, so that the paradigm must necessarily be incomplete and needs to be embedded into a self-consistent cosmological model. Nevertheless, although details will of course depend on the

specific parent theory, we can study some general consequences of a MOND-like force in cosmology. First, because of the stronger force, dynamical measures of the cosmic mass density will be overestimated, just as in galaxies. Also, MOND forms structures more rapidly than Newtonian gravity with the same source perturbation (Llinares et al. 2008). In the early universe, perturbations would (in principle) not grow because the baryons are coupled to the photon fluid. Since the mass density is lower in MOND, matter domination would occur later than in ΛCDM. Consequently, MOND structure formation initially would lag behind ΛCDM. However, as the influence of the photon field would then decline and perturbations begin to enter the MOND regime, structure formation would rapidly accelerate. Large galaxies could form as early as $z \approx 10$ from the collapse of gas clouds, which we model hereafter. Note that this means that there would be fewer mergers than in ΛCDM at intermediate redshifts, and in those rare mergers and much more numerous encounters between spirals, tidal dwarf galaxies would be formed and survive more easily than in ΛCDM (Sanders 1998). This may indeed plausibly be the source of satellite galaxies and their observed phase-space correlation (Okazaki & Taniguchi 2000; Banik et al. 2018b; Bílek et al. 2018; Javanmardi & Kroupa 2020).

At this point we want to focus on one important difference between simulations done in ΛCDM and MOND. Simulations done in the ΛCDM framework depend strongly on the feedback processes that are implemented in the code that is used, which are assumed by the community to be very uncertain. For example, supernovae, active galactic nucleus activity, and ram pressure introduce baryonic outflows from the galaxy that lead to a lower value of the baryon-to-DM fraction within the galaxy compared to the cosmic mean value. Therefore, the implementation of these processes and their fine-tuning are very important for the shape and properties of the galaxy. However, in MOND this should not be the case, because there is no large gravitating dark matter halo with its own dynamical history and no reason to eject material to obtain the right baryon fraction. We will show that in MOND it is indeed not necessary to introduce complex baryonic physics, as the general properties of the galaxies will not be significantly different. Furthermore, in ΛCDM simulations, angular momentum can be exchanged between the baryons and the dark halo. Again, this is not the case in simulations done in MOND.

A short introduction into the MOND framework and the numerical methods that were used is given in Section 2, and afterward the models are described in Section 3. Section 4 contains a detailed analysis of the formed galaxies. In Section 5 the computations are compared with observed scaling relations (MDAR and MASR), and in Section 6 the results are discussed and further work will be outlined that should be done in the future.

## 2. Theory and Numerical Methods

### 2.1. Modified Newtonian Dynamics (MOND)

Since galaxies show systematic and tightly correlated dynamical deviations from Newtonian behavior, here we assume this to be the result of effective gravitation deviating from the classical Newtonian law. This is plausible because the law of universal gravitation was empirically motivated only by solar system objects, in the case of both Newton (Newton 1687) and Einstein (Einstein 1916). We follow the empirical finding
that galaxies follow scale-invariant dynamics in the low-acceleration limit. Scale invariance is the symmetry under which the gravitational acceleration and dynamical (actual) acceleration transform in the same way under scaling of lengths and time by a factor $\lambda$. Thus, if the actual acceleration $g$ (which scales as $\lambda^{-1}$) is a function of the Newtonian gravitational acceleration, $g_N$ (which scales as $\lambda^{-2}$), their relation has to be of the form shown in Equation (1). As scale invariance only leads to proportionality in Equation (1), $a_0$ needs to be normalized such that equality holds (Milgrom 1983, 2009; Wu & Kroupa 2015). Thus, when scale invariance applies, the true gravitational acceleration, $g$, is given by

$$g = \sqrt{a_0 g_N}, \qquad (1)$$

where $a_0 \approx 10^{-10} \text{ m s}^{-2} \approx 3.7 \text{ pc Myr}^{-2}$ is Milgrom's constant and $g_N$ is the Newtonian acceleration.

Equation (1) is the basic equation underlying the MOND paradigm, and it is only valid in the low-acceleration deep-MOND limit $g \ll a_0 \approx 10^{-10} \text{ m s}^{-2}$. To also encapsulate the symmetry breaking in the Newtonian regime, Milgrom's law is formulated as follows:

$$g = \nu \left( \frac{g_N}{a_0} \right) g_N, \qquad (2)$$

with $\nu(y)$, $y = g_N/a_0$ being the transition function, which is defined by its limits:

$$\nu(y) \rightarrow 1 \text{ for } y \gg 1 \text{ and } \nu(y) \rightarrow y^{-1/2} \text{ for } y \ll 1. \qquad (3)$$

Equation (1) correctly predicts galaxy scaling relations, e.g., the MASR (BTFR), and led to the discovery of new scaling relations like the MDAR (Sanders 1990; McGaugh 2004, 2005, 2012). Additionally, it is the best-fitting model to the aforementioned RAR, which can be considered to be the most general galaxy scaling relation to date. MOND can be formulated as a classical gravitation theory with a fully developed Lagrangian formalism (Bekenstein & Milgrom 1984; Famaey & McGaugh 2012), and it may be related to the quantum physics of the vacuum (Milgrom 1999; Smolin 2017).

## 2.2. Simulation Code

Ramses is a publicly available hydrodynamical N-body code based on the adaptive mesh refinement technique; hence, it is a grid-based code with a tree-like data structure so that the grid can be recursively refined on a cell-by-cell basis (Teyssier 2002). Simulations start with the coarsest grid, and the cells in this grid will split up into $2^{\text{dim}}$ child cells (octs) if certain conditions are fulfilled, e.g., if the density or particle number threshold for the cell is reached. This also works vice versa, so that octs can be destroyed if, e.g., the density becomes too low. Therefore, the resolution increases in areas of interest.

The basic steps for grid-based codes to update the positions and velocities for the particles are the following (Teyssier 2002):

1. The mass density $\rho$ on the mesh is calculated using a “cloud-in-cell” (CIC) interpolation scheme.

2. The potential $\phi$ on the mesh is calculated by solving the Poisson equation (see the description of POR below for the MOND case).

3. The acceleration on the mesh is calculated using a standard finite-difference approximation of the gradient.

4. The particles’ accelerations are computed using an inverse CIC interpolation scheme.

5. The particles’ velocities are updated according to their accelerations.

6. The particles’ positions are updated according to their velocities.

To be more precise, the Poisson equation is solved in Ramses by minimizing the residual $\Delta\phi - 4\pi G\rho$ using the Gauss-Seidel method, and the particles’ positions and velocities are updated using the leap-frog scheme with adaptive time steps.

It is important to note that Ramses uses a “one-way interface” scheme, i.e., the finer-level solutions are updated using the coarser-level solutions and not vice versa, which would be a “two-way interface” scheme. But this basic scheme is not exactly true for all simulations done in this work, as most of the simulations were computed with MOND (see Section 2.2.3).

As mentioned before, not only can the trajectories and accelerations of particles be computed with Ramses, but it also includes a hydrodynamical solver. Therefore, star formation and cooling/heating processes of a given gas content are available, as well as dynamical phenomena like shocks. In the code the Euler equations in their conservative form are solved using a second-order Godunov scheme with a Riemann solver (see Teyssier 2002 and references therein).

### 2.2.1. Star Formation and Sink Particles

Star formation is very important for the simulations done in this work, because all particles are born from the initial gas content.

Star formation in general occurs as a result of the fragmentation of a gas cloud and the formation of dense clumps. If these clumps are heavier than the Jeans mass, they should collapse under self-gravity and form stars. However, this mechanism is not implemented in codes that are made for cosmological and galactic scales for several reasons. The resolution of the finest grid would have to be smaller than the size of a star, which would increase the computation time drastically. Further, if the resolution is on that level, all physical processes that act on this scale would have to be implemented, e.g., a correct initial mass function (Kroupa et al. 2013; Jeřábková et al. 2018), the evolution of single stars, and, most importantly, star-star interactions and binary stars (e.g., Oh & Kroupa 2016), altogether amounting to a prohibitive computational cost.

In Ramses, star formation is implemented as follows. Stellar particles are only created if two criteria are met. First, the gas mass-volume-density, $\rho$, needs to exceed a user-defined threshold, and second, a dimensionless integer number, $n_*$, drawn from a Poisson distribution, has to be unequal to zero. The distribution has a mean value of $\rho_{SFR} d_x^3 dt / M_*$ and depends on the local star formation rate (SFR) density, $\rho_{SFR}$; the length of the time step, $dt$; the 1D size of the cell, $d_x$, in which the density threshold is exceeded; and $M_*$, which is the mass scale of newly born stars (either user defined or given by $M_* = \rho_{\text{thres}} d_x^3$). The SFR follows the Schmidt law, $\rho_{SFR} = \rho/t_* \propto \rho^{3/2}$, where $t_* = t_{\text{ff}} (\rho/\rho_{\text{thres}})^{-1/2}$ is the star formation timescale, which is proportional to the local freefall time, $t_{\text{ff}} = \sqrt{3\pi/32G\rho}$. The code checks at every time step whether any cell exceeds the density threshold, then $n_*$ is drawn for the respective cell, and if it is nonzero, a stellar
particle is created by converting the amount of gas, $n_{\star} M_{\star}$, into the mass of the stellar particle. In our work the mass of these particles exceeds $10^4 M_{\odot}$ owing to our resolution, so they may be viewed as clusters of stars rather than single stars.

All particles produced this way are disconnected from the hydrodynamical evolution and only interact with the gas through gravity and feedback, which starts after a certain delay or directly depending on the user's choice (here a delay of 10 Myr is applied).

Also, a different kind of particle can be produced within the simulation, namely, sink particles (sinks), which are meant to stabilize the simulation.

In galaxy-scale simulations it is nearly impossible to resolve single stars as mentioned before, so following the gravitational collapse of a gas cloud has to be artificially stopped at some point owing to the maximum resolution. But not resolving the Jeans length and mass can lead to artificial fragmentation of the gas in collapsing regions (Truelove et al. 1998). This can be avoided using different techniques. One is to introduce a barotropic equation of state, so that the gas heats up if a certain density threshold is met. The downside of this implementation is that the gas formations stay extended and are more vulnerable to disruption (e.g., through shocks or tidal stripping). Also because of this, A. Bleuler and R. Teyssier invented and implemented a new sink particle scheme into Ramses. The idea is that collapsing gas clouds are artificially stopped if a certain density is exceeded and the gas is condensed into a point mass, which is decoupled from the hydrodynamical evolution as mentioned above. So the user can choose the density threshold according to a maximum resolution or certain physical scale (for further details about the implementation and tests of the sink particle scheme see Bleuler & Teyssier 2014).

### 2.2.2. Cooling/Heating and Radiative Transfer

As explained before, Ramses uses a second-order Godunov scheme to solve the Euler equations in their conservative form. There is also a radiative transfer (RT) addition to Ramses called Ramses-RT (Rosdahl et al. 2013), which uses a first-order Godunov solver, which will only briefly be explained at the end of this section.

In the code itself the function, which contains radiative cooling and heating, is computed separately; therefore, the thermochemistry scheme can be changed without changing the whole hydrodynamical solver. It depends on the gas density, temperature, ionization states of the gas (e.g., helium, hydrogen), and metallicity, but in simulations without the addition of the RT scheme, collisional ionization equilibrium (CIE) is assumed. That means that the ionization states can be calculated with the temperature and the density alone, so they do not need to be tracked in the code. Cooling and heating of the gas are therefore computed by using tables that are included in the code, which describe the cooling and heating rates due to several physical processes. The heating term includes photo-ionization processes, and the radiative cooling processes are collisional excitation, collisional ionization, recombination, bremsstrahlung, and Compton scattering. The Sutherland & Dopita (1993) cooling model is used in Ramses with look-up tables in the temperature and metallicity plane, and also the tables from Courty & Alimi (2004) are used for the different cooling/heating processes.

In fact, there is no rigorous treatment of RT by using this scheme, because the ionization states are computed by assuming CIE so that the cooling and heating rates can be computed with the temperature and density as described above. In the end the temperature is updated, and with it the total energy density, by also taking the rest of the Euler equations into account.

As mentioned at the beginning of this section, there is an option called RT that can be enabled for simulations. By doing so, the simulation is run with a different hydrodynamical solver. With this addition photon fluxes are introduced and the code keeps track of the ionization states, because they are now correctly computed by using photons, collisions, and none- equilibrium thermochemistry. So radiative transfer between sources like stellar particles is implemented.

Using the RT option increases the memory requirement massively, because the conservative state vector $U = (\rho, \rho u, E, \rho Z)$, which stores the hydrodynamical properties of each cell (gas density $\rho$, momentum density $\rho u$, total energy density $E$, and metal mass density $\rho Z$), becomes $U = (\rho, \rho u, E, \rho Z, \rho x_{HII}, \rho x_{HeII}, \rho x_{HeIII}, N_i, F_i)$, where new variables, connected to the photons, are introduced (ionization fraction densities of hydrogen and helium $\rho x_{HII}, \rho x_{HeII}, \rho x_{HeIII}$ with, e.g., $x_{HII} = n_{HII}/n_{H}$, photon number density $N_i$, and photon Fluxes $F_i$), and this roughly increases the memory requirement by a factor of 3.5.

### 2.2.3. Phantom of Ramses (POR)

Every computation done in this work used a customized version of Ramses, called POR. POR (Lüghausen et al. 2015) is the only publicly available code, which is capable of running full hydrodynamical simulations with a MONDian Poisson solver. Part of the reason why this is the only publicly available code is that the generalization of the Poisson equation in MONDian theories is nonlinear, which makes it very challenging or impossible to implement into an N-body code. However, there is a formulation of MOND that can be implemented without changing the Poisson solver of the code. It is called the quasi-linear formulation of MOND (QUMOND; Milgrom 2010), which is derived from an action and obeys the standard conservation laws, with the generalized Poisson equation,

$$ \Delta\Phi(\mathbf{x}) = 4\pi G\rho_b(\mathbf{x}) + \nabla \cdot [\tilde{\nu}(|\nabla\phi|/a_0)\nabla\phi(\mathbf{x})], \quad (4) $$

$$ \Delta\Phi(\mathbf{x}) = 4\pi G(\rho_b(\mathbf{x}) + \rho_{ph}(\mathbf{x})). \quad (5) $$

Here $\rho_b(\mathbf{x})$ is the baryonic density; $\phi(\mathbf{x})$ is the Newtonian potential, which fulfills the standard Poisson equation $\Delta\phi(\mathbf{x}) = 4\pi G\rho_b(\mathbf{x})$; $\Phi(\mathbf{x})$ is the total gravitational potential; and $\tilde{\nu}$ is the transition function between the Newtonian and the MOND regime, which is $\nu$ from Equation (2) minus 1 with the limits: $\tilde{\nu}(y) \to 0$ if $y \gg 1$ (Newtonian regime) and $\tilde{\nu}(y) \to y^{-1/2}$ if $y \ll 1$ (MOND regime) with $y = g_N/a_0$. Several functions that fulfill this criterion have been used in the literature (see Lüghausen et al. 2015), but here

$$ \tilde{\nu}(y) = -\frac{1}{2} + \frac{1}{2} \sqrt{1 + \frac{4}{y}} \quad (6) $$

is used.<fnref n="5" />

<footnote n="5">⁵ Note that Equation (6) corresponds to the “simple” interpolating function; see, e.g., Section 6.2 of Famaey & McGaugh (2012).</footnote>
Note that the second term on the right-hand side of Equation (4) was condensed into

$$ \rho_{\text{ph}}(\mathbf{x}) = \frac{\nabla \cdot [\tilde{\nu}(\nabla\phi/a_0)\nabla\phi(\mathbf{x})]}{4\pi G}, \quad (7) $$

so the generalized Poisson equation also visualizes its quasi-linearity. Equation (5) shows that the total gravitational potential in MOND depends on the baryonic density, and an additional term, which is called $\rho_{ph}$, also has the units of matter density and depends on the Newtonian potential alone. Therefore, the total gravitational potential can be written as $\Phi = \phi + \Phi_{ph}$, where it is divided into a Newtonian and a MONDian part. $\rho_{ph}$ is called “phantom dark matter” (PDM) density, but it is not a real matter distribution. It is a mathematical formulation to compare more easily MOND with Newtonian dynamics. It is also exactly the density that would be interpreted as dark matter in the standard cosmology framework.

As mentioned before, the generalized Poisson equation is not linear (through the computation of $\rho_{ph}$), so the Poisson solver in POR is different from the one in Ramses. However, the QUMOND formulation makes the implementation much easier, because the already-existing standard Poisson solver can be used and no new solver needs to be implemented. Therefore, to solve the complete Poisson equation for QUMOND, three steps were implemented in the code:

1. After the smoothed matter density distribution is calculated in a given cell from gas and particles, the standard Poisson equation,

$$ \Delta\phi(\mathbf{x}) = 4\pi G \rho_{\text{b}}(\mathbf{x}), \quad (8) $$

is solved to compute the Newtonian potential $\phi$ and its gradient.

2. The PDM density is calculated using Equation (7) (for the detailed scheme see Lüghausen et al. 2015).

3. With the matter density and the PDM density, the whole Poisson equation in QUMOND (Equation (5)) is solved to compute the total gravitational potential. The gradient of this potential at location $\mathbf{x}$ then yields the acceleration, $\mathbf{a}$, at $\mathbf{x}$, $\mathbf{a}(\mathbf{x}) = -\nabla\Phi(\mathbf{x})$.

We stress here that only the calculation of the potential is slightly changed, such that the total gravitational potential in MOND is calculated. Therefore, the hydrodynamical solver is unchanged but uses the total gravitational potential from the MONDian gravity solver.

### 3. The Models

The computations start from pure gas clouds, setting up identical models, except for the initial rotation velocity, radius, and mass. This was done to see whether it is possible to reproduce roughly the Hubble sequence (Hubble 1926) by changing the initial rotation velocity from 0 to a certain value. According to Disney et al. (2008), galaxies are simple objects, and their properties only depend on a single parameter, so if MOND contains the right description of gravitational dynamics,<fnref n="6" /> it should be possible to reproduce the Hubble sequence in simulations, that is, if all members of the sequence

can form in isolation from one gas cloud without further accretion of gas.

The initial distribution of gas is a uniform sphere with the same mass density for all simulations, which leads to a surface mass density distribution, $\Sigma_{\text{init}}(r_{\text{cyl}})$, that only depends on the initial spherical radius, $r_{\text{init}}$, and on the cylindrical radius, $r_{\text{cyl}}$,

$$ \Sigma_{\text{init}}(r_{\text{cyl}}) = 2\rho_{\text{init}}\sqrt{r_{\text{init}}^2 - r_{\text{cyl}}^2}, \quad (9) $$

with $\rho_{\text{init}}$ being the constant initial mass density. Hence, the models are initially morphologically unrelated to late-type galaxies, which are rotationally supported exponential disks. So if a disk galaxy forms, it forms because of the dynamics and not because of the initial distribution of gas.

#### 3.1. Initial Conditions

The initial conditions used may be viewed to be a first rough approximation of an early-universe gas cloud, which is gravitationally unstable. Later work will investigate initially turbulent gas clouds. We emphasize that the initial conditions used here are not based on standard cosmology, as it would be unlogical and unphysical to use initial conditions from the LCDM model in a model of galaxy formation without dark matter. In particular, since possible parent theories of MOND might be inherently nonlocal, using, e.g., the high-$z$ power spectrum in a MOND context is practically impossible, and we need to use some initial conditions, which are likely to not be very unphysical (i.e., pure gas clouds). Ultimately, by studying which initial conditions lead to which types of galaxy, we should be able to constrain MONDian cosmological theory.

Here it is assumed that neutral gas clouds form by a redshift of about 15–50 (Barkana & Loeb 2001). The calculations thus begin with a neutral cloud with the following initial properties (see also Table 1). The angular momentum of the gas cloud is expected to result from its internal motions as it cools and begins to collapse and from tidal torquing. In the past, the angular momentum of galaxies has been associated with tidal torquing (Efstathiou & Jones 1979; Wesson 1985; Voglis 1994; Catelan & Theuns 1996). In the models explored here a very simple initial law, which significantly differs from that of galaxies, is assumed in order to ascertain that the properties of the final galaxies are not assured through the specific initial conditions.

The models start as a constant-density sphere of gas, where the initial rotation velocity has a radial dependence,

$$ v_i(\mathbf{r}) = \eta \cdot \mathbf{r} \times \mathbf{I}_z, \quad (10) $$

with $\eta$ being an angular velocity parameter, $\mathbf{r}$ the radius vector, $v_i$ the initial rotation velocity vector, and $\mathbf{I}_z$ the unit vector in the direction of the $z$-axis. This is motivated by the first Larson relation (Larson 1981), which relates the velocity dispersion of a molecular gas cloud proportionally to the radius. Only three parameters are available to change the model: $\eta$ for the initial rotation velocity distribution, the initial 3D radius, $r_{\text{init}}$, and mass, $M_{\text{init}}$, of the sphere. Important to note here is that the initial velocity and density distribution are vastly different from those of observed disk galaxies.

The starting temperature of the gas is $T = 10^4$ K, and the density threshold for star formation to take place is set to $\rho_{\text{star}} = 0.1$ Hcm$^{-3}$, with a star formation efficiency of 5% (Dubois & Teyssier 2008). Also, all models evolve in isolation,

<footnote n="6">⁶ Right in the sense of reproducing the observed properties of galaxies.</footnote>
<table><caption>**Table 1**
Initial Conditions of All Models</caption><thead><tr><th>Model<br>Name/No.</th><th><i>M</i><sub>init</sub><br>(10<sup>9</sup><i>M</i><sub>☉</sub>)</th><th><i>r</i><sub>init</sub><br>(kpc)</th><th>η<br>(km s<sup>-1</sup> kpc<sup>-1</sup>)</th><th>Additions</th><th>Poisson<br>Solver</th></tr></thead><tbody><tr><td>M1/1</td><td>6.4</td><td>20</td><td>1.44</td><td>...</td><td>MOND</td></tr><tr><td>M1sn/2</td><td>||</td><td>||</td><td>||</td><td>sink, rt, sn</td><td>||</td></tr><tr><td>M1N/3</td><td>||</td><td>||</td><td>||</td><td>...</td><td>Newton</td></tr><tr><td>M1const/4</td><td>||</td><td>||</td><td>6.56 kpc</td><td>...</td><td>MOND</td></tr><tr><td>M1Zpoor/5</td><td>||</td><td>||</td><td>1.44</td><td>Z = 10<sup>-4</sup><i>z</i><sub>☉</sub></td><td>||</td></tr><tr><td>M1Zpoorsn/6</td><td>||</td><td>||</td><td>||</td><td>Z = 10<sup>-4</sup><i>z</i><sub>☉</sub><br>sink,rt,sn</td><td>||</td></tr><tr><td>M1111/7</td><td>||</td><td>||</td><td>||</td><td>res<sub>max</sub><br>= 468.75pc</td><td>||</td></tr><tr><td>M1113/8</td><td>||</td><td>||</td><td>||</td><td>res<sub>max</sub><br>= 117.19pc</td><td>||</td></tr><tr><td>M2/9</td><td>21.6</td><td>30</td><td>0.39</td><td>...</td><td>||</td></tr><tr><td>M2sn/10</td><td>||</td><td>||</td><td>||</td><td>sink, rt, sn</td><td>||</td></tr><tr><td>M2N/11</td><td>||</td><td>||</td><td>||</td><td>...</td><td>Newton</td></tr><tr><td>M3/12</td><td>||</td><td>||</td><td>0.96</td><td>...</td><td>MOND</td></tr><tr><td>M3sn/13</td><td>||</td><td>||</td><td>||</td><td>sink, rt, sn</td><td>||</td></tr><tr><td>M4/14</td><td>100.0</td><td>50</td><td>0.58</td><td>...</td><td>||</td></tr><tr><td>M4sn/15</td><td>||</td><td>||</td><td>||</td><td>sink, rt, sn</td><td>||</td></tr></tbody></table>

**Note.** Column (1): name and number of the model. Column (2): initial mass. Column (3): initial radius of the sphere. Column (4): initial velocity parameter $\eta$. Column (5): which options are added/changed (sn = supernova, sink = sink particles, rt = radiative transfer, Z = metallicity, res<sub>max</sub> = maximum resolution). Column (6): information about which Poisson solver is used. Note that $\eta$ for model M1const has the units km s<sup>-1</sup>, because it starts with a constant rotation velocity throughout the whole sphere of gas.

and there is no UV background radiation. The initial density is equal for all simulations with a value of $\rho_{init} = 7.83 \times 10^{-3}$ Hcm$^{-3}$. To mimic isolation, the density of the intergalactic medium outside the gas sphere is set to $\rho_{IGM} = 10^{-5} \times \rho_{init}$ such that after the collapse of the gas cloud further accretion from the ambient intergalactic medium onto the formed galaxy is not significant. The minimum mass of a stellar particle is $M_{\star} \approx 3 \times 10^4 M_{\odot}$. The initial metallicity is solar, except for M1Zpoor and M1Zpoorsn, and constant throughout the simulations. The maximum resolution for most simulations is 234.375 pc (for M1111 and M1113 see Table 1), which is the length of the smallest grid cell, and the size of the coarsest grid cell (the simulation box itself) is set to 960 kpc. Milgrom's constant is set to $a_0 = 1.12 \times 10^{-10} \text{ m s}^{-2}$ within the code. Supernovae and sink particles will be discussed later, when such models are explained (Section 4.1).

Fifteen different models were calculated for this work with the following properties (see also Table 1):

1. M1 assumes simple cooling/heating and no additional or more complex baryonic physics, such as supernovae. It starts with $M_{init} = 6.4 \times 10^9 M_{\odot}$, $r_{init} = 20$ kpc, and $\eta = 1.44 \text{ km s}^{-1} \text{ kpc}^{-1}$ and is simulated with the MON-Dian Poisson solver.

2. M1sn has the same initial conditions as M1 with the addition of supernovae, sink particles, and explicit radiative transfer.

3. M1N also is initially identical to M1, but it is simulated with the Newtonian Poisson solver.

4. M1const is set up with a different (constant) initial rotation law with $\eta = 0.1(\text{ kpc/r}) \times 65.6 \text{ km s}^{-1}$ and is otherwise identical to M1.

5. M1Zpoor and M1Zpoorsn are M1 and M1sn with lower metallicity $Z = 10^{-4} \times Z_{\odot}$.

6. M1113 and M1111 are M1 with a higher and smaller maximum refinement level, respectively, and therefore different maximum resolution, res<sub>max</sub> = 117.1875 pc and res<sub>max</sub> = 468.75 pc.

7. M2: $M_{init} = 21.6 \times 10^9 M_{\odot}$, $r_{init} = 30$ kpc, $\eta = 0.39 \text{ km s}^{-1} \text{ kpc}^{-1}$.

8. M3: $M_{init} = 21.6 \times 10^9 M_{\odot}$, $r_{init} = 30$ kpc, $\eta = 0.96 \text{ km s}^{-1} \text{ kpc}^{-1}$.

9. M4: $M_{init} = 100 \times 10^9 M_{\odot}$, $r_{init} = 50$ kpc, $\eta = 0.58 \text{ km s}^{-1} \text{ kpc}^{-1}$.

10. The name convention for M2, M3, and M4 is identical to M1, so sn means with sinks, supernovae, and radiative transfer, and N means with the Newtonian Poisson solver. Only $M_{init}$, $r_{init}$, and $\eta$ are different.

## 3.2. Galaxy Models

The initial collapse time is similar for all models, except for the Newtonian ones (this will be further discussed in Section 4.7), $\approx 0.5$ Gyr. So after approximately 0.5 Gyr the spheres collapse and form a rotating dense, thin disk with less dense gas surrounding it. All figures and videos of the simulations are made by projecting the relevant property of the galaxy onto the respective plane (edge-on/face-on view), i.e., a figure consists of equally sized cells ($0.46875 \times 0.46875 \text{ kpc}^2$), and for every figure cell the mass-weighted average of the stellar surface mass density and the gas density, along the axis that is not shown, is computed.

During the collapse, stellar particles form a thin, dense disk with a stellar halo surrounding the center and a less dense component in the outskirts. Very shortly after the formation of the disk-like structure, a spiral pattern arises, which is especially visible in the gaseous component (face-on).
<figure><img src="image_2.png" /><figcaption>**Figure 1.** M1sn after 8.45 Gyr shown face-on (top panels) and edge-on (bottom panels). Left panels: gas density, $\rho$, of M1sn color-coded as indicated by the scale on the right of the plot. Also, a void-like structure can be seen near the center, which is a “recent” supernova explosion (top left panel). Right panels: stellar surface mass density, $\Sigma$, plotted with the same color scheme as $\rho$. For movies of the formation of these galaxies, see YouTube entry “Formation of disk galaxies in MOND by Nils Wittenburg.”</figcaption></figure>

Until the end of the simulation (10 Gyr), the appearance of the galaxy does not change much, but it spreads radially, especially the less dense gaseous component. A spiral pattern in the face-on view of the gas remains, while the stellar component shows a very dense core (face-on) and a thin, dense disk surrounded by a flattened stellar halo (edge-on); see, e.g., Figure 1.

The left panel of Figure 2 shows the rotation curve of M1 (red crosses) and the “pristine” MOND rotation curve (Equation (2), purple squares) after 10 Gyr. The rotation curve is calculated by averaging the cylindrically radial accelerations of the stellar particles in radial bins of 500 pc and then computing the circular velocity according to $v_c = \sqrt{a\bar{r}}$, with $a$ being the acceleration and $r$ being the radius of the radial bin. The uncertainties are calculated from the scatter around the

average value of the acceleration in the bins and then propagated by Gaussian error propagation.

The “pristine” rotation curve is calculated by using Equation (2), assuming circular motion and calculating the Newtonian circular velocity beforehand. In general, the relation between $g$ and $g_N$ is equivalent for QUMOND and Milgrom’s law up to a curl field correction (Brada & Milgrom 1995), so every difference between the two rotation curves shows directly the effect of this curl field. Therefore, comparing the rotation curve of the simulation (QUMOND) with the “pristine” rotation curve (approximation that stems from Milgrom’s law) shows not only the effect of the curl field but also how significant this curl field correction is.

As is evident from Figure 2, the rotation curve rises steeply in the center ($r < 1$ kpc), decreases afterward for approximately
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Rotation curves of both models with $r_{init} = 20$ kpc shown in red, $v_{c,acc}$, and the “pristine” rotation curve of the respective model in purple. $v_{c,acc}$ shows the rotation velocity calculated from the mean radial acceleration in the respective radial bin. Left panel: M1; right panel: M1sn. See the text for further details on the computation and analysis.</figcaption></figure>

2 kpc, this being the Newtonian regime, and then becomes flat for the majority of the galaxy as indicated by the fit. To calculate the asymptotically flat rotation velocity of the galaxy, $v_{rot,flat}$, a constant fit function,

$$c(r) = v_{rot,flat}, \qquad (11)$$

is used here, which results in $v_{rot,flat} = 103.03 \pm 0.19$ km s$^{-1}$ for M1. The fit is performed for $3 < r/kpc < 20$ (all fits with Equation (11) are calculated from the beginning of the flat part of the rotation curve until the end of the stellar disk, which changes depending on the model).

The differences between the “pristine” curve and the measured curve are most prominent in two regions of the galaxy. The first is around 3 kpc away from the center, where the rotation velocity decreases significantly compared to the initial peak but is still outside the flat part of the curve, and the second region is the end of the stellar disk, where the measured rotation velocity increases again slightly in contrast to the behavior of the “pristine” curve.

The left panel of Figure 3 shows the surface mass density distribution, $\Sigma$, of stars, gas, and the sum of both after 10 Gyr. $\Sigma$ is calculated by adding the mass of every star or gas cell within a radial bin and dividing the total mass of the bin by its area. The uncertainties are Poisson uncertainties. Both the stellar and gaseous surface mass distributions are fitted by a simple exponential function to verify whether exponential disks form during the simulation. The function

$$\Sigma_{\text{exp}}(r) = n \times \exp(-r/r_e) \qquad (12)$$

is used, with $n$ being the normalization and $r_e$ the radial exponential scale length. For $\Sigma_{stars}$ only the stellar disk is used for the fit, meaning that the central part ($r < 2$ kpc) and radial bins with less than 10 stellar particles are neglected. $\Sigma_{gas}$ is fitted outside of the stellar disk to trace the outer mass distribution of the galaxy. Indeed, both components show a radial exponential profile, so the total surface mass density is essentially the sum of two exponential profiles. The inner part

of the galaxy is dominated by the stellar particles; therefore, $\Sigma_{tot}$ decreases according to $\Sigma_{stars}$ until it becomes smaller than the surface mass density of the gas. At that point the decrease becomes shallower and $\Sigma_{tot}$ follows $\Sigma_{gas}$. This is also reflected in the different exponential scale lengths of the stellar and gaseous distribution: $r_{e,stars} = 1.61 \pm 0.04$ kpc and $r_{e,gas} = 5.75 \pm 0.06$ kpc. The exponential scale lengths of every model are discussed in more detail in Figure 18 with additional properties of the models, while the numerical values of all fit parameters are shown in the Appendix (Table 4). The more extended gas disks result from the star formation condition not being fulfilled at larger radii such that the gas is not consumed there, while in the inner regions (the region occupied by the stellar particles) star formation has consumed the vast majority of the gas.

It should be noted here that observed galaxies show exponential stellar surface mass density profiles, so it is a major result that this very simple simulation without complex baryonic physics is showing a similar behavior. What is more, it is also noteworthy that all galaxy models form an exponential disk shortly after their collapse (roughly 0.5–1 Gyr after the computation begins). Additionally, we want to stress that the occurrence of exponential radial profiles in star-forming rotationally supported galaxies is still an unsolved problem. There are several different approaches attempting to account for the observed exponential surface density profiles: inclusion of scattering of stars in idealized models by Struck & Elmegreen (2017) generates thick and warm-to-hot stellar disks, while a phenomenological model of cloud disruption launching gas to large distances suggests the settling of the gas into an exponential disk that may form stars (Struck & Elmegreen 2018). Although potentially promising, such models of baryonic processes remain, to some degree, ad hoc.

After the discussion of M1, the differences between models computed with simple and more complex baryonic physics will be shown in the following.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Surface mass density profiles of both models with $r_{\text{init}} = 20$ kpc. Left panel: M1; right panel: M1sn. The red points correspond to the stellar surface mass density, the green ones correspond to the surface mass density distribution of the gas, and the blue ones show the total baryonic surface mass density distribution. The lines show the respective exponential fit to data.</figcaption></figure>

# 4. Results

In this section the resulting galaxies are analyzed.

## 4.1. Impact of Different Physical Processes

In order to understand the role of including more complex baryonic physics in galaxy formation and evolution within the MONDian framework, identical simulations were carried out with different physical processes involved.

The “simple” model (M1) was calculated with star formation and simple cooling/heating, while for the “complex” model (M1sn) star formation, radiative transfer, sink particle formation, and supernova explosions are enabled.

The density threshold for the sink particle formation was set to $2.77 \text{ Hcm}^{-3}$, such that the Jeans length is resolved by at least four grid cells. The parameters for the implementation of the supernovae were taken from Dubois & Teyssier (2008) and Teyssier et al. (2013). The supernova scheme in Ramses injects the kinetic part of the energy as spherical blast waves with the size of galactic superbubbles (here $r_{\text{SN}} = 150$ pc). If the energy would be injected only as thermal energy, then most of it would be radiated away very quickly and the effect of the supernovae on the kinematics of the stars and the gas would be minor.

After the collapse, both show a dense disk in the gaseous and stellar component with a low-density environment and a dense stellar core. Differences can be seen in the extent of the low-density environment in the gaseous component, and a spiral pattern is evident in the gas disk of M1 seen face-on, while M1sn does not show this pattern shortly after the collapse. The differences between both models most likely stem from supernova explosions, as they make the gas more turbulent, so that early on no spiral pattern can form in model M1sn. The gas is also spread wider, so low-density gas can be seen at larger radii.

After 10 Gyr, both models are still morphologically quite similar. Both show a spiral pattern in the gaseous component face-on, and the stellar components look nearly identical. The only major difference is the vertical extent of the gas, because

supernovae push the gas in every direction in M1sn, while for M1 the gas stays within the disk.

Figures 2–3 show the comparison of the analysis of the two simulations. The inner nonflat part of the rotation curves depicts the major difference between the two models, which can be explained by the effects of supernova explosions on the gas in M1sn, which heat the gas, thereby reducing the central deep collapse of the baryons. At a radius near 4 kpc the surface mass density increases for model M1sn, in contrast to the slow decrease that can be seen for model M1. Most of the supernovae take place in the center of the model, so the gas density is decreased significantly in this region, and therefore the surface mass density rises again to a certain degree outside of the center. This also leads to a difference in the gravitational potential and so to a difference in the rotational velocity. Also, this is the region where both models deviate most from their “pristine” curve. In general, the transition region between the flat part of the rotation curve and the decreasing flank after the peak is mostly influenced by the curl field correction, which will be underlined further in the following sections.

Apart from the different distributions of gas in the center, the surface mass density distributions are quantitatively in agreement. Both models show a steep decrease at the center for the first few kiloparsecs, after which $\Sigma_{\text{stars}}$ exponentially decreases. In the outer parts of the galaxy $\Sigma_{\text{gas}}$ becomes greater than $\Sigma_{\text{stars}}$, and it decreases also exponentially. So both components for both models show exponentially decreasing profiles except in the center. What is more, the exponential scale lengths are also slightly, but not substantially, different, because of the redistribution of gas due to supernovae, leading to $r_{\text{e,stars,M1}} = 1.61 \pm 0.04$ kpc, $r_{\text{e,stars,M1 sn}} = 1.82 \pm 0.06$ kpc and $r_{\text{e,gas,M1}} = 5.37 \pm 0.10$ kpc, $r_{\text{e,gas,M1 sn}} = 6.28 \pm 0.09$ kpc.

Although the supernovae have an obvious impact on the simulation, the values of the rotation velocity in the flat part of the two rotation curves are very comparable, $v_{\text{rot,flat,M1}} = 103.03 \pm 0.19 \text{ km s}^{-1}$ and $v_{\text{rot,flat,M1 sn}} = 99.98 \pm 0.09 \text{ km s}^{-1}$, the fits (Equation (11)) being obtained over $3 < r/\text{kpc} < 20$ and $6 < r/\text{kpc} < 15.5$, respectively.
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Comparison of the rotation curves between models with different rotation laws (top) and different values of η (bottom). From top left to bottom right: M1, M1const, M2, and M3. See the text for further details.</figcaption></figure>

All in all, both models are very similar, despite the
differences that a simulation with and without supernova
explosions has, e.g., supernovae lead to a larger extent of the
gaseous component, a higher gas fraction, and therefore lower
mass in stellar particles.

The next step is to see whether different initial values of the
rotation velocity parameter, η, or a different initial rotation law
change the properties and the morphology of a model.

## 4.2. *Comparison between Different Initial Rotation Laws*

In this section the impact of different initial rotation laws and
η will be examined.

The comparison is shown in Figures 4 and 5. Four models
are compared (M1, M1const, M2, M3), and the differences
between M1 and M1const and between M2 and M3 are
discussed. All initial parameters can be seen in Table 1. M2 and
M3 assume a more extended and more massive gas cloud than
M1, whereby M3 rotates faster than M2. M1const has a
constant rotation velocity throughout the cloud.

Differences can be seen already after the collapse of the
gaseous sphere. The collapse of M3 is not as smooth as the one
for M2, because three dense clumps form prior to the formation

of the whole galaxy, due to the higher rotation velocity than for
model M2. These clumps encounter each other and merge
asymmetrically, so an asymmetric, rotating gas distribution
forms around the dense disk. Also, the whole galaxy is more
extended (especially the stellar component), which is shown
clearly in the rotation curve and the surface mass density
distribution (Figures 4 and 5, respectively).

Despite these initial differences, the evolution of both
models is very similar. Both form a thin, rotating, dense disk
in both components.

M1const is initially different from all other models, because
the initial rotation velocity distribution is constant, $v_i = 6$
$6 \text{ km s}^{-1}$. Therefore, the collapse of the sphere at the beginning
of the simulation is more violent, because especially the infall
of gas from the outermost radii is not slowed down by a higher
rotation velocity. So, compared to M1, M1const is overall
smaller, which is in agreement with the findings from the
comparison between M2 and M3.

Not only is the initial collapse different depending on the
initial rotation velocity, but also the radial extent of star
formation differs, as the rotation curve of M3 reaches radii that
are 10 kpc larger than the one for M2, while the stellar disk of
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Comparison of the surface mass density distributions between models with different rotation laws (top) and different values of $\eta$ (bottom). From top left to bottom right: M1, M1const, M2, and M3. See the text for further details.</figcaption></figure>

M1 reaches out approximately 5 kpc farther than the disk of M1const. M1const and M2 have higher peak rotation velocities compared to M1 and M3, respectively, and the effect of the curl field correction is most prominent in the region directly after the peak and before the flat part, although the rotation curve of M3 is nearly identical to its “pristine” curve. The flat rotation velocities are very comparable between M1 and M1const and between M2 and M3, but not identical, as the mass within the stellar disk varies slightly owing to different extents and gas densities, $v_{\text{rot,flat,M2}} = 142.41 \pm 0.46$ km s$^{-1}$, $v_{\text{rot,flat,M3}} = 140.35 \pm 0.21$ km s$^{-1}$, $v_{\text{rot,flat,M1}} = 103.03 \pm 0.19$ km s$^{-1}$, and $v_{\text{rot,flat,M1const}} = 104.50 \pm 0.14$ km s$^{-1}$. The corresponding ranges for the fits (Equation (11)) are $3 < r/kpc < 20$ (M1), $10 < r/kpc < 15.25$ (M1const), $10 < r/kpc < 18.25$ (M2), and $20 < r/kpc < 26.75$ (M3).

Similar to the models shown in the last section, the surface mass density profiles for the stellar and the gaseous component are exponentially decreasing. Also, the difference in the extent of the stellar component is visible via the slope of the fit of this component, while the fit of $\Sigma_{\text{gas}}$ for every model shows that the surface mass density distribution of the gas varies as well, which is shown by the exponential scale lengths in Table 2.

<table><caption>**Table 2**
Exponential Scale Lengths of M1, M1const, M2, and M3 for the Stellar and Gaseous Surface Mass Distributions</caption><thead><tr><td>Model Name</td><td>r<sub>e,stars</sub>(kpc)</td><td>r<sub>e,gas</sub>(kpc)</td></tr></thead><tbody><tr><td>M1</td><td>1.61 ± 0.04</td><td>5.75 ± 0.06</td></tr><tr><td>M1const</td><td>1.19 ± 0.01</td><td>6.44 ± 0.07</td></tr><tr><td>M2</td><td>1.35 ± 0.01</td><td>7.24 ± 0.16</td></tr><tr><td>M3</td><td>1.98 ± 0.02</td><td>7.79 ± 0.10</td></tr></tbody></table>

Hence, changing the initial rotation velocity or the rotation law to a constant one does not alter the fact that a flat rotation curve and an exponentially decreasing disk are obtained. However, the initial collapse is different, which changes the morphology (especially the radial extent/scale length) and accordingly the density distribution. This suggests that low surface brightness galaxies result from collapsing gas clouds with high specific angular momentum. As a side remark, note that, initially, we have constant temperature and density, and hence constant pressure. Different initial temperatures and/or different profiles might also lead to slightly different collapse patterns.
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Comparison between rotation curves of MONDian and Newtonian simulations for two different initial conditions. From top left to bottom right: M1, M1N, M2, and M2N. See the text for further details.</figcaption></figure>

### 4.3. Comparison between Newtonian and MONDian Models

Due to the fact that exponential stellar surface mass density profiles occur in all the MONDian galaxy models calculated here, and for the purpose of testing the correct behavior of the code and the role of the gravitational law, M1 and M2 were simulated again, while only using the Newtonian Poisson solver, but still without dark matter.

Figures 6 and 7 show, respectively, the surface mass density profiles and the rotation curves of the equivalent Newtonian and MONDian models. First of all, the rotation curves of the Newtonian models show no flatness but a Keplerian decrease, which is expected and demonstrates that the code works properly in this regard. On the other hand, the majority of the rotation curves of the MONDian models are flat, as seen before. Also, the measured rotation curves of the Newtonian models are very close to the theoretical spherical-equivalent curves, $v_{c,N}$ (purple boxes), which are calculated by using the mass inside the respective radius, $M(<r)$. Comparing $v_{c,N}$ with the MONDian curves reveals that the center of M2 lies within the Newtonian acceleration regime of our interpolating function, because the curves overlap in this region. Contrary to that, the accelerations in M1 do not reach values well above

$10^{-10}$ m s$^{-2}$, as the curves never meet. This overlapping region depends on the depth of the potential and extends further in radius the more massive a model is. Of course, the gap between the MONDian rotation curve and the theoretical Newtonian one increases with radius, because the rotation curves become asymptotically flat for M1 and M2, while the theoretical Newtonian curves show a Keplerian decrease.

Qualitatively the comparison of the surface mass densities shows an encouraging result, as the stellar surface mass density profiles of the Newtonian models do not show an exponential profile but are rather curved. However, $\Sigma_{gas}$ shows an exponentially decreasing profile similar to the simulations done with MOND, but with a slower radial falloff, therefore leading to different exponential scale lengths, $r_{e,gas,M1} = 5.37 \pm 0.10$ kpc, $r_{e,gas,M1N} = 9.62 \pm 0.11$ kpc and $r_{e,gas,M2} = 7.24 \pm 0.16$ kpc, $r_{e,gas,M2N} = 8.55 \pm 0.13$ kpc.

Although the forces for the two gravitation laws are different, when using the same mass content, it is still noteworthy that $\Sigma_{stars}$ is decreasing exponentially for all models simulated with the MONDian Poisson solver, while the Newtonian models do not show this behavior. This indicates that the natural appearance of exponentially decreasing stellar surface mass density profiles seems to be a feature of
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Surface mass density profiles of MONDian and Newtonian simulations for two different initial conditions. From top left to bottom right: M1, M1N, M2, and M2N. See the text for further details.</figcaption></figure>

the MONDian framework and not a feature of the initial conditions.

### 4.4. Initial Size/Mass Limit

As seen in Figures 4 and 5, the slope, $\eta$, of the rotation law used here changes the size of the computed galaxy noticeably, in contrast to the subtle differences that occur by using more complex baryonic physics. Furthermore, the models are sensitive not only to the slope of the initial rotation law but also to $r_{init}$ and $M_{init}$. The models M4 and M4sn show significantly different stellar surface mass density profiles and formation histories compared to all other simulations.

Again, the initial collapse is very important for the evolution of these models. M4 and M4sn have nearly five times more mass than even the intermediately massive models, but they are collapsing into a comparably sized region. Therefore, the collapse itself is more violent and dense, and more stars are formed that have an elongated orbit. Also, many more stars are outside of the dense disk, which leads to a more populated and more extended lower-density stellar component. Furthermore, like for M3, several dense clumps form very shortly before the overall collapse happens, but in addition to that, in M4 one dense spiral-arm-like structure drifts away from the

collapse and small satellites form out of it, which collide with the galaxy at about 6 Gyr after the start of the simulation. Supernovae seem to suppress the formation of these elongated dense clouds and also weaken the collapse overall, but still a small satellite forms during the collapse and merges very shortly afterward with the main galaxy in M4sn. Considering the stability criterion for disks in MOND (Banik et al. 2018a),

$$ \begin{aligned} c_s &\ge \frac{\pi G \nu \Sigma \left(1 + \frac{K_0}{2}\right)}{\Omega_{r}} &\Leftrightarrow 1 \\ &\le \frac{c_s \Omega_{r}}{\pi G \nu \Sigma \left(1 + \frac{K_0}{2}\right)} &\equiv Q_{\text{MOND}}, \end{aligned} \quad (13) $$

where $c_s$ is the sound speed of the gas, $\nu \propto \sqrt{a_0/|g_N|}$ is the transition function, $\Sigma$ is the surface mass density, $K_0 = d \ln[\nu(y)]/d \ln(y)$ at $y = g_{N0}/a_0$, where the subscript 0 for $g_{N0}$ and $K_0$ indicates that these are unperturbed values (see Banik et al. 2018a), and $\Omega_r^2 = -3g_r/r - \partial g_r/\partial r$, the formation of clumps for the more massive models, in contrast to M1, for example, becomes clearer. Equation (13) can be approximated using the analytical proportionalities of $\nu$, $\Sigma$ and $\Omega_r$, as
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Rotation curves of both models with $r_{\text{init}} = 50$ kpc. Left panel: M4; right panel: M4sn. See the text for further details.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Stellar surface mass density profiles of both models with $r_{\text{init}} = 50$ kpc. Left panel: M4; right panel: M4sn. The stellar surface mass density distribution does not follow a single exponential profile, as the formation of these models is dominated by mergers of smaller satellites during and shortly after the initial collapse. The surface density profile may be approximated by more than one exponential profile (not attempted here though). See the text for further details.</figcaption></figure>

the rest is constant. With $v \propto r/\sqrt{M}$, $\Sigma \propto M/r^2$, and $\Omega_r \propto v_{\text{rot}}/r \propto \sqrt[4]{M}/r$, $Q_{\text{MOND}}$ scales as $1/\sqrt[4]{M}$, which shows that there is a mass limit above which multiple fragments are likely to form for a given sound speed, $c_s$. In reality, this is more complicated owing to the explicit form of $v$, $\Sigma$ and $\Omega_r$, which, for example, can then lead to central instabilities and the formation of a bulge, while the outer disk is stable.

The comparison between M4 and M4sn is shown in Figures 8 and 9. Although the initial collapse in M4 and M4sn gives rise to the formation of small satellites and a very extended stellar density profile, it does not affect the rotation curve. As for every other model simulated with the MONDian Poisson solver, the rotation curve flattens in the outer part of

the galaxies. The rotation curve of M4 is nearly identical to the “pristine” curve except in the very center, while deviations are also apparent in the transition region for M4sn.

The major difference of M4/M4sn from all other models is the surface mass density distribution. The decrease for $\Sigma_{\text{stars}}$ is shallower than for the less massive models, because the low-density stellar component is very populated and extended to a degree where $\Sigma_{\text{stars}}$ nearly becomes flat at about 50 kpc. Furthermore, $\Sigma_{\text{stars}}$ is higher than $\Sigma_{\text{gas}}$ within a radius of 50 kpc, which can be explained by the low gas fractions of the two models, $f_{\text{g,M4}} = 0.04$ and $f_{\text{g,M4sn}} = 0.05$. The central galaxy (not the satellites) forms partially from mergers of the smaller satellites during and shortly after the initial collapse, and the resulting stellar surface mass density profile ceases to
<figure><img src="image_11.png" /><figcaption>**Figure 10.** The $r_{\text{eff}}-M_s$ relation (Equation (14)) compared with our models and the MW. Our models are colored according to the key, with darker color indicating sn and lighter color Newtonian models (e.g., M1 = green, M1sn = dark green, and M1N = light green). The data for the MW are taken from Banik & Zhao (2018; mass) and Mor et al. (2018; exp. scale length). The upper distribution of points shows the relation between stellar mass and effective radius, where the latter is based on the fit to the exponential stellar disks, $r_{\text{eff}} = 1.678r_e$. The lower distribution shows the same for the stellar half-mass radius, $r_h$. Missing points correspond to nonexistent exponential stellar surface mass density profiles, e.g., M4 (plum).</figcaption></figure>

be a simple exponential profile. This suggests that galaxies, which do not have a single exponential stellar surface mass density profile, may have been partially formed by baryonic mergers.

Increasing the initial mass further will lead to more satellites and an even more violent collapse, so there seems to be a limit for the formation of isolated galaxies from a single rotating gas cloud with insignificant later accretion of additional gas. Moreover, the formation of dense clumps prior to the collapse can be enhanced by increasing the initial mass and also by increasing $\eta$. If both are done, it is possible to produce not only satellites but also groups of galaxies, which will be shown in a forthcoming contribution (N. Wittenburg et al., in preparation).

To test whether our models that show an exponentially decreasing stellar surface mass density profile are compatible with observations, we compared our models with the size-mass relation (Lange et al. 2015),

$$r_{\text{eff}} = 0.13(M_s)^{0.14} \left( 1.0 + \frac{M_s}{14.03 \times 10^{10} M_{\odot}} \right)^{0.77}, \quad (14)$$

where $r_{\text{eff}}$ is the effective radius, with $r_{\text{eff}} = 1.678r_e$ for exponential disks (where $r_e$ is the scale length), and $M_s$ is the stellar mass.

Figure 10 depicts the $r_{\text{eff}}-M_s$ relation from Lange et al. (2015; black line), with the uncertainty shown as the gray area; our models are color-coded as indicated by the key (the sn models are shown with the respective darker color and the Newtonian with lighter color, e.g., M1 = green, M1sn = dark green, and M1N = light green), and the MW is depicted with a mass estimate from Banik & Zhao (2018) and the scale length measurement from Mor et al. (2018).

Indeed, most of our models follow the observed relation within the uncertainties, when using the effective radius, which is based on the exponential fit to the stellar disk. When using the stellar half-mass radius, $r_h$, the distribution looks similar, but it is offset to too low radii. There are two important results that emerge from these different distributions. First, the surface density profiles of the stellar disks of our models are compatible with observed ones. Second, because $r_{\text{eff}}$ is defined as the half-mass radius of a purely exponential disk, the discrepancy between $r_{\text{eff}}$ and $r_h$ should be much smaller if at all evident in disk-dominated galaxies. Because $r_{\text{eff}} > r_h$ here, our models have a somewhat massive central region, which may stem from the formation process, as the galaxies do not grow, but form in practically one monolithic collapse. These compact bulges may also be related to the observed red nuggets (de la Rosa et al. 2016). This feature of our models will also be important in Section 4.6, where we show that they slightly lack stellar angular momentum.

Important to note here is that it is a priori not clear that MONDian models should follow the size-mass relation, but this appears to be the natural result in MOND.

### 4.5. Star Formation History and Gas Depletion Time

Although the structure and morphology of the present models are very comparable to those of real galaxy disks (see also the comparison between observational scaling relations and the data of the models in Section 5), the simplicity of the initial conditions, i.e., a collapsing gas sphere with insignificant further accretion of gas, gives rise to a different star formation history (SFH) compared to observations of real galaxies.

The left panel of Figures 11 and 12 shows the SFH of M1 and M1N, respectively. The SFR is calculated by separating all stellar particles in bins of 100 Myr according to their age, summing up the stellar mass in every time bin and dividing this by the length of the time bin.

The SFR increases sharply at the time of the collapse of the sphere until the maximum is reached (the halo forms during this time), and then it decreases exponentially as shown by the fit. The function

$$h(t) = e \times \exp(-t/t_e), \quad (15)$$

with $e$ and $f$ being the fit parameters, is used to fit the first (green line, 1) and second (blue line, 2) part of the decline for M1, and for M1N the second part is fitted by a constant, $c$. This results in $t_{e,M1,1} = 0.55 \pm 0.02$ Gyr, $e_{M1,1} = 14.67 \pm 0.56 M_{\odot}yr^{-1}$, $t_{e,M1,2} = 2.53 \pm 0.06$ Gyr, $e_{M1,2} = 0.47 \pm 0.02 M_{\odot}yr^{-1}$ for M1 and in $t_{e,M1N,1} = 0.57 \pm 0.02$ Gyr, $e_{M1N,1} = 33.46 \pm 4.10 M_{\odot}yr^{-1}$, $c_{M1N} = 0.02 \pm 4 \times 10^{-4} M_{\odot}yr^{-1}$ for M1N. The fit parameters for every model are shown in Table 5.

The right panel of Figures 11 and 12 shows the gas depletion time of M1 and M1N at every output time,

$$t_{\text{gas,d}}(t) = \frac{M_{\text{gas}}(t)}{\text{SFR}(t)}. \quad (16)$$

Here $t_{\text{gas,d}}(t)$ corresponds to the time it would take to convert the remaining gas mass, $M_{\text{gas}}$, at time $t$ with the SFR at that time into stellar particles.

The evolution of the gas depletion time is very similar to the SFH (just inverted), so the initial collapse corresponds to the minimum at the beginning of the curve, after that $t_{\text{gas,d}}$ rises exponentially until most of the gas has been consumed, and
<figure><img src="image_12.png" /><figcaption>Figure 11. SFH and evolution of the gas depletion time for M1. The red circles show the model data; the green line shows the exponential fit for the first part of the decaying SFR/increasing gas depletion time directly after the initial collapse. The blue line shows the exponential fit for the shallower part after most of the gas was converted into stellar particles during and shortly after the collapse.</figcaption></figure>

then the increase becomes shallower or, in the case of M1N, $t_{gas,d}$ becomes constant. The exponential parts (directly after the collapse: green line; shallower/constant part: blue line) are again fitted by a simple exponential function, $g(t)$,

$$g(t) = i \times \exp(t/j), \qquad (17)$$

where i and j are the fit parameters. The constant part of M1N is also again fitted by a constant, c. The results are $j_{M1,1} = 0.85 \pm 0.02$ Gyr, $i_{M1,1} = 0.36 \pm 0.02$ Gyr, $j_{M1,2} = 3.17 \pm 0.08$ Gyr, $i_{M1,2} = 4.24 \pm 0.21$ Gyr for M1 and $j_{M1N,1} = 0.84 \pm 0.02$ Gyr, $i_{M1N,1} = 0.32 \pm 0.04$ Gyr, $c_{M1N} = 117.92 \pm 2.27$ Gyr for M1N. The fit parameters for every model are shown in Table 6. Again, indices 1 and 2 correspond to the steeper (green) and shallower/constant (blue) part, respectively (the SFH and the evolution of $t_{gas,d}$ for every model can be seen in the Appendix of the arXiv version of this paper).

Observed galaxies show a mildly decreasing SFR with cosmic time (Speagle et al. 2014), such that the present-day gas depletion times are independent of galaxy mass and about 2.8 Gyr (Pflamm-Altenburg & Kroupa 2009). Galaxies therefore need to accrete gas accordingly. This is different for the models computed here, because there is no reservoir of gas surrounding the model galaxies such that accretion of gas is insignificant, and so the gas density in the galaxy model decreases continuously. Therefore, the SFR decreases and the gas depletion time increases. In the case of M1N the gas spreads further, due to the weaker potential compared to MOND, so the peak of the SFH is smaller but also more extended in time. After approximately 5 Gyr, the SFR and gas depletion time become constant, in contrast to the behavior of M1, because the gas density is higher for M1N compared to M1 owing to the initial collapse not consuming as much gas. This also leads to different gas fractions as can be seen in Figure 18.

Figure 13 shows the Kennicutt-Schmidt diagram (Kennicutt 1998) for all simulations. We adopted the layout of Figure 15 of Bigiel et al. (2008), as they further investigated the relation between the SFR density, $\Sigma_{SFR}$, and the gas surface mass density, $\Sigma_{gas}$, and to compare their results with our data. $\Sigma_{gas}$ is calculated by averaging the gas surface mass density over the stellar disk, and $\Sigma_{SFR}$ shows the SFR from, e.g., Figure 11 divided by the area of the stellar disk.

The diagonal dotted lines in Figure 13 correspond to star formation efficiencies of 1%, 10%, and 100% as indicated in the plot, meaning that these lines show the SFR needed to turn 1%, 10%, or 100% of the remaining gas mass into stars within $10^8$ yr. The two vertical black and dotted lines show the limits Bigiel et al. (2008) established to divide three distinctly different regimes in the $\Sigma_{SFR}-\Sigma_{gas}$ plane. The left vertical line corresponds to the saturation of H I at a surface density of $\Sigma_{gas} \approx 9 M_{\odot} \text{ pc}^{-2}$, while the second line at $\Sigma_{gas} \approx 200 M_{\odot} \text{ pc}^{-2}$ is the transition from “normal” galactic SFRs to starburst galaxies (see Bigiel et al. 2008). The green line shows the best fit, Equation (4), of Kennicutt (1998),

$$\Sigma_{SFR} = (2.5 \pm 0.7) \times 10^{-4} \times \left( \frac{\Sigma_{gas}}{1 M_{\odot} \text{ pc}^{-2}} \right)^{1.4 \pm 0.15} M_{\odot} \text{ yr}^{-1} \text{ kpc}^{-2}. \qquad (18)$$

The black diamonds are the two Newtonian models, M1N and M2N, while the colored crosses depict all MONDian simulations. The evolution of each model is plotted in steps of 100 Myr, while the age of every point is shown by its color. Similar to Figure 15 of Bigiel et al. (2008), our simulations show a shallower increase of $\Sigma_{SFR}$ in the central regime compared to the points with $\Sigma_{gas} \le 9 M_{\odot} \text{ pc}^{-2}$. Due to the fact that very few of our models reach gas surface mass densities well within the starburst regime, we cannot make any reasonable statement for MONDian models in this regime. However, we can also confirm that the fit from
<figure><img src="image_13.png" /><figcaption>Figure 12. Same as Figure 11, but for M1N.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 13. KS diagram for all 15 models shown as a time evolution, while the age is indicated by color. The Newtonian models are emphasized as black diamonds, though they show the same trend as all other models. The green line corresponds to the best fit to observational data by Kennicutt (1998), while the vertical dotted lines correspond to the dividing lines between different regimes according to Bigiel et al. (2008; see the text). The diagonal dotted lines show star formation efficiencies of 1%, 10%, and 100%, respectively, meaning that 1%, 10%, and 100% of the remaining gas mass is consumed within 10$^8$ yr.</figcaption></figure>

Kennicutt (1998) is offset to higher SFRs, because the author used starburst and nonstarburst galaxies in one fit, as Bigiel et al. (2008) stated.

In general, our models are qualitatively similar to the observational data from Bigiel et al. (2008), showing two distinct regimes with a transition close to the H I saturation line. On the other hand, the rapid evolution of our models is also evident, as the slope in the intermediate region, $\Sigma_{gas} \ge 9 M_{\odot} pc^{-2}$, is steeper than observed. Important to note is that also the Newtonian models show a similar but slightly

shallower behavior than the MONDian models, which suggests that the critical law here is the description of star formation itself, rather than gravity. Additionally, we stress that we show the time evolution of our models, while Bigiel et al. (2008) plot the SFR of nearby galaxies.

To investigate the evolution of the SFRs of our models from another angle and to test whether they lie near the observed main sequence of galaxies (Speagle et al. 2014), we also plotted the SFR versus stellar mass ($M_s$) dependency. We use the function $SFR_{ms}^7$

$$ \log(SFR_{ms}(M_s)) = \log\left(\frac{a}{1 M_{\odot} yr^{-1}}\right) + b \times \log\left(\frac{M_s}{M_{\odot}}\right), \quad (19) $$

with $a$ and $b$ being the fit parameters, to fit our data.

The left panel of Figure 14 shows our models directly at their peak SFR (upper part of the plot) and after ≈5 Gyr compared to the observed main sequence of galaxies of Speagle et al. (2014) at $t = 0.5$ Gyr (red area) and $t = 5.0$ Gyr (violet area). We only fitted the upper part, resulting in $a = (3.66 \pm 0.35) \times 10^{-9} M_{\odot} yr^{-1}$ and $b = 1.041 \pm 0.004$, while the distribution of the models after ≈5 Gyr suggests that the scatter of the correlation increases with time. At $t = 0.5$ Gyr the fit to our models suggests a steeper slope than Equation (20), but the models lie within the uncertainties of the observed relation. After $t = 5.0$ Gyr, the picture is different, as the scatter is strongly increased and the SFR of our models is offset by approximately two orders of magnitude compared to the observed relation. The right panel of Figure 14 shows the evolution of M1, M3, and M4 with time in the SFR-$M_s$ plane and the evolution of the linear fit to all models at various times. The thin lines are the fits from the peak SFR at ≈500 Myr to ≈1.5 Gyr (the first thick, dotted line) in 100 Myr steps. From that time onward, the step size is increased to 1 Gyr, and the last line shows the fit at 9.5 Gyr. The slopes, $b$, of the fits range from 1 to 0.25, while the uncertainty from the fit increases with

<footnote n="7">7 Note that $\log \equiv \log_{10}$ throughout.</footnote>
<figure><img src="image_15.png" /><figcaption>Figure 14. Left panel: SFR–$M_s$ plane, with all models at their respective peak SFR at the top and the distribution of the models 5 Gyr afterward at the bottom. The black line shows the linear fit to the peak SFR data, resulting in $a = (3.66 \pm 0.35) \times 10^{-9} M_{\odot}yr^{-1}$ and $b = 1.041 \pm 0.004$. The colors are identical to Figure 10, with the addition of M4 = plum, M4sn = dark plum, and the Newtonian models indicated by light colors (e.g., M2 = blue and M2N = light blue). The red and violet areas show the observed main sequence of galaxies of Speagle et al. (2014) at $t = 0.5$ and 5.0 Gyr, respectively (Equation (20)). Right panel: evolution of the linear fit to all models shown as the dashed lines at various times (see the text) and the complete time evolution of M1 (green), M3 (yellow), and M4 (plum).</figcaption></figure>

time up to 85% of the fitted value. We note that the observed main sequence of galaxies follows

$$ \log(\text{SFR}(M_s, t)) = (0.84 \pm 0.02 - 0.026 \pm 0.003 \times t) \times \log M_s - (6.51 \pm 0.24 - 0.11 \pm 0.03 \times t), \quad (20) $$

with $t$ being the age of the universe (Speagle et al. 2014), such that the model result appears to be promising in terms of possibly helping to understand the origin and evolution of the main sequence. Clearly much more work is required, in particular to address the issue of continued gas accretion.

Different models react differently to the depletion of gas, due to the insignificant accretion, as indicated by the left panel of Figure 14. For the same $\eta$ and initial mass, there is a difference between models of approximately one order of magnitude in SFR after ≈5 Gyr (points at the left end). For example, the highest SFR is evident in M1113, which is shown by the dark-turquoise diamond. This is expected, as the density threshold for a star-forming event to take place needs to be met in a smaller volume compared to all other simulations. The comparison with the Speagle et al. (2014) relations underlines the fact that the SFR decreases more rapidly than observed, which we would expect owing to the lack of further accretion. The most prominent effect of the absence of accretion can be seen in the right panel of Figure 14. The evolutionary track of M4 in the SFR–$M_s$ plane becomes nearly vertical ≈300 Myr after its peak SFR, while the track of M1 shows that behavior more than ≈1 Gyr later. This is a direct consequence of the main sequence of galaxies, because the more massive a galaxy is, the more gas it has to accrete to sustain its SFR.

## 4.6. Stellar Angular Momentum and Disk versus Halo Populations, Thin and Thick Disks

After analyzing the SFH and showing that there might be a mass and initial cloud angular momentum limit to the single collapsing gas sphere approach, the question remains how the

angular momentum distribution of the stars is shaped by the simple initial conditions used here. We note that a series of simulations of the same model (M1) with different resolution (M1111 and M1113) yield the same total angular momentum in the stellar particles within 5%.

Figure 15 shows the ratio of the specific stellar angular momentum, $j_{star}$, and its z-component, $j_{z,star}$, for all stars (left panel) and only for old (>8 Gyr) halo/bulge stars (right panel) against the spherical radius, $r_{sph}$, for M1 (the plots for every other model can be seen in the Appendix of the arXiv version of this paper). The age of every stellar particle is shown via color from ≈10 Gyr (bright yellow) to a few megayears (black) as indicated by the colored age scale. In order to obtain the contribution of Population II particles (older than 8 Gyr) to the mass of the bulge and halo, the disk particles were sorted out. These are assumed to have $j_{z,star}/j_{star} \ge 0.98$ and an age smaller than 8 Gyr. The green line in the right panel corresponds to the dividing line between halo and bulge, which is defined here as the radius at which the smallest ratio of densities in neighboring radial bins occurs, i.e., $n_{star}(i)/n_{star}(i-1)$, where i is the number of the radial bin, is computed for all bins, and the radius of bin i with the smallest ratio within a reasonable distance to the center, $r_{sph} < 10$ kpc, is chosen to be the radius of the bulge, $r_{bulge}$.

The left panel of Figure 15 reveals that the youngest stars are concentrated in the center of the galaxy and most of them move within the disk, as the contribution of $j_{z,star}$ to $j_{star}$ is greater than or equal to 0.98. The rest of the young stars can be found within the bulge, but not within the halo. This is very different from observed late-type galaxies, because especially the most massive ones show a higher SFR in the disk compared to the center. On the other hand, the discrepancy between the models and real galaxies is expected, as the models computed here do not accrete gas during their evolution, and hence the center is the only region where star formation continues after the galaxy has formed.
<figure><img src="image_16.png" /><figcaption>**Figure 15.** Radial distribution of the ratio of the specific stellar angular momentum and its z-component for M1 after 10 Gyr, with $r_{\text{sph}}$ being the spherical radius, $j_{\text{star}}$ the specific angular momentum, and $j_{\text{star}}$ its z-component. The age of every star is shown via color from 10 Gyr (bright yellow) to a few megayears (black). Left panel: all stellar particles. Right panel: only bulge/halo stars that are older than 8 Gyr (see the text for details).</figcaption></figure>

<table><caption>**Table 3**

Stellar Masses and Mass Fractions of the Bulge and the Halo and also the Total Stellar and Gas Mass of Every Model after 10 Gyr</caption><thead><tr><th>Model Name</th><th><i>M</i><sub>tot,star</sub>(10<sup>9</sup><i>M</i><sub>☉</sub>)</th><th><i>M</i><sub>tot,gas</sub>(10<sup>9</sup><i>M</i><sub>☉</sub>)</th><th><i>M</i><sub>B,star</sub>(10<sup>9</sup><i>M</i><sub>☉</sub>)</th><th><i>M</i><sub>B,star</sub>/<i>M</i><sub>tot,star</sub></th><th><i>M</i><sub>B,star</sub>(10<sup>9</sup><i>M</i><sub>☉</sub>)</th><th><i>M</i><sub>B,star</sub>/<i>M</i><sub>tot,star</sub></th></tr></thead><tbody><tr><td>M1</td><td>5.21</td><td>1.19</td><td>2.00</td><td>0.38</td><td>0.98</td><td>0.19</td></tr><tr><td>M1sn</td><td>4.33</td><td>1.69</td><td>1.75</td><td>0.40</td><td>1.07</td><td>0.25</td></tr><tr><td>M1N</td><td>3.72</td><td>2.66</td><td>1.34</td><td>0.36</td><td>0.75</td><td>0.20</td></tr><tr><td>M1const</td><td>5.88</td><td>0.52</td><td>3.19</td><td>0.54</td><td>0.60</td><td>0.10</td></tr><tr><td>M1Zpoor</td><td>5.12</td><td>1.28</td><td>2.11</td><td>0.41</td><td>0.86</td><td>0.17</td></tr><tr><td>M1Zpoorsn</td><td>4.83</td><td>1.44</td><td>1.85</td><td>0.38</td><td>0.89</td><td>0.18</td></tr><tr><td>M1111</td><td>4.90</td><td>1.50</td><td>2.27</td><td>0.46</td><td>0.68</td><td>0.14</td></tr><tr><td>M1113</td><td>5.36</td><td>1.04</td><td>2.53</td><td>0.47</td><td>1.16</td><td>0.22</td></tr><tr><td>M2</td><td>20.48</td><td>1.12</td><td>10.75</td><td>0.52</td><td>2.03</td><td>0.10</td></tr><tr><td>M2sn</td><td>19.47</td><td>1.82</td><td>10.25</td><td>0.53</td><td>2.26</td><td>0.12</td></tr><tr><td>M2N</td><td>20.39</td><td>1.19</td><td>9.79</td><td>0.48</td><td>3.55</td><td>0.17</td></tr><tr><td>M3</td><td>19.23</td><td>2.37</td><td>7.47</td><td>0.39</td><td>3.44</td><td>0.18</td></tr><tr><td>M3sn</td><td>18.38</td><td>3.03</td><td>7.27</td><td>0.40</td><td>3.35</td><td>0.18</td></tr><tr><td>M4</td><td>95.62</td><td>4.07</td><td>36.49</td><td>0.38</td><td>20.78</td><td>0.22</td></tr><tr><td>M4sn</td><td>94.90</td><td>4.91</td><td>35.51</td><td>0.37</td><td>21.73</td><td>0.23</td></tr></tbody></table>

**Note.** Note that the mass fraction of the bulge decreases the weaker the initial collapse is.

The right panel of Figure 15 shows the radial distribution of $j_{\text{star}}/j_{\text{star}}$ for all stars older than 8 Gyr. The bulge and the halo are separated by the green line, although the bulge is also visible by eye, as it is the region where a significant number of stars are not rotating around the z-axis. Comparing the left panel with the right one displays that there are many stars (especially older ones) that do not move within the disk, which already hints that there is a slight lack of stellar angular momentum. This is further emphasized by Table 3, where the masses and mass fractions of the halo and bulge for all models are shown. The result is that more than 50% of the total stellar mass is not in the disk, which will be put into perspective in Figure 17.

Figure 16 shows the distribution of mass along the z-axis, $\Sigma_z$, versus the vertical distance (height not thickness), z, for M1. In contrast to the distribution of mass outside of the disk, $\Sigma_z$ is comparable to observations, as the distribution can be

fitted by two simple exponential functions, which would correspond to the profile of a thin and thick disk. For example, the MW also shows a vertical stellar mass distribution with a thin and a thick disk. The fits are obtained with a similar function to that for the surface mass profiles,

$$ \Sigma_z(z) = n_z \times \exp(-z/z_e), \qquad (21) $$

resulting in $n_{z,\text{thin,M1}}/2 \times 10^7 M_{\odot} \text{ pc}^{-1} = 0.316 \pm 0.009$, $z_{e,\text{thin,M1}} = 0.305 \pm 0.007 \text{ kpc}$ and $n_{z,\text{thick,M1}}/2 \times 10^7 M_{\odot} \text{ pc}^{-1} = 0.096 \pm 0.004$, $z_{e,\text{thick,M1}} = 0.606 \pm 0.011 \text{ kpc}$.

There are a few differences between the models computed here. Models based on M1 show a very smooth transition between thick disk and halo, where the halo can only be recognized in Figure 16 by the increasing scatter around the blue line of the thick-disk fit. All other MONDian models (i.e., based on M2, M3, and M4) show a shallower halo profile. The
<figure><img src="image_17.png" /><figcaption>Figure 16. Mass distribution along the rotation axis (z-axis), Σz. The green line corresponds to the exponential fit for the thin disk, and the blue line corresponds to the thick disk, fitted in the ranges 0–1 kpc and 1–3 kpc, respectively (see the text for details).</figcaption></figure>

fit parameters for every model are shown in Table 7. All Newtonian models do not show a thick disk, but they decrease according to a power law beyond z ≈ 1 kpc,

$$ \Sigma_z(z) = az^{-b}, \qquad (22) $$

where a and b are the fit parameters for the Newtonian models for the power-law part (all vertical profiles can be seen in the Appendix of the arXiv version of this paper).

Figure 17 depicts the Fall relation, which compares the total specific stellar angular momentum of a galaxy, jstar, with its stellar mass, Mstar. Note that, for this plot, we do not compute the true specific stellar angular momentum from individual stellar particles, but rather estimate it exactly as for observations (Posti et al. 2018). The magenta area with the dashed line shows the best fit of the observational data from Posti et al. (2018), while the black line shows the best fit for the data of the models. The same functions as used by Posti et al. (2018) are here applied to calculate jstar and fit the data,

$$ j_{\text{star}}(<R) = \frac{\int_{0}^{R} dR'R'^{2}\Sigma_{\text{star}}(R')V_{\text{star,rot}}(R')}{\int_{0}^{R} dR'R'\Sigma_{\text{star}}(R')}, \qquad (23) $$

$$ \log j_{\text{star}} = \alpha[\log(M_{\text{star}}/M_{\odot}) - 11] + \beta. \qquad (24) $$

Equations (23) and (24) are Equations (1) and (4) of Posti et al. (2018), where Σstar is the stellar surface mass density, Vstar,rot is the stellar rotation curve of the galaxy, and α, β are the fit parameters. Posti et al. (2018) find their best fit with α = 0.55 ± 0.02 and β = 3.34 ± 0.03, while the best fit for the models here is α = 0.51 ± 0.04 and β = 3.05 ± 0.04. Figure 17, together with these results, shows that the simulations done in this work follow the Fall relation closely (α is compatible with observations), but they are slightly offset to lower values of jstar. Thus, the slope α of the specific stellar angular momentum is very similar to observed galaxies, but jstar is slightly too small. A large fraction of the initial angular

<figure><img src="image_18.png" /><figcaption>Figure 17. Relation between the specific stellar angular momentum, jstar, of all models with their respective stellar mass, Mstar (Fall relation, computed here as for observations). The magenta area with the dashed line shows the observed data from Posti et al. (2018), the symbols for the different models are shown in the legend, and the black line (Equation (24)) shows the best fit to the models. See the text for further details.</figcaption></figure>

momentum of the original cloud remains in the gas component. This is connected to the simplistic initial conditions, especially the major collapse at the beginning and the absence of accretion, resulting in mass-concentrated cores of the models (see again Figure 10 and the corresponding explanation), and hence jstar is offset from the observed value. Such compact spheroidal components may have been observed as red nuggets (de la Rosa et al. 2016). This may also be partially related to the detailed subgrid prescriptions used to form stars, which are fully applicable in the ΛCDM framework that they were calibrated on but are not necessarily applicable in MOND. Nevertheless, that the model and observed values of α are so similar is unexpected and suggests again that the formation of galaxies in Milgromian gravitation may contain nontrivial aspects of reality.

## 4.7. Overview of All Modeled Galaxies

Figure 18 shows an overview of the analysis of all models at the final time (10 Gyr after the start of the simulation). We also include additional simulations with lower metallicity, Z = 10⁻⁴ × Z☉, and higher and lower resolution (117.1875 and 468.75 pc) to confirm that these parameters have only a minor effect on our results. Panels (a) to (c) present properties of the SFH, panels (d) and (e) properties of the rotation curve, panels (j) and (k) properties of the surface mass density profiles, and panel (l) shows the gas fraction, fg,

$$ f_g = \frac{M_{\text{gas}}}{M_{\text{tot}}}, \qquad (25) $$

with Mgas being the gas mass and Mtot the total baryonic mass.
<figure><img src="image_19.png" /><figcaption>**Figure 18.** Properties of all galaxy models at 10 Gyr. The color for every point is the same as in Figure 10 (see Table 1 for the number and description of the models), but the Newtonian models are shown as circles instead of diamonds. Panels (a)–(c) are linked to the SFH: (a) exponential star formation decay time (see Equation (15)), (b) time of the highest SFR, (c) value of the highest SFR. Panels (d)–(f) show properties of the rotation curves: (d) flat rotation velocity (flat part of the rotation curve, Equation (11)), (e) highest rotation velocity, (f) distance from the center at which the highest $v_{\text{rot}}$ occurs. Panels (j) and (k) show the exponential scale lengths for the stellar and gas profiles, respectively (see Equation (12)), and panel (l) shows the gas fraction of every model (Equation (25)). Missing points imply that the model does not have the respective property, e.g., Newtonian models do not have an asymptotic rotation velocity.</figcaption></figure>

Panel (a) shows the exponential star formation decay time of every model, which is decreasing with increasing $M_{\text{init}}$, but increasing when using a steeper slope (larger $\eta$) of the initial rotation law. This emphasizes the explanations made earlier about how violent or strong the collapse is. As the initial mass increases, not only the peak height of the SFR increases but also the flanks of the peak become steeper, and therefore the exponential decay time decreases. On the other hand, if the initial rotation velocity is increased, the collapse is weakened and a high SFR is evident over a longer period of time, so the peak broadens and the exponential decay time increases. Panel (b) shows that for all MONDian models the collapse takes place at the same time, $t = 350$ Myr, regardless of $M_{\text{init}}$ or $\eta$,

while the two Newtonian models collapse later (M1N: $t = 950$ Myr; M2N: $t = 650$ Myr). Panel (c) shows the aforementioned trend that higher initial mass corresponds to a higher peak SFR, but also that supernovae and higher initial rotation velocities lead to lower peak SFRs and therefore weaken the initial collapse.

Panel (e) is related to the findings about the star formation histories, as it also shows that the peak rotation velocity increases with initial mass but decreases with the addition of supernovae and steeper slopes, $\eta$, of the rotation law used. This comes about because $v_{\text{rot}}$ is obtained from the acceleration at the respective radii and is therefore closely linked to the gravitational potential, which depends on the mass distribution.
The mass distribution is shaped by the initial collapse and rotation velocity as seen before. Panel (f) shows that the highest rotation velocity is reached close to the center of the models, and panel (d) shows that regardless of the complexity of the baryonic physics or $\eta$, the flat rotation velocity is nearly equal for the same initial mass, which is an alternative way of describing the MASR (BTFR).

The exponential scale lengths shown in panels (j) and (k) for stars and gas, respectively, demonstrate that changes of the initial mass or rotation velocity and the addition of supernovae affect the stellar and gaseous distribution differently. While the stellar surface mass density distribution steepens with increasing mass and becomes shallower with the addition of supernovae and higher $\eta$, the steepest gas distributions, i.e., with the smallest radial exponential scale lengths, for models with the same $M_{init}$ are evident for the ones with supernovae enabled, except for M1Zpoorsn. Supernovae seem to push considerable amounts of gas outside of the stellar disk such that $\Sigma_{gas}$ is larger directly beyond the stellar disk for the supernova models compared to the simple ones. However, outside of the gaseous disk, both the sn and the simple feedback models show comparable surface mass densities, so the slope of the supernova models is steeper and therefore $r_{e, gas}$ is smaller. Moreover, a higher initial rotation velocity at the same $M_{init}$ does not change the slope of the gas distribution significantly for the supernova models (nos. 10 and 13), but it increases the exponential scale length for model M3 (no. 12) compared to M2 (no. 9). The shallowest surface mass density profiles for the gaseous component are produced by the Newtonian models, but this is expected, because the gravitational force is weaker in these models, so overall the mass distribution is more extended.

Panel (l) depicts the gas fraction of every model, and this again shows the expected behavior already seen in panels (c) and (e). The gas fraction decreases with higher initial mass, because the collapse is more violent and more gas is turned into stars, and it increases again with higher initial rotation velocity and the addition of supernovae, as these weaken the initial collapse. Real galaxies are well known to have an increasing gas fraction with decreasing baryonic mass, so this aspect is qualitatively reproduced by the present models.

We also want to emphasize the differences between the Newtonian and MONDian models, as the panels closely related to the dynamics of the system (panels (d), (e), (j), and (k)) show a significant difference between the two. In the Newtonian models there is no flat rotation curve and the peak rotation velocity is smaller than that of the corresponding MOND model, the stellar surface mass density profile is not described by an exponential profile, while the profile for the gas is significantly shallower. Additionally, star formation starts later and the evolution of the SFR is also shallower compared to the MONDian models.

Although M1const is simulated with a different (constant) rotation law, it underlines the aforementioned explanations about the initial collapse, because the major difference between the evolution of M1 and M1const is the more violent collapse for M1const. As can be seen in Figure 18, M1const shows the expected behavior compared to M1 in every panel that is related to the initial collapse.

At this point it is important to reiterate what has been shown for the non-Newtonian models, before proceeding to the next section:

1. All stellar surface mass density profiles decrease exponentially within the respective stellar disk (except M4/M4sn).

2. Changing the initial rotation velocity does not alter the resulting galaxy regarding the formation of a dense disk and an exponentially decreasing disk.

3. Adding more complex baryonic physics also does not change the fact that an exponentially decreasing disk and the same shape for the SFH come about. However, star formation is slightly suppressed, because supernovae push gas out of dense regions.

More complex baryonic physics therefore does not yield significantly different results for the same combination of $M_{init}$, $\eta$, and $r_{init}$.

4. All vertical mass density profiles show a thin and thick disk, while distinguishable disk, halo, and bulge components within the stellar distribution are evident.

5. The models slightly lack stellar angular momentum, which is linked to the initial conditions and the absence of accretion and thus the formation of a massive Population II halo and bulge. This may also be related to the detailed prescriptions used to form stars, which are fully applicable in the standard cosmological framework they were calibrated on, but may not necessarily be applicable in MOND. There is substantial observational evidence that galactic disks can grow around bulges (i.e., early-formed elliptical galaxies or red nuggets) through accretion (de la Rosa et al. 2016; Mancini et al. 2019).

These simulations with very simple initial conditions seem to match major properties of real galaxies. Galaxies with a stellar halo and nearly completely flat rotation curves form, and exponentially decreasing disks also seem to appear naturally within MOND.

## 5. The MASR and MDAR/RAR

After analyzing the differences between simple and more complex baryonic physics and different initial rotation laws and initial masses, we finally study, as a consistency check, whether the models comply with the MASR and MDAR/RAR, which arise from analytical MOND formulations. Of course, it is well known that these two relations should always be obeyed in MOND, as they are natural laws in this paradigm. But it is reassuring to see that our numerical models do comply with them, since they were initialized as pure gas clouds. The MDAR/RAR combines several scaling relations, including the MASR (BTFR), and has, for real galaxies, very little to no intrinsic scatter (Lelli et al. 2017). It relates the observed acceleration, $g_{obs}$, which is obtained from the rotation curve, to the baryonic acceleration, $g_{bar}$, which is the gravitational acceleration coming from the baryonic potential using the Newtonian gravitational law.

To compare the simulations with observations, the same fit function was used, and the best fit with uncertainties of the observational data is plotted together with the data from all models,

$$g_{\text{obs}} = \mathcal{F}(g_{\text{bar}}) = \frac{g_{\text{bar}}}{1 - e^{-\sqrt{g_{\text{bar}}/g_{\text{\textit{f}}}}}}, \quad (26)$$

Equation (26) here being the same as Equation (11) in Lelli et al. (2017), where only one free parameter, $g_{\text{\textit{f}}}$, is available to fit the data. This parameter is in MOND equivalent to $a_0$, and
<figure><img src="image_20.png" /><figcaption>Figure 19. Comparison of the simulations with observed galaxy scaling relations. Panel (a) shows the radial acceleration relation of the data obtained directly from the accelerations in the simulations and the observational data from Lelli et al. (2017; red area), the best fit to the data from all models (black solid line), and the purely Newtonian line, where $g_{\text{obs}} = g_{\text{bar}}$ (black dotted line). The colors of the models are shown in the legend. Panel (b) shows the MASR (BTFR) of the models (red diamonds), the best fit (black solid line) with uncertainties (gray area), the MASR (BTFR) in MOND assuming spherical symmetry (blue line), and the observational data from Lelli et al. (2016a; magenta area).</figcaption></figure>

$F/g_{\text{bar}}$ to the interpolation function $\nu$. The fit function could be precisely transformed into an interpolating function of the theory, which we could have implemented into PoR to ensure perfect agreement between the models and data. However, the simple interpolating function, which we chose to use in this work, is a good approximation to the Lelli et al. fit function (Famaey & McGaugh 2012); we demonstrate this here once again.

The observational fit is shown as the red area in Figure 19 with $g_{\ddagger} = 1.20 \pm 0.02 \times 10^{-10} \text{ m s}^{-2} = 3.87 \pm 0.06 \text{ pc Myr}^{-2}$ and a systematic error of $0.24 \times 10^{-10} \text{ m s}^{-2} = 0.77 \text{ pc Myr}^{-2}$. The best fit to the simulation data is the black line in panel (a) with $g_{\ddagger} = (1.217 \pm 0.006) \times 10^{-10} \text{ m s}^{-2} = 3.92 \pm 0.02 \text{ pc Myr}^{-2}$, which is well within the uncertainties of the observational data. Note, however, that we do not have galaxies that are deeply in the low-acceleration regime in their central parts, where the difference between the MDAR/RAR and QUMOND would be expected to be more important. The dashed black line shows $g_{\text{bar}} = g_{\text{obs}}$, which corresponds to a purely Newtonian gravitational acceleration even at the lowest accelerations. As can be seen, all Newtonian models lie on this line, which again tests the correct behavior of the code.

Panel (b) shows the MASR (BTFR), where $v_{\text{rot,flat}}$ is the flat rotation velocity of the models (red points); the blue line shows the theoretical MASR (BTFR) in MOND within the spherical approximation ($v_{\text{rot}} = (GM a_0)^{1/4}$) with $a_0 = 1.12 \times 10^{-10} \text{ m s}^{-2}$, which is well known to deviate somewhat from the true MOND predicted MASR (BTFR), as we confirm here. This is simply due to the fact that a flattened mass distribution spins slightly faster than the equivalent spherical mass distribution (McGaugh 2011) and that we use the flat rotation velocity and not the true asymptotic velocity. The black line

shows the best fit to the model data, with uncertainties shown as the gray area, and the magenta area shows the observational data from Lelli et al. (2016a). Again, the same fit function was used as in Lelli et al. (2016a; the baryonic mass, $M_{\text{bar}}$, is in units of $M_{\odot}$, $v_{\text{rot,flat}}$ is in units of km s$^{-1}$, and A is in units of $M_{\odot} (\text{sk m}^{-1})^{\gamma}$),

$$ \log_{10}(M_{\text{bar}}) = \gamma \times \log_{10}(v_{\text{rot,flat}}) + \log_{10}(A), \quad (27) $$

where $v_{\text{rot,flat}}$ is the rotation velocity of the flat part of the rotation curve, $M_{\text{bar}}$ the corresponding baryonic mass, $\log_{10}(A)$ the normalization, and $\gamma$ the slope of the relation. Lelli et al. find $\gamma = 3.71 \pm 0.08$ and $\log_{10}(A) = 2.27 \pm 0.18$, and the best fit to the simulation data is $\gamma = 3.97 \pm 0.04$ and $\log_{10}(A) = 1.78 \pm 0.08$. The best fit including uncertainties lies within the observational constraints, while the theoretical MASR (BTFR) in MOND is also covered by the fit.

After analyzing both relations separately, we will focus on the difference between the inputted value of $a_0$, $a_0 = 1.12 \times 10^{-10} \text{ m s}^{-2}$, and the one from the RAR, $a_0 = (1.217 \pm 0.006) \times 10^{-10} \text{ m s}^{-2}$. There are several reasons why this discrepancy may occur: first, the two interpolation functions used (the fit for the RAR and the one implemented in POR) can differ up to 5%, and we also do not expect an exact MDAR/RAR in QUMOND, meaning that we also do not expect the best-fit value for $a_0$ to be exactly the same as the inputted one. Future work performed at significantly higher numerical resolution will allow us to revisit this problem.

## 6. Discussion and Future Work

In this work we carried out the first full hydrodynamical simulations in a Milgromian framework to study the formation and evolution of single galaxies. Despite using very simple
initial conditions and a computationally constrained resolution,
our results are remarkably close to observations:

1. The work shows that late-type galaxies with exponen-
tially decreasing surface mass density profiles, thin and
thick disks, and bulge and Population II halos form
naturally in the MOND framework.

2. We compared simulations with simple and more complex baryonic physics to analyze whether the additions lead to major differences between the models. Surprisingly, the properties of simulations with and without more complex baryonic physics do not change significantly. This result underlines the hypothesis that the evolution of galaxies seems to follow a strict and simple law, where complex but realistic feedback processes play a very minor role (Kroupa 2015).

3. We set up simulations with gas spheres having different initial rotation laws, which resulted in comparable flat rotation velocities but different spatial extents and therefore different density distributions.

4. To check whether the occurrence of exponentially decreasing stellar surface mass density distributions is connected to our initial conditions, we computed several models again with the Newtonian Poisson solver. Indeed, only the Milgromian models show this feature, so the different gravitational theory and not the initial setup leads to an exponentially decreasing stellar surface mass density profile.

5. We also showed that the formation of a single galaxy in isolation has its limits with our method. Above a certain mass and rotation velocity scale, the initial collapse produces satellites in addition to the main galaxy, and, as mentioned in Section 4.4, this can also lead to the formation of a group of galaxies (N. Wittenburg et al., in preparation). Additionally, during the collapse of such models, a few protogalaxies form, which merge owing to gas dissipative processes. The resulting surface mass density profiles deviate somewhat from the single exponential form (Equation (12)).

6. The specific stellar angular momentum is slightly too small for the models as a whole, which is connected to insignificant later accretion and the simple initial conditions and may also be connected to the detailed prescriptions used to form stars as mentioned in Section 4.7. The models presented here assume an extremely small gas density beyond the primordial cloud such that accretion onto the collapsed and evolving galaxy models is negligible. There is substantial observational evidence that galactic disks can grow around bulges (i.e., early-formed elliptical galaxies or red nuggets) through accretion (de la Rosa et al. 2016; Mancini et al. 2019).

Overall our simulations suggest that major properties of galaxies, e.g., exponentially decreasing stellar surface mass density profiles, are a result of the dynamics produced by the gravitation potential, rather than the complexity of the baryonic processes in the Milgromian framework. It will be interesting to investigate whether the disk internal matter and angular momentum redistribution mechanism suggested by Herpich et al. (2017) is at work here or a completely different one. Important to note here is that these findings are somewhat surprising, as the effectively stronger gravity in MOND could have led to compact, pressure-supported systems that are well

within the Newtonian regime owing to a stronger collapse. We
also began running simulations with initial fluctuation of
the velocity and density distribution of up to 20% of the
unperturbed value. The first impression is that the properties of
the resulting galaxies are not changed significantly.

The recent observation of the gas kinematics of galaxies at a redshift of z = 6.8 revealed that these young objects are rotation dominated and not as turbulent as suggested before (Smit 2018). This supports our findings that baryonic processes like supernova explosions have a minor effect on the dynamics and morphology of a galaxy as a whole (even for the earliest population of galaxies).

We want to stress that all simulations done in the MOND framework show a tight MDAR/RAR, which compares with the observed MDAR/RAR and disks compatible with the size-mass relation. Hence, full hydrodynamical simulations in MOND produce objects that show dynamics comparable to real galaxies. The models, however, match neither the observed present-day SFHs nor the observed gas consumption timescales and have systematically slightly too small specific stellar angular momentum, but they do reproduce the slope of the Fall relation of observed galaxies and show that more massive galaxies have a smaller gas fraction, while also reproducing qualitatively the same picture in the KS diagram as real galaxies and initially matching the main sequence of star-forming galaxies.

Although the simulations done in this work are surprisingly
close to real galaxies in many aspects, we only made the first
steps in analyzing the formation and evolution of galaxies in
the MOND framework. Therefore, more simulations have to be
done with higher resolution, initially turbulent gas clouds, gas
accretion, an external field, and more baryonic processes (e.g.,
UV background radiation) in order to compare more observa-
tional results with our simulations. The calculations presented
here show that we can hope to constrain the initial conditions
needed in a MONDian cosmology, which in turn constrain
such a cosmology, and these conditions will get more precise
with more detailed simulations.

What is more, the formation of massive early-type galaxies is
also not understood today, as only 3%–4% of the observed
galaxies with a baryonic mass larger than about 10¹⁰ M☉ are
ellipticals (Delgado-Serrano et al. 2010). Among others, this
problem will be addressed in an upcoming publication (N.
Wittenburg et al., in preparation), where we will show how
early-type galaxies might form in the Milgromian framework.
Ellipticals could be results of rare mergers, or collisions that
strip the outer disks. But we can already conclude that elliptical
galaxies are an exception in MOND, since every collapsing
protogalactic gas cloud has angular momentum. The typical
outcome of galaxy formation in MOND is therefore a
rotationally supported disk galaxy, as is also observed in the
real universe (Delgado-Serrano et al. 2010).

B.F. acknowledges funding from the Agence Nationale de la Recherche (ANR project ANR-18-CE31-0006) and from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation program (grant agreement No. 834148).

**Appendix**

Tables 4–7 present the fit parameters for Equations, 12, 15, 17, 21, and 22 for every model.
<table><caption>**Table 4**
Fit Parameters for Equation (12) for All Models</caption><thead><tr><th>Model<br>Name/No.</th><th>n<sub>stars</sub> (M<sub>☉</sub>pc<sup>-2</sup>)</th><th>n<sub>gas</sub> (M<sub>☉</sub>pc<sup>-2</sup>)</th><th>r<sub>e,stars</sub> (kpc)</th><th>r<sub>e,gas</sub> (kpc)</th></tr></thead><tbody><tr><td>M1/1</td><td>168.61 ± 9.43</td><td>12.47 ± 0.93</td><td>1.61 ± 0.04</td><td>5.37 ± 0.10</td></tr><tr><td>M1sn/2</td><td>119.06 ± 7.89</td><td>8.32 ± 0.40</td><td>1.82 ± 0.06</td><td>6.28 ± 0.09</td></tr><tr><td>M1N/3</td><td>...</td><td>3.36 ± 0.12</td><td>...</td><td>9.62 ± 0.11</td></tr><tr><td>M1const/4</td><td>283.35 ± 7.53</td><td>2.05 ± 0.10</td><td>1.19 ± 0.01</td><td>6.44 ± 0.07</td></tr><tr><td>M2/5</td><td>845.68 ± 16.50</td><td>3.32 ± 0.31</td><td>1.35 ± 0.01</td><td>7.24 ± 0.16</td></tr><tr><td>M2sn/6</td><td>822.86 ± 17.59</td><td>10.55 ± 0.66</td><td>1.35 ± 0.01</td><td>5.80 ± 0.07</td></tr><tr><td>M2N/7</td><td>...</td><td>2.30 ± 0.12</td><td>...</td><td>8.55 ± 0.13</td></tr><tr><td>M3/8</td><td>443.41 ± 10.04</td><td>4.31 ± 0.23</td><td>1.98 ± 0.02</td><td>7.79 ± 0.10</td></tr><tr><td>M3sn/9</td><td>323.24 ± 16.04</td><td>24.02 ± 1.34</td><td>2.23 ± 0.06</td><td>5.72 ± 0.06</td></tr><tr><td>M4/10</td><td>...</td><td>20.39 ± 1.22</td><td>...</td><td>7.15 ± 0.09</td></tr><tr><td>M4sn/11</td><td>...</td><td>32.27 ± 1.54</td><td>...</td><td>6.40 ± 0.06</td></tr></tbody></table>

Note. Missing entries correspond to nonexponential behavior.

<table><caption>**Table 5**
Fit Parameters for Equation (15) for All Models</caption><thead><tr><th>Model<br>Name/No.</th><th>e<sub>1</sub> (M<sub>☉</sub>Myr<sup>-1</sup>)</th><th>f<sub>1</sub> (Gyr)</th><th>e<sub>2</sub> (M<sub>☉</sub>Myr<sup>-1</sup>)</th><th>f<sub>2</sub> (Gyr)</th></tr></thead><tbody><tr><td>M1/1</td><td>14.67 ± 0.56</td><td>0.554 ± 0.018</td><td>0.47 ± 0.02</td><td>2.53 ± 0.06</td></tr><tr><td>M1sn/2</td><td>8.73 ± 0.44</td><td>0.722 ± 0.035</td><td>0.46 ± 0.17</td><td>1.97 ± 0.18</td></tr><tr><td>M1N/3</td><td>33.46 ± 4.10</td><td>0.566 ± 0.023</td><td>0.0226 ± 0.0004</td><td>...</td></tr><tr><td>M1const/4</td><td>102.69 ± 9.12</td><td>0.207 ± 0.009</td><td>1.60 ± 0.09</td><td>1.03 ± 0.02</td></tr><tr><td>M2/5</td><td>2683.14 ± 273.4</td><td>0.108 ± 0.003</td><td>0.85 ± 0.07</td><td>1.56 ± 0.05</td></tr><tr><td>M2sn/6</td><td>1661.98 ± 160.30</td><td>0.122 ± 0.003</td><td>3.37 ± 0.10</td><td>0.97 ± 0.01</td></tr><tr><td>M2N/7</td><td>14522.20 ± 3331.00</td><td>0.128 ± 0.006</td><td>8.05 ± 0.60</td><td>0.85 ± 0.02</td></tr><tr><td>M3/8</td><td>426.77 ± 32.97</td><td>0.205 ± 0.006</td><td>2.14 ± 0.08</td><td>1.52 ± 0.03</td></tr><tr><td>M3sn/9</td><td>337.91 ± 29.49</td><td>0.222 ± 0.008</td><td>2.65 ± 0.05</td><td>1.55 ± 0.02</td></tr><tr><td>M4/10</td><td>8740.40 ± 1832.00</td><td>0.122 ± 0.008</td><td>7.13 ± 0.20</td><td>1.38 ± 0.02</td></tr><tr><td>M4sn/11</td><td>5019.00 ± 1828.00</td><td>0.146 ± 0.020</td><td>2.38 ± 0.06</td><td>2.67 ± 0.03</td></tr></tbody></table>

Note. The index 1 corresponds to the first exponential part (green line) and 2 to the second part (blue line). The second part for model M1N is constant and therefore fit not with Equation (15) but with a constant, c. e₂ represents c in this case.

<table><caption>**Table 6**
Fit Parameters for Equation (17) for All Models</caption><thead><tr><th>Model<br>Name/No.</th><th>i<sub>1</sub> (Gyr)</th><th>j<sub>1</sub> (Gyr)</th><th>i<sub>2</sub> (Gyr)</th><th>j<sub>2</sub> (Gyr)</th></tr></thead><tbody><tr><td>M1/1</td><td>0.29 ± 0.03</td><td>0.83 ± 0.04</td><td>4.11 ± 0.21</td><td>3.17 ± 0.07</td></tr><tr><td>M1sn/2</td><td>0.45 ± 0.02</td><td>0.89 ± 0.01</td><td>11.04 ± 4.03</td><td>2.63 ± 0.30</td></tr><tr><td>M1N/3</td><td>0.26 ± 0.03</td><td>0.83 ± 0.02</td><td>116.36 ± 2.21</td><td>...</td></tr><tr><td>M1const/4</td><td>0.07 ± 0.02</td><td>0.51 ± 0.04</td><td>0.76 ± 0.66</td><td>1.39 ± 0.03</td></tr><tr><td>M2/5</td><td>0.017 ± 0.007</td><td>0.31 ± 0.04</td><td>2.74 ± 0.26</td><td>2.10 ± 0.06</td></tr><tr><td>M2sn/6</td><td>0.031 ± 0.008</td><td>0.36 ± 0.03</td><td>0.74 ± 0.08</td><td>1.03 ± 0.03</td></tr><tr><td>M2N/7</td><td>0.013 ± 0.005</td><td>0.37 ± 0.04</td><td>0.44 ± 0.04</td><td>1.19 ± 0.03</td></tr><tr><td>M3/8</td><td>0.024 ± 0.006</td><td>0.30 ± 0.02</td><td>1.97 ± 0.95</td><td>1.94 ± 0.03</td></tr><tr><td>M3sn/9</td><td>0.025 ± 0.007</td><td>0.28 ± 0.02</td><td>1.52 ± 0.10</td><td>1.65 ± 0.03</td></tr><tr><td>M4/10</td><td>0.007 ± 0.003</td><td>0.24 ± 0.02</td><td>0.89 ± 0.03</td><td>1.67 ± 0.02</td></tr><tr><td>M4sn/11</td><td>0.009 ± 0.004</td><td>0.25 ± 0.03</td><td>2.35 ± 0.10</td><td>2.62 ± 0.04</td></tr></tbody></table>

Note. The index 1 corresponds to the first exponential part (green line) and 2 to the second part (blue line). The second part for model M1N is constant and therefore fit not with Equation (17) but with a constant, c. e₂ represents c in this case.
<table><caption>**Table 7**
Parameters for the Σ<sub>z</sub> Profile Fitted with Equation (21) for All Models</caption><thead><tr><th>Model<br>Name/No.</th><th>n<sub>z,thin</sub>(2 × 10<sup>7</sup> M<sub>☉</sub>pc<sup>-1</sup>)</th><th>n<sub>z,thick</sub>(2 × 10<sup>7</sup> M<sub>☉</sub>pc<sup>-1</sup>)<br>N: a(2 × 10<sup>7</sup> M<sub>☉</sub>pc<sup>-1</sup>)</th><th>z<sub>e,thin</sub>(kpc)</th><th>z<sub>e,thick</sub>(kpc)<br>N: b</th></tr></thead><tbody><tr><td>M1/1</td><td>1.610 ± 0.165</td><td>0.117 ± 0.005</td><td>0.144 ± 0.011</td><td>0.592 ± 0.009</td></tr><tr><td>M1sn/2</td><td>1.047 ± 0.067</td><td>0.231 ± 0.006</td><td>0.175 ± 0.011</td><td>0.495 ± 0.006</td></tr><tr><td>M1N/3</td><td>0.832 ± 0.034</td><td>0.028 ± 0.001</td><td>0.188 ± 0.007</td><td>2.304 ± 0.029</td></tr><tr><td>M1const/4</td><td>2.442 ± 0.190</td><td>0.248 ± 0.013</td><td>0.115 ± 0.004</td><td>0.315 ± 0.006</td></tr><tr><td>M2/5</td><td>8.488 ± 0.453</td><td>0.295 ± 0.026</td><td>0.117 ± 0.004</td><td>0.420 ± 0.014</td></tr><tr><td>M2sn/6</td><td>7.039 ± 0.407</td><td>0.237 ± 0.015</td><td>0.134 ± 0.005</td><td>0.479 ± 0.010</td></tr><tr><td>M2N/7</td><td>5.950 ± 0.241</td><td>0.075 ± 0.001</td><td>0.158 ± 0.005</td><td>2.929 ± 0.029</td></tr><tr><td>M3/8</td><td>5.874 ± 0.412</td><td>0.407 ± 0.022</td><td>0.149 ± 0.008</td><td>0.574 ± 0.015</td></tr><tr><td>M3sn/9</td><td>5.561 ± 0.346</td><td>0.344 ± 0.017</td><td>0.145 ± 0.008</td><td>0.639 ± 0.013</td></tr><tr><td>M4/10</td><td>27.290 ± 2.229</td><td>0.415 ± 0.013</td><td>0.139 ± 0.009</td><td>2.636 ± 0.037</td></tr><tr><td>M4sn/11</td><td>50.490 ± 4.984</td><td>0.389 ± 0.013</td><td>0.152 ± 0.012</td><td>2.664 ± 0.039</td></tr></tbody></table>

Note. The Newtonian models do not show an exponential profile after the thin disk; therefore, these regions are fitted with Equation (22).

## ORCID iDs

Nils Wittenburg https://orcid.org/0000-0001-9332-0000

Pavel Kroupa https://orcid.org/0000-0002-7301-3377

Benoit Famaey https://orcid.org/0000-0003-3180-9825

## References

Banik, I., Milgrom, M., & Zhao, H. 2018a, arXiv:1808.10545

Banik, I., O’Ryan, D., & Zhao, H. 2018b, MNRAS, 477, 4768

Banik, I., & Zhao, H. 2018, MNRAS, 473, 419

Barkana, R., & Loeb, A. 2001, PhR, 349, 125

Bekenstein, J., & Milgrom, M. 1984, ApJ, 286, 7

Bertone, G., & Tait, T. M. P. 2018, Natur, 562, 51

Bigiel, F., Leroy, A., Walter, F., et al. 2008, AJ, 136, 2846

Bílek, M., Thies, I., Kroupa, P., & Famaey, B. 2018, A&A, 614, A59

Bleuler, A., & Teyssier, R. 2014, MNRAS, 445, 4015

Brada, R., & Milgrom, M. 1995, MNRAS, 276, 453

Brada, R., & Milgrom, M. 1999, ApJ, 519, 590

Candlish, G. N., Smith, R., & Fellhauer, M. 2015, MNRAS, 446, 1060

Catelan, P., & Theuns, T. 1996, MNRAS, 282, 436

Combes, F. 2014, A&A, 571, A82

Courty, S., & Alimi, J. M. 2004, A&A, 416, 875

Dalcanton, J. J., Spergel, D. N., & Summers, F. J. 1997, ApJ, 482, 659

de Blok, W. J. G., & McGaugh, S. S. 1997, MNRAS, 290, 533

de la Rosa, I. G., La Barbera, F., Ferreras, I., et al. 2016, MNRAS, 457, 1916

Delgado-Serrano, R., Hammer, F., Yang, Y. B., et al. 2010, A&A, 509, A78

Desmond, H. 2017a, MNRAS, 464, 4160

Desmond, H. 2017b, MNRAS, 472, L35

Di Cintio, A., & Lelli, F. 2016, MNRAS, 456, L127

Disney, M. J., Romano, J. D., Garcia-Appadoo, D. A., et al. 2008, Natur, 455, 1082

Dubois, Y., & Teyssier, R. 2008, A&A, 477, 79

Dutton, A. A. 2009, MNRAS, 396, 121

Efstathiou, G., & Jones, B. J. T. 1979, MNRAS, 186, 133

Einstein, A. 1916, AnP, 354, 769

Faber, S. M., & Jackson, R. E. 1976, ApJ, 204, 668

Fall, S. M., & Efstathiou, G. 1980, MNRAS, 193, 189

Famaey, B., & McGaugh, S. S. 2012, LRR, 15, 10

Freeman, K. C. 1970, ApJ, 160, 811

Herpich, J., Tremaine, S., & Rix, H.-W. 2017, MNRAS, 467, 5022

Hubble, E. P. 1926, ApJ, 64, 321

Ibata, R. A., Lewis, G. F., Conn, A. R., et al. 2013, Natur, 493, 62

Javanmardi, B., & Kroupa, P. 2020, MNRAS, 493, 44

Jeřábková, T., Hasani Zonoozi, A., Kroupa, P., et al. 2018, A&A, 620, A39

Keller, B. W., & Wadsley, J. W. 2017, ApJL, 835, L17

Kennicutt, R. C. J. 1998, ApJ, 498, 541

Kroupa, P. 2012, PASA, 29, 395

Kroupa, P. 2015, CaJPh, 93, 169

Kroupa, P., Famaey, B., de Boer, K. S., et al. 2010, A&A, 523, A32

Kroupa, P., Theis, C., & Boily, C. M. 2005, A&A, 431, 517

Kroupa, P., Weidner, C., Pflamm-Altenburg, J., et al. 2013, in The Stellar and Sub-stellar Initial Mass Function of Simple and Composite Populations, ed. T. D. Oswalt & G. Gilmore (Dordrecht: Springer), 115

Lange, R., Driver, S. P., Robotham, A. S. G., et al. 2015, MNRAS, 447, 2603

Larson, R. B. 1981, MNRAS, 194, 809

Lelli, F., Fraternali, F., & Verheijen, M. 2013, MNRAS, 433, L30

Lelli, F., McGaugh, S. S., & Schombert, J. M. 2016a, ApJL, 816, L14

Lelli, F., McGaugh, S. S., Schombert, J. M., & Pawlowski, M. S. 2016b, ApJL, 827, L19

Lelli, F., McGaugh, S. S., Schombert, J. M., & Pawlowski, M. S. 2017, ApJ, 836, 152

Llinares, C., Knebe, A., & Zhao, H. 2008, MNRAS, 391, 1778

Ludlow, A. D., Benítez-Llambay, A., Schaller, M., et al. 2017, PhRvL, 118, 161103

Lüghausen, F., Famaey, B., & Kroupa, P. 2015, CaJPh, 93, 232

Mancini, C., Daddi, E., Juneau, S., et al. 2019, MNRAS, 489, 1265

McGaugh, S. 2011, arXiv:1109.1599

McGaugh, S. 2014, Galax, 2, 601

McGaugh, S. S. 2004, ApJ, 609, 652

McGaugh, S. S. 2005, ApJ, 632, 859

McGaugh, S. S. 2012, AJ, 143, 40

McGaugh, S. S., Schombert, J. M., Bothun, G. D., & de Blok, W. J. G. 2000, ApJL, 533, L99

Metz, M., Kroupa, P., & Jerjen, H. 2007, MNRAS, 374, 1125

Milgrom, M. 1983, ApJ, 270, 365

Milgrom, M. 1999, PhLA, 253, 273

Milgrom, M. 2009, ApJ, 698, 1630

Milgrom, M. 2010, MNRAS, 403, 886

Milgrom, M. 2014, SchpJ, 9, 31410

Mor, R., Robin, A. C., Figueras, F., & Antoja, T. 2018, A&A, 620, A79

Müller, O., Pawlowski, M. S., Jerjen, H., & Lelli, F. 2018, Sci, 359, 534

Navarro, J. F., Benítez-Llambay, A., Fattahi, A., et al. 2017, MNRAS, 471, 1841

Newton, I. 1687, Philosophiae Naturalis Principia Mathematica. Auctore Js. Newton

Oh, S., & Kroupa, P. 2016, A&A, 590, A107

Okazaki, T., & Taniguchi, Y. 2000, ApJ, 543, 149

Pawlowski, M. S. 2018, MPLA, 33, 1830004

Pawlowski, M. S., & Kroupa, P. 2020, MNRAS, 491, 3042

Pawlowski, M. S., Pflamm-Altenburg, J., & Kroupa, P. 2012, MNRAS, 423, 1109

Pflamm-Altenburg, J., & Kroupa, P. 2009, ApJ, 706, 516

Posti, L., Fraternali, F., Di Teodoro, E. M., & Pezzulli, G. 2018, A&A, 612, L6

Renaud, F., Famaey, B., & Kroupa, P. 2016, MNRAS, 463, 3637

Rosdahl, J., Blaizot, J., Aubert, D., Stranex, T., & Teyssier, R. 2013, MNRAS, 436, 2188

Sancisi, R. 2004, in IAU Symp. 220, Dark Matter in Galaxies, ed. S. Ryder (San Francisco, CA: ASP), 233

Sanders, R. H. 1990, A&ARv, 2, 1

Sanders, R. H. 1998, MNRAS, 296, 1009

Serra, P., Oosterloo, T., Cappellari, M., den Heijer, M., & Józsa, G. I. G. 2016, MNRAS, 460, 1382

Smit, R. 2018, AAS Meeting Abstracts, 231, 454.03

Smolin, L. 2017, PhRvD, 96, 083523
Speagle, J. S., Steinhardt, C. L., Capak, P. L., & Silverman, J. D. 2014, ApJS, 214, 15

Struck, C., & Elmegreen, B. G. 2017, MNRAS, 469, 1157

Struck, C., & Elmegreen, B. G. 2018, ApJL, 868, L15

Sutherland, R. S., & Dopita, M. A. 1993, ApJS, 88, 253

Tamburri, S., Saracco, P., Longhetti, M., et al. 2014, A&A, 570, A102

Teyssier, R. 2002, A&A, 385, 337

Teyssier, R., Pontzen, A., Dubois, Y., & Read, J. I. 2013, MNRAS, 429, 3068

Thomas, G. F., Famaey, B., Ibata, R., et al. 2018, A&A, 609, A44

Tiret, O., & Combes, F. 2007, A&A, 464, 517

Tiret, O., & Combes, F. 2008, A&A, 483, 719

Truelove, J. K., Klein, R. I., McKee, C. F., et al. 1998, ApJ, 495, 821

Tully, R. B., & Verheijen, M. A. W. 1997, ApJ, 484, 145

van Albada, T. S., & Sancisi, R. 1986, RSPTA, 320, 447

Voglis, N. 1994, in Lecture Notes in Physics, Galactic Dynamics and N-Body Simulations, ed. G. Contopoulos, N. K. Spyrou, & L. Vlahos (Berlin: Springer), 365

Wesson, P. S. 1985, A&A, 151, 105

Wu, X., & Kroupa, P. 2015, MNRAS, 446, 330