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

HAL
open science

# Globular cluster orbital decay in dwarf galaxies with MOND and CDM: Impact of supernova feedback

Michal Bílek, Françoise Combes, Srikanth T. Nagesh, Michael Hilker

► To cite this version:

Michal Bílek, Françoise Combes, Srikanth T. Nagesh, Michael Hilker. Globular cluster orbital decay in dwarf galaxies with MOND and CDM: Impact of supernova feedback. Astronomy & Astrophysics - A&A, 2024, 690, pp.A119. 10.1051/0004-6361/202449666. insu-04726360

HAL Id: insu-04726360

https://insu.hal.science/insu-04726360v1

Submitted on 8 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Globular cluster orbital decay in dwarf galaxies with MOND and CDM: Impact of supernova feedback

Michal Bílek<fnref n="1" />,<fnref n="2" />,*, Françoise Combes<fnref n="1" />,*, Srikanth T. Nagesh<fnref n="3" />, and Michael Hilker<fnref n="4" />

Received 19 February 2024 / Accepted 11 June 2024

## ABSTRACT

Dynamical friction works very differently for Newtonian gravity with dark matter and in modified Newtonian dynamics (MOND). While the absence of dark matter considerably reduces the friction in major galaxy mergers, analytic calculations indicate the opposite for very small perturbations, such as globular clusters (GCs) sinking in dwarf galaxies. Here, we study the decay of GCs in isolated gas-rich dwarf galaxies using simulations with the Phantom of Ramses code, which enables both the Newtonian and the QUMOND MOND gravity. We modeled the GCs as point masses, and we simulated the full hydrodynamics, with star formation and supernovae feedback. We explored whether the fluctuations in gravitational potential caused by the supernovae can prevent GCs from sinking toward the nucleus. For GCs of typical mass or lighter, we find that this indeed works in both Newtonian and MOND simulations. The GC can even make a random walk. However, we find that supernovae cannot prevent massive GCs ($M \ge 4 \times 10^5 M_\odot$) from sinking in MOND. The resulting object looks similar to a galaxy with an offset core, which embeds the sunk GC. The problem is much milder in the Newtonian simulations. This result thus favors Newtonian over QUMOND gravity, but we note that it relies on the correctness of the difficult modeling of baryonic feedback. We propose that the fluctuations in the gravitational potential could be responsible for the thickness of the stellar disks of dwarf galaxies and that strong supernova winds in modified gravity can transform dwarf galaxies into ultra-diffuse galaxies.

**Key words.** galaxies: dwarf – galaxies: evolution – galaxies: kinematics and dynamics – galaxies: star clusters: general – dark matter

## 1. Introduction

The missing mass problem at multi-scale from galaxies to large-scale structures in the Universe is one of the most stubborn problems in astrophysics. Since its discovery in galaxy clusters (Zwicky 1937), then in galaxies through rotation curves (Rubin et al. 1980), and its confirmation through gravitational lenses and cosmic background radiation (e.g. Planck Collaboration XIII 2016), many candidates for a dark matter particle have been proposed. Perhaps the most favored candidate for cold dark matter (CDM) has been neutralino, the most stable particle from supersymmetry (e.g. Roszkowski et al. 2018). However, no supersymmetric particles have been discovered yet at the Large Hadron Collider of CERN; therefore, other candidates have been seriously considered, such as the axion or axion-like particles (ALPs; Hui et al. 2017).

Alternatives have been proposed in terms of modified gravity, and one of the most successful is modified Newtonian dynamics (MOND), which was designed 40 years ago by Milgrom (1983a,b,c). MOND has several possible forms and can also be modified inertia. In our subsequent analysis, we chose modified gravity in its special flavor: quasi-linear modified Newtonian dynamics (QUMOND). There exist well-known problems of MOND in galaxy clusters (e.g., Famaey & McGaugh 2012) and in the early universe, where the acoustic peaks of the cosmic microwave background (CMB) cannot be repro-

duced, as well as in the power spectrum of large-scale structures. However, these issues might be in the process of being solved (Skordis & Złośnik 2021). The predictions from MOND have been tested many times and compared to those from the standard CDM model for galaxies in the hope of choosing the model that better fits the observations. For example, the existence and frequency of bars (Tiret & Combes 2007), galaxy interactions (Tiret & Combes 2008), and dynamical friction (Kroupa 2015; Bílek et al. 2019, 2021) have been compared. Recently, Freundlich et al. (2022) reported a failure of MOND in the Coma cluster, where ultra-diffuse galaxies (UDGs) reveal kinematics in accord with MOND as isolated galaxies, while the external field effect (EFE) should exist for galaxies in a cluster and reduce their expected internal velocity, with respect to the isolated case. But the problem of MOND in clusters has been known for a long time, and it awaits a new physical phenomena, such as dark baryons and/or screening effects, due to the introduction of another characteristic scale (Milgrom 2023a) expected to solve both problems.

Another test for MOND is based on the observation of globular clusters (GCs) in dwarf galaxies, where the dynamical friction should be high and drive the GCs quickly into the nucleus (Nipoti et al. 2008). Dynamical friction is a phenomenon very different in MOND and CDM. In the deep MOND regime, for very light galaxies and very small perturbations, such as a GC in a dwarf, analytical formulations have shown that the friction could be higher in MOND than in CDM (Ciotti & Binney 2004). However, for large perturbations, such as those in mergers of two

<footnote n="1">¹ Observatoire de Paris, LERMA, Collège de France, CNRS, PSL University, Sorbonne University, F-75014 Paris</footnote>

<footnote n="2">² FZU – Institute of Physics of the Czech Academy of Sciences, Na Slovance 1999/2, Prague 182 21 Czech Republic</footnote>

<footnote n="3">³ Strasbourg University, CNRS, Observatoire astronomique de Strasbourg, F-67000 Strasbourg, France</footnote>

<footnote n="4">⁴ European Southern Observatory, Karl-Schwarzschild-Strasse 2, 85748 Garching bei München, Germany</footnote>

<footnote>* Corresponding author; michal.bilek@obspm.fr</footnote>
Milky Way-type galaxies, dynamical friction is much reduced compared to Newtonian gravity with dark matter, as shown by Tiret & Combes (2008) and Combes (2014). A simple analytical formula has been proposed for the deep MOND dynamical friction on small perturbers in Sánchez-Salcedo et al. (2006) and tested through simulations by Bílek et al. (2021). Such analytic calculations indicate that GCs should sink from their currently observed positions into the centers of low-mass, low-surface-brightness galaxies on the timescale of one gigayear. This contradicts the observations of 10-gigayear-old GCs in these galaxies, unless fine-tuning is invoked.

Bílek et al. (2021) have considered the decay of single or multiple GCs in an isolated UDG in the deep MOND regime. Using the Phantom of Ramses code (POR; Lüghausen et al. 2015), they simulated the sinking GCs modeled as point masses to check whether all GCs could sink to the center of the system. What happened is that the GC sinking stopped at a fraction of a kiloparsec due to core stalling (Hernandez & Gilmore 1998). The apocenter of the GC first reduced rather quickly, and the decay slowed down when the apocenter reachd 0.5 kpc. This phenomenon is due to the existence of a core in the UDG density profile, implying a harmonic potential, with particles not being able to absorb the energy of the GC.

In these simulations, the interaction between GCs was not taken into account, neither the possible stripping nor the destruction of the sinking GC through the UDG tidal forces. Also, UDGs were simulated with collisionless particles without gas. In the present paper, we consider sinking GCs in isolated dwarfs, which are smaller and lighter gas-rich galaxies, with effective radii as low as 0.5 kpc, and the galaxies rotate. The results of Bílek et al. (2021) indicate that the dynamical friction timescale decreases with the mass of the galaxy, and thus we tested the survivability of GCs for the least massive galaxies that are known to have GCs. We show that the introduction of gas dynamics, star formation, and feedback can have a substantial effect on the dynamics of GCs. Medium- and low-mass GCs are totally prevented from sinking. The feedback, however, is not strong enough, at least in our simulations, to solve the problem of fast sinking for the massive GCs.

In Sect. 2, we compile observations of dwarf galaxies and their GC populations in order to define the initial conditions of the simulations. The technical details of the simulations are described in Sect. 3. The results regarding the survivability of the GCs are displayed in Sect. 4, and in Sect. 5, we discuss other interesting phenomena seen in the simulations. Sect. 6 presents a summary of our work.

## 2. Properties of observed isolated dwarf galaxies and of their GCs

In order to make simulations of relevant objects, we first compiled information about the characteristics of the observed isolated dwarfs that are known to have GCs. It turns out that there are not many studies on GCs of isolated dwarfs. We took our data from the works of Sharina et al. (2005) and Georgiev et al. (2009). The authors of each work were looking for GC candidates in Hubble Space Telescope (HST) images of nearby dwarf galaxies ($2 < D < 11$ Mpc with a few exceptions).

Regarding the selection criteria, Sharina et al. (2005) and Georgiev et al. (2009) chose their GC candidates on the basis of their absolute magnitudes, colors, and structural photometric parameters since they appear spatially resolved in the HST images. We included in our galaxy sample only the galaxies marked by the authors to be located in the field. Sharina et al.

(2005) estimated the contaminant fraction to be at most 10% among their GC candidates. Georgiev et al. (2009) stated that their sample contains at most two contaminants per field. Nevertheless, the GCs of these nearby galaxies appear as resolved objects, which is one of the important selection criterion. Therefore, the contaminants are rather expected among the low-luminosity GC candidates.

It was then necessary to obtain the characteristics of the galaxies. Their HI masses and distances were taken from the online database of the Catalog and Atlas of the Local Volume galaxies<fnref n="1" /> (the LV database hereafter, Karachentsev & Kaisina 2019) in its form from April 1, 2022. Both Sharina et al. (2005) and Georgiev et al. (2009) provided absolute magnitudes of the galaxies in the V band. We converted the absolute magnitudes to stellar masses, assuming the absolute magnitude of the Sun in V of 4.83 (Binney & Merrifield 1998) and the stellar mass-to-light ratio of one. We neglected the mass of the possible molecular gas, helium, and metals. Because the distances of the galaxies are generally different in Sharina et al. (2005), Georgiev et al. (2009), and the LV database, we adopted for all galaxies the distances from the most recent source, namely, from the LV database. Accordingly, we corrected the stellar masses of the galaxies and the projected galactocentric distances of the GCs. The distance in the respective sources and in the LV database are given in Table B.1. The median deviation of the Sharina et al. (2005) distances from the LV distances is 0.03 dex (ca. 10%), and for the Georgiev et al. (2009) data, this deviation is 0.006 dex (ca. 1%).

The properties of the isolated dwarf galaxies are shown in Fig. 1 as a function of the stellar plus HI mass of the galaxy. The figure shows that GCs are common in isolated dwarf galaxies more massive than ~$10^8 M_\odot$. The least massive galaxy that has a GC candidate has a mass of $1 \times 10^7 M_\odot$. The GCs have projected galactocentric distances of a few kiloparsecs. Sharina et al. (2005) and Georgiev et al. (2009) do not list the stellar disk scale lengths. However, after using Eq. 1 derived below to estimate the stellar disk scale lengths, we found that the projected galactocentric radial distances of the GCs are comparable to stellar disk scale lengths. Sharina et al. (2005) nevertheless warns that their list of GC candidates can omit GCs that are far from their hosts because of the limited field of view of the HST camera. In the fourth panel of Fig. 1, we show the estimated mass of the GC candidates, assuming a V-band mass-to-light ratio of 2.2 (Spitler & Forbes 2009). We note that this might be an overestimate because dwarf galaxies can also contain young GCs with lower mass-to-light ratios (Pace et al. 2021). The last panel tells us that isolated dwarfs are usually gas-rich objects. We found that the median HI mass fraction of the galaxies is 0.49. The median HI mass fraction of the galaxies with GCs is the same.

It was also necessary to estimate the sizes of the galaxies because Sharina et al. (2005) and Georgiev et al. (2009) do not list them. We estimated them from the mass-size relation. The relation for galaxies of this type was extracted from the database of the LITTLE THINGS survey of gas-bearing dwarf irregular galaxies (Hunter et al. 2012). Hunter et al. (2021) listed for these galaxies the half-mass radii of the HI distribution from Sérsic fits. We converted them to the estimates of scale-lengths by dividing them by 1.68, as if the HI gas were organized into exponential disks. Fitted exponential stellar disk scale-lengths were listed by Hunter et al. (2021) directly. The paper, however, does not list the stellar and HI masses of the galaxies. To estimate them, we converted the V band magnitudes of the galaxies that they give to stellar masses, assuming the stellar mass-to-light

<footnote n="1">¹ https://relay.sao.ru/lv/lvgdb/</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 1. Properties of observed isolated dwarf galaxies. From top to bottom: (1) Number of GCs in each galaxy, (2) projected galactocentric radius of each GC in each galaxy ($R_{GC}$) in the units of kiloparsecs, (3) same but in units of the scale lengths of the stellar disks of the galaxies ($R_{D,est}$) estimated from Eq. 1, (4) stellar mass of each GC ($M_{GC}$) in each galaxy, and (5) HI mass fraction in the total baryonic mass of each galaxy. The quantity on the horizontal axis is the total mass in stars and gas of each galaxy. The colors in the middle three tiles help distinguish between the GCs of different galaxies.</figcaption></figure>

ratio of one. We assumed, inspired by our main galaxy sample, an HI mass fraction of 0.5.

The results are shown in Fig. 2. The scale lengths of the HI and stellar disks differ, in median, by 0.2 dex; that is, the HI disk

<figure><img src="image_4.png" /><figcaption>Fig. 2. Properties of observed isolated dwarf galaxies. From top to bottom: (1) Ratio of the scale lengths of the HI and stellar disks and (2) scale length of the stellar disk. These quantities are plotted as functions of the total mass in stars and gas of the galaxy. The dashed line in the bottom panel shows the best linear fit given by Eq. 1.</figcaption></figure>

size is typically 1.6 times the stellar scale length. The galaxies follow a well-defined mass-size relation. The best linear fit of the logarithms of the baryonic mass and radius reads

$$ \log_{10} R_{D,*} = -2.74 + 0.314 \log_{10}(M_{*} + M_{HI}). \quad (1) $$

The stellar size-mass relation is rather independent of redshift for low-mass star-forming galaxies (Nedkova et al. 2021). The axial ratio of HI disks of dwarf galaxies is about 0.6 (Roychowdhury et al. 2010).

To put the GCs of the considered dwarfs in a context of galaxies with known, well-sampled GC systems, we note that the absolute magnitudes of the GCs of a given galaxy generally follow a Gaussian distribution. While the width of the distribution can vary with the mass of the galaxy, the distribution usually peaks near the absolute V-magnitude of -7.5 (Rejkuba 2012). Using the same assumptions as before, this corresponds to the mass of the GC of $2 \times 10^5 M_{\odot}$. This is close to the typical mass of the GCs in our sample.

Figure 3 shows the relation between the mass of the GC and its projected distance to the galaxy center. This distance is expressed in the units of the estimated stellar effective radius of each host galaxy. The line with error bars in the figure indicates averages in bins and the uncertainties. If the dynamical friction strongly affects the GCs, then we expect the more massive GCs to be on average at lower galactocentric distances than the less massive GCs. Indeed, there is a tendency for such a correlation. The Spearman test gives that the probability that there is no correlation is just 0.2%.
<figure><img src="image_5.png" /><figcaption>Fig. 3. Mass versus projected galactocentric radius for the observed GCs. The GCs have the same colors as in Fig. 1.</figcaption></figure>

Our host galaxies hold up to 3% (10%) of their baryonic (stellar) mass in GCs.

### 3. Numerical methods

We used the POR code to run all the simulations. The code is a patch to the publicly available adaptive mesh refinement (AMR) grid-based code RAMSES (Teyssier 2002). The POR patch was added to the 2015 version of RAMSES by Lüghausen et al. (2015), and it implements the quasi linear formalism of MOND (QUMOND; Milgrom 2010). The governing equation of QUMOND is

$$ \nabla \cdot \mathbf{g} = \nabla \cdot (\mathbf{v}g_N), \qquad (2) $$

where $g$ and $g_N$ are the true and Newtonian acceleration vectors. At each step, the Newtonian acceleration $g_N$ is obtained using only the baryonic distribution, and subsequently, the interpolating function $v$ is determined at every step, which enables the code to compute the source term $g$. The “simple” form of the interpolating function used in POR is

$$ v = \frac{1}{2} + \sqrt{\frac{1}{4} + \frac{a_0}{g_N}}. \qquad (3) $$

It relates the MOND and Newtonian regimes (Famaey & Binney 2005), and it has also been shown to work observationally (Gentile et al. 2011; Iocco et al. 2015; Banik et al. 2018; Chae et al. 2018).

#### 3.1. Hydrodynamics, star formation, and feedback

As mentioned earlier, POR only modifies the Poisson gravity solver from a Newtonian one to a MOND one, but it inherits the hydrodynamical solver and star formation and feedback recipes from the 2015 version of RAMSES. A second-order Godunov scheme with a Riemann solver for the conservative Euler equations are used (Teyssier 2002).

Star formation was modeled according to the Schmidt law (see, e.g., Shi et al. 2011 for a discussion). Schmidt (1959) proposed that the stellar mass formed per unit time in unit volume, $\dot{\rho}_*$, is related to the local density of interstellar material, $\rho$, as

a power law $\dot{\rho}_* \propto \rho^\nu$. As it was found observationally, the projected versions of these quantities are related as $\dot{\Sigma}_* \propto \Sigma^\nu$, where $\nu \approx 1.4$ (Kennicutt 1989). This makes sense because, from the dynamical point of view in Newtonian gravity, we expect a relation of the type

$$ \dot{\rho}_* = \epsilon_* \frac{\rho}{t_{\text{ff}}}, \qquad (4) $$

where the free fall time is $t_{\text{ff}} = \sqrt{3\pi/32G\rho}$ (e.g., Katz 1992). In POR, it is possible to prescribe star formation either by using Eq. 4 (Teyssier et al. 2010, 2013) or the equivalent form (Rasera & Teyssier 2006; Dubois & Teyssier 2008)

$$ \dot{\rho}_* = \frac{\rho}{t_{\text{dep}}}, \qquad (5) $$

which we used in our study. Introducing $n$ as the number density of hydrogen atoms corresponding to $\rho$, the gas depletion timescale in Eq. 5, $t_{\text{dep}}$, is calculated from the user-supplied parameters $t_*$ and $n_*$ as $t_{\text{dep}} = t_* (\frac{n}{n_*})^{-1/2}$. At the same time, the parameter $n_*$, with the dimension of the number density of hydrogen atoms, also plays the role of the threshold above which star formation is enabled. This is motivated by the observed fact that at low gas densities, the star formation rate is much lower than expected from extrapolating the Schmidt relation derived for denser gas. Both $t_*$ and $n_*$ are theoretically poorly constrained, depend on the resolution of the simulation, and have to be tuned for the specific purpose. To derive the expression that POR uses to model the star formation, both sides of Eq. 5 have to be multiplied by the volume of the computational cell $V_{\text{cell}}$ and by the timestep $dt$, and divided by the user-supplied mass of a star particle, $m_*$. In turn, the ideal number of stellar particles formed in a given time step would be

$$ dN_* \equiv \frac{\dot{\rho}_* V_{\text{cell}} dt}{m_*} = \frac{m_{\text{cell}}}{m_*} \left(\frac{n_{\text{cell}}}{n_*}\right)^{1/2} \frac{dt}{t_*}, \qquad (6) $$

where $m_{\text{cell}} = \rho_{\text{cell}} V_{\text{cell}}$ is the gas mass in the cell and $n_{\text{cell}}$ is the number density of hydrogen atoms in the cell. The quantity $dN_*$ is generally not and integer. Therefore, the actual number of stars formed in the given cell is drawn from a Poisson distribution that has the mean of $dN_*$. In POR, Eq. 5 is activated by setting $t_* > 0$; otherwise, Eq. 4 is used.

Supernova (SN) feedback plays a vital role in regulating star formation in a galaxy. Different types of feedback prescriptions available in POR are described in detail in Sect. 2.2 of Nagesh et al. (2023), of which the intermediate feedback (Sect. 2.2.1 of Nagesh et al. 2023, originally implemented by Dubois & Teyssier 2008) is used in the present work. The intermediate feedback prescription allows the user to specify the fraction of energy from the SNe explosion to be injected back into the interstellar medium as kinetic energy, which is carried by a SN blast wave with a user-specified blast radius, $r_{\text{bubble}}$, where velocity of the wave is computed using a local Sedov blast wave solution. Every time a stellar particle forms, metallicity, energy, and momentum are released into the respective cells. As far as the mass is concerned, whenever a stellar particle with mass $m_*$ is formed, $m_*(1 + \eta_{sn})$ is removed from the gas cell, with $\eta_{sn}$ being a user specified fraction of mass of the new stellar particle that goes into the SNe. Subsequently, $m_* \eta_{sn}$ is injected back into the interstellar medium after a delay of $t_{sn}$ (Dubois & Teyssier 2008). The feedback parameters used are listed in Appendix D.

Since its development, POR has been applied on diverse scenarios using N-body-only polar ring galaxies (Lüghausen et al.
<table><caption>**Table 1.** Parameters of the basic MOND models used to generate the initial conditions in MOND DICE.</caption><thead><tr><th>Name</th><th><i>M</i><sub>bar</sub><br>[10<sup>8</sup> <i>M</i><sub>☉</sub>]</th><th><i>r</i><sub>D</sub><br>[kpc]</th><th><i>q</i></th><th><i>f</i><sub>gas</sub></th><th><i>T</i><br>[K]</th><th><i>Q</i></th><th><i>N</i></th><th><i>r</i><sub>GC</sub><br>[kpc]</th></tr></thead><tbody><tr><td>GDw1e8</td><td>1</td><td>0.8</td><td>0.4</td><td>0.9</td><td>25000</td><td>0.9<sup>(f)</sup></td><td>10<sup>5</sup></td><td>1</td></tr><tr><td>GDw2e7</td><td>0.2</td><td>0.5</td><td>0.4</td><td>0.9</td><td>25000</td><td>0.9<sup>(f)</sup></td><td>10<sup>5</sup></td><td>0.5</td></tr><tr><td>GUDG1e8</td><td>1</td><td>2</td><td>0.5</td><td>0.9</td><td>25000</td><td>1.25<sup>(m)</sup></td><td>10<sup>5</sup></td><td>2</td></tr><tr><td>GUDG2e7</td><td>0.2</td><td>2</td><td>0.5</td><td>0.9</td><td>21000</td><td>1.25<sup>(m)</sup></td><td>10<sup>5</sup></td><td>2</td></tr></tbody></table>

**Notes.** *M*<sub>bar</sub>: Total baryonic mass. *r*<sub>D</sub>: Disk scale length (the same for the stellar and gas disk). *q*: Axial ratio of the disk. *f*<sub>gas</sub>: Gas fraction. *T*: Gas temperature. *Q*: Toomre parameter (type: (f)fixed, (m)minimal). *N*: Number of particles in the stellar disk. *r*<sub>GC</sub>: Initial distance of the GC from the galaxy center.

2013), stellar streams (Thomas et al. 2017), the Local Group (Bílek et al. 2018), and galaxies with shells (Bílek et al. 2022), as well as hydrodynamical simulations, the antennae pair of interacting galaxies (Renaud et al. 2016), disk stability of M33 in MOND (Banik et al. 2020), the UDG AGC 114905 with star formation (Banik et al. 2022a), and the satellite plane (Banik et al. 2022b). The code POR has also been used to run a cosmological hydrodynamical simulation in the MONDian framework (Wittenburg et al. 2023). Nagesh et al. (2023) compared the effect of two types of feedback on galaxy models with stellar mass ranging between $10^7 - 10^{11} M_☉$ and have found that the star formation rates and gas consumption timescale are in agreement with observations. Similarly, the effect of these feedback prescriptions on a monolithic cloud collapse of a rotating (non-rotating) gas cloud into disk (elliptical) galaxies in MOND have been investigated, too (Wittenburg et al. 2020; Eappen et al. 2022). Simulations of the formation of disk galaxies in MOND automatically lead to the correct scale lengths and exponential disks (Wittenburg et al. 2020). The monolithic collapse of non-rotating gas clouds can naturally explain the emergence of compact massive relic galaxies within a MOND-based universe (Eappen & Kroupa 2024). A user manual to set up N-body and hydrodynamical isolated disk galaxies in MOND using POR (Nagesh et al. 2021), as well as with POR and other relevant packages, is available here<fnref n="2" />.

## 3.2. Generating initial conditions

The Disk Initial Condition Environment (DICE; Perret et al. (2014)), adapted to MOND gravity by Banik et al. (2020), was used to generate the initial conditions for the MOND simulations<fnref n="2" />. The details of the implementation of the MOND gravity in DICE are discussed in Banik et al. (2020). For the prescribed properties of the galaxy to be prepared, the code returns a list of positions, velocities, and masses of the initial stellar particles and a table with the rotation curve. Based on the rotation curve provided by the MOND DICE code, the POR_HYDRO patch in the BonnPoR package<fnref n="2" /> calculates and generates the necessary distribution of the gas in the galaxy (Banik et al. 2022b). We used the MAGI code (Miki & Umemura 2018) to generate the initial conditions for the Newtonian models. We used the code to generate three components in the galaxy: the stellar disk, the gas disk, and a dark NFW halo (Navarro et al. 1996). The characteristic density and scale radius of the halo were chosen such that the galaxies follow the mean stellar-to-halo mass relation (Behroozi et al. 2013) and the halo mass-concentration relation (Diemer & Kravtsov 2015). The stellar mass to generate these values was assumed to be half of the total initial baryonic mass

of the galaxies because this is the typical value observed for our galaxies (Sect. 2). The galaxies created by MAGI consist only of particles. In order to include gas, we calculated the rotation curve of the galaxy, removed the gas-disk particles, and let POR deduce and generate the distribution of gas, as in the MOND case. The stability of a halo created by MAGI is demonstrated in Appendix C. We modeled the GCs as point masses. In most of this work, the GCs had a mass of $10^5 M_☉$, which is about the typical mass of GCs in real dwarfs (Sect. 2). In Sect. 4.4, we explore how the results change for different GC masses.

## 3.3. Description of the galaxy models

We were interested in whether GCs can survive in isolated dwarf galaxies without sinking to their centers for 10 Gyr. This time was chosen because this is the typical age of GCs of the Milky Way (Massari et al. 2023). We initiated the simulated galaxies to be mostly gaseous because real galaxies were like this at the epoch of GC formation (Tacconi et al. 2010).

We focused on four MOND models, GDw1e8, GDw2e7, GUDG1e8, and GUDG2e7, and for comparison, we considered two Newtonian models, GDw1e8N and GDw2e7N (see Tables 1 and 2). According to the scaling relations derived in Bílek et al. (2021), the sinking of GCs should be the fastest for the smallest dwarfs. The baryonic masses and radii of the models GDw1e8 and GDw1e8N were thus chosen according to the parameters of the observed dwarfs with a mass of $10^8 M_☉$ because this is the mass above which GCs in dwarfs start to be common (Sect. 2). Similarly, the parameters of the models GDw2e7 and GDw2e7N were chosen to resemble those of the observed dwarfs with a mass of $2 \times 10^7 M_☉$ because this is the mass of some of the least massive galaxies that are known to have GCs. Actually, we could not have set a lower mass for the galaxy because the software we had at our disposal could not generate such galaxies and have them be stable. The initial scale lengths of the stellar and gas disks were the same. Mostly to explore the parameter space, we also considered the MOND models GUDG1e8 and GUDG2e7. They are more spatially extended than their counterparts GDw1e8, GDw2e7, and real isolated dwarf galaxies (Sect. 2). They rather resemble observed isolated gas-rich UDGs (Mancera Piña et al. 2019). These galaxies may host globular clusters (Jones et al. 2023). The main characteristics of all the MOND and Newtonian models are presented in Tables 1 and 2, respectively. Some experimentation and compromises were necessary in order to make the models relatively stable. This is why some of the parameters are not the same for all the models. In the MOND DICE, we found, for example, that the Toomre *Q* parameter has very little influence on the stability of the galaxy. The used computational parameters of the PoR code for the different galaxy models are listed in Table 3.

<footnote n="2">² https://bitbucket.org/SrikanthTN/bonnPoR/src/master/</footnote>
<table><caption>**Table 2.** Parameters of the basic Newtonian models used to generate the initial conditions in PoR and MAGI.</caption><thead><tr><th>Name</th><th><i>M</i><sub>bar</sub><br>[10<sup>8</sup> <i>M</i><sub>☉</sub>]</th><th><i>r</i><sub>D</sub><br>[kpc]</th><th><i>q</i></th><th><i>f</i><sub>gas</sub></th><th><i>T</i><br>[K]</th><th>σ<sub>0</sub><br>[kpc/Myr]</th><th><i>Q</i><sub>R<sub>D</sub></sub></th></tr></thead><tbody><tr><td>GDw1e8N</td><td>1</td><td>0.8</td><td>0.4</td><td>0.9</td><td>23 000</td><td>5</td><td>1</td></tr><tr><td>GDw2e7N</td><td>0.2</td><td>0.5</td><td>0.4</td><td>0.9</td><td>23 000</td><td>0.3</td><td>1</td></tr></tbody></table><table><thead><tr><th>Name</th><th><i>N</i><sub>*</sub><br>[10<sup>6</sup>]</th><th><i>M</i><sub>vir</sub><br>[10<sup>9</sup> <i>M</i><sub>☉</sub>]</th><th><i>r</i><sub>s</sub><br>[kpc]</th><th><i>r</i><sub>trunc</sub><br>[kpc]</th><th><i>w</i><sub>trunc</sub><br>[kpc]</th><th><i>M</i><sub>h</sub><br>[10<sup>9</sup> <i>M</i><sub>☉</sub>]</th><th><i>N</i><sub>h</sub><br>[10<sup>6</sup>]</th></tr></thead><tbody><tr><td>GDw1e8N</td><td>0.1</td><td>21</td><td>4.83</td><td>15</td><td>2</td><td>7.5</td><td>8</td></tr><tr><td>GDw2e7N</td><td>0.1</td><td>6.6</td><td>2.98</td><td>15</td><td>2</td><td>3.3</td><td>8</td></tr></tbody></table>

**Notes.** *M*<sub>bar</sub>: Total baryonic mass. *r*<sub>D</sub>: Disk scale length (the same for the stellar and gas disk). *q*: Axial ratio of the disk. *f*<sub>gas</sub>: Gas fraction. *T*: Gas temperature. *σ*<sub>0</sub> and *Q*<sub>R<sub>D</sub></sub>: Special parameters of the MAGI code (see Miki & Umemura 2018). *N*<sub>s</sub>: Number of particles representing the stellar disk. *M*<sub>vir</sub>: Virial mass of the NFW dark halo (if the halo were not truncated). *r*<sub>s</sub>: Scale radius of the NFW dark halo. *r*<sub>trunc</sub>: Truncation radius of the halo. *w*<sub>trunc</sub>: Width of the truncation. *M*<sub>h</sub>: Actual mass of the halo (after applying the truncation). *N*<sub>h</sub>: Number of particles representing the halo.

<table><caption>**Table 3.** Setup of the PoR code for the different galaxy models.</caption><thead><tr><th>Parameter</th><th>Value (GDw1e8, GDw2e7, GUDG1e8, GUDG2e7, GDw1e8N, GDw2e7N)</th></tr></thead><tbody><tr><td>levelmin</td><td>6, 6, 7, 7, 6, 6</td></tr><tr><td>levelmax</td><td>12, 15, 18, 18, 12, 13</td></tr><tr><td>boxlen</td><td>150 kpc</td></tr><tr><td>mass_sph</td><td>1e-9</td></tr><tr><td>m_refine</td><td>2e3</td></tr></tbody></table>

Current-day simulations of galaxies do not have sufficient resolution to resolve all relevant processes. To model the gas physics, one would have to precisely model processes from the scale of astronomical units (~10<sup>11</sup> m) at which gas collapses into individual stars to the scale of kiloparsecs (~10<sup>19</sup> m) at which galactic tidal fields trigger the collapse of the gas clouds into stars. One thus resorts to approximate recipes for the subgrid physics. They are not only approximate, but they also contain parameters whose values are not known. Therefore, each of our models had several variants with different recipes for star formation and SN feedback in order to explore the possibilities. The variants of the models and their parameters can be found in the tables in Appendix D.

We limited the masses of the newly formed stellar particles to be at least 500 $M_☉$. This roughly corresponds to the minimum masses of OB associations, which can reach up to 10 000 $M_☉$ (Massey & Hunter 1998). This was necessary because RAMSES assumes that a certain fraction of mass of each of the newly formed particles explodes as in a SNe. All stars in the OB association that are destined to explode as SNe are assumed to die immediately after the stellar particle is created because the lifetimes of the very massive stars (5×10<sup>-3</sup> Myr) are shorter than the timestep of the simulation. Their explosion then occurs after the delay of 10 Myr mentioned in Sect. 3.1. If the masses of the new stellar particles in the simulation were too low, the SN feedback would be unrealistically spread in space and time.

Unless stated otherwise, we always started the simulations with GCs at the apocenter and at the distance of one scale length of the galaxy. This roughly agrees with observed positions (Sect. 2). In reality, the GCs could have been formed at larger distances. For the purposes of the current paper, where we are interested in whether the SN feedback can prevent GC sinking, this assumption seems appropriate. Moreover, the fact that the distances of all the observed GCs are on the order of

the stellar scale length of the galaxy indicates that it has always been the case; otherwise, we would encounter a fine-tuning problem. Moreover, as we show below, in the absence of star formation, the sinking times depend dramatically on the mass of the clusters, while only a mild correlation is observed between the masses and galactocentric radii (Fig. 3). This would make the fine-tuning problem even worse.

We modeled the GCs by point masses. We could not model them as resolved objects for the following reason. To model them as such, the simulation would have to have a spatial resolution of a fraction of a parsec so that we could model the internal dynamics of the GCs. At such a high spatial resolution, the mass in each gas cell would be so low that it would not cross our limit of 500 $M_☉$ for the formation of new stellar particle, and no new stellar particles would be formed.

The barycenter of the galaxy was always set as the origin of a right-handed Cartesian coordinate system at the beginning of the simulation. The spin vector of the galaxy pointed in the *z* direction, and the midplane of the disk was in the *x*–*y* plane.

# 4. Results

## 4.1. MOND models without star formation

We first explored whether the GCs can survive in the galaxies without sinking for 10 Gyr if star formation is disabled. In each of the four MOND models, we considered five orbits of the GC: (1) the GC moves on a circular orbit in the plane of the disk, corotating with it and starting at the *x*-axis (prograde orbit hereafter); (2) the same but the GC moves against the direction of the rotation of the disk (retrograde orbit hereafter); (3) the GC is initiated at the *x*-axis and is assigned the circular velocity in the direction of the *z*-axis (polar orbit hereafter); (4) the GC is dropped with a zero velocity with respect to the galaxy center from the *x*-axis (radial orbit hereafter); (5) the GC is dropped with a zero velocity with respect to the galaxy center from the *z* axis (axial orbit hereafter). In all cases, the initial distance of the GC from the galaxy center was about one scale length of the disk (see Table 1 for exact values). The GC in our simulations had a mass of 10<sup>5</sup> $M_☉$. This value is roughly the typical mass of a GC (Sect. 2).

After some time, strong bars formed in most of the models. This time depended on the type of the orbit of the GC. The bars then affected the motions of the GCs substantially. The bar could form either before or after the sinking of the GC was completed.
<figure><img src="image_6.png" /></figure>

**Fig. 4.** Evolution of the distance of the GC from the stellar barycenter of the galaxy for the MOND models without star formation. Each panel corresponds to the galaxy model and each curve to a specific orbit of the GC. The dotted parts of the lines indicate unreliable parts: either the galaxy develops a bar that affects the motion of the GC, or (for GUDG2e7) the GC moves less than three resolution elements from the center.

The results are presented in Fig. 4. The dotted parts of the lines indicate the period when the galaxy has a bar or, for GUDG2e7 (which does not develop a bar), the period when the GC moves less than three resolution elements from the galaxy center. In all models, the GC on the prograde orbit sinks in a few gigayears<fnref n="3" />. In contrast to the isolated gas-free UDGs explored in Bílek et al. (2021), there is no core stalling phase. The curve of the orbital decay of the GCs in the model GUDG2e7 show a slowdown of the decay once the GC approaches the core of the galaxy. We note, however, that during this period, the GC is separated just by about two resolution elements such that the sinking is not modeled properly. If the GC is put on the radial orbit, it sinks in almost the same amount of time, albeit mostly somewhat later. Snapshots from the simulations show that the GC feels the drag of the surrounding gas such that the GC starts corotating with the galaxy. The retrograde orbit is the least prone to orbital decay. In the small dwarf GDw2e7, the sinking of the GCs is the fastest. For all the explored orbits, the sinking takes at most 1 Gyr. In all simulations, images show trailing waves behind the GCs (see the example in Fig. 5).

In Bílek et al. (2021), GCs of UDGs were showing core stalling because of the harmonic potential in the center of the galaxies. Here, we did not observe the stalling. Indeed, the upper panel of Fig. 6 shows that the gravitational potential is not harmonic (the acceleration of a harmonic potential is proportional to the galactocentric radius). This plot was constructed from particles in the simulation GDw2e8_e3 before the onset of SN explo-

sions. Only the particles within 10° from the disk plane were used.

We explored whether the presence of two GCs can affect the survivability of at least one of them. The idea behind this was that one of the GCs would serve as a reservoir of energy for the other. One of the GCs of the model GUDG1e8 (we call it GC1) was placed on the prograde orbit with the initial radius of 2 kpc. The second cluster, GC2, had the same mass and orbital plane but was on a retrograde circular orbit with a radius of 1 kpc. We tried four different initial positions for GC2, offset from each other by ninety degrees. It turned out that GC2 does not affect the sinking time of GC1 appreciably. The reason is that the region in which the gravitational field of a GC dominates over the gravitational field of the galaxy is too small, and therefore the GCs influence each other too rarely to have an effect.

## 4.2. MOND models with star formation

For the simulations without star formation, we found that the orbital decay is usually the strongest if the GC is on the prograde orbit. For the simulations with star formation, we therefore explored only GCs initialized on the prograde orbit to see whether the SN explosion can prevent GCs from sinking even in this least favorable case.

Figures 7–A.3 show how the models evolve<fnref n="4" />. The SN explosions can affect the trajectories of the GC substantially. Depending on how we set the baryonic physics parameters, the GC makes either larger or smaller oscillations in the galaxy.

<footnote n="3">³ A movie (GUDG1e8 prog.mp4) of the GC orbiting on the initially prograde orbit in the model GUDG1e8 can be seen online</footnote>

<footnote n="4">⁴ A movie (GUDG1e8_e1.mp4) of the GC orbiting on the initially prograde orbit in the model GUDG1e8_e1 can be seen online</footnote>
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Map of gas density in the midplane of the disk for the model GDw1e8 without star formation and with a GC on the prograde orbit. There is a density wave trailing behind the cluster. Both the GC and the material of the disk orbit the center of the galaxy counterclockwise.</figcaption></figure>

However, as long as there are SNe exploding (seen as a non-zero star formation rate in the figures), the GC does not sink completely. The SN explosions cause fluctuations in the density of the gas (Fig. 8). These, in turn, give rise to fluctuations in the gravitational potential (bottom panel of Fig. 6). The GC experiences pushes from random directions, which prevent it from staying at any particular place in the galaxy. While the GCs initially move on prograde orbits, the random pushes eventually cause them to make random walks.

However, the SN feedback must have the right strength. If it is too weak, the GC stays very close to the center of the galaxy. This contradicts observations, where we see that GCs are located at around one scale length of the stellar disk (Sect. 2). On the other hand, if the feedback is too strong, the galaxy explodes – the gas is expelled, and since the gravitational field decreases, the stars spread as well. In MOND, for isolated objects, the escape velocity is infinite, so the galaxy cannot dissolve completely. It just becomes large in size and does not follow the observed mass-size relation (Sect. 2). It is possible that some UDGs were formed from dwarf galaxies by such intensive star-bursts (Sect. 5).

The survival of the GC is not possible in the model GDw2e7. This is the model in which the sinking was the fastest in the simulations without star formation. We were not able to find star formation parameters that would lead to sustained star formation. The GC sunk shortly after the SNe had ceased. After the first SNe exploded, gas was heated up and diluted such that it stopped forming any further stars. It did not cool enough even after several gigayears. We suspect that this could be due to insufficient treatment of gas cooling in RAMSES. The code determines the cooling rate of a computational cell only on the basis of its temperature, density, and metallicity. There are no parameters to tune it. However, in a real galaxy, gas has a fractal structure (Elmegreen & Falgarone 1996). There are

<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Effect of the SN bubbles on the profile of the radial acceleration $a_{\text{rad}}$. A positive (negative) radial acceleration indicates that the vector points away from (toward) the galaxy center. At the beginning of the simulation, before SNe appear, the profile is smooth (top). Later, the radial acceleration depends on the position within the galaxy (bottom). The images are 2D histograms of individual stellar particles in the simulation GDw1e8_e3. Only particles near the disk plane are shown here. The upper panel also shows that the gravitational potential of the galaxy is not harmonic (i.e., the acceleration is not $\propto r$).</figcaption></figure>

compact dense cores of the size of a few astronomical units (Pfenniger & Combes 1994), much smaller than the cells in simulations, in which cooling is much more efficient. Therefore it is not sufficient to consider just the average gas properties in the computational cell. In addition, from the old stellar population the explosions of the type Ia SNe occur even without any local gas concentration; however, they were not taken into account in the present simulations. The impact of SNe on the neighboring gas, and therefore the gravitational potential, might be underestimated. In any case, this example underlines the finding that SNe are essential for the survival of GCs in isolated dwarfs.

At first glance, it might appear surprising that SNe affect the trajectory of a GC. The SN forms a bubble that is a spherically symmetric structure. Outside of a sphere enclosing such a bubble, the gravitational field does not depend on whether gas inside the sphere is distributed homogeneously or if it is concentrated toward the surface of the sphere. Similarly, the gravitational force inside this spherical shell is zero. However, one needs to take into account that the galaxy is a disk, and the SN
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Different star formation variants of the model GDw1e8. The values of the star formation parameters are stated in Appendix D. From top to bottom, the panels show: (1) Evolution of the distance of the GC from the stellar barycenter of the galaxy projected to the disk plane, $R_{xy}$, in the units of kiloparsecs. (2) Same but in the units of stellar half-mass radius. (3) Evolution of the stellar half-mass radius, $R_{50}$. (4) Evolution of the star formation rate. In the legend, we indicate just the suffix of the variant of the model indicated in the title of the plot (see the tables in Appendix D).</figcaption></figure>

bubble forms a hollow ring with an increased density near its border. It is known, at least for Newtonian gravity, that the gravitational field inside a ring points outside the center of the ring. When viewing the simulation movies, it really looks as though the GC is pushed away from the centers of the SN bubbles. Moreover, if two or more SN bubbles overlap, as they often do in our simulations, the density disturbances are then not spherically symmetric (Fig. 8).

It is interesting that the galaxies in our simulation have relatively low star formation rates. If the galaxies are to convert 50% of their gas mass into stars in 10 Gyr to reach the current observed gas fraction, then the average star formation rate must be, in decadic logarithm, $-2.3 (-3.0) M_{\odot} \text{yr}^{-1}$ for the $1 \times 10^8 M_{\odot}$ ($2 \times 10^7 M_{\odot}$) galaxy. It turned out that achieving such a high star formation rate is difficult. The models with the highest star formation rates have star formation parameters that are quite different from the fiducial values (Sect. 3.1). While we could manage an intensive starburst right after the star of the simulation, the

resulting SNe diluted the gas such that the subsequent star formation was strongly decreased. We were able to achieve higher star formation rates easily in the Newtonian simulations, as the massive dark halos prevent gas from being spread. Nevertheless, the consequence of this disparity between the observed and simulated star formation indicates that in the real Universe, there are more SNe than in our simulation. Therefore, the real effect of SNe on the motion of GCs and stars must be larger than in the simulations.

### 4.3. Newtonian simulations

We were interested in whether the problem of fast sinking of GCs is encountered in Newtonian gravity with CDM. In the CDM scenario, dwarf galaxies are dominated in mass by dark matter, even in the very center. Although the NFW profile predicts a DM cusp in the center, it is possible that it is mitigated by the star formation feedback, although the weak level
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Map of gas density in the midplane of the disk for the model GDw1e8_e3.</figcaption></figure>

of star formation might not be able to transform the cusp into a core.

In the model GDw1e8N, we put the GC on a radial orbit with an initial distance of 1 kpc. We expected the strongest dynamical friction when the GC goes through the center of the halo, where it has an infinite density (Navarro et al. 1996). The result is shown in Fig. 9 by the green curve. The GC sinks in 5 Gyr. This indicates GCs have the potential to sink in the centers of their hosts even in the ΛCDM cosmology.

Next, we were interested in whether SNe can prevent such a cluster from sinking. The star formation parameters are in Table D.4. The simulations turned out to run much slower than the MOND simulations. This is why we continued the simulation only until ca. 3 Gyr and had to adjust the standard star formation parameters (Sect. 3.1) in order to decrease the star formation rate. The result is shown by the thin red curve in Fig. 9. Unlike in the case without star formation, there is no hint of sinking within the simulation time, meaning that the SNe really prevent the GC from sinking. We then ran another simulation, placing the GC at an initial distance of just 0.4 kpc from the galaxy center while keeping the star formation parameters the same. It again had a zero relative velocity with respect to the galaxy. The result is shown in Fig. 9 by the thick red curve. In this case, the GC sinks in ca. 1.5 Gyr. We repeated the exercise for the model GDw2e7N. The GC started from 0.5 kpc. This time, we used the standard star formation parameters. The GC did not show any hints of sinking for the whole simulation time until ca. 2 Gyr.

We conclude from these simulations that SN explosions can but do not always prevent GCs from sinking in the context of ΛCDM cosmology. The initial position and velocity of the GC is an important factor.

It is worth noting that while we explored just three situations, the GCs approximately kept their initial apocentric distances in two of them. In most MOND simulations, the apocentric distance oscillates either around a larger or lower value. Moreover,

unlike the MOND simulations, the star formation rate remained close to or higher than the value required by observations without any ad hoc tuning.

### 4.4. Different masses of the GCs

We also briefly explored whether our conclusions for the MOND models would change if we changed the mass of the GCs. According to the Chandrasekhar formula for dynamical friction and the results of the MOND simulations of Bílek et al. (2021), the sinking time should decrease when increasing the mass of the GC.

Indeed, when the GC in the model GUDG1e8 was decreased to $10^4 M_{\odot}$ on the prograde orbit, it did not have a problem surviving in the galaxy for 10 Gyr, even without star formation. The galaxy, however, developed a strong bar in the simulation time of about 4 Gyr, which influenced the GC substantially.

We then explored the effect of increasing the mass of the GC to $10^6 M_{\odot}$. This is the maximum we encountered in the sample of the observed isolated dwarf galaxies (Sect. 2), even if only in one galaxy. We entered this GC in the model GDw1e8, first on the prograde orbit, without star formation enabled<fnref n="5" />. The sinking progressed differently than usual. The GC first accumulated a massive cocoon of gas and stars around itself, which effectively made its mass even larger. After finishing a single revolution around the galaxy barycenter, the GC captured the central density peak of the galaxy. This large overdensity kept orbiting the galaxy pericenter without approaching the barycenter much more. If observed, the galaxy would be classified as having lopsided outer parts of the disk, and the GC would be right in the middle of the galaxy center. This process took just about 300 Myr.

If the $10^6 M_{\odot}$ GC is put on a retrograde orbit, the situation is different<fnref n="6" />. Instead of capturing gas, the GC induced a V-shaped wave trailing behind it – like the one that can be seen on the water behind a ship. The sinking was slower than in the prograde case. Before the GC reached the center of the galaxy, the wave formed an overdensity near the galaxy center, orbiting it in the direction of the disk rotation. Once the GC sunk close enough to this overdensity, the GC suddenly switched the direction of its orbit and stuck to the central density peak of the galaxy. Again, observationally, we would classify the object as a lopsided galaxy with a GC in its central density peak. Reaching this state took just about 700 Myr.

We then enabled star formation. The GC was put on the prograde orbit with an initial radius of 1 kpc<fnref n="7" />. The feedback parameters were set as in the simulation GDw1e8_e3 because in that simulation, the GC was kept at large galactocentric distances the easiest when it had a mass of $10^5 M_{\odot}$. In spite of that setup, with the massive GC, the sinking went virtually the same as without star formation. The SNe were too weak to affect the process substantially. The GC harbored in the center of the central density peak in 300 Myr. Exchanging four random stellar particles of the galaxy by GCs of this mass did not change the outcome much. All four GCs and the galaxy density peak merged into a single off-centered object in 540 Myr. Then, back with a single GC in the simulation, we reduced the mass of the GC to $4 \times 10^5 M_{\odot}$. This is more appropriate for a galaxy of this mass (Fig. 1). The GC merged with the central density peak in 500 Myr. When this GC was put on the retrograde orbit, which is usually the orbit

<footnote n="5">⁵ A movie (GDw1e8_GC1e6_program.mp4) can be seen online</footnote>

<footnote n="6">⁶ A movie (GDw1e8_GC1e6_retro.mp4) can be seen online</footnote>

<footnote n="7">⁷ A movie (GDw1e8_e3_GC1e6.mp4) can be seen online</footnote>
<figure><img src="image_11.png" /></figure>

**Fig. 9.** Newtonian models. The top panel shows the evolution of the distance of the GC from the particle barycenter of the galaxy. The bottom panel shows the evolution of the star formation rate.

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

**Fig. 10.** Map of gas density in the midplane of the disk of the model GDw1e8_e3. In this case, the mass of the GC was $4 \times 10^5 M_{\odot}$ and moved on the retrograde orbit. The frame shows the final state of the simulation, when the GC, marked by the dot, settles in the offset central density peak of the galaxy. The stellar barycenter of the galaxy is nearly the center of the coordinate system and is marked by the asterisk.

that is the least prone to sinking, the GC merged with the density peak in 1300 Myr, as shown in Fig. 10. Just as without star formation, 20 Myr before merging, the GC suddenly reversed the direction of its motion. It stuck to the stellar and gas density peaks, being offset with respect to the galaxy stellar barycenter. We also put the GC of the mass of $10^6 M_{\odot}$ on a circular retrograde orbit inclined by 45° with respect to the galactic disk. In this case, the GC captured the central density peak after 540 Myr.

In total, in our MOND simulations, the SNe cannot prevent massive GCs from sinking into the center of the galaxy. More precisely, the central density peak of the galaxy and the GC stick together.

We then explored whether the situation is different in the Newtonian case. A GC with a mass of $4 \times 10^5 M_{\odot}$ was entered in the model GDw1e8N with star formation and initiated to move on a retrograde orbit with a radius of 1 kpc. Due to the slow progress of the simulation, we ran it only for 1800 Myr so that we could compare it with the MOND case. The sinking was quite mild, just to 0.82 kpc. Indeed, unless our simulations are wrong, Newtonian gravity with NFW dark halos can explain the observations of massive GCs in dwarf galaxies more easily than MOND.

**5. Other interesting phenomena seen in the MOND simulations**

The SNe bubbles randomly push not only GCs but also stars, increasing their velocity dispersion. We saw in some of our MOND simulations that strong SN winds can disperse the gas, which holds most of the baryonic mass of the galaxy. Once the gas is expelled from the center of the galaxy, the gravitational field becomes weaker in this region such that the stellar component of the galaxy grows in radius too. This can be seen most prominently in the model GDw1e8_e4, where the stellar half mass radius grew from 1 kpc to 3 kpc. Such a mechanism could be a formation channel of UDGs, at least for their small-mass end. Other formation mechanisms of UDGs proposed in the literature for Newtonian gravity include the scenario of gas outflows associated with star formation in several burst episodes within dwarf-sized halos Di Cintio et al. (2017), which is similar to ours. Also, tidal interactions and ram-pressure stripping have been invoked (Carleton et al. 2019).

In addition, the random pushing of stars by SNe bubbles makes the stellar disks thicker. This phenomenon plays a role in both stars and gas and particularly in the smallest mass galaxies. In agreement with this, stellar and gas disks of dwarf galaxies are relatively thick compared to more massive galaxies (e.g., Roychowdhury et al. 2010, 2013).
# 6. Summary and conclusions

Analytic arguments indicate that GCs of low-mass, low-surface-brightness galaxies experience strong dynamical friction in MOND. It should cause the GC to sink in the centers of their host galaxies on timescales on the order of 1 Gyr. Observations of 10-gigayear-old GCs in such galaxies, however, have raised doubts about the validity of MOND – or the correctness of the analytic calculations. Indeed, the simulations of GCs of isolated gasless UDGs by Bílek et al. (2021) showed that the GCs sink much slower than it initially appeared because of the core-stalling mechanism that was not taken into account in the analytic calculations. In the present work, we focused on GCs in isolated dwarf galaxies and UDGs with MOND (more precisely QUMOND) and Newtonian gravity. Our main aim was to see whether SN explosions can prevent GCs from sinking.

We first explored the observational properties of isolated dwarf galaxies and their GCs on the basis of literature data (we are not aware of any observational studies of GCs of isolated UDGs; see Sect. 2). Our main findings are the following: (1) Isolated dwarfs are usually gas-rich objects. The median gas fraction is 50%. Given that the SNe only directly affect gas, this observation justifies the expectation that SNe can cause fluctuations of the gravitational potential and, as a consequence, dynamically heat the orbits of the GCs. (2) The least massive dwarfs do not contain GCs. The least massive dwarf in our sample has a mass of $10^7 M_\odot$. Only dwarfs with masses over about $10^8 M_\odot$ usually have GCs. This defined the galaxies of interest for our study. (3) The distribution of GCs does not seem to differ much from the distribution of the stars of the host galaxy. This indicates that they roughly formed together and are subject to similar forces; otherwise, we would run into a fine-tuning problem. (4) More massive GCs tend to be located closer to the center of the galaxy, with the possible exception of the most massive GCs (mass over about $5 \times 10^5 M_\odot$), in agreement with the expectation that more massive GCs are more affected by dynamical friction.

In the simulations, we initiated the galaxies to follow roughly the observed mass-size relations. They were initiated to be 90% gaseous in order to mimic the conditions at the birth of the GCs 10 Gyr ago. In most of the simulations we performed, the GC had a mass of $10^5 M_\odot$, which is roughly the typical observed mass. The GC was modeled by a point mass. The GCs always started at the apocenter, at about one scale length of the disk.

We started with MOND simulations without star formation. There were four galaxy models with more and less massive dwarfs and more and less massive UDGs. Regardless of the direction of the initial velocity vector, the GCs sunk in the centers of their host galaxies on the timescale of 1 Gyr. Unlike in the simulations of GCs of UDGs (Bílek et al. 2021), no core stalling occurred. The sinking was usually the fastest for the radial and prograde orbit, and it was the slowest for the retrograde orbit. However, in some cases it was not possible to study the orbital decay in its completeness because the galaxies became unstable and started reorganizing themselves. This process was halted or modified when the SN feedback was included later. Nevertheless, this demonstrated that if a bar occurs in the galaxy, it can push the GC and extend its lifetime.

We then introduced star formation and SN feedback. Because star formation and SN explosions can be modeled only approximately in current simulations, we explored various subgrid parameters. The GC was always initiated on the prograde orbit, which is the most prone to sinking. We indeed found that SNe can prevent the $10^5 M_\odot$ GCs from sinking. The GC can for-

get the direction of its initial orbital motion and make a random walk. However, it was quite difficult, or impossible, to choose the baryonic parameters such that the galaxies form an appreciable number of stars. In the low-mass dwarf, the star formation ceased shortly after the simulations started. In the result, the GC sunk quickly after that because there were no SNe. It is not clear at the moment whether the inefficient star formation indicates an improper numerical treatment of the gas physics in PoR or a failure of the QUMOND gravity.

We then explored a few additional situations less thoroughly. We ran a few Newtonian simulations with cuspy NWF halos. Even in them, without star formation, the GC can sink due to dynamical friction in a few gigayears. Star formation can but does not have to prevent the GC from sinking. It depends on the initial position and velocity of the GC. The star formation rate was reasonable with the default baryonic parameters of the code.

We then explored how our conclusions change if we change the mass of the GC. In the MOND simulations, if it is reduced to $10^4 M_\odot$, then it often does not have difficulties surviving, even without star formation. When we increased the mass of the GC to $4 \times 10^5 M_\odot$ or $10^6 M_\odot$, which are the largest observed GC masses, the situation became the opposite. The SNe were too weak to prevent the GC from sinking, regardless of the initial tangential velocity of the GC. More precisely, after getting close to the galaxy center, the GC attracted and captured the central density peak of the galaxy. In the result, the galaxy would be classified as having an offset center – or lopsided outskirts – and the GC would be classified as a nuclear star cluster. When several massive GCs were put in the simulation, they all ended in the offset central density peak of the galaxy. A single simulation demonstrated that the sinking of massive GCs is by far not as fierce a problem in Newtonian simulations as it is in the MOND ones.

The simulations we have performed therefore favor Newtonian gravity with NFW dark matter halos over the QUMOND gravity. However, the test is not entirely conclusive because hydrodynamics, star formation, and SN feedback are notoriously difficult to model. We identified two potential problems in PoR: (1) The code does not take into account the fact that there can be subgrid compact gas clouds that enhance gas cooling and thus star formation, and (2) PoR currently does not change the prescription for star formation from the Newtonian gravity to QUMOND gravity. While this choice was motivated by the results of Zonoozi et al. (2021), namely, that the free fall times of gas clouds are almost equal in the two types of gravity (Nagesh et al. 2023), we note that Zonoozi et al. (2021) derived their results for gas clouds of specific masses and sizes and not for a general computational cell. Furthermore, it is important to stress that in this work, we have investigated only the QUMOND version of MOND. It is not possible to conclude that our simulations disfavor MOND as a paradigm because the results can be different in other MOND theories. Indeed, MOND theories differing substantially from QUMOND have been discussed recently (Milgrom 2022, 2023b,a). It will be interesting to see whether improvements in feedback modeling will eventually be able to at least rule out QUMOND. Another caveat is that our simulations have not been able to reproduce the core that is often observed in isolated gas-rich dwarfs (Swaters et al. 2009). In such galaxies, their rotation curves suggest, assuming Newtonian gravity, cored dark matter halos with a fraction of a kiloparsec in size (Carignan & Beaulieu 1989; Gentile et al. 2007; Elson et al. 2010). This should be able to stall the GC at 0.3–0.5 kpc in both the CDM and MOND cases. Taken together,
further research is necessary before adopting the CDM hypoth-
esis since multiple phenomena, different from those investigated
here, raise serious doubts about it (Kroupa 2015; Banik & Zhao
2022; Kroupa et al. 2023).

Data availability

Movies are available at https://www.aanda.org

Acknowledgements. We thank Jokić Rosdahl for helpful advice. This work was granted access to the HPC resources of MesoPSL financed by the Region Ile de France and the project EquipMeso (reference ANR-10-EQPX-29-01) of the programme Investissements d'Avenir supervised by the Agence Nationale pour la Recherche. FC acknowledges support from the Programme National de Galaxies et Cosmologie (PNCG), du CNRS/INSU.

References

Banik, I., & Zhao, H. 2022, Symmetry, 14, 1331
Banik, I., O'Ryan, D., & Zhao, H. 2018, MNRAS, 477, 4768
Banik, I., Thies, I., Candlish, G., et al. 2020, ApJ, 905, 135
Banik, I., Nagesh, S. T., Haghi, H., Kroupa, P., & Zhao, H. 2022a, MNRAS, 513, 3541
Banik, I., Thies, I., Truelove, R., et al. 2022b, MNRAS, 513, 129
Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57
Bílek, M., Thies, I., Kroupa, P., & Famaey, B. 2018, A&A, 614, A59
Bílek, M., Samurović, S., & Renaud, F. 2019, A&A, 629, L5
Bílek, M., Zhao, H., Famaey, B., et al. 2021, A&A, 653, A170
Bílek, M., Fensch, J., Ebrová, I., et al. 2022, A&A, 660, A28
Binney, J., & Merrifield, M. 1998, Galactic Astronomy (Princeton, NJ: Princeton University Press)
Carignan, C., & Beaulieu, S. 1989, ApJ, 347, 760
Carleton, T., Errani, R., Cooper, M., et al. 2019, MNRAS, 485, 382
Chae, K.-H., Bernardi, M., & Sheth, R. K. 2018, ApJ, 860, 81
Ciotti, L., & Binney, J. 2004, MNRAS, 351, 285
Combes, F. 2014, A&A, 571, A82
Di Cintio, A., Brook, C. B., Dutton, A. A., et al. 2017, MNRAS, 466, L1
Diemer, B., & Kravtsov, A. V. 2015, ApJ, 799, 108
Dubois, Y., & Teyssier, R. 2008, A&A, 477, 79
Eappen, R., & Kroupa, P. 2024, MNRAS, 528, 4264
Eappen, R., Kroupa, P., Wittenburg, N., Haslbauer, M., & Famaey, B. 2022, MNRAS, 516, 1081
Elmegreen, B. G., & Falgarone, E. 1996, ApJ, 471, 816
Elson, E. C., de Blok, W. J. G., & Kraan-Korteweg, R. C. 2010, MNRAS, 404, 2061
Famaey, B., & Binney, J. 2005, MNRAS, 363, 603
Famaey, B., & McGaugh, S. S. 2012, Liv. Rev. Relat., 15, 10
Freundlich, J., Famaey, B., Oria, P.-A., et al. 2022, A&A, 658, A26
Gentile, G., Salucci, P., Klein, U., & Granato, G. L. 2007, MNRAS, 375, 199
Gentile, G., Famaey, B., & de Blok, W. J. G. 2011, A&A, 527, A76
Georgiev, I. Y., Puzia, T. H., Hilker, M., & Goudfrooij, P. 2009, MNRAS, 392, 879
Hernandez, X., & Gilmore, G. 1998, MNRAS, 297, 517
Hui, L., Ostriker, J. P., Tremaine, S., & Witten, E. 2017, Phys. Rev. D, 95, 043541
Hunter, D. A., Ficut-Vicas, D., Ashley, T., et al. 2012, AJ, 144, 134
Hunter, D. A., Elmegreen, B. G., Goldberger, E., et al. 2021, AJ, 161, 71
Iocco, F., Pato, M., & Bertone, G. 2015, Phys. Rev. D, 92, 084046

Jones, M. G., Karunakaran, A., Bennet, P., et al. 2023, ApJ, 942, L5
Karachentsev, I. D., & Kaisina, E. I. 2019, Astrophys. Bull., 74, 111
Katz, N. 1992, ApJ, 391, 502
Kennicutt, R. C., Jr 1989, ApJ, 344, 685
Kroupa, P. 2015, Can. J. Phys., 93, 169
Kroupa, P., Gjergo, E., Asencio, E., et al. 2023, ArXiv e-prints [arXiv: 2309.11552]
Lüghausen, F., Famaey, B., Kroupa, P., et al. 2013, MNRAS, 432, 2846
Lüghausen, F., Famaey, B., & Kroupa, P. 2015, Can. J. Phys., 93, 232
Mancera Piña, P. E., Fraternali, F., Adams, E. A. K., et al. 2019, ApJ, 883, L33
Massari, D., Aguado-Agelet, F., Monelli, M., et al. 2023, A&A, 680, A20
Massey, P., & Hunter, D. A. 1998, ApJ, 493, 180
Miki, Y., & Umemura, M. 2018, MNRAS, 475, 2269
Milgrom, M. 1983a, ApJ, 270, 371
Milgrom, M. 1983b, ApJ, 270, 384
Milgrom, M. 1983c, ApJ, 270, 365
Milgrom, M. 2010, MNRAS, 403, 886
Milgrom, M. 2022, Phys. Rev. D, 106, 064060
Milgrom, M. 2023a, Phys. Rev. D, 108, 084005
Milgrom, M. 2023b, Phys. Rev. D, 108, 063009
Nagesh, S. T., Banik, I., Thies, I., et al. 2021, Can. J. Phys., 99, 607
Nagesh, S. T., Kroupa, P., Banik, I., et al. 2023, MNRAS, 519, 5128
Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563
Nedkova, K. V., Häußler, B., Marchesini, D., et al. 2021, MNRAS, 506, 928
Nipoti, C., Ciotti, L., Binney, J., & Londrillo, P. 2008, MNRAS, 386, 2194
Pace, A. B., Walker, M. G., Koposov, S. E., et al. 2021, ApJ, 923, 77
Perret, V., Renaud, F., Epinat, B., et al. 2014, A&A, 562, A1
Pfenniger, D., & Combes, F. 1994, A&A, 285, 94
Planck Collaboration XIII. 2016, A&A, 594, A13
Rasera, Y., & Teyssier, R. 2006, A&A, 445, 1
Rejkuba, M. 2012, Ap&SS, 341, 195
Renaud, F., Famaey, B., & Kroupa, P. 2016, MNRAS, 463, 3637
Roszkowski, L., Sessolo, E. M., & Trojanowski, S. 2018, Rep. Prog. Phys., 81, 066201D
Roychowdhury, S., Chengalur, J. N., Begum, A., & Karachentsev, I. D. 2010, MNRAS, 404, L60
Roychowdhury, S., Chengalur, J. N., Karachentsev, I. D., & Kaisina, E. I. 2013, MNRAS, 436, L104
Rubin, V. C., Ford, W. K., Jr, & Thonnard, N. 1980, ApJ, 238, 471
Sánchez-Salcedo, F. J., Reyes-Iturbide, J., & Hernandez, X. 2006, MNRAS, 370, 1829
Schmidt, M. 1959, ApJ, 129, 243
Sharina, M. E., Puzia, T. H., & Makarov, D. I. 2005, A&A, 442, 85
Shi, Y., Helou, G., Yan, L., et al. 2011, ApJ, 733, 87
Skordis, C., & Zlošnik, T. 2021, Phys. Rev. Lett., 127, 161302
Spitler, L. R., & Forbes, D. A. 2009, MNRAS, 392, L1
Swaters, R. A., Sancisi, R., van Albada, T. S., & van der Hulst, J. M. 2009, A&A, 493, 871
Tacconi, L. J., Genzel, R., Neri, R., et al. 2010, Nature, 463, 781
Teyssier, R. 2002, A&A, 385, 337
Teyssier, R., Chapon, D., & Bournaud, F. 2010, ApJ, 720, L149
Teyssier, R., Pontzen, A., Dubois, Y., & Read, J. I. 2013, MNRAS, 429, 3068
Thomas, G. F., Famaey, B., Ibata, R., Lüghausen, F., & Kroupa, P. 2017, A&A, 603, A65
Tiret, O., & Combes, F. 2007, A&A, 464, 517
Tiret, O., & Combes, F. 2008, ASP Conf. Ser., 396, 259
Wittenburg, N., Kroupa, P., & Famaey, B. 2020, ApJ, 890, 173
Wittenburg, N., Kroupa, P., Banik, I., Candlish, G., & Samaras, N. 2023, MNRAS, 523, 453
Zonoozi, A. H., Lieberz, P., Banik, I., Haghi, H., & Kroupa, P. 2021, MNRAS, 506, 5468
Zwicky, F. 1937, ApJ, 86, 217
**Appendix A: Additional figures**

<figure><img src="image_13.png" /><figcaption>**Fig. A.1.** Different star formation variants of the model GDw2e7. (See Fig. 7 for a description.)</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Fig. A.2.** Different star formation variants of the model GUDG1e8. (See Fig. 7 for a description.)</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Fig. A.3.** Different star formation variants of the model GUDG2e7. (See Fig. 7 for a description.)</figcaption></figure>
## Appendix B: Distances to observed galaxies

<table><caption>**Table B.1.** Distances to the observed galaxies according to different sources.</caption><thead><tr><th>Name in LV</th><th>Name in source</th><th>Distance in LV [Mpc]</th><th>Distance in source [Mpc]</th><th>Source</th></tr></thead><tbody><tr><td>KKR25</td><td>KKR25</td><td>1.91</td><td>1.86</td><td>S</td></tr><tr><td>UGC08508</td><td>U8508</td><td>2.67</td><td>2.56</td><td>S</td></tr><tr><td>DDO190</td><td>DDO190</td><td>2.83</td><td>2.79</td><td>S</td></tr><tr><td>ESO379-007</td><td>E379-07</td><td>5.45</td><td>5.22</td><td>S</td></tr><tr><td>ESO321-014</td><td>E321-014</td><td>3.33</td><td>3.19</td><td>S</td></tr><tr><td>KKH5</td><td>KKH5</td><td>5.42</td><td>4.26</td><td>S</td></tr><tr><td>KKH34</td><td>KKH34</td><td>7.28</td><td>4.61</td><td>S</td></tr><tr><td>KKH98</td><td>KKH98</td><td>2.58</td><td>2.45</td><td>S</td></tr><tr><td>KK16</td><td>KK16</td><td>5.62</td><td>4.74</td><td>S</td></tr><tr><td>KK17</td><td>KK17</td><td>5.01</td><td>4.72</td><td>S</td></tr><tr><td>KKH18</td><td>KKH18</td><td>4.79</td><td>4.43</td><td>S</td></tr><tr><td>ESO489-056</td><td>E489-56</td><td>6.31</td><td>4.99</td><td>S</td></tr><tr><td>ESO490-017</td><td>E490-17</td><td>6.34</td><td>4.23</td><td>S</td></tr><tr><td>UGC03755</td><td>U3755</td><td>7.69</td><td>5.22</td><td>S</td></tr><tr><td>KK65</td><td>KK65</td><td>7.98</td><td>4.51</td><td>S</td></tr><tr><td>UGC04115</td><td>U4115</td><td>7.87</td><td>5.49</td><td>S</td></tr><tr><td>KKH86</td><td>KKH86</td><td>2.61</td><td>2.61</td><td>S</td></tr><tr><td>UGC438</td><td>UA438</td><td>2.22</td><td>2.23</td><td>S</td></tr><tr><td>LSBCD634-03</td><td>D634-03</td><td>9.59</td><td>9.46</td><td>G</td></tr><tr><td>DDO052</td><td>DDO52</td><td>9.86</td><td>10.28</td><td>G</td></tr><tr><td>ESO121-020</td><td>ESO121-20</td><td>6.08</td><td>6.05</td><td>G</td></tr><tr><td>HIPASSJ1247-77</td><td>HIPASSJ1247-77</td><td>3.47</td><td>3.16</td><td>G</td></tr><tr><td>HS117</td><td>HS117</td><td>3.96</td><td>3.96</td><td>G</td></tr><tr><td>IC4662</td><td>IC4662</td><td>2.55</td><td>2.44</td><td>G</td></tr><tr><td>KK182</td><td>KK182</td><td>5.94</td><td>5.78</td><td>G</td></tr><tr><td>KK230</td><td>KK230</td><td>2.21</td><td>1.92</td><td>G</td></tr><tr><td>KK246</td><td>KK246</td><td>6.85</td><td>7.83</td><td>G</td></tr><tr><td>UGC07242</td><td>KKH77</td><td>5.45</td><td>5.42</td><td>G</td></tr><tr><td>NGC4605</td><td>NGC4605</td><td>5.55</td><td>5.47</td><td>G</td></tr><tr><td>IC779</td><td>UGC7369</td><td>16.67</td><td>11.60</td><td>G</td></tr></tbody></table>

**Notes.** “LV” stands for the Local Volume database (Karachentsev & Kaisina 2019) (Sect. 2). “Source” refers to the publication from which we took the parameters of the GCs: “S” stands for Sharina et al. (2005) and “G” for Georgiev et al. (2009).

## Appendix C: Stability of the simulated NFW halo

<figure><img src="image_16.png" /><figcaption>**Fig. C.1.** Demonstration that the MAGI code generates stable dark matter halos. Blue curve: the cumulative mass profile of the dark matter particles at the beginning of the simulation GDw1e8N without star formation. Orange curve: the same profile at the end of the simulation, 6300 Myr later.</figcaption></figure>

## Appendix D: Star formation parameters

<table><caption>**Table D.1.** Variants of the model GDw1e8.</caption><thead><tr><th>Variant</th><th>n_star</th><th>t_star</th><th>eta_sn</th><th>f_ek</th></tr></thead><tbody><tr><td>GDw1e8_e1</td><td>0.5</td><td>3.0</td><td>0.1</td><td>0.5</td></tr><tr><td>GDw1e8_e2</td><td>0.1</td><td>3.0</td><td>0.1</td><td>0.5</td></tr><tr><td>GDw1e8_e3</td><td>0.1</td><td>6.0</td><td>0.1</td><td>0.5</td></tr><tr><td>GDw1e8_e4</td><td>0.01</td><td>3.0</td><td>0.1</td><td>0.5</td></tr><tr><td>GDw1e8_e6</td><td>2</td><td>3.0</td><td>0.1</td><td>0.5</td></tr><tr><td>GDw1e8_e7</td><td>0.1</td><td>10.0</td><td>0.1</td><td>0.5</td></tr></tbody></table>

<table><caption>**Table D.2.** Variants of the model GDw2e7.</caption><thead><tr><th>Variant</th><th>n_star</th><th>t_star</th><th>eta_sn</th><th>f_ek</th></tr></thead><tbody><tr><td>GDw2e7_e1</td><td>0.5</td><td>3.0</td><td>0.1</td><td>0.5</td></tr><tr><td>GDw2e7_e2</td><td>0.1</td><td>3.0</td><td>0.1</td><td>0.5</td></tr><tr><td>GDw2e7_e3</td><td>0.1</td><td>2.0</td><td>0.1</td><td>0.5</td></tr><tr><td>GDw2e7_e4</td><td>0.1</td><td>1.0</td><td>0.1</td><td>0.5</td></tr><tr><td>GDw2e7_e5</td><td>0.5</td><td>3.0</td><td>0.1</td><td>1</td></tr><tr><td>GDw2e7_e6</td><td>0.5</td><td>3.0</td><td>0.1</td><td>0.5, Z=0.001</td></tr><tr><td>GDw2e7_e7</td><td>0.5</td><td>1.0</td><td>0.1</td><td>0.5, r<sub>bubble</sub>=50pc</td></tr></tbody></table>

**Notes.** For the GDw2e7_e6 model, the default gas metallicity of 0.1 was decreased to 0.001. For the GDw2e7_e7 model, the SN bubble radius $r_{bubble}$ was changed from the default value of 150 pc to 50 pc.
<table><caption>**Table D.3.** Variants of the models GUDG1e8 and GUDG2e7.</caption><thead><tr><th>Variant</th><th>n_star</th><th>t_star</th><th>eta_sn</th><th>f_ek</th></tr></thead><tbody><tr><td>GUDG1e8_e1</td><td>0.1</td><td>3.0</td><td>0.1</td><td>0.5</td></tr><tr><td>GUDG1e8_e2</td><td>0.1</td><td>0.3</td><td>0.1</td><td>0.5</td></tr><tr><td>GUDG1e8_e3</td><td>0.05</td><td>0.3</td><td>0.1</td><td>0.5</td></tr><tr><td>GUDG1e8_e4</td><td>0.05</td><td>3.0</td><td>0.1</td><td>0.5</td></tr><tr><td>GUDG1e8_e5</td><td>0.1</td><td>6.0</td><td>0.1</td><td>0.5</td></tr><tr><td>GUDG1e8_e6</td><td>0.01</td><td>3.0</td><td>0.1</td><td>0.5</td></tr><tr><td>GUDG1e8_e8</td><td>0.001</td><td>3.0</td><td>0.1</td><td>0.5</td></tr><tr><td>GUDG1e8_e9</td><td>0.01</td><td>3.0</td><td>0.01</td><td>0.5</td></tr><tr><td>GUDG1e8_e10</td><td>0.001</td><td>3.0</td><td>0.01</td><td>0.5</td></tr><tr><td>GUDG1e8_e12</td><td>0.01</td><td>3.0</td><td>0.1</td><td>0.0</td></tr><tr><td>GUDG1e8_e13</td><td>0.01</td><td>3.0</td><td>0.1</td><td>1.0</td></tr></tbody></table>

**Notes.** The parameters are given only for the GUDG1e8 models. For the GUDG2e7 models, we used the same sets of parameters.

<table><caption>**Table D.4.** Star formation and feedback parameters of the Newtonian models.</caption><thead><tr><th>Model</th><th>n_star</th><th>t_star</th><th>eta_sn</th><th>f_ek</th></tr></thead><tbody><tr><td>GDw1e8N</td><td>0.1</td><td>6.0</td><td>0.1</td><td>0.5</td></tr><tr><td>GDw2e7N</td><td>0.1</td><td>3.0</td><td>0.1</td><td>0.5</td></tr></tbody></table>

**Notes.** Only one variant was considered for each model.