© M. Bílek et al. 2021

1. Introduction

Modified Newtonian dynamics, or Milgromian dynamics (MOND; Milgrom 1983a, b) was introduced four decades ago as a possible explanation for the missing mass or gravity in galaxies not relying on an elusive dark matter component in these stellar systems (see Famaey & McGaugh 2012, for a review). Under MOND theory, rather than adding dark matter, the laws of gravity or inertia have to be adjusted in the regime of small acceleration a ≪ a₀, where a₀ is on the order of 10⁻¹⁰ m s⁻² and would constitute a new natural constant. This constant a₀ would represent the transition from Newtonian into the so-called deep-MOND regime and could even be related to the cosmological constant (Milgrom 1999). The main MOND implication under deep-MOND and spherical symmetry is that the true acceleration norm a should become $a=\sqrt{a_{\mathrm{N}}{a}_0}$, where a_N is the norm of the Newtonian gravitational acceleration. As an extension of classical gravity, MOND requires a new Poisson equation derived from a Lagrangian. Two such classical MOND Lagrangians have been proposed (Bekenstein & Milgrom 1984; Milgrom 2010), while several Lorentz covariant extensions have been proposed over the last two decades; the most recent were proposed by Skordis & Złośnik (2019, 2020), which are consistent with the gravitational wave signals.

The MOND theory has celebrated some major achievements on galaxy scales during the last decades. A key prediction of MOND is the tight coupling between the acceleration generated by baryons inside an object and its intrinsic acceleration (Milgrom 1983a). About 40 years after its proposition, this coupling has recently been empirically confirmed by high-quality observations of rotation curves of galaxies with stellar masses spanning several orders of magnitude (Lelli et al. 2016, 2017). This coupling extends into the regime of low surface brightness galaxies, which were not known at the time of the development of MOND. Strong signs of the external field effect, a unique feature of MOND, have been detected (McGaugh & Milgrom 2013a,b; McGaugh 2016; Caldwell et al. 2017; Haghi et al. 2016; Hees et al. 2016; Chae et al. 2020). Simulated galaxy formation under MOND easily produces exponential galactic disks (Wittenburg et al. 2020). Furthermore, MOND may provide a solution (Zhao et al. 2013; Bílek et al. 2018; Banik & Zhao 2018) to 100 kpc scale phenomena such as the planes-of-satellite problem (Pawlowski 2018; Pawlowski & Kroupa 2020; Müller et al. 2021a). It could also for instance help in shedding light on the origin of large density fluctuations on a 600 Mpc scale (Haslbauer et al. 2020) or on the rapid formation of galaxy clusters (Asencio et al. 2021).

While MOND is able to account for many aspects of galactic phenomenology – in particular the detailed shape of the HI rotation curves of spirals, but also X-ray temperature profiles, galaxy-galaxy lensing, and globular cluster (GC) kinematics of most ellipticals (Milgrom 2012, 2013; Samurović 2014; Bílek et al. 2019), this phenomenology can also be understood with appropriate, although fine-tuned, dark matter distributions in the Newtonian context. A profound fundamental difference, however, between MOND and Newtonian dynamics with dark matter halos is the drag experienced by a massive satellite moving through the bath of particles of its host galaxy, which is comprised mostly of stars in the MOND case and stars and dark matter in the Newtonian case. This drag, due to the back-reaction of the wake and angular momentum transfer generated by the massive satellite, is called dynamical friction, and has been not been explored very much in the MOND context. Under MOND, this was pioneered by Ciotti & Binney (2004), who found with analytical arguments that the dynamical friction timescale is shorter in the deep-MOND regime than with Newtonian gravity including an equivalent dark matter distribution. However, later simulations showed this to be valid only in the regime of small perturbations, while on the contrary, dynamical friction in MOND is rather inefficient when the perturber has a large mass (e.g., Nipoti et al. 2008; Combes 2014). Going further than the mere timescale of dynamical friction, Sánchez-Salcedo et al. (2006) proposed a heuristic formula for the drag force, which we set out in this work to test by N-body simulations for the first time.

We focus on the case of GCs moving in and around isolated ultra-diffuse galaxies (UDGs), that is, in galaxies with very low surface brightnesses and very large radial extents (Sandage & Binggeli 1984; van Dokkum et al. 2015a). Whilst this type of galaxies has been identified for a long time, interest in these galaxies has grown again over the last few years owing to their ubiquity in galaxy groups and galaxy clusters (e.g., Crnojević et al. 2014; Koda et al. 2015; Muñoz et al. 2015; van Dokkum et al. 2015b,a; Mihos et al. 2015; Yagi et al. 2016; van der Burg et al. 2017; Venhola et al. 2017; Müller et al. 2018; Zaritsky et al. 2019; Iodice et al. 2020), but also in the field (Román & Trujillo 2017; Leisman et al. 2017). In the standard context, many hypotheses have been put forward regarding their possible origin, such as feedback mechanisms that might have made them lose their gas early or expand through dynamical heating (e.g., Freundlich et al. 2020) or tidal debris from mergers or tidally disrupted dwarfs (e.g., Greco et al. 2018; Toloba et al. 2018; Jiang et al. 2019; Carleton et al. 2019). As extreme systems, they are ideal testbeds for various models of galaxy formation. In the context of MOND, irrespective of their formation scenario, such galaxies in the field are very well suited to the study of dynamical friction because they are deep in the MOND regime all the way down to their center. Of course, this is not the case for galaxies residing in an environment in which the external field should dominate their dynamics. The exercise of studying dynamical friction in isolated UDGs is also interesting from the point of view of testing MOND itself. Analytic arguments suggested that GCs of galaxies in the deep MOND regime experience strong dynamical friction and sink quickly to the centers of the galaxies (Ciotti & Binney 2004; Nipoti et al. 2008). Many UDGs contain many GCs (Lim et al. 2018, 2020), which might even be exceptionally massive (van Dokkum et al. 2018a) and therefore experiencing exceptionally strong dynamical friction. This brings up the question of whether GCs could have survived orbiting in (relatively isolated) UDGs since their formation 10 Gyr ago in MOND, or if MOND predicts that such GCs should have sunk to the centers of UDGs and formed nuclear star clusters. It is thus important to test the analytic formulas by simulations. The analytic formulas also neglect the gravitational interactions between the individual GCs in the galaxy, which might have slowed down the sinking, as is known from simulations with Newtonian gravity.

This paper is organized as follows. In Sect. 2 we review the current knowledge of dynamical friction in MOND in more detail and identify the open questions that we later answer by N-body self-consistent simulations. The technical details of the simulations and the characteristics of the simulated UDG are described in Sect. 3. In Sect. 4 we demonstrate that a single GC does not sink right to the center of the UDG but it experiences core stalling. We further compare the simulation of one GC orbiting the UDG with the prediction of a MOND analog of the Chandrasekhar formula in Sect. 5. The evolution of a whole GC system incorporated into our simulated UDG is studied in Sect. 6, where we find the GC system shrinks but the shrinking does not progress to a zero radius. We then use scaling arguments to show that the same is expected for a vast majority of the structural parameters of observed UDGs (as long as they are not external field dominated) in Sect. 7. In Sect. 8 we investigate analytically how much energy is absorbed internally in GCs during the GC-GC interaction, which is a process that is not present in our simulations with GCs modeled as single particles. We find that that our simulations are a good approximation for usual GC masses but less good for the GC masses such as those observed in the NGC 1052-DF2 galaxy. The results are summarized in Sect. 9.

2. Dynamical friction

We first introduce dynamical friction in the context of Newtonian dynamics. Dynamical friction is a process that transfers the orbital energy and angular momentum of a pair of bodies orbiting each other into the energy and angular momentum of their constituents. For example, when a galaxy is orbited by a satellite, the orbital energy of the satellite is transferred to the orbital energy of the stars or dark matter particles of the galaxy. Dynamical friction causes the bodies to merge. It has two dominant sources (see for instance Mo et al. 2010 or Just & Peñarrubia 2005 for a review). The first is a density wake of the constituent particles behind the object being decelerated. The wake backreacts gravitationally on the orbiting body against the direction of its motion. Supposing Newtonian gravity, the frictional deceleration can in most situations be approximated by the Chandrasekhar formula (Chandrasekhar 1943) as follows:

$$\begin{array}{c}\hfill a_{\mathrm{DF},\mathrm{NWT}}=\frac{2\pi ln\mathrm{\Lambda }{G}^2\rho m}{{\sigma }^2{X}^2}[\mathrm{erf}(X)-\frac{2X}{\sqrt{\pi }}exp(-X^2)]·\end{array}$$(1)

In the above equation, ρ stands for the density of the environment of the decelerated body (the enhancement of the density due to the wake is negligible) and σ the local velocity dispersion of the particles of the environment. Next, m denotes the mass of the decelerated body and

$$\begin{array}{c}\hfill X=\frac{v}{\sqrt{2}\sigma },\end{array}$$(2)

where v is the velocity of the decelerated body with respect to its environment. Finally, lnΛ is the Coulomb logarithm. In principle Λ is defined as Λ = b_max/b_min, where b_min is the impact parameter of the background particle whose trajectory is bent owing to the decelerated body by 90° and b_max the maximum impact parameter that the background particle can have, that is, in principle the extent of the system. The definition of the Coulomb logarithm is obviously ambiguous and simulations show that its value depends somewhat on the particular problem under consideration (see, e.g., Hamilton et al. 2018, for an in-depth criticism of the Coulomb logarithm approach). It is possible to roughly estimate lnΛ from its definition, but if a more precise estimate is needed, a calibration of Eq. (1) by numerical simulations is required. It turns out that if the Coulomb logarithm is chosen correctly, Eq. (1) works very well. The suitable value of lnΛ differs from the default one by a factor of a few. The situations under which the Chandrasekhar formula works well were summarized in Mo et al. (2010). While the exact value of Coulomb logarithm has yet to be found by simulations, the Chandrasekhar formula is still highly useful without its exact knowledge. The Chandrasekhar formula gives us a basic understanding of dynamical friction and allows order-of-magnitude estimates of its importance. Formulas based on the Chandrasekhar formula are used for example as a basis when we look for fitting formulas for galaxy merger rates and timescales in simulations (Lacey & Cole 1993; Jiang et al. 2008, 2014; Solanes et al. 2018).

The other source of dynamical friction is called resonant coupling. It happens because the angular momentum of the decelerated body is transferred to the angular momentum of the particles whose orbits are in resonance with the orbit of the decelerated object. For this reason, a satellite can experience dynamical friction even if it moves outside of the galaxy (Lin & Tremaine 1983; Tremaine & Weinberg 1984). In such a case no density wake is formed behind the body and Eq. (1) incorrectly predicts a zero deceleration. In agreement with this second mechanism of dynamical friction, it was found that the effective value of the Coulomb logarithm increases with the circularity of the orbit (Chan et al. 1997).

In the context of our paper, one more effect is important. It again goes against expectations based on the Chandrasekhar formula. It transpires that the satellite can stop its sinking to the center of the host galaxy at a certain radius before reaching the center of the galaxy. This effect is known as core stalling (Hernandez & Gilmore 1998; Read et al. 2006; Petts et al. 2015). It occurs for galaxies whose gravitational potential is nearly harmonic in the center, that is, having a cored density profile for the case of Newtonian dynamics. In such a potential there are no background particles that would have suitable periods to absorb the energy of the decelerated body. If the density profile is only close to harmonic, the sinking satellite can dynamically heat the central particles and form a core that prevents the sinking again.

Dynamical friction is far less explored in MOND. An important analytic result was obtained by Ciotti & Binney (2004), assuming the aquadratic formulation of MOND (Bekenstein & Milgrom 1984). These authros considered a galaxy governed by MOND gravity and the same galaxy in Newtonian gravity with an extra auxiliary rigid gravitational field. The auxiliary field is such that the total gravitational field of the Newtonian galaxy is as in MOND. They found that the ratio of the dynamical friction timescales of the two systems are as follows:

$$\begin{array}{c}\hfill \frac{t_{\mathrm{DF},\mathrm{MOND}}}{t_{\mathrm{DF},\mathrm{N}+\mathrm{AF}}}=\frac{a_0^2}{\sqrt{2}{a}^2},\end{array}$$(3)

where a denotes the typical value of gravitational acceleration in the system. Making use of this and other related equations, they concluded that dynamical friction is more effective in MOND than in Newtonian gravity with live dark matter halos. In spite of this analytic finding, comparative simulations of major mergers of galaxies show the opposite tendency: the merging of galaxies is much faster in Newtonian simulations (Nipoti et al. 2007; Tiret & Combes 2007a,b; Renaud et al. 2017). Combes (2014) simulated the orbital evolution of massive gas clumps in young galaxies (clumps containing about 20% of the baryonic mass of the galaxies) and found that dynamical friction was stronger with Newtonian gravity than in MOND. A resolution of the inconsistency of the analytic calculations and simulations was proposed by Nipoti et al. (2008). These authors verified with a simulation that the analytic estimates work for a small bar rotating in the center of a galaxy. But realistic bars are much larger and take up a significant fraction of the baryonic mass in the central zone. This means that the reservoir of particles to interact with, assumed infinite in the case of the analytic treatment of Ciotti & Binney (2004), is in reality insufficient to slow down the bar pattern speeds in MOND. In conclusion, Ciotti & Binney (2004) suggested that dynamical friction in MOND is not effective if the mass of the perturbers. such as the massive clumps of Combes (2014), is comparable to the mass of the whole system because in such a case it is difficult for the remaining “background” particles to absorb a large amount of energy.

The MOND results we mentioned so far pertained only to the timescales of dynamical friction. It would be desirable to have a tool to describe the effects dynamical friction in more detail, for example, for predicting the position of the decelerated satellite in time. To this end, Sánchez-Salcedo et al. (2006) proposed an analog to the Chandrasekhar formula. They heuristically suggested that dynamical friction in MOND can be evaluated by multiplying the Newtonian Chandrasekhar formula by the ratio of the dynamical friction timescales given by Eq. (3),

$$\begin{array}{c}\hfill a_{\mathrm{DF},\mathrm{MOND}}=a_{\mathrm{DF},\mathrm{NWT}}\frac{a_0^2}{\sqrt{2}{a}^2}·\end{array}$$(4)

This last equation, which we call the Sánchez-Salcedo formula, has never been proven analytically or verified by simulations. We note that even Eq. (3) has only been verified for small galactic bars and not for other orbital configurations of the interacting objects.

The Sánchez-Salcedo formula has been used for investigating the orbital evolution of the GCs of the Fornax dwarf galaxy (Sánchez-Salcedo et al. 2006; Angus & Diaferio 2009). Analytic estimates predicted a fast sinking of GCs to the centers of dwarf galaxies, where the GCs should merge and form a nuclear star cluster (Ciotti & Binney 2004; Sánchez-Salcedo et al. 2006; Nipoti et al. 2008). The presence of GCs and the absence of a nuclear star cluster in Fornax were initially claimed to contradict MOND, but it was later found that the problem can be solved if we assume that the GCs were formed further away from the galaxy than where they are observed today (Angus & Diaferio 2009). Moreover, the applicability of the analytic formula has never been verified by simulations for the case of GCs orbiting a galaxy.

The discovery of a large number of UDGs (e.g., Sandage & Binggeli 1984; Impey et al. 1988; Crnojević et al. 2014; Koda et al. 2015; van der Burg et al. 2017; Müller et al. 2018; Habas et al. 2020; Iodice et al. 2020) brings up the question of survival of GC systems again. These diffuse objects, with a surface brightness fainter than ∼24 in the g band and an effective radius larger than 1.5 kpc, can contain a large number of GCs (Beasley et al. 2016; van Dokkum et al. 2016, 2018a; Müller et al. 2020). The fact that UDGs are deep-MOND objects (Müller et al. 2019) and the GCs of some of these objects seem to be exceptionally massive (van Dokkum et al. 2018a) indicates that the GC systems of UDGs can be affected strongly by dynamical friction.

Newtonian simulations have already been performed to investigate the evolution of the GC system of the UDG NGC 1052-DF2 (DF2 hereafter) by Dutta Chowdhury et al. (2019). This galaxy has been claimed to have a system of unusually massive GCs (van Dokkum et al. 2018a) if it lies at the distance of 20 Mpc (Danieli et al. 2020), which we assume hereafter (see e.g., Trujillo et al. 2019 or Haghi et al. 2019 for a discussion of the effect of the distance on the properties of DF2). The galaxy hosts ten confirmed GCs, all of which have stellar masses between 10⁵ and 10⁶ M_⊙. They contain a non-negligible fraction of the total stellar mass of the galaxy. The simulations by Dutta Chowdhury et al. (2019) captured several effects that do not readily follow from the Chandrasekhar formula or its modified-gravity analogs – the core stalling and the interactions of GCs among each other (the so-called dynamical buoyancy). Such effects have never been included in the studies of the orbital decay of GCs in MOND.

There are many open questions regarding dynamical friction in MOND. In this paper, we wanted to investigate several of them using self-consistent N-body simulations. First, we wanted to explore the precision of the Sánchez-Salcedo formula by studying the motion of a single GC in a UDG. We wanted to learn whether MOND really predicts that the GC will sink in the center of the galaxy over a few gigayears, as following the Sánchez-Salcedo formula. Then we were interested whether the answer will change if we include in the simulation several GCs that interact with each other. The arising dynamical buoyancy could possibly slow down the sinking. These simulation would eventually address whether MOND predicts that GCs of UDGs should have merged and formed nuclear stars clusters in a few gigayears after the formation of the GCs, if the GCs were initially distributed around the galaxy at similar distances as observed today.

Our simulations were inspired by the DF2 galaxy and its possibly overmassive GC system. It is a typical UDG in terms of stellar mass, size and Sérsic index (Habas et al. 2020). We stress that our simulations are not intended to describe DF2 itself because in MOND this galaxy is subjected to a strong external field effect (Famaey et al. 2018; Kroupa et al. 2018; Haghi et al. 2019). Thus, MOND predicts DF2 to behave in the Newtonian way. Our present investigations pertain to the study of DF2 analogs that would be far from a massive neighbor and hence fully unaffected by the external field effect.

3. Setup of the simulations

In this work, we assumed the QUMOND formulation of MOND (Milgrom 2010). If the distribution of matter ρ is given, then the gravitational field can be obtained by solving the QUMOND equivalent of the Poisson equation:

$$\begin{array}{c}\hfill \mathrm{\Delta }\varphi =\mathbf{\nabla }·[\nu (|\mathbf{\nabla }{\varphi }_{\mathrm{N}}|/a_0)\mathbf{\nabla }{\varphi }_{\mathrm{N}}],\end{array}$$(5)

where

$$\begin{array}{c}\hfill \mathrm{\Delta }{\varphi }_{\mathrm{N}}=4\pi G\rho \end{array}$$(6)

is the standard Poisson equation for the Newtonian potential ϕ_N. We employed the so-called simple form of the interpolation function ν,

$$\begin{array}{c}\hfill \nu (x)=\frac{1+{(1+4x^{-1})}^{1/2}}{2},\end{array}$$(7)

and the value of the MOND acceleration constant of a₀ = 1.12 × 10⁻¹⁰ m s⁻², which is consistent with the measurements of McGaugh et al. (2016).

The simulations were performed using the Phantom of RAMSES code (Lüghausen et al. 2015), which is a QUMOND adaptation of the popular adaptive-mesh-refinement code RAMSES (Teyssier 2002). For the simulation described below, unless specified otherwise we used the computational parameters summarized in Table 1. We simulated the galaxy with 10⁷ particles of 20 M_⊙. Each GC was represented by a point mass. For reference, the overmassive GCs of DF2 were estimated to have effective radii of around 8 pc (van Dokkum et al. 2018a). The grid was refined if the mass in the cell exceeded 200 M_⊙, that is, if it contained at least ten stellar particles or one GC. The maximal spatial resolution was 49 pc.

We initialized the density of the simulated galaxy as a Sérsic sphere that has the parameters observed for DF2, that is, a stellar mass of 2 × 10⁸ M_⊙, effective radius of 2 kpc, and Sérsic index 0.6 (van Dokkum et al. 2018b). The galaxy was modeled with 10⁷ particles, each of a mass of 20 M_⊙. Each particle was given a velocity drawn from a three-dimensional isotropic Gaussian distribution. The local velocity dispersion was determined by solving the spherically symmetric Jeans equation. Before introducing any GCs into the simulation, the galaxy was let to virialize for 10 Gyr (over 70 free-fall times from the distance of two effective radii to the center of the galaxy). After this period, a Sérsic fit of the density yields an effective radius of 1.9 kpc and a Sérsic index of 0.7. The radial profiles of density, velocity dispersion, and gravitational acceleration after the virialization are plotted in Fig. 1. The inner ∼7 kpc of the galaxy became dominated by stars on tangential orbits, while the more distant regions are inhabited by stars preferentially on radial orbits.

Fig. 1. Properties of the simulated galaxy just before entering GCs. Top: radial profile of density. Middle: radial profile of velocity dispersion, separately for the full velocity dispersion (full line), velocity dispersion in the radial direction (dotted line), and in the tangential direction (dashed line). Bottom: radial profile of gravitational acceleration in the units of the MOND constant a0.

4. Core stalling: Demonstration on the fiducial model

We begin presenting our results with the simulation that we hereafter refer to as the fiducial. It starts with one GC placed at the distance of 5 kpc from the center of the galaxy such that the two objects are at rest with respect to each other. The GC is modeled as a single particle with a mass of 10⁶ M_⊙, mimicking some of the most massive GCs of DF2. Figure 2 shows the time evolution of the distance of the GC from the center of the galaxy. We can note that the apocentric distances initially decrease at an approximately constant rate, but at a time of about 3 Gyr the orbital decay slows down dramatically. This effect has also been noted in the Newtonian simulations of DF2 of Dutta Chowdhury et al. (2019), who attribute the effect to the core-stalling mechanism.

Fig. 2. Orbital decay of the GC in the fiducial model. The decay slows down substantially once the apocentric distance of the GC reaches about 0.5−1 kpc.

Nevertheless, our simulation is not Newtonian. It is important to verify that the observed core stalling is not due to numerical errors. We can immediately check that the apocentric radius of the GC at the beginning of the stalling phase, at about 0.5 kpc, is an order of magnitude higher than the maximum spatial resolution of the simulation, 0.049 kpc. In addition, we run extra simulations with a higher and lower spatial resolution to demonstrate the convergence of the simulation for our default parameters. The result is shown in Fig. 3. The orbits of the GCs in the fiducial model and that with quadruple resolution nearly coincide. The agreement is not that tight when comparing the default and half resolution simulations, but they still agree well on the apocentric distance of the GC in the stalling phase. This demonstrates that the default resolution is a good balance between precision and computing demands. In total, we found no hints that the observed core stalling would be caused by an insufficient resolution of the simulation.

Fig. 3. Demonstration of the effect of changing of the spatial resolution in the fiducial model on the orbital decay of the GC. The core stalling phase is present for all three choices of resolution.

5. Testing the Sánchez-Salcedo formula for dynamical friction in MOND

The Sánchez-Salcedo formula (Eq. (4)) is supposed to be applicable to objects in the deep-MOND regime (accelerations lower than a₀), which is the case of our UDG (see Fig. 1). We investigated whether the orbital decay of a GC can be modeled analytically by solving the following equation of motion:

$$\begin{array}{c}\hfill \mathit{a}={\mathit{a}}_{\mathrm{grav}}(|\mathit{r}|)+{\mathit{a}}_{\mathrm{DF}}(|\mathit{r}|,|\mathit{v}|;\mathrm{\Lambda }),\end{array}$$(8)

for a suitable value of the Coulomb logarithm Λ, where a_grav is the gravitational acceleration given by Eq. (9) and a_DF the acceleration caused by dynamical friction according to the Sánchez-Salcedo formula. We assumed that the mass of the GC is negligible compared to the mass of the galaxy, such that we can treat the GC as a test particle. A similar approach was used in Sánchez-Salcedo et al. (2006) or Angus & Diaferio (2009). Solving such an equation is of course much faster than running a simulation.

When solving Eq. (8), we made use of the fact that the radial acceleration of a test particle in the gravitational field of a spherically symmetric object in QUMOND can be evaluated as

$$\begin{array}{c}\hfill a_{\mathrm{grav}}=a_{\mathrm{grav},\mathrm{N}}\times \nu \left(\frac{\left|a_{\mathrm{grav},\mathrm{N}}\right|}{a_0}\right),\end{array}$$(9)

where a_grav, N is the radial acceleration expected from Newtonian dynamics (Milgrom 2010). In order to evaluate a_grav, N, we modeled the mass profile of the galaxy analytically as a Sérsic sphere with the post virialization parameters of our stimulated galaxy described in Sect. 3. The density of the Sérsic sphere was calculated by exploiting the analytic approximations by Lima Neto et al. (1999), updated in Márquez et al. (2000). When solving numerically the equation of motion, we also added 1 pc to the distance of the GC from the center of the UDG to avoid numerical problems close to the center of the galaxy.

5.1. Varying the mass of the GC

First, we explored the role of the mass of the GC, keeping the initial orbital configuration fixed. At the beginning of each simulation, we again placed the GC at the distance of 5 kpc from the center of the UDG such that the two objects had a zero relative velocity. We only varied the mass of the GC in each simulation. The resulting evolution of distance of the GCs from the galaxy center in each simulation are shown in Fig. 4 with the blue lines. The mass of the GC in each simulation is indicated next to the corresponding panel of the figure.

Fig. 4. Testing the precision of the Sánchez-Salcedo formula for dynamical friction in MOND (Eq. (4)). Blue lines: Evolution of the GC-galaxy distances in N-body simulations. Orange lines: The same but calculated analytically employing the Sánchez-Salcedo formula. Each panel corresponds to another mass of the GC under consideration. The GCs in all simulations or calculations started with the same initial conditions, that is, starting from rest with respect to the galaxy and at a distance of 5 kpc from the galaxy center. To the right of each figure, the corresponding mass of the GC and the value of the Coulomb logarithm that is needed to be entered in the Sánchez-Salcedo formula to mimic the result of the simulation are stated.

The simulations were compared to the analytic solutions of the model (Eq. (8)). They are indicated in Fig. 4 by the orange lines. The values of the Coulomb logarithm were fitted manually to obtain a balance between matching the instants of apocenters and matching of the overall slopes of the curves. The fitted values are again shown next to the corresponding panels of the figure. All realistic values for the mass of the GC require a value of lnΛ ≃ 3 for these radial orbits. This value has already been proposed by Sánchez-Salcedo et al. (2006).

5.2. Varying the orbital eccentricity of the GC

We also explored how the Sánchez-Salcedo formula works for non-radial orbits. The mass of the GC was fixed at 10⁶ M_⊙. The GC was always initiated as moving in the tangential direction with respect to the center of the galaxy. We specified the orbital initial conditions by the initial distance of the two objects and their relative tangential velocity in units of the circular velocity at the given radius evaluated from Eq. (9).

The results are presented in Fig. 5, where each panel represents one simulation. The text lines on the right of each panel give the initial tangential velocity, initial distance of the GC from the galaxy, and the fitted value of the Coulomb logarithm, respectively. The main conclusion of these experiments is the finding that the effective value of the Coulomb logarithm is higher than in the case of the radial orbits as it ranges between the values of 6.6 and 15. In the case of the GC moving on a circular orbit with a radius of 10 kpc = 5.2 R_e, there is no clear hint of orbital decay in the simulation because of the low density of background particles. The analytic model gives good agreement with this simulation for any reasonable value of the Coulomb logarithm; we tested the range between 0 and 50.

Fig. 5. As in Fig. 4, but varying orbital parameters of the GC. The GC always had the mass of 106 M⊙. The first lines of the notes on the right of the panels indicate the initial tangential velocity of the GC with respect to the galaxy (the radial velocity was zero). It is given in the units of the local circular velocity. The second row indicates the initial galactocentric distance and the third the value of the Coulomb logarithm in the Sánchez-Salcedo formula that provides the best match of the simulation and analytic calculation.

6. Evolution of the density profile of the GC system of a UDG

The properties of the GC systems of UDGs are still a matter of debate. We considered here two types of GC systems. In one, the GCs had the high masses of the GCs of DF2. In the other the masses of GCs were consistent with the standard globular cluster luminosity function (GCLF) of dwarf galaxies. The standard GCLF mass function seems to be more common for observed UDGs at least in the Fornax and Virgo clusters (Prole et al. 2019; Lim et al. 2020). We generally obtained different results for these two types of GC systems.

6.1. DF2-like GC mass function

We added into the galaxy 10 GCs that have masses as inferred for DF2 (Dutta Chowdhury et al. 2019), see Table 2. Their spatial distribution was inspired by the distribution of the GCs of DF2. Dutta Chowdhury et al. (2019) found that the GC system of this galaxy can be described by a Sérsic sphere of the index 1 whose effective radius is 1.3 times the effective radius of the galaxy.

The initial conditions for the GC system were prepared in the following way. We initially drew the positions of the GCs randomly from a distribution matching the observed Sérsic profile. The GCs were first assigned velocities in the same way as we did for the initial velocities of the particles of the galaxy (see Sect. 3). It was then necessary to let the initial conditions for the GCs virialize while avoiding GC-GC interactions and the dynamical friction on the particles of the galaxy. To achieve this, we extracted the virialized mass profile of the galaxy (see Sect. 3), calculated its spherically symmetric gravitational potential, and integrated the equation of motion of each GC individually in this rigid potential for 10 Gyr. The resulting positions and velocities of the GCs were then used as the initial conditions for the GCs in the self-consistent simulation. Contrary to Dutta Chowdhury et al. (2019), we did not strive to assign the GCs the observed radial velocities or the observed projected galactocentric distances because we aimed to investigate only a general UDG similar to DF2.

We ran ten simulations, each for a different random realization of the initial conditions for the GCs. The time evolution of galactocentric distances of the GCs in one of those simulations is plotted in Fig. 6. This plot shows that the motion of individual GCs are influenced by the others, particularly those near the center of the galaxy. The positions of the GCs were recorded in the simulations approximately every 20 Myr. The recorded galactocentric distances of the GCs in all simulations are plotted together in Fig. 7. This figure shows how the radial profile of the GC density evolves on average. The GCs tend to gather in the central 1 kpc of the galaxy after a few gigayears. The profile of GC density gets initially depleted at around a radius of 3 kpc (1.5 R_e of the galaxy), but later the depleted region starts extending toward lower and higher radii. In order to get more quantitative results, we inspected the half-number radius of the GCs. We divided the sample of GCs in halves according to the masses of GCs. The full line in the figure indicates the half-number radius of the high-mass group of GCs and the dashed line is the same for the low-mass group. The figure shows that the high-mass GCs settled faster because of mass segregation. The GCs did not sink right to the very center, likely again because of the core stalling mechanism. The horizontal dashed line in the figure indicates the highest spatial resolution of the simulation. It is again substantially smaller than the radius where the stalling occurs. It can be noted that the core stalling phase establishes after about half of the typical age of a GC. The fact that the spatial distribution of GCs of the DF2 galaxy does not have a core suggests that it has spent the majority of its life under the influence of a strong EFE imposed by its massive neighbor NGC 1052.

Fig. 6. Example of evolution of galactocentric radii of GCs in a simulation with 10 GCs. The masses of the GCs are set as observed for the GCs of the DF2 galaxy (first column of Table 2).

Fig. 7. Evolution of the distances of GCs from the center of the UDG. The black points show combined data from 10 random realizations of a UDG orbited by 10 GCs. The masses of the GCs in each simulation mimic the measured masses of the GCs of the DF2 galaxy (see the first column of Table 2). The horizontal black dashed line indicates the maximum resolution of the simulation. The red lines indicate the median galactocentric distance of the GCs at the given time. This was derived separately for the high-mass half of the GCs (full red line) and for the low-mass half of the GCs (dashed red line).

We divided this dataset into 1 Gyr time intervals. The data points in each interval were converted to radial profiles of GC number density, surface number density, and cumulative fraction of GCs inside a given radius. The results are shown in the left column of Fig. 8. Each line in these plots corresponds to a given time interval. The centers of the time intervals are indicated in the legend according to the color key. These plots can be compared to observations.

Fig. 8. Evolution of GC systems. Left column: simulations with 10 GCs that have the masses of those observed in DF2 (left column of Table 2). Simulations of 10 random realizations of the system were combined in these plots. Galactocentric distance is given in units of the effective radius of the galaxy (Re = 1.9 kpc). Right column: same but for GC masses set according to the standard GCLF (right column of Table 2). Top row: profiles of number density of GCs as a function of galactocentric distance. The individual lines correspond to different times since the start of the simulation. Middle row: evolution of the profile of the surface number density of GCs. Bottom row: evolution of the radial profile of the cumulative fraction of GCs.

In order to compare a galaxy to our simulation with a single parameter, we propose hereafter an observational concentration parameter (not yet used elsewhere to the best of our knowledge) defined as

$$\begin{array}{c}\hfill R=\frac{N(<1R_{\mathrm{e}})-N(<0.5R_{\mathrm{e}})}{N(<1R_{\mathrm{e}})},\end{array}$$(10)

where the symbol N(<x) denotes the number of GCs projected inside a circle of radius x. The time evolution of the R parameter in our simulations is shown in Fig. 9. For example, for the real DF2 galaxy, we have N(<1 R_e) = 6 and N(<0.5 R_e)  =  2 (Dutta Chowdhury et al. 2019) and therefore R = 0.67, such that this galaxy does not show any hint of a central concentration.

Fig. 9. Time evolution of the concentration parameter R of the GC system defined by Eq. (10). The two lines correspond to the two considered distributions of the masses of GCs (see Table 2).

The distribution of velocities of GCs also evolves. Figure 10 shows that the line-of-sight velocity dispersion of the whole GC system decreases by about one-third during the simulated period of 10 Gyr. The line-of-sight velocity dispersion measured from the velocities of the five GCs projected the closest to the center, decreases even more, dropping from about 23 km s⁻¹ at the beginning of the simulation to 8 km s⁻¹ at its end. It is interesting that the velocity dispersion of this inner part of the GC system stabilizes after the first 6 Gyr of the simulation, likely due to the core stalling. The radial profiles of the line-of-sight velocity dispersion of the GC system at the beginning and end of the simulation are compared in the top panel of Fig. 11. This shows that the drop of velocity dispersion happened mostly in the inner part of the GC system. The same plot also represents the line-of-sight velocity dispersion profiles of the stars of the galaxy at the beginning and end of the simulation. In contrast with the GCs, the central velocity dispersion of the stars increased with time. There was also a mild increase in velocity dispersion of the stars in the outer parts of the galaxy. We attributed these changes to energy equipartition of stars and GCs. The radial profiles of the line-of-sight velocity dispersions of stars and GCs are not the same even at the beginning of the simulation. This is because of a different density profile of the stars and GCs. A different profile of the anisotropy parameter can also play a role.

Fig. 10. Line-of-sight velocity dispersions of the GC systems in our simulations as functions of time.

Fig. 11. Radial profile of the velocity dispersion profile of the stars and of the GC system of the simulated galaxy. The curves correspond to the beginning and end of the simulations, as indicated in the legend. The vertical lines indicate the 1σ uncertainty limits.

6.2. Standard GC mass function

We repeated the same procedure as above but substituting the masses of the GCs with values according to a GCLF typical for dwarf galaxies (Durrell et al. 1996; Georgiev et al. 2009). Other than the GC masses, the initial conditions of the GCs remained the same. The distribution of absolute magnitudes of GCs in a galaxy can generally be described by a Gaussian distribution. It turns out that the mean of the distribution is common for nearly all galaxies. Therefore the peak of the GC luminosity distribution is being used as a distance indicator (Rejkuba 2012). The width of the Gaussian distribution depends on the mass of the galaxy (Jordán et al. 2007; Villegas et al. 2010). For UDG galaxies in the Virgo Cluster, Lim et al. (2020) found that their g-band GCLF can be described by a Gaussian with a mean of μ = −7.5 mag and a width of σ = 1.0 mag, which is consistent with similar dwarf galaxies in the same galaxy cluster.

We selected the absolute magnitudes of the GCs such that they follow the mentioned Gaussian distribution, written as $M_k=\sqrt{2}\sigma {\mathrm{erf}}^{-1}[2(\frac{1}{2n}+\frac{k}{n})-1]+\mu$ for k = 0, 1, …, n − 1. These magnitudes were then converted to luminosities adopting the g-band absolute magnitude of the Sun of 5.12 (Sparke & Gallagher 2000). Finally the luminosities were converted to masses by multiplying by 2.2, following Spitler & Forbes (2009), who derived this as the typical value of the mass-to-light ratio of GCs in the V band, a photometric band similar to g.

Similarly to the previous section, we plotted in Fig. 12 the positions of all GCs for every time step for all ten simulations with this more typical GC mass function. The GC system becomes more concentrated with time again and we can again note a mass segregation. The evolution of the GC system is nevertheless notably slower than with the DF2-like GC mass function. The profiles of the GC system at different times after the beginning of the simulation are shown in the right column of Fig. 8.

Fig. 12. Same as in Fig. 7 but with GC masses set according to the standard GCLF.

Figure 10 shows that the line-of-sight velocity dispersion of the GC system again decreases with time, but the decrease is weaker than in the simulation with the GCs with the DF2-like masses. As Fig. 11 shows, there is nevertheless still a strong drop of the central velocity dispersion of the GC system. The central velocity dispersion of the stars of the galaxy again increased but less prominently than in the simulation with the massive GCs. It can be noted from Figs. 10 and 11 that our simulations predict that the observed line-of-sight velocity dispersion of a GC system of a UDG should be lower than that of the stars of the galaxy. There might be exceptions from this rule caused by the low number of GCs.

7. Extension of the results to other UDGs

Finally, we explore the scaling relations of the evolution of GC systems. This is useful for applying the results of our simulations to galaxies that have properties different from those of the galaxy in our simulations. We did this via a dimensional analysis. We were interested in the scaling of the cumulative fraction of GCs f within a galactocentric radius r at time t after the creation of the GC system. We had to assume that GC systems are formed with properties similar to the initial conditions of our simulations. Apart from r and t, f has to depend on the mass of the galaxy, M, its effective radius, R_e, the typical mass of the GC in the GC system, m, and on the constants G and a₀. For simplicity we assumed that all GCs in the GC system have the same mass, galaxies and the GC systems always have the Sérsic indices as in our simulation, the GC systems and galaxies always have a fixed ratio of effective radii, GC systems always contain the same number of GCs, and that the whole system is in the deep-MOND regime. The last assumption implies that in the expression for f, G and a₀ and all masses M_i can appear only together in the products of the form of Ga₀M_i (Milgrom 2001, 2008, 2009). Dimensional analysis then implies that f has to have the form of

$$\begin{array}{c}\hfill f=f(\frac{m}{M},\frac{r}{R_{\mathrm{e}}},\frac{t\sqrt[4]{GM{a}_0}}{R_{\mathrm{e}}})·\end{array}$$(11)

The cumulative fraction of GCs within radius r is invariant for all configurations of the problem that keep the arguments of the function on the right-hand side invariant.

We now illustrate the use of this formula on an example. Let us assume that we want to compare our simulation to an observed galaxy that has a mass of M_o = 0.7 × 10¹⁰ M_⊙, an effective radius of R_e, o = 3 kpc, a Sérsic index of one and whose GCs follow the GCLF, such that their average mass is m_o = 2.6 × 10⁵ M_⊙, and they are t_o = 10 Gyr old. Let the symbols indexed by “s” denote the analogous quantities in the simulation. Equation (11) indicates that we have to compare the observed galaxy to our simulation with the DF2-like GC mass function because the first argument of the function in Eq. (11) is invariant only if

$$\begin{array}{c}\hfill m_{\mathrm{s}}=m_{\mathrm{o}}\frac{M_{\mathrm{s}}}{M_{\mathrm{o}}},\end{array}$$(12)

which equals 7.4 × 10⁵ M_⊙. The third argument of the function in Eq. (11) is invariant if

$$\begin{array}{c}\hfill t_{\mathrm{s}}=t_{\mathrm{o}}\sqrt[4]{\frac{M_{\mathrm{o}}}{M_{\mathrm{s}}}}\frac{R_{\mathrm{e},\mathrm{s}}}{R_{\mathrm{e},\mathrm{o}}},\end{array}$$(13)

or 5.8 Gyr in our case, which is the simulation time to be compared to the observation.

We now applied this scaling relation to investigate the survivability of the GC systems of all UDGs. In Fig. 13 we plotted the time scaling factor from Eq. (13), that is, the product of the last two terms, as a function of the absolute V-band magnitude of the galaxy and of its effective radius. We assumed the mass-to-light ratio of the galaxy to be two. The displayed range of the parameters on the axes of this plot covers most of the observed UDGs. If a UDG has a standard GCLF mass function, then our simulations, in combination with the scaling relations, show that the GC system does not form a nuclear star cluster if the parameters of the UDG lie above the white dashed line in Fig. 13 and to the right of the black dashed line. Nearly all UDGs satisfy this condition; see the data in Muñoz et al. (2015), van Dokkum et al. (2015b), Román & Trujillo (2017), Habas et al. (2020). Thus we deduce, that for the majority of UDGs, GC systems do not merge into a nuclear star cluster in MOND.

Fig. 13. Time scaling factor given by Eq. (13) as a function of the absolute magnitude and effective radius of a galaxy. The white point indicates our simulated galaxy. The white dashed line indicates where the scaling factor equals 1. The black dashed line denotes where the GC-mass scaling factor, the second factor in Eq. (12), equals 2.66, which is the value where an observed UDG with GCs following the standard GCLF mass function has to be compared to our simulations with the DF2-like GC mass function. Supposing that the galaxy has a GC system following the standard GC mass function, our results imply that the galaxies to the right of the black line and to the top the white line do not create a nuclear star cluster by inspiralling of GCs.

8. Influence of the internal structure of the GCs

We simulated the GCs by point masses. The real GCs consist of many stars. Therefore, we might miss some important phenomena. For example, GCs can merge or evaporate because of their mutual encounters. The GC system could potentially also sink faster because its orbital energy would be transferred into the internal energy of GCs during the encounters. Similar questions were already studied for the DF2 galaxy in the Newtonian context by Dutta Chowdhury et al. (2020) with simulations and analytically. They found that, in this context, hardly any GC mergers can be expected in 10 Gyr. Similarly, the internal structure of GCs evolved very weakly. In this section, we discuss the situation for GCs of an isolated UDG in MOND. In this context, because of the creation of the central dense core of the GC system, the interactions between the clusters are stronger and more frequent than in the Newtonian case. It would be ideal to investigate the questions above by simulations with resolved GC. Since this is highly computationally demanding, we progress in this initial study by coupling the many GC simulations described in Sect. 6 with analytic estimates.

8.1. Destruction of GCs by encounters

We investigated the possibility that the evaporation of GCs is caused by tidal shocks during GC-GC encounters. We took all pairs of GCs in every simulation and searched for the times of the relative pericenters of the GC pairs. At every pericenter, we determined the relative distance and tangential velocity of the GCs. We took these as the input for calculating the energy gain of each of the GCs through the impulse approximation. The total gain of energy of a particular GC by encounters ΔE_enc was obtained by summing over all pericenters with respect to all other GCs in the system during the course of the whole simulation. We approximated the GCs as Newtonian Plummer spheres because most GCs are high-acceleration objects.

In particular, the gain of the internal energy of a given GC in one encounter, ΔE_1, enc, was calculated using the formula for the impulse approximation provided by Banik & van den Bosch (2021). Unlike the classical formula (Spitzer 1958; Mo et al. 2010), it takes into account the actual density profiles of the interacting bodies. It is therefore suitable even for the close encounters of GCs that occurred in our simulations. The impulse approximation, nevertheless, is designed only for fast encounters, the so-called tidal shocks. This means that we can consider the stars of the investigated clusters stationary with respect to the cluster center during the characteristic duration of the encounter. In our simulation, this is often not satisfied because the velocity dispersion of the GC system (Fig. 10) is not much higher than the velocity dispersion of the stars inside the GCs, which is around 11 km s⁻¹ (Dutta Chowdhury et al. 2020). When the encounter is very slow (when there is adiabatic variation of the gravitational potential), the energy gain of the cluster after the encounter is zero. This effect can be accounted for by multiplying the energy gain of a GC predicted by the impulse approximation formula by the “adiabatic correction” factor A = [1+(ω/τ)²]^−γ (Gnedin & Ostriker 1999). The typical angular velocity of a star in the cluster ω was evaluated as the ratio of the velocity dispersion at the half-mass radius of the cluster, r_h and r_h itself. The variable τ stands for the timescale of the encounter, which was evaluated as the pericentric distance of the encounter divided by the pericentric velocity. Following Banik & van den Bosch (2021), we put γ = 2−0.5 erf[(τ−2.5 t_dyn)/(0.7 t_dyn)], where $t_{\mathrm{dyn}}=\sqrt{{\pi }^2{r}_{\mathrm{h}}^3/(2GM)}$ is the dynamical time at the half-mass radius of the cluster of the mass M. In order to decide whether the gain of energy of a GC by the encounters can cause a substantial evaporation of GCs or a creation of tidal tails, we compared ΔE_enc to the Newtonian gravitational binding energy of a Plummer sphere, 3πGM²/(32r_h). The effective radius of GCs of DF2 is around 8 pc (Dutta Chowdhury et al. 2020), while the MOND transitional radius, where deviations from Newtonian gravity become important, is 34 pc. If a GC became that big, it would not be classified as a GC, and this is why we consider the Newtonian binding energy. The radii and masses of the GCs of DF2 do not show an obvious correlation and therefore we adopted a universal radius for all GCs.

For the simulation with the DF2-like GC masses, it is suitable to accept the characteristic radius for the GCs as the typical radius of the GCs of DF2, 8 pc (Dutta Chowdhury et al. 2020). For this value we got a total relative change of internal energy greater than (0.1,0.5,1) for (66,61,57) GCs out of the total 100 in the whole ensemble of our 10 simulations with the DF2-like GC masses. Adopting the characteristic radius of the GCs instead as the maximum scale radius of the GCs of the DF2 galaxy (12 pc) (Dutta Chowdhury et al. 2020), we got a total relative change of internal energy greater than (0.1,0.5,1) for (67,62,58) cases. This suggests that a substantial fraction of GCs can be affected by evaporation of GCs due to the encounters. Most of the interactions leading to big gains of internal energy happened during the interactions involving the three most massive GCs. For the simulation with the standard GC masses, we first adopted the characteristic radius of a GC in a typical galaxy (3 pc) (Masters et al. 2010). For this value we got a total relative change of internal energy greater than (0.1,0.5,1) for (14,4,4) GCs. Adopting instead a twice larger characteristic radius for the GCs (6 pc), we got a total relative change of internal energy greater than (0.1,0.5,1) in (17,12,6) cases in the whole ensemble of 10 simulations. In other words, for these lower GC masses, the structure of GCs by GC-GC encounters can affect a small fraction of GCs.

8.2. Mergers between GCs

In order for two GCs to merge, they first have to become bound. After that, they have to lose relative orbital kinetic and potential energy. We therefore estimated the frequency of mergers in two steps that we describe in detail below: First, we detected pairs of GCs in the simulations in relative pericenters and decided the GCs would be bound. Second, we checked whether they had dispositions to transfer the orbital energy into the internal energy.

Let us define the function U(M₁, M₂, b, r) as the binding potential energy of two Plummer spheres with masses M₁, M₂, the characteristic radius b being separated by the distance r. The energy U consists of three contributions: the two-body energy U_2b necessary to separate the two Plummer spheres to infinity and the binding energy of the two Plummer spheres themselves. Let us further assume that the centers of the Plummer spheres follow the average positions and velocities of the stars that originally belonged to the two interacting GCs. Let K be the relative kinetic energy of the spheres. Before the merger, neglecting the influence of all the other bodies in the system, the two clusters oscillate in the U − K plane, keeping the sum E_ini = U + K constant. After the completion of the merger, K = 0 and U is minimized at U(r = 0). In order for the two isolated GCs to merge, the energy E_merg = E_ini − U(r = 0) has to be transferred to the internal energy of the GCs.

For every pair of GCs in a simulation we identified the instances of pericenters and calculated the energy loss by tidal shock as above. A GC pair was considered bound if the difference of the energy of the GC-GC system, U_2b + K, and of the internal energy of the GCs gained by the shock was less than zero. In order for a bound GC pair to merge, a substantial amount of orbital energy has to be transferred to the internal energy of the GCs. If the orbital energy is not lost, encounters with other GCs can make the pair unbound or the pair does not merge during the whole life of the galaxy. We thus determined the candidates for merging clusters such that the loss of energy by the current encounter by tidal shock, ΔE_1, enc, was at least 0.1 E_merg, for U and K evaluated at the pericenter.

In MOND, the binding energy depends on the value of the external field. For isolated objects, it is infinite because isolated objects produce infinitely deep potential wells in MOND. For our purposes, we assumed that the GC pairs are embedded in an external field of the intensity 0.1a₀, which is the typical gravitational acceleration imposed by the galaxy itself on the GCs close to the center of the galaxy (Fig. 1), where we expect most of the interactions between the GCs. The binding energy of two point masses was obtained by integrating the MOND approximate two-body force formula in the presence of an external field by Zhao et al. (2013). Following Bílek et al. (2018), we included a softening parameter equal to the accepted characteristic radius of the GCs to account for the non-zero size of the GCs. The potential binding energy of a Plummer sphere, necessary for evaluating the energy U, was calculated as

$$\begin{array}{c}\hfill U_{\mathrm{Plummer}}={\displaystyle {\int }_0^{\infty }}4\pi G{r}_1^2\rho (r_1){\displaystyle {\int }_{r_1}^{\infty }}\mathcal{T}\left[\frac{-GM(<r_1)}{r_2^2}\right]\mathrm{d}{r}_2\mathrm{d}{r}_1,\end{array}$$(14)

where ρ and M(< r) stand for the standard density and cumulative mass of a Plummer sphere, respectively, and 𝒯 for a function transforming the Newtonian acceleration to the MOND acceleration through the approximate “1D” formula for the EFE in QUMOND (Famaey & McGaugh 2012). The post-merger energy U(r = 0) was evaluated from Eq. (14), when considering the two Plummer spheres overlaid on each other.

In this way, we obtained for the simulation with the DF2-like GC masses 15 candidate merging GC pairs in the whole set of 10 simulations for the characteristic radius of the GCs of 8 pc. Adopting the maximum GC characteristic radius in DF2, 12 pc, we found 8 GC merger candidates pairs. This suggests that for isolated UDGs with GCs similar to DF2, GC mergers can occasionally happen. For the standard GC masses, we obtained no mergers regardless of whether we assumed the GC radii of 3 or 6 pc.

8.3. Loss of orbital energy of GCs from GC interactions

We also explored the possibility that the GCs might sink toward the center of the UDG faster because their orbital energy with respect to the galaxy would be transferred to the internal energy of the GCs because of their interactions with each other. We therefore again detected the times of pericenters of all GC pairs in the simulations and calculated the energy absorbed by the GCs via the tidal approximation approach described above. We then compared the sum of all of these energy losses to the total kinetic energy of all GCs at the end of the simulation.

For the DF2-like GC masses and a characteristic radius of 8 pc, we found that the energy loss is 2.8 times higher than the kinetic energy at the end of the simulation. Adopting the GC radius of 12 pc, 2.3 times the final kinetic energy of GCs is transferred to the internal energy of GCs. This indicates that the GC-GC interactions might also have an effect on the spatial structure of the GC system. This is not the case for the standard GC masses. If their characteristic radius is assumed to be 3 pc, the absorbed energy makes only 0.05 of the kinetic energy of the GC system at the end of the simulation. Adopting the double GC radius, the GC system transferred a fraction of 0.07 of its kinetic energy to the internal energy of the GCs.

9. Summary and conclusions

Dynamical friction in MOND gravity is far less explored than in Newtonian gravity. In this contribution, we aimed to improve this deficiency. We first sought to test for the first time the validity of the proposed MOND analog of the Chandrasekhar formula using high-resolution, self-consistent simulations. We started by studying dynamical friction in an example of a single GC, represented by a point mass, moving in the gravitational field of a UDG with a mass of 2 × 10⁸ M_⊙ similar to the DF2 galaxy (but in isolation). We found that the Sánchez-Salcedo formula (Eq. (4)) works excellently for a suitable choice of the effective value of the Coulomb logarithm, as long as the GC does not have its apocenter closer than about 0.75 kpc to the center of the UDG. For radial orbits and realistic masses of the GC (i.e., less than about 2 × 10⁶ M_⊙), the best value of the Coulomb logarithm is around lnΛ ≃ 3. For higher masses of the GC, between 2 × 10⁶ and 2 × 10⁷ M_⊙ (one-tenth of the mass of the galaxy), the value of the Coulomb logarithm decreases. This is in line with the findings of previous works that dynamical friction in MOND becomes ineffective for interactions of objects with comparable masses. On the other hand, the effective value of the Coulomb logarithm is higher for non-radial orbits. This agrees with the results of Newtonian simulations. It can likely be explained by a higher contribution of the resonant coupling mechanism to dynamical friction. The effective value of the Coulomb logarithm depends not only on the circularity of the orbit, but also on its size: the highest values were encountered for orbits close to the center of the galaxy. The highest effective value of the Coulomb logarithm that we encountered in our simulations was lnΛ = 15.

The Sánchez-Salcedo formula stopped being applicable to our problem when the apocentric distance of the GC decreased below about 0.75 kpc (0.4 R_e) with respect to the center of the galaxy. Dynamical friction becomes much less effective in this central region. The GCs in our simulations, even the most massive ones, never sink right to the center of the galaxy, even if predicted to be the case by the Sánchez-Salcedo formula. We verified that this is not because of insufficient resolution of the simulations. A similar effect, known from Newtonian simulations, is known as core stalling. Whether core stalling occurs in Newtonian simulations depends on the inner density profile of the galaxy, but we did not explore this aspect in MOND. For a GC of 10⁶ M_⊙ falling from 5 kpc (2.5 R_e), it takes about 3 Gyr to reach the stalling phase, but the same takes more than 10 Gyr for a GC with a mass of 10⁵ M_⊙. These inspiralling times depends on the particular density distribution of the galaxy. In total, we have shown that the Sánchez-Salcedo formula can serve as a rule of thumb to estimate the order of dynamical friction. It has a limited range of applicability, but in certain circumstances it is a very good approximation.

We also explored the evolution of a whole GC system consisting of ten GCs around the UDG. We considered two types of GC systems: the GCs were either “overmassive” as deduced in DF2 (if it is at a distance of 20 Mpc), or the GCs had masses following the standard GC mass function. It is not established yet how often anomalous GC mass functions occur in UDGs. Follow-up observations of another UDG with a suspected DF2-like mass function (Müller et al. 2020) have shown that it rather follows a standard GCLF mass function (Müller et al. 2021b).

The more massive GC system undertakes a dramatic structural change during the first 6 Gyr of its evolution, when its half-number radius changes from 3.5 kpc to about 0.6 kpc. The most massive GCs sink more quickly. After that time, the half-number radius remains nearly constant up to the end of the simulation at 10 Gyr, at which point most of the GCs have reached their core-stalling phase. Even in this phase, the GCs would likely be too far from each other to merge together to form a nuclear star cluster, since the average radius of the massive GCs of UDGs is only 6−8 pc (van Dokkum et al. 2018a). Hence we need to be careful when arguing against MOND on the basis of the survival of the GC system with analytic formulas such as Eq. (4). We note that we simulated a galaxy that is in isolation. While its parameters were inspired by the DF2 galaxy, the real DF2 galaxy, if indeed at the measured distance, has to be close to its neighboring galaxy NGC 1052 dominated by its external field and therefore its internal dynamics is effectively Newtonian. The GC system of DF2 in Newtonian dynamics has been studied by Dutta Chowdhury et al. (2019), who found that it is consistent with observations if it was somewhat more extended in the past. The GC system of DF2 is thus consistent with Newtonian and MOND gravity. The fact that the GC system of DF2 does not have a core (see Sect. 6.1) then suggests that DF2 has spend most of its life dominated by the external field of NGC 1052.

The simulated GC system following the standard GC mass function experiences a weaker evolution, decreasing its 3D half-number radius from 3.5 kpc to about 2 kpc in 10 Gyr. Extrapolating this result backward in time, we expect that GC systems of isolated UDGs were larger by a factor of almost 2 at their formation 10 Gyr ago with the standard GCLF.

For both choices of the mass function of the GC system, we detected a mass segregation of the GCs. This can be taken as a prediction for observations. We found that the GC system of the real DF2 galaxy has no signature of the central core. This indicates that in MOND, if DF2 is a member of the NGC 1052 group, it has to have spent a large fraction of its life dominated by the EFE imposed by this galaxy. Further, we were able to show by scaling our simulations for most GC systems of UDGs that they are not expected to merge and form nuclear star clusters.

Not only the spatial distribution of GCs evolved in our simulations, but also the central velocity dispersion of the GC system decreased and that of the stars of the galaxy increased, likely because of energy equipartition (see Figs. 10 and 11). We can thus expect that in observations, UDGs usually have a higher velocity dispersion of stars than of GC systems.

Our simulations were limited in the sense that the GCs were modeled as point masses instead of bodies consisting of many particles. The many GC simulations were thus not able to capture the effects related to absorbing the orbital energy of the GC system in the internal energy of the individual GCs. We partly resolved this issue, at least for a UDG of the properties of the simulated galaxy, by joining the simulations with analytic estimates. For the simulations with a GC system inspired by DF2, we found that GC-GC interactions are expected to play a substantial role in the evolution of the GC system. About one-half of the GCs are expected to receive energy that is comparable to their internal Newtonian binding energy. Therefore, not only the GCs might lose a large fraction of their stars because of the encounters, but the GC system might sink substantially faster to the center of the galaxy because of the loss of orbital energy. For the DF2-like GCs, we therefore still cannot reliably exclude the possibility that some of the GCs would form a nuclear star cluster in an isolated UDG. This is to be investigated by dedicated simulations with GCs simulated as many-body systems. For the standard GCs, which seem to be the most common even in UDGs, our analytic estimates indicate a much less dramatic influence of mutual GC interactions. The increase in internal energy of GCs is expected to be comparable to the Newtonian binding energy of GCs for only about one GC out of the ten in the GC system. We estimated that the GC-GC mergers are nearly non-existent, and there is no acceleration of the sinking of the GC system to the center of the UDG. We thus do not expect that the GCs would form a nuclear star cluster if the GCs were simulated as many-particle systems.

The question of survivability of GC systems in MOND remains open for dwarf galaxies. Such galaxies have lower effective radii and the same or lower brightnesses than UDGs. For example, for a dwarf galaxy similar to the Fornax dwarf, having a stellar mass of 3 × 10⁷ M_⊙ and an effective radius of 0.5 kpc, we deduce a GC-mass scaling factor of 6.7 (Eq. (12)) and a time scaling factor of 2.6 (Eq. (13)). The evolution of the GC system due to dynamical friction is thus much faster than in our simulation of the UDG with the DF2-like GC mass function. It is difficult to predict from our present simulations whether a nuclear star cluster would form. Observations of chemical properties of some dwarf galaxies suggest that their nuclear star clusters formed by dynamical sinking of GCs (Fahrion et al. 2020); the galaxies had much smaller radii and brightnesses than the UDGs investigated in this work. Future work should address in which galaxies MOND predicts the formation of a nuclear star cluster through the merging of the GC system.

Previous simulations showed that major mergers are much less effective in MOND than in Newtonian simulations with dark matter. Analytic arguments however indicated that the situation is the opposite for minor mergers around spheroidal galaxies. Our simulations confirmed, for the first time, the analytic prediction for situations similar to small galaxies being accreted by much more massive spheroidal galaxies (previous works confirmed them only for the rotation of galactic bars). This forces us to consider seriously the potential importance of minor mergers for galaxy evolution in MOND. Nevertheless, minor merging will require a more careful examination. The first reason is that Newtonian simulations show that the size of the satellite can also influence dynamical friction (Mo et al. 2010), while our simulated satellite was a point mass. It should also be explored better how dynamical friction works in MOND when the satellite moves outside of the main galaxy, where the density of stars is negligible. In Newtonian gravity, in these regions we still have dark matter halos up to hundreds of kiloparsec from the galaxy that can absorb the energy of the satellite. These are questions are to be investigated by future works.

Acknowledgments

The authors thank the anonymous referee for a constructive report which helped to significantly improve the paper. MB and HSZ are thankful for the financial support by Cercle Gutenberg. MB acknowledges the support from the Polish National Science Centre under the grant 2017/26/D/ST9/00449. BF and RI acknowledge financial support from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation programme (grant agreement No. 834148). O.M. is grateful to the Swiss National Science Foundation for financial support.

References

All Tables

All Figures

Fig. 1. Properties of the simulated galaxy just before entering GCs. Top: radial profile of density. Middle: radial profile of velocity dispersion, separately for the full velocity dispersion (full line), velocity dispersion in the radial direction (dotted line), and in the tangential direction (dashed line). Bottom: radial profile of gravitational acceleration in the units of the MOND constant a0.

In the text

	Fig. 2. Orbital decay of the GC in the fiducial model. The decay slows down substantially once the apocentric distance of the GC reaches about 0.5−1 kpc.
In the text

	Fig. 3. Demonstration of the effect of changing of the spatial resolution in the fiducial model on the orbital decay of the GC. The core stalling phase is present for all three choices of resolution.
In the text

Fig. 4. Testing the precision of the Sánchez-Salcedo formula for dynamical friction in MOND (Eq. (4)). Blue lines: Evolution of the GC-galaxy distances in N-body simulations. Orange lines: The same but calculated analytically employing the Sánchez-Salcedo formula. Each panel corresponds to another mass of the GC under consideration. The GCs in all simulations or calculations started with the same initial conditions, that is, starting from rest with respect to the galaxy and at a distance of 5 kpc from the galaxy center. To the right of each figure, the corresponding mass of the GC and the value of the Coulomb logarithm that is needed to be entered in the Sánchez-Salcedo formula to mimic the result of the simulation are stated.

In the text

Fig. 5. As in Fig. 4, but varying orbital parameters of the GC. The GC always had the mass of 106 M⊙. The first lines of the notes on the right of the panels indicate the initial tangential velocity of the GC with respect to the galaxy (the radial velocity was zero). It is given in the units of the local circular velocity. The second row indicates the initial galactocentric distance and the third the value of the Coulomb logarithm in the Sánchez-Salcedo formula that provides the best match of the simulation and analytic calculation.

In the text

	Fig. 6. Example of evolution of galactocentric radii of GCs in a simulation with 10 GCs. The masses of the GCs are set as observed for the GCs of the DF2 galaxy (first column of Table 2).
In the text

Fig. 7. Evolution of the distances of GCs from the center of the UDG. The black points show combined data from 10 random realizations of a UDG orbited by 10 GCs. The masses of the GCs in each simulation mimic the measured masses of the GCs of the DF2 galaxy (see the first column of Table 2). The horizontal black dashed line indicates the maximum resolution of the simulation. The red lines indicate the median galactocentric distance of the GCs at the given time. This was derived separately for the high-mass half of the GCs (full red line) and for the low-mass half of the GCs (dashed red line).

In the text

Fig. 8. Evolution of GC systems. Left column: simulations with 10 GCs that have the masses of those observed in DF2 (left column of Table 2). Simulations of 10 random realizations of the system were combined in these plots. Galactocentric distance is given in units of the effective radius of the galaxy (Re = 1.9 kpc). Right column: same but for GC masses set according to the standard GCLF (right column of Table 2). Top row: profiles of number density of GCs as a function of galactocentric distance. The individual lines correspond to different times since the start of the simulation. Middle row: evolution of the profile of the surface number density of GCs. Bottom row: evolution of the radial profile of the cumulative fraction of GCs.

In the text

	Fig. 9. Time evolution of the concentration parameter R of the GC system defined by Eq. (10). The two lines correspond to the two considered distributions of the masses of GCs (see Table 2).
In the text

	Fig. 10. Line-of-sight velocity dispersions of the GC systems in our simulations as functions of time.
In the text

	Fig. 11. Radial profile of the velocity dispersion profile of the stars and of the GC system of the simulated galaxy. The curves correspond to the beginning and end of the simulations, as indicated in the legend. The vertical lines indicate the 1σ uncertainty limits.
In the text

	Fig. 12. Same as in Fig. 7 but with GC masses set according to the standard GCLF.
In the text

Fig. 13. Time scaling factor given by Eq. (13) as a function of the absolute magnitude and effective radius of a galaxy. The white point indicates our simulated galaxy. The white dashed line indicates where the scaling factor equals 1. The black dashed line denotes where the GC-mass scaling factor, the second factor in Eq. (12), equals 2.66, which is the value where an observed UDG with GCs following the standard GCLF mass function has to be compared to our simulations with the DF2-like GC mass function. Supposing that the galaxy has a GC system following the standard GC mass function, our results imply that the galaxies to the right of the black line and to the top the white line do not create a nuclear star cluster by inspiralling of GCs.

In the text