© E. Vitral and P. Boldrini 2022

1. Introduction

There is no doubt that the nature of dark matter (DM) is one of the most elusive concepts in modern-day physics. However, the existence of this astrophysical component has been used and requested to explain a vast range of phenomena for a considerable amount of time, including the following: Back when Zwicky (1933, 1937) proposed that some sort of non-luminous matter could compose the amount of mass needed to explain the discrepancy between mass measurements of the Coma cluster based on the Virial theorem (e.g. Binney & Tremaine 2011) and based on the brightness and number of galaxies, up to recent measurements from the Planck Collaboration XIII (2016), Planck Collaboration VI (2020) which yield very robust fits of the cosmic microwave background (CMB) using a Λ-cold dark matter (ΛCDM) model accounting for the existence of DM. For the period in between, other important confirmations of this mysterious dark component were provided: Rubin & Ford (1970) and Rubin et al. (1980) showed that the rotation curves of outer stars in nearby galaxies needed an extra amount of mass (compared to observable, luminous matter) to explain their high velocity values and finally, gravitational lensing studies (e.g. Taylor et al. 1998) have also confirmed that the total amount of mass in many galaxy clusters corresponds to the dynamical measurements accounting for DM.

Such findings point to this DM as a fundamental component in galaxy formation, which is present in most galaxies as an enveloping halo, from the smallest to the highest scales. Thus, in principle, it seems curious that dense collections of stars such as globular star clusters (GCs), that span from ∼10⁵ to 10⁷ M_⊙, with some of them thought to be accreted dwarf galaxies (e.g. Majewski et al. 2000; Bekki & Freeman 2003; Pechetti et al. 2022), do not seem to require any significant amount of DM to explain their dynamical mass (e.g. Shin et al. 2013; Conroy et al. 2011; Ibata et al. 2013; Moore 1996; Baumgardt et al. 2009; Lane et al. 2010; Feng et al. 2012; Hurst et al. 2015), although some recent works seem to suggest otherwise for a few particular sources (Carlberg & Grillmair 2022; Errani et al. 2022). In fact, Peebles (1984) proposed a formation scenario where GCs are formed inside their own DM mini-halo¹, and further studies defended that if formed before the end of re-ionisation², GCs could be smaller counterparts of galaxies (e.g. Bromm & Clarke 2002, Fig. 2 from Mamon et al. 2012; Silk & Mamon 2012 for a review on galaxy formation).

On the other hand, different formation scenarios, where GCs are not necessarily embedded in DM mini-halos also exist. For instance, GCs could be formed as bound gas clouds (Peebles & Dicke 1968), as galaxy DM-free fragments (e.g. Searle & Zinn 1978; Abadi et al. 2006), as relics of young massive clusters (YMCs, Portegies Zwart et al. 2010; Longmore et al. 2014) that formed in the high-redshift Universe (Kruijssen 2014, 2015), or simply that formed in situ along with its host galaxy (e.g. Harris 1991; Forbes et al. 1997). Moreover, recent cosmological simulations indicate realistic mechanisms through which all these scenarios (i.e. DM and DM-free) can actually happen (Trenti et al. 2015; Kimm et al. 2016; Ricotti et al. 2016; Keller et al. 2020; Ma et al. 2020), making it reasonable to argue that GCs likely originate from more than a single formation channel³.

Another reason that makes the understanding and further confirmation of different formation channels difficult is that much of the main consequences of these channels are better observable in the outskirts of GCs, where one often lacks good quality data. For example, the detection of tidal tails or stellar streams that might relate to an accretion event can be complicated by an observational bias where the stars in the stream are less luminous than the ones in the central and dense regions of the cluster, and thus more difficult to observe, comparatively (Balbinot & Gieles 2018). Similarly, a potential DM mini-halo could present a much more diffuse structure than the stellar component, so that its dynamical detectability might only be possible beyond several scale radii (Peñarrubia et al. 2017), where GC stars are usually confused with galactic field stars.

With the astrometric revolution brought by the Gaia mission (Gaia Collaboration 2018a,b, 2021), and the promising future discoveries of the James Webb Space Telescope (JWST, Gardner et al. 2006) and the Euclid mission (Laureijs et al. 2011; Lançon et al. 2021), the need to better constrain the expected differences between the multiple GC formation scenarios is a priority, so that these rich data sets can be fully exploited to better understand the many long-sought questions regarding GC formation. As a matter of fact, although many robust attempts to better model the observational implications of the DM mini-halo scenario have been made, the high computational cost of simulating a GC+DM system in a Milky Way (MW) type of galaxy forced these attempts to be placed in idealised scenarios: for instance, isolated GCs not experiencing tidal forces (Peñarrubia et al. 2017), or orbiting GCs in a static potential (Mashchenko & Sills 2005a). Preferably, simulations with clusters experiencing tidal forces, in a host galaxy fully composed of particles (and not just a static potential) would allow to take into account more correctly dynamical friction, tidal effects between globular clusters and their host galaxy, and dynamical shocks with larger structures of the host.

In this work, we aim to clearly separate the observational behaviours of GCs that are formed inside a DM mini-halo and which are devoid of one, both of them orbiting a host galaxy and thus permanently experiencing a tidal field. We do it by performing N-body simulations of a GC system with and without a DM embedding mini-halo, alongside with a host galaxy. In order to bypass the high computational cost mentioned above while still keeping a high resolution, we place our GCs in a dwarf spheroidal (dSph) galaxy, following the same prescriptions as the Fornax dSph, in a similar manner than done in Boldrini et al. (2020a) ⁴. Such setup allows one to consider much less stars than it would be needed if the satellites were orbiting a MW-like galaxy and keep a satisfactory resolution, and also to avoid using a potential instead of particles. Besides, by taking a real galaxy (i.e. Fornax) as our model, we are likely better exploring the dynamics and orbital evolution of different parameters, and thus reaching more realistic conclusions.

We describe our methods in Sect. 2, while our results and main DM signatures are presented in Sect. 3. We discuss the implications of our work and conclude it in Sects. 4 and 5. Throughout the rest of the paper, we reference the stellar component of the system of GC stars without DM as ⋆, and the stellar component of the system of GCs formed inside DM mini-halos as ⋆•. The dark matter is labelled as •.

2. Methods

2.1. Initial conditions

The initial conditions for the Fornax system were taken from Boldrini et al. (2020a), and although a complete description can be seen in their Sect. 2, we remind readers here of the most salient points, which are summarised in Table 1. We consider the scenario from that paper where the GCs were accreted recently, at z = 0.36 (i.e. 4 Gyr ago) by Fornax. It ensures that at z = 0, the GCs embedded in DM are still orbiting and no star clusters form in the centre of Fornax (in accordance with observations), and also that at 3 Gyr, the cluster’s positions relative to Fornax are consistent with their observed projected distances.

The Fornax dSph is modelled with a stellar component following the Plummer (Plummer 1911) profile, with total mass of 3.82 × 10⁷ M_⊙ (de Boer & Fraser 2016) and Plummer scale radius of 1.038 kpc (Strigari et al. 2010). Its DM halo follows a NFW (Navarro et al. 1996) profile with total mass 10⁹ M_⊙. Given the halo mass and redshift, the scale radius of the Fornax DM halo was estimated from cosmological N-body simulations (Prada et al. 2012).

The GCs on the other hand, started with a King (1962) stellar distribution, with total mass of 10⁶ M_⊙ and a King radius of 1 pc, lower than the observed radius, since it is susceptible to increase through dynamical processes such as mass loss (Mackey & Gilmore 2003). When assigning their DM mini-halo, we prescribed a mass of 2 × 10⁷ M_⊙, also following a NFW profile, and computed the respective scale radius in a similar fashion than for the dSph. The GC initial positions in phase space were set after analysing 7395 isolated dwarf galaxies with mass ∼10⁹ M_⊙ in the Illustris TNG-100 simulation (Pillepich et al. 2018), whose mass resolution is 7.5 × 10⁶ M_⊙, hence similar to the DM mini-halo mass. The positions were then selected by picking the cases where the sub-halos projected distances to the dwarf, at z = 0, were in the same range as the ones observed for the Fornax dSph GCs. From this subset, a new set of positions was obtained by going back 4 Gyr in time, and finally, the final values were taken by computing the maximum weights in each dimension.

Our choice of ratio for the mass of the GC mini-halo and stellar component (i.e. M_DM/M_stars = 20) is consistent with recent works targeting this scenario. For instance, Mashchenko & Sills (2005a) simulated GCs with a stellar mass of 8.8 × 10⁴ M_⊙ and a mini-halo virial mass of 10⁷ M_⊙, while Peñarrubia et al. (2017) tested models with mass ratios M_DM/M_stars = {0, 10, 100}. Indeed, GC systems with a DM mini-halo are very similar to orbiting (or accreted) DM sub-halos, which sit in a specific mass range for dwarf galaxies, predicted by the extended Press-Schechter (EPS) formalism (Bond et al. 1991; Lacey & Cole 1993). A Fornax-like galaxy is considered to have accreted roughly twelve DM sub-halos (see Fig. 1 from Boldrini et al. 2020b) with a mass ratio M_Fornax/M_sub-halo ∼ 10 − 10² (Zavala & Frenk 2019), consistent with our values.

2.2. Simulations

To generate our N-body objects, we use the initial condition code MAGI (Miki & Umemura 2018). Adopting a distribution-function-based method, it ensures that the final realisation of the galaxy is in dynamical equilibrium (Miki & Umemura 2018). We perform our simulations with the high performance collisionless N-body code GOTHIC (Miki & Umemura 2017). This gravitational octree code runs entirely on a Graphics Processing Unit (GPU) and is accelerated by the use of hierarchical time steps in which a group of particles has the same time step (Miki & Umemura 2017). In order for the simulated GCs to fully relax prior to introducing the Fornax tidal field, they were previously left to evolve in isolation for 2 Gyr.

We evolve the Fornax-GC system over 4 Gyr in each scenario, i.e. for GCs with and without a DM mini-halo. We set the particle resolution of all the live objects to 50 M_⊙⁵ and the gravitational softening length to 0.1 pc. Numerical convergence tests have been performed in Boldrini et al. (2020a). The use of a host composed of live particles rather than just a static potential provides a better handling of dynamical friction, which is the deceleration of massive particles (here, the GC) after interacting with less massive particles (here, the stellar and DM particles from the host). This is necessary to interpret the evolution of orbital radii of GCs, both with and without a DM mini-halo, and to robustly quantify their differences. Furthermore, the pruning of tidal structures (Binney & Tremaine 2011, Sect. 7.3) due to the shocking of such structures with the host field tracers is an effect completely neglected when using static potentials. Hence, the study of stellar stream’s properties and projected ellipticity are also better addressed when using live particles, such as we do.

2.3. Data analysis

Using the instantaneous orbital radius r_orbit and the instantaneous satellite mass M_sat, we calculate the theoretical tidal radius of each GC (or GC+DM system) at each snapshot from the simulation as derived by Bertin & Varri (2008), and used in many GC-related studies (e.g. Khalaj & Baumgardt 2016; Daniel et al. 2017; Webb et al. 2019). For the system composed of GC stars plus DM particles (i.e. GC+DM), we measure the total tidal radius in our simulations by taking into account the mass of both DM and stars in our calculation. In order to determine the bound particles of the DM and stellar components of GCs, we follow the procedure described in van den Bosch & Ogiya (2018, see their Sect. 2.3). The bound mass is determined by computing iteratively the binding energy of each particle. Initially, the method assumes that each particle is bound. We use the publicly available code galpy⁶ (Bovy 2015) to compute the quantities E and L_Z, for the integrals of motion space.

We performed most of the data analysis with internal routines of the astrometry modelling code BALROGO⁷ (Vitral 2021). This includes fits of the surface density (and ellipticity), construction of velocity anisotropy and dispersion profiles, and sky projection. Whenever targeting these themes, we mention in detail how we proceeded to perform the analysis.

3. Results: Dark matter signatures

In this section, we analyse the main implications of the presence of a DM mini-halo in the overall GC dynamics and morphology, according to our simulations. We remind that our analysis is modelled on the specific case of the Fornax dSph, in order to probe a realistic scenario.

3.1. Dark matter shield

In the beginning of our simulations with DM mini-halos, all GCs presented a DM envelope massive and concentrated enough, so their total tidal radius was considerably greater than the case without a DM mini-halo (see Fig. 1). Whenever we had such tidal radius increase, the DM particles are stripped beforehand GC stars, so the DM envelope works effectively as a shield against tidal stripping of the stellar component, being gradually removed as the system experiences stronger tidal forces from the host galaxy.

Fig. 1. Tidal radius growth due to an extra dark mass: Scheme illustrating the increase of tidal radius (blue) in the globular cluster system (orange) embedded in a dark matter mini-halo (grey), in the left, in comparison with the globular cluster system without dark matter, in the right, at a same orbital radius. This tidal radial expansion is mostly explained by its dependence on the satellite’s mass, which is greater for the case with a DM mini-halo.

In order to visualise this effect, we created a velocity dispersion map (see Appendix B.1 for details) of DM particles, which are much more spread than GC stars, and thus provide better spatial completeness to the map. Such map helps us to spot the transitory region where the DM particles start to effectively feel tidal effects of the host galaxy, which is characterised by a steep increase of the velocity dispersion produced by tidal heating of the system.

In Fig. 2, we show this velocity dispersion map for GC5, one of the clusters where this DM shield seemed most effective. We display the last six pericentres of its orbit, when the tidal effects are stronger. The extension of bound GC stars and bound DM particles are highlighted as dotted and dashed green lines, respectively, while the theoretical tidal radius, calculated according to Sect. 2.3, is displayed as a solid green circle. The centres of Fornax and of the GC are represented as a thick green cross and a plus sign, respectively.

Fig. 2. Dark matter shield: Velocity dispersion map of dark matter particles for GC5, projected in the X vs. Y plane and centred in the centre of mass of the globular cluster system. We display the last six pericentres of its orbit, where the tidal effects are stronger. The extension of bound globular cluster stars and bound dark matter particles are highlighted as dotted and dashed green lines, respectively, while the theoretical tidal radius, calculated according to Sect. 2.3, is displayed as a solid green circle. The maps are colour-coded logarithmically from blue (lower dispersion) to red (higher dispersion). The centres of Fornax and of the globular cluster are represented as a thick green cross and a plus sign, respectively. For this cluster, we notice that the empirical tidal radius, well traced by the blue region, remains always larger than the bound stars radii. This argues in favour of the dynamical presence of a dark matter shield, as illustrated in Fig. 1.

We also observe, still in Fig. 2, a gradual decrease in the extension of the blue region, characterised by a low velocity dispersion, where the DM shield is effective. This means that at first, the DM shield is highly protective, and with time, as DM particles are stripped from the cluster, the system mass decreases, and the shield becomes weaker. The red colours, on the other hand, point to a region of high dynamical heating, which becomes more intense in the centre of Fornax and favours tidal stripping of the DM shield.

This protective DM shield is a key mechanism to explain most of the DM signatures we observed in GCs embedded in DM. Although this phenomenon is ubiquitous in our GCs formed in DM mini-halos, the survival of such shield depends mostly on the orbital parameters of each cluster, and a specific discussion on this matter is addressed in Sect. 4.1. In the following, we focus on describing the main impacts of such shield with respect to the case where GCs are not embedded in DM mini-halos.

3.2. Orbital decay

The exchange of energy between the satellite and the field stars from its host galaxy will lead to the drag of the satellite, also referred to as dynamical friction. The use of live particles for both the satellite and its host, such as we do, allows to provide a better handling of this effect. The dynamical friction time of the cluster is the timescale needed for the satellite to reach the centre of mass of the host galaxy. It has been defined in Binney & Tremaine (2011, Eqs. (7)–(26)), and follows the relation t_fric ∝ 1/M, where M the satellite’s total mass.

Applying this relation to our simulations, where in one case we have a GC system alone of 10⁶ M_⊙ and in the other case where the total system mass is that of the GC plus the DM mini-halo (i.e. 2 × 10⁷ M_⊙), we find that the system with DM is supposed to sink to the centre twenty times faster than the system without the DM mini-halo. Indeed, when looking at the upper panels of Fig. 3, one can notice that the systems with DM (solid green) occupy much shorter orbital radii than the systems without DM (dashed red) throughout their orbits.

Fig. 3. Velocity anisotropy: Evolution of the orbital radius (i.e. rorbit, upper plots) and of the outer velocity anisotropy (i.e. β∞, lower plots), for the five simulated globular clusters. The clusters originally embedded in a dark matter mini-halo are depicted in solid green, and those without it are displayed in solid red. For β∞, we display the running mean over the five closest snapshots (time-wise), for better visualisation, which explains why the initial β∞ values seem to be different for the two scenarios. We observe typically lower rorbit and more radial (i.e. higher) β∞ for clusters with a dark matter mini-halo.

As a consequence, systems with a DM mini-halo tend to be located closer to their host centre sooner, and thus feel a stronger dynamical heating imposed by the host galaxy. This dynamical heating, on its turn, can potentially work to remove more GC stars (along with their DM envelope) than in the case without DM. Hence, our simulations help to answer whether the presence of a DM mini-halo is rather protective or disruptive with respect to the cluster’s stellar component, a point that we discuss in detail in Sect. 4.1.

3.3. Velocity anisotropy

One of the points not addressed in previous works targeting a DM mini-halo in GCs (e.g. Mashchenko & Sills 2005a; Peñarrubia et al. 2017) is the evolution of the velocity anisotropy of the stellar component, in both cases with and without a DM mini-halo. The velocity anisotropy (‘anisotropy’ for short) is defined as in (Binney 1980):

$$\begin{array}{c}\hfill {\displaystyle \beta (r)=1-\frac{{\sigma }_{\theta }^2(r)+{\sigma }_{\varphi }^2(r)}{2{\sigma }_{\mathrm{r}}^2(r)},}\end{array}$$(1)

where θ and ϕ are the tangential components of the coordinate system, while ${\sigma }_i^2$ stands for the velocity dispersion of the component i of the coordinate system. In spherical symmetry, σ_ϕ = σ_θ.

We fit the anisotropy for each cluster and snapshot of our simulations in a Bayesian frame detailed in Appendix A. For that, we consider the generalised parametrisation from Osipkov (1979), Merritt (1985), which has been applied in many works modelling the dynamics of spherical stellar systems (e.g. Walker & Peñarrubia 2011; Mamon et al. 2013; Read et al. 2021; Vitral & Mamon 2021; Vitral et al. 2022)

$$\begin{array}{c}\hfill {\displaystyle {\beta }_{\mathrm{gOM}}(r)={\beta }_0+({\beta }_{\infty }-{\beta }_0)\frac{r^2}{r^2+r_{\beta }^2},}\end{array}$$(2)

where r_β is the anisotropy radius, and (β₀, β_∞) are the anisotropy values at the centre of the system, and at infinity, respectively. We keep in mind that the initial β(r) in our simulations is not necessarily isotropic as a result of the prior 2 Gyr relaxation run we mention in Sect. 2.2, which can change the anisotropy shape, specially in the cluster outskirts (β₀ still lies closer to zero).

As we are most interested in the differences pertaining to the external dynamics of the GC systems, we chose to plot in Fig. 3 (lower plots) the evolution⁸ of β_∞ for the GCs with a DM mini-halo (solid green) and without it (solid red). This allows us to quantitatively differentiate the impact of the Fornax tidal field on the cluster outer regions. Indeed, the external tidal field has been shown to play a major role on the outer anisotropy of GCs, by removing stars in radial orbits in the cluster’s outskirts, which are more easily pruned by tides (e.g. Takahashi 1995; Tiongco et al. 2016a; Zocchi et al. 2016; Bianchini et al. 2017). As a result, one should expect more tangential (or equivalently, less radial) anisotropy profiles in systems that undergo more severe tidal interactions.

The variations of β_∞ are more intense for clusters without the DM mini-halo, since they cover a broader range of orbital radii. Furthermore, an important diagnosis is that the outer anisotropy of clusters with the DM mini-halo is, in general, more radial, even though those clusters occupy regions closer to the host galaxy. We thus associate this tidal resilience with the dynamical presence of the DM shield discussed in Sect. 3.1.

We also notice that the initial value of β_∞ (not necessarily zero) does not really impact these conclusions. This is because the values of β_∞ not only depend considerably on the orbital radius, but their variation speed is almost identical to the respective orbital radius variation. This means that regardless of the initial anisotropy conditions, the values of β_∞ adapt rapidly to the ongoing tidal field and its further interpretation is more straightforward.

Hence, the outer anisotropy is an important parameter that allows to differentiate both scenarios. While clusters with more radial orbits (i.e. β_∞ ≳ 0.4) could belong to both scenarios, clusters with tangential anisotropy, or less radial orbits (i.e. β_∞ ≲ 0.4) were mostly obtained in our simulations without a DM mini-halo. This does not necessarily mean that clusters with very tangential orbits did not have a DM mini-halo in the past, but rather that searches for a current mini-halo surrounding them would likely be in vain.

3.4. Velocity dispersion profile

Peñarrubia et al. (2017) simulated isolated GCs embedded in DM mini-halos and showed that due to the extra DM mass, an inflation of the radial velocity dispersion profile towards outer radii⁹ is to be expected. In the current work, we are able to test such predictions for a more realistic scenario, where the GC+DM subsystem is orbiting a host galaxy, and therefore experiencing tidal effects. For that, we measured the radial velocity dispersion profiles (see Appendix B.2 for details) of the bound stars in our clusters in both the cases with and without a DM mini-halo, through the course of their evolution in the Fornax tidal field.

In Fig. 4 (right plots), we show the values of the maximum radial velocity dispersion as a function of time, colour-coded by the distance to the centre of Fornax (i.e. r_orbit), in kpc. One of the most important realisations of our analysis is the fact that the overall shape of the velocity dispersion is much more impacted by the host galaxy’s tidal field than by an eventual DM mini-halo: The amplitude of the dispersion, traced by its value at the peak follows a periodic variation, with same period than the GC orbit, and has values almost uniquely dependent on the ongoing tidal forces.

Fig. 4. Stellar velocity dispersion: Series of plots for GC5. The upper plots relate to the simulations devoid of dark matter mini-halos (⋆), while the lower plots indicate the results for the globular clusters formed in such mini-halos (⋆•). The two columns on the left display hand-picked snapshots where the radial velocity dispersion profile (i.e, σr(r)) resembled better to an isolated case without a dark matter mini-halo (left), and with a massive dark matter mini-halo (right), according to Fig. 2 from Peñarrubia et al. (2017). The column on the right presents the evolution of the maximum value of the radial velocity dispersion profile (i.e, σrmax), for each scenario concerning the dark matter mini-halo, colour-coded according to the distance of the cluster to the centre of the host galaxy (i.e. rorbit), with two vertical dashed lines corresponding to the instants from the two columns on the left. These plots argue that the tidal field from the host galaxy tends to have a much greater impact on inflating the velocity dispersion than the presence of a dark matter mini-halo. In fact, such mini-halos help to protect the cluster from tidal effects, rather than contributing to it.

In the plot, we can clearly observe that at pericentres (blueish), the velocity dispersion inflates as a whole: the tidal heating from the host galaxy is effectively felt more intensively, leading to a higher velocity dispersion. In contrast, at apocentres (reddish or blackish), the cluster is closer from an ideal isolated scenario, and tidal heating is less effective, leading to low velocity dispersion peaks.

To illustrate the much stronger dependence on tidal forces than on possible DM mini-halos, we selected two snapshots for GC5 for the case with and without a DM mini-halo (lower and upper plots, respectively). In these snapshots (Fig. 4, four left plots), we can verify that the radial velocity dispersion profile of both the DM and DM-free scenarios assume forms similar to both the isolated cases with no DM mini-halo (left), and with a massive DM mini-halo (right), as presented in Fig. 2 of Peñarrubia et al. (2017).

As a general trend, all the clusters had an increasing velocity dispersion close to pericentre, with multiple points of velocity dispersion inflation throughout the GC radial extension. In apocentres, as mentioned above, the clusters resembled better to an isolated case (with one or two inflation points), specially for the case with DM mini-halo, where the shapes retrieved by Peñarrubia et al. (2017) could be better spotted. The reason behind the best resemblance in apocentre for the case with DM mini-halos is directly connected to the DM shield, which protects the GC stars and approaches better the ideal isolated framework.

3.5. Rotation and ellipticity (flattening)

Satellites that are more affected by tidal effects should exchange more internal angular momentum from energy transfer with their host, and this will directly impact the degrees of internal rotation of the analysed system (e.g. Tiongco et al. 2016b, 2018). Hence, to quantify the differences pertaining to rotation between GCs formed in DM mini-halos and those who were not, we plot in the upper panel of Fig. 5, the evolution of their mean line-of-sight (LOS) velocity dispersion¹⁰ (with 1-σ Poisson error bars), since this quantity is used in many observational studies to constrain the degree of rotation of GCs (e.g. Bianchini et al. 2013; Boberg et al. 2017; Lanzoni et al. 2018). In order to recreate a two-dimensional projected distribution of stars, such as we observe in true data, and extract the respective LOS velocity, we placed, for each snapshot, the GC centre in the position (α,  δ)=(0,  0), at five kpc away from the observer¹¹. The dispersion in LOS velocity of all bound GC stars is neatly higher for clusters without a DM mini-halo, meaning that they suffered more prominent tidal interactions. In other words, GCs embedded in dark matter mini-halos exchange less internal angular momentum from tidal interactions, thus presenting smaller degrees of rotation.

Fig. 5. Rotation: On the top, we show the mean evolution of the line-of-sight (LOS) velocity dispersion of bound globular cluster particles, for each simulated cluster, while the bottom plot displays the mean ellipticity (or flattening) of the stellar distribution. We display the running mean over the ten closest snapshots (time-wise), for better visualisation. Error bars are represent 1-σ Poisson uncertainties. Clusters originally embedded in a dark matter mini-halo (i.e. ⋆•) are displayed in green, while those devoid of the mini-halo (i.e. ⋆) are represented in red. As a general trend, globular clusters embedded in dark matter mini-halos exchange less internal angular momentum from tidal interactions, thus presenting smaller degrees of rotation, which in turn is associated with smaller ellipticities. This again attests the efficacy of the dark matter shield illustrated in Fig. 2.

An increased amount of rotation can also contribute to larger values of ellipticities (or flattening) for the stellar distribution from GCs (Fabricius et al. 2014; Kamann et al. 2018). We thus fitted the distribution of bound stars with a Sérsic (Sérsic 1963; Sersic 1968) asymmetric model, using the recipe described in Appendix C. The fits yielded a semi-major (a) and semi-minor axis (b), from which the ellipticity (or flattening) could be calculated as

$$\begin{array}{c}\hfill e=\sqrt{1-{\left(\frac{b}{a}\right)}^2}·\end{array}$$(3)

The lower panel from Fig. 5 was then constructed by computing, at each instant, the mean ellipticity¹² of the five simulated clusters with and without DM. The error bars were calculated as $\sqrt{⟨ϵ⟩{}^2+{\sigma }_{\mathrm{e}}^2/n}$, where ⟨ϵ⟩ is the mean uncertainty of the ellipticity fit, σ_e is the standard deviation of the values from the n = 5 clusters.

First, we observe a decreasing pattern in the mean ellipticities of both scenarios, which is expected from orbiting satellites, due to the stripping of outer stars by the tidal field (Akiyama 1991)¹³. Nonetheless, one notices that clusters embedded in a DM mini-halo show considerably smaller ellipticities (i.e. e ≲ 0.35) than the clusters devoid of DM, which presented rather e ≳ 0.35. Such relatively small ellipticities are indeed predictable for systems less affected by tidal forces (van den Bergh 2008), attesting the efficacy of the DM shield illustrated in Fig. 2.

3.6. Tidal tails

The study of tidal tails in GCs has been revolutionised by Gaia data and simulations. For instance, a troubling question that arises when simulating GCs on tidal fields is to understand why these simulations usually predict much more prominent tidal tails (e.g. Boldrini & Vitral 2021; Montuori et al. 2007) than what is observed for the majority of MW clusters. The answer to this timing problem is partially given by Balbinot & Gieles (2018), who showed with simulations that there is a preferential bias towards the escaping of low-mass stars, specially in denser clusters. Such trend reduces considerably the visibility of the tails. In addition, Gieles et al. (2021) recently defended this trend by showing that the visible and extended tails of Palomar 5 are well explained by a supra-massive population of stellar-mass black holes, which is a characteristic associated with less dense GCs (Kremer et al. 2020). On the other hand, the presence of a DM mini-halo could also reduce the prominence of tidal tails, since the mini-halo is expected to be stripped beforehand GC stars, thus delaying tail formation (e.g. Bromm & Clarke 2002; Mashchenko & Sills 2005a; Saitoh et al. 2006; Bekki & Yong 2012; Boldrini & Vitral 2021).

In any case, the increase of better quality data such as Gaia EDR3 has allowed to go deeper into this question: Although some GCs, such as NGC 1851 and NGC 7089 (M2), were thought not to have tails based on ground-based imaging (Kuzma et al. 2016, 2018), further Gaia studies revealed long tails associated to them (Ibata et al. 2021). However, many clusters with no tidal features, or with only extended envelopes (without tails, see Piatti & Carballo-Bello 2020) are still present¹⁴, and there is still no consensus on which mechanisms are behind the lack of very extended tails in these GCs. On a similar note, Martin et al. (2022) recently used Gaia EDR3 and CFHT data to constrain a stellar stream, C-19, whose metallicity is consistent with it being a remnant of the oldest GC known in our Galaxy, even though the stream still presents a coherent structure to this day, thus requiring some sort of shielding mechanism to protect it from tidal stripping.

We decide to explore the consequences of a DM mini-halo on tidal tails, and we illustrate it in Fig. 6, with the stellar distribution of GC5 (one of the cases where the DM shield was most effective), for the case with (right) and without (left) an embedding DM mini-halo. The chosen snapshot was such that GC5 was located at roughly the same orbital radius for both scenarios, for a fairer comparison. We also provide online a video of the evolution of the tails in the two kind of clusters, with a respective zoomed-in version too. We observe much more prominent and obvious tails in the case without DM, with thick and well-defined streams measuring up to ≳10 kpc long. However, in the case where the GCs are formed inside DM mini-halos, the stellar distribution remains roughly spherical. We connect this difference between two kind of GCs to the protection of the DM shield (as explained in Sect. 3.1, and in Figs. 1, 2), which reduces and delays tidal effects on the GCs where it is present.

Fig. 6. Dark matter influence on tidal tails: Stellar distribution (and respective zoom in) of GC5 for the case where it was originally embedded in a dark matter mini-halo (right) and where it was devoid of it (left), both centred in the GC centre of mass. The chosen snapshot was such that GC5 was located at roughly the same orbital radius for both scenarios, for a fairer comparison. This plot highlights more prominent tails in the dark matter-free case, while clusters formed inside dark matter mini-halos present a more compact stellar envelope. Hence, the dark matter shield (Sect. 3.1) has the effect of delaying the formation of tidal tails. Movies are available online.

Hence, although GCs embedded in DM mini-halos can still form stellar streams eventually, the ones devoid of DM seem to develop much longer, thicker and well-defined tails. Evidently, once the DM shield is destroyed, which can happen much before a Hubble time for clusters having specific orbital parameters (see Sect. 4.1), the development of tidal tails can be considered similar to the case without DM, and thus even clusters originally embedded in DM could present extended tails by present time, in those conditions. Having said that, the DM shield has the effect of delaying the formation of tidal tails, by making them much milder while the shield is present, and thus providing a possible answer to the timing problem mentioned in the beginning of this section.

3.7. Integrals of motion

By analysing the orbital evolution of simulated satellites in a static MW potential, Helmi & de Zeeuw (2000) proposed that the space of integrals of motion (IOM) was an optimal tool to recover the orbital history of accreted and disrupted satellites. In particular, by analysing the diagram of energy versus angular momentum in the Z galactic direction, they showed that GCs originating from an accreted satellite should remain in a similar position than its respective progenitor in IOM space, regardless of their phase and velocities being no longer the same. Even though energy is not conserved throughout orbital evolution of the satellite, the authors argued that those IOM positions should not move significantly from their original positions. Indeed, by using such method, Massari et al. (2019) related the origin of the galactic GC population to different past major accretions events from the MW.

We decided to analyse the impact of a DM mini-halo in the evolution of a GC in IOM space by plotting their respective position in the energy-angular momentum (i.e. E × L_Z) diagram for our five analysed clusters, at each snapshot. Figure 7 displays this evolution for both scenarios: with a DM mini-halo in green, and without it in red, both beginning in the same position, marked by a light grey colour, and ending at their respective darker crosses.

Fig. 7. Integrals of motion: Evolution of our simulated clusters in the integrals of motion (IOM) space, in particular, the E × LZ diagram. Their positions in this diagram are scattered in a sequential colour-map, starting from the same point in light grey at t = 0, and ending in darker tones (green for globular clusters originally formed in dark matter mini-halos, and red for clusters devoid of dark matter). The last snapshot is marked by a cross, for each cluster. This plot highlights a bimodal evolutionary distribution of clusters originally embedded in dark matter (moving towards lower energies) and those who were not (moving only slightly towards higher energies), and more importantly, it shows that clusters originally embedded in dark matter move significantly in IOM space, such that their association with past merger events through this diagram is not reliable.

One observes that clusters devoid of DM (red) do not significantly move in the diagram throughout their evolution, in agreement with Helmi & de Zeeuw (2000), which in turn serves as a validity check for our simulations. However, the same cannot be said for GCs originally embedded in DM, which move considerably during their orbital evolution. This is mainly due to the massive loss of energy from dynamical friction (see Sect. 3.2), and further changes on the host potential, which are better handled by our use of an evolving host composed of particles, instead of static potentials. Furthermore, this energy loss shifts the GC population formed in DM mini-halos to lower regions of the E × L_Z diagram, creating a clear bimodal evolutionary distribution in Fig. 7, where green evolution tracks move downwards (lower energies), while the red ones tend to move only slightly upwards (higher energies). This does not imply though, that at a given time, a cluster with low energy necessarily had a DM mini-halo, since a DM-free GC could naturally form in such regions and not move significantly thereafter.

Given that there is no reason why these conclusions should change qualitatively when considering the case of the MW, Fig. 7 has important implications. Not only it shines a light on the bimodal evolutionary distribution commented above, which can be used to guide the search of GCs candidates that have been formed in a DM mini-halo, but most of all, it shows that GCs formed in DM mini-halos move significantly in IOM space due to their more intense energy loss from dynamical friction. Hence, the association of such clusters with past merger events by using the IOM method is not reliable. Indeed, Boldrini & Bovy (2022) recently showed, by performing orbital integrations, that there was no clear association between many MW GCs and other MW satellites, even though a large fraction of them are believed to be accreted from previous IOM analyses.

4. Discussion

In the previous section, we have consistently showed that, at least during the first orbits of a GC embedded in a DM mini-halo, the DM behaves as a protecting shield that prevents the development of prominent tidal features. Now, we discuss a different aspect, which is under which circumstances this DM shield is effective.

4.1. Survival of the dark matter mini-halo

The DM mini-halo of our simulated GCs was shown to effectively work as a shield against tidal effects. However, not only this mini-halo is gradually removed by being itself stripped beforehand GC stars, as the increase of mass in the GC+DM system (with respect to the GCs devoid of DM) also brings the system closer to the centre of the host galaxy¹⁵, where tidal effects are stronger. As a result, the DM mini-halo is expected to become negligible with time, so that the question of interest becomes: When is a DM mini-halo effectively disrupted?

The answering of this question starts when comparing the five GCs we simulated, and their different orbital parameters. First, we can separate these five GCs into two categories: (i) The GCs where the DM shield is not effective by the end of our simulation and (ii) GCs whose shield still manages to protect them against strong tidal effects. Looking at Fig. 8, we can clearly assign GCs 1 and 2 to category (i), while GCs 3, 4 and 5 are better suited to category (ii), with GC 4 approaching to category (i). This assignment is due to the fact that the limit of bound stars in GCs 1 and 2 is comparable or greater than the respective tidal radius of the GC+DM system, meaning that the DM shield illustrated in Fig. 1 is no longer well observed. We also managed to observe that such disrupted mini-halos presented r_t >  r_1/2, • at some point of their orbital evolution (Fig. D.1), where r_t is the tidal radius of the GC+DM system, and r_1/2, • the half-mass radius of the bound DM particles. This criterion is similar to the relation presented in Hayashi et al. (2003) for the disruption of DM halos.

Fig. 8. Survival of the dark matter mini-halo: Similarly to Fig. 2, we show the velocity dispersion map of dark matter particles for the five globular clusters we simulated, projected in the X vs. Y plane (with exception of GC4, projected on the Z vs. Y plane for better visualisation of its orbit) and centred in the centre of mass of the globular cluster system. We display the last pericentre passages of their respective orbits, where the tidal effects are stronger. One can notice that while GCs 3, 4 and 5 manage to keep the dark matter mini-halo (and thus the dark matter shield effect depicted in Fig. 2), GCs 1 and 2 do not. The extent of the red region for GCs 1 and 2 is related to their different orbital parameters. For GCs 3, 4 and 5, we again notice that the empirical tidal radius is well traced by the internal blue region.

We argue that the DM loss in GCs 1 and 2 is accelerated by the many passages in pericentre (i.e. relatively small orbital periods), as well as the relatively small orbital radii (r_orbit) throughout the clusters’ orbits. Those factors force these GCs to stay longer next to the centre of the host galaxy, where the tidal interactions are more severe. As a consequence, the DM mini-halos in GCs 1 and 2 suffer from stronger tidal effects and are stripped faster. In Table 2, we display the mean of the structural parameters of stars and DM particles, over the last ten snapshots of the simulation, and by analysing the difference between Cols. 5 and 8 (i.e. the half-mass radius of the DM mini-halo, and the tidal radius of the respective system), one can infer that the DM mini-halo survives better to worse, in the following order: GC5 and GC3 similarly, followed by GC4, GC1 and GC2, respectively.

4.2. Detectability of dark matter

4.2.1. On the amount of detectable dark matter

If some GCs manage to preserve part of their original DM mini-halo, it is important that we understand the observational limitations that may allow us (or not) to detect this DM component. For instance, as the mini-halo mass distribution is more diffuse than the GC stellar component, one might ask if the DM amount in the GC inner regions is significant. For that purpose, we plot in Fig. 9 the mass profile of different tracers (i.e. DM and GCs, for both scenarios) as a function of distance from the cluster’s centre, at the last snapshot of our simulations.

Fig. 9. Cumulative mass profiles: We display the cumulative mass profiles at the last snapshot of our simulations, for the dark matter particles (black), their embedded stellar component (green) and the stellar component from the simulations without dark matter (dashed red). The blue line depicts the MCMC fit of a Zhao (1996)αβγ model to the dark matter mini-halo (see Sect. 4.2.2), with a fainter blue region encompassing the [2.5, 97.5] percentiles. The dotted lines represent local slopes of the stellar mass distribution at the respective regions shown in the plot.

Figure 9 reveals that up to ∼10 pc¹⁶, the mass ratio between the DM component (solid black) and the stellar component of the GCs (solid green) is of the order of 1%, meaning that the influence of DM in the internal dynamics of these clusters is mostly negligible¹⁷. In the past, many studies (Shin et al. 2013; Conroy et al. 2011; Ibata et al. 2013; Moore 1996; Baumgardt et al. 2009; Lane et al. 2010; Feng et al. 2012; Hurst et al. 2015) measured very low mass-to-light ratios in GCs, and sometimes used that as an argument to rule out the presence of dark matter in them. We show, in agreement with what has been previously pointed out in Bromm & Clarke (2002), Mashchenko & Sills (2005a), Saitoh et al. (2006), Bekki & Yong (2012), that GCs originally embedded in DM lose much of their initial DM content due to tidal interactions. This naturally leads to mass ratios (stellar to total mass) very close to unity, up to distances as far as available data usually goes.

Figure 9 also compares the mass profile of the GCs formed in DM mini-halos and those devoid of it (dashed red). Clearly, clusters with DM have an inner mass much higher than GCs devoid of DM, with GCs 3 and 5, which managed to retain most of their initial DM mini-halo, having inner masses up to ten times higher than their DM-free counterparts. Hence, in addition to denser clusters producing an observational bias where tidal structures are more difficultly detected (Balbinot & Gieles 2018), we point that the ones that might be formed in DM mini-halos also tend to have at the same time higher inner densities and less prominent tidal structures (see Sect. 3.6), which can render the search for stellar streams in these GCs particularly hard.

Finally, it is worth mentioning that the outer slopes of the stellar mass profile (dotted lines in Fig. 9) do not seem to show a necessary correlation with the analysed scenario (with or without DM), differently than what is shown in Peñarrubia et al. (2017, Fig. 2). We associate it once again, as in Sect. 3.4 for the velocity dispersion profile, with the presence of an intense tidal field (neglected in the former work) that shapes the outer regions of the stellar distribution in a more impactful way than an eventual DM mini-halo.

4.2.2. Analytical description

One of the main ways to look for a DM mini-halo in GCs is to perform mass modelling of observed data, which can be done by solving the Jeans equation (Binney 1980) for spherical systems with no streaming motions

$$\begin{array}{c}\hfill \frac{\mathrm{d}\left(\rho {\sigma }_r^2\right)}{\mathrm{d}r}+2\frac{\beta (r)}{r}\rho (r){\sigma }_r^2(r)=-\rho (r)\frac{G{M}_{\mathrm{T}}(r)}{r^2},\end{array}$$(4)

where we assume a total mass profile M_T(r), a anisotropy profile β(r), and a previously determined mass density profile ρ(r) for the kinematic tracers (here, stars). The term $\rho {\sigma }_r^2$ is the dynamical pressure that counteracts gravity. Since the observable tracers are stars, we only need a parametrisation for the DM mass profile in order to model it. This parametric form is then incorporated in the total mass profile M_T(r)=M_⋆•(r)+M_•(r). To obtain a suitable mass profile for the DM component, we tested two models, namely the Kazantzidis et al. (2004) model¹⁸,

$$\begin{array}{c}\hfill {\rho }_{\bullet }(r)={\rho }_0{r}^{-\gamma }{e}^{-r/a},\end{array}$$(5)

and the Zhao (1996)αβγ model,

$$\begin{array}{c}\hfill {\displaystyle {\rho }_{\bullet }(r)={\rho }_0{r}^{-\gamma }{[1+{\left(\frac{r}{a}\right)}^{\alpha }]}^{(\gamma -\beta )/\alpha },}\end{array}$$(6)

where ρ₀ is a scaling factor, a the scale radius, γ ∈ (0, 3) an inner slope, β an outer slope, which we take to be 5, such as in the Plummer profile, and α a scaling exponent, taken as 2, according to the Plummer profile as well¹⁹ Then, we converted those densities to mass profiles by performing $M_{\bullet }(r)\equiv {\displaystyle {\int }_0^r}4\pi x^2{\rho }_{\bullet }(x)\mathrm{d}x$ with MATHEMATICA12.

Similarly to what is shown in Mashchenko & Sills (2005a, Fig. 7), we noticed that the stripped DM particles from the original mini-halo had a tendency to form a second shell surrounding the remaining mini-halo. Such effect is only noticeable in the DM mass profile beyond a few hundred pc (see Fig. 9), but since we choose to analyse the mass profiles up to the outer regions of the clusters, we took that into account by assigning M_•(r)=M_•, 1(r)+M_•, 2(r), where M_•, 1(r) is an inner component, most susceptible to be detected or to influence the dynamics of the embedded star cluster, and M_•, 2(r) is an outer shell, modelled only for the sake of better constraining the outer slope of M_•, 1(r). This yields, for the Kazantzidis et al. (2004) profile,

$$\begin{array}{c}\hfill {\displaystyle M_{\bullet ,i}(r)=4\pi a^{3-\gamma }{\rho }_0\mathrm{\Gamma }(3-\gamma ,\frac{r}{a})}\end{array}$$(7)

and for the Zhao (1996)αβγ profile,

$$\begin{array}{c}\hfill {\displaystyle M_{\bullet ,i}(r)=\frac{4\pi a^{3-\gamma }{\rho }_0}{3-\gamma }{[1+{\left(\frac{a}{r}\right)}^2]}^{(\gamma -3)/2}}\end{array}$$(8)

where we note $\mathrm{\Gamma }(b,x)={\displaystyle {\int }_0^x}{t}^{b-1}{e}^{-t}\mathrm{d}t$, which is the lower incomplete gamma function, and i = {1, 2}.

Since we fix the αβ parameters when using the Zhao (1996)αβγ profile, we end up with the same amount of free parameters than the Kazantzidis et al. (2004) model, which allows us to compare the fits from the two models by simply choosing the one with highest likelihood as our preferred model. Both models performed well, and the preference of likelihoods depended not only on the simulated cluster, but on the considered radial range. In Fig. 9, we display in blue the fits of the Zhao (1996)αβγ model, which had an overall better performance for radii smaller than 1 kpc. Indeed, going beyond such distances is not of practical interest when modelling observed data, since at such lengths, the confusion with MW interlopers is considerably damaging for the data set. Hence, for data sets going up to a few tens of pc, Jeans modellers could simply use the Zhao (1996)αβγ model, with fixed α = 2 and β = 5, and a single shell (i.e. M_• = M_•, 1) to probe the existence of a DM mini-halo in GCs. We provide in Table D.1, the best fit parameters of our simulated clusters (at z = 0, for M_•, 1) in each mass model.

4.3. Comparison with the literature

We now compare our study with previous attempts to model GCs in DM mini-halos. Beginning with Mashchenko & Sills (2005a), the authors simulated GCs inside a tidal field, using a M_DM/M_stars ratio of nearly a hundred. Their simulations used static potentials for the host, which tend to underestimate the impacts of dynamical friction and the trimming of outer structures due to passages in high density regions. Their work showed nevertheless that GCs embedded in DM mini-halos are more resilient to severe tidal stripping, in contrast to DM-free GCs, and also that higher mass to light ratios are found in in the clusters outskirts, while the inner regions keep mass to light ratios similar to the purely baryonic case. They finally proposed that presence of obvious stellar streams is probably the only observational evidence that can reliably rule out a presence of significant amounts of DM in GCs.

Our conclusions are in general very similar to those from Mashchenko & Sills (2005a), with some important caveats. Indeed, we find that the presence of a DM mini-halo tends to initially protect the GC from tidal effects, but our more suitable handling of dynamical friction from a time-changing potential with live particles suggests that clusters with a DM mini-halo sink faster to the clusters centre, where the tidal field is more intense. As a result, clusters whose orbital parameters tend to naturally place them closer to the host will have this effect strengthened, resulting in a faster disruption of their mini-halo, and further cluster dissolution. Thus, the question if GCs embedded in DM are always more resilient to tidal effects is not as simple, and some of them might actually be disrupted faster, if originally orbiting near the host’s centre (see Sect. 4.1).

As in Mashchenko & Sills (2005a), we also find that the DM mass profile is negligible, with respect to the stellar component, up to at least roughly a half-mass radius from the clusters centre (Fig. 9). However, we do not associate obvious stellar streams as the only indicator of a DM-free GC. In fact, it is important to mention that although systems with DM do not form obvious stellar streams, they can still form milder and thinner ones before the DM mini-halo is entirely disrupted. Another important DM-free indicators could be a very tangential velocity anisotropy in the outskirts and much higher degrees of rotation, which indirectly leads to higher ellipticities as well. These indicators suggest that the GC stellar component is highly impacted by the host’s tidal field, hinting that a DM mini-halo is likely not present at current time²⁰. We also propose to use IOM space locations to evaluate if the existence of a DM mini-halo is likely (Fig. 7), since that clusters formed in this scenario tend move towards regions of lower energy in the E × L_Z diagram, due to dynamical friction, while still keeping in mind that DM-free clusters can eventually form in these regions and remain there during their evolution. Our new analysis of the IOM evolution additionally shows that associations of GCs and accreted satellites through nearby locations in IOM space are not reliable for clusters that were embedded in a DM mini-halo.

Recently, Peñarrubia et al. (2017) analysed the effects of isolated GCs embedded in DM mini-halos, with DM to stellar mass ratios of ten and a hundred. Their analyses predicted a flattening or increase of the velocity dispersion profile at outer radii for GCs with DM, along with shallower surface density profiles. They also provided analytical forms for the DM mass profile, which can be used in Bayesian approaches. Although such predictions seem indeed valid for isolated systems, when placing them inside a realistic tidal field, such signatures are mostly erased in the cluster outskirts. We have showed in Sect. 3.4 that an eventual inflation in the velocity dispersion profile can also happen in DM-free clusters due to tidal heating, sometimes in even more pronounced proportions (see Fig. 4). When at apocentres though, where tidal interactions are milder, we retrieved similar velocity dispersion shapes than Peñarrubia et al. (2017) for GCs embedded in DM mini-halos, since those are better shielded from tides and resemble better the isolated case. We did not find a clear relation between the outer mass slopes of GC with and without DM, likely due to the same arguments above.

As in Peñarrubia et al. (2017), we provide an analytical shape for the mass profiles of DM mini-halos, that can be used in Bayesian analyses, specially in Jeans modelling routines (Read et al. 2021). In our case though, our mass profiles were derived also for partially disrupted mini-halos, such as it is expected from systems evolving in a tidal field. We thus consider our parametric forms better suited to the realistic cases of GCs orbiting the MW or other Local Group galaxies.

Another recent simulation of GCs pertaining to the DM formation scenario was done by Carlberg & Keating (2022). However, the authors do not consider a single mini-halo per cluster such as in our work (or such as in Mashchenko & Sills 2005a; Peñarrubia et al. 2017), but they consider the scenario where many GCs are formed inside higher mass sub-halos²¹. Although that analysis is also interesting, a direct comparison with our work would be unfair, given such very different initial conditions.

4.4. General caveats

As the simulations we ran required a substantial computer power and running time, we needed to constrain some of our assumptions and exploratory analyses, in order to meet all the numerical convergence tests performed in Boldrini et al. (2020a). In this section, we briefly discuss novel aspects that can be envisaged in future works and how they could affect some of our results, by changing a few of our initial assumptions.

4.4.1. Mass profiles

Aiming for a more homogeneous treatment and interpretation, we assume all five GCs and respective mini-halos to have the same mass distributions. Although such assumption is reasonable, given the observed stellar masses of the Fornax dSph GCs (de Boer & Fraser 2016), the more general scenario of Galactic GCs for example encompasses a much broader mass and scale radii range (Baumgardt & Hilker 2018), which can deliver more different results. Another way to look into this variable is that all of our GCs have different filling factors²², and one could eventually envisage a numerical simulation where different GCs are formed with a similar filling factor.

In such scenario, clusters and respective mini-halos that lie closer to the host would be more compact than in our simulations, which does seem to happen to some extent in the Milky Way (see GC half-number radii from Vitral 2021 and orbital radii from Baumgardt et al. 2019, for example). By being more compact, inner clusters could have their DM shield effective for longer timescales, and thus present milder or no stellar streams at present age, as well as smaller degrees of rotation, ellipticity and higher velocity anisotropies, according to our findings. Similarly, if innermost GCs and respective mini-halos happen to be less massive, dynamical friction would be less prominent, and the GC+DM system would sink more slowly to regions of stronger tides. Again, the shielding consequences mentioned above would then last for longer in such clusters. We thus propose, for instance, that future works could explore a greater variety of GC and mini-halo mass profiles for clusters starting at similar positions with respect to their host.

4.4.2. Particle masses

In our simulation, all live objects (i.e. stars and DM particles) are set up with the same mass of 50 M_⊙. This numerical recipe helps to avoid possible heating and segregation effects that might affect the long-term internal evolution of the GCs embedded in DM. We thus assume the simplest case of equal particle mass in order not to bias our results in the case of heavier or lighter DM particles, and choose to study rather the external regions of the GCs. Indeed, mass segregation effects are less severe in external cluster regions due to a lower density of particles, which in turn limits the amount of dynamical interactions that lead to energy equipartition.

Nonetheless, if usual DM particles or sub-structures have masses considerably lower or higher than stars, our simulations could then fail to spot possible mass segregation effects specially related to the cluster internal dynamics. Meanwhile, the outer dynamics should not be considerably affected (as explained above), such that our results would remain valid up to a reasonable extent.

Mashchenko & Sills (2005b) studied the aspects of cluster relaxation, using a particle mass ratio of m_⋆•/m_• = {0.044, 0.44, 0.88}. They found that because the inner regions of the GC+DM system are dominated by stars (as we confirm in this study), the dynamics in the cluster’s core is not severely changed, and the impact of DM only starts to be visible at the radius where the amount of DM particles and stars are equivalent. At this point, energy equipartition leads the more massive DM particles to scatter the surrounding stars, creating a cutoff radius similar to what is observed in the density profiles of many Milky Way GCs (e.g. Trager et al. 1995). Equivalently, one could expect that less massive DM particles would themselves be scattered and deplete the mini-halo faster, thus decreasing the effectiveness of the shielding mechanism we demonstrated.

4.4.3. Different orbits

When comparing simulated GCs with and without a DM mini-halo in a tidal field, the important mass contribution from the mini-halo is such that the orbits of these two scenarios are bound to be different, as a consequence of dynamical friction (e.g. upper panels of Fig. 3). In such setup, it can be sometimes difficult to disentangle which dynamical effects are due to the DM mini-halo, or to the ongoing tidal interactions at a given time step. For that reason, we choose to analyse the evolution of these two scenarios starting from the same distance relative to the host, in an attempt to provide a fairer comparison. This forcibly leads clusters formed inside a DM mini-halo to inner regions of the host’s potential, where tidal forces are stronger.

Hence, all indications of milder tidal effects in GCs embedded in DM add robustness to our conclusion of a DM shielding mechanism. This is particularly true for the DM consequences we highlighted on stellar streams, velocity anisotropy, rotation and ellipticity. However, populating our simulations with more GCs would undeniably provide a more complete picture on the distribution of GC positions in IOM space at a given time, for instance. Also, by simulating clusters at higher distances from their host, they would suffer less dynamical friction²³, which could allow us to compare velocity dispersion profiles of clusters at more similar orbital radii. Due to the high cost of our simulations, we restrained such study to the five observed clusters in the Fornax dSph, so as to approach a realistic case, but adding similar clusters at a broader range of orbital radii seems to be a tantalising direction for future studies.

5. Summary and conclusions

We have used full N-body simulations on GPU of the globular cluster system in the Fornax dwarf spheroidal galaxy to probe the differences from globular clusters originally embedded in a dark matter mini-halo, and those devoid of it. For that, we simulated the two cases, where globular clusters have or not an enveloping dark matter mini-halo, using the initial conditions from Boldrini et al. (2020a), where all tracers (stars and dark matter) are assumed as N-body particles, for both the globular clusters and the host galaxy. Although the extra mass of the system containing dark matter accelerates orbital decay, sinking the clusters to regions of stronger tidal fields, the same extra mass triggers an increase of the tidal radius, which generally goes well beyond the limiting region of globular cluster bound stars (see Fig. 1). As a result, this tidal radius growth causes the dark matter mini-halo to work as a protective shield (see Fig. 2), which is itself stripped beforehand globular cluster stars. Below, we summarise some of the main diagnostics from our work.

• Globular clusters with a dark matter mini-halo should have, in general, more radial outer velocity anisotropy profiles, throughout all their orbit.

• Inflations on the outer velocity dispersion profile are expected in both clusters with and without a dark matter mini-halo, as a result of the host’s tidal field.

• Globular clusters with a dark matter mini-halo should have smaller degrees of internal rotation, and as a consequence, smaller ellipticities for their stellar distribution.

• Globular clusters without a dark matter mini-halo develop more prominent stellar streams, and present a more diffuse stellar distribution.

• Due to stronger dynamical friction, clusters originally embedded in a dark matter mini-halo evolve towards considerably lower energy regions of the integrals of motion space (differently from their dark matter-free counterparts), and cannot be reliably associated to previous accretion events due to their fast displacements in such space.

• Even clusters that retain large amounts of dark matter tend to have inner densities vastly dominated by baryonic components (up to a factor of a hundred in mass), and dark matter searches should target rather the cluster’s outskirts.

• We provide parametric mass profiles for dark matter mini-halos that evolved in a tidal field, which is of much usefulness for dark matter mini-halo searches employing Jeans mass-modelling (e.g. Carlberg & Grillmair 2022).

All in all, our work helps to confirm previous predictions (e.g. the shielding mechanism from dark matter mini-halos), but also challenges a few preceding points. For instance, we can mention the shape of stellar velocity dispersion profiles, which we show to be more dependent on the tidal field rather than the presence of a mini-halo. In addition, we also point to the long-term survival of clusters embedded in dark matter, which in some scenarios can be disrupted faster due to pronounced dynamical friction, since the latter leads them to regions of stronger tidal forces, where the shielding mechanism will vanish in shorter timescales, all while leaving the cluster closer to its host. Such new predictions are due mostly to our better handling of tidal interactions and dynamical friction, respectively. We also provide the first direct analysis pertaining to the evolution of the stellar velocity anisotropy, rotation, ellipticity and position in integrals of motion space for the dark matter mini-halo scenario. Finally, we evaluate some of the caveats of our analyses, and point towards novel aspects that can be further explored in more detail with the help of numerical simulations.

Movies

Movie 1 associated with Fig. 6 (movie) Access here

Movie 2 associated with Fig. 6 (movie_zoom) Access here

Acknowledgments

We thank the anonymous referee for the constructive report, with many insightful comments that have helped us to improve the quality of our results and clarify some descriptions in the manuscript. We also thank Geraint Lewis and Amaury Micheli for useful tips concerning the 3D visualisation of stellar streams. Eduardo Vitral was funded by an AMX doctoral grant from École Polytechnique. This work was supported by the EXPLORAGRAM Inria AeX grant. We greatly benefited from the public software PYTHON (Van Rossum & Drake 2009) packages BALROGO (Vitral 2021), SCIPY (Jones et al. 2001), NUMPY (van der Walt et al. 2011) and MATPLOTLIB (Hunter 2007). We also used GEOGEBRA (Hohenwarter 2002), as well as the SPYDER Integrated Development Environment (Raybaut 2009).

References

Appendix A: Velocity anisotropy

The Bayesian computation of the velocity anisotropy in BALROGO follows the procedure described in the following. First, we consider the local velocity distribution h(v, r) as a Gaussian product over all spherical coordinates (i.e. r, θ, ϕ), in a similar fashion than done in MAMPOSST-PM (Mamon et al. 2013; Read et al. 2021)

$$\begin{array}{c}\hfill h(\mathit{v},r)=\mathcal{G}[v_r,{\mu }_r,{\sigma }_r(r)]\mathcal{G}[v_{\theta },{\mu }_{\theta },{\sigma }_{\theta }(r)]\mathcal{G}[v_{\varphi },{\mu }_{\varphi },{\sigma }_{\varphi }(r)].\end{array}$$(A.1)

Differently than in MAMPOSST-PM, which solves the Jeans equation (Binney & Mamon 1982) to derive σ_r(r), we fit σ_r(r) with a generalisation of the parametric Plummer form from Dejonghe (1987, similarly used in Vasiliev 2019).

$$\begin{array}{c}\hfill {\displaystyle {\sigma }_r(r)=\frac{{\sigma }_0}{{[1+{(r/r_{\sigma })}^{\alpha }]}^{\gamma }},}\end{array}$$(A.2)

where σ₀ is the radial velocity dispersion at the cluster’s centre, r_σ a characteristic radius, α and γ characteristic exponents, which we allow to vary between {1, 8}, and {0.05, 2}, respectively. Then, σ_θ(r) and σ_ϕ(r) can be derived, by assuming spherical symmetry and using the definition of the velocity anisotropy (Eq. [1]), such that ${\sigma }_{\theta }^2(r)={\sigma }_{\varphi }^2(r)={\sigma }_r^2(r)[1-\beta (r)]$. We then use the gOM expression for the anisotropy (Eq. [2]) in order to write β(r). Since we do not use the Jeans equation, we can also account for the cluster’s systemic motions, and so we assign μ_i as the mean of the velocities in the i coordinate (differently from MAMPOSST-PM, which fixes μ_i = 0).

Finally, we minimise the logarithm of the likelihood, which is basically the expression from Eq. A.1, multiplied over each tracer. We do it by using a Markov chain Monte Carlo (MCMC) approach from the PYTHON package EMCEE (Foreman-Mackey et al. 2013), with 10⁴ steps²⁴ and 15 chains. In Figure A.1, we provide a goodness of fit of the velocity anisotropy of GC5 for the last snapshot from our simulations, for the case with and without a DM mini-halo.

Fig. A.1. Goodness of fit: Bayesian fit of the velocity anisotropy with BALROGO. The upper plot shows the case without a dark matter mini-halo, while the lower plot shows the case with a dark matter mini-halo. The black dots are the binned values of the empirical velocity anisotropy, logarithmically spaced and with error bars representing the 3-σ Poisson error. The solid coloured line (red or green) represent the best likelihood fit of the Osipkov (1979), Merritt (1985) model (Eq. [2]), with respective shaded regions encompassing the [2.5, 97.5] percentiles of the fit.

Appendix B: Velocity dispersion

We describe here how we constructed velocity dispersion profiles for the output of our simulations.

B.1. Dispersion map

The velocity dispersion maps displayed in Figures 2 and 8 was constructed by first binning the projected X vs. Y map with PYTHON’s HEXBIN routine, setting the argument gridsize = 100. For each bin, we then computed the 3D velocity dispersion of the DM component by summing quadratically the velocity dispersion on X, Y and Z directions. Whenever there were more than hundred particles inside the bin, we used all the bin’s particles, otherwise we completed the sample by picking the closest particles to the bin’s centre of mass, until the threshold of a hundred particles was attained.

The process above assured that each bin had a statistically significant number of tracers, which helped us to reach a better spatial resolution, eventually. However, as the results still presented an important amount of statistical noise, we decided to smooth our maps with PYTHON’s GAUSSIAN_FILTER routine, with the argument sigma = 3. We finally displayed the outcome of this procedure in a map colour-coded logarithmically from blue (lower dispersion) to red (higher dispersion).

B.2. Dispersion radial profiles

To build the radial dispersion profiles from Figure 4, we used the routines ANGLE.CART_TO_SPH and DYNAMICS.DISPERSION from the BALROGO PYTHON package, to first convert the data into spherical coordinates and then compute the velocity dispersion of the radial component, as a function of the distance to the cluster’s centre. The way we set BALROGO to compute this dispersion is by first dividing the radial extent into thirty equally spaced logarithmic bins and then calculating the respective dispersion and Poisson error associated to it.

Next, it smooths the profile in order to remove statistical noise, by fitting a 10-th order polynomial to it with the NUMPY.POLYFIT routine. 1-σ regions were obtained by considering the covariance matrix returned by the same algorithm. We finally chose to consider the results inside a more restrict spatial extent in order to neglect the potentially bad fits of the polynomial on the borders of our data.

Appendix C: Asymmetric surface density

Here we describe how the BALROGO’s POSITION.ELLIPSE_LIKELIHOOD routine performs asymmetric Sérsic fits with (α,  δ) data. First, it projects the (α,  δ) data according to classical spherical trigonometry relations, translating it so it is centred at the origin:

$$\begin{array}{cc}\hfill x_{\mathrm{p}}& =cos\delta sin(\alpha -{\alpha }_0),\hfill \end{array}$$(C.1a)

$$\begin{array}{cc}\hfill y_{\mathrm{p}}& =sin\delta cos{\delta }_0-cos\delta sin{\delta }_{0}cos(\alpha -{\alpha }_0),\hfill \end{array}$$(C.1b)

where (α₀,  δ₀) is the centre of the cluster. Next, it rotates the axis so the data can be easily handled:

$$\begin{array}{cc}\hfill x& =x_{\mathrm{p}}cos\theta +y_{\mathrm{p}}sin\theta ,\hfill \end{array}$$(C.2a)

$$\begin{array}{cc}\hfill y& =-x_{\mathrm{p}}sin\theta +y_{\mathrm{p}}cos\theta ,\hfill \end{array}$$(C.2b)

where θ is the angle between the original reference frame and the new one. With this new set, we are able to define a likelihood function of the stellar distribution as:

$$\begin{array}{c}\hfill \mathcal{L}={\displaystyle \prod _i}\frac{\mathrm{\Sigma }(m)}{N_{\mathrm{tot}}}·\end{array}$$(C.3)

with $m=\sqrt{{(x/a)}^2+{(y/b)}^2}$, and where (a,  b) are the semi-axis of the ellipse. The surface density Σ(m) and the number of tracers at infinity N_tot are defined such as in van de Ven & van der Wel (2021), so that we have:

$$\begin{array}{c}\hfill \frac{\mathrm{\Sigma }(m)}{N_{\mathrm{tot}}}=\frac{b_n^{2n}\exp [-b_n{m}^{1/n}]}{2\pi abn\mathrm{\Gamma }(2n)}·\end{array}$$(C.4)

In equation C.4, n is the Sérsic index, Γ(x) is the gamma function of the variable x and the b_n is computed with the precise approximation from Ciotti & Bertin (1999). Thus, the fit finds the parameters which maximise ℒ.

Finally, to derive statistical errors of our Bayesian estimates, we use PYTHON’s NUMDIFFTOOLS.HESSIAN method to compute the Hessian matrix of the probability distribution function (i.e. Eq. [C.4]). After, we assign the uncertainties of each parameter as the square root of the respective diagonal position of the inverted Hessian matrix.

Appendix D: Extra material

Fig. D.1. Evolution of scale radii: Evolution of the GC+DM system tidal radius (rt, blue) and of the half-mass radius of the bound dark matter particles (r1/2, •, black), for each simulated cluster. We observe that clusters who lost most of their dark matter mini-halo (i.e. GC1 and GC2) presented rt >  r1/2, • at some point of their orbital evolution, similarly to the relation for sub-halo disruption from Hayashi et al. (2003).

Fig. D.2. Evolution of the dark mass: Evolution of the total mass of bound dark matter particles, for each simulated cluster. Bound particles were defined according to Section 2.3.

Fig. D.3. Evolution of the stellar mass: Evolution of the total stellar mass of bound globular cluster particles, for each simulated cluster. Clusters originally embedded in a dark matter mini-halo (i.e. ⋆•) are displayed in green, while those devoid of the mini-halo (i.e. ⋆) are represented in red. Bound particles were defined according to Section 2.3.

Fig. D.4. Evolution of the half mass radius: Evolution of the half mass radius of bound globular cluster particles, for each simulated cluster. Clusters originally embedded in a dark matter mini-halo (i.e. ⋆•) are displayed in green, while those devoid of the mini-halo (i.e. ⋆) are represented in red. Bound particles were defined according to Section 2.3.

All Tables

All Figures

Fig. 1. Tidal radius growth due to an extra dark mass: Scheme illustrating the increase of tidal radius (blue) in the globular cluster system (orange) embedded in a dark matter mini-halo (grey), in the left, in comparison with the globular cluster system without dark matter, in the right, at a same orbital radius. This tidal radial expansion is mostly explained by its dependence on the satellite’s mass, which is greater for the case with a DM mini-halo.

In the text

Fig. 2. Dark matter shield: Velocity dispersion map of dark matter particles for GC5, projected in the X vs. Y plane and centred in the centre of mass of the globular cluster system. We display the last six pericentres of its orbit, where the tidal effects are stronger. The extension of bound globular cluster stars and bound dark matter particles are highlighted as dotted and dashed green lines, respectively, while the theoretical tidal radius, calculated according to Sect. 2.3, is displayed as a solid green circle. The maps are colour-coded logarithmically from blue (lower dispersion) to red (higher dispersion). The centres of Fornax and of the globular cluster are represented as a thick green cross and a plus sign, respectively. For this cluster, we notice that the empirical tidal radius, well traced by the blue region, remains always larger than the bound stars radii. This argues in favour of the dynamical presence of a dark matter shield, as illustrated in Fig. 1.

In the text

Fig. 3. Velocity anisotropy: Evolution of the orbital radius (i.e. rorbit, upper plots) and of the outer velocity anisotropy (i.e. β∞, lower plots), for the five simulated globular clusters. The clusters originally embedded in a dark matter mini-halo are depicted in solid green, and those without it are displayed in solid red. For β∞, we display the running mean over the five closest snapshots (time-wise), for better visualisation, which explains why the initial β∞ values seem to be different for the two scenarios. We observe typically lower rorbit and more radial (i.e. higher) β∞ for clusters with a dark matter mini-halo.

In the text

Fig. 4. Stellar velocity dispersion: Series of plots for GC5. The upper plots relate to the simulations devoid of dark matter mini-halos (⋆), while the lower plots indicate the results for the globular clusters formed in such mini-halos (⋆•). The two columns on the left display hand-picked snapshots where the radial velocity dispersion profile (i.e, σr(r)) resembled better to an isolated case without a dark matter mini-halo (left), and with a massive dark matter mini-halo (right), according to Fig. 2 from Peñarrubia et al. (2017). The column on the right presents the evolution of the maximum value of the radial velocity dispersion profile (i.e, σrmax), for each scenario concerning the dark matter mini-halo, colour-coded according to the distance of the cluster to the centre of the host galaxy (i.e. rorbit), with two vertical dashed lines corresponding to the instants from the two columns on the left. These plots argue that the tidal field from the host galaxy tends to have a much greater impact on inflating the velocity dispersion than the presence of a dark matter mini-halo. In fact, such mini-halos help to protect the cluster from tidal effects, rather than contributing to it.

In the text

Fig. 5. Rotation: On the top, we show the mean evolution of the line-of-sight (LOS) velocity dispersion of bound globular cluster particles, for each simulated cluster, while the bottom plot displays the mean ellipticity (or flattening) of the stellar distribution. We display the running mean over the ten closest snapshots (time-wise), for better visualisation. Error bars are represent 1-σ Poisson uncertainties. Clusters originally embedded in a dark matter mini-halo (i.e. ⋆•) are displayed in green, while those devoid of the mini-halo (i.e. ⋆) are represented in red. As a general trend, globular clusters embedded in dark matter mini-halos exchange less internal angular momentum from tidal interactions, thus presenting smaller degrees of rotation, which in turn is associated with smaller ellipticities. This again attests the efficacy of the dark matter shield illustrated in Fig. 2.

In the text

Fig. 6. Dark matter influence on tidal tails: Stellar distribution (and respective zoom in) of GC5 for the case where it was originally embedded in a dark matter mini-halo (right) and where it was devoid of it (left), both centred in the GC centre of mass. The chosen snapshot was such that GC5 was located at roughly the same orbital radius for both scenarios, for a fairer comparison. This plot highlights more prominent tails in the dark matter-free case, while clusters formed inside dark matter mini-halos present a more compact stellar envelope. Hence, the dark matter shield (Sect. 3.1) has the effect of delaying the formation of tidal tails. Movies are available online.

In the text

Fig. 7. Integrals of motion: Evolution of our simulated clusters in the integrals of motion (IOM) space, in particular, the E × LZ diagram. Their positions in this diagram are scattered in a sequential colour-map, starting from the same point in light grey at t = 0, and ending in darker tones (green for globular clusters originally formed in dark matter mini-halos, and red for clusters devoid of dark matter). The last snapshot is marked by a cross, for each cluster. This plot highlights a bimodal evolutionary distribution of clusters originally embedded in dark matter (moving towards lower energies) and those who were not (moving only slightly towards higher energies), and more importantly, it shows that clusters originally embedded in dark matter move significantly in IOM space, such that their association with past merger events through this diagram is not reliable.

In the text

Fig. 8. Survival of the dark matter mini-halo: Similarly to Fig. 2, we show the velocity dispersion map of dark matter particles for the five globular clusters we simulated, projected in the X vs. Y plane (with exception of GC4, projected on the Z vs. Y plane for better visualisation of its orbit) and centred in the centre of mass of the globular cluster system. We display the last pericentre passages of their respective orbits, where the tidal effects are stronger. One can notice that while GCs 3, 4 and 5 manage to keep the dark matter mini-halo (and thus the dark matter shield effect depicted in Fig. 2), GCs 1 and 2 do not. The extent of the red region for GCs 1 and 2 is related to their different orbital parameters. For GCs 3, 4 and 5, we again notice that the empirical tidal radius is well traced by the internal blue region.

In the text

Fig. 9. Cumulative mass profiles: We display the cumulative mass profiles at the last snapshot of our simulations, for the dark matter particles (black), their embedded stellar component (green) and the stellar component from the simulations without dark matter (dashed red). The blue line depicts the MCMC fit of a Zhao (1996)αβγ model to the dark matter mini-halo (see Sect. 4.2.2), with a fainter blue region encompassing the [2.5, 97.5] percentiles. The dotted lines represent local slopes of the stellar mass distribution at the respective regions shown in the plot.

In the text

Fig. A.1. Goodness of fit: Bayesian fit of the velocity anisotropy with BALROGO. The upper plot shows the case without a dark matter mini-halo, while the lower plot shows the case with a dark matter mini-halo. The black dots are the binned values of the empirical velocity anisotropy, logarithmically spaced and with error bars representing the 3-σ Poisson error. The solid coloured line (red or green) represent the best likelihood fit of the Osipkov (1979), Merritt (1985) model (Eq. [2]), with respective shaded regions encompassing the [2.5, 97.5] percentiles of the fit.

In the text

Fig. D.1. Evolution of scale radii: Evolution of the GC+DM system tidal radius (rt, blue) and of the half-mass radius of the bound dark matter particles (r1/2, •, black), for each simulated cluster. We observe that clusters who lost most of their dark matter mini-halo (i.e. GC1 and GC2) presented rt >  r1/2, • at some point of their orbital evolution, similarly to the relation for sub-halo disruption from Hayashi et al. (2003).

In the text

	Fig. D.2. Evolution of the dark mass: Evolution of the total mass of bound dark matter particles, for each simulated cluster. Bound particles were defined according to Section 2.3.
In the text

	Fig. D.3. Evolution of the stellar mass: Evolution of the total stellar mass of bound globular cluster particles, for each simulated cluster. Clusters originally embedded in a dark matter mini-halo (i.e. ⋆•) are displayed in green, while those devoid of the mini-halo (i.e. ⋆) are represented in red. Bound particles were defined according to Section 2.3.
In the text

	Fig. D.4. Evolution of the half mass radius: Evolution of the half mass radius of bound globular cluster particles, for each simulated cluster. Clusters originally embedded in a dark matter mini-halo (i.e. ⋆•) are displayed in green, while those devoid of the mini-halo (i.e. ⋆) are represented in red. Bound particles were defined according to Section 2.3.
In the text