© The Authors 2024

1. Introduction

Globular clusters (GCs) are compact (a few parsecs), massive (10⁴ − 10⁶ M_⊙) star systems found in nearly all galaxies. A galaxy similar to the Milky Way has a few hundreds of them, while giant ellipticals can have more than ten thousand GCs. The colors of many galaxies form a bimodal distribution, with rather universal positions for the two peaks. Therefore, GCs are divided into two types – the metal-poor “blue GCs” and the metal-rich “red GCs” (Brodie & Strader 2006; Cantiello et al. 2020). Red GCs generally follow the kinematics of the stars in a galaxy, with a similar rotational velocity and velocity dispersion. In contrast, blue GCs often show complex kinematics (Schuberth et al. 2010; Coccato et al. 2013; Chaturvedi et al. 2022). The spatial distribution of GCs around galaxies is more centrally concentrated for the red GCs than for the blue GCs. The distinct properties of the red and blue GCs point toward their different formation pathways (Ashman et al. 1995; Peng et al. 2006; Brodie & Strader 2006). It seems that the blue GCs are added to massive galaxies via accretion of low-mass galaxies, while most red GCs form in situ, together with the stars of the host galaxy (Côté et al. 1998; Harris 2001; Tonini 2013; Renaud et al. 2017). Globular cluster systems thus carry information about the assembly history of their host galaxies (Peng et al. 2008; Brodie et al. 2014; Harris et al. 2016). The number of GCs that a galaxy hosts is proportional to the expected mass of its dark matter halo (Spitler & Forbes 2009; Harris et al. 2015). The kinematics and distribution of GC systems (GCSs) reflect the profiles of the gravitational fields of their host galaxies (Samurović 2014, 2016; Alabi et al. 2017; Bílek et al. 2019b).

In the paper by Bílek et al. (2019a; BSR19 hereafter), an interesting new property of GCSs of early-type galaxies was noted. They parametrized¹ the volume number density of GCs in a GCS, ρ by a broken power law as

$$\begin{array}{cc}& \rho (r)={\rho }_0{r}^a\text{for}r<r_{\mathrm{br}},\hfill \\ \hfill & \rho (r)={\rho }_0{r}_{\mathrm{br}}^{a-b}{r}^b\text{for}r\ge r_{\mathrm{br}},\hfill \end{array}$$(1)

where r is the galactocentric radius. The parameter r_br was called the break radius. The authors found that the break radius coincides well with the a₀ radius, that is the radius at which the expected gravitational acceleration generated by the baryons of the galaxy equals the galactic acceleration scale a₀ = 1.2 × 10⁻¹⁰ m s⁻². The values of the break radii did not agree with the values of the other characteristic lengths of the galaxies, such as stellar effective radii or dark halo scale radii. These other lengths were either several times bigger or smaller than the break radii, at least for a substantial fraction of the galaxy sample (see their Table 1).

The galactic acceleration scale is known best from the behavior of the observed gravitational fields of galaxies (e.g., Lelli et al. 2017; Li et al. 2017). In the regions of galaxies where the gravitational acceleration expected by Newtonian gravity from the distribution of baryons g_N is greater than a₀, the observed gravitation acceleration equals g_N, meaning that the Newtonian dynamics does not require dark matter. On the other hand, in the regions where g_N is lower than a₀, the observed gravitational acceleration is very close to $\sqrt{g_{\mathrm{N}}{a}_0}$. The same rules apply even for many, or perhaps all, GCs (Scarpa & Falomo 2010; Scarpa et al. 2011; Ibata et al. 2011; Sanders 2012; Hernandez et al. 2012; Hernandez & Lara-D I 2020).

This behavior was initially predicted by the modified Newtonian dynamics (MOND), which is a class of modified gravity and inertia theories (Milgrom 1983c). Here we assume MOND to be a modified gravity theory. It predicts that the gravitational acceleration in spherical isolated objects is (Milgrom 1983c, 2010; Famaey & McGaugh 2012)

$$\begin{array}{c}\hfill g_{\mathrm{M}}=g_{\mathrm{N}}\nu (g_{\mathrm{N}}/a_0).\end{array}$$(2)

The function ν is not known exactly, but it must have the limit behavior ν(x)∼x^−1/2 for x ≪ 1, and ν(x)∼1 for x ≫ 1. This gives rise to two regimes of a gravitational field around a galaxy: the strong field, the so-called Newtonian regime, and the weak field, the so-called deep-MOND regime. The observed counterpart of Eq. (2) is known as the radial acceleration relation (McGaugh et al. 2016).

Apart from the radial acceleration relation, MOND predicted or explained many other observational laws (Milgrom 1983a,b,c), all of which contain the constant a₀. This is the case for the baryonic Tully-Fisher relation (McGaugh et al. 2000; Lelli et al. 2019), the Faber-Jackson relation (Faber & Jackson 1976; Famaey & McGaugh 2012), and the radial acceleration relation, which connect the mass or mass distribution of galaxies to the velocities of stars and gas in them. The Fish law (Fish 1964; Allen & Shu 1979) and Freeman limit (Freeman 1970; McGaugh et al. 1995; Fathi 2010; Famaey & McGaugh 2012) give upper limits on the surface brightness for elliptical and spiral galaxies, respectively, above which galaxies are rare. Recently, there appeared a MOND explanation (Milgrom 2021) of the Fall relation (Fall 1983; Posti et al. 2018), which connects the mass and specific angular momentum of galaxies. The law of the universal surface density of the cores of the putative dark matter halos (Kormendy & Freeman 2004; Donato et al. 2009; Salucci et al. 2012) can be explained by MOND as well (Milgrom 2009).

Finally, there are interesting numerical coincidences of a₀ with the constants of cosmology (Milgrom 1983c, 2020). If we denote H₀ the Hubble constant, c the speed of light, G the gravitational constant, R_H the size of the cosmic horizon, and M_H the total mass inside the cosmic horizon, then we find the order-of-magnitude equalities a₀ ≈ cH₀ ≈ c²Λ^1/2 ≈ c²/R_H ≈ c⁴/GM_H. There is no clear explanation of these coincidences yet (Navarro et al. 2017; Milgrom 2020).

The finding of BSR19, of the equality of the break and a₀ radii, is thus another case of the many occurrences of the constant a₀ in extragalactic astronomy. More precisely, in this work we consider two types of a₀ radii: i) the one where the acceleration calculated from the distribution of baryons and Newtonian gravity equals a₀ and ii) where the acceleration calculated for MOND gravity via Eq. (2) equals a₀. For most galaxies, the two a₀ radii are numerically similar. Therefore, in this paper, if we do not specify whether we are referring to Newtonian or MOND a₀, then we mean that the statement is valid for both options.

The theoretical explanation for why the a₀ radii coincide with the break radii of GCSs has not been clarified yet, even if some initial proposal were given in BSR19. Importantly, according to one of the proposed explanations that involves the Newtonian gravity and dark matter, the match of the a₀ and break radii is of practical importance. It had been found before that the number of GCs that a galaxy has is proportional to the mass of its dark matter halo (Spitler & Forbes 2009; Harris et al. 2015). The new finding allows one to estimate the scale radius of the halo: the break radius should be located at the radius where the gravitational attraction of the stars of the galaxy equals that of the dark matter halo. One can thus solve the equation of the equity of the accelerations to obtain the scale radius of the halo.

The paper BRS19 left several important questions open, which we aim to answer here. We investigate the distribution of GCs primarily using photometric data for early-type galaxies in the Fornax galaxy cluster, but we also analyze new spectroscopic data for two galaxies. The galaxy sample of BSR19 spanned only about one order of magnitude in stellar mass. The data investigated here allow us to verify that the match between the a₀ and break radii holds true for early-type galaxies spanning three orders of magnitude – from dwarfs of the mass of the Magellanic clouds to the brightest cluster galaxies. The distribution of GCs was investigated in BSR19 on the basis of the catalog of spectroscopically confirmed GCs. This approach can lead to distorted results, because spectroscopic surveys usually are spatially incomplete. This could have affected the derived parameters of the broken power-law profile. The verification of the match of the a₀ and break radii in the new data thus removes the shade of doubt that was left about the results of the previous work. The data in BSR19 do not allow the density profiles of the red and blue GC subpopulations to be investigated. We do this here and find that there are no statistically significant trends of the profile parameters with the color of the GCs. We exploit the new data for a further exploration of the profiles of GCSs. In particular, we investigate whether the parameters of the broken power-law profiles correlate with each other and with the parameters of the host galaxy. We also investigate in detail the GCS of NGC 1399, that is the central galaxy of the Fornax cluster. We find that its break radius depends on the relative velocity of the GCs with respect to the center of the galaxy, and that the break radius varies as a function of the position angle. In this paper, we also aim to explain the reason of why the GCSs of our galaxies have the broken power-law density profiles and why the break radii coincide with the a₀ radii. Several explanations were proposed in BSR19, and here we add a few more. Then we make the first steps toward finding which of them is correct. None seem perfect at this point.

This paper is organized as follows. In Sect. 2 we describe the observational data that we analyze here. The methods to extract and fit the radial profiles of volume number densities of GCs in the GCSs of the investigated galaxies are detailed in Sect. 3. The derivation of the a₀ radii for our sample galaxies is explained in Sect. 4. We present our results in Sect. 5. In particular, in Sect. 5.1 we compare the structural parameters of GCSs for the total GC population and the red and blue subpopulations. In Sect. 5.2 we investigate the correlations of the structural parameters between each other and with the characteristics of the host galaxies and, finally, we explore the relation between the break and a₀ radii in Sect. 5.3. We explore the details of the distribution of GCs of NGC 1399 in Sect. 6. We explore the credibility of some potential explanations of the approximate coincidence of the a₀ and break radii in Sect. 7. Finally we synthesize and summarize our findings in Sect. 8. In this work, we denote the natural logarithm by log and the logarithm of the base m by log_m. We assume the distance of the Fornax cluster and all the investigated galaxies to be 20.0 Mpc (Blakeslee et al. 2009). This corresponds to the angular scale of 5.8 kpc per arcmin.

2. Data analyzed

The galaxies investigated in this paper are member galaxies of the Fornax galaxy cluster. They are listed in Table 1. Low-mass Fornax cluster galaxies (M_* < 10^9.5 M_⊙) have too few GCs for constructing the density profiles of their GCSs individually. Therefore, we stacked the GC candidatess of many low-mass galaxies in three mass bins in order to get their average GCS density profiles. These are the three “Stack” entries in Table 1. The details of the stacking procedure are explained in Sect. 2.3.

All structural parameters of the galaxies in Table 1 were taken from Su et al. (2021). They are based on GALFIT profile fitting (Peng et al. 2002), using Sérsic functions, to photometric g′-band data of the Fornax Deep Survey (FDS; Iodice et al. 2016; Venhola et al. 2018). Stellar masses also come from that work. They were derived from empirical relations between colors and stellar mass-to-light ratio. We use archival photometric and spectroscopic catalogs of GC candidates for investigating their spatial distribution. Here follows a brief description of the datasets investigated in our work.

2.1. ACS Fornax cluster survey data

The ACS Fornax cluster survey (ACSFCS), taken using the Advanced Camera for Surveys (ACS) of the Hubble Space Telescope (HST), imaged 43 early-type galaxies of the Fornax cluster. Full details of the ACSFCS, scientific motivations and data reduction techniques, are given in Jordán et al. (2007). Each galaxy in ACSFCS was imaged in the F475W (g) and F850LP (z) bands. For studying the GCs, each image was sufficiently deep such that 90% of the GCs within the ACS FOV can be detected. The selection and identification of bonafide GCs of the ACSFCS galaxies are performed in the size-magnitude plane (Jordán et al. 2015). The resulting catalog of GCs in the ACSFCS provides the probability of an object being a GC, denoted by Pgcs, where Pgcs is a function of half-light radius, apparent magnitude and local background. For our analysis, we selected GCs with Pgcs larger than 0.5. This leaves the faintest GC candidates of m_g = 26.3 mag. To separate the GCs into red and blue subpopulations, we adopted a dividing g − z color of 1.1 mag (Fahrion et al. 2020).

2.2. Spectroscopically confirmed GCs

We studied the spectroscopically confirmed sample of GCs of the Fornax cluster from the recent catalog produced by Chaturvedi et al. (2022). Reanalyzing data of the Fornax cluster taken using the Visible Multi-Object Spectrograph (VIMOS) at the Very Large Telescope (VLT; Pota et al. 2018) and adding literature work, they have produced the most extensive GC radial velocity catalog of the Fornax cluster (see Chaturvedi et al. 2022, for details), comprising more than 2300 confirmed GCs. The faintest GC has m_g = 24.2 mag. They used a Gaussian modeling mixture technique to divide the GC population into red and blue GCs, with g − i ∼ 1.0 mag as separating color, which we adopt in our analysis.

2.3. Fornax Deep Survey data

The Fornax Deep Survey (FDS) is a joint project based on a guaranteed time observation of the FOCUS (P.I. R. Peletier) and VEGAS (P.I. E. Iodice, Capaccioli et al. 2015) surveys. It consists of deep multiband (u, g, r and i ) imaging data from OmegaCAM at the VST (Kuijken 2011; Schipani et al. 2012) and covers an area of 30 square degrees out to the virial radius of the Fornax cluster.

We applied morphological and photometric selection criteria to the photometric compact sources catalog of the FDS (Cantiello et al. 2020) to decrease the fraction of contaminant objects that are not GCs. The criteria on the colors of the selected objects g − i > 0.5 and g − i < 1.6 were chosen according to the colors of spectroscopically confirmed GCs around the central cluster galaxy NGC 1399. Further criteria were inspired by cross-matching the FDS and ACS catalogs of GCs of NGC 1399. They were chosen such that we do not exclude too many real GCs and, at the same time, exclude as many contaminants as possible. We found a good balance when using the following criteria: CLASS_STAR > 0.031, m_g > 20 mag, and Elongation < 3. The meaning of the parameters is explained in Cantiello et al. (2020). After applying the selection criteria, the faintest GC candidate in this calatog has m_g = 27.0 mag.

It turned out that the FDS catalog, after applying the GC candidate selection criteria, shows systematic variations of the surface density of sources that have a tile-like pattern, see Fig. 1. The tiles corresponds to the OmegaCAM imaging tiles of the FDS survey. They probably result from varying observing conditions, like seeing, sky transparency, etc. When fitting the surface density profiles of GCSs of individual galaxies, we had to be careful that the tile borders do not introduce any kinks in the profiles. We did this first by visual inspection of plots of positions of the GC candidates in the wide neighborhood of the investigated galaxies and second by visual inspection of the plots of the radial profiles of density of the GC candidates around the target galaxies.

Fig. 1. Positions of GC candidates in the FDS catalog, after applying selection criteria. The coordinate system is centered on NGC 1399. North is up, and west is to the right.

As mentioned above, we decided to stack the faint galaxies of similar stellar mass in order to have enough sources to extract the density profiles of their GCSs. We used only the FDS data for this, since most faint galaxies were not covered by the ACSFCS. If the hypothesis of this paper, that the break radii coincide with the acceleration radii, is correct, then the break radii of all GCSs stacked in this way should be roughly equal in each mass bin (an acceleration radius nevertheless depends also on the particular distribution of mass in the given galaxy). We then treated the stacks as single galaxies. We created three such artificial objects. Their logarithmic masses were centered on the values of 8.0, 8.5 and 9.0 M_⊙. The widths of all logarithmic mass bins were 0.5 M_⊙. Stacking of galaxies of even lower masses did not provide sufficiently clear profiles of the projected density of GCs.

We stacked all galaxies from the Su et al. (2021) catalog that met the following criteria. First we excluded spiral galaxies, because our data do not allow to distinguish between GCs and star forming clumps or young star clusters in their disks. Early-type galaxies were identified by requiring their asymmetry parameter, stated in Su et al. (2021), to be lower than 0.06. We further excluded galaxies that are located close to the borders of the tiles of the FDS mosaic, and galaxies whose GCSs surface density profiles might be affected by interloping GCSs of other nearby galaxies. These were identified by visual inspection of the positions of the GC candidates in the neighborhood of the galaxies being stacked. The list of stacked objects in every mass bin is stated in Table A.1.

We assigned to every stacked galaxy a stellar mass and Sérsic parameters. These were obtained as the median values for all galaxies included in the corresponding stack.

3. Fitting the surface density profiles of the GCSs by analytic functions

Here we describe how we determined for the investigated galaxies the parameters of the GCSs density profiles.

3.1. Extracting the observed profiles of surface density

For a given galaxy, we divided the GC candidates in radial bins according to their galactocentric distance, and for each bin, we calculated the surface density of sources in it. The bins had the shape of circular annuli, that is we ignored the possible ellipticity of the GCSs. We demonstrate in Appendix C that this simplification has no appreciable impact on the derived profile of the GCS density. We chose the widths of the radial bins such that

$$\begin{array}{c}\hfill N=n-A\overline{\gamma }\end{array}$$(3)

was constant in each bin. Here n stands for the total number of sources falling in the annulus, A the area of the annulus and $\overline{\gamma }$ an estimate of the surface density of contaminants. The number N is thus the estimate of the number of GCs in the aperture. The parameter N was chosen to be big enough so that the profile did not appear too noisy (as judged by eye) and small enough so that the two straight parts of the broken power law profile were resolved by at least three data points. The purpose of subtracting the expected number of contaminants was to increase the signal to noise in the outer bins. This condition was used only for choosing the bin widths; the final profile parameters of the GCSs (including the contaminants level) were derived from fitting the surface density of all sources, that is n.

For the datasets that are not expected to contain many contaminants, that is the spectroscopic and ACS catalogs, we used simply $\overline{\gamma }=0$. Otherwise the value of $\overline{\gamma }$ was found iteratively in the following way. We first chose a low value of $\overline{\gamma }$ and constructed the observed surface density profile. That was then fitted by some of the analytic model profiles described below. One of the parameters of the fitted profile was the density of the background sources, γ. In the next iteration we used a value of $\overline{\gamma }$ that lied between the value of $\overline{\gamma }$ used in the previous iteration and the fitted value of γ. We had to stop increasing $\overline{\gamma }$ at some point, because when it is too large, the extracted surface density profile would be truncated or distorted, because the expression on the right-hand-side of the condition Eq. (3) would be lower than the assigned value of N at large radii. In some cases, the final $\overline{\gamma }$ was chosen such that the extracted profile did not show any small bumps arising from small clusters of contaminants (e.g., distant galaxy clusters in the background). We note that if the definitive $\overline{\gamma }$ was chosen to be somewhat below the true value, the resulting extracted profile would not be affected substantially; only the bins would not have the optimal widths.

The contamination by the light of the host galaxy makes it difficult to detect GCs near its center. In order to detect them, we have to subtract a model of the light of the galaxy from the image. The model is typically not perfect, and the fit residuals in the difference image still complicate the detection of the GCs. Furthermore, the light of the galaxy introduces photon noise, that decreases the signal-to-noise ratio of the GCs. Fainter GCs are affected more. This can deform the observed GCS density profile. The ground-based FDS data are more sensitive to these problems than the HST data: the GCs in HST images appear sharper because of no atmospeheric seeing effects, and therefore reach a higher signal-to-noise ratio. We indeed found the signature of this in our data: the inner slopes of the profiles of surface density of GC candidates were shallower in the FDS data than in the HST data. Examples can be seen in Fig. A.14 (NGC 1427) or Fig. A.8 (NGC 1380B).

We were using two strategies to mitigate the flattening of the inner GC density profile because of the problems with the contamination by the light of the host galaxy. If there were HST data available in the central region, we adopted the inner slope a derived from these data. For galaxies without HST data, we constructed the GCS surface density profiles only from bright-enough sources in FDS. We determined the magnitude cut using a plot magnitude – projected galactocentric distance of the sources. The limiting magnitude was set as that of the faintest sources that were still detectable in the very center of the galaxy. This approach has the downside that it reduces substantially the number of GC candidates, such that it becomes harder to trace the GC density profile, particularly at large galactocentric distances. In Appendix B we demonstrate on the example of NGC 1399 that the break radius does not depend significantly on the magnitude cut.

The surface density profiles of the GCSs were extracted and analyzed only within a restricted range of the projected galactocentric radius defined by the limits r_min and r_max. The upper limit, r_max, was usually enforced either by the proximity to the GCSs of other galaxies, or because of the fluctuations of the density of contaminants, or because there was a border of tiles of the FDS mosaic. The lower limit, r_min, was usually taken to be zero, but in some galaxies, it was used to reduce the problems with the contamination by the light of the host galaxy. The radii r_min and r_max were determined by a visual inspection of both the map of the GC candidates near the inspected galaxy and the surface density profile of the GC candidates.

When constructing the GCS surface density profiles for a galaxy, we had to exclude the areas occupied by the GCSs of neighboring galaxies. We excluded sources close to the following major galaxies (unless studying the GCSs of these galaxies themselves): ESO 358-33, NGC 1396, NGC 1317, NGC 1373, NGC 1374, NGC 1375, NGC 1379, NGC 1380, NGC 1381, NGC 1382, NGC 1386, NGC 1387, NGC 1389, NGC 1404, NGC 1427 and NGC 1427A. We marked all sources closer than 3′ to these galaxies as excluded. This distance was chosen by visual inspection of the map of positions of all sources. We note that we did not exclude the area occupied by the GCS of NGC 1399 because it is very extended. If we did so, it would not be possible to construct the observed surface density profiles of several smaller galaxies in the vicinity of NGC 1399. Instead, whenever reasonable, we just assumed that the GCs of NGC 1399 are distributed with a constant surface density in the vicinity of the investigated smaller galaxies, so that we could treat them as additional contaminants. This approach was not applied to NGC 1404, which has a rather extended GCS and is located very close in projection to NGC 1399. Due to the changing number density of NGC 1399 GCs across the area of NGC 1404, we could not consider them as uniformly distributed contaminants and we had to use another method, see Sect. 3.1.

For a few galaxies, we had to mark as excluded sources in the regions of blemishes of the FDS survey. Most often these were regions around bright stars that appear as holes in the maps of the positions of sources from the FDS catalog.

In general, when calculating surface density in a given annulus, we did not consider sectors in which there was at least one excluded source. In total, the observed surface density of sources in a given annulus was calculated as

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathrm{obs}}=\frac{2\pi }{(2\pi -\alpha )}\frac{n}{A},\end{array}$$(4)

where α stands for the sum of angular extents of all excluded sectors, n for the number of sources in the not-excluded sectors, and A for the total area of the inspected annulus. Given that the number of sources in a given area follows the Poisson distribution, we estimated the uncertainty of the surface density of sources as

$$\begin{array}{c}\hfill \mathrm{\Delta }{\mathrm{\Sigma }}_{\mathrm{obs}}=\frac{2\pi }{(2\pi -\alpha )}\frac{\sqrt{n}}{A}.\end{array}$$(5)

The galaxy NGC 1404 lies in projection close to the central cluster galaxy NGC 1399. Their GCSs overlap. The GCS of NGC 1399 has a strong gradient of surface density at the position of NGC 1404. The method described above would not allow us to produce a reliable GC density profile of this galaxy. Instead, we use the catalog of spectroscopic GCs by Chaturvedi et al. (2022) for NGC 1404 and assume that it covers the galaxy homogeneously. We could then make use of not only the information about the spatial positions of the sources, but also the information about their radial velocities, since the radial velocity of NGC 1404 is ∼520 km s⁻¹ larger than that of NGC 1399. We first applied spatial criteria on the sources to be used for constructing the observed profile of the GCS. They are depicted in the upper-right panel of Fig. A.13. We avoided regions that are closer than 8′ to NGC 1399. We also avoided the region that is closer than 10′ to the point with the J2000 coordinates (54.86017, −35.75988), because this region seemed to suffer from geometrical incompleteness of the spectroscopic survey because the density of sources in this region was lower than in its surroundings. Just as for all other galaxies, we applied a limit for the maximum distance of the used sources from the galaxy. We also applied a radial-velocity limit: all sources that had radial velocities lower than the center of NGC 1404, 1947 km s⁻¹, were excluded, as shown in the bottom-right panel of Fig. A.13. This helped us to reduce substantially the contamination by the GCs of NGC 1399. While the observed number density profile was constructed from only 23 sources and there were only 2–3 sources per bin, the break in the profile is clearly visible and the break radius follows the same correlations as the break radii of the other galaxies.

3.2. Models of density and surface density profiles of GCSs

The extracted profiles of surface density of GCs candidates were fitted by one of the analytic functions described in this section. We made use of the fact that for a spherically symmetric GCS, the 3-dimensional density profile ρ corresponds to the projected surface density profile Σ given by the Abel transform:

$$\begin{array}{c}\hfill \mathrm{\Sigma }(R)=2{\int }_R^{\infty }\frac{\rho (r)r}{\sqrt{r^2-R^2}}\mathrm{d}r.\end{array}$$(6)

Here r stands for actual galactocentric distance and R for the projected galactocentric distance.

We considered several types of volume density profiles. The first was a pure broken power law given by Eq. (1). Its corresponding surface density profile is given by the equations

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathrm{BPL}}(R)=\{\begin{array}{cc}2{\rho }_0\{f_0(u_{\mathrm{br}},a)-f_0(0,a)+\hfill & \\ r_{\mathrm{br}}^{a-b}[f_0(\infty ,b)-f_0(u_{\mathrm{br}},b)]\}\hfill & \text{for}R<r_{\mathrm{br}}\hfill \\ \hfill \\ 2{\rho }_0{r}_{\mathrm{br}}^{a-b}[f_0(\infty ,b)-f_0(0,b)]\hfill & \text{for}R\ge r_{\mathrm{br}},\hfill \end{array}\end{array}$$(7)

where

$$\begin{array}{c}\hfill f_0(u,\alpha )=R^{\alpha +1}u{}_2{F}_1(\frac{1}{2},\frac{-\alpha }{2};\frac{3}{2};-u^2),\end{array}$$(8)

$$\begin{array}{c}\hfill u=\sqrt{{\left(\frac{r}{R}\right)}^2-1},\end{array}$$(9)

the symbol ₂F₁ denotes the Gaussian hypergeometric function, and u_br = u(r_br). This form of surface density profile was used to fit only the datasets that contain a negligible fraction of contaminants, that is usually for the GC candidates coming from the ACS catalogs.

For the datasets that contain contaminants, we had to add one more free parameter in the model profile, γ, that expresses the surface density of contaminants:

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathrm{BPLB}}(R)={\mathrm{\Sigma }}_{\mathrm{BPL}}(R)+\gamma .\end{array}$$(10)

The contaminants were thus assumed to be distributed homogeneously, that is with no density gradient. This profile was mostly used for fitting of the profiles extracted from the FDS catalog.

For some galaxies, it was not possible to fit the surface density profile of GCS by a broken power law. This happened either when the GCS was so extended that the ACS field captured only the inner part of the broken power law profile, or the broken power law would be only poorly constrained because the galaxy had too few GCs. In these situations, we only made a fit by a simple power-law density profile. The density profile of a single power law was parametrized as

$$\begin{array}{c}\hfill \rho (r)={\rho }_0{r}^a,\end{array}$$(11)

which corresponds to the surface density profile

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathrm{SPL}}(R)=\sqrt{\pi }{\rho }_0{R}^{a+1}\frac{\mathrm{\Gamma }\left(\frac{-a-1}{2}\right)}{\mathrm{\Gamma }\left(\frac{-a}{2}\right)}.\end{array}$$(12)

When contaminants were expected to contribute to the surface density profile substantially, we added the background term to the single power law:

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathrm{SPLB}}(R)={\mathrm{\Sigma }}_{\mathrm{SPL}}(R)+\gamma .\end{array}$$(13)

3.3. The fitting method

We found the best-fit parameters of the models by the maximum likelihood method. Here we assumed the uncertainty of the number of sources in each bin follows the Gaussian distribution. Then the likelihood reads:

$$\begin{array}{c}\hfill ln\mathcal{L}(\mathit{p})={\displaystyle \sum _i}-0.5ln\left(2\pi \right)-ln\left(\mathrm{\Delta }{\mathrm{\Sigma }}_{\mathrm{obs},i}\right)-\frac{{[{\mathrm{\Sigma }}_{\mathrm{m}}(r_i,\mathit{p})-{\mathrm{\Sigma }}_{\mathrm{obs},i}]}^2}{2\mathrm{\Delta }{\mathrm{\Sigma }}_{\mathrm{obs},i}^2}.\end{array}$$(14)

Here p = (p₁,p₂,p₃,…,p_ν) denotes the vector of the free parameters of the model, and r_i the central radius of the ith bin, that is the arithmetic average of the inner and outer radius of the given annulus, Σ_m denotes the surface density predicted by the fitted model and Σ_obs, i the observed surface density of GC candidates in the given bin.

We used the following method to estimate the uncertainty of the fitted parameters. It is sometimes called the method of support. It is based on the statistical likelihood-ratio test. Let p_max denote the vector of the best-fit free parameters and ℒ_max = ℒ(p_max) the maximum value of the likelihood function. Then the upper (lower) limit on the value of the jth free parameter, p_j, was estimated by maximizing (minimizing) p_j over the region of the parameter space satisfying the condition

$$\begin{array}{c}\hfill ln{\mathcal{L}}_{\max }-ln\mathcal{L}(\mathit{p})<0.5.\end{array}$$(15)

The examples in Appendix C and Sect. 6.2 demonstrate on artificial data that our methods are able to recover correctly the intrinsic parameters of the density of the GCSs.

3.4. Resulting parameters of the fits

For every galaxy, we eventually obtained density profiles of their GCSs and their fits from up to three types of data (ACS, FDS, spectroscopic, see Sect. 2). All the fitted parameters are provided in Table A.1. In Figs. A.1–A.19 we show plots of the observed density profiles together with the fitted models. For all galaxies, we show the background-subtracted profiles – this means that the fitted value of the surface density of contaminant sources was subtracted from both the measured and fitted profiles. We also show profiles without background subtraction for a few galaxies, namely for NGC 1352 (Fig. A.5), NGC 1387 (Fig. A.11) and NGC 1427 (Fig. A.15). In all plots showing the background-subtracted FDS profiles we also indicated by the dashed horizontal lines the fitted values of the background.

It was necessary to accept for every galaxy one final set of the parameters of the density profile of the GCS. We had to choose for every galaxy the values of the inner slope a, outer slope b and the break radius r_br. We did not make any final choice of the surface density of contaminant sources and the normalization of the density profile of the GCS – they are strongly influenced by the method of observation and were not necessary for the subsequent analysis.

The way we selected the final set of parameters of the GCS profile density was different for different galaxies. We chose the parameters derived from the spectroscopic data only if no other data were available for the given galaxy, because spectroscopic data can easily be degraded by geometric incompleteness of the survey. Whenever possible, we based the final parameters on the FDS or ACS data. If data from both surveys were available for a given galaxy, the strategy of accepting the final parameters was decided on the basis of a visual inspection of the surface density profiles. If the inner slope, a, was different in the FDS and ACS data, we accounted the difference to the problems with contamination by the light of the host galaxy in the FDS data (see Sect. 3.1) and preferred the ACS value. For many galaxies, the ACS data did not fully cover the outer parts of the broken power law profiles. Therefore, if the outer slopes of the profile came out differently, we adopted b from the FDS data. Regarding the break radii, if the FDS profiles appeared affected by the contamination by the light of the host galaxy (i.e., the inner slope was shallower for the FDS data than for the ACS data), we preferred the break radius estimated from the ACS data. If the break was close to the outer limit of the ACS data, we preferred r_br from the FDS data.

In the cases that a given parameter p_j appeared consistent between the two datasets for a given galaxy, we combined the measurements using the following form of the weighed average. It can account for the fact that our estimates of the uncertainty limits on p_j were asymmetric. Let Δ⁺p_j and Δ⁻p_j denote the upper and lower errorbars of the parameter p_j, respectively. Next, let us introduce the “joint Gaussian distribution”

$$\begin{array}{c}\hfill g(x,\overline{x},{\sigma }_{+},{\sigma }_{-})=\{\begin{array}{c}\frac{\sqrt{2}}{\sqrt{\pi }({\sigma }_{+}+{\sigma }_{-})}exp[-\frac{{(x-\overline{x})}^2}{2{\sigma }_{+}^2}]\text{for}x\ge \overline{x}\hfill \\ \hfill \\ \frac{\sqrt{2}}{\sqrt{\pi }({\sigma }_{+}+{\sigma }_{-})}exp[-\frac{{(x-\overline{x})}^2}{2{\sigma }_{-}^2}]\text{for}x<\overline{x},\hfill \end{array}\end{array}$$(16)

which is subsequently used on several occasions in this paper. We approximated the probability distribution of p_j by $g(p_j,\overline{p_j},{\mathrm{\Delta }}^{+}{p}_j,{\mathrm{\Delta }}^{+}{p}_j)$. Here $\overline{p_j}$ stands for the best-fit value of p_j. If g_FDS denotes the probability distribution for the FDS data and g_ACS that for the ACS data, then the function

$$\begin{array}{c}\hfill g_{\mathrm{tot}}(p_j)=g_{\mathrm{FDS}}(p_j)g_{\mathrm{ACS}}(p_j)\end{array}$$(17)

represents the total likelihood function of p_j. We got the final estimate of p_j as the argument p_{j, jointmax} maximizing the function g_tot(p_j). The uncertainty limits were obtained through the method of support, that is as the border values of the interval of p_j meeting

$$\begin{array}{c}\hfill ln{g}_{\mathrm{tot}}(p_{j,\mathrm{jointmax}})-ln{g}_{\mathrm{tot}}(p_j)<0.5.\end{array}$$(18)

For some galaxies, the GCS break radius was close to the border of the ACS data and the inner reliable limit of the FDS data. For such galaxies, we had only the fits of the inner slope from ACS data and the outer slope from the FDS data. In such cases, we opted for a “manual” method when we estimated the break radii and their uncertainties by visual inspection of both profiles in the background-subtracted surface density plots. That was, for example, the case for the profile of the blue subsample of the sources around NGC 1336, as shown in Fig. A.4. In these cases, when we accepted this subjective method, we stated very conservative estimates of the uncertainty limits.

The final parameters are listed in Table 2. The notes indicate the data or joining method that was used to obtain the value: f indicates the FDS data, h the HST ACS data, s the spectroscopic data, a the weighed average, and m the manual joining method. There are also some values marked by an asterisk. They indicate suspicious measurements. Those were identified by visual inspections of the observed and fitted profiles of surface density in Appendix A. The measurements were usually deemed suspicious if one of the parts of the broken power law profile was resolved by only one or two data points. This corresponds to the galaxies that have few GCs. Next, we marked all parameters of the density profile of NGC 1404 suspicious, because this galaxy is known to be undergoing a tidal interaction with NGC 1399, their GCSs overlap and we used spectroscopic data to fit the profile, that can suffer from incomplete spatial coverage of the spectroscopic surveys. Finally, we deemed suspicious the fitted values of b for the two galaxies with the most extended GCS, that is NGC 1399 and NGC 1316, because they might have been affected by the large-scale sensitivity variations of the FDS survey (Sect. 2.3). This is discussed in detail in Appendix D. It turned out that all of the suspicious measurements actually do not deviate noticeably from the scaling relations followed by the trusted measurements, as it will be shown in Sect. 5.

4. Calculation of a₀ radii

The a₀ radii were calculated on the basis of the Sérsic fits of the galaxies, and the estimates of their stellar masses (Su et al. 2021, see Sect. 2 for details). We were interested in two types of a₀ radii: those predicted by Newtonian gravity and those predicted by MOND gravity. Assuming that the galaxies are spherically symmetric, we could use the approximate analytic formulas for the density Sérsic spheres of Lima Neto et al. (1999; as updated by Márquez et al. 2000) to calculate the profiles of Newtonian gravitational acceleration, g_N(r). The profiles of MOND acceleration, g_M, were obtained by transforming g_N using the formula Eq. (2). We adopted the observationally motivated interpolation function ν:

$$\begin{array}{c}\hfill \nu (x)={[1-exp(-\sqrt{x})]}^{-1}\end{array}$$(19)

and value of a₀ = 1.2 × 10⁻¹⁰ m s⁻² (McGaugh et al. 2016; Li et al. 2018). Once the radial profiles of gravitational accelerations were known, we could find the a₀ radii numerically.

Finding the a₀ radii for the stacked galaxies required a more elaborate approach. First, we calculated the acceleration profile for each of the stacked galaxies individually. We then assigned to each stack a final acceleration profile that was calculated as the median acceleration profile of all objects contributing to the stack. Then we could solve for the a₀ radii.

In the following, we denote by r_a0, N the Newtonian a₀ radius and by r_a0, M the MOND one. In the case that we do not need to distinguish between them or the statements hold true for both, we denote the a₀ radius by r_a0.

Uncertainties on the a₀ radii were derived only from the uncertainties of stellar masses. These were tabulated by Su et al. (2021) for several stellar mass bins. We used linear interpolation to obtain the uncertainties of the stellar mass of individual galaxies. We neglected the uncertainties in distance. The total line-of-sight distance scatter for the galaxies in the Fornax cluster is about 0.5 Mpc (Blakeslee et al. 2009), which at the assumed distance of the center of the cluster of 20 Mpc, would result in the difference of stellar mass of 0.02 dex. This is negligible compared to the uncertainty in masses caused by the uncertainty in the mass-to-light ratio (Table 1).

For NGC 1399 we considered including the contribution of the hot intergalactic gas to the gravitational field of the galaxy. This could potentially influence the a₀ radius of the galaxy. The profile of cumulative mass of the hot gas for this galaxy was presented in Paolillo et al. (2002) and Samurović & Danziger (2006). It turned out that including the gas mass has only very little effect on the position of the a₀ radii. At the position of the MOND a₀ radius the cumulative mass of hot gas is just ∼10⁹ M_⊙, which is negligible compared to the mass of the galaxy (see Table 1). For this reason, we neglected the contribution of the gas mass in the gravitational field when estimating the a₀ radii.

The resulting a₀ radii are tabulated in Table 1. Some galaxies at the low-mass end of our sample, do not have a₀ radii. These are low-surface-brightness galaxies inside which the gravitational acceleration does not exceed a₀. Interestingly, the GCSs of such galaxies still can show broken power-law profiles as we see below.

5. Results

5.1. Structural parameters of GCSs for different GC subpopulations

We explored how the fitted values of the structural parameters of GCSs, a, b, a − b and r_br, differ for the total, red and blue GC populations. Inspection of plots of these parameters against the stellar masses of the galaxies revealed that the suspicious measurements identified in Sect. 3.4 often lied far from the reliable measurements, with the exception of the break radii. Therefore, all suspicious measurements were excluded from the subsequent analysis in this section.

For each structural parameter, we counted the number of galaxies for which the parameter is greater for the red subpopulation of GCs than for the blue subpopulation, the number of galaxies for which the situation is opposite, and the number of galaxies for which the parameter is consistent for the blue and red GC populations. The consistency means that the uncertainty intervals of the measurements overlap. The results are listed in Table 3. Given that we do not have all parameters for all galaxies, the last column of this table indicates for each structural parameter the total number of galaxies that this comparison is based on.

The inner slope of the GCS density profile, a, shows the most variate behavior. It can be both smaller or greater for the red GCs than for the blue GCs, or the two GCs populations can have it consistent. All these cases occur approximately in the same number of galaxies. This result could have been influenced by the problems with the contamination by the light of the host galaxies. The b parameter is usually consistent for the red and blue GC populations, but in some cases the parameter is higher for the blue GCs. The difference of the slopes, a − b, is consistent for the red and blue GCs for all but one galaxy, for which a − b is greater for the red GCs (i.e., the break is more pronounced). The values of r_br are in a roughly equal number of galaxies either consistent for the red and blue GC populations or they are lower for the blue population.

Next, we estimated the mean and intrinsic scatter of the structural parameters for all GC populations. We assumed that the intrinsic distribution of every parameter p_j is Gaussian:

$$\begin{array}{c}\hfill f(p_j)=\frac{1}{\sqrt{2\pi }{\sigma }_{\mathrm{int},j}}exp[-\frac{{(p_j-\overline{p_j})}^2}{2{\sigma }_{\mathrm{int},j}^2}],\end{array}$$(20)

where $\overline{p_j}$ is the mean and σ_int, j the intrinsic scatter of the distribution. We took into account the asymmetric uncertainty intervals of the structural parameters. The probability distribution of the parameter p_j of the ith galaxy, p_j, i, was modeled as $g_{\mathrm{mes},j,i}(p_{j,i},\overline{p_{j,i}},{\sigma }_{p_{j,i}}^{+},{\sigma }_{p_{j,i}}^{-})$, where the function g was defined by Eq. (16). The best-fit values of $\overline{p_j}$ and σ_int, j were found by maximizing the likelihood function:

$$\begin{array}{c}\hfill \mathcal{L}(\overline{p_j},{\sigma }_{\mathrm{int},j})={\displaystyle \prod _{i=1}}[f(p_{j,i})\ast g_{\mathrm{mes},j,i}(p_{j,i})](p_{j,i}),\end{array}$$(21)

where the symbol denotes convolution. The uncertainty limits of $\overline{p_j}$ and σ_int, j were found through the method of support. In order to avoid numerical difficulties, we required the intrinsic scatter to be at least 0.01. The results are stated in Table 4.

First, the table reveals the typical values of the structural parameters: a = −1.7, b = −3.4 and a − b = 1.7. The parameters a, b and their difference a − b do not differ substantially between all GCs, red GCs and blue GCs. The values are consistent within 1σ. The intrinsic scatters of the parameters are consistent with each other for the different GC populations too. It is worth noting that the intrinsic scatter of the prominence of the break, that is a − b, is also consistent with being zero for the total and blue GC populations. This means that the intrinsic scatter is smaller than the measurement uncertainties of the individual data points.

We also inspected the statistical distribution of the differences of the structural parameters for red and blue GCs in individual galaxies, that is p_{j, i, red} − p_{j, i, blue}. We derived the mean and intrinsic scatter as before. The results are summarized in Table 5. It confirms that the structural parameters of the red and blue GC populations are in average the same within 1-2σ uncertainty limits. The outer slope b is at the border of the 2σ uncertainty limit, which suggests that the outer slope might be systematically steeper (i.e., more negative) for the red GCs. The break radii for the red and blue GC populations seem to be remarkably consistent, differing typically just by 0.3 kpc.

In summary, we found at most marginal evidence that the structural parameters of the GCSs would depend on the color of the GCs. Following the rule of Occam’s razor, we consider hereafter no difference between the structural parameters of the GC subpopulations.

5.2. Correlations of the structural parameters of GCSs

We were interested in how the fitted structural parameters of GCSs, a, b, a − b and r_br, correlate with the stellar mass, effective radius, Sérsic index and the Newtonian and MOND a₀ radii of the host galaxy. We exploited that through Spearman’s correlation coefficient and its p-value. The p-value expresses the probability that the two quantities under consideration actually do not correlate and the observed amount of correlation is there because of a coincidence. We comment here only on the pairs of quantities which correlate at least at the 5% confidence level (i.e., their p ≤ 0.05).

We found that the parameter b correlates significantly with both Newtonian and MOND a₀ radii. The correlations are shown in Fig. 2. The red numbers in the corners of the tiles of the figure indicate the p-value of the correlations. The open symbols indicate the suspicious measurements.

Fig. 2. Statistically significant correlations of the outer slope of the density profile of the GCSs b with the properties of their host galaxies. Top: the MOND a0 radius. Bottom: the Newtonian a0 radius. The open symbol indicates a suspicious measurements. The red numbers in the corners indicate the p-value of the Spearman’s correlations.

It is noteworthy that the parameter b correlates better with the a₀ radii than with the galaxy stellar mass. The p-value of the correlation with the latter is 0.03, after removing of the galaxies for which the a₀ radii do not exist. In contrast, the p-value of the correlation of b with the a₀ radii is 0.008 (0.009) for the Newtonian (MOND) case. This suggests that there is a connection between the distribution of stars in the galaxy and the distribution of GCs beyond the break radius, which typically exceeds the effective radius of the galaxy in our sample twice.

The break radii correlate significantly with all the considered characteristics of the galaxies, as shown in Fig. 3. We note, however, that only the correlation with the a₀ radii is close to a one-to-one relation, as we describe in more detail in Sect. 5.3, in agreement with the finding of BSR19. The mean of the ratio of the break radius and the effective radius r_br/R_e = 2.4 with a root-mean-square scatter of 0.7. The ratios r_br/r_a0 are discussed in Sect. 5.3.

Fig. 3. Statistically significant correlations of the break radii of GCSs with the properties of their host galaxies. From top to bottom: stellar mass, effective radius, Sérsic index, MOND a0 radius, Newtonian a0 radius. The open symbols indicate suspicious measurements. The red numbers in the corners indicate the p-value of the Spearman’s correlations.

We also explored whether the structural parameters a, b, a  −  b and r_br correlate with each other. We found only two significant correlations, that are not surprising in the light of what has been said above. The break radius correlates with the b parameter, as it can be expected because we already found above that b correlates with the a₀ radius, which in turn correlates with the break radius. Next, we found that the difference a − b correlates with b. This is again not too surprising given that the absolute value of b is typically larger than the absolute value of a (Table 4).

5.3. Equality of break radii and a₀ radii

Here we come to the main part of the paper – the comparison of the break radii of GCSs with the a₀ radii of their host galaxies. The work of BSR19 pointed out their approximate match, nevertheless their sample contained only galaxies with a relatively narrow range of stellar masses and therefore also of the a₀ radii. Moreover, they investigated only the total populations of GCs, such that the shapes of profiles of the red and blue GCSs, including the break radii, remained unclear. The sample investigated here rectifies these insufficiencies.

The break and a₀ radii are compared in Fig. 4. The top row corresponds to the a₀ radii calculated from the MOND gravity, and the bottom row to those calculated assuming Newtonian gravity. The three columns, from left to right, correspond to the population of all, red and blue GCs, respectively. The squares represent the data from the sample from the present paper, while the circles those from BSR19. The open squares indicate the suspicious measurements. The diagonal black dashed lines indicate the one-to-one relation. The break radii of the galaxies that do not have a₀ radii are shown by the vertical dashed lines. In these galaxies, the acceleration generated by their stars, either as predicted by the Newtonian or MOND gravity, is less than a₀ everywhere in their extents. The galaxy NGC 1399 is plotted twice in each panel – one measurement comes from the data analyzed in this paper, the other comes from the paper BSR19, where it was derived from spectroscopic data. The match of the break and a₀ radii for NGC 1399 became worse in this paper (see Sect. 6 for more details).

Fig. 4. Relation of the a0 radius versus the break radius of the GCS for each galaxy. The upper row corresponds to the MOND a0 radius, the bottom row to the Newtonian a0 radius. The three columns in the horizontal direction correspond to the whole, red, and blue GC (sub)population, respectively. The squares in the first column mark the galaxy sample inspected in this paper, the circles mark those from BSR19. The empty squares denote the suspicious measurements. The diagonal dashed black lines represent the one-to-one relation. The vertical gray dashed lines mark the galaxies that do not have a0 radii.

The figure demonstrates a good match of the a₀ radii and break radii. The match is better for the a₀ radii calculated from MOND. This holds true regardless of the separation of the GCs in red and blue subpopulations. It should be pointed out that an exact match cannot be expected because of the tidal interactions between the galaxies and galaxy mergers in the dense environment of a galaxy cluster. The tidal interactions can even lead to a loss of GCs from the galaxies or to a transfer of GCs from one galaxy to another (Bekki et al. 2003). Our data indeed provide an observational indication that galaxy interactions influence the break radius, see Sect. 6.3.

Figure 5 gives us an alternative view of the same data. It shows the ratios of the break and a₀ radii plotted against the break radii. The estimates of the expected values and uncertainty limits of the ratios were based on the fact that the distribution h_z(z) of the variable z, that is a ratio of two independent variables x and y, z = x/y, can be calculated from the known distribution functions of x, h_x(x) and of y, h_y(y) using the so-called ratio distribution formula:

$$\begin{array}{c}\hfill h_z(z)={\displaystyle {\int }_{-\infty }^{\infty }}|y|h_x(zy)h_y(y)\mathrm{d}y.\end{array}$$(22)

Fig. 5. Demonstration of the approximate equality of the break radii of the GCSs and the a0 radii of their host galaxies. First row: the ratio of the break radius and the MOND a0 radius as a function of the break radius. The vertical axis is in the base-2 logarithmic scale, while the horizontal axis in the base-10 logarithmic scale. The horizontal line indicates the ratio of one (the break and a0 radii are equal). Second row: the same as the first row but for the a0 radii calculated for the Newtonian gravity. Columns: from left to right, the rows show the data derived for all GCs and for the red and blue subpopulations of GCs.

The distributions of the break and a₀ radii for every individual galaxy were again approximated by Eq. (16). We estimated the expected value and the 1σ uncertainty limits of r_br/r_a0 as the 18, 50 and 84th percentiles of the ratio distribution.

The top row of Fig. 5 refers to the a₀ radii calculated in the MOND way, the bottom row to those calculated in the Newtonian way. The points that lie at the top border were shifted downward because they lie out of the displayed radial range, but their error bars are displayed correctly. We could make several observations from these plots. 1) The break and MOND a₀ radii agree with each other within a factor of two for most galaxies – either in terms of the most likely values or within 1-2 uncertainty limits. This applies to all GC populations. 2) For the total GC population, there is a hint of a correlation of the ratio r_br/r_a0 with r_br if the break radius is greater than about 20 kpc. In this region, however, most data points come from the spectroscopic data analyzed in BSR19, which could be biased by systematic errors. 3) For the blue GC population, there is a correlation of the ratio r_br/r_a0 with r_br for the whole galaxy sample. 4) The Newtonian ratios r_br/r_a0 seem to be systematically more offset from unity than the MOND ratios. We inspect the points 1) and 4) more in detail below and the points 2) and 3) in Appendix E.

Let us denote ζ = r_br/r_a0. We aimed to fit the distribution of log₁₀ζ by a normal distribution with the mean of log₁₀μ and a scatter of σ_int. Specifically, we wanted to fit the distribution of ζ by the lognormal distribution:

$$\begin{array}{c}\hfill g_{\mathrm{LN}}(\zeta )=\frac{1}{\sqrt{2\pi }ln(10){\sigma }_{\mathrm{int}}\zeta }exp\{-\frac{{[{log}_{10}(\zeta )-{log}_{10}(\mu )]}^2}{2{\sigma }_{\mathrm{int}}^2}\}.\end{array}$$(23)

The best-fit parameters and their uncertainties were found through a likelihood function analogous to Eq. (21). This was applied to the sets of the total population and the red and blue subpopulations of the Fornax cluster galaxy sample investigated here. In addition, we created a union sample of the galaxy set of BSR19, with NGC 1399 excluded, and the sample of total GC populations of the Fornax galaxies. This was repeated for the a₀ radii calculated in the Newtonian and MOND ways. The results are listed in Table 6.

The table indeed shows that the MOND r_br/r_a0 ratios are indeed close to one for all GC (sub)populations. The intrinsic scatter is around 0.28 dex (i.e., the factor of 1.9). The Newtonian r_br/r_a0 ratios are somewhat offset from unity, with a mean of about 1.4. The intrinsic scatter is similar to the MOND case.

In Appendix E, we fit the relation of the break radius and a₀ radius by a linear function with a lognormal intrinsic scatter. The fitted value of the slope is not consistent with one. Nevertheless, it is possible to argue that the deviation of the slope from one is caused by a few outliers. The data inspected in this work thus indicate just the fact that the break and a₀ radii agree within a factor of two.

6. The break in the GCS of NGC 1399 under scrutiny

Here we explore in detail the profile of the GCS of NGC 1399, the central galaxy of the Fornax Cluster. Apart from that, we chose this galaxy because is has the richest GC system in our sample, and we moreover have multiple datasets for it (FDS, ACS and spectroscopic data). Its GC density profile has been fitted by a broken power law already in BSR19.

6.1. Comparison of profiles extracted from different datasets

We analyzed the positions of GCs for this galaxy from three sources: the photometric catalogs of FDS and ACS and the spectroscopic data of Chaturvedi et al. (2022) and Fahrion et al. (2020). In addition, even another spectroscopic dataset (Schuberth et al. 2010) was fitted by a broken power law in BSR19 (and without the uncertainty limits already in Bílek et al. 2019b). The comparison of all the profiles is shown in Fig. A.12.

The inner slope of the broken power law is virtually the same for all the datasets analyzed here. This suggests that the presence of NGC 1399 in the center of the GC system did not affect our ability to detect the GCs. This might be counterintuitive, given that this is the brightest galaxy analyzed here. This is probably because the relation between magnitude and surface brightness of elliptical galaxies shows a peak at intermediate magnitudes (Graham & Guzmán 2003). The match of the inner slope of the profile derived from the spectroscopic GCs of Chaturvedi et al. (2022) and Fahrion et al. (2020) with the photometric samples indicates that it has a good spatial coverage near the center of the galaxy.

The outer slope and the break radius of the spectroscopic data do not agree with the FDS data that well. While the outer slopes for the total population of GCs agree within the 2σ uncertainty limits, the break radii do not (Table A.1). This might indicate either issues with a geometric incompleteness of the spectroscopic sample or an imprecise estimation of the surface density of contaminant sources when extracting the profile from the photometric data. Indeed, the GCS of NGC 1399 is extended and the extraction of the data might have be affected by the variations of the sensitivity of OmegaCAM near the edges of its FOV (Sect. 2.3). The outer slope of the FDS profile is rather bumpy, whereas that of the spectroscopic sample looks cleaner. This is probably because of contaminants in the FDS sample, such as background galaxy clusters. These would be excluded from the spectrocopic sample on the basis of radial velocity.

The profile analyzed in BSR19 deviates most from the other profiles. We attribute this to a geometric incompleteness of the survey (see Fig. 1 in Schuberth et al. 2010). The position of the break radius is however not affected too much: BSR19 found $7.3_{-0.4}^{+2^{\prime }}$, while we accepted here 11.0′±0.6 (Table 2). The match of the break radius with the a₀ radius, located at about 2.4′, got worse in this paper. We note that the surface density profile of the GCs from the catalog of Schuberth et al. (2010) was fitted by a broken power law already in Samurović & Danziger (2006). Their results are in good agreement with those of BSR19.

6.2. Dependence of the profile parameters on radial velocity cuts

The large number of GCs around NGC 1399 allowed us to explore how the radial profile of the density of the GCs depends on the line-of-sight velocity of the GCs with respect to the center of the galaxy. We assumed a radial velocity of NGC 1399 to be 1424.9 km s⁻¹. The profiles for different velocity cuts are displayed in Fig. 6. The fitted values of the parameters of the broken power-law profiles are presented in Table A.3.

Fig. 6. Radial profiles of surface density of GCs of NGC 1399 in several bins of radial velocity with respect to the center of the galaxy.

The figure shows that the break radius shifts toward the center of the galaxy for the GCs that have larger velocities with respect to the galaxy. The outer slope, b, remains constant. The inner slope, a, becomes steeper for GCs having lower radial velocities with respect to the galaxy center. As a consequence, the GC density profile is almost a simple, that is unbroken, power law for the lowest radial velocity bin.

As the first step toward the understanding of this observation, we constructed simple dynamical models of the GCS of NGC 1399 and explored how the profile of GC density depends on the chosen radial velocity bin. To this end, we solved the spherical Jeans equation, which gave us the profile of velocity dispersion of the GCS, see Bílek et al. (2019b) for details. We assumed the dark matter halo and the stellar mass-to-light ratio that were derived in that paper from the best-fit isotropic Jeans model of the GCS of this galaxy. For our present model, we assumed the same broken-power-law profile of the density of the GC system as in Bílek et al. (2019b). We solved the Jeans equation for three different choices of the anisotropy parameter: β = 0 (isotropic), β = 1 (purely radial) and β = −3 (highly tangential). The β parameter specifies the typical shapes of orbits of the GCs around the galaxy. Having solved the Jeans equation, we generated a three-dimensional model of the GC system. Each GC was randomly assigned a position and velocity according to the assumptions described above. The velocities were drawn from an ellipsoidal Gaussian distribution, according to the solution of the Jeans equation and the assumed anisotropy parameter. There were 10⁴ GCs in each model. We eventually created a catalog of the projected positions and line-of-sight velocities of the artificial GCs. That was analyzed like the real data.

The results are shown in Fig. 7. For the isotropic and radial models, the position of the break indeed depends on the radial velocity of the GCs with respect to the galaxy, such that the GCs that are slower have their break radius further away from the galaxy, in agreement with the observed data. The tangential model does not show any obvious dependence of the break radius on the GC velocity range. This hints at a radial or isotropic anisotropy of the real GCS of NGC 1399. Next, the slopes of the broken power law depend on the inspected range of radial velocities of GCs in all models. The difference is most pronounced for the inner slopes of the radial and tangential models. The inner slopes of the radial model become more shallow toward the low radial velocities of GCs. The trend is the opposite for the tangential model. The tangential model thus resembles the real NGC 1399 in this regard.

Fig. 7. Same as in Fig. 6, but for modeled data. The models are Jeans spherical models with the indicated types of anisotropy parameters. The models were based on the gravitational potential of NGC 1399 derived in Bílek et al. (2019b) and the parameters of the GCS used in that work. The vertical dashed line at 30′ is an eye-guide to facilitate the reading of the variations of the break radius.

We explored many other values of the anisotropy parameters, allowing them even to be a function of the galactocentric distance. We were never able to fully reproduce the trends observed in NGC 1399. Most notably, we never got close to equalizing the inner and outer slopes of the GC density profile that is observed for the GCs that have low radial velocities with respect to the center of NGC 1399. It might be necessary to model the GC system by several GC populations with different spatial distributions and anisotropy parameters. They would correspond to GCs formed in situ and accreted from possibly several galaxies. This more realistic modeling is beyond the scope of this paper. It should also be pointed out that in all of our models, including the non-isotropic ones, we used a gravitational potential that was derived assuming an isotropic and spherical GCS.

As a side note, the values of the parameters recovered from the projected positions of the artificial GCs agreed, within the uncertainty limits with the parameters that were used to generate the three-dimensional positions of the GCs. This demonstrates the correctness of our methods.

6.3. Dependence of the break radius on azimuth: elongation toward the interacting neighbor NGC 1404

We inspected also how the break radius varies in azimuth in NGC 1399. We defined eight sectors of the same angular widths centered on the galaxy. The first sector was centered on the western photometric major semi-axis of NGC 1399. This semi-axis points to the west almost exactly. We initially experimented with the GC candidates selected from the FDS data for this exercise, but it turned out that the profiles contained spikes. By inspecting images of the sky from the Digitized Sky Survey at the position of the spikes in the Aladin software (Bonnarel et al. 2000; Boch & Fernique 2014), it seemed that the spikes were caused by clusters of galaxies in the background, that is by contaminating sources. This is why we eventually decided to utilize the spectroscopic data. This also is the reason why we could not inspect in detail the second largest GCS in the Fornax cluster that is possessed by NGC 1316.

The fitted parameters of the individual sectors are presented in Table A.2. The found break radii are plotted as a function of azimuth in Fig. 8. The azimuth 0°coincides with the photometric western major semi-axis of NGC 1399; the azimuth 90°coincides with the northern minor semi-axis. We can see that the break radii form approximately an ellipse. It is interesting that the ellipse is inclined substantially with respect to the stellar body of the galaxy. Instead, one of the major semi-axes of the ellipse formed by the break radii points toward the galaxy NGC 1404. It is well known that both galaxies are tidally interacting (Bekki et al. 2003; Sheardown et al. 2018). This suggests that galaxy interactions can influence the distribution of GCs, and thus the values of break radii.

Fig. 8. Break radii of the GCS of NGC 1399 in different sectors centered on the galaxy. The radial coordinate is in arcminutes. The red point indicates the position of NGC 1404.

7. Interpretation

BSR19 proposed several potential explanations for the existence of the breaks in the profiles of the densities of GCS and of the coincidence of the break radii with the a₀ radii. In the following subsections, we develop these ideas one step further and discuss them in the light of the new data. New ideas are added. Some of the interpretations are specific for a MOND universe, because in MOND, the constant a₀ enters naturally in many phenomena. It is in part because a₀ marks the boundary between the two regimes of the law of gravity or inertia. Some of our interpretations are applicable also in the ΛCDM cosmology. Some of the explanations are based on the validity of the radial acceleration relation, which is an empirical fact. In MOND, the radial acceleration is a trivial implication of the theory, while the ΛCDM cosmology still is finding its way toward its full explanation through the galaxy formation theory (e.g., Di Cintio & Lelli 2016; Santos-Santos et al. 2016; Navarro et al. 2017).

7.1. Consequence of two regimes of gravitational potential and of the accretion of GCs in mergers

Supposing Newtonian gravity, the gravitational field in the inner region of a high-surface-brightness galaxy is dominated by the contribution of baryons. On the other hand, far from the galaxy center the gravitational field is dominated by the contribution of the dark matter halo. Given that according to the radial acceleration significant amounts of dark matter are needed only beyond the a₀ radius, the strengths of the gravitational fields generated by the baryons and the dark matter halo have to be equal roughly at the a₀ radius. In the inner region of the galaxy, the gravitational potential is steep, since it can roughly be approximated by that of a point mass that represents the baryonic component. In the outer region, the potential is more shallow. According to the radial acceleration relation it can be approximated by a logarithmic potential.

A large fraction of the baryonic mass of massive galaxies is expected to be gained by accretion of smaller galaxies. These bring also their GCs into the system. There are multiple pieces of evidence that a large fraction of blue GCs gets into giant galaxies through the accretion of dwarf galaxies (Côté et al. 1998; Hilker et al. 1999).

The GCs that are brought in with the satellites are tidally stripped from the satellites preferentially when the satellites are close to their orbital pericenter with respect to the host. Various accreted satellites have various pericentric velocities. If the stripped GCs have a large enough velocity at the moment of stripping, they reach the outer shallow part of the gravitational potential of the host and spread over a large range of apocentric distances. In contrast, if the GCs have a low radial velocity with respect to the host when stripped, they reach in the apocenters of their orbits only the inner steep part of the gravitational potential. Therefore their apocentric distances span only a narrow range (illustrated in Fig. 2 of BSR19). A break in the profile of the GC system is then expected at the border between the steep an shallow part of the gravitational potential of the host, that is near the a₀ radius.

The same reasoning can be applied to MOND. The only difference is that the change of the slope of the gravitational potential in the inner and outer regions of the galaxy is not because of a dark matter halo, but because of the strong and weak field regimes of MOND.

If such an interpretation was true, then the breaks in the density of GCSs would be very useful for investigation of dark matter halos under the assumption of the ΛCDM cosmology. The breaks would mark the radius at which the gravitational field changes from the baryon-dominated to the dark-mater dominated regime. This would allow estimating the effective radius of the dark matter halos. It is already known that the masses of dark matter halos of galaxies can be estimated from the number of GCs that the galaxies have (Spitler & Forbes 2009; Harris et al. 2015).

It would be ideal to explore if this mechanism of formation of breaks in GCS density profiles actually works through simulations. This is beyond the scope of this paper. Here we instead resort to looking for observational signatures of this scenario.

For giant galaxies, most blue GCs are deemed to be accreted while many of the red GCs are deemed to be formed in situ (Côté et al. 1998; Hilker et al. 1999). Under this scenario of the origin of the breaks, we would thus expect that the breaks would be more pronounced for the blue GCs. Our data do not show this – the broken number density profiles of GCSs are observed for both red and blue subpopulations and the magnitude of the break, a − b, does not seem to be different for the two subpopulations. A quantification through simulations however remains desirable.

Next, the galaxy sample presented here contains, unlike the sample of BRS19, also relatively low-mass galaxies, with the masses of the Small Magellanic Cloud. Low-mass galaxies are expected to form mostly by in situ star formation, without a substantial growth through mergers. From this point of view, it is unexpected that we detected breaks also in the low-mass galaxies. On the other hand, the early-type galaxies of our sample, that are supported mainly by velocity dispersion, might have richer merging histories than the typical dwarf galaxies, which are rotating. It is known, however, that dwarf galaxies in cluster environments can transform from disky, rotating dwarfs to spheroidal, dispersion dominated dwarfs via tidal stirring (e.g., Mayer et al. 2001a,b). In total, we found tentative evidence against the origin of the breaks of GCSs through this mechanism.

7.2. MOND external field effect and stripping of dark halos

In MOND, the dynamics of an object follows the radial acceleration relation only if the object is isolated. If the object is located in an external gravitational field (for example a galaxy in a galaxy cluster), the apparent enhancement of gravity compared to the Newtonian gravity diminishes because of the nonlinearity of the theory. Once the strength of the external field surpasses the constant a₀ substantially, then the dynamics of the object behaves as in Newtonian dynamics without dark matter. This is called the external field effect (Milgrom 1983c; Bekenstein & Milgrom 1984). If interpreted in the Newtonian way, an object that is exerted to a gradually stronger external field behaves as it was losing its dark matter halo. Observational evidence for the external field effect has been reported (McGaugh & Milgrom 2013; McGaugh 2016; Caldwell et al. 2017; Chae et al. 2020), even if the galaxies in clusters might be an exception (Freundlich et al. 2022).

Galaxy clusters assemble by accreting individual galaxies from less dense environments. The external field effect then reduces the gravitational fields of the galaxies beyond their a₀ radii. At smaller radii, the gravitational field remains Newtonian, as it was when the galaxies were far from the cluster. Therefore, one expects that a break will develop at the a₀ radius.

We made a simple model to explore the impact of this process. The GCs were initiated as orbiting a point source with a mass of 10¹⁰ M_⊙ on circular orbits and a zero external field was assumed. Then the model was evolved for 10 Gyr, increasing the magnitude of the external field linearly with time, such that the external field eventually reached the value of 2a₀, a value typical for the cores of galaxy clusters (Milgrom 2008). The magnitude of the gravitational force in the presence of an external field was calculated using the so-called 1-D approximation of QUMOND (Famaey & McGaugh 2012).

The trajectories of the modeled GCs are shown in Fig. 9. The horizontal dotted black line indicates the MOND a₀ radius. While the GCs below the a₀ radius stay at their initial orbits, the distant GCs recede as the consequence of the external field effect. This causes a dilution of the GCS beyond the a₀ radius, such that the profile bends down.

Fig. 9. Impact of the increasing of the external field on the trajectories of GCs that were originally on circular orbits of different radii, assuming the MOND gravity. The horizontal line indicates the MOND a0 radius.

It is possible to estimate the impact of the external field effect on the distribution of GCs analytically. If the gravitational potential of the host remains spherical when it is changing, then there is no tidal torque acting on a GC moving on a circular orbit. The angular momentum of the GC is then conserved. The angular momentum of a GC orbiting the host galaxy well beyond the a₀ radius when the galaxy is far from the cluster, is $r_0\sqrt[4]{GM{a}_0}$. Once the GC and its host appear in a strong external field, compared to a₀, the gravitational field of the host becomes Newtonian and thus the angular momentum of a GC on a circular orbit becomes $r_1\sqrt{\frac{\mathit{GM}}{r_1}}$. From the conservation of angular momentum we get:

$$\begin{array}{c}\hfill r_1=r_0^2\sqrt{\frac{a_0}{\mathit{GM}}}.\end{array}$$(24)

If the Newtonian gravity is assumed, galaxies that enter galaxy clusters can reduce their gravitational fields by stripping of their dark matter halos in the consequence of tidal interactions with the other galaxies (Lee et al. 2018; Mitrašinović 2022). This can mimic the external field effect to a certain degree.

Even this mechanism of creating breaks in GCS density profiles through the reduction of the gravitational fields of galaxies is not perfect. The central galaxies of clusters, that spent their whole lives in the clusters, are not expected to reduce their gravitation fields. Yet NGC 1399 shows a break in the GCs density profile. Moreover, given that the intensity of the external field or tidal stripping is probably different for every galaxy, it is not clear why the galaxies would end up with a relatively narrow range of the external slope of the GCSs, b. Finally, some of the galaxies investigated in BSR19 were isolated (e.g., NGC 821 and NGC 3115), and thus never had an opportunity to develop a break through the external field effect or the tidal stripping of their halo.

7.3. Two regimes of dynamical friction

When a GC orbits its host galaxy, it attracts the stars, gas, or dark matter particles of the host galaxy gravitationally and gives them kinetic energy and angular momentum. As consequence of the laws of conservation of these quantities, the orbital angular momentum and energy of the GCs decrease. This is called dynamical friction. Dynamical friction manifests itself as a force acting on the GC against the direction of the velocity of the GC with respect to its environment. In Newtonian gravity, the magnitude of the dynamical friction force can be estimated by Chandrasehkar’s formula (Chandrasekhar 1943):

$$\begin{array}{cc}& F_{\mathrm{DF},\mathrm{NWT}}=\frac{2\pi ln\mathrm{\Lambda }{G}^2\rho m^2}{{\sigma }^2{X}^2}[\mathrm{erf}(X)-\frac{2X}{\sqrt{\pi }}exp(-X^2)],\hfill \\ \hfill & X=\frac{v}{\sqrt{2}\sigma },\hfill \end{array}$$(25)

where m stands for the mass of the GC, v its velocity with respect to its local environment, ρ the density of the local environment and σ the velocity dispersion of the environment. The expression lnΛ is called the Coulomb logarithm. Its value depends on the exact configuration of the problem under consideration, but it is on the order of a few.

A MOND analog of Chandrasekhar’s formula was proposed by Sánchez-Salcedo et al. (2006) on the basis of heuristic arguments and theoretical results by Ciotti & Binney (2004):

$$\begin{array}{c}\hfill F_{\mathrm{DF},\mathrm{MOND}}^{\mathrm{WF}}=\frac{a_0^2}{\sqrt{2}{a}^2}{F}_{\mathrm{DF},\mathrm{NWT}}.\end{array}$$(26)

Here a denotes the gravitational acceleration exerted by the host galaxy on the GC. Sánchez-Salcedo’s formula has recently been verified by simulations of GCs orbiting ultra-diffuse galaxies by Bílek et al. (2021). It is however supposed to work only in the weak-field regime of MOND, that is for a ≪ a₀. If a ≫ a₀, the dynamical friction force should reduce back to the one given by Chandrasekhar’s formula.

Therefore here we heuristically propose a universal formula for dynamical friction in MOND:

$$\begin{array}{c}\hfill F_{\mathrm{DF},\mathrm{MOND}}=\xi \left(\frac{a}{a_0}\right)F_{\mathrm{DF},\mathrm{NWT}},\end{array}$$(27)

where ξ(x) is a function satisfying the limit behavior:

$$\begin{array}{cc}& \xi (x)\to 1\text{for}x\gg 1,\hfill \\ \hfill & \xi (x)\to \frac{1}{\sqrt{2}{x}^2}\text{for}x\ll 1.\hfill \end{array}$$(28)

It was proposed in BSR19 that these two regimes of dynamical friction in MOND might give rise to the observed breaks in the GCS density profiles at the a₀ radii. For the purpose of the exercise below, we arbitrarily chose

$$\begin{array}{c}\hfill \xi (x)=1+\frac{1}{\sqrt{2}{x}^2}.\end{array}$$(29)

We explored if dynamical friction can influence the distribution of GCs noticeably through simplistic models of the dynamics of the GCs of NGC 1399 and NGC 1373. These two are representatives of the most and the least massive galaxies in our sample, respectively. We consider models assuming either Newtonian gravity with dark matter or the MOND gravity.

The distribution of stars of the galaxies was modeled by Sérsic spheres described by the parameters stated in Table 1. It also was important to include the circumgalactic media in our models given the appearance of ρ in Eq. (25). The break radii of the GCSs in our galaxy sample exceed 2-3 times the effective radii of the host galaxies and thus the contribution of the circumgalactic gas to the local baryonic density might be substantial. The masses of circumgalactic media are generally difficult to determine. They could possibly be comparable to the stellar masses of galaxies (Gupta et al. 2012; Werk et al. 2014; Keeney et al. 2017; Das et al. 2019). Here we modeled the distribution of gas around NGC 1399 according to the observed density of hot gas presented in Paolillo et al. (2002) in their Fig. 17. As central cluster galaxy, NGC 1399 is not expected to contain much cold gas (Yoon & Putman 2017; Burchett et al. 2018). The extended gaseous component of NGC 1373 was modeled by a Sérsic sphere whose mass and Sérsic index are the same as of the stellar component but the effective radius is four times larger.

For the Newtonian models, we assumed the NFW dark matter halo derived for NGC 1399 by Bílek et al. (2019b). It was deduced from an isotropic model of the kinematics of the GCS of the galaxy. For NGC 1373, we assumed a NFW dark matter halo of the virial mass of 10^11.35 M_⊙ and a scale radius 13.4 kpc according to the stellar-to-halo mass relation (Behroozi et al. 2013) and the halo mass-concentration relation (Diemer & Kravtsov 2015).

In the MOND models, it was necessary to consider the external field effect coming from the Fornax cluster. For the MOND model of NGC 1373, we included the external field effect such that the galaxy was assigned the Newtonian gravitational field produced solely by the stars of the galaxy. For NGC 1399, as for a central galaxy of the cluster, a zero external field was assumed.

The velocity dispersion that appears in Eq. (25) was obtained by solving the spherical Jeans equation with a zero anisotropy parameter. All material was approximated as consisting of a single population of collisionless particles. Bílek et al. (2021) found that for GCs on circular orbits, the value of the Coulomb potential is around 10 and it is around 3 for radial orbits. In the following calculations, we assumed the value of 10. The GCs in our models had a mass of 10⁵ M_⊙. Their orbits were integrated for 10 Gyr.

For NGC 1399, we found that neither for Newtonian or MOND gravity dynamical friction is able to affect the orbits of GCs noticeably. We explored various orbital eccentricities of the GCs and different apocentric distances in the range between 10 and 200 kpc. Dynamical friction started to influence the orbit of the modeled GC noticeably only after the mass of the GC was increased to at least about 10^8 − 9 M_⊙, depending on the orbit shape. Dynamical friction is thus rather expected to affect the spatial distribution of the satellite galaxies. With the mass of the GC of 10⁵ M_⊙, the situation did not change for the MOND model even if we multiplied the density of matter of the galaxy by a factor of 30 to model the hypothesized cluster baryonic dark matter (Milgrom 2008).

For NGC 1373, we explored various orbital eccentricities of the GCs and different apocentric distances in the range between 0.6 to 8 kpc. The friction came out to be the strongest near the center of the galaxy. Radial orbits were affected more than circular. The GCs that moved within about one effective radius of the galaxy eventually settled in the center of the galaxy, both in the Newtonian and MOND models. The central deficit of GCs was immediately balanced by sinking of the GCs that originally moved on bigger orbits. As result, generally no substantial breaks in the GCS density profiles were induced. We show in Fig. 10 an example of orbits that create a very strong break, compared to the other orbits we explored. The break is interestingly located near the a₀ radius. The figure shows the initial versus the final apocentric distances of GCs that move in the Newtonian model and were launched with velocities that are equal to 0.1 of the local circular velocity. The sinking of the GCs caused a mild bend of the density profile at around 2 kpc, which is where the break in the real galaxy is observed. The black dotted line indicates the one-to-one relation. Once the GCs were put on circular orbits, the break stopped being noticeable. In a real galaxy, where the GCs orbits have various eccentricities, the break induced by dynamical friction would thus be even milder than in Fig. 10. We verified that there is no change in these conclusions even for the intermediate-mass galaxy NGC 1379.

Fig. 10. Effect of dynamical friction on the apocentric radii of GCs that initially move in apocenter with a velocity that is equal to 0.1 of the local circular velocity.

This demonstrates that dynamical friction on a GC is alone probably not responsible for the observed breaks in GCS density profiles. In low-mass galaxies, it is important for the dynamics of GCs within the stellar body of the galaxy, but the contribution to the formation of the breaks seems to be insignificant.

7.4. Change of mass or a₀

Galaxies can increase their baryonic mass, either by forming stars in situ or by accreting other galaxies. Early-type galaxies, that are investigated in this paper, are also suspected to loose their baryonic mass through baryonic outflows (Fan et al. 2008; Damjanov et al. 2009; Ragone-Figueroa & Granato 2011). Let us approximate the galaxy by a point source that has a mass of M initially and ϵM after the mass change. According to the radial acceleration relation, the gravitational acceleration in the strong field regime will change by a factor of ϵ, while in the weak field regime by a factor of $\sqrt{ϵ}$. Using the argument with the conservation of angular momentum from Sect. 7.2, we found that a circular orbit with an initial radius of r₀ will change, after the mass change, to r₁ = r₀/ϵ in the strong field regime and to r₁ = r₀/ϵ^1/4 in the weak field regime. This suggests that the change of the baryonic mass could introduce only a different normalization of the GCS density profile inside and beyond the a₀ radius, but not the observed change of slope. On the other hand, GCs in GCSs move on trajectories with different eccentricities, which might possibly change the situation. In total, the change of mass does not seem to be responsible for the formation of the observed breaks, but simulations are necessary to exclude this option definitively.

It was proposed that in a MOND universe, the value of a₀ could vary with time (Milgrom 1983c). Milgrom (2015) showed that the variation would affect the orbits of objects in a similar fashion as the change of mass described above, but the orbits would be expanding. Again, the above arguments can be used against this mechanism to be responsible for the formation of the observed breaks of density of GCSs.

7.5. Other influences shaping the density profile of a GCS

The radial profiles of number density of GCSs can be also shaped by other processes than those listed above. Here we mention some examples, even if it is not currently clear why they would introduce GCS density breaks near the a₀ radii. They should be considered in more complex models of the formation of the density profiles of GCSs.

Early-type galaxies are known to grow their effective radius with time (Daddi et al. 2005; Trujillo et al. 2006; van Dokkum et al. 2009), even if they do not form new stars or do not increase their average stellar masses substantially. A promising explanation of the growth of their radii are repeated minor mergers (Naab et al. 2009). The premerger potential energies of the accreted satellites are transformed into kinetic energy of the stars and dark matter particles of the merger remnant. It is plausible that GCs would absorb some of the energy of the satellites too. This would make the distribution of the GCs more extended, as already explored for major mergers in Bekki & Forbes (2006).

The distribution of the GCs can further be influenced by the tidal destruction of the GCs. Many GCs of the Milky Way are known to have stellar streams (e.g., Ibata et al. 2021). It is possible that the mechanisms listed above will direct some GCs on radial orbits and once the GCs appear close to the center of the galaxy, they would be destroyed by tidal forces (Brockamp et al. 2014). This process would cause the a decrease of the density of the GCS near its center.

The density profiles of dark halos derived from observations under the assumption of Newtonian dynamics tend to show central cores (Kormendy & Freeman 2004; Oh et al. 2008; Donato et al. 2009; Salucci et al. 2012). The proposed explanations include the baryonic feedback (e.g., Governato et al. 2010; Di Cintio et al. 2014). Given that the feedback influences the motion and distribution of the dark matter particles, it will most probably influence the distribution of GCs in the same way. This would explain the observed flattening of the GC density profile toward the center of the galaxy. It is however unclear whether the core would have the a₀ radius. A mere presence of the dark matter core, formed by any mechanism, could affect the distribution of the GCs in the way described in Sect. 7.1.

8. Summary and conclusions

It was found in BSR19 that the number density profiles of GCSs follow broken power laws (Eq. (1)) whose break radii coincide with the a₀ radii of their host galaxies. An a₀ radius is defined as the radius at which the acceleration generated by the baryons of the galaxy equals the galactic acceleration constant a₀. It was shown in BSR19 that the a₀ radii coincide with the break radii better than other characteristic length scales of galaxies, such as the effective radii or the dark matter halo characteristic radii. The galaxy sample of BRS19 nevertheless spans only a relatively narrow range of baryonic masses of the galaxies, namely about one decade. They investigated the distribution of GCs only on the basis of spectroscopic catalogs that can suffer from geometrical incompleteness. They investigated only the total population of GCs. It was unclear, for example, whether the density profiles of the blue and red GC subpopulations follow the broken power-law as well and, if so, whether the break radii of the subpopulations are located at the a₀ radius. In the current contribution, we aimed to overcome these deficiencies. We analyzed two catalogs of photometric GC candidates in the Fornax cluster, one based on the ground-based Fornax Deep Survey data and the other on the ACS Fornax cluster data. The profiles of density of GCSs for the lowest-mass galaxies were derived from stacks of the GC candidates over many galaxies of a similar stellar mass. Additionally, we inspected a new spectroscopic catalog of GCs in the vicinity of the central galaxy of the cluster, NGC 1399. We were investigating only the GCSs of ETGs, since for LTGs, it is not possible to distinguish GCs from the numerous compact regions of star formation in the disks. The galaxy sample studied here spans logarithmic stellar masses from 8.0 to 11.4 M_⊙. The fitted parameters of the GCS density profiles are listed in Table 2.

Our observational findings can be summarized as follows:

1. We were able to detect breaks in the GCS profiles of virtually all galaxies (Appendix A). The exceptions are only the galaxies that have too few GCs. The breaks were found in the entire GC population as well as in the blue and red GC subpopulations (see the figures in Appendix A).

2. In the cases when the outer part of the broken power-law profile was observed in both the FDS and ACS catalogs, the outer slopes, b, generally agreed well, just as the break radii (Appendix A). The inner slope in the FDS data was usually shallower than in the ACS data. We attribute this difference to the greater difficulty of detecting faint GCs near galactic centers in the ground-based data because of the contamination by the light of the host galaxy (Sect. 3).

3. The break radii of the total GC population and the red and blue subpopulations are rather similar (Sect. 5.1). There is a marginal tendency for the blue GCs to have systematically higher break radii than the red GCs at the 1.5σ confidence level, namely by 0.3 ± 0.2 kpc, on average.

4. We calculated the a₀ radii in two ways: assuming the Newtonian gravity and the MONDian gravity, which provides a clear theoretical understanding of the existence of the acceleration scale a₀. We found that the break and a₀ radii agree typically within a factor of two (Sect. 5.3). The a₀ radii calculated from the MOND gravity are less offset from the break radii than the a₀ radii calculated from the Newtonian gravity.

5. The break radii also show significant correlations with the stellar masses of the galaxies, their effective radii, and Sérsic indices (Sect. 5.2). None of these correlations is however close to the one-to-one relation.

6. Gravitational fields of some galaxies are weaker than the constant a₀ in the whole extents of the galaxies (Sect. 4). Such galaxies thus do not have a₀ radii. They still show broken power-law density profiles for their GCSs.

7. The outer slope of the GCSs profiles, b, strongly correlates with the a₀ radii (both Newtonian and MOND, Sect. 5.2). The correlations of b with the a₀ radii are more significant than the correlation of b with the stellar mass of the host galaxy. This suggests that the mechanism that sets the profile of the GCS is causally linked rather with the spatial distribution of mass in the galaxy than with its total stellar mass.

8. The parameter b for the blue GCs is higher (i.e., the profile is less steep) than for the red GCs at the 2σ confidence level (Sect. 5.1).

9. We inspected in more detail the galaxy with the highest number of GCs, which is NGC 1399, the central galaxy of the Fornax cluster. We used the new catalog of spectroscopic GCs (Chaturvedi et al. 2022). We divided the GCs into groups according to a similar absolute value of the radial velocity with respect to the center of the galaxy. The shape of the profile shows systematic trends with the mean velocity of the selected group (Sect. 6.2, Fig. 6).

10. We divided the spectroscopic GCs around NGC 1399 in azimuthal sectors centered on the galaxy and derived the radial profiles of the density of GCs in each of the sectors. When the break radii are plotted in polar coordinates according to the angle of the middle line of the sector, they form an ellipse in the plane of the sky (Sect. 6.3, Fig. 8). The major axis of the ellipse points toward the neighboring galaxy NGC 1404. It is well known that these two galaxies are undergoing a tidal interaction. This demonstrates that break radii are influenced by galaxy interactions.

BSR19 proposed several explanations for the approximate match of the break and a₀ radii. We explored those and a few others in more detail (Sect. 7). We made use of simple models and observational arguments. None of them explains our findings completely satisfactory. More elaborate models and simulations are desirable. More data, ideally coming from a larger variety of environments and galaxy morphological types, could give us hints as to why the a₀ and break radii are similar.

Acknowledgments

We thank the referees for their comments that helped to improve the quality of the manuscript. M.B. acknowledges the support by the ESO SSDF grant 21/10. F.R. acknowledges support from the Knut and Alice Wallenberg Foundation, and from the University of Strasbourg Institute for Advanced Study (USIAS) within the French national program Investment for the Future (Excellence Initiative) IdEx-Unistra. S.S. acknowledges the financial support of the Ministry of Education, Science and Technological Development of the Republic of Serbia through the contract No. 451-03-68/2022-14/200002. The Digitized Sky Surveys were produced at the Space Telescope Science Institute under U.S. Government grant NAG W-2166. The images of these surveys are based on photographic data obtained using the Oschin Schmidt Telescope on Palomar Mountain and the UK Schmidt Telescope. The plates were processed into the present compressed digital form with the permission of these institutions. The National Geographic Society – Palomar Observatory Sky Atlas (POSS-I) was made by the California Institute of Technology with grants from the National Geographic Society. The Second Palomar Observatory Sky Survey (POSS-II) was made by the California Institute of Technology with funds from the National Science Foundation, the National Geographic Society, the Sloan Foundation, the Samuel Oschin Foundation, and the Eastman Kodak Corporation. The Oschin Schmidt Telescope is operated by the California Institute of Technology and Palomar Observatory. The UK Schmidt Telescope was operated by the Royal Observatory Edinburgh, with funding from the UK Science and Engineering Research Council (later the UK Particle Physics and Astronomy Research Council), until 1988 June, and thereafter by the Anglo-Australian Observatory. The blue plates of the southern Sky Atlas and its Equatorial Extension (together known as the SERC-J), as well as the Equatorial Red (ER), and the Second Epoch [red] Survey (SES) were all taken with the UK Schmidt. Supplemental funding for sky-survey work at the ST ScI is provided by the European Southern Observatory.

References

Appendix A: Detailed results of the GC surface density fitting

Fig. A.1. Fits of the surface density profile of the sources centered on ESO 358-006.

Fig. A.2. Fits of the surface density profile of the sources centered on ESO 358-050.

Fig. A.3. Fits of the surface density profile of the sources centered on NGC 1316.

Fig. A.4. Fits of the surface density profile of the sources centered on NGC 1336.

Fig. A.5. Fits of the surface density profile of the sources centered on NGC 1351.

Fig. A.6. Fits of the surface density profile of the sources centered on NGC 1373.

Fig. A.7. Fits of the surface density profile of the sources centered on NGC 1379.

Fig. A.8. Fits of the surface density profile of the sources centered on NGC 1380.

Fig. A.9. Fits of the surface density profile of the sources centered on NGC 1380B.

Fig. A.10. Fits of the surface density profile of the sources centered on NGC 1381.

Fig. A.11. Fits of the surface density profile of the sources centered on NGC 1387.

Fig. A.12. Fits of the surface density profile of the sources centered on NGC 1399.

Fig. A.13. Fits of the surface density profile of the sources centered on NGC 1404.

Fig. A.14. Fits of the surface density profile of the sources centered on NGC 1419.

Fig. A.15. Fits of the surface density profile of the sources centered on NGC 1427.

Fig. A.16. Fits of the surface density profile of the sources centered on NGC 1428.

Fig. A.17. Fits of the surface density stacked profile of sources centered on galaxies with stellar masses log10M*/M⊙ = 8 ± 0.25.

Fig. A.18. Fits of the surface density stacked profile of sources centered on galaxies with stellar masses log10M*/M⊙ = 8.5 ± 0.25.

Fig. A.19. Fits of the surface density stacked profile of sources centered on galaxies with stellar masses log10M*/M⊙ = 9 ± 0.25.

Appendix B: The influence of magnitude cut on the measurement of the break radius

Before fitting the GC density profiles, we sometimes had to remove from the analyzed sample of GC candidates the faint objects, such that we mitigate the problems with the contamination by the light of the host galaxy. In this appendix, we investigate whether imposing the magnitude cut has any impact on the break radius.

We used the galaxy NGC 1399 for this test. Since it hosts the highest number of GC candidates, it is possible to measure the break radius relatively precisely even if most of them are excluded from the fit.

The results are shown in Fig. B.1. The values of the break radius agree with each other within the uncertainty limits. The faintest GC candidates in the sample are around 26 mag.

Fig. B.1. Demonstration that the break radius does not depend on the applied magnitude cut. The data for the galaxy NGC 1399 were used.

Appendix C: The effect of ellipticity of the GCS on the extracted profile of its surface density

In our work, we used circular annuli to extract the radial profiles of surface density of GC candidates. We neglected the fact that the distribution of GCs around a galaxy can be flattened. It has been found that the ellipticity of GCs, particularly of the red ones, follows the ellipticity of the host galaxy (see Brodie & Strader 2006 for a review). For example the galaxy NGC 1380, one of the galaxies of our sample, has a rather high ellipticity of 0.5 (Kissler-Patig et al. 1997). We thus explore here how the profile of the surface denstity of GC candidates, extracted in our way, would change if the GC system would actually be elliptical.

The real data for the galaxy NGC 1380 are not suitable for this task. The break radius is located at the border of the ACS FOV. In the FDS data the break is located inside the region that was excluded because of the problems with the contamination by the light of the host galaxy, see Fig. A.8. We therefore decided to use artificial data for the test.

We constructed a circular and an elliptical GCS. The positions of 500 GCs were initially generated in three dimensions. The positions were drawn from a spherical distribution whose volume density followed a broken power law described by the parameters a = −1.7, b = −3.4, r_br = 1′. These are the typical parameters for the real galaxies (Table 4). For constructing the circular GCS, the Cartesian coordinates x_c and y_c of the modeled GCs were then directly considered the on-sky positions of the GCs. For the elliptical GCS, the on-sky coordinates of the GCs were constructed as $x_e=\sqrt{ϵ}{x}_c$ and $y_c/\sqrt{ϵ}$, where ϵ is the ellipticity of the GCS, chosen to be ϵ = 0.5. Finally, we added contaminants in the list of the modeled GCs candidates. The contaminants were assumed to have a uniform distribution on the sky with a mean density of 1 arcsec⁻².

Both the circular and elliptical GCSs were then analyzed in the same way as the observational data. Figure C.1 shows the comparison of the extracted profiles of surface density of the modeled GC candidates and fits by the broken power laws (Eq. 10). Both the extracted surface density profiles and the best-fit curves are very similar. This is confirmed by consistency of the values of the fitted parameters within their uncertainty limits, as shown in Table C.1. The last column of the table, listing the true values of parameters of the GCS density profile and the density of contaminants, confirms that our method is able to recover these parameters correctly.

Fig. C.1. Demonstration that the ellipticity of the distribution of the GCs in the GCS does not influence the extracted profile of surface density of GC candidates substantially. Artificial data. The flattened GCS has the on-sky ellipticity of 0.5.

Appendix D: Inhomogeneity of the sensitivity of the FDS survey

As we mentioned in Sect. 2.3, the catalog of FDS GC candidates shows spatial inhomogeities that form a tile-like pattern, see Fig. 1. The pattern seems to be arising because of a varying sensitivity of the survey both between the individual tiles of the mosaic and inside of the individual tiles. We were not able to remove the tile pattern even if excluding all but the brightest sources. Here we investigate whether the sensitivity variations could affect our measurements of the GCS density profiles. We focused on the cases of the galaxy with the richest and most extended GCS, that is NGC 1399, because it would be affected most by the large-scale sensitivity variations. This galaxy is located close to the center of one of the tiles.

To this end, we constructed radial profiles of the density of sources in four tiles close to NGC 1399, centered approximately on the centers of the respective tiles. In particular, we used the two tiles adjacent from the west to the tile containing the galaxy, and the two tiles adjacent from the east, along the line of a constant declination, see Fig. 1. These tiles do not contain any galaxies with substantial GCSs. The profiles were extracted up to the distance of 30′, which is the distance that was used for constructing the GCS profile of NGC 1399. The regions occupied by the GCSs of intermediate galaxies were masked in the way described in Sect. 3.1 before extracting the profiles. The radial bins were chosen such that each contains 800 sources. The measured profiles are shown in Fig. D.1 as the thin jiggles lines. The figure also shows by the thicker lines linear fits to the extracted profiles. The figure shows that the sensitivity at the outermost radius is different by around 10% from the central sensitivity. The sensitivity can both increase and decrease with the distance from the center of a given tile. The central surface densities of sources can differ by several tens of percent between the individual tiles.

Fig. D.1. Variation of the density of sources with distance from the centers of four tiles of the FDS mosaic. These are the two tiles to the west of the tile containing NGC 1399 and the two tiles to the east.

We further investigated quantitatively what is the effect of the sensitivity variations on our estimates of the parameters of the broken power-law profiles of the GCSs. We assumed that the sensitivity variations visible in Fig. 1 affect GCs and contaminating sources in the same way. As a first approximation, we assumed that the sensitivity is a linear function of the distance from the center of the tile. We thus took the measured profile of the source density of NGC 1399 Σ(R) and transformed it to:

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_2(R)=\mathrm{\Sigma }(R)(1+vR/R_{\max }),\end{array}$$(D.1)

where R_max = 30′ is the radius of the last bin measured for NGC 1399. The parameter v quantifies the magnitude of the variation of sensitivity. It is the ratio of the sensitivities in the center of the tile (which is the same as the center of NGC 1399) and at R_max. The profile Σ₂ was then fitted by the projected broken power law, in the same way we did with the real data.

The results are presented in Table D.1. The table reveals that the position of the break radius is not very sensitive to this type of spatial variation of sensitivity, even for very strong variations of sensitivity. On the other hand, a realistic variation of the sensitivity of plus or minus 10% already has a substantial effect on the measured outer slope b.

This sheds doubt on the measured value of the slopes b for the two galaxies with the most extended GCSs, that is NGC 1399 and NGC 1316. This is why we marked in Table 2 the values of these parameters for these two galaxies as suspicious. The GCSs of the other galaxies in our sample are much smaller (Fig. 1 and the values of break radii in Table 2). The profiles of these galaxies would not be affected much.

Appendix E: Linear fits of the relation between the break radii and a₀ radii

Fig. E.1. Fit of the relation of the observed break radius and the MOND a0 radius by a law in the form of log10rbr = C1log10ra0 + C0. The black crosses represent the uncertainty intervals of the data. The full red line show the fit and the dotted red lines the fitted intrinsic scatter of the relation. The green dashed line indicates the one-to-one relation. The deviation of the best-fit curve from the one-to-one relation is driven by a few data points.

We made linear fits of the relation between the a₀ and break radii, and recover the intrinsic scatter of the correlation for several datasets. More specifically, we assumed that at a given r_br, log₁₀r_a0 follows a Gaussian distribution with a mean of C₁log₁₀r_br + C₀ and a standard deviation of σ_int. This implies that at a given r_br, r_a0 follows the lognormal distribution

$$\begin{array}{cc}\hfill G(r_{\mathrm{br}},r_{a_0})& =\frac{1}{\sqrt{2\pi }ln(10){\sigma }_{\mathrm{int}}{r}_{\mathrm{br}}}\hfill \\ \hfill & exp[-\frac{{\{{log}_{10}[r_{a_0}]-(C_1{log}_{10}[r_{\mathrm{br}}]+C_0)\}}^2}{2{\sigma }_{\mathrm{int}}^2}].\hfill \end{array}$$(E.1)

We had to take into account the uncertainty in the estimates of the break and a₀ radii. For the i-th galaxy, we assumed for the pair of variables (r_br, r_a0) the distribution function:

$$\begin{array}{c}\hfill g_{\mathrm{mes},i}(r_{\mathrm{br}},r_{a_0})=g(r_{\mathrm{br}},\overline{r_{\mathrm{br},i}},{\sigma }_{r_{\mathrm{br},i}}^{+},{\sigma }_{r_{\mathrm{br},i}}^{-})g(r_{a_0},\overline{r_{\mathrm{eq},i}},{\sigma }_{r_{\mathrm{eq},i}}^{+},{\sigma }_{r_{\mathrm{eq},i}}^{-}),\end{array}$$(E.2)

where g was defined in Eq. 16. The estimated slope, intercept and intrinsic scatter of the correlation between r_br and r_a0 were found by maximizing the likelihood function:

$$\begin{array}{c}\hfill \mathcal{L}(C_1,C_0,{\sigma }_{\mathrm{int}})={\displaystyle \prod _i}[G(r_{\mathrm{br}},r_{a_0})\ast g_{\mathrm{mes},i}(r_{\mathrm{br}},r_{a_0})](\overline{r_{\mathrm{br},i}},\overline{r_{\mathrm{eq},i}}),\end{array}$$(E.3)

where the asterisk denotes convolution over both variables. The intrinsic scatter was restricted to be at least 0.01 dex, which is much less than the measurement uncertainties of r_br and r_a0, in order to avoid numerical problems.

This procedure was applied to several datasets. First, to the Fornax Cluster data presented in this paper, for the total, red and blue GC populations, then it was applied to the union of the data from this paper and from BSR19 for the total GC populations. For NGC 1399, that appears in both of the samples of this paper and BSR19, we used only the datapoint from this paper. Finally, we divided the union dataset in two parts according to whether the GCS has a break radius greater or lower than 20 kpc, since, as we noted in Sect. 5.3, the trend might be different in these two regimes. All of this was done for the Newtonian and MOND a₀ radii. The results are stated in Table E.1. The values are in many cases not consistent with the one-to-one relation between the break and a₀ radii.

We demonstrate the origin of this disagreement on the example of the blue GC population and the MOND a₀ radii, for which the fit indicated a zero slope. The fitted data and the best-fit curve are shown in Fig. E.1. The full red line indicates the best-fit linear fit of the relation between log₁₀r_a0 and log₁₀r_br and the dotted lines the fitted intrinsic scatter. The green dashed line indicates the one-to-one relation. A likelihood ratio test indeed indicates that the best-fit model is significantly better than the green line. It seems that the deviation of the fitted curve from the one-to-one relation is caused by a few datapoints. At this moment, it is not clear whether these points deviate because they are outliers (either because of observational errors or because the galaxies experienced an unusual event) or the true correlation between the break and a₀ radii is not a one-to-one correlation exactly. This will become clearer when more data points become available. For the moment, we were able to prove the match between the break and a₀ radii only in the sense that they agree within a factor of a about two, as demonstrated in Sect. 5.3.

All Tables

All Figures

	Fig. 1. Positions of GC candidates in the FDS catalog, after applying selection criteria. The coordinate system is centered on NGC 1399. North is up, and west is to the right.
In the text

	Fig. 2. Statistically significant correlations of the outer slope of the density profile of the GCSs b with the properties of their host galaxies. Top: the MOND a0 radius. Bottom: the Newtonian a0 radius. The open symbol indicates a suspicious measurements. The red numbers in the corners indicate the p-value of the Spearman’s correlations.
In the text

	Fig. 3. Statistically significant correlations of the break radii of GCSs with the properties of their host galaxies. From top to bottom: stellar mass, effective radius, Sérsic index, MOND a0 radius, Newtonian a0 radius. The open symbols indicate suspicious measurements. The red numbers in the corners indicate the p-value of the Spearman’s correlations.
In the text

Fig. 4. Relation of the a0 radius versus the break radius of the GCS for each galaxy. The upper row corresponds to the MOND a0 radius, the bottom row to the Newtonian a0 radius. The three columns in the horizontal direction correspond to the whole, red, and blue GC (sub)population, respectively. The squares in the first column mark the galaxy sample inspected in this paper, the circles mark those from BSR19. The empty squares denote the suspicious measurements. The diagonal dashed black lines represent the one-to-one relation. The vertical gray dashed lines mark the galaxies that do not have a0 radii.

In the text

Fig. 5. Demonstration of the approximate equality of the break radii of the GCSs and the a0 radii of their host galaxies. First row: the ratio of the break radius and the MOND a0 radius as a function of the break radius. The vertical axis is in the base-2 logarithmic scale, while the horizontal axis in the base-10 logarithmic scale. The horizontal line indicates the ratio of one (the break and a0 radii are equal). Second row: the same as the first row but for the a0 radii calculated for the Newtonian gravity. Columns: from left to right, the rows show the data derived for all GCs and for the red and blue subpopulations of GCs.

In the text

	Fig. 6. Radial profiles of surface density of GCs of NGC 1399 in several bins of radial velocity with respect to the center of the galaxy.
In the text

	Fig. 7. Same as in Fig. 6, but for modeled data. The models are Jeans spherical models with the indicated types of anisotropy parameters. The models were based on the gravitational potential of NGC 1399 derived in Bílek et al. (2019b) and the parameters of the GCS used in that work. The vertical dashed line at 30′ is an eye-guide to facilitate the reading of the variations of the break radius.
In the text

	Fig. 8. Break radii of the GCS of NGC 1399 in different sectors centered on the galaxy. The radial coordinate is in arcminutes. The red point indicates the position of NGC 1404.
In the text

	Fig. 9. Impact of the increasing of the external field on the trajectories of GCs that were originally on circular orbits of different radii, assuming the MOND gravity. The horizontal line indicates the MOND a0 radius.
In the text

	Fig. 10. Effect of dynamical friction on the apocentric radii of GCs that initially move in apocenter with a velocity that is equal to 0.1 of the local circular velocity.
In the text

	Fig. A.1. Fits of the surface density profile of the sources centered on ESO 358-006.
In the text

	Fig. A.2. Fits of the surface density profile of the sources centered on ESO 358-050.
In the text

	Fig. A.3. Fits of the surface density profile of the sources centered on NGC 1316.
In the text

	Fig. A.4. Fits of the surface density profile of the sources centered on NGC 1336.
In the text

	Fig. A.5. Fits of the surface density profile of the sources centered on NGC 1351.
In the text

	Fig. A.6. Fits of the surface density profile of the sources centered on NGC 1373.
In the text

	Fig. A.7. Fits of the surface density profile of the sources centered on NGC 1379.
In the text

	Fig. A.8. Fits of the surface density profile of the sources centered on NGC 1380.
In the text

	Fig. A.9. Fits of the surface density profile of the sources centered on NGC 1380B.
In the text

	Fig. A.10. Fits of the surface density profile of the sources centered on NGC 1381.
In the text

	Fig. A.11. Fits of the surface density profile of the sources centered on NGC 1387.
In the text

	Fig. A.12. Fits of the surface density profile of the sources centered on NGC 1399.
In the text

	Fig. A.13. Fits of the surface density profile of the sources centered on NGC 1404.
In the text

	Fig. A.14. Fits of the surface density profile of the sources centered on NGC 1419.
In the text

	Fig. A.15. Fits of the surface density profile of the sources centered on NGC 1427.
In the text

	Fig. A.16. Fits of the surface density profile of the sources centered on NGC 1428.
In the text

	Fig. A.17. Fits of the surface density stacked profile of sources centered on galaxies with stellar masses log10M*/M⊙ = 8 ± 0.25.
In the text

	Fig. A.18. Fits of the surface density stacked profile of sources centered on galaxies with stellar masses log10M*/M⊙ = 8.5 ± 0.25.
In the text

	Fig. A.19. Fits of the surface density stacked profile of sources centered on galaxies with stellar masses log10M*/M⊙ = 9 ± 0.25.
In the text

	Fig. B.1. Demonstration that the break radius does not depend on the applied magnitude cut. The data for the galaxy NGC 1399 were used.
In the text

	Fig. C.1. Demonstration that the ellipticity of the distribution of the GCs in the GCS does not influence the extracted profile of surface density of GC candidates substantially. Artificial data. The flattened GCS has the on-sky ellipticity of 0.5.
In the text

	Fig. D.1. Variation of the density of sources with distance from the centers of four tiles of the FDS mosaic. These are the two tiles to the west of the tile containing NGC 1399 and the two tiles to the east.
In the text

Fig. E.1. Fit of the relation of the observed break radius and the MOND a0 radius by a law in the form of log10rbr = C1log10ra0 + C0. The black crosses represent the uncertainty intervals of the data. The full red line show the fit and the dotted red lines the fitted intrinsic scatter of the relation. The green dashed line indicates the one-to-one relation. The deviation of the best-fit curve from the one-to-one relation is driven by a few data points.

In the text