# The survival of globular clusters in a cuspy Fornax

Shi Shao1, Marius Cautun2,1, Carlos S. Frenk1, Marta Reina-Campos3,4,5, Alis
J. Deason1, Robert A. Crain6, J. M. Diederik Kruijssen3 and Joel Pfeffer6
1 Institute for Computational Cosmology, Department of Physics, Durham
University, South Road Durham DH1 3LE, UK
2 Leiden Observatory, Leiden University, PO Box 9513, NL-2300 RA Leiden, the
Netherlands
3 Astronomisches Rechen-Institut, Zentrum für Astronomie der Universität
Heidelberg, Mönchhofstrasse 12-14, D-69120 Heidelberg, Germany
4 Department of Physics & Astronomy, McMaster University, 1280 Main Street
West, Hamilton, ON, L8S 4M1, Canada
5 Canadian Institute for Theoretical Astrophysics (CITA), University of
Toronto, 60 St George St, Toronto, ON M5S 3H8, Canada
6 Astrophysics Research Institute, Liverpool John Moores University, 146
Brownlow Hill, Liverpool L3 5RF, UK
E-mail<EMAIL_ADDRESS>
###### Abstract

It has long been argued that the radial distribution of globular clusters
(GCs) in the Fornax dwarf galaxy requires its dark matter halo to have a core
of size $\sim 1$ kpc. We revisit this argument by investigating analogues of
Fornax formed in E-MOSAICS, a cosmological hydrodynamical simulation that
self-consistently follows the formation and evolution of GCs in the EAGLE
galaxy formation model. In EAGLE, Fornax-mass haloes are cuspy and well
described by the Navarro-Frenk-White profile. We post-process the E-MOSAICS to
account for GC orbital decay by dynamical friction, which is not included in
the original model. Dynamical friction causes $33~{}$per cent of GCs with
masses $M_{\rm GC}\geq 4\times 10^{4}{~{}\rm M_{\odot}}$ to sink to the centre
of their host where they are tidally disrupted. Fornax has a total of five
GCs, an exceptionally large number compared to other galaxies of similar
stellar mass. In the simulations, we find that only 3 per cent of the Fornax
analogues have five or more GCs, while 30 per cent have only one and 35 per
cent have none. We find that GC systems in satellites are more centrally
concentrated than in field dwarfs, and that those formed in situ (45 per cent)
are more concentrated than those that were accreted. The survival probability
of a GC increases rapidly with the radial distance at which it formed ($r_{\rm
init}$): it is 37 per cent for GCs with $r_{\rm init}\leq 1~{}\rm kpc$ and 92
per cent for GCs with $r_{\rm init}\geq 1~{}\rm kpc$. The present-day radial
distribution of GCs in E-MOSAICS turns out to be indistinguishable from that
in Fornax, demonstrating that, contrary to claims in the literature, the
presence of five GCs in the central kiloparsec of Fornax does not exclude a
cuspy DM halo.

###### keywords:

globular clusters: general – galaxies: dwarfs – methods: numerical

††pubyear: 2020††pagerange: The survival of globular clusters in a cuspy
Fornax–C

##  1 Introduction

One of the fundamental predictions of the standard cosmological model
($\Lambda$CDM) is that the density profiles of dark matter (DM) haloes have
cuspy profiles of the Navarro-Frenk-White (NFW) form, where the inner DM
density profile follows $\rho\propto r^{-1}$ (Navarro et al., 1996b, 1997). In
contrast, some studies of stellar motions and rotation curves of faint
galaxies often appear to indicate constant-density profiles at the centre,
$\rho\propto r^{0}$ (e.g. Moore, 1994; Flores & Primack, 1994; Battaglia et
al., 2008; Walker & Peñarrubia, 2011). The apparent discrepancy between
theoretical predictions and observations has become known as the core-cusp
problem and several hypotheses have been proposed to explain it. These include
modifications of the assumed nature of DM (e.g. Spergel & Steinhardt, 2000;
Rocha et al., 2013; Shao et al., 2013; Kaplinghat et al., 2016; Schneider et
al., 2017), galaxy formation processes that alter the structure of the halo,
(e.g. Navarro et al., 1996a; Read & Gilmore, 2005; Pontzen & Governato, 2012;
Di Cintio et al., 2014; Benítez-Llambay et al., 2019), and systematic biases
in the interpretation of the observational data (e.g. Oman et al., 2019).

Dwarf galaxies are important cosmological probes since their inner structure
is particularly sensitive to the nature of DM; this has served to stimulate
numerous studies of the rotation curve and kinematics of dwarfs. The
kinematics of stars are tracers of the potential and multiple analysis methods
have been employed to infer the mass profile of the DM halo in which they
orbit, such as Jeans, and more complex, modelling (e.g. Gilmore et al., 2007;
Walker et al., 2009; Amorisco et al., 2013; Pascale et al., 2018; Read et al.,
2019). However, due to the relatively low number of stars with precise radial
velocities measurements and the lack of accurate proper motion measurements,
the current data are unable to distinguish robustly between halo profiles with
central cores or cusps in dwarfs (e.g. Strigari et al., 2010, 2017; Genina et
al., 2018, 2019a). Dwarfs contains not only stars but many of them also host
globular clusters (GCs) (e.g. Mackey & Gilmore, 2003; Law & Majewski, 2010;
McConnachie, 2012; Forbes et al., 2018), and these can also be used to
constrain the profile of DM haloes.

Even though the scatter in the number of GCs associated with dwarf galaxies is
large (e.g. Forbes et al., 2018), their presence in dwarf spheroidals like
Fornax has long been considered a useful probe of the DM distribution in these
galaxies. In particular, the GC system in Fornax has received considerable
attention. It consists of five GCs that are found at an average projected
distance of $\sim 1.1~{}\rm kpc$ from the centre (Mackey & Gilmore, 2003).
They are typically very old, with a lookback formation time, $t_{\rm
age}\approx 12$ Gyrs (Fornax 4 is an exception with $t_{\rm age}\approx 10$
Gyrs; see Table 1). The interest in the Fornax GCs stems from comparing their
ages to the orbital decay times (the time required for the GC to sink to the
centre) due to dynamical friction. If the GCs formed at their present-day
positions, then the orbital decay time in a cuspy DM halo would be shorter
than their present age (e.g. Hernandez & Gilmore, 1998; Goerdt et al., 2006);
this discrepancy has been called the “GC timing problem”. One simple solution
is to increase the decay timescales by assuming larger initial radii for the
GCs rather than assuming that they formed at their present-day positions (e.g.
Angus & Diaferio, 2009; Boldrini et al., 2019; Meadows et al., 2020). However,
the initial positions where the GCs were born are difficult to constrain due
to the long timescales involved.

An alternative solution is to assume that the inner density profile of the
Fornax halo has a kiloparsec-sized core. Goerdt et al. (2006) used idealized
N-body simulations to investigate the orbits of GCs under different
assumptions for the DM profile. They argued that when dynamical friction is
taken into account, the GCs would sink in a halo with a flat inner density
profile and stall at some radius for a few gigayears. In contrast, in a cuspy
halo, the GCs would gradually sink all the way to the centre and be disrupted
quickly by the central tidal field or perhaps end up as a nuclear star cluster
(e.g. Tremaine et al., 1975). Cole et al. (2012) also reported that if the
core is relatively large, GCs could be pushed out by the so-called “dynamical
buoyancy” mechanism. Based on the argument that GCs would stall in a halo with
a central core, several studies have estimated the core radius of Fornax from
the present-day positions of its GCs (e.g. Read et al., 2006; Inoue, 2009;
Cole et al., 2012; Petts et al., 2015; Kaur & Sridhar, 2018; Boldrini et al.,
2019).

Another possibility has been proposed by Boldrini et al. (2020) who have shown
that if the GCs are surrounded by their own DM minihaloes, then there is no
“GC timing problem” even in cuspy haloes. Finally, it is worth mentioning that
tidal streams produced by the disruption of GCs stripped from dwarf galaxies
can perhaps be used to distinguish between a central cusp or a core in the
halo of the dwarf from which the GC was stripped (Malhan et al., 2020).

In this paper, we investigate if the present-day number and radial
distribution of the GCs in the Fornax dSph are consistent with $\Lambda$CDM
predictions in a cuspy DM halo. For this, we use one of the E-MOSAICS
simulations that self-consistently model galaxy and GC formation and evolution
in a fully cosmological context (Pfeffer et al., 2018; Kruijssen et al.,
2019a). The presence of GCs does not affect the DM halo profile and so the
E-MOSAICS dwarfs have the same DM profiles as in the EAGLE simulation, namely
they are cuspy and well fitted by the NFW parametrisation (e.g. Schaller et
al., 2015; Benítez-Llambay et al., 2019; Bose et al., 2019).

We proceed by selecting a sample of Fornax analogues in E-MOSAICS, that is
galaxies with similar stellar mass to the Fornax dSph, and identify the GCs
associated to each analogue. We cannot directly use the number and positions
of those GCs since E-MOSAICS does not include an accurate calculation of GC
orbital decay due to dynamical friction. Instead we account for this effect in
post-processing. Firstly we trace back the GCs to either their formation time
(for in-situ GCs) or their accretion time (for accreted GCs). Then, starting
from these initial positions and velocities, we follow the orbital decay of
the GCs by fitting an NFW profile to the host halo and analytically
integrating the orbits subject to dynamical friction down to the present day.
We account for the tidal disruption of the GCs that sink to the centre of the
host galaxy. Finally, we compare the resulting present-day population of GCs
in the simulation against the distribution of GCs in the Fornax dSph.

This paper is organised as follows. In Section 2 we describe the E-MOSAICS
simulations and, in Section 3, we introduce our Fornax analogue sample and
examine the properties of their DM haloes and GCs. In Section 4 we describe
our method for modelling dynamical friction and tidal disruption of the GCs.
In Section 5 we present our results on the number and radial distribution of
GCs in Fornax-mass dwarfs. We conclude with a discussion of our results in
Section 6 and a short summary of our main findings in Section 7.

Table 1: Selected properties of the Fornax GCs. The columns give: GC designation, projected distance from the centre of the Fornax dSph, stellar mass, age and metallicity. The projected distances are taken from Mackey & Gilmore (2003) and updated for the current distance of the Fornax dSph of $147~{}\rm kpc$ (McConnachie, 2012). The remaining properties are taken from de Boer & Fraser (2016). Name | $r_{p}\;[\rm{kpc}]$ | $M_{\rm GC}\;[10^{5}\;M_{\odot}]$ | Age [Gyrs] | [Fe/H]
---|---|---|---|---
Fornax 1 | 1.72 | 0.42$\pm$0.10 | 12.1$\pm$0.8 | -2.5$\pm$0.3
Fornax 2 | 1.13 | 1.54$\pm$0.28 | 12.2$\pm$1.0 | -2.5$\pm$0.3
Fornax 3 | 0.46 | 4.98$\pm$0.84 | 12.3$\pm$1.4 | -2.5$\pm$0.2
Fornax 4 | 0.26 | 0.76$\pm$0.15 | 10.2$\pm$1.2 | -1.2$\pm$0.2
Fornax 5 | 1.54 | 1.86$\pm$0.24 | 11.5$\pm$1.5 | -1.7$\pm$0.3
Fornax 6† | 0.030 | $\sim 0.29$ | – | –
† This is a recently discovered faint and highly elliptical GC, which is
possibly undergoing tidal disruption (Wang et al., 2019). We estimated its
mass by assuming that it has the same mass-to-light ratio as Fornax 1, which
is the Fornax GC with the closest V-band magnitude to Fornax 6.

##  2 Simulations

The E-MOSAICS (MOdelling Star cluster population Assembly In Cosmological
Simulations within EAGLE) suite of cosmological hydrodynamical simulations is
an extension of the ‘Evolution and Assembly of GaLaxies and their
Environments’ (EAGLE) simulations (Crain et al., 2015; Schaye et al., 2015)
that includes a subgrid model of stellar cluster formation, evolution, and
disruption (Kruijssen et al., 2011; Pfeffer et al., 2018). For a detailed
description of the physical ingredients of the model, we refer the reader to
(Pfeffer et al., 2018; Kruijssen et al., 2019a).

Currently, the E-MOSAICS project consists of two groups of simulations. The
first is a suite of 25 cosmological “zoom” simulations of MW-mass haloes
(Pfeffer et al., 2018; Kruijssen et al., 2019a). These contain only a small
number of Fornax analogues and are not used in this analysis. The second
group, which is the one we use, is a simulation of a full cosmological volume
in a periodic cube of side-length 34.4 (comoving) Mpc (Crain et al. in prep.).
The volume of the simulation is 2.6 times larger than the EAGLE “high-
resolution” simulation which has the same resolution (labelled
Recal-L025N0752). The simulation follows the evolution of $1034^{3}$ DM
particles and an initially equal number of gas particles. The DM particle mass
is $1.2\times 10^{6}{~{}\rm M_{\odot}}$, and the initial gas particle mass is
$2.3\times 10^{5}{~{}\rm M_{\odot}}$. E-MOSAICS assumes the Planck-1 cosmology
(Planck Collaboration XVI, 2014) with cosmological parameters: $\Omega_{\rm
m}=0.307,\Omega_{\rm
b}=0.04825,\Omega_{\Lambda}=0.693,h=0.6777,\sigma_{8}=0.8288$ and $n_{\rm
s}=0.9611$, which are those used by the EAGLE project.

As in the EAGLE project, the simulation we analyze here was performed with a
modified version of the Gadget code (Springel, 2005), which includes state-of-
the-art smoothed particle hydrodynamics (Dalla Vecchia & Schaye, 2012;
Hopkins, 2013; Schaller et al., 2015) and subgrid models, such as element-by-
element gas cooling, star formation, metal production, stellar winds, and
stellar and black hole feedback (Wiersma et al., 2009; Springel et al., 2005;
Booth & Schaye, 2009; Schaye et al., 2015). The parameters were calibrated so
as to reproduce three present-day observables: the stellar mass function, the
galaxy size–mass relation and the normalization of the relation between the
masses of supermassive black holes and the stellar mass of their host
galaxies. For a more detailed description we refer the reader to Schaye et al.
(2015).

The semi-analytic model for MOSAICS was coupled with EAGLE to track the
formation and evolution of star clusters. The model is calculated on-the-fly
within EAGLE since the time resolution required to resolve the rapid change of
the tidal field experienced by the GCs ($<1$ Myr), which drives most of their
disruption, is much finer than the time interval between simulation snapshots
(roughly $70$ Myrs). The formation of stellar particles in the simulation
triggers the formation of a subgrid population of star clusters which inherit
properties of their host stellar particles such as position, velocity, age,
and metallicity.

The number of star clusters and their masses at formation are determined by
two parameters, the cluster formation efficiency (CFE, Bastian 2008) and an
upper truncation mass scale in the Schechter (1976) initial cluster mass
function ($M_{c,\star}$), with a power-law index of -2 at low masses. These
parameters are described in terms of the local natal properties of the GC,
such as local ambient gas properties (pressure, density and mass) and stellar
velocity dispersion. Environments with higher gas pressure lead to the
formation of more star clusters which can also be more massive (Kruijssen,
2012; Reina-Campos & Kruijssen, 2017).

Four GC models are applied to the simulation in parallel. The fiducial model
allows the CFE and $M_{c,\star}$ to vary as a function of properties of the
local environment in which the stars are formed (Kruijssen, 2012; Reina-Campos
& Kruijssen, 2017). The other three keep either, or both of the CFE and
$M_{c,\star}$ fixed. For more details of the cluster formation models, we
refer the reader to Pfeffer et al. (2018) and Reina-Campos et al. (2019).
After the clusters form, they lose mass due to stellar evolution (according to
the fractional mass-loss of the parent stellar particle) and by dynamical
effects such as two-body relaxation and tidal shocks which are based on the
strength and change of the local tidal field respectively (Kruijssen et al.,
2011).

The E-MOSAICS simulations have been able to reproduce a broad range of
observational properties such as the deficit of massive metal-poor GCs (i.e.
the “blue tilt”) in galaxies across a wide range of environments (Usher et
al., 2018) and the diversity of age-metallicity relations of GCs in different
galaxies (Kruijssen et al., 2019a). The $M_{\rm GC}/M_{\rm halo}$ and $M_{\rm
GC}/M_{\star}$ relations predicted by the simulation are also in good
agreement with observations (Bastian et al., 2020). Additionally, the large
sample of MW-mass galaxies allows the assembly history of our MW and the GC
formation history to be probed by reference to the observed present-day GC
populations (Kruijssen et al., 2019b; Kruijssen et al., 2020; Reina-Campos et
al., 2019; Pfeffer et al., 2020; Trujillo-Gomez et al., 2020). The simulations
also have strong implications for the origin of the stellar bulge and stellar
halo whose masses are made up, in part, of GC remnants (Hughes et al., 2020;
Reina-Campos et al., 2020). In addition, the simulations can be used to make
predictions for the properties of the recently observed GCs in the M31 system
(Hughes et al., 2019).

The halo and galaxy catalogues in E-MOSAICS have been built using the tools
described by Schaye et al. (2015). Haloes are initially identified by the
friends-of-friends (FOF) algorithm (Davis et al., 1985) with a linking length
$0.2$ times the mean interparticle separation. The resulting FOF groups are
further split into gravitationally bound substructures using the subfind code
(Springel et al., 2001; Dolag et al., 2009), applied to the total matter
distribution (dark matter, gas and stars) associated with each FOF group. The
central subhalo is defined as the subhalo that contains the most bound
particle, while the remaining subhaloes are classified as satellites. The
stellar distribution associated with the main subhalo is identified as the
central galaxy. The central haloes are characterized by the mass, $M_{200}$,
and radius, $R_{200}$, that define an enclosed spherical overdensity of $200$
times the critical density. The position of each galaxy, for both centrals and
satellites, is given by their most bound particle.

Fig. 1 shows the relation between the stellar mass of central galaxies and the
mass of their host haloes in E-MOSAICS. The satellite galaxies are not shown
since the total mass associated with the subhalo changes as it orbits in the
main halo. The scatter in the stellar-to-halo mass relation for Fornax-mass
dwarfs is larger than for more massive galaxies, but is significantly smaller
than for lower mass dwarfs (Sawala et al., 2015).

Figure 1:  The relation between stellar mass, $M_{\star}$, and total halo
mass, $M_{200}$, for central galaxies in the simulation. The colours indicate
the number of galaxies in each halo and stellar mass bin (see legend). The
grey shaded region shows galaxies with stellar masses in the range $2-8\times
10^{7}{~{}\rm M_{\odot}}$, which corresponds to our sample of field Fornax-
mass dwarfs (we also select Fornax-mass satellites, which are not shown in
this diagram).

##  3 Sample selection and Methods

Figure 2:  The properties of the DM haloes that host a Fornax-mass dwarf
galaxy. The plot shows the PDF of the halo mass (left panel), concentration
(centre panel), and scale radius (right panel). The results are for present-
day field galaxies, since the same properties for satellite galaxies are not
directly comparable (see text for discussion on this). However, we do show the
distribution of scale radii for satellites as the solid red line in the right
panel. The vertical arrows indicate the median of each PDF.

###  3.1 The Fornax analogue sample selection and their DM halo properties

We select Fornax analogues, to which we also refer as Fornax-mass dwarfs, in
the simulation by requiring that they have a stellar mass in the range,
$M_{\star}\in[2,8]\times 10^{7}{~{}\rm M_{\odot}}$. This results in 1154
objects, of which 599 are field galaxies and 555 are satellites. The mass
range used for the selection corresponds to a factor of two variation around
the Fornax dSph stellar mass, which we take to be $4\times 10^{7}{~{}\rm
M_{\odot}}$ (de Boer et al., 2012). The mass range is relatively wide since we
need a large sample of Fornax-mass dwarfs for good statistics and a factor of
2 is representative of the uncertainty in Fornax’s mass estimates. A typical
Fornax-mass dwarf in E-MOSAICS is resolved with ${\sim}2\times 10^{4}$ DM
particles, and hundreds of star particles, which allows for a robust
determination of its present day properties as well as its formation history.

We study the general properties, such as halo mass, concentration, and scale
radius, $R_{s}$, of the Fornax-mass dwarf sample in Fig. 2. The left panel
shows the distribution of the host halo mass, $M_{\rm 200}$, for the field
Fornax-mass dwarfs. The satellites are not shown since they are tidally
stripped by the hosts after their infall. Nevertheless, for satellite dwarfs,
we refer to the $M_{\rm 200}$ of their main progenitors, before they became a
satellite of a more massive host, as their host halo mass. The distribution
for the field sample peaks at a value of $2.8^{+1.1}_{-1.0}\times
10^{10}{~{}\rm M_{\odot}}$ (68 per cent confidence limit), with sharp drop-
offs on both sides; this is in agreement with abundance matching predictions
and results from other hydrodynamical simulations (e.g. Moster et al., 2013;
Sawala et al., 2015; Oñorbe et al., 2015). The width of the distribution
largely reflects the fact that haloes of a given mass can have a range of
different concentrations; higher concentration haloes, which typically form
early, have more time to form stars and experience less efficient feedback.

The concentration and scale radius of the host haloes are calculated by
fitting the spherical DM density within a radial distance of $50~{}\rm kpc$
from the centre (roughly the median $R_{\rm 200}$ of the field Fornax-mass
haloes) to an NFW profile (Navarro et al., 1996b, 1997):

$\rho(r)=\dfrac{\rho_{0}}{\dfrac{r}{r_{\rm s}}\left(1+\dfrac{r}{r_{\rm
s}}\right)^{2}}\;,$ (1)

where $\rho_{0}$ and the scale radius, $r_{\rm s}$, are parameters (see
appendix A). For the satellites, we fit the DM profile of the main progenitor
just before infall. The concentration is given by the ratio $c=r_{\rm
s}/R_{200}$. The distribution of concentration values is shown in the middle
panel of Fig. 2; it has a median value of $c=10.5$ for the field Fornax-mass
dwarfs. The result is in good agreement with previous high-resolution
cosmological simulations (e.g. Hellwing et al., 2016).

The right panel of Fig. 2 shows the distribution of the scale radius, $r_{\rm
s}$, for both the field and satellite samples. The field sample has a median
value of $r_{\rm s}=6.2~{}\rm kpc$, slightly larger than that of the satellite
sample, $r_{\rm s}=5.9~{}\rm kpc$. The two samples have similar $z=0$ stellar
mass by construction; thus the difference in $r_{\rm s}$ between the two
samples is due to satellites having experienced tidal stripping and also,
potentially, to small differences in the assembly history of satellite versus
central galaxies.

###  3.2 GCs sample selection and their formation time

Here we compare the E-MOSAICS predictions with the five most massive GCs in
the Fornax dSph. These have been known for a long time and have been
thoroughly studied (e.g. de Boer & Fraser, 2016); we give a few selected
properties of the Fornax GCs in Table 1. More recently, Wang et al. (2019)
discovered a sixth GC in Fornax. It has a lower mass than the other five and
is possibly undergoing tidal disruption. For simplicity, and better to compare
with previous works, we limit our comparison to the five well-known Fornax
GCs.

To identify analogues of the Fornax GC population, we proceed by selecting all
GCs that at $z=0$ are associated with a Fornax analogue and have a stellar
mass above $4\times 10^{4}{~{}\rm M_{\odot}}$. The GC stellar mass threshold
corresponds to the mass of the Fornax-1 GC, the least massive of the five
Fornax GCs we study here. Our selection results in 2133 GCs, out of which only
1439 survive to the present day (the others were tidally disrupted after
dynamical friction dragged them to the centre – we describe this in detail in
Section 4.1).

GCs that are associated with their $z=0$ host may have formed inside another
galaxy that was subsequently accreted into the present-day host. Such GCs can
have different properties (e.g. age, metallicity, radial distribution) from
those born in situ. We study this by splitting our sample into in situ and
accreted GCs according to whether they were formed in the main (in situ) or in
a sub- (accreted) branch of the progenitor of the $z=0$ host galaxy of each
GC. We find that around half ($45~{}$per cent) of the GCs formed in situ.

Fig. 3 shows the probability distribution function (PDF) of GC formation times
for both in situ and accreted objects that survive to $z=0$. Overall, GCs in
E-MOSAICS have very early formation times, with more than half of the full
sample having formed within the first $2$ Gyrs after the Big Bang; this is
consistent with the inferred ages of the GCs in Fornax. The only exception is
Fornax 4 which is $\sim 10$ Gyrs old and is relatively ‘younger’ than the
other four GCs, as indicated by its more metal-rich stellar population (de
Boer & Fraser, 2016).

More interestingly, when comparing the accreted and in situ GC populations, we
find that $77$ per cent of accreted GCs formed early ($\sim 12$ Gyrs old),
while this fraction is only $23$ per cent for the in situ ones. The variation
in formation times between the two GC populations is expected, since there are
at least two processes that act differently in the two populations. The
accreted population is brought in by mergers, and most such mergers take place
at $z\gtrsim 2$; by construction, those GCs must have formed before the
merger. Secondly, as we shall discuss later, in situ GCs are more radially
concentrated and thus dynamical friction is more efficient at causing them
migrate to the centre, where they are tidally disrupted. This would suggest
that recently formed in situ GCs are more likely to survive to $z=0$ than
their older siblings. Both of these processes lead to an excess of younger in
situ GCs.

It is also worth noting that Fornax has two or more distinctive populations of
stars: a metal poor and a younger, more compact, metal rich population (e.g.
Battaglia et al., 2006). One way to obtain such configurations in simulations
is through a major merger between two dwarfs (e.g. Benítez-Llambay et al.,
2016; Genina et al., 2019b), raising the intriguing possibility that Fornax
might have had one or more such mergers (see also Yozin & Bekki, 2012).

Figure 3:  The distribution of formation times for GCs found in Fornax-mass
dwarfs at $z=0$. The plot shows the PDF for all GCs (black dotted line), in
situ (blue dashed line) and accreted (red solid line) GCs. Arrows indicate the
formation time of Fornax’s GCs (de Boer & Fraser, 2016). The left-hand side of
the plot corresponds to the Big Bang and $t=13.8$ Gyrs to present day.

##  4 GC orbital evolution

In the E-MOSAICS model GCs are associated with stellar particles whose
dynamics they follow. Like the stars, GCs are collisionless objects and thus
this is a reasonable approximation, except for one aspect: GCs can experience
significant dynamical friction as they orbit in their host galaxy. This
process is central to this study and thus we need to supplement the E-MOSAICS
model with a treatment of dynamical friction. We do this by post-processing
the simulated GC populations from E-MOSAICS as follows. We trace back each GC
associated with a present-day Fornax analogue to the simulation snapshot
closest to the time when the GC formed (for in situ objects) or to the time
when it was accreted into its $z=0$ host (for accreted GCs). The corresponding
positions and velocities, which are given by E-MOSAICS, are taken as the
starting position of each GC and its orbit is modelled in the potential of its
host galaxy accounting for dynamical friction. Finally, we follow the orbit of
each GC until the present day keeping track of whether it got close enough to
the galaxy centre to experience tidal disruption. In this section, we present
the model used to follow the GC orbit evolution.

###  4.1 Dynamical friction

We follow the evolution of each GC separately in the potential of the DM halo
and stellar component of its host galaxy. The DM halo is modelled as an
analytic Navarro, Frenk & White (Navarro et al., 1996b, 1997, hereafter NFW)
profile with total mass, $M_{200}$, and concentration, $c$, whose potential is
given by:

$\Phi_{\rm halo}=-\frac{GM_{200}}{r}\frac{\ln\left(1+r/r_{\rm
s}\right)}{\ln(1+c)-c/(1+c)}\;,$ (2)

where $r_{\rm s}$ denotes the scale radius of the DM halo. The halo
parameters, $c$, $r_{\rm s}$ and $M_{200}$, are inferred by fitting an NFW
profile to the host DM halo at $z=0$ for field Fornax-mass galaxies and at
infall for the Fornax-mass galaxies that are satellites at the present day. In
the latter case, we assume that the GC orbit is determined only by its Fornax-
mass host, and we ignore additional forces coming from the more massive system
to which the Fornax-mass satellite belongs. This approximation is valid as
long as the inner region, $r\lower 3.22916pt\hbox{$\sim$}\hbox
to0.0pt{\hss\raise 1.1625pt\hbox{$<$}}{}r_{\rm s}$, of the Fornax-mass
satellite does not undergo severe tidal disruption.

The GC dynamics are dominated by the DM halo potential, but, for completeness,
we also include the stellar potential. We model the stellar component of the
Fornax-mass galaxy as a Plummer profile with potential:

$\Phi_{\rm stars}=-\frac{GM_{\star}}{\sqrt{r^{2}+b^{2}}}\;,$ (3)

where $M_{\star}$ denotes the stellar mass and the parameter, $b$, is given by
$b=R_{1/2}/1.3$, with $R_{1/2}$ the 3D stellar half-mass radius. The Plummer
profile gives a good match to the stellar density profile of Fornax and other
dwarf spheroidal galaxies around the Milky Way (e.g. see Wang et al., 2018).

We implement dynamical friction as a deceleration experienced by the GC while
orbiting within the host halo of its galaxy. We assume that the deceleration
is given by Chandrasekhar’s formula,

$\frac{{\rm d}{\bf v}}{{\rm d}t}=-\frac{4\pi G^{2}M_{\rm
GC}\rho\ln\Lambda}{v^{2}}\left[{\rm
erf}(X)-\frac{2X}{\sqrt{\pi}}e^{-X^{2}}\right]\frac{{\bf v}}{v}\;,$ (4)

(Binney & Tremaine, 2008), where $G$ is the gravitational constant, $M_{\rm
GC}$ is the GC mass, ${\bf v}$ is the relative velocity of the GC with respect
to its galaxy, $\rho$ is the density of the DM halo at the GC’s position, and
$X=v/(\sqrt{2}\sigma_{v})$, with $\sigma_{v}$ the local 1D velocity dispersion
of the DM halo. We take the Coulomb logarithm as:

$\ln\Lambda=\ln\frac{b_{\rm{max}}\sigma_{v}^{2}}{GM_{\rm GC}}\;,$ (5)

(Goerdt et al., 2006), where $b_{\rm max}$ is the largest impact parameter to
be considered. By comparing with high resolution simulations of GC orbit
evolution, Goerdt et al. (2006) found that for cuspy haloes, $b_{\rm
max}=0.25~{}\rm kpc{}$ is the best fitting value. We include only the
dynamical friction arising from the DM halo, which is the dominant effect, and
neglect any contribution from the stellar distribution.

When integrating the orbits of the GCs we keep the potential of the DM halo
fixed at all times. This is a reasonable approximation since, while the halo
mass can grow by factors of several since the time most GCs formed, ${\sim}12$
Gyrs ago, the growth takes place by adding new mass to the outskirts of the
halo while leaving the inner region mostly unchanged (Wang et al., 2011). Most
stars and GCs orbit in the inner few kiloparsecs of the halo and thus their
orbits will not be affected by mass growth at the halo outskirts. However,
this approximation ceases to be valid when a halo undergoes major mergers. In
that case, the GC orbits can also be affected and, on average, they are pushed
slightly towards higher energy and more extended orbits (e.g. Benítez-Llambay
et al., 2016). This would reduce the effect of dynamical friction and slow
down the orbital decay of GCs. Thus, by not accounting for the effect of major
mergers, we are likely overestimating the number of GCs that sink to the
centre of their host.

###  4.2 Tidal disruption of GCs

To account for the disruption of GCs by the tidal field of their host galaxy
we calculate their tidal radius. For an NFW profile the tidal radius for a GC
on a circular orbit at distance, $r$, from the galaxy centre is given by:

$r_{\rm tidal}=(r+r_{\rm s})\left(\frac{M_{\rm GC}}{M_{\rm
200}}\frac{r}{3r+r_{\rm s}}\right)^{\frac{1}{3}}\;,$ (6)

(Renaud et al., 2011; Orkney et al., 2019). This expression neglects the tides
arising from the stellar distribution of the host galaxy, which are much
smaller than the tidal field of the DM halo. We calculate the tidal radius at
each point along the GC’s orbit and consider the GC to be destroyed by the
tidal field of its host when the tidal radius is comparable to the half-mass
radius of the GC. As GCs lose mass due to two-body relaxation as well as tidal
stripping, their half-mass radius increases. According to the N-body
simulations of Orkney et al. (2019), GCs evolving in cuspy profiles can reach
a half-mass radius of ${\sim}6~{}\rm pc$, after which the mass loss rate
increases rapidly and the GC is disrupted shortly thereafter. Based on these
results, we assume that GCs are fully disrupted if at any point along their
orbit the tidal radius becomes smaller than $6~{}\rm pc$.

###  4.3 Decay of the GC orbits

Figure 4:  The time evolution of the radial distance for two GCs on circular
orbits whose starting positions are at $0.5~{}\rm kpc$ and $1~{}\rm kpc$ from
the centre of its host. The top panel shows how the orbits change when varying
the host halo mass (here, the GCs have the same mass, $M_{\rm GC}=2\times
10^{5}{~{}\rm M_{\odot}}$). The bottom panel shows how the orbits change when
varying the GC mass (here, the host halo has the same mass, $M_{\rm
200}=3\times 10^{10}{~{}\rm M_{\odot}}$).

The orbital decay time of a GC is mainly determined by two factors: 1) the
initial distance to the center of its host, and 2) the mass ratio between the
GC and the host halo. We illustrate this in Fig. 4, where we show GC orbits
for two starting positions and a selection of GC and host halo masses. In all
cases, we start with GCs on circular orbits; in reality, E-MOSAICS predicts
highly elliptical orbits even for the in situ GCs of our Fornax-analogue
sample, with a median ellipticity at birth of 0.64 (and 16 to 84 per centile
of the ellipticity distribution of 0.36 and 0.84 respectively). However, here
we want to illustrate the systematic effect of dynamical friction, which is
most clearly seen for circular orbits.

In the first example, in the top panel in Fig. 4, we fix the GC mass to
$M_{\rm GC}=2\times 10^{5}{~{}\rm M_{\odot}}$, which is the median mass of
Fornax’s GCs, and vary the mass of the host halo. We select a reference halo
mass, $M_{200}=3\times 10^{10}{~{}\rm M_{\odot}}$, which corresponds to the
median halo mass in our Fornax analogue sample, as well as a much lower mass,
$M_{200}=2\times 10^{9}{~{}\rm M_{\odot}}$, which was the value used by Goerdt
et al. (2006). We also present orbits for host halo masses $M_{200}=1$ and
$5\times 10^{10}{~{}\rm M_{\odot}}$, which correspond to roughly the 10 and 90
per centile of the distribution of halo masses for our sample of Fornax-mass
dwarfs.

We show the orbital evolution for two starting distances from the host centre:
$0.5~{}\rm kpc$ and $1~{}\rm kpc$. As expected, GCs that are initially closer
to the centre sink more quickly than those which start further away. GCs with
$r_{\rm init}=0.5~{}\rm kpc$ sink to the center within $\sim 2$ Gyr in the
least massive halo ($M_{\rm 200}=2\times 10^{9}{~{}\rm M_{\odot}}$), and $\sim
5$ Gyr in the most massive host halo ($M_{\rm 200}=5\times 10^{10}{~{}\rm
M_{\odot}}$). The sinking time of GCs increases when increasing $r_{\rm
init}$. For $r_{\rm init}=1~{}\rm kpc$, the GC can survive for a Hubble time
in the two most massive hosts, but not in the lower mass haloes.

In view of the diversity of the mass of the Fornax GCs, in the bottom panel we
study the variation in orbits with GC mass. In this case, we fix the halo mass
to $M_{200}=3\times 10^{10}{~{}\rm M_{\odot}}$, the median value for Fornax-
mass hosts in the E-MOSAICS simulation (see Fig. 3). We show orbits for four
GC masses, from $4$ to $50\times 10^{4}{~{}\rm M_{\odot}}$. The lowest mass GC
experiences the weakest orbital decay, and survives for a Hubble time even
when starting at $r_{\rm init}=0.5~{}\rm kpc$. In contrast, the most massive
GC sinks rapidly to the centre, with a sinking time of $2$ and $7$ Gyrs
depending on whether it starts at 0.5 or $1~{}\rm kpc$ from the centre. This
suggests that the most massive GC in Fornax, Fornax-3, must have formed
significantly farther than its present-day position ($0.46~{}\rm kpc$
projected distance); we discuss this point in greater detail in Sections 6 and
in Fig. 10.

The oldest GCs can lose up to 60 per cent of their initial mass due to
dynamical effects, such as two-body relaxation, tidal stripping and stellar
evolution (see Appendix C). We account for this effect by having a time-
varying mass for our GCs, with the mass decreasing in time according to the
prescription detailed in Appendix C.

##  5 Results

We now compare the number and radial distribution of GCs in the Fornax dSph
with the E-MOSAICS predictions. These predictions have been post-processed to
include dynamical friction and exclude GCs that were tidally disrupted after
sinking to the centre of their host, as described in the previous section.

Figure 5:  Top panel: PDF of the number of GCs at $z=0$ in Fornax analogues.
We only show GCs with present-day masses above $4\times 10^{4}{~{}\rm
M_{\odot}}$, the lowest mass GC in Fornax. The three lines show the
distribution for Fornax-mass dwarfs in the field (dotted line), for satellites
around all hosts (dashed line) and for satellites around MW-mass hosts (solid
line). Bottom: same as the top panel but counting also the GCs more massive
than $4\times 10^{4}{~{}\rm M_{\odot}}$ that were tidally disrupted after
sinking to the centre of their host galaxy.

###  5.1 The number of GCs in Fornax-mass dwarfs

The first question we address is how frequently do the simulations produce
dwarfs with the number of GCs observed in the Fornax dSph, that is at least
five with mass, $M_{\rm GC}>4\times 10^{4}{~{}\rm M_{\odot}}$? We answer this
question in the top panel of Fig. 5, which shows the distribution of the
number of GCs that survive to $z=0$ in our Fornax analogue sample. Regardless
of whether they are satellites or field galaxies, nearly $35~{}$per cent of
the Fornax-mass hosts do not have any GCs at all and $30~{}$per cent have only
one object. Thus, E-MOSAICS predicts a low number of GCs in Fornax-mass
dwarfs. Field dwarfs tend to have more GCs than satellite galaxies, but the
effect is rather small. The difference could be due to: i) satellite dwarfs
having stopped forming stars, and hence GCs, after falling into their host
(although infall can lead to at least a temporary increase in star formation,
e.g. Shao et al., 2018; Genina et al., 2019b; Hughes et al., 2019), and ii)
having had some of the GCs that were on extended orbits tidally stripped by
their more massive host.

When comparing with the Fornax dSph we find that only $3~{}$per cent of the
simulated analogues have five or more GCs. Thus, systems of GCs as rich as the
one observed in Fornax are predicted to be rather rare. As we mentioned in the
previous section, around $33~{}$per cent of the GCs that formed within, or
were accreted into, our Fornax analogues sunk to the centre of these systems
and were tidally destroyed. So, could the rarity of GC-rich dwarfs in our
model be due to tidal destruction of GCs? We answer this question in the
bottom panel of Fig. 5 where we again show the PDF of the GC count for each
Fornax-mass dwarf, but now we include also the tidally destroyed GCs. Thus,
the panel shows the E-MOSAICS prediction for the PDF of all the GCs that ever
formed within, or fell into, a present-day Fornax-mass dwarf.

As expected, we find an increase in the number of GCs per dwarf galaxy. In
particular, we find a decrease in the fraction of dwarfs with none or one GC,
and an increase in systems with two or more GCs. However, the tail of the
distribution corresponding to five or more GCs changes only slightly, with the
chance of having a GC system as rich as the Fornax dSph increasing to only
$6.5~{}$per cent. This value is still rather low and indicates that the tidal
destruction of GCs is not the main factor responsible for the low prevalence
in E-MOSAICS of GC systems as rich as the one in Fornax. The main difference
between haloes with a central core and those with a cusp is the fraction of
tidally disrupted GCs, which is expected to be lower for the former case (e.g.
Goerdt et al., 2006; Meadows et al., 2020). Thus, the number of Fornax GCs is
not a reflection of whether or not the dark matter halo profile has a central
core or cusp.

The number of GCs in a galaxy depends strongly on the mass cut used for the GC
selection. Here, we only select GCs with a present-day stellar mass above
$4\times 10^{4}{~{}\rm M_{\odot}}{}$, which is the mass of the lightest GC in
Fornax. The uncertainty in the mass estimate of that GC is ${\sim}1\times
10^{4}{~{}\rm M_{\odot}}$ (25 per cent fractional error), as shown in Table 1.
To have a fair comparison between our model prediction and the data we should
include this measurement error in our predictions. We recalculate the fraction
of hosts that have five or more GCs by adding a random error to the GC mass,
given by a Gaussian distribution with zero mean and standard deviation,
$1\times 10^{4}{~{}\rm M_{\odot}}$. We find that the inclusion of errors leads
to a small increase in the fraction of GC systems as rich as the one in
Fornax. For our fiducial case that includes GC tidal disruption, the
prevalence of Fornax GC systems grows from 3 per cent, when no mass error is
included, to 4 per cent when modelling mass measurement errors.

The E-MOSAICS predictions have been shown to agree well with observational
data (Kruijssen et al., 2019a) such as the total mass in GCs at a given host
halo mass or host galaxy stellar mass (for details see Fig. 1 in Bastian et
al., 2020). This suggests that the paucity of rich Fornax-like GC systems is
unlikely to be due to a failure of the GC formation and evolution model and
that, most likely, suggests that Fornax has an excess of GCs for its stellar
mass.

A census of GCs in galaxies as faint as Fornax is currently lacking and the
available observations are rather heterogeneous, which makes it difficult to
perform a robust comparison between Fornax and other equal-mass dwarfs.
However, from the currently available data in Forbes et al. (2018), which
provides a table of all dwarf galaxies with one or more GCs, we conclude that
our predictions in Fig. 5 are broadly consistent with observations. For
example, the table of Forbes et al. contains 9 dwarfs with luminosity within
$\pm 1~{}\rm{mag}$ of the Fornax dSph. Of those, the majority (7 out of 9)
have either one or two GCs, and only two systems, Fornax and UGC685, have more
than two (interestingly, both these galaxies have 5 GCs).

Figure 6: The 3D radial distribution of GCs in Fornax-mass dwarfs. The dashed
line shows the distribution of GCs at birth, if they formed in situ, or at
infall into their Fornax analogue host, if they were accreted. Dynamical
friction leads to some GC sinking to the centre of their host and the
resulting $z=0$ distribution is shown by the dotted line. Once a GC gets close
enough to the centre, it is disrupted by the tides of the host galaxy. The
resulting $z=0$ radial distribution is shown by the solid line and includes
the effect of both dynamical friction and tidal disruption.

###  5.2 The radial distribution of GCs

Since, as we have seen, the high number of GCs in Fornax is not in itself a
reliable way to infer if the DM halo is cuspy or has a core. In this section
we examine the constraints that can be placed on the inner structure of Fornax
by considering instead the radial distribution of its GCs. Most dwarfs have
one or two GCs, so it is difficult and not very meaningful to show the radial
distribution for each host. Instead, we stack all our Fornax-mass dwarfs and
study the mean radial profile of GCs obtained by stacking the distances to
their host’s centre.

We start with Fig. 6 which shows how the 3D radial distribution of GCs changes
from their initial position to the present day. The dashed line shows the
initial positions of GCs which correspond to either their birth location (for
in situ GCs) or to their position at the time they were accreted into their
$z=0$ host (for accreted GCs). There is a large range of initial positions,
with a median value of 1.7 kpc and 10 and 90 per centiles of 0.3 and 8.5 kpc
respectively. Although not shown, we have checked that the in situ population
is more concentrated than the accreted sample.

The subsequent orbital evolution leads to the GCs moving inwards on average,
as illustrated by the dotted line in Fig. 6, which shows the $z=0$ positions
after including the effect of dynamical friction but without removing tidally
destroyed GCs. We find that around one third of the GC sample sinks to the
centre of their galaxies, i.e. to $r<0.1~{}\rm kpc$, and thus are likely to be
tidally disrupted. When accounting for tidal disruption, we find that the
majority of GCs with $r\lesssim 0.2~{}\rm kpc$ are removed and the
distribution of surviving GCs, which is shown by the solid line in Fig. 6, is
now somewhat less concentrated than the initial GC positions.

Figure 7: The dependence of the 3D radial distribution of GCs on the number of
GCs in a Fornax-mass dwarf. It shows that there is no correlation between the
radial distribution of GCs and the number of GCs in a galaxy.

In Fig. 7 we investigate if GC systems in all Fornax-mass dwarfs have the same
average radial profile. In particular, we study if the radial profile depends
on the number of GCs in a given galaxy. This test is motivated by our previous
result that Fornax has a large number of GCs for its stellar mass and we wish
to investigate if such an excess also impacts the GC radial distribution. The
figure shows that there is no significant correlation between the number of
GCs and their radial distribution. (While the dashed line shows some
deviations from the mean trend, these are consistent with random scatter given
the small sample size.) In practice, this shows that no systematic bias is
introduced when comparing the Fornax GC distribution with the average profile
of the full Fornax analogue sample.

Figure 8: Comparison of the present-day radial distribution of GCs in field
(solid lines) and satellite (dashed lines) Fornax-mass hosts. For each
subsample we show the distribution of all GCs (in black), those formed in situ
(in blue) and those accreted (in red).

Fig. 8 investigates how the radial distribution of GCs differs in satellite
and field dwarfs. GCs in satellite dwarfs (black-dashed line) are more
radially concentrated than those in field dwarfs (black-solid line): nearly
half of the GCs in satellites are located within $1$ kpc from the centre while
only 25 per cent in the field dwarfs are. The difference between satellite and
field dwarfs is due to tidal stripping of GCs in satellite dwarfs, which
preferentially removes objects on extended orbits and thus leads to more
concentrated distributions. While not shown, we have checked this by tracing
back the satellite and field dwarfs to redshift, $z=2$, and comparing their GC
distributions at that time. To make a meaningful comparison, for the $z=0$
satellites, we used only the $z=2$ progenitors that were central galaxies,
which represent the bulk of the population. We find that at redshift $z=2$
there is no significant difference in the GC distributions of progenitors of
satellite and field dwarfs. Thus, the differences seen in Fig. 8 are mainly
due to tidal stripping from the outskirts of satellite dwarfs.

In Fig. 8 we further split the GCs into two additional subsamples: accreted
and in situ. In both satellite and field Fornax analogues, the in situ GCs are
more radially concentrated than the accreted population. The segregation
between the in situ and accreted components is largest for field dwarfs. The
fraction of accreted GCs in satellite dwarfs is $50~{}$per cent, lower that
the $60~{}$per cent fraction in field dwarfs. This difference could be due to
two effects. First, satellite galaxies experience mergers mostly before infall
into a more massive host halo (see e.g. Angulo et al., 2009; Deason et al.,
2014); they thus have less time to devour other galaxies and grow their
accreted GC population. Secondly, as we have just seen, accreted GCs have a
more radially-extended distribution than the in situ GCs and so they are more
easily tidally stripped from their hosts.

Figure 9: The distribution of projected 2D distances of GCs in Fornax
analogues. The solid and dashed lines show the results for respectively field
and satellite Fornax-mass dwarfs. The observed projected distances of GCs in
Fornax are shown by the dotted line (each triangle symbol corresponds to a
Fornax GC). The distribution of GCs in Fornax matches very well that of its
analogues in the simulation.

Finally, we compare the GC radial distribution of our Fornax analogues with
that of the real Fornax dSph. We take the projected distances of Fornax GCs
from Mackey & Gilmore (2003), which we recalculate using a more recent
distance estimate for Fornax of $147\pm 4$ kpc (McConnachie, 2012); this new
distance is about $7~{}$per cent larger than the 137 kpc used by Mackey &
Gilmore. The updated GC distances are given in Table 1 and their radial
distribution is shown in Fig. 9, where we also show the projected distance
distribution predicted by our model for satellite and field dwarfs.
Qualitatively, we find good agreement between observations and the
simulations. To quantify the degree of agreement, we use the Kolmogorov-
Smirnov (KS) test. The KS test returns a p-value, $p=0.98$, when comparing the
Fornax GCs with those of satellite analogues in the E-MOSAICS simulation, and
a lower value, $p=0.55$, when comparing with the field analogues. This
indicates that any differences seen in Fig. 9 are consistent with random noise
and thus not significant.

##  6 Discussion

In this section we discuss in detail the implications of our results and what
they can tell us about the formation history of the Fornax dSph.

###  6.1 The DM halo of the Fornax dSph

Our stellar-mass selected Fornax analogues are found in haloes of present-day
mass, $M_{200}=2.8^{+1.1}_{-1.0}\times 10^{10}{~{}\rm M_{\odot}}$ (68 per cent
confidence limit), and concentration, $c=10.5^{+3.6}_{-2.6}$. These values
apply to field galaxies, i.e. to Fornax-mass dwarfs that are central galaxies.
The $M_{200}$ and $c$ values can be used to infer the halo scale radius,
$r_{s}$, and characteristic density, $\rho_{0}$ (see Eq. 1). The resulting
values of $r_{s}$ and $\rho_{0}$ provide a good description also of the inner
region ($r\lesssim r_{s}$) of our Fornax-mass satellites, as long as these
satellites have not experienced stripping of their inner DM profile. This is
because DM halo growth proceeds primarily by the deposition of newly accreted
material in the outer parts of the halo with the inner parts remaining largely
unchanged (Wang et al., 2011). Once a galaxy becomes a satellite, it stops
accreting mass and, in fact, can lose mass by tidal stripping, with most of
the decrease taking place at the outskirts of the halo while the inner region
changes only slowly (e.g. Peñarrubia et al., 2008; Errani & Peñarrubia, 2020).
The profile of the inner halo changes appreciably only once a large fraction
($\gtrsim 50~{}$per cent) of the mass has been lost, and this will be manifest
in changes in the orbits of stars and GCs (Peñarrubia et al., 2008; Fattahi et
al., 2016). Thus, as long as the satellites have not suffered a large degree
of tidal stripping, we expect that they should have similar inner DM profiles
as their field counterparts (we explicitly show this for the halo scale
radius, $r_{s}$, in Fig. 2).

We have found that the Fornax dSph has an atypically large number of GCs for
its stellar mass. Observations have shown that the number and, in particular,
the total mass in GCs correlates with halo mass (e.g. Forbes et al., 2018).
The same correlation is well reproduced in the E-MOSAICS simulation (Bastian
et al., 2020). This raises the question of whether we can take into account
the richness of the Fornax GC system to better determine properties of its
dark halo. We investigate this possibility in Appendix B, where we show that
our Fornax analogues that have at least three GCs (the fraction of systems
with at least five GCs is too low to obtain robust conclusions) reside in DM
haloes that are, on average, 20 per cent more massive than the full Fornax-
mass sample. Interestingly, the DM halo concentration does not show a trend
with GC abundance.

Our average, our DM halo mass estimates are an order of magnitude larger than
previous determinations based on the stellar kinematics of the Fornax dSph,
which suggest a DM halo mass of $\rm{few}\times 10^{9}{~{}\rm M_{\odot}}$
(Goerdt et al., 2006; Boldrini et al., 2019; Meadows et al., 2020). However,
our estimates are in good agreement with the recent non-parametric Jeans
modelling of Read et al. (2019), who found a total mass of ${\sim}2\times
10^{10}{~{}\rm M_{\odot}}$. The discrepancy could be due to the Fornax halo
having experienced a large amount of tidal stripping, as suggested by Fattahi
et al. (2016), Wang et al. (2017) and Genina et al. (2020). These studies used
hydrodynamical simulations to select Fornax counterparts by matching the
observed stellar velocity-dispersion and stellar half-mass radius in the first
two of these papers and the stellar mass in the last one. If stripping is
important then a ${\sim}10^{9}{~{}\rm M_{\odot}}$ halo is not representative
of the halo in which the Fornax GCs formed and evolved for most of their
lifetime before the galaxy fell into the Milky Way.

Figure 10: The orbital evolution of the GCs in Fornax assuming the median
mass, $M_{200}=3\times 10^{10}{~{}\rm M_{\odot}}$, and the median
concentration, $c=10$, for the DM halo of Fornax’s. Negative time values
correspond to past orbits and positive ones to future orbits. Each colour
corresponds to one of the 5 GCs in Fornax (see bottom-right legend). In the
past the GCs were farther from the galaxy centre, especially Fornax 3, which
is the most massive, and Fornax 4, which today is the closest to the centre.

###  6.2 What determines the sinking time scale for GCs?

Our result agrees with earlier work that showed that in a cuspy DM halo, the
orbits of GCs decay because of dynamical friction (e.g. Goerdt et al., 2006;
Meadows et al., 2020). However, we estimate a longer decay time for the Fornax
GCs that in previous work. The differences arise from two novel aspects
introduced into this analysis. Firstly, we show that previous measurements of
the Fornax DM halo based on simple stellar kinematic analyses give a mass that
is too low (see discussion in Section 6.1); this leads to an estimate of the
sinking timescale that is correspondingly too low (see Section 4). On average,
the sinking timescale doubles in an $M_{200}=3\times 10^{10}{~{}\rm
M_{\odot}}$ halo compared to that in a halo of the previously assumed value of
$M_{200}=2\times 10^{9}{~{}\rm M_{\odot}}$.

Secondly, as pointed out in previous work, we know only the present-day
positions of the GCs, not their birth positions (e.g. Angus & Diaferio, 2009;
Boldrini et al., 2019; Meadows et al., 2020). The sinking time increases
rapidly with increasing distance, and even small changes in distance can have
a considerable effect. To take this into account in our analysis, we used the
GC positions and velocities at birth from the E-MOSAICS simulation which self-
consistently models the formation and evolution of GCs in cuspy DM haloes. Our
calculation shows that the majority (67 per cent) of GCs in Fornax-mass dwarfs
survive to the present day and thus only a third are expected to sink to the
centre of their galaxy in a Hubble time. The survival chance depends on the
starting position of the GC: it is as low as 37 per cent for GCs formed within
1 kpc from the centre and as large as 92 per cent for GCs formed beyond 1 kpc
from the centre.

###  6.3 Past and future orbits of the GCs in Fornax

Given the updated GC sinking times discussed in Section 6.2, we recalculate,
in Fig. 10, the possible orbits of the Fornax GCs. The figure shows both the
past (to the left of the vertical dotted line) and the future orbits assuming
the median halo mass, $M_{\rm 200}=3\times 10^{10}{~{}\rm M_{\odot}}$, and the
median concentration, $c=10$, of our Fornax analogue sample. The distances of
GCs are in 3D and the orbit integration was done assuming circular orbits with
present-day radii of $\sqrt{3/2}$ times their observed projected radial
distance.

This simple model suggests that only one GC, Frn 4, was born at a small radial
distance, $r_{\rm init}=700~{}\rm pc$, which is approximately the half-light
radius of the present-day Fornax (as given by McConnachie 2012). For the
remaining four GCs, the initial distances are larger than $1.5$ kpc, which
suggests that at least some of them were accreted (see also Boldrini et al.,
2020). This result is in good agreement with our model prediction that around
half of the present-day surviving GCs are of accreted origin. The sinking time
varies from GC to GC depending on the initial distances and masses; Frn 1, Frn
2, and Frn 5 could still survive for another Hubble time mainly because of
their large initial distances. Interestingly, the time for Frn 3 to sink is
the shortest amongst the five even though it started from $r_{\rm
init}=1.5~{}\rm kpc$, which is twice as far as Frn 4. This is due to the fact
that Frn 3 is the most massive GC of the five, with a mass, $M_{\star}=5\times
10^{5}{~{}\rm M_{\odot}}$.

We note that the orbital evolution in Fig. 10 is highly simplified and does
not take into account that the GCs are likely to be on eccentric orbits. For
example, E-MOSAICS predicts that at birth the GCs in Fornax-mass dwarfs have
an orbital ellipticity of $0.64^{+0.20}_{-0.28}$ (68th per centile).
Furthermore, there is indirect evidence that the Fornax dSph might have had a
major merger ${\sim}10$ Gyrs ago (Yozin & Bekki, 2012; Benítez-Llambay et al.,
2016) and that it is potentially undergoing severe tidal stripping (Wang et
al., 2017; Genina et al., 2020). The latter process, if it takes place rapidly
enough, could lead to the dissolution of the GC system in the outskirts of the
Milky Way halo before Frn 3 can sink to the centre.

###  6.4 GCs and the core-cusp problem

The existence of GCs close to the centre of Fornax has often been adduced as
evidence that this galaxy’s DM halo has a constant density core. Here we have
shown, to the contrary, that cuspy DM haloes with galaxies of simliar stellar
mass to Fornax are expected to have a similar GC radial distribution: the GCs
in Fornax cannot be used to rule out a cuspy halo (see also Angus & Diaferio,
2009; Meadows et al., 2020). In fact, the absence of an extended core, of
radius $\sim 1~{}\rm kpc$, in Fornax is supported by other probes. Read et al.
(2019) used higher-order moments of the velocity distribution to break the
degeneracy between mass and velocity anisotropy and found that the stellar
kinematics in Fornax favour an increase of an order of magnitude in DM density
between radii of 1 and 0.1 kpc, with a core, if present, having a radius
$\lesssim 0.2~{}\rm kpc$. Such a core radius is too small to lead to a
stalling of the orbital decay of GCs, which takes place at around ${\sim}1/3$
of the core radius (Meadows et al., 2020).

Our model predicts that a third of the GCs associated with Fornax-mass dwarfs
have sunk to the centre of their host galaxies where they may have been
tidally disrupted. It remains to be seen what is the corresponding fraction of
destroyed GCs if the DM profile had a core, but it is likely lower. Thus,
identifying in observations the fraction of GCs that were destroyed is, in
principle, a promising avenue for distinguishing between a cusp and a core in
Fornax. The extent of tidal stripping experienced by GC is sensitive to the
central DM density profile: in cuspy haloes GCs are more extended and have
lower mass loss that in haloes with cores (e.g. Amorisco, 2017; Webb &
Vesperini, 2018; Orkney et al., 2019) but it is unclear if this difference can
be used as a test of the inner structure of the halo.

##  7 Conclusions

We have investigated the number and radial distribution of GCs in the Fornax
dSph with the goal of testing whether the GC population in this galaxy is
consistent with a cuspy halo, the simplest prediction of the $\Lambda$CDM
cosmological model (Navarro et al., 1996b, 1997). Our study has been motivated
in part by the ongoing debate on whether the Fornax GCs are a robust signature
of a kiloparsec-sized core at the centre of this galaxy’s halo (e.g. see
discussion in Meadows et al., 2020). To this end we have analyzed the
E-MOSAICS simulation (Pfeffer et al., 2018; Kruijssen et al., 2019a), which
includes a subgrid prescription for self-consistently following galaxy and GC
formation and evolution within the EAGLE cosmological hydrodynamics framework
(Schaye et al., 2015; Crain et al., 2015).

We have proceeded by identifying a sample of Fornax analogues, selected to
have a similar stellar mass, within a factor of two, as the Fornax dSph, i.e.
in the range $[2,8]\times 10^{7}{~{}\rm M_{\odot}}$. E-MOSAICS is well suited
for this study because it contains a large sample (1154) of Fornax-mass
dwarfs, within which GC formation and evolution is followed in a full
cosmological context. To account for the GC orbital decay due to dynamical
friction, we have post-processed the orbits of GCs starting from their initial
positions as determined by E-MOSAICS, and assuming that they evolve in the NFW
profile that best fits each host halo in the simulation. The initial
conditions are given by the GC positions and velocities at birth for in situ
GCs, and at infall for accreted GCs. We applied a tides model to account for
the disruption of GC by their host halo. We then investigated the number and
radial distribution of GCs of stellar mass, $M_{\star}\geq 4\times
10^{4}{~{}\rm M_{\odot}}$, corresponding to the smallest of the five well-
studied GCs in Fornax. Our main conclusions may be summarised as follows:

  1. 1.

Field Fornax-mass dwarfs reside in haloes of median mass, $M_{\rm 200}\approx
2.8\times 10^{10}{~{}\rm M_{\odot}}$, and concentration, $c\approx 10$ (see
Fig. 1).

  2. 2.

The population of GCs of Fornax analogues consists, in nearly equal amounts,
of in situ (45 per cent) and accreted objects (55 per cent).

  3. 3.

The GCs of Fornax are typically old, with a median age, $t_{\rm age}\approx
12$ Gyrs. The fraction of GCs younger than $12$ Gyrs is 77 per cent for in
situ GCs and 23 per cent for accreted GCs. This is because in situ GCs, which
form close to their host centre, are more likely to be destroyed, while the
accreted GCs started orbital decay at larger initial distances (see Fig. 3).

  4. 4.

Orbital decay leads to 33 per cent of GCs sinking to the centre where they can
be tidally disrupted in a cuspy halo (see Fig. 6).

  5. 5.

Most of the Fornax-mass dwarfs have one ($\approx$ 30 per cent) or no
($\approx$ 35 per cent) GCs, with very few having several. The results are
similar for both satellite and field dwarfs (see Fig. 5).

  6. 6.

The median radial distance of surviving GCs is $1.7~{}\rm kpc$, with 10 and 90
percentiles of $0.3$ and $8.5~{}\rm kpc$ respectively (see Fig. 6).

  7. 7.

Satellite galaxies have more concentrated GC distributions than their field
analogues. In situ GCs are also more radially concentrated than their accreted
counterparts (see Fig. 8).

Our model predicts that only a small fraction (3 per cent) of Fornax-mass
dwarfs have at least as many GCs as are observed in Fornax. This result does
not pose a challenge to our model since observations similarly show that
Fornax has a surprisingly large number of GCs for its stellar mass (e.g. see
Figure C1 in Forbes et al. 2018). In fact, most observed dwarfs with similar
masses to Fornax have, at most, only one or two GCs, in agreement with our
model predictions.

The observed radial distribution of the GCs in Fornax agrees very well with
our model prediction for satellite Fornax analogues (KS test p-value of 0.98).
This indicates that, contrary to previous claims, the distribution of GCs in
Fornax can be reproduced in a DM halo that has a central cusp. Furthermore,
our model predicts that the $z=0$ surviving GCs represent ${\approx}$ 67 per
cent of the GCs ever associated with Fornax, ie. born within, or accreted
onto, Fornax and thus suggests that Fornax may have had an additional
${\sim}2$ GCs that sunk to its centre and could have been destroyed. One such
candidate has recently been discovered: a low mass GC ${\sim}30$ pc away in
projection from the centre of Fornax which has an irregular shape that could
indicate that it is undergoing tidal disruption (Wang et al., 2019).

By modelling GC formation and evolution in a cuspy DM halo, our work has
demonstrated that the Fornax GC system cannot be used to rule out a central
cusp in this galaxy’s halo. The next steps involve investigating if the number
and radial distribution of GCs differs between DM haloes with cusps and cores.
One way to do so would be to run E-MOSAICS with an increased star formation
gas density threshold, which has been shown to lead to the formation of cores
within the EAGLE galaxy formation model (Benítez-Llambay et al., 2019). Such a
simulation would predict the GC distribution for profiles with cores and, by
comparing directly against our results, would help identify which GC
statistics are best suited to distinguish between a central core or cusp in
the DM profile.

## Acknowledgements

SS and CSF were supported by the European Research Council through ERC
Advanced Investigator grant, DMIDAS [GA 786910], to CSF and by the Science and
Technology Facilities Council (STFC) [grant number ST/F001166/1, ST/I00162X/1,
ST/P000541/1]. MC acknowledges support by the EU Horizon 2020 research and
innovation programme under a Marie Skłodowska-Curie grant agreement 794474
(DancingGalaxies). MRC acknowledges financial support from a CITA National
Fellowship. AD and RAC are supported by Royal Society University Research
Fellowships. JP gratefully acknowledges funding from a European Research
Council consolidator grant (ERC-CoG-646928-Multi-Pop). MRC and JMDK gratefully
acknowledge funding from the European Research Council (ERC) under the
European Union’s Horizon 2020 research and innovation programme via the ERC
Starting Grant MUSTANG (grant agreement number 714907). JMDK gratefully
acknowledges funding from the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) through an Emmy Noether Research Group (grant number
KR4801/1-1) and the DFG Sachbeihilfe (grant number KR4801/2-1). This work used
the DiRAC Data Centric system at Durham University, operated by ICC on behalf
of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by
BIS National E-infrastructure capital grant ST/K00042X/1, STFC capital grant
ST/H008519/1, and STFC DiRAC Operations grant ST/K003267/1 and Durham
University. DiRAC is part of the National E-Infrastructure.

## Data availability

The used and data produced in this paper are available upon reasonable request
to the corresponding author.

## References

  * Amorisco (2017) Amorisco N. C., 2017, ApJ, 844, 64
  * Amorisco et al. (2013) Amorisco N. C., Agnello A., Evans N. W., 2013, MNRAS, 429, L89
  * Angulo et al. (2009) Angulo R. E., Lacey C. G., Baugh C. M., Frenk C. S., 2009, MNRAS, 399, 983
  * Angus & Diaferio (2009) Angus G. W., Diaferio A., 2009, MNRAS, 396, 887
  * Bastian (2008) Bastian N., 2008, MNRAS, 390, 759
  * Bastian et al. (2020) Bastian N., Pfeffer J., Kruijssen J. M. D., Crain R. A., Trujillo-Gomez S., Reina-Campos M., 2020, preprint,  (arXiv:2005.05991)
  * Battaglia et al. (2006) Battaglia G., et al., 2006, A&A, 459, 423
  * Battaglia et al. (2008) Battaglia G., Helmi A., Tolstoy E., Irwin M., Hill V., Jablonka P., 2008, ApJ, 681, L13
  * Benítez-Llambay et al. (2016) Benítez-Llambay A., Navarro J. F., Abadi M. G., Gottlöber S., Yepes G., Hoffman Y., Steinmetz M., 2016, MNRAS, 456, 1185
  * Benítez-Llambay et al. (2019) Benítez-Llambay A., Frenk C. S., Ludlow A. D., Navarro J. F., 2019, MNRAS, 488, 2387
  * Binney & Tremaine (2008) Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edition. Princeton University Press
  * Boldrini et al. (2019) Boldrini P., Mohayaee R., Silk J., 2019, MNRAS, 485, 2546
  * Boldrini et al. (2020) Boldrini P., Mohayaee R., Silk J., 2020, MNRAS, 492, 3169
  * Booth & Schaye (2009) Booth C. M., Schaye J., 2009, MNRAS, 398, 53
  * Bose et al. (2019) Bose S., et al., 2019, MNRAS, 486, 4790
  * Cole et al. (2012) Cole D. R., Dehnen W., Read J. I., Wilkinson M. I., 2012, MNRAS, 426, 601
  * Crain et al. (2015) Crain R. A., et al., 2015, MNRAS, 450, 1937
  * Dalla Vecchia & Schaye (2012) Dalla Vecchia C., Schaye J., 2012, MNRAS, 426, 140
  * Davis et al. (1985) Davis M., Efstathiou G., Frenk C. S., White S. D. M., 1985, ApJ, 292, 371
  * Deason et al. (2014) Deason A., Wetzel A., Garrison-Kimmel S., 2014, ApJ, 794, 115
  * Di Cintio et al. (2014) Di Cintio A., Brook C. B., Macciò A. V., Stinson G. S., Knebe A., Dutton A. A., Wadsley J., 2014, MNRAS, 437, 415
  * Dolag et al. (2009) Dolag K., Borgani S., Murante G., Springel V., 2009, MNRAS, 399, 497
  * Errani & Peñarrubia (2020) Errani R., Peñarrubia J., 2020, MNRAS, 491, 4591
  * Fattahi et al. (2016) Fattahi A., et al., 2016, MNRAS, 457, 844
  * Flores & Primack (1994) Flores R. A., Primack J. R., 1994, ApJ, 427, L1
  * Forbes et al. (2018) Forbes D. A., Read J. I., Gieles M., Collins M. L. M., 2018, MNRAS, 481, 5592
  * Genina et al. (2018) Genina A., et al., 2018, MNRAS, 474, 1398
  * Genina et al. (2019a) Genina A., et al., 2019a, preprint,  (arXiv:1911.09124)
  * Genina et al. (2019b) Genina A., Frenk C. S., Benítez-Llambay A. r., Cole S., Navarro J. F., Oman K. A., Fattahi A., 2019b, MNRAS, 488, 2312
  * Genina et al. (2020) Genina A., Read J. I., Fattahi A., Frenk C. S., 2020, preprint, p. arXiv:2011.09482 (arXiv:2011.09482)
  * Gilmore et al. (2007) Gilmore G., Wilkinson M. I., Wyse R. F. G., Kleyna J. T., Koch A., Evans N. W., Grebel E. K., 2007, ApJ, 663, 948
  * Goerdt et al. (2006) Goerdt T., Moore B., Read J. I., Stadel J., Zemp M., 2006, MNRAS, 368, 1073
  * Harris et al. (2017) Harris W. E., Blakeslee J. P., Harris G. L. H., 2017, ApJ, 836, 67
  * Hellwing et al. (2016) Hellwing W. A., Frenk C. S., Cautun M., Bose S., Helly J., Jenkins A., Sawala T., Cytowski M., 2016, MNRAS, 457, 3492
  * Hernandez & Gilmore (1998) Hernandez X., Gilmore G., 1998, MNRAS, 297, 517
  * Hopkins (2013) Hopkins P. F., 2013, MNRAS, 428, 2840
  * Hughes et al. (2019) Hughes M. E., Pfeffer J., Martig M., Bastian N., Crain R. A., Kruijssen J. M. D., Reina-Campos M., 2019, MNRAS, 482, 2795
  * Hughes et al. (2020) Hughes M. E., Pfeffer J. L., Martig M., Reina-Campos M., Bastian N., Crain R. A., Kruijssen J. M. D., 2020, MNRAS, 491, 4012
  * Inoue (2009) Inoue S., 2009, MNRAS, 397, 709
  * Kaplinghat et al. (2016) Kaplinghat M., Tulin S., Yu H.-B., 2016, Phys. Rev. Lett., 116, 041302
  * Kaur & Sridhar (2018) Kaur K., Sridhar S., 2018, ApJ, 868, 134
  * Kruijssen (2012) Kruijssen J. M. D., 2012, MNRAS, 426, 3008
  * Kruijssen et al. (2011) Kruijssen J. M. D., Pelupessy F. I., Lamers H. J. G. L. M., Portegies Zwart S. F., Icke V., 2011, MNRAS, 414, 1339
  * Kruijssen et al. (2019a) Kruijssen J. M. D., Pfeffer J. L., Crain R. A., Bastian N., 2019a, MNRAS, 486, 3134
  * Kruijssen et al. (2019b) Kruijssen J. M. D., Pfeffer J. L., Reina-Campos M., Crain R. A., Bastian N., 2019b, MNRAS, 486, 3180
  * Kruijssen et al. (2020) Kruijssen J. M. D., et al., 2020, preprint,  (arXiv:2003.01119)
  * Law & Majewski (2010) Law D. R., Majewski S. R., 2010, ApJ, 714, 229
  * Mackey & Gilmore (2003) Mackey A. D., Gilmore G. F., 2003, MNRAS, 340, 175
  * Malhan et al. (2020) Malhan K., Valluri M., Freese K., 2020, preprint,  (arXiv:2005.12919)
  * McConnachie (2012) McConnachie A. W., 2012, AJ, 144, 4
  * Meadows et al. (2020) Meadows N., Navarro J. F., Santos-Santos I., Benítez-Llambay A., Frenk C., 2020, MNRAS, 491, 3336
  * Moore (1994) Moore B., 1994, Nature, 370, 629
  * Moster et al. (2013) Moster B. P., Naab T., White S. D. M., 2013, MNRAS, 428, 3121
  * Navarro et al. (1996a) Navarro J. F., Eke V. R., Frenk C. S., 1996a, MNRAS, 283, L72
  * Navarro et al. (1996b) Navarro J. F., Frenk C. S., White S. D. M., 1996b, ApJ, 462, 563
  * Navarro et al. (1997) Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493
  * Oñorbe et al. (2015) Oñorbe J., Boylan-Kolchin M., Bullock J. S., Hopkins P. F., Kereš D., Faucher-Giguère C.-A., Quataert E., Murray N., 2015, MNRAS, 454, 2092
  * Oman et al. (2019) Oman K. A., Marasco A., Navarro J. F., Frenk C. S., Schaye J., Benítez-Llambay A. r., 2019, MNRAS, 482, 821
  * Orkney et al. (2019) Orkney M. D. A., Read J. I., Petts J. A., Gieles M., 2019, preprint,  (arXiv:1906.04759)
  * Pascale et al. (2018) Pascale R., Posti L., Nipoti C., Binney J., 2018, MNRAS, 480, 927
  * Peñarrubia et al. (2008) Peñarrubia J., Navarro J. F., McConnachie A. W., 2008, ApJ, 673, 226
  * Petts et al. (2015) Petts J. A., Gualandris A., Read J. I., 2015, MNRAS, 454, 3778
  * Pfeffer et al. (2018) Pfeffer J., Kruijssen J. M. D., Crain R. A., Bastian N., 2018, MNRAS, 475, 4309
  * Pfeffer et al. (2020) Pfeffer J. L., Trujillo-Gomez S., Kruijssen J. M. D., Crain R. A., Hughes M. E., Reina-Campos M., Bastian N., 2020, MNRAS, 499, 4863
  * Planck Collaboration XVI (2014) Planck Collaboration XVI 2014, A&A, 571, A16
  * Pontzen & Governato (2012) Pontzen A., Governato F., 2012, MNRAS, 421, 3464
  * Read & Gilmore (2005) Read J. I., Gilmore G., 2005, MNRAS, 356, 107
  * Read et al. (2006) Read J. I., Goerdt T., Moore B., Pontzen A. P., Stadel J., Lake G., 2006, MNRAS, 373, 1451
  * Read et al. (2019) Read J. I., Walker M. G., Steger P., 2019, MNRAS, 484, 1401
  * Reina-Campos & Kruijssen (2017) Reina-Campos M., Kruijssen J. M. D., 2017, MNRAS, 469, 1282
  * Reina-Campos et al. (2019) Reina-Campos M., Kruijssen J. M. D., Pfeffer J. L., Bastian N., Crain R. A., 2019, MNRAS, 486, 5838
  * Reina-Campos et al. (2020) Reina-Campos M., Hughes M. E., Kruijssen J. M. D., Pfeffer J. L., Bastian N., Crain R. A., Koch A., Grebel E. K., 2020, MNRAS, 493, 3422
  * Renaud et al. (2011) Renaud F., Gieles M., Boily C. M., 2011, MNRAS, 418, 759
  * Rocha et al. (2013) Rocha M., Peter A. H. G., Bullock J. S., Kaplinghat M., Garrison-Kimmel S., Oñorbe J., Moustakas L. A., 2013, MNRAS, 430, 81
  * Sawala et al. (2015) Sawala T., et al., 2015, MNRAS, 448, 2941
  * Schaller et al. (2015) Schaller M., Dalla Vecchia C., Schaye J., Bower R. G., Theuns T., Crain R. A., Furlong M., McCarthy I. G., 2015, MNRAS, 454, 2277
  * Schaye et al. (2015) Schaye J., et al., 2015, MNRAS, 446, 521
  * Schechter (1976) Schechter P., 1976, ApJ, 203, 297
  * Schneider et al. (2017) Schneider A., Trujillo-Gomez S., Papastergis E., Reed D. S., Lake G., 2017, MNRAS, 470, 1542
  * Shao et al. (2013) Shao S., Gao L., Theuns T., Frenk C. S., 2013, MNRAS, 430, 2346
  * Shao et al. (2018) Shao S., Cautun M., Deason A. J., Frenk C. S., Theuns T., 2018, MNRAS, 479, 284
  * Spergel & Steinhardt (2000) Spergel D. N., Steinhardt P. J., 2000, Phys. Rev. Lett., 84, 3760
  * Springel (2005) Springel V., 2005, MNRAS, 364, 1105
  * Springel et al. (2001) Springel V., Yoshida N., White S. D. M., 2001, New A, 6, 79
  * Springel et al. (2005) Springel V., Di Matteo T., Hernquist L., 2005, MNRAS, 361, 776
  * Strigari et al. (2010) Strigari L. E., Frenk C. S., White S. D. M., 2010, MNRAS, 408, 2364
  * Strigari et al. (2017) Strigari L. E., Frenk C. S., White S. D. M., 2017, ApJ, 838, 123
  * Tremaine et al. (1975) Tremaine S. D., Ostriker J. P., Spitzer L. J., 1975, ApJ, 196, 407
  * Trujillo-Gomez et al. (2020) Trujillo-Gomez S., Kruijssen J. M. D., Reina-Campos M., Pfeffer J. L., Keller B. W., Crain R. A., Bastian N., Hughes M. E., 2020, preprint, p. arXiv:2005.02401 (arXiv:2005.02401)
  * Usher et al. (2018) Usher C., Pfeffer J., Bastian N., Kruijssen J. M. D., Crain R. A., Reina-Campos M., 2018, MNRAS, 480, 3279
  * Walker & Peñarrubia (2011) Walker M. G., Peñarrubia J., 2011, ApJ, 742, 20
  * Walker et al. (2009) Walker M. G., Mateo M., Olszewski E. W., Peñarrubia J., Evans N. W., Gilmore G., 2009, ApJ, 704, 1274
  * Wang et al. (2011) Wang J., et al., 2011, MNRAS, 413, 1373
  * Wang et al. (2017) Wang M. Y., et al., 2017, MNRAS, 468, 4887
  * Wang et al. (2018) Wang M.-Y., et al., 2018, preprint, p. arXiv:1809.07801 (arXiv:1809.07801)
  * Wang et al. (2019) Wang M. Y., et al., 2019, ApJ, 875, L13
  * Webb & Vesperini (2018) Webb J. J., Vesperini E., 2018, MNRAS, 479, 3708
  * Wiersma et al. (2009) Wiersma R. P. C., Schaye J., Smith B. D., 2009, MNRAS, 393, 99
  * Yozin & Bekki (2012) Yozin C., Bekki K., 2012, ApJ, 756, L18
  * de Boer & Fraser (2016) de Boer T. J. L., Fraser M., 2016, A&A, 590, A35
  * de Boer et al. (2012) de Boer T. J. L., et al., 2012, A&A, 544, A73

##  Appendix A The DM density profile of Fornax-mass dwarfs

Figure 11: Mean, radially averaged density profile, $\rho(r)$, for the field
Fornax-mass haloes in the simulation. The red line shows the best-fit NFW
profile.

In agreement with previous results from Eagle (Schaller et al., 2015), the
spherically averaged density profiles of our Fornax-mass haloes are well fit
by the NFW profile. This is shown in Fig. 11.

##  Appendix B The halo properties of GC rich dwarfs

Figure 12:  Comparison of the mass and concentration of DM haloes of two
Fornax-analogue populations. The black dotted line shows the distribution for
our reference population of all Fornax-mass dwarfs. The red solid line shows
the subset of systems that have 3 or more GCs. The vertical arrows indicate
the median of each PDF.

In Fig. 12 we show how the mass and concentration of the field Fornax
analogues varies when, in addition to selecting the analogues according to
stellar mass, we also require that they have a large number of GCs. To obtain
a reasonable sample of dwarfs with many GCs, we select the subsample that has
three or more GCs. The left-hand panel of Fig. 12 shows the PDF of $M_{\rm
200}$ for all the field dwarfs (dotted line) that were part of our initial,
stellar-mass only selection and field dwarfs that have three or more GCs
(solid line). The latter reside in slightly more massive haloes, with a
median, $M_{\rm 200}=3.4\times 10^{10}{~{}\rm M_{\odot}}$, that is 1.2 times
larger than the median of the full Fornax-analogue sample. This agrees with
the fact that more massive galaxies host more GCs (e.g. Harris et al., 2017).
The difference in the concentrations between the two samples is rather small,
as shown in the right-hand panel.

##  Appendix C The mass evolution of GCs in Fornax-mass dwarfs

Figure 13: The average mass evolution of our GC sample. The plot shows the
ratio of the mass at a given time to the initial GC mass as a function of the
time since the GC formed. The result for all the GCs is shown by the black
solid line with the shaded region showing the 68 percentile object-to-object
scatter. The various coloured lines show the results for GCs in different
initial mass ranges (see legend). Each of the GC subsamples contains roughly
half of the full sample.

In Section 4 we discussed the effect of GC mass on the decay of its orbit. To
make our predictions as realistic as possible, we accounted for the decrease
in GC mass as a function of time. The average trend, quantified as the ratio
of the GC mass at time $t$ and the initial mass, as a function of the GC age,
obtained from E-MOSAICS, is shown in Fig. 13. As described in Section 2,
E-MOSAICS accounts for the GC mass lost by stellar evolution, two-body
relaxation and tidal shocks. When calculating dynamical friction in post-
processing, we assumed that the GC mass follows the relation shown in the
figure (solid black line). The mass evolution is slightly affected by the
initial mass of the GCs but, for simplicity, we neglect this effect.