Constraining the shape of dark matter haloes with globular clusters and diffuse stellar light in the E-MOSAICS simulations

ABSTRACT

We explore how diffuse stellar light and globular clusters (GCs) can be used to trace the matter distribution of their host halo using an observational methodology. For this, we use 117 simulated dark matter (DM) haloes from the |$(34.4~\rm cMpc)^3$| periodic volume of the E-MOSAICS project. For each halo, we compare the stellar surface brightness and GC projected number density maps to the surface density of DM. We find that the dominant structures identified in the stellar light and GCs correspond closely with those from the DM. Our method is unaffected by the presence of satellites and its precision improves with fainter GC samples. We recover tight relations between the dimensionless profiles of stellar-to-DM surface density and GC-to-DM surface density, suggesting that the profile of DM can be accurately recovered from the stars and GCs (σ ≤ 0.5 dex). We quantify the projected morphology of DM, stars, and GCs and find that the stars and GCs are more flattened than the DM. Additionally, the semimajor axes of the distribution of stars and GCs are typically misaligned by ∼10 degrees from that of DM. We demonstrate that deep imaging of diffuse stellar light and GCs can place constraints on the shape, profile, and orientation of their host halo. These results extend down to haloes with central galaxies M_⋆ ≥ 10¹⁰ M_⊙, and the analysis will be applicable to future data from the Euclid, Roman, and the Rubin observatories

1 INTRODUCTION

Most of the mass in dark matter (DM) haloes lies at large distances from the centre, so probing the galactic outskirts beyond the extent of the visible galaxy is critical to trace the structure of the halo. In our current Lambda Cold Dark Matter paradigm, galaxies are built via a combination of in situ star formation and hierarchical accretion of satellite galaxies. The relative role of these mechanisms changes with galaxy mass: the fraction of accreted mass in galaxies increases strongly with their mass (e.g. Abadi, Navarro & Steinmetz 2006; Genzel et al. 2010; Qu et al. 2017; Behroozi et al. 2019). In particular, the stellar component of galaxy clusters pre-dominantly assembles via the accretion of satellites. An outcome of this hierarchical build-up is the presence of diffuse stellar light surrounding the central galaxy, which permeates the space between the galaxies within the cluster (the ‘intracluster light’; e.g. Contini 2021; Montes 2022).

Since both the DM and stars are collisionless, the large spatial extent covered by the intracluster light has led to it being posited as a viable tracer of the outer matter distribution of the host halo (e.g. Sampaio-Santos et al. 2021). Following this idea, Montes & Trujillo (2019) demonstrate that the diffuse stellar light in six galaxy clusters follows more closely the mass distribution estimated from gravitational lensing compared to that inferred from X-ray observations. This striking result suggests that deep imaging might be sufficient to characterize the distribution of DM in galaxy clusters. Expanding on this work, Alonso Asensio et al. (2020) use the Cluster-Evolution and Assembly of GaLaxies and their Environments (EAGLE) simulations to suggest that the diffuse stellar light is an even a better tracer of the matter distribution than observations had previously suggested and propose an indirect method to obtain the halo mass and its mass density profile from the radial profile of the intracluster mass density.

Along with the presence of diffuse stellar light, observations have also confirmed the presence of populations of intergalactic globular clusters (GCs) out to large distances in a variety of galaxy systems, e.g. the galaxy cluster Abell 1689 (e.g. Alamo-Martínez et al. 2013; Alamo-Martínez & Blakeslee 2017), Virgo cluster (e.g. Durrell et al. 2014), Perseus cluster (e.g. Harris et al. 2020b), Fornax cluster (e.g. Chaturvedi et al. 2022), and even the M81/M82/NGC3077 system of galaxies (Chies-Santos et al. 2022). When compared to the X-ray surface brightness distribution, the GCs in the Abell 1689 cluster show more substructure (Alamo-Martínez et al. 2013) but those in Virgo show a similar distribution to the hot gas (Durrell et al. 2014). These intergalactic GCs are also likely an outcome of the hierarchical build-up of their current halo: these objects were stripped from their host galaxies as they were accreted, and so they now populate the outskirts of the halo. Together with the intracluster light, these GC populations are a relic of the formation and assembly history of their host cluster.

It has been established observationally that more massive haloes host a larger number of GCs (e.g. Blakeslee 1997; Peng et al. 2008; Spitler & Forbes 2009; Georgiev et al. 2010; Hudson, Harris & Harris 2014; Harris, Harris & Hudson 2015; Harris, Blakeslee & Harris 2017; Forbes et al. 2018). Together with the increasingly larger fraction of accreted mass in more massive haloes (e.g. Behroozi et al. 2019), one could expect that GCs can hold as much information about the matter distribution as the diffuse stellar light in galactic systems in which GCs are very numerous.

Additionally, GCs are interesting objects from an observational perspective for a number of reasons. First, their high luminosities (i.e. the peak being at M_V ≃ −7.3, e.g. Harris et al. 2014) and compact sizes (∼3–|$4~\mbox{${\rm pc}$}$|⁠, e.g. Krumholz, McKee & Bland -Hawthorn 2019) imply that their surface brightnesses are higher than that of the diffuse stellar light. Secondly, their compact nature makes them appear as bright point sources in extragalactic surveys, and with colour information they are easier to distinguish from other objects such as background galaxies (e.g. Muñoz et al. 2014; Brito-Silva et al. 2021). Lastly, their high luminosities mean that they can be observed further out from galaxy centres, i.e. observations typically extend ∼5–20 times the stellar effective radius of the galaxy, depending on its mass (e.g. Alabi et al. 2016). At those large radii, GCs are more likely to be tracing the outer matter distribution. Overall, GC populations can be easier and quicker to observe than the diffuse stellar light in their host haloes. For these reasons, we examine both stellar and GC populations as possible tracers of the matter distribution of their host halo.

In addition to determining the shape of the DM halo, prior work has also explored correlations between the DM and baryonic morphologies, and their internal misalignments using numerical simulations of galaxy formation and evolution (e.g. Deason et al. 2011; Tenneti et al. 2014; Velliscig et al. 2015a, b; Pillepich et al. 2019; Thob et al. 2019; Hill et al. 2021). These studies aim to provide a theoretical interpretation for upcoming measurements of the apparent alignment of galaxy shapes due to the gravitational lensing effect caused by the underlying matter distribution (e.g. from the Euclid and Rubin observatories). Using the eagle simulations, the stellar component of galaxies is found to be well aligned with the local mass distribution, but it is often misaligned significantly with respect to the global halo (Velliscig et al. 2015a; Hill et al. 2021). This indicates that stars follow the local DM distribution, and that the orientation of the DM distribution changes from the inner to the outer halo. In contrast, the hot and diffuse gas in the halo (taken as a proxy for X-ray emitting gas) is found to be significantly misaligned relative to the local matter distribution. This give credence to the suggestion that the diffuse stellar light in haloes can be a good tracer of their matter distribution.

Building upon these previous studies, we explore whether diffuse stellar light and GC populations correlate with the outer matter distribution of their host halo. We also examine if the DM surface density profile can be recovered from these observational tracers, and we extend the analysis of the projected morphologies of DM and stars to the GC populations. The goal of the current study is to test these correlations using hydrodynamical numerical simulations of galaxy formation that include star clusters. In the local Universe (⁠|$D\lt 50~\mbox{${\rm Mpc}$}$|⁠), GC populations do not require deep observations to be characterized and might therefore provide a more efficient method to track their host DM distribution.

In future analysis, we will compare the predictions of this method to observations for which studies of the intracluster GCs, diffuse stellar light, and DM distributions exist. Some of these include the galaxy cluster Abell 1689 (Alamo-Martínez et al. 2013), Perseus cluster (e.g. Harris et al. 2020b), Fornax cluster (e.g. Chaturvedi et al. 2022), and, recently observed with the James Webb Space Telescope, the galaxy cluster SMACS0723 (Diego et al. 2023). Once the method developed using simulations is tested against direct DM observations, it can be applied to the very large samples of galaxies obtained by upcoming deep and wide surveys with the Euclid, Roman, and Rubin observatories.

In this work, we use the |$(34.4~\rm cMpc)^3$| periodic volume from the MOdelling Star Clusters Populations Assembly In Cosmological Simulations within eagle project (E-MOSAICS project; Pfeffer et al. 2018; Kruijssen et al. 2019a). The combination of a subgrid description for stellar cluster and formation with the eagle galaxy formation model (Crain et al. 2015; Schaye et al. 2015) allows us to study the formation and assembly of GC populations together with their host galaxies. For this, we select the 117 DM haloes with central galaxies more massive than M_⋆ ≥ 10¹⁰ M_⊙, which correspond to halo masses M₂₀₀ ≳ 4 × 10¹¹ M_⊙. By considering lower-mass haloes that do not host galaxy clusters, we aim to investigate whether the agreement between the diffuse light and the matter distribution still holds at those masses.

The simulation setup is described in Section 2. We follow the methodology of Montes & Trujillo (2019) and explore whether the 2D structures that can be identified in stellar surface brightness and GC number density maps correspond to the underlying structures in the DM distribution. The generation of the projected maps is discussed in Section 3. We explore the similarity of structures identified in the DM, stellar, and GC maps in Section 4 and examine how the DM surface density radial profile can be recovered from the stellar surface brightness and GC number density maps in Section 5. We quantify the projected morphology of DM, stars, and GCs in Section 6, and the summary is presented in Section 7.

2 SIMULATIONS

2.1 The E-MOSAICS project

The E-MOSAICS project (Pfeffer et al. 2018; Kruijssen et al. 2019a) combines a subgrid description for stellar cluster formation and evolution with the eagle galaxy formation model (Crain et al. 2015; Schaye et al. 2015). The combination of these models allows us to study the formation and evolution of stellar cluster populations alongside their host galaxies. We briefly summarize the models describing the subgrid star clusters, and we refer the reader to Pfeffer et al. (2018) and Kruijssen et al. (2019a) for further details.

Every time a star particle is formed in our simulations, a subgrid stellar cluster population can be spawned within it. The formation of the subgrid clusters is regulated by two physical ingredients: the cluster formation efficiency and the upper mass scale truncation of the Schechter (1976) initial cluster mass function. We model these two ingredients with environmentally dependent descriptions, in which natal environments with higher gas pressures lead to more mass forming in clusters and more massive objects (Kruijssen 2012; Reina-Campos & Kruijssen 2017). Once the subgrid clusters have formed, they evolve alongside their parent star particle. We model their mass evolution due to stellar evolution following the EAGLE model (Wiersma et al. 2009), and we model three additional dynamical disruption mechanisms. In our simulations, clusters lose mass due to sudden tidal shocks with their surrounding environment (Kruijssen et al. 2011), due to two-body interactions within the cluster, and due to their in-spiral towards the centre of galaxies caused by dynamical friction. Given the subgrid nature of the stellar cluster populations, the latter mechanism is only applied in post-processing.

The eagle galaxy formation model reproduces many key properties of the observed galaxy population. Of particular interest for this project, the simulations reproduce the evolution of galaxy sizes across time (Furlong et al. 2017), the present-day luminosities and colours of galaxies (Trayford et al. 2015), and the evolution of the galaxy stellar masses and their star formation rates (Furlong et al. 2015). Additionally, the simulations have been extensively used to characterize the distributions of DM, hot gas and stars, and their relationship to one another (e.g. Velliscig et al. 2015a, b; Thob et al. 2019; Hill et al. 2021, 2022).

The E-MOSAICS simulations have been successful at reproducing the properties of old and young star clusters in the local Universe (Kruijssen et al. 2019a; Pfeffer et al. 2019b; Hughes et al. 2020). Among other properties, these simulated cluster populations reproduce the ‘blue tilt’ (i.e. the overabundance of massive red clusters; Usher et al. 2018), the high-mass end of the GC mass function (Pfeffer et al. 2018; Kruijssen et al. 2019a), and the metallicity distributions across galaxy stellar mass (Pfeffer et al. 2023). In particular, the number density radial profiles of GCs are also in very good agreement to observations for a wide range of galaxy stellar masses (Reina-Campos et al. 2021). These simulations have also informed predictions for the high-redshift conditions for cluster formation (Reina-Campos et al. 2019; Pfeffer et al. 2019a; Keller et al. 2020), and they have showcased the potential of using GCs to trace the assembly history of their host galaxy (Hughes et al. 2019; Kruijssen et al. 2019b, 2020; Pfeffer et al. 2020; Trujillo-Gomez et al. 2021). Thus, this data set provides an excellent framework to explore whether diffuse stellar and GC populations correlate with the overall matter distribution of their host halo.

2.2 Sample selection

We study the stellar clusters and their host DM haloes from the |$(34.4~\rm cMpc)^3$| periodic volume evolved within the E-MOSAICS project, as featured in Bastian et al. (2020). We identify DM haloes in the periodic volume using the Friends-of-Friends algorithm (FoF; Davis et al. 1985) with a linking length of 0.2 times the mean particle separation. Next, we associate gas and stellar particles to the nearest DM particle, and we identify gravitationally bound susbtructures within each halo with the SUBFIND algorithm (Springel, Yoshida & White 2001; Dolag et al. 2009). We select the DM haloes that contain central galaxies more massive than M_⋆ ≥ 10¹⁰ M_⊙, which correspond to halo masses M₂₀₀ ≳ 4 × 10¹¹ M_⊙.¹ This selection yields a sample of 117 haloes.

In our fiducial selection of particles, we include all the particles that are bound and unbound to the main halo (i.e. FOF group), as determined by the FoF algorithm. This selection includes all the particles in the central galaxy and in the diffuse unbound component, as well as those locked in the satellite galaxies within the halo.

3 PROJECTED DISTRIBUTIONS OF DM, STARS, AND GCS

We present the projected distributions of DM, stars, and GCs in our sample of haloes, and we explore correlations between the observational tracers (i.e. stars and GCs) and the matter component (i.e. DM). We produce mock images of our tracers to mimic observational biases and provide a method that can be used directly on observational data.

3.1 Converting physical into observational quantities

To mimic observational selection criteria, we calculate the luminosity emitted by stars and stellar clusters in different filters. To do so, we assume that our stars and clusters correspond to single stellar populations and we use the stellar and cluster ages, metallicities, and masses together with the Flexible Stellar Population Synthesis (FSPS) model (Conroy, Gunn & White 2009; Conroy, White & Gunn 2010; Conroy & Gunn 2010a, b) to calculate the absolute magnitudes. We further use the MILES spectral library (Sánchez-Blázquez et al. 2006), as well as Padova isochrones (Girardi et al. 2000; Marigo & Girardi 2007; Marigo et al. 2008), a Chabrier (2003) initial stellar mass function, and the default FSPS parameters.

Since stellar populations in the galactic outskirts are more prominent in redder bands, we calculate the stellar luminosities in the SDSS r-band. In the case of the stellar clusters, we calculate their absolute magnitudes in the F336W, F475W, and F814W bands of the Hubble Space Telescope (HST) ACS filter system. We use these HST bands as they have been widely used to explore extragalactic GC populations in different galactic environments (e.g. Alamo-Martínez et al. 2013; Harris et al. 2014; Harris 2016; Alamo-Martínez & Blakeslee 2017). These HST bands roughly correspond to the u, g, and i + z filters, respectively, in the SDSS band system. All of these magnitudes are calculated in the AB system.

3.2 Producing the projected maps

We produce mock images of our tracers by doing projections of their halo. These simulated projected maps do not contain foreground nor background stars or galaxies, as would be the case in real observations.

First, we create the maps that represent the true matter distribution within each halo. We restrict our analysis to the DM distribution, as it represents ∼84 per cent of the matter in the Universe (Planck Collaboration 2020). We create these maps by selecting the DM particles in the host halo, and using the pynbody.plot.image.² routine to produce the projected surface density images. This routine uses a k-dimensional tree to interpolate and smooth the projected densities.

Next, we create the maps corresponding to the observational tracers. We start by producing the stellar surface brightness maps. We do so by converting the absolute magnitudes in the SDSS r–band to luminosities using a solar absolute magnitude of |${\rm M}_{\odot }^{r-{\rm band}} = 4.65$| (Willmer 2018). We then use the pynbody.plot.image routine to produce projected luminosity maps, which we convert to surface brightness maps by accounting for the angular area of a pixel in the image. Lastly, we keep all the pixels with surface brightness brighter than |$\mu _{r-{\rm band}} \lt 28~\rm mag\, arcsec^{-2}$| to mimic observational constraints. We explore the effect of the surface brightness limit in Section 4.3.

Finally, we create the number density maps of GC populations. To facilitate direct comparisons with observations, we select the GCs using observational constraints. First, we select objects more metal-rich than |$\mbox{$[{\rm Fe}/{\rm H}]$}\gt -3$|⁠.³ Then, we apply an absolute magnitude cut in the F475W band to select objects brighter than the peak of the GC mass function (M ∼ 10⁵ M_⊙). We use a magnitude cut of M_F475W < −6.2, which corresponds to an |$8~\mbox{${\rm Gyr}$}$| old stellar cluster of solar metallicity (for a solar absolute magnitude of |${\rm M}_{\odot }^{475W} = 5.09$|⁠; Willmer 2018).⁴ The resulting cluster sample also includes lower mass clusters than the peak of the mass function due to the metallicity dependence of the mass-to-light ratio. We examine the effect of the absolute magnitude cut in Section 4.4.

Lastly, we also apply cuts to our cluster sample in colour combinations with a UV filter (F336W). As previous studies have shown, the lack of a description for the cold gas phase of the interstellar medium in the eagle galaxy formation model implies that the low-mass, young, and metal-rich stellar clusters disrupt too slowly in the E-MOSAICS simulations (see appendix D of Kruijssen et al. 2019a). Thus, the selection in magnitude and colour space allows us to efficiently detect the contaminant objects and remove them from our samples. Identifying GCs in colour space leads to samples with more complex demographics than if simple age and metallicity cuts were applied (e.g. Brito-Silva et al. 2021). We illustrate the procedure in Fig. 1. First, we calculate the colours F336W − F475W and F475W − F814W. Then, we motivate the colour cuts by examining the ranges spanned by those clusters older than |$\gt 8~\mbox{${\rm Gyr}$}$|⁠. Using 10 bins in the F475W − F814W colour, we fit linear regressions to the F336W − F475W colours within each bin, and we calculate the standard deviation around the fits (see small inset in the top panel of Fig. 1). We then keep those clusters whose F336W − F475W colour lies within |$2\, \sigma$| from the linear fit in their bin. We have verified that these selection criteria lead to realistic cluster populations in all the haloes in our sample. Finally, we assume that GCs are point sources and create the projected number density maps of GC populations using the numpy.histogram2d routine.

The requirement of reaching the peak of the GC luminosity function limits the distances of the DM haloes for which we will be able to apply this analysis in observational data. We can estimate these distances by considering the apparent magnitude limits expected to be reached by the Euclid Wide Survey (m = 26.2 in the visible band; Scaramella et al. 2021), and the Rubin observatory (m = 27.4 in the g–band after 10 yr of integration; Ivezić et al. 2019) as these provide bracketing limits for upcoming observational facilities. For our fiducial absolute magnitude cut of M_F475W < −6.2, these observatories will be able to observe GC populations at the required completeness level in haloes up to |$30.4~\mbox{${\rm Mpc}$}$| and |$52.9~\mbox{${\rm Mpc}$}$| away, respectively. Given that within our simulated (34.4 cMpc)³ volume there are 117 haloes of interest, one could naively expect to find ∼340 and ∼1800 such haloes within spheres of radius |$30.4~\mbox{${\rm Mpc}$}$| and |$52.9~\mbox{${\rm Mpc}$}$|⁠, respectively.

As a last step in the production of the maps, we smooth them with a Gaussian filter of kernel size |$\sigma =0.02\, r_{200}$|⁠. This scale roughly corresponds to the effective radii of the central galaxy (Kravtsov 2013), and it removes the small-scale noise from the images while preserving the large-scale structures that we aim to identify. We show in Appendix  A the influence of the kernel in the structures that can be identified in the images.

The resulting smoothed projected surface density map of DM, the stellar surface brightness map, and the number density map of GCs for FOF000 are shown in Fig. 2. This is the most massive halo in our sample, with a mass log₁₀(M₂₀₀/M_⊙) = 13.7 and a virial radius of |$r_{200} = 783.8~\mbox{${\rm kpc}$}$|⁠. Its central galaxy is a giant elliptical of stellar mass log₁₀(M_⋆/M_⊙) = 11.3 that seems to be undergoing a major merger at z = 0. Encouragingly, the DM halo substructures are clearly seen in the observational tracers. This halo also exhibits a diffuse distribution of DM, stars, and GCs around the galaxies. This diffuse component is a relic of the assembly history: stars and GCs were stripped during the accretion of their galaxy onto the central galaxy and populate the outskirts of the halo. We find similar diffuse features in all the haloes in our sample.

3.3 Measuring the correlation between projected maps

As a first test, we examine whether the structures present in pairs of images are analogous. To do so, we measure the Pearson correlation coefficient of pixel intensities among pairs of smoothed maps.⁵ We only consider those pixels in the images with significant signal, and we take the logarithm of the maps of DM and GCs.

We calculate these coefficients among our sample of 117 DM haloes, and we show the resulting values in Fig. 3. Stellar light shows the highest correlation with the underlying DM distribution, with a median coefficient of r_P = 0.92 across our haloes. In comparison, number density maps of GCs show a lower degree of correlation with DM (although still statistically significant), and this correlation decreases with decreasing galaxy mass. The median correlation coefficient is r_P = 0.54. The correlation between GCs and stars is slightly stronger, with r_P = 0.63. The p-values of all the correlations are consistent with being null. Thus, the projected structure of the stellar and GC components shows a strong correlation with the DM.

3.4 Identifying the isodensity contours

To compare the projected maps of stellar light and GCs against the DM maps, we identify structures defined by isodensity contours. We do this following the methodology described in Montes & Trujillo (2019), which we illustrate with GC maps in Fig. 4.

Dividing the images in 32 radial bins that are logarithmically spaced, we calculate the median value of the density (or surface brightness) within each bin. We interpolate the measured radial profile to obtain its value at seven radial distances, which correspond to 1 per cent, 5 per cent, 10 per cent, 20 per cent, 40 per cent, 70 per cent, and 90 per cent of r₂₀₀. Then, we identify the structures in each map by joining pixels with a value equal to the interpolated value from the radial profile using the matplotlib.contour routine. At each radial distance, there can be several identified structures, and we use the longest contour at each level as the dominant structure. As we explore below, the presence of ongoing mergers can cause the algorithm to choose different structures in two maps for the same contour.⁶

4 RECONSTRUCTING THE SHAPE OF THE DM HALO

Here, we explore whether the dominant isodensity contours identified in the stellar and GCs maps (i.e. observational distributions) correspond to the isodensity structures identified in the DM maps (i.e. true matter distribution). A qualitative comparison of the isodensity contours in three representative haloes is done in Fig. 5.

In some haloes (e.g. FOF000), the presence of massive subhaloes within the halo implies that different structures are identified depending on the tracer. Based on the stellar-light contour at 0.1r₂₀₀, the central galaxy in FOF000 dominates, but from the GC distribution, we can identify a diffuse extended population associated with the prominent satellite (at x/r₂₀₀, y/r₂₀₀ ≃ 0.4, 0). Ongoing major mergers can heavily distort the identification of the isodensity contours, as the incoming massive subhalo can break the assumption of smooth contours. Interestingly, these cases are very rare among our sample of simulated haloes at z = 0.

Interestingly, the isodensity contours identified from the stellar and GC maps follow closely that of the DM distribution in the majority of haloes that we examine (e.g. FOF042 and FOF166). As we consider lower mass haloes, we find that the stellar and GC distributions are restricted to the inner regions, i.e. the isodensity contours can only be drawn for the inner 3–5 radial distances. This is related to these haloes having undergone lower levels of satellite accretion in their assembly histories than more massive haloes, and so have formed less extended haloes.

4.1 The Modified Hausdorff Distance

We now perform a quantitative analysis of the similarity between the isodensity contours identified in the projected maps of each component. For that, we use the Modified Hausdorff Distance (MHD; Huttenlocher, Klanderman & Rucklidge 1993; Dubuisson & Jain 1994), which is a measure of how far apart two subsets are. This metric is a modification of the Hausdorff Distance that prevents the contamination from outliers at large distances, and it is commonly used in shape matching. For two sets of points, X = {x₁, x₂,..., x_N} and Y = {y₁, y₂,..., y_N}, the MHD is defined as

The distance d(X, Y) is calculated as

and it corresponds to the closest distance from a point on the contour X to the contour Y, averaged over all points on the contour X.

We calculate the MHD between pairs of isodensity contours identified in observational maps (i.e. stars and GCs) and the DM distributions among our sample of galaxies, which we show in Fig. 6. We focus on the contours drawn at 0.05 and |$0.2\, r_{200}$| to have sufficient haloes with both stellar and GC contours to enable meaningful comparisons. We quantify the agreement between the observational and true isodensity contours with the median and 25–75th percentiles, which we summarize in Table 1. A smaller median MHD means that the observational isodensity contour is more accurate, but a smaller interquartile range indicates that the observational tracer is more reliable overall.

Comparing the isodensity contours drawn at R/r₂₀₀ = 0.05, we find that the stellar light contours differ from the underlying DM contours by a median |$0.003\, r_{200}$| (with percentiles 0.002–|$0.005\, r_{200}$|⁠) among our sample of haloes. In contrast, the isodensity contours from GCs differ from the DM distribution by a median |$0.002\, r_{200}$| (with percentiles 0.001–|$0.003\, r_{200}$|⁠). The smaller median and scatter indicates that the projected distribution of GCs correlate more closely to the underlying DM distribution than stars in the inner region of the halo. We also measure the distance between the isodensity contours of GCs and stellar light (green points in Fig. 6), and find that the median is a factor of two larger than to the DM, but the percentile range is comparable.

Among our haloes, there are three outliers in which the measured MHDs are larger than the radius at which the contour is drawn (red shaded region in Fig. 6). There is one system in the distances measured to the stellar maps, and there are two outliers in the GCs maps. This indicates that identifying structures in the stellar surface brightness maps or in the number density maps of GCs is equally likely to be adversely affected by the presence of satellites. Overall, caution is necessary when identifying isodensity structures in haloes that are undergoing major mergers.

At larger radii (R/r₂₀₀ = 0.2), the stellar-light contours are a median |$0.02\, r_{200}$| (with percentiles 0.01–|$0.03\, r_{200}$|⁠) apart from the DM structure, whereas the isodensity contours of GCs differ by a median |$0.04\, r_{200}$| (with percentiles 0.03–|$0.05\, r_{200}$|⁠). This implies that GC populations are better tracers of the matter distribution within the inner halo, whereas stars are more accurate on average at larger radii. However, from a practical perspective, observational limitations must also be considered. GCs are easily detected at large radii and colour and shape cuts with space-based imaging can eliminate much of the contamination, whereas diffuse stellar light is much more difficult to measure over large areas at low surface brightness levels (e.g. Galactic cirrus can severely contaminate measurements at these low surface brightness levels; Rudick et al. 2010; Mihos et al. 2017). Together, it implies that the distribution of GC populations can yield the distribution of DM more accurately.

To examine why stars and GCs track the matter distribution of their host halo, we explore the radial profiles of the fractions of accreted stars and GCs. We show these together with the profile of in situ objects in Fig. 7. At |$0.05\, r_{200}$|⁠, the fraction of accreted GCs is larger than that of stars across the three galaxy mass bins, whereas both fractions have comparable values at |$0.2\, r_{200}$|⁠. At the radial distances probed, the accreted fractions are larger than their in situ counterparts in galaxies more massive than log₁₀(M_⋆/M_⊙) ≥ 10.5. The accreted stars and GCs dominate the stellar and GC populations, and thus, they are the strongest tracers of the DM halo.

In the lowest galaxy mass bin, the fractions of in situ and accreted objects become comparable at the radial distances probed. The populations of GCs are similar at |$0.05\, r_{200}$|⁠, and the fractions of stars have similar values at |$0.2\, r_{200}$|⁠. The signal from the accreted component is thus diluted and there is more scatter in the MHDs with the DM isodensity contours.

As a proof of concept for this methodology, in the next subsections, we test whether our results are affected by three different aspects: the presence of satellites in the images, the stellar surface brightness limit, and the absolute magnitude cut applied to the GC populations.

4.2 The effect of the presence of satellites

We evaluate the effect of satellite galaxies on the computed MHDs. If a halo hosts a prominent satellite (e.g. it is undergoing a major merger), the isodensity contours drawn at small radii might be more prone to pick up the satellite rather than the central galaxy in the different maps. In that case, the longest contour used to calculate the MHD can be drawn around two different objects (the central and the massive satellite) in the different tracers. Hence, the measured MHD would be the distance between the two galaxies, and it would be larger than expected (red regions in Fig. 6). Our halo FOF000 is an excellent example of this issue (see Fig. 5).

To explore this issue, we prepare new smoothed mock images considering all the particles in the main halo and discard those that are locked in satellite galaxies. This selection retains the populations that are bound to the central galaxy, as well as the diffuse objects that lie in the main halo as a result of its accretion history. We then compute the MHDs between the isodensity contours in the DM and the observational maps at R/r₂₀₀ = 0.2, and we compare them against the values obtained from the fiducial images. We show this comparison in Fig. 8 as a function of the mass fraction in satellite galaxies within the halo.

Examining haloes with increasing stellar mass locked in satellite galaxies, there is an increasing trend: the contours drawn in the images with satellites are increasingly further apart from those without satellites in them. This trend is especially noticeable in haloes with more than 20 per cent of their stellar mass in satellite galaxies.

In the vast majority of cases, we find that the presence of substructure in our images has an effect smaller than |$0.01\, r_{200}$|⁠. There is only one instance in which the distances measured from the images with and without subhaloes differ significantly (i.e. by more than 10 per cent of r₂₀₀). This halo has more than 50 per cent of its stellar mass in satellite galaxies, possibly indicating the presence of one or more massive satellites within the halo. Thus, this methodology can be applied to observational data without the need of removing satellite galaxies from the images – unless the halo is undergoing a major merger or it contains a large number of satellites.

4.3 Modifying the stellar surface brightness limit

We examine the effect of the surface brightness limit of the stellar maps on the measured MHDs. For this, we produce smoothed stellar surface brightness images at different limiting depths. We consider the stellar surface brightness limits of SDSS, (⁠|$\mu _{r-{\rm band}} = 26.7~\rm mag\, arcsec^{-2}$|⁠; York et al. 2000), the Euclid Wide Survey, (⁠|$\mu _{r-{\rm band}} = 29.4~\rm mag\, arcsec^{-2}$|⁠; Euclid Collaboration 2022), and the limit after 10 yr of full survey integration from the Rubin r–band, (⁠|$\mu _{r-{\rm band}} = 30.6~\rm mag\, arcsec^{-2}$|⁠; Ivezić et al. 2019), in addition to our fiducial limit of |$\mu _{r-{\rm band}} = 28~\rm mag\, arcsec^{-2}$|⁠.

Deepening the stellar surface brightness limit allows us to measure the stellar distribution at larger radii for more haloes. Thus, we calculate the MHD at three radial distances that probe the outer regions of the halo: |$0.2\, r_{200}$|⁠, |$0.4\, r_{200}$| and |$0.7\, r_{200}$|⁠, which we show in Fig. 9. At a radius of |$0.4\, r_{200}$|⁠, there are 57 haloes in our fiducial images (⁠|$\mu _{r-{\rm band}}=28~\rm mag\, arcsec^{-2}$|⁠) for which we can measure the MHD, whereas we can do so for the entire sample of 117 haloes when lowering the cut to the end-of-mission Rubin limit. Pushing outwards to a distance of |$0.7\, r_{200}$|⁠, there are 3 and 33 haloes for which we can measure their MHDs in our fiducial and Rubin images, respectively. We summarize the number of galaxies for each observational limit in Table 2.

Interestingly, deeper observations barely increase the precision of the MHD compared to shallower images of the same halo. Although deeper data would allow us to apply this methodology to lower-mass haloes and further out in the galaxies, diffuse stellar light is very difficult to measure over large areas (e.g. Mihos et al. 2017). Thus, present-day shallower measurements of diffuse stellar light are sufficient to reconstruct the inner distribution of DM in haloes with central galaxies more massive than M_⋆ ≥ 10¹⁰ M_⊙.

4.4 Modifying the GC absolute magnitude limit

We examine the effect of the absolute magnitude limit in the F475W band applied to GC populations on the measured MHDs. The fiducial limit of M_F475W = −6.2 corresponds to the peak of the GC mass function (∼10⁵ M_⊙) for a |$8~\mbox{${\rm Gyr}$}$| old cluster of solar metallicity. Thus, increasing and decreasing this value, we select GC populations that contain, respectively, lower and higher-mass clusters than the peak of the mass function. Regardless of the absolute magnitude limit used, we also select GCs based on colour combinations with a UV filter.

Using these different absolute magnitude limits, we produce the corresponding smoothed projected number density maps of the GC populations. We compute the MHDs between them and the DM contours, which we show in Fig. 10.

Observing fainter GC populations (M_F475W < −4.2) allows us to compute the MHDs to further radial distances for more systems. Using the faintest GC samples, we can measure the MHD to the DM distribution for 114 and 57 haloes at |$0.2\, r_{200}$| and |$0.4\, r_{200}$|⁠, respectively. In comparison, there are 13 and 3 haloes, respectively, for which we can do the same analysis with only the brightest GC populations. Additionally, we find that observing fainter GC populations leads to a more accurate recovery of the DM distribution (i.e. smaller MHDs). This is in contrast with the result obtained varying the stellar surface brightness limit, in which deeper observations only added information towards the outskirts of the haloes.

At |$0.4\, r_{200}$| in the faintest GC populations, there is an increasing trend of MHDs towards lower-mass galaxies. This trend is related to the mass range of the satellite galaxies accreted by these haloes; lower-mass satellites host fewer GCs (e.g. Harris et al. 2017) and therefore contribute fewer accreted GCs to populate the outskirts during their host galaxy accretion. Lastly, we find that the MHDs can only be computed at a distance of |$0.7\, r_{200}$| for the faintest GC populations in the six most massive haloes. For those haloes, the measured MHDs from the GC contours are |$\sim 0.1\, r_{200}$| away from the DM distribution. These results suggest that GC populations can provide accurate estimates of the shape of the matter distribution out to |$0.4\, r_{200}$| of their host DM halo for haloes with central galaxy stellar masses of M_⋆ ≳ 10¹⁰ M_⊙.

Dynamical modelling of large samples of GCs has been used to recover the distribution of DM around galaxies (e.g. Bellstedt et al. 2018). However, measuring the velocities of large samples of tracers can quickly become prohibitive. In the absence of velocities, our method opens up the possibility to recover the profile of DM in large photometric surveys.

5 RECOVERING THE DM SURFACE DENSITY PROFILE

We examine whether the spatial distributions of stars or GCs can be calibrated to recover the DM radial profile. We explore this calibration using physical (i.e. surface densities) and observational (i.e. surface brightness and number density) quantities. We first determine the calibration of the physical quantities. To do this, we calculate the dimensionless ratios of the stellar-to-DM surface density profiles and GCs-to-DM surface density profiles. We show these ratios in Fig. 11.

There are very tight relations between both the stellar and GC surface density profiles and those of the DM. To quantify these relationships, we fit linear regressions of the form

where y = Σ_stars and y = Σ_GCs for stars and GCs, respectively. We perform the fits for the radial range 0.05 ≤ R/r₂₀₀ ≤ 0.4, and we provide the best-fitting parameters in Table 3, alongside the standard deviation around the linear fit.

For both the stellar and GC profiles, the slopes of the fits become shallower for more massive haloes. This trend is driven by the increasingly larger fraction of accreted stars and GCs in more massive haloes, as this material is deposited further out. Additionally, massive galaxies accrete more objects that have high stellar-to-halo mass ratios, so there is more accreted DM in the outskirts.

Across the three galaxy mass bins, the relation between the stellar surface density profiles and those of the DM is tighter (σ = 0.2–0.3 dex) than the relation between the surface density profiles of GCs and those of the DM (σ = 0.3–0.5 dex). In the most massive galaxy bin, the scatter around the GC-to-DM relation is similar to that of stars, and it increases for decreasing galaxy mass. This indicates that the stellar distribution would yield a more accurate reconstruction of the DM surface density profile across our galaxy mass range. However, observational constraints that limit our ability to measure large-scale low surface brightness structures (e.g. Galactic cirrus; Mihos et al. 2017) can hinder their practical use. Since GC populations do not suffer from these observational constraints, they would be a better tracer overall for the reconstruction of the DM profile.

A similar relation between the stellar and the total mass surface density radial profiles is found by Alonso Asensio et al. (2020) (red dashed line in Fig. 11). Compared to our fits, their linear relationship is shallower than the one obtained in the most massive galaxy bin. The authors derive the relation from simulated clusters of galaxies of masses M₂₀₀ = 6 × 10¹³–2 × 10¹⁵ M_⊙ (Cluster-eagle project; Bahé et al. 2017; Barnes et al. 2017), which are a factor of 1.2–40 times larger than our most massive halo. Thus, their shallower slope confirms the trend observed among a sample of haloes: as haloes built larger fractions of their mass via accretion, more material is deposited in the outer regions.

Similarly, we can do a qualitative comparison to the stellar-to-DM and GC-to-DM ratios from the galaxy cluster Abell 1689 (fig. 10 of Alamo-Martínez et al. 2013, purple lines in Fig. 11). This galaxy cluster has a mass of 10¹⁵ M_⊙ and a factor of 20 times larger than our most massive halo. Both the observational ratios of the stellar and GC profiles are shallower than our simulated profiles for the most massive galaxy bin, and the stellar-to-DM slope is comparable to the fit from Alonso Asensio et al. (2020). This is further evidence that the calibration to recover the DM surface density profile is sensitive to the mass of the halo.

Using these physical calibrations requires being able to determine the average stellar and GC surface density within a radial bin. Our goal is to obtain calibrations that can be readily applied to observations. Thus, we have repeated this analysis using mock observational quantities (i.e. stellar surface brightness and GC number density). There are similarly tight relations between these observational quantities and the DM surface density profile. We provide the best-fitting coefficients for linear fits in Table 3, which is adapted to be

in the case of the stellar light. The small scatter around the observational linear fits suggests that either the stellar or GC number density profiles are sufficient to recover the DM profile with high accuracy, i.e. within factor of ∼2 for stars, and ∼3 for GCs. Hence, it opens a novel avenue to recover the DM surface density profile for a large sample of galaxies with upcoming surveys such as Euclid or the Rubin observatory.

6 PROJECTED MORPHOLOGY OF DM, STARS AND GCS

We provide quantitative descriptions of the overall morphology of the DM, stellar, and GC distributions. To do so, we model the projected spatial distributions (as described by the particles) with ellipses following the methods outlined by Thob et al. (2019) and Hill et al. (2021, 2022).⁷

We use the ellipticity parameter to describe the projected spatial distributions of DM, stars, and GCs. The ellipticity is defined as ϵ = (a − b)/(a + b), where a and b are the moduli of the semimajor and semiminor axes of the ellipses, respectively. Given a particle distribution, the eigenvalues of a matrix describing its 2D mass distribution (λ_i for i = 1, 2) describe the modulus of the major and minor axes of the corresponding ellipse: |$a = \sqrt{\lambda _1}$| and |$b = \sqrt{\lambda _2}$| with λ₁ > λ₂. Similarly, the associated eigenvectors (e_i for i = 1, 2) correspond to the semimajor and semiminor axes of the ellipse.

We use the reduced mass distribution tensor (see also Davis et al. 1985; Dubinski & Carlberg 1991; Cole & Lacey 1996; Bett 2012; Schneider, Frenk & Cole 2012) in an iterative scheme.⁸ By considering the reduced version of the mass distribution tensor, we suppress possible strong influences from features in the outer regions of the distributions, and the iterative approach enables the scheme to adapt to particle distributions that deviate from the initial selection. The reduced inertia tensor can be described as,

where |$\tilde{r}_{p}$| is the elliptical radius and r_p, i is the i-th component of the coordinate vector of the particle p with respect to the centre of the halo.⁹ In the first iteration, the scheme considers all the particles enclosed within an initial circular aperture. For consistency with the previous sections, we quantify the morphologies within two initial apertures: |$0.05\, r_{200}$| and |$0.2\, r_{200}$|⁠. This provides a first estimate of the axes lengths a and b. In the next iterations, the scheme only considers those particles that satisfy that their elliptical radius is

where |$\tilde{a}$| and |$\tilde{b}$| are re-scaled axes lengths such that |$\tilde{a}=a\, r_{p}\, (ab)^{-1/2}$|⁠. The scheme continues iterating until the fractional change in b/a falls below 1 per cent.

The iterative scheme requires at least 100 particles within the initial circular aperture in order to converge. This constraint limits the number of haloes for which we can measure the ellipticity and orientation of GC populations. Our sample of 117 haloes is thus reduced to 88 and 93 haloes, respectively, for the |$0.05\, r_{200}$| and |$0.2\, r_{200}$| initial apertures. This requirement does not affect the number of haloes for which we can measure the morphology of DM or stars.

Finally, we define the orientation of each component as the unit vector parallel to the major axis of the best-fitting ellipse (e₁^x with |$x=\,$|DM, stars or GCs). We determine the relative alignment between components x and y as the angle between these units vectors,

The misalignment angle is invariant under a rotation of 180 degrees, and it is thus confined between [0,90] degrees. It has been shown that misalignment angles between projected distributions tend to be smaller (i.e. more aligned) at all radii and halo masses than when measured in three dimensions (Tenneti et al. 2014; Velliscig et al. 2015a).

We show the distribution of the ellipticities describing the spatial distributions of DM, stars, and GCs in Fig. 12, and we summarize the median and the 25–75th percentiles in Table 4. Overall, we find that the DM distributions are close to being circular, whereas the stellar and GC distributions are more flattened. For the smaller initial aperture (⁠|$0.05\, r_{200}$|⁠), the median ellipticity of DM is 0.09, whereas the median ellipticities of stars and GCs are 0.23 and 0.20, respectively. The distribution of ellipticities of the stellar and GC distributions are very similar, with the 25–75th percentiles spanning the same range of values, ∼0.13–0.35. The flattened stellar and GC distributions might be explained by the flatness of galaxy discs, which reside within |$r\lt 0.05\, r_{200}$|⁠, as a result of the dissipative collapse of the gas.

The distribution of ellipticities of the stellar and GC populations barely change when considering the outer initial aperture (⁠|$0.2\, r_{200}$|⁠). We relate this to most of their mass already being enclosed within the smaller initial aperture. In contrast, the DM distributions appear slightly more flattened with a median of 0.11. One way to improve the agreement between the DM and the stellar and GC distributions would be to remove the central galaxy from the distribution being fitted. Another method would be to remove the mass weighting in the reduced inertia tensor, as most of the mass in stars and GCs mass resides in the inner region of the halo whereas most of the DM mass sits in the outskirts.

The relative alignment between the semimajor axes of the best-fitting ellipses to the DM, stellar, and GC distributions are shown in Fig. 13. Regardless of the initial aperture, all the distributions are preferentially aligned, with the orientation of ∼70 per cent of the GC systems being less than 20 degrees away from their host DM halo. The DM and stellar distributions show the closest alignment with a median difference of ∼8–9 degrees at the two radial distances considered (also see fig. 13 of Hill et al. 2021). GCs are more aligned with the stellar than the DM distribution of their host halo, although the median misalignment angles among all components only differ by ∼2–3 degrees. Previous studies have already found that galaxies in the eagle simulations are well aligned with the local mass distribution but are often misaligned with respect to the global halo. This indicates that the stellar component follows the local DM distribution, which is the dominant matter component, and that the structure of the DM halo changes from the inner to the outer halo (Velliscig et al. 2015a; Hill et al. 2021).

It is worth noting that some severe misalignments are simply cases where the two components are actually aligned, just not along the same axes. This happens most often in very prolate haloes (e.g. see fig. 7 of Hill et al. 2021). We have also examined the ellipticities and the misalignment for different galaxy stellar mass bins, but no trends could be identified. Thus, we find that characterizing the orientation of the spatial distribution of stars and GCs is a very accurate probe of the orientation of the DM halo.

7 CONCLUSIONS

The aim of this study is to evaluate whether the diffuse stellar and GC populations correlate with the overall structure of the matter distribution of their host halo and thus, whether the DM distribution can accurately be recovered using observations of GCs and diffuse stellar light. For this, we use the |$(34.4~\rm cMpc)^3$| periodic volume from the E-MOSAICS project (Pfeffer et al. 2018; Kruijssen et al. 2019a), which is a unique data set to study stellar cluster populations and their host galaxies after a Hubble time of evolution. We select the 117 DM haloes containing central galaxies more massive than M_⋆ > 10¹⁰ M_⊙, which corresponds to halo masses M₂₀₀ ≳ 4 × 10¹¹ M_⊙. For each of these haloes, we create smoothed surface density maps of DM and stellar surface brightness maps and number density maps for the GC populations. The stellar light shows the highest degree of correlation to the DM distribution across our sample, whereas GC populations are increasingly uncorrelated to both the stellar and the DM distributions in lower mass systems (Fig. 3).

To compare the distributions of DM, stars and GCs, we identify isodensity contours for a given set of radial distances (Fig. 4). In haloes with prominent satellites (e.g. ongoing major mergers), the algorithm can identify different structures in different tracers (i.e. the central and the satellite), and thus the comparison of the contours becomes more challenging. In most of our haloes, the structures identified in the diffuse stellar light and in the GC populations qualitatively agree with those in the DM maps (Fig. 5).

We measure the MHD between pairs of mock observational maps (i.e. stars and GCs) and the DM maps (Fig. 6). This metric provides a quantitative measure of how close together two spatial distributions are. For our fiducial limits (⁠|$\mu _{r-{\rm band}} \lt 28~\rm mag\, arcsec^{-2}$| in the stellar maps and M_F475W < −6.2 for the GCs), the GC populations follow closely the matter distribution in the inner region of the halo (⁠|$0.05\, r_{200}$|⁠), whereas the diffuse stellar light is more accurate at larger radii. However, observing GC populations is less challenging than mapping the diffuse stellar light, and thus it can be a more effective way of probing their host halo.

The accuracy of each observational tracer is related to the dominance of accreted objects (Fig. 7), i.e. the material deposited during the assembly of the host halo contains the information about the matter distribution. At the two radial distances considered, the fraction of accreted GCs is greater than or similar to that of stars. In low mass haloes (10 < log₁₀(M_⋆/M_⊙) < 10.5), the accreted and in situ fractions become comparable, and thus the signal tracing the matter distribution dilutes and it is more difficult to recover.

As a proof of concept, we examine how much the measured MHDs are affected by three different aspects. First, we consider the presence of substructure within each halo. Removing the galactic substructure from our maps has an effect smaller than |$0.01\, r_{200}$| in the measured MHDs relative to the DM contours (Fig. 8). This is a promising result for future applications of this analysis to observational data since removing substructure in observations is difficult.

Secondly, we consider the influence of the stellar surface brightness limit (Fig. 9). We consider the limits from SDSS, the Euclid Wide Survey and the 10-yr integration limit of the Rubin observatory and compare the results to our fiducial limit. Deeper observations do not lead to smaller MHDs but allow us to apply this analysis to lower-mass haloes and to larger radial distances within the halo.

Thirdly, we examine the role of the GC absolute magnitude limit by selecting populations that reach above and below the fiducial value (i.e. the peak of the GC luminosity function, Fig. 10). Observing fainter GCs allows the measurement of the MHDs to be extended to lower-mass haloes and to larger radial distances within the halo. Additionally, we find that the contours identified in the deeper GC observations are closer to the DM isodensity contours than those identified in brighter populations.

We examine whether the spatial distributions of stars and GCs can be calibrated to recover the DM surface density radial profile (Fig. 11). More massive haloes have shallower ratios of stars to DM and GCs to DM, which is related to having built up a larger fraction of their mass from accreted satellites that deposit stars and GCs further out in the halo. We also find very tight relations between the stellar-to-DM surface density profiles, and the GC-to-DM profiles (the scatter around the fit is σ ∼ 0.2–0.5 dex). Similar tight relations are found for the mock observational quantities (i.e. stellar surface brightness and GC number density), suggesting that the DM surface density profile can be recovered within a factor of ∼2 and ∼3 from deep observations of stars and GCs, respectively.

We quantify the projected ellipticity and orientation of DM, stars, and GCs by modelling the spatial distribution of their particles using ellipses (see also e.g. Tenneti et al. 2014; Velliscig et al. 2015a; Thob et al. 2019; Hill et al. 2021). DM shows the most circular distributions, whereas those of stars and GCs are more flattened (Fig. 12). All distributions are preferentially aligned (Fig. 13), with median misalignment values of ∼10 degrees and the 75th percentiles being less than 30 degrees (Table 4). The stellar distribution shows a closer alignment to DM than GCs, but the difference is only ∼2 degrees. Stars have previously been shown to follow closely the local matter distribution (Velliscig et al. 2015a; Hill et al. 2021), and our results suggest GC populations are also effective tracers. Thus, by characterizing the orientation of diffuse stellar light and GC populations, we can determine the orientation of their host DM halo within a few degrees.

Our results suggest that the host DM haloes of observed galaxies can be mapped out to distances of nearly half the virial radius using the projected distribution of GCs, and to even further distances using very deep stellar observations. The accuracy when using GCs is highest at |$0.05\, r_{200}$|⁠, and for galaxies with the most abundant GC systems (including massive ellipticals). In the local Universe (⁠|$D \lt 50~\mbox{${\rm Mpc}$}$|⁠), observations of GC populations do not require deep observations, in contrast with diffuse stellar haloes. Therefore, maps of GC systems provide a more efficient way to trace the DM distribution. In future works, we will compare these results to observations for which studies of GC populations, diffuse stellar light, and mass distribution exist. Environments of particular interest include the galaxy cluster Abell 1689 (Alamo-Martínez et al. 2013), Fornax cluster (e.g. Chaturvedi et al. 2022), and, recently observed with the James Webb Space Telescope, the galaxy cluster SMACS0723 (Diego et al. 2023). Upcoming deep and wide surveys with the Euclid, Roman, and Rubin observatories will eventually allow the reconstruction of the host DM halo of very large samples of galaxies using this method. Comparing those results to DM maps of hundreds of galaxy clusters obtained using strong lensing analyses in next-generation surveys could yield new constraints on the nature of DM (e.g. Despali et al. 2022).

ACKNOWLEDGEMENTS

MRC thanks Bill Harris for valuable discussions on mimicking observations and Alex Hill for providing his code to measure the shape of particle distributions. MRC gratefully acknowledges the Canadian Institute for Theoretical Astrophysics (CITA) National Fellowship for partial support, and this work was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC). STG 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). JLP is supported by the Australian government through the Australian Research Council’s Discovery Projects funding scheme (DP200102574). AD is supported by a Royal Society University Research Fellowship. AD also acknowledges support from the Leverhulme Trust and the Science and Technology Facilities Council (STFC; grant numbers ST/P000541/1, ST/T000244/1). RAC is a Royal Society University Research Fellow. JMDK gratefully acknowledges funding from the German Research Foundation (DFG) in the form of an Emmy Noether Research Group (grant number KR4801/1-1). COOL Research DAO is a Decentralised Autonomous Organisation supporting research in astrophysics aimed at uncovering our cosmic origins.

This work used the DiRAC Data Centric system at Durham University, operated by the Institute for Computational Cosmology 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 grants ST/H008519/1 and ST/K00087X/1, STFC DiRAC Operations grant ST/K003267/1 and Durham University. DiRAC is part of the National E-Infrastructure. The work also made use of high performance computing facilities at Liverpool John Moores University, partly funded by the Royal Society and LJMU’s Faculty of Engineering and Technology.

Software: This work made use of the following python packages: h5py (Collette et al. 2021), jupyter notebooks (Kluyver et al. 2016), numpy (Harris et al. 2020a), pynbody (Pontzen et al. 2013), and scipy (Jones et al. 2001), and all figures have been produced with the library matplotlib (Hunter 2007). The comparison to observational data was done more reliably with the help of the WEBPLOT-DIGITIZER¹⁰ webtool.

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

Footnotes

References

APPENDIX A: INFLUENCE OF THE KERNEL SIZE

In this appendix, we explore the influence of the size of the Gaussian kernel on the structures that can be identified from the smoothed maps. We consider five different kernel sizes, |$\sigma =0,0.01,0.02,0.05,0.10\, r_{200}$|⁠, with the former implying that no smoothing is performed. We show the resulting GC number density maps of FOF000 in Fig. A1. As we increase the size of the kernel, we find that the small-scale noise is smoothed out. Once the kernel is larger than the effective radius of the central galaxy (⁠|$\sim 0.015\, r_{200}$|⁠; Kravtsov 2013), the galactic-scale signal is washed out and the isodensity contours are very smooth. For these reasons, we decide to use a kernel size of |$\sigma =0.02\, r_{200}$| in our main analysis. Using a kernel of size |$\sigma =0.02\, r_{200}$|⁠, we smooth out the small-scale noise while retaining the galactic-scale perturbations in the contours.