methods: numerical – galaxies: clusters: general – galaxies: evolution

1 Introduction↩︎

The light produced by the diffuse ensemble of stars that permeate the intergalactic space in a galaxy cluster – not belonging to any particular galaxy but still bound to the cluster potential as a whole – is termed intracluster light (ICL; see [1], [2], [3] for relevant reviews). The existence of this diffuse stellar component was first theorised and observed by [4], [5]. Though much of the ICL in a typical cluster will be centrally concentrated – composing an extended diffuse stellar envelope around the brightest cluster galaxy (BCG; e.g. [6], [7], [8]) – this low surface-brightness cluster component (\(\mu_{r} \gtrsim 26\,\textrm{mag arcsec}^{-2}\); [9]) can extend to the far fringes of a cluster, with observational studies tracing the ICL out to \(\gtrsim1\) Mpc ([10], [11], [12]).

Typically over 10 per cent of the stellar mass in a cluster may be associated with the ICL, though with substantial scatter between individual clusters (e.g. [13], [14], [15]). Significant scatter also occurs between different definitions and methodologies for extracting the ICL (e.g. [16], [17], [18], [9]). Likewise unresolved is the question of how this typical ICL mass fraction evolves with redshift, though a mounting body of both observational (e.g. [19], [20], but see also e.g. [21]) and theoretical work (e.g. [22], [23], [24], but see also e.g. [25]) suggests this diffuse stellar component should remain significant to \(z\sim1\) (and plausibly to \(z\gtrsim2\); e.g. [26], [24], [27]).

As hierarchically formed structures, galaxy clusters are thought to grow through a protracted procession of accretion and merger events ([28], [29]) and so the ICL – as the stellar debris generated by these and related processes – is expected to encode a record of each cluster’s merger history as well as its past and present dynamical state (e.g. [30], [24], [31], and references therein). Given that the ICL can span the entire cluster halo, with observational detections made even beyond the splashback radius [10], and that these intergalactic stars can be considered collisionless particles with motions governed only by the overall cluster potential, it has also been proposed that the ICL could be used as a visible tracer of the total mass distribution in galaxy clusters (e.g. [32], [33], [34]; see also [35]). However, exploiting observations of the ICL for these purposes requires both a detailed understanding of the processes that contribute stars to the ICL and the objects from which these constituent stars are sourced ([36]).

The currently prevailing paradigm is that the bulk of ICL stars should originate from cluster satellite galaxies, chiefly through a combination of tidal stripping (e.g. [37], [38], [39]) and ejection by violent relaxation processes during galaxy mergers (chiefly between the BCG and massive satellites; e.g. [40], [41], [42]). There is an emerging consensus that the main progenitor objects from which these ICL stars are drawn should be massive galaxies (stellar mass, \(M_{*} \gtrsim10^{10.5}\) M). A developing body of theoretical work supports this perspective (e.g. [43], [44], [45], [46], [47], [39], [48], [15], [49]; but see also e.g. [50], [51]). The alternative, that ICL stars are predominately sourced from less massive satellites, requires in some combination either significant modifications to the galaxy stellar mass function or a considerable suppression of stripping rates at higher satellite masses compared to current expectations ([3], [39], [48], [36]).

Observational studies remain more divided. Though many agree with the prevailing view from theoretical studies, that the ICL should be largely sourced from more massive galaxies (e.g. [52], [53], [54], [55], [56]), some instead suggest that less massive galaxies (\(M_{*}\lesssim10^{10}\) M) are the main progenitors of ICL stars, on the basis of e.g. observed ICL colours and metallicities (e.g. [57], [58], [59], [7]). However, among this latter body of work there is an acknowledged degeneracy which might ease this tension: declining stellar metallicity with increasing galactocentric radius (e.g. [60], [61], [62], [63]), allowing a population of stars tidally stripped from the outskirts of a more massive galaxy to have a typical metallicity mimicking that characteristic of a less massive galaxy ([59], [7]).

Numerous prior studies have already investigated the ICL in individual hydrodynamical simulations, but the findings of these studies can significantly diverge concerning various key ICL properties and statistics such as the typical ICL stellar mass fraction (e.g. see table 1 in [18]). Some disagreement is generally expected between different simulations, from (for example) the differing numerical methodologies implemented to model gravity and baryonic physics ([64], [65], [66], [67], [68]) as well as the simulation resolution achieved ([39], [69], [70]). Conceivably even the broad simulation strategy adopted, such as uniform-resolution vs “zoom-in” simulation, may have an effect (e.g. [71], [72]). However, an understanding of these essential differences is obscured by these studies also generally employing bespoke and often incompatible methods to identify the ICL. In order that any differences can be understood and reconciled separately from the impact of differing ICL extraction methodologies, comparison studies employing a uniform approach to identifying the ICL are necessary.

Here, and in a forthcoming paper (Brown et al. in prep.), we investigate the assembly of the ICL in a number of simulated galaxy clusters drawn from four state-of-the-art hydrodynamical simulations (\(z\approx0\) halo masses \(\in[10^{14},10^{15}]\) M; \(10-16\) clusters per simulation) which leverage a diverse selection of simulation practices and numerical methods, upon which we have imposed a uniform ICL definition. In this paper we consider the division of stellar mass between galaxies and the ICL in these simulated clusters. We do this both at \(z\approx0\) and in the progenitors of these structures back to \(z\sim3\), and probe how our findings are influenced by the specific ICL definition implemented. We also investigate the relative contributions to the ICL from satellite galaxies with differing stellar masses preceding cluster infall and so assess the main progenitor galaxies of ICL stars, extending the work described in [48]. We leave a comprehensive evaluation of the ICL contribution from assembly channels beyond the liberation of satellite galaxy stars as well as an investigation into various observed ICL trends with cluster-centric radius and an interrogation of plausible origins for these trends to the forthcoming companion paper.

This paper is organized as follows. In Section 2 we introduce and briefly detail the four simulations (2.1), our fiducial homogenized ICL definition (2.2), how we associate \(z\approx0\) ICL stars to progenitor galaxies (2.3), and the cluster samples we employ (2.4). In Section 3 we present our results for the division of stellar mass between galaxies and the ICL at \(z\approx0\) (3.1), and how this has evolved since \(z\sim3\) (3.2), as well as probe how our findings might be influenced by the imposed BCG-ICL separation (3.3). In Section 4 we investigate the ICL contribution from a cluster’s own satellite galaxies (4.1) as well as the differing contributions made to the ICL by satellite galaxies with different masses (4.2). In Section 5 we discuss the physical origins of noted trends (5.1), apparent inter-simulation discrepancies (5.2), as well as important caveats and limitations of our analysis (5.3). We conclude by summarising the main results of our investigation in Section 6. For the following analysis, we adopt the respective cosmologies of the simulations being considered; distances should be assumed proper (not comoving); and presented quantities that depend on the dimensionless Hubble parameter have been evaluated using the value from the corresponding cosmology and do not contain additional factors of “little \(h\)” unless explicitly indicated otherwise.

2 sec:Methods↩︎

2.1 Simulations↩︎

We use data from four cosmological hydrodynamical simulation projects: the Horizon-AGN simulation (hereafter Hz-AGN); the TNG100-1 simulation of the IllustrisTNG project (hereafter TNG100); the Gizmo-Simba 7K galaxy cluster simulations of the The Three Hundred project (hereafter The300 GS-7K); and the Hydrangea galaxy cluster simulations. We summarise below the properties, models, and methods employed by each of the simulations. Some key properties are also presented for comparison in Table 1. Further details can be found in the respective introductory papers of each simulation.

2.1.1 Horizon-AGN↩︎

The Horizon-AGN simulation² ([73], [74], [75]) is a uniform-resolution cosmological-volume hydrodynamical simulation employing the Adaptive Mesh Refinement (AMR) Eulerian hydrodynamics code Ramses [76], with a cubic simulation volume of side-length \(100\,h^{-1}\,\)Mpc (comoving). A \(\Lambda\)CDM cosmology consistent with the 7-year Wilkinson Microwave Anisotropy Probe data (WMAP7; [77]) is adopted with \((\Omega_{m,0}, \Omega_{\Lambda,0}, \Omega_{b,0}, \sigma_{8}, n_{s}, h) = (0.272, 0.728, 0.045, 0.810, 0.967, 0.704)\).

The simulation follows \(1024^3\) dark matter (DM) particles (mass, \(m_{\textrm{DM}}\sim8 \times 10^{7}\,\)M) and an initially uniform \(1024^3\) cell gas grid (initial gas resolution \(m_{\textrm{gas,ini}}\sim1\times10^{7}\,\)M) dynamically refined according to a quasi-Lagrangian criterion (based on either the total baryonic or DM mass in a cell exceeding eight times the mass of a DM particle) down to a minimum cell size of \(1\) kpc after seven levels of refinement. Star formation occurs beyond a threshold gas density of \(n_{\textrm{H}}=0.1\,\textrm{cm}^{-3}\) following the Kennicutt–Schmidt law ([78]) with a constant star efficiency of 2 per cent per free-fall time. Stellar particles have a mass of \(m_{*}\sim3 \times 10^{6}\,\)M. Gas heating by a uniform UV background begins after \(z = 10\); H and He cooling with a contribution from metals using a [79] model can cool gas to \(10^{4}\,\textrm{K}\). Continuous stellar feedback that includes momentum, mechanical energy and metals from stellar winds and supernovae is included, in addition to dual-mode AGN feedback switched according to gas accretion rate following the [80] model.

2.1.2 TNG100↩︎

The TNG100 simulation of the IllustrisTNG project³ ([81]; see also [82], [83], [84], [85], [86]) is a uniform-resolution cosmological-volume (magneto-)hydrodynamical simulation employing the moving-mesh (magneto-)hydrodynamics code Arepo ([87], [88], [89]) with a cubic simulation volume of side-length \(75\,h^{-1}\,\)Mpc (comoving). From the three flagship IllustrisTNG volumes (TNG50, TNG100, and TNG300) we adopt the intermediate TNG100 simulation for this study as a compromise between resolution and volume. A \(\Lambda\)CDM cosmology consistent with the [90] results is assumed with \((\Omega_{m,0}, \Omega_{\Lambda,0}, \Omega_{b,0}, \sigma_{8}, n_{s}, h) = (0.309, 0.691, 0.049, 0.816, 0.967, 0.677)\).

TNG100 follows \(1820^3\) DM particles (\(m_{\textrm{DM}}\approx7.5\times10^6\,\)M) and initially has \(1820^3\) gas cells (minimum cell radius \(14\) pc; target gas cell mass, \(m_{\textrm{gas}}\approx1.4\times10^6\,\)M). Stochastic star formation occurs above a threshold gas density of \(n_{\textrm{H}}=0.1\,\textrm{cm}^{-3}\) following the [91] model at rates that empirically reproduce the Kennicutt-Schmidt law [78], with the typical stellar particle formation mass approximately matching the target gas cell mass (i.e. \(m_{*}\sim1\times10^6\,\)M). The Plummer equivalent gravitational softening of collisionless particles is \(0.5\,h^{-1}\,\)kpc at \(z\leq1\) and \(h^{-1}\,\)comoving kpc at \(z\geq1\); adaptive gravitational softening is used for gas, set at \(2.5\) times the effective cell radius with an enforced minimum of \(0.125\,h^{-1}\,\)kpc (comoving). Radiative gas cooling (primordial and metal-line) occurs in the presence of a redshift-dependent UV background switched on at \(z=6\) ([92]). Stellar feedback with galactic-scale, energy-driven kinetic winds is implemented following the [92] model, and dual-mode AGN feedback switched according to accretion rate is implemented following the [93] model.

2.1.3 The Three Hundred Gizmo-Simba 7K↩︎

The The Three Hundred simulation suite⁴ ([94]) consists of hydrodynamical zoom-in simulations of the 324 most massive haloes (masses \(\in[10^{14.8},10^{15.5}]\,\)M) drawn from the MultiDark Planck 2 cosmological N-body simulation (side length \(1\,h^{-1}\,\textrm{Gpc}\); [95]). In this work we use data from the new GIZMO-SIMBA 7K simulation runs (Cui et al. in prep.), a higher-resolution and otherwise improved follow-up of the original GIZMO-SIMBA (3K) simulations [96]. The meshless finite-mass Lagrangian hydrodynamics solver GIZMO ([97], [98]) is employed, coupled with the SIMBA-C galaxy formation model ([99], [100]). This fork of the SIMBA model [101] is differentiated chiefly by adopting the state-of-the-art Chem5 cosmic chemical enrichment model ([102] and references therein). The simulation regions are each centred on a cluster-scale halo, with the high-resolution region extending at least \(15\,h^{-1}\,\textrm{Mpc}\) from the cluster centre at \(z=0\). A \(\Lambda\)CDM cosmology consistent with the [103] results is assumed with \((\Omega_{m,0}, \Omega_{\Lambda,0}, \Omega_{b,0}, \sigma_{8}, n_{s}, h) = (0.307,0.693,0.048, 0.829, 0.961, 0.678)\).

In the high-resolution regions, DM particles have \(m_\text{DM}\sim2\times10^{8}\,\)M and gas particles \(m_\text{gas}\sim4\times10^{7}\,\)M. Star formation occurs above a threshold gas density of \(n_{\textrm{H}}=0.1\,\textrm{cm}^{-3}\) following a [104] derived model which scales with gas molecular hydrogen content, and stellar particles form with \(m_{*}\sim4\times10^{7}\,\)M. The gravitational softening length for collisionless particles is \(1.25\,h^{-1}\) kpc at \(z\leq1\) and \(2.5\,h^{-1}\) comoving kpc at \(z\geq1\); adaptive gravitational softening is used for gas [97] with an enforced minimum softening length of \(0.1\,h^{-1}\) kpc at \(z\leq1\) and \(0.2\,h^{-1}\) comoving kpc at \(z\geq1\). Radiative gas cooling (including metal cooling) and photoionization is incorporated via the GRACKLE-3.3 library [105]. SIMBA-C largely follows the dual-mode AGN feedback implementation of SIMBA (as described in [101], [106]; see [99] for minor updates) and a recalibration of SIMBA’s star-formation-driven two-phase galactic winds model is also used, but the stellar feedback treatment in SIMBA-C has otherwise been substantially overhauled as compared to SIMBA (see [99] for details).

2.1.4 Hydrangea↩︎

The Hydrangea simulations ([107], [108]) are a suite of 24 hydrodynamical galaxy cluster zoom-in simulations run using a variant of the EAGLE galaxy formation and evolution code (variant “AGN-dT9”; [109]) and form part of the C-EAGLE project [110]. The EAGLE code is a substantial modification of the Gadget-3 Smoothed Particle Hydrodynamics (SPH) code [111]. Each of the simulated regions is selected from an N-body simulation (MACSIS; [112]) with side length \(3.2\) Gpc, and are each centred on a cluster-scale halo (mass \(\in[10^{14.0},10^{15.4}]\,\)M; subject to the additional selection criteria of there being no more massive halo nearby at \(z=0\)). The high-resolution regions extend to at least \(10\times r_{200c}\)⁵ from the centre of the central cluster at \(z=0\). A \(\Lambda\)CDM cosmology consistent with the [103] results is assumed with \((\Omega_{m,0}, \Omega_{\Lambda,0}, \Omega_{b,0}, \sigma_{8}, n_{s}, h) = (0.307,0.693,0.048, 0.829, 0.961, 0.678)\).

In the high-resolution regions, DM particles have \(m_\text{DM}\approx9.8\times10^{6}\,\)M, and gas particles have an initial mass of \(m_\text{gas,ini}\approx1.8\times10^{6}\,\)M. Star formation occurs above a metallicity-dependent gas density threshold of \(n_{\textrm{H}}=(Z/0.002)^{-0.64}\times0.1\) cm\(^{-3}\) [113] following the pressure implementation of the Kennicutt–Schmidt law [114] from [115], and stellar particles have \(m_{*}\sim1\times10^{6}\,\)M. The Plummer-equivalent gravitational softening length is 0.7 kpc at \(z<2.8\). Additions atop the base hydrodynamics scheme of Gadget-3 include radiative cooling, photoheating, and reionization (using [116] models); mass and metal enrichment of gas due to stellar outflows (based on [117]); and thermal feedback from both star formation and supermassive black holes ([118], [119], [109]). For further detail on the Eagle code and its “ANARCHY” hydrodynamics scheme update see [109] and [120].

2.2 Homogenized ICL definition↩︎

2.2.1 Distinguishing galactic and intergalactic stars↩︎

The four simulations we employ for this study by default adopt different codes for structure identification. Hz-AGN uses the AdaptaHOP structure finder ([121], [122]). TNG100 uses Subfind ([123], [124]). The300 GS-7K uses the Amiga Halo Finder (AHF; [125], [126]). We employ the Cantor structure catalogue (Bahé et al. in prep.) for Hydrangea. We omit detailed general descriptions of each structure finder for concisenesses, referring the interested reader to the introductory papers of each code (see also [127], [128] for a comparison of the different approaches).

Subfind, AHF, and Cantor do not innately separate “galaxies” from the DM (sub-)structures hosting them. The stars and gas of a galaxy, the DM of its host (sub-)halo, as well as any circumgalactic medium or diffuse stellar envelope within that (sub-)halo are all combined together into one amalgamated structure: a “sub-halo”⁶. In the context of ICL studies, this means that there is no separation between the BCG and the ICL of a cluster imposed by these structure finders, and the implementation of any such division is left to the end-user. Conversely, AdaptaHOP separately identifies DM structures (“DM haloes”) and stellar structures (“galaxies”) and so does intrinsically impose an explicit distinction between galactic and intergalactic stars: the local stellar density of each stellar particle is estimated using 20 nearest neighbour particles [122], and then a density threshold applied to split the two populations of stars. No unbinding procedure is employed. In Hz-AGN, this threshold density is \(178\times\) the cosmic average total matter density [73]. We adopt this local stellar density definition for (inter-)galactic stars, along with an arbitrary cluster border at \(r_{178c}\) (centred on the central galaxy), as our fiducial ICL definition for this study – calculating local stellar density following the methodology of [122] and [73].

2.2.2 Galaxy catalogue construction↩︎

Under our fiducial ICL definition all stellar particles with local stellar densities below \(178\times\) the cosmic average total matter density are considered intergalactic. For the purpose of then identifying “galaxies” among the remaining “galactic” stellar particles in TNG100, The300 GS-7K, and Hydrangea sub-haloes (referring to the mixed structures of DM and baryons identified by the structure finder), we assume each sub-halo to contain at most one galaxy, to which the sub-halo’s merger tree is associated (see Section 2.3.1). However, we do not simply assign all stellar particles with sufficient local density in each sub-halo to this galaxy. In massive haloes, the dynamically unbound but dense stellar envelopes around infalling satellites are reassigned to the central structure by the unbinding procedure (for further discussion see e.g. [129], [47]). We seek to exclude these spatially distinct stellar envelopes when identifying the central galaxies of massive haloes (associating these envelopes instead with the satellite galaxies they surround). For this purpose, we perform a particle clustering search on the galactic-density stellar particles of each sub-halo with the DBSCAN clustering algorithm ([130]). For this search we target neighbourhood densities of \(\geq178\times\) the cosmic average matter density and consider a neighbourhood size equal to the half mass radius of all the sub-halo’s “galactic” stellar particles, centred on the density maximum as found by a “shrinking sphere” approach. The largest structure identified by this clustering search is then used as a seed for the main stellar structure in each sub-halo and – once these searches have been carried out on all sub-haloes – leftover galactic density stars are allocated to these structure seeds by a density walk through their 20 nearest neighbour stellar particles, emulating how stars are assigned to local density peaks to construct galaxies in AdaptaHOP.

Subfind and Cantor sometimes fragment individual galaxies into several distinct objects. This behaviour is undesirable for our use case, and our procedure for reconsolidating these fragments into a single galaxy is described in Appendix 7. AHF follows the inclusive particle convention (in which the particle membership lists of structures can overlap) unlike the other structure finders, which assign particles to (at most) one structure each. Before carrying out our galaxy identification procedure, we obtain exclusive membership lists for AHF objects by assigning each particle to the lowest particle count structure it is a member of.

In Hz-AGN a minimum particle count of 51 for stellar structures to be kept is also imposed. For consistency, we impose this same minimum stellar particle count uniformly across the simulations: sub-threshold sub-structures of stars identified within a larger stellar structure are absorbed into this host galaxy, and the stellar particles of isolated small objects are instead recategorised as “intergalactic”. We confirm that imposing this minimum stellar particle count has no qualitative impact on any of our findings (though see Sections 5.2.1 and 5.2.3).

2.2.3 BCG determination and fiducial definitions summary↩︎

In Hz-AGN we identify the cluster BCG as the most massive galaxy within \(0.1\times r_{178c}\) of the cluster halo centre at \(z\approx0\), and assume the main progenitor of this galaxy to always be the BCG of the cluster (progenitor structure) in all earlier snapshots (see Section 2.3.1). In TNG100, The300 GS-7K, and Hydrangea we assume the cluster BCG to be the galaxy identified within the cluster’s central structure at \(z=0\) and that the main progenitor of this structure always contains the central galaxy of the cluster (progenitor) in all earlier snapshots.

Our fiducial BCG definition can thus be summarised as the BCG being the stellar structure associated with the cluster’s main halo, composed exclusively of stellar particles with local stellar densities greater than \(178\times\) the cosmic average total matter density (with no unbinding criteria employed). Under our fiducial ICL definition, all stellar particles within an outer cluster border at \(r_{178c}\) that fall below this same local stellar density threshold are considered part of the cluster’s ICL component (with a small additional contribution from unresolved satellites; Section 2.2.2).

2.2.4 BCG-ICL separation arbitrariness↩︎

Arguably the most disputed aspect of different ICL extraction methodologies is the separation of the BCG and the ICL: determining where to place the border between these two cluster components is non-trivial and somewhat subjective, and different selections can yield considerable differences in the determined properties of both components (e.g. [9], [24], [14]). As a result, rather than implement an arbitrary border, several prior ICL studies have opted not to separate the BCG from the ICL and consider the combined system only (e.g. [131], [45], [132]; see also [82] and [14] for further relevant discussion). We probe the consequences of our fiducial approach for separating the BCG and the ICL on determined ICL stellar mass fractions in Section 3.3, but also provide alternative versions of our analysis where applicable considering the combined BCG + ICL rather than the ICL alone.

2.3 Linking ICL stars to progenitor galaxies↩︎

We follow a very similar scheme for assigning \(z\approx0\) ICL stars to progenitor satellite galaxies as was implemented in . We first trace the main progenitor structures of \(z\approx0\) clusters back to a lookback time of approximately 12 Gyr (\(z\sim3\)). We then track the stellar particles of satellites falling into these structures after \(z\sim3\) and so generate lists of stellar particles associated with particular satellite galaxies, which can then be used to look up the progenitor objects of \(z\approx0\) ICL stars.

2.3.1 Snapshot spacing and merger trees↩︎

We perform our analysis using simulation snapshots with a coarse time resolution of \(\sim1\) Gyr (henceforth referred to as “coarsely spaced snapshots”) – selected to limit the computational expense of our analysis while still employing a spacing smaller than one cluster crossing time, allowing each orbit of an infalling satellite to be sampled multiple times. We follow structures between snapshots in Hz-AGN using merger trees generated by the Treemaker algorithm ([133], [122]) built using a time resolution of \(\sim0.05\) Gyr between \(z\approx3\) and \(z\approx0\). In TNG100 we use the provided merger trees output by the Sublink algorithm (we specifically use the baryonic-based variant trees, “Sublink_gal”; [129]). In The300 GS-7K we use the provided outputs of the MergerTree algorithm incorporated into AHF (we specifically use the “skipping+stars” variant trees first described in [134] that incorporate tree break patching and consider both DM and stars in the utilised merit function). In Hydrangea we use the provided merger trees for the Cantor catalogue objects as generated by the Spiderweb algorithm [108].

It is known that during a near-binary merger between a central galaxy and a massive satellite, structure finders can sometimes behave erratically – with significant “mass sloshing” between the central and satellite sub-haloes, and with which structure is considered the central able to flip back-and-forth between consecutive snapshots (the so-called “sub-halo switching problem”; for further discussion see e.g. [135], [136], [137]). To limit the potential impact of this behaviour on our analysis, for some clusters we slightly stray (for a single snapshot) from our usual snapshot spacing of \(\sim1\) Gyr to avoid using snapshots when the cluster central is in the midst of such a merger.

2.3.2 Tracking infalling satellites↩︎

We consider galaxies to become satellites of clusters once they pass \(r_{178c}\), and track the stellar particles of all galaxies that become cluster satellites since a lookback time of \(\sim12\) Gyr in order to associate \(z\approx0\) ICL stars with progenitor objects. Infalling galaxies are identified as galaxies with centres-of-mass outside \(r_{178c}\) whose main descendants in the next coarsely spaced snapshot are cluster satellites. We also label as an infalling satellite any galaxy outside \(r_{178c}\) which has no detected descendant in the next coarsely spaced snapshot (or that merges into the cluster’s central structure) provided we find the centre-of-mass of all stellar particles formerly belonging to that galaxy within \(r_{178c}\). To avoid double counting from “splashback” galaxies that resurface above \(r_{178c}\) before falling back into the cluster, this infalling satellite identification is performed chronologically and each galaxy’s first observed \(r_{178c}\) crossing considered its original infall.

All stellar particles belonging to an infalling galaxy when it is first identified as such are tagged as having that galaxy as their progenitor object. We then follow the main descendant of that progenitor galaxy through all remaining coarsely spaced snapshots, tagging to the same progenitor stellar particles formed within that galaxy after first cluster infall. In order to avoid double counting due to galaxy mergers, we stop tagging new stellar particles to a progenitor once the original infalling galaxy is no longer the main progenitor of the descendant galaxy those stellar particles form in.

It is known that merger tree algorithms can sometimes struggle in the extreme environments of clusters, occasionally leading to tracking failures or inappropriate links between objects (for further discussion, see e.g. [138], [139], [137]). When a satellite either shares zero stellar particles with its supposed descendant or has no apparent descendant in the next coarsely spaced snapshot, and yet another galaxy exists containing \(>50\) per cent of the stellar particles formerly belonging to the satellite, we overrule the merger tree and connect the tracked satellite to this assumed “correct” descendant.

We do not attempt to separately quantify the ICL contributions from tidal stripping and violent relaxation during galaxy mergers as part of this study, and refer to the combination of these processes as the liberation of stars from satellite galaxies. Disentangling the contributions from these two processes during a major galaxy merger is non-trivial, and the proportional ICL contributions of the two channels highly sensitive to the arbitrary demarcation imposed (see [23] and references therein for further relevant discussion).

2.4 Cluster sample selection↩︎

As no structures more massive than \(10^{15}\,\)M are present in either Hz-AGN or TNG100, we restrict ourselves to The300 GS-7K and Hydrangea clusters with \(z=0\) halo masses less than \(10^{15}\,\)M to limit the mass discrepancy between our cluster samples. We also exclude from our samples any clusters with a central merger tree that cannot be followed back to a lookback time of \(\sim12\) Gyr. As described in Pearce et al. (in prep.), in some The300 GS-7K clusters the main merger tree branch of the cluster’s DM halo clearly decouples from that of the associated \(z=0\) central galaxy (due to a “stellar core switch” occurring during a close encounter between two massive systems). Our assumption that the main progenitor of the \(z=0\) central structure always contains the main progenitor of the identified \(z=0\) BCG in TNG100, The300 GS-7K, and Hydrangea breaks down if such a “core switch” occurs. We thus follow the methodology of Pearce et al. (in prep.; using our local stellar density defined central galaxies and coarsely spaced snapshots) to identify clusters that experience such a decoupling at \(z\lesssim2\) and exclude these from our cluster samples as well. Additional details on all the simulated clusters we consider in this study are provided in Appendix 8, including individual cluster identifiers.

In the latest-time Hz-AGN simulation snapshot we use for this work (snapshot 761, \(z=0.0556\)), there are \(14\) cluster-scale haloes (\(M_{178c}\gtrsim10^{14}\,\)M). We consider only the sub-set of 10 for which we are able to follow the BCG back to a lookback time of \(\sim12\) Gyrs. These clusters have \(M_{178c}\) values in the range \([1.00,7.52]\times10^{14}\) M.

Our TNG100 sample consists of the 15 most massive Friends-of-Friends (FoF) groups in the TNG100 simulation box at \(z=0\) (snapshot 99; we use \(M_{200c}\) masses for this selection). These clusters have \(M_{178c}\) values in the range \([0.94,3.72]\times10^{14}\) M. Four of these clusters (snapshot 99 FoF IDs 0, 4, 5, and 10; see Appendix 8) are within 5 Mpc of the simulation boundaries at \(z=0\), with two of these (IDs 0 and 10) crossing the boundaries during our period of study. Addressing concerns about edge-effects in TNG100 raised by [140] and [141], we confirm that the inclusion of these clusters does not meaningfully alter any of our results.

We restrict ourselves to The300 GS-7K clusters less massive than \(10^{15}\) M at \(z=0\) (using \(M_{200c}\) masses for this selection), and also exclude clusters with central structure merger trees that either cannot be followed back to a lookback time of \(\sim12\) Gyr or that are otherwise problematic (see above). This yielded a The300 GS-7K sample consisting of 15 clusters with \(M_{178c}\) values in the range \([8.56,10.71]\times10^{14}\) M.

Of the 24 Hydrangea clusters, we restrict ourselves to the sub-set of 18 clusters with \(M_{178c}(z=0)\lesssim10^{15}\,\)M. From this sub-set, we exclude two further clusters: one due to issues with the merger tree of its central structure (see above) and another for an unphysical AGN outburst at \(z\sim1.5\) (identifiers CE-11 and CE-10 respectively; [107]). This yielded a sample of 16 clusters, which have \(M_{178c}\) values in the range \([1.16,6.86]\times10^{14}\) M.

3 Stellar mass fractions↩︎

3.1 Distribution of cluster stellar mass at \(z\approx0\)↩︎

We begin by examining the \(z\approx0\) ICL (and BCG + ICL) stellar mass fractions of the cluster samples drawn from each of the four simulations, homogenized to our fiducial ICL definition based on local stellar density (Section 2.2), to probe for significant intrinsic inter-simulation differences. We present in Figure 1 violin plots showing the distribution of cluster stellar mass (within \(r_{178c}\)) between satellite galaxies, the BCG, and the ICL at \(z\approx0\) separately for the cluster sample from each simulation. Combined BCG + ICL fractions are also included. Stellar mass fractions are given for each cluster individually in Appendix 8.

The ICL fractions of the four cluster samples are broadly consistent, all lying in the range \(\sim0.1\) to \(\sim0.2\), which is far less scatter than the wide range of values reported for various different ICL extraction methodologies by prior theoretical studies (between \(\sim0.05\) and \(>0.5\); [18]). We note a slight dispersion in typical ICL fraction between the four simulations, which we discuss in Section 5.2.1. The ICL fractions we obtain lie towards the lower end of the broad range of previously reported values, but are markedly similar to the findings of past studies employing comparable ICL definitions, such as [25], who reported \(z\approx0\) ICL fractions between \(\sim0.09\) and \(\sim0.15\) for a sample of similarly massive simulated clusters, and [22], who reported a typical ICL fraction of \(\sim0.12\) for a large sample of simulated haloes (mass \(>10^{13}\,\)M) at \(z\sim0.625\).

While the BCG + ICL stellar mass fractions of the four simulations span very similar ranges (\(\sim0.3\) to \(\sim0.7\)), there is noticeable scatter between the average BCG + ICL fractions, driven chiefly by a corresponding scatter in typical BCG mass fractions. Most apparent is the comparatively low average BCG fraction for The300 GS-7K (of only \(\sim0.2\)), though the distribution of values in TNG100 also appears mildly top-heavy, consistent with the findings of [9] and [24], who both reported slightly higher average BCG + ICL fractions in TNG100 clusters compared to those from Hz-AGN. We discuss possible factors contributing to these differing typical BCG (and hence also BCG + ICL) mass fractions in Section 5.2.2 (see also Appendix 8) but in brief we largely attribute the usually lower BCG fractions of The300 GS-7K to these generally more massive clusters typically assembling more of their stellar mass closer to \(z\approx0\).

3.2 ICL mass fraction evolution↩︎

We examine next the evolution of ICL (and BCG + ICL) stellar mass fractions with redshift. Prior theoretical works have reported conflicting findings for how the ICL component evolves with redshift (e.g. [25], [24]). We probe for significant inter-simulation discrepancies which might contribute to this tension, unobstructed by inconsistent ICL extraction methodologies. Figure 2 depicts median ICL stellar mass fractions as a function of lookback time for the samples of \(z\approx0\) galaxy clusters and their \(z>0\) progenitor structures⁷ from the four simulations, homogenized to our fiducial ICL definition (top panels; the bottom panels indicate typical progenitor structure mass as a function of lookback time). The shaded regions are bounded by 16th and 84th percentiles as an indication of the cluster-to-cluster scatter. Combined BCG + ICL fractions are also included. We present an alternative version of this analysis confined to cluster mass (\(\gtrsim10^{14}\,\)M) structures at all times in Appendix 9.

Although the ICL and BCG + ICL stellar mass fractions of individual structures can fluctuate over time, we find no significant evolution in either the average ICL or BCG + ICL mass fraction in any simulation until at least \(z\sim2\). At \(z\gtrsim2\) there are hints of a common downturn in typical ICL fraction and the average BCG + ICL fractions all begin to rise, which we consider in both cases symptomatic of progenitor structures falling below the group scale as opposed to intrinsic redshift evolution.

Our finding of no evolution in the average ICL fraction (at least for \(z\lesssim2\)) is in agreement with the results of some prior theoretical studies – [23] predicted no trend in ICL fraction with redshift (for \(\gtrsim10^{13}\,\)M haloes at \(0<z<3\) based on semi-analytical models) and [22] found a nearly constant median ICL fraction across various redshifts (for \(\gtrsim10^{13}\,\)M haloes at \(z=0.625\) and their progenitors back to \(z\sim3\) using the Horizon Run 5 hydrodynamical simulation) – but disagrees with the findings of [25], who reported ICL fractions that steadily fall with increasing lookback time, approaching zero by \(z\sim2\). Although [25] employed a density-based ICL definition broadly similar to our fiducial ICL definition, their methodology differs by imposing a fixed stellar density threshold to distinguish ICL stars (e.g. \(\rho_{\textrm{thresh}} = 10^{-5}\,\)M pc^-3), whereas for this study we emulate AdaptaHOP and employ a stellar density threshold that scales with the redshift-dependent cosmic average matter density (e.g. \(\rho_{\textrm{thresh}}(z=0)\sim7\times10^{-6}\,\)M pc^-3, \(\rho_{\textrm{thresh}}(z=2)\sim2\times10^{-4}\,\)M pc^-3) for our fiducial ICL definition.

3.3 Influence of BCG-ICL separation↩︎

In order to probe the influence that the fiducial ICL extraction methodology we employ for this study has on our findings for the evolution of the ICL fraction with lookback time, we repeat our analysis from Section 3.2 for a selection of common alternative approaches for separating the BCG from the ICL. The first alternative BCG-ICL separation we consider is a fixed local stellar density threshold, emulating [25] and using the same fixed stellar density threshold (of \(\rho_{*}=10^{-5}\) M pc\(^{-3}\)). We additionally consider three alternative methodologies all based on different spherical aperture cuts: splitting BCG from ICL at 10 per cent of \(r_{178c}\); splitting at three times the combined BCG + ICL stellar half-mass radius, \(r_{*,1/2}\); and splitting at an arbitrary fixed radius. Lastly, we also consider a kinetic BCG-ICL separation employing a double Maxwellian fit to the velocities of BCG + ICL stellar particles (following the methodology of [142], [143] as implemented in [9]). For these repeat analyses, we retain both our arbitrary cluster border at \(r_{178c}\) and our fiducial procedure for demarcating satellite galaxies from the ICL.

The results of these alternative analyses are depicted in Figure 3, which shows how the median ICL fractions for each of the different BCG-ICL separation methodologies evolve as a function of lookback time for each of the four simulated cluster samples. We present an alternative version of this analysis restricted to cluster mass (\(\gtrsim10^{14}\,\)M) structures at all times in Appendix 9.

Fixed apertures with radii between 30 and 200 kpc were considered, but we include only the results for a 100 kpc aperture in Figure 3 as a representative example; smaller apertures result in larger characteristic ICL fractions (e.g. \(\sim0.3-0.4\) at \(z\approx0\) for a 30 kpc aperture compared to \(\sim0.1-0.2\) for a 200 kpc aperture) but follow the same overall trends with lookback time (see [14] for further relevant discussion). Similarly, separating the BCG from the ICL at \(2\times r_{*,1/2}\) (as in e.g. [82]) rather than \(3\times r_{*,1/2}\) does not alter the overall trends with lookback time shown in the bottom middle panel of Figure 3 and just raises the typical ICL mass fractions found slightly (such that the median ICL fractions at \(z\approx0\) fall in the range \(\sim0.12-0.20\)). Projected circular apertures with fixed radii were also examined, which yield ICL fractions typically just below their spherical cut equivalent (e.g. within 4 percentage points at \(z\approx0\) for 100 kpc) and reproduce qualitatively the same trends with lookback time.

Considerable scatter is apparent in Figure 3 between the typical ICL fractions that the different BCG-ICL split methodologies yield, even for the same sample of clusters from the same simulation (a phenomenon also reported in prior studies e.g. [9], [14]). Particularly conspicuous is the kinematic methodology yielding typical ICL fractions significantly larger than any other BCG-ICL separation considered (with typical \(z\approx0\) ICL fractions in the range \(\sim0.2-0.3\); approximately twice those obtained by the other methods), qualitatively similar to the findings of [9].

How the typical ICL fraction evolves as a function of lookback time is also not consistent between the different BCG-ICL separation methodologies. Under the fixed aperture and fixed density cut methodologies, all four simulations yield ICL fractions that decrease with lookback time, approaching zero by \(z\sim3\), reproducing the findings of [25]. In contrast, the kinematic separation methodology as well as BCG-ICL cuts scaled by \(r_{178c}\) or \(r_{*,1/2}\) yield no significant evolution in the typical ICL fraction from \(z\approx0\) to \(z\sim2\), just as with our fiducial ICL extraction methodology. However, it merits highlighting that when homogenized to the same BCG-ICL split methodology, there are no significant discrepancies between the different simulations for the evolution of the ICL fraction as a function of lookback time. Essentially the same trends are reproduced by every considered simulation for each considered BCG-ICL separation methodology, and so the different trends with lookback time appear to just stem from the different ICL definitions being imposed.

[24] previously investigated how the fraction of stellar mass in different cluster components evolved with time for galaxy clusters from the Magneticum simulation for various BCG-ICL splits. They used effectively the same selection of BCG-ICL separations we trial in the three bottom panels of Figure 3 (i.e. spherical aperture cuts based on halo over-density radii, stellar half-mass radii, and fixed apertures; see their figure B.2) and recovered the same trends we observe: little-to-no evolution of typical ICL mass fractions with lookback time (back to \(z\sim2\)) except for fixed aperture approaches, which yield fractions that fall with increasing redshift.

For the analysis presented in Figure 3, we do not control for evolving halo or central-galaxy mass when applying the fixed aperture methodology to progenitor structures at \(z>0\). This is the primary cause behind the consistent redshift evolution trends seen for this methodology, which are largely curtailed in the alternative version of this analysis limited to cluster mass objects which we present in Appendix 9. However, it also bears highlighting that for a fixed mass both galaxies and haloes are typically more centrally concentrated and physically smaller at higher redshifts (e.g. [144], [145], [146]). Consequently, any aperture of fixed physical size would be expected to enclose more of a typical cluster mass halo (and combined BCG + ICL system) at \(z\gtrsim2\) as compared to \(z\approx0\). Not controlling for evolving typical halo or galaxy size could thus also be a secondary driver of the evolution towards lower ICL fractions at higher redshift seen with the fixed aperture methodology (traces of which are still present in the alternative analysis employing a mass-cut shown in Appendix 9).

4 Assembly of ICL from cluster satellites↩︎

Previous theoretical works have occasionally disagreed on whether the main progenitor galaxies of ICL stars should be more or less massive galaxies, with differing predictions sometimes attributed to the different methodologies implemented to separate galactic and intergalactic stars (e.g. [51]). In this section, we investigate which satellites are the dominant contributors of stars to the \(z\approx0\) ICL as predicted by each of the four simulations (using an approach similar to that employed in ) – probing for significant discrepancies between the simulations for the predicted peak progenitor masses despite implementing a uniform scheme for separating galactic and intergalactic stars. For this analysis, we return to our fiducial ICL extraction methodology based on local stellar density (as compared to the cosmic average matter density; Section 2.2).

4.1 The ICL contribution from a cluster’s own satellites↩︎

Our scheme for linking ICL stars to progenitor galaxies (described in Section 2.3) only considers stars liberated on or after cluster infall from galaxies we track falling into our sample clusters after \(z\sim3\), and is not recursive for the pre-existing diffuse stellar component of smaller clusters or massive groups accreted during cluster assembly. Though tidal stripping is generally acknowledged to be the dominant assembly channel for ICL stars by \(z=0\), this usually requires considering this “pre-processing” of ICL stars in accreted groups to be an indirect channel (as in e.g. [23]) – with prior studies often finding the contribution of this pre-processing channel to be substantial when it is considered separately.

In Figure 4 we present for each of the four simulated cluster samples the mean fraction of \(z=0\) ICL (and BCG + ICL) stellar mass assembled by the liberation of stars within a cluster’s own halo from satellite galaxies that join the cluster after \(z\sim3\). These mean ICL mass fractions span the range \(\sim0.3-0.6\) (\(\sim0.4-0.7\) for BCG + ICL), always falling significantly short of unity, in qualitative agreement with the findings of prior studies: [23] reported pre-processing to typically contribute between 20 and 40 per cent of the \(z=0\) ICL stellar mass in \(10^{14}-10^{15}\) M halo mass clusters (based on semi-analytic models); and [132] reported only \(\sim55\) per cent of the combined BCG + ICL stellar mass in a \(z=0.79\) cluster (from the NewCluster simulation; \(z=0.79\) halo mass \(\sim1\times10^{14}\) M) to originate from satellite galaxies, with ICL pre-processing already contributing \(>10\) per cent of the combined BCG + ICL mass even at this early redshift.

Though we defer a more comprehensive analysis of ICL assembly via channels other than the liberation of satellite galaxy stars to the forthcoming companion paper (Brown et al. in prep.), we confirm that in all four considered simulations most of the \(z\approx0\) ICL mass not liberated within the cluster halo from a cluster’s own satellite galaxies can be attributed to ICL pre-processing – with stars formed within the central galaxy since \(z\sim3\) also typically making a substantial contribution to the combined BCG + ICL (see Section 5.2.2). We defer discussion of the inter-simulation differences in Figure 4 to Section 5.2.3. Though for the remainder of our analysis we consider only the sub-set of ICL stars liberated within the cluster halo from a cluster’s own satellite galaxies, we highlight that under hierarchical assembly it is broadly expected that this pre-existing diffuse stellar component of massive groups should be assembled in a self-similar fashion to the ICL of clusters (see Section 5.3 for further discussion).

4.2 Progenitor galaxies of ICL (and BCG + ICL) stars↩︎

4.2.1 Liberated fraction of satellite stellar mass↩︎

To begin our investigation into which satellites are the largest contributors of ICL stars overall, we first probe how the typical ICL contribution of an individual satellite galaxy varies as a function of cluster infall stellar mass. We refer to the fraction of the stellar particles associated with a satellite – meaning those stellar particles composing the galaxy just proceeding first cluster infall or that were born in the descendant of that galaxy thereafter (see Section 2.3.2) – that become part of the \(z\approx0\) ICL component as the liberated fraction, \(f_\text{lib-ICL}\), of that satellite.

For each of the four simulations individually we calculate the mean value of \(f_\text{lib-ICL}\) for 0.5 dex wide rolling bins (bin step 0.15 dex) in infall stellar mass, \(M_{*}\), for satellite galaxies that enter any of the clusters for the first time at least \(1\,\)Gyr before \(z\approx0\). We then fit a cubic spline function to these mean \(f_\text{lib-ICL}\) values as a function of infall stellar mass, ignoring both bins containing poorly-resolved galaxies (infall masses corresponding to fewer than 100 stellar particles) and also poorly sampled bins (fewer than 5 objects per cluster). We linearly extrapolate this fit beyond the range of masses used for fitting (with the output constrained to \([0,1]\)).

We plot the resulting fitted curves for each of the four simulations in the top panel of Figure 5. The bottom panel of Figure 5 shows the equivalent analysis for the combined BCG + ICL. The shaded regions indicate the dispersion (16th-84th percentile) of the fit lines generated in the same fashion for \(10^{4}\) bootstrap resamples of the aggregate infalling satellite population of each simulation. The data each fit is based on are presented in Appendix 10. Though we include stars formed post-infall when calculating \(f_{\text{lib-ICL}}\), we confirm that this has no qualitative impact on the presented findings. The fitted functions shown in Figure 5 are not meaningfully different if reproduced as a function of total stellar mass associated to an infaller by \(z=0\) rather than infall mass. Liberated satellite stars found at \(r>r_{178c}\) at \(z\approx0\) are not considered \(z\approx0\) ICL stars when we calculate \(f_{\text{lib-ICL}}\); we also exclude satellites just crossing \(r_{178c}\) for the first time at \(z\approx0\) as well as those galaxies which only “skim” clusters, entering then exiting \(r_{178c}\) in consecutive coarsely spaced snapshots and never returning before \(z\approx0\). We confirm that none of these exclusions significantly influence our findings.

Qualitatively the same trend between expected \(f_\text{lib-ICL}\) value and infall mass is reproduced by all four simulations: comparatively higher average \(f_\text{lib-ICL}\) fractions at lower infall masses, which decrease with increasing infall mass until \(M_{*}\sim10^{10}\,\)M, at which point they level out around \(f_\text{lib-ICL}\sim0.05-0.1\). An initially similar general trend is also seen for the equivalent analyses with the combined BCG + ICL, differing chiefly through \(f_\text{lib-BCG+ICL}\) beginning to rise again at high masses (\(M_{*}\gtrsim10^{10.5}\,\)M; from a rapidly rising typical BCG contribution) rather than level out. Essentially the same qualitative trends were previously reported in , and we discuss the physical origins of these trends in Section 5.1. We defer discussion of the inter-simulation differences to Section 5.2.4.

4.2.2 Infalling satellite mass function↩︎

In addition to understanding how the typical ICL mass contribution per galaxy varies with satellite infall mass, determining which population of satellites contribute the bulk of ICL stars also requires quantifying the relative population sizes of differingly massive satellites falling into assembling clusters. We present in Figure 6 infalling satellite mass function fits for each simulation for galaxies that enter any of the clusters for the first time at least \(1\,\)Gyr before \(z\approx0\). These fits adopt the form of the [147] function, \[\Phi(M_{*})\cdot dM_{*} = \Phi_{n} \left( \frac{M_{*}}{M_{k}} \right)^{\alpha} \exp\left({-M_{*} / M_{k}} \right) \cdot dM_{*} \label{eq:Schechter}\tag{1}\] where \(\alpha\) is the low-mass-end slope, \(\Phi_{n}\) a normalization parameter, and \(M_{k}\) a characteristic mass (corresponding to the “knee” of the function i.e. the mass when the function exhibits a rapid change in slope). The data each fit is based on are presented in Appendix 10, and fitted parameter values are given in Table 2. To compensate for resolution effects and enable predictions for the ICL contribution made by low-mass satellites, only stellar masses corresponding to well-resolved galaxies (\(\geq100\) stellar particles) are considered for fitting and the fitted functions then extrapolated to lower masses. To facilitate comparison between the low-mass slopes and “knee” positions of each fit, the mass functions depicted in Figure 6 have all been arbitrarily renormalised to each correspond to a total accreted stellar mass of \(5\times10^{12}\,\)M.

We assume both parameter fit uncertainty and the dispersion in fitted parameters values across \(10^{4}\) bootstrap resamples of the satellite population to be independent sources of Gaussian error, and combine these in quadrature to obtain the estimated errors on fit parameter value given in Table 2. The shaded regions in Figure 6 then show 16th-84th percentile confidence intervals on the fitted functions for this estimated uncertainty. Within estimated uncertainty the infalling satellite mass function fits for the different simulations appear in reasonable agreement (see Section 5.2.5 for a brief discussion of the inter-simulation differences).

4.2.3 Star formation after cluster infall↩︎

A complete accounting of the ICL contribution by satellite galaxies requires acknowledging the population of stars formed in some satellites after cluster infall (rather than solely the initial infall mass of these galaxies). In , this post-infall star formation was accounted for with a constant factor equal to the mean ratio between total associated stellar mass by \(z\approx0\) and infall mass across the entire population of tracked satellite galaxies. Here we instead calculate the mean ratio between total associated stellar mass by \(z\approx0\) and infall stellar mass for the same rolling bins used in the analysis depicted in Figure 5. We fit a cubic spline function to these mean mass ratios in the same manner as our fitted \(f_{\mathrm{lib}}(M_{*})\) functions, including linear extrapolation beyond the range of masses used for fitting (though with the output constrained to \([1,\infty]\) rather than \([0,1]\) and no longer disregarding low-mass satellites when fitting). The resulting fitted functions can be found in Appendix 11 but all adopt the same essential shape, with a peak at intermediate infall masses (\(10^{9.5}-10^{10}\,\)M), falling back to unity at lower and higher masses. We confirm none of our remaining findings are significantly altered by substituting this approach for that used in .

4.2.4 Total ICL contributions from galaxies of different mass↩︎

The product of our fitted functions for liberated fraction of stellar mass (\(f_{\text{lib-ICL}}(M_{*})\); Figure 5) and total associated stellar mass by \(z\approx0\) (\(M^{z=0}_{*,\textrm{Tot}}(M_{*})\); Appendix 11) yield a prediction for the typical ICL mass contribution per satellite for a given infall stellar mass. Combining this with our infalling satellite mass function fit (\(\Phi(M_{*})\); Figure 6) then allows us to make a prediction for the typical mass fraction of satellite galaxy sourced ICL stars at \(z\approx0\) contributed by satellites with infall mass \(M_{*}\), i.e. \[f_\text{ICL}(M_{*})=\frac{f_{\mathrm{lib}}(M_{*})M^{z=0}_\text{*,\textrm{Tot}}(M_{*})\Phi(M_{*}) }{ \int f_{\mathrm{lib}}(M_{*}) M^{z=0}_\text{*,\textrm{Tot}}(M_{*})\Phi(M_{*}) \cdot dM_{*} }. \label{eq:final95fit}\tag{2}\] The results for each of the four simulations are shown in Figure 7 (top panel; equivalent analysis for the combined BCG + ICL shown in the bottom panel), normalized to per dex in infall mass, with dotted lines indicating any reliance on extrapolation beyond the mass ranges used to fit the constituent functions. The shaded regions indicate estimated uncertainty (16th-84th percentile confidence interval) from propagating the uncertainties of the component fitted functions. For convenience, we state the specific infall stellar masses corresponding to the peak contributors of ICL (and BCG + ICL) stars as per Figure 7 in Table 3 (with uncertainties estimated from the 16th-84th percentile scatter of the peak masses from the \(10^{4}\) bootstrap resample repeat analyses). We re-emphasise that the progenitor analysis presented in Figure 7 considers only stars liberated within the cluster halo from a cluster’s own satellite galaxies – and does not include the “pre-processed” diffuse stellar component of smaller clusters and groups accreted during cluster assembly (Section 4.1; see also Section 5.3 for relevant discussion).

Though there are some noticeable differences between the four predicted \(f_\text{ICL}(M_{*})\) curves (which we discuss in Section 5.2.6), it can be seen in Figure 7 (and Table 3) that all four simulations predict the dominant progenitors of ICL stars to be massive galaxies, with the peak contributors being roughly Milky Way mass galaxies (\(M_{*}\sim10^{11}\,\)M), in good agreement with the findings of several prior theoretical studies (e.g. [43], [44], [45], [46], [47], [148], [15]). The peak contributors of BCG + ICL stars are then predicted by every simulation to be very slightly more massive (by \(\sim0.05-0.2\) dex) than for the ICL alone – consistent with the findings of previous works (e.g. [49], [15]).

We present the same analysis reframed as the predicted cumulative fractional contribution of satellite galaxy sourced ICL stars against satellite infall mass in Figure 8. A comparison between these fitted cumulative fractional ICL contribution curves and the same obtained directly from the simulations for the aggregate ICL of each cluster sample are presented in Appendix 12.

From Figure 8 it can be seen that all four simulations predict \(\gtrsim50\) per cent of the satellite sourced ICL mass of an average \(z\approx0\) cluster to be contributed by \(M_{*}\gtrsim10^{10}\,\)M satellites (\(M_{*}\gtrsim10^{10.5}\,\)M for every simulation besides Hydrangea), with \(\sim80\) per cent from \(M_{*}\gtrsim10^{9}\,\)M satellites (\(M_{*}\gtrsim10^{9.5}\,\)M for every simulation besides Hydrangea). The corresponding minimum satellite infall mass thresholds for the combined BCG + ICL are \(M_{*}\gtrsim10^{10.5}\,\)M and \(M_{*}\gtrsim10^{10}\,\)M for \(\sim50\) per cent and \(\sim80\) per cent (of the satellite sourced BCG + ICL stellar mass), respectively.

In a prior study employing TNG100, [15] reported the mean ICL mass fraction attributable to progenitor objects with stellar masses \(\gtrsim10^{10}\,\)M to be \(\sim60\) per cent. This is broadly in line with our findings, with our slightly higher equivalent prediction of \(\sim70\) per cent using TNG100 thought chiefly due to [15] also including in their sample galaxy groups (which should be assembled from on average slightly less massive progenitors than clusters). The minimum progenitor stellar mass required to recover at least \(50\) per cent of the BCG + ICL in TNG100 clusters was previously reported by [82] to be \(\sim10^{10.5}\,\)M (requiring \(\gtrsim10^{9.5}\,\)M progenitors for \(\geq90\) per cent; see their figure 13), also in reasonably good agreement with our findings.

Using the lower-resolution predecessor simulations of The300 GS-7K, [148] reported \(M_{*}\geq10^{11}\,\)M progenitors to contribute \(65-80\) per cent of the ICL in the \(14.8 < \text{log}_{10}(M_{200}/\text{M}_{\sun}) < 15.6\) clusters they considered, in qualitative agreement with our result that massive galaxies appear the dominant contributors of ICL stars, though we do not find quite so large a contribution from these very massive galaxies (instead predicting a typical contribution \(\sim30-40\) per cent). This discrepancy is likely due in part to [148] considering typically more massive clusters, and that they only consider the ICL contribution from objects explicitly resolved in their (lower-resolution) simulations, whereas we extrapolate to predict the contribution from unresolved low mass galaxies as well. Additionally, the definitions they use include the pre-existing diffuse stellar component of groups that merge into their clusters as part of the ICL contribution by the former central galaxies of those groups (see Section 5.3 for further relevant discussion).

5 sec:Discussion↩︎

5.1 Physical origins of liberated fraction trends↩︎

All four simulations reproduce the same general trends in Figure 5: decreasing \(f_\text{lib-ICL}\) values with increasing satellite infall mass, and \(f_\text{lib-BCG+ICL}\) values that initially fall with increasing satellite mass before reversing and beginning to rapidly rise at \(M_{*}\gtrsim10^{10.5}\) M. These are qualitatively the same trends noted and discussed in , and we consider these to emerge chiefly from the convolution of two phenomena: more massive satellites generally experience more rapid orbital decay in clusters due to dynamical friction (e.g. [149], [150], [151]); and more massive galaxies also have deeper potential wells and proportionally larger DM haloes (once past the peak in the stellar-to-halo mass relation), so are less readily stripped of stars by gravitational interactions (e.g. [152], [153], [39]).

The orbits of less massive satellites typically taking longer to decay provides ample opportunity for a significant fraction of their more readily stripped stellar mass to be gradually siphoned into the ICL as they spiral slowly inwards towards the cluster centre, yielding higher average \(f_\text{lib-ICL}\) values but in general not meaningfully contributing to the BCG. By contrast, the very massive former-central galaxies accreted when a cluster subsumes a smaller cluster or massive group will not only be highly resistant to stripping, but also typically experience only a brief window for this stripping to occur between cluster infall and the inevitable merger with the cluster central galaxy relatively soon thereafter, yielding low average \(f_\text{lib-ICL}\) values but high average \(f_\text{lib-BCG+ICL}\) values at the highest satellite masses (see Pearce et al. in prep. for further relevant discussion). We then consider the levelling out of \(f_\text{lib-ICL}\) at high-masses indicative of violent relaxation processes during BCG mergers becoming an increasingly significant mechanism for liberating stars to the ICL for especially massive galaxies ().

Despite the scheme described above, it can be seen (by comparison between the two panels of Figure 5) that the average BCG stellar mass contribution is generally non-zero even for the lowest infall masses. This is driven chiefly by low mass satellites which join to-be clusters very early (first infall at \(z\gtrsim2\)) which, having been afforded \(\gtrsim10\) Gy of orbital decay after joining the much smaller, early progenitor of a \(z\approx0\) cluster (and having joined with a more significant infall mass ratio compared to a similarly massive galaxy at lower redshifts; Pearce et al. in prep.) are able to eventually merge into the BCG, resulting in the underlying distribution being mildly bimodal. We note that earlier infall times generally correspond to higher average \(f_\text{lib-ICL}\) and \(f_\text{lib-BCG+ICL}\) values in all four simulations for all infall masses, but qualitatively the same trends seen in Figure 5 persist even when considering only satellites with a specific infall time.

5.2 Inter-simulation differences↩︎

5.2.1 ICL mass fractions↩︎

A slight dispersion between the typical ICL stellar mass fractions of the four simulations at \(z\approx0\) can be seen in Figure 1, with the typical values for Hz-AGN and TNG100 (median 0.13 and 0.12 respectively) falling slightly below those for The300 GS-7K and Hydrangea (median 0.18 and 0.17 respectively). It is initially curious that these typical fractions appear stratified by simulation approach (uniform-volume vs zoom-in) but, while it is conceivable that such a stratification could stem from some bias associated with simulation approach (e.g. the Hz-AGN and TNG100 clusters occupying the very top of the mass hierarchy in those simulations and so being biased towards rapid recent mass growth; [154]), that this stratification is usually not retained under the various alternative BCG-ICL split methodologies trialled in Figure 3 prompts us to consider it somewhat coincidental.

The fiducial ICL definition we implement is based on instantaneous local stellar density, thus dynamically unbound stars transiently located within a region of sufficiently high density are in that instance regarded as galactic and not as ICL stars. This is pertinent for the Hz-AGN ICL fractions we determine, as Hz-AGN is known to have central galaxies with larger than expected effective radii (as compared to observations; e.g. [155], [156], see also figure 1 in ). In combination with the ICL being highly centrally concentrated, these distended central galaxies may be decreasing the ICL mass fractions we find in Hz-AGN somewhat, with a population of stars that might have otherwise contributed to the ICL occluded beneath the extended BCG.

We consider the comparatively high ICL fractions of the The300 GS-7K sample shown in Figure 1 to at least partially be a consequence of the lower resolution of this simulation (and also of our adopted galaxy identification methodology). Reduced simulation resolution artificially enhances the efficiency of tidal stripping ([39], [69]; see also [43], [157]). One might therefore anticipate that were the resolution of The300 GS-7K improved to be on par with the other three simulations, the typical ICL stellar mass fractions yielded may fall slightly. It is however also true that stellar stripping rates are not yet fully converged even at the higher resolutions of Hz-AGN, TNG100, and Hydrangea, thus it could equally be said that the ICL fractions predicted by these simulations may be slight overestimates as well.

Additionally, we impose a minimum stellar particle count of 51 when identifying galaxies in TNG100, The300 GS-7K, and Hydrangea (for consistency with the implementation of AdaptaHOP employed in Hz-AGN; Section 2.2), corresponding to a minimum galaxy stellar mass of approximately \(1\times10^{8}\) M in Hz-AGN, \(5\times10^{7}\) M in TNG100 and Hydrangea, and \(2\times10^{9}\) M in The300 GS-7K. If this constraint is relaxed, and a minimum particle count of unity instead employed, the median ICL fraction of The300 GS-7K under our fiducial ICL methodology falls to \(\sim0.16\) – indicating that the ICL fractions in The300 GS-7K from our main analysis include a small contribution from the stellar mass of neglected, poorly-resolved galaxies (of stellar mass \(\lesssim2\times10^9\,\)M; the median fractions of TNG100 and Hydrangea are unaltered by this change). The typical ICL fractions of the The300 GS-7K sample under the various alternative methodologies considered in Figure 3 likewise fall slightly when this constraint is relaxed (usually by \(\lesssim0.02\); with TNG100 and Hydrangea again broadly unaffected), though none of the redshift evolution trends shown in Figure 3 are qualitatively altered by relaxing this constraint.

We discuss the consequences of relaxing this minimum stellar particle threshold on the apparent fractional ICL contribution from a cluster’s own satellites in Section 5.2.3, but confirm the progenitor analysis presented in Section 4.2 to otherwise be insensitive to the exact threshold employed. The predicted peak ICL (and BCG + ICL) progenitor masses obtained by repeating the analysis presented in Figure 7 for TNG100, The300 GS-7K, and Hydrangea with a minimum galaxy stellar particle threshold of unity are functionally identical to those presented in Table 3.

5.2.2 BCG (and BCG + ICL) mass fractions↩︎

We note some scatter between the median BCG stellar mass fractions of the four simulations in Figure 1, most conspicuously the comparatively low median BCG fraction found for The300 GS-7K (of 0.19, compared to 0.35, 0.43, and 0.29 for Hz-AGN, TNG100, and Hydrangea respectively). Though less prominent, the distribution of BCG mass fractions in TNG100 also appears mildly skewed towards larger values – with the median BCG fraction in TNG100 coinciding with the maximum of the The300 GS-7K sample. This BCG fraction scatter is largely responsible for the similar scatter also seen in the typical BCG + ICL mass fractions, as the ICL fractions of the different simulations are broadly consistent. We consider this scatter chiefly the joint consequence of differing typical cluster assembly times and differing dependencies on accretion for assembling stellar mass between the different cluster samples.

The The300 GS-7K sample clusters typically have lower DM halo assembly redshifts compared to those from the other samples (with a median \(z_{50}\) value of 0.53, compared to 0.62, 0.82, and 0.77 for Hz-AGN, TNG100, and Hydrangea respectively; see Appendix 8). A more recent assembly redshift suggests a larger fraction of the total cluster mass to be newly introduced – brought in with satellite galaxies that have yet had less opportunity for disruption – and so a larger fraction of the cluster’s total stellar mass should be retained in the satellite component. As such, a correlation between assembly redshift and BCG + ICL mass fraction appears intuitive, and prior studies have noted such trends already (e.g. [158], [24]; see also Appendix 8). Furthermore, we highlight that the The300 GS-7K clusters are systematically more massive at \(z\approx0\) than those from the other samples and prior studies have also reported broad trends towards lower central stellar mass fractions in increasingly massive haloes (e.g. [82], [14]).

We also attribute part of the dispersion in typical BCG stellar mass fractions between the simulations to differing amounts of late-time central star formation. We note significantly more “in-situ” BCG + ICL star formation occurring at \(z<3\) in TNG100 and Hydrangea as compared to Hz-AGN or The300 GS-7K. This star formation – occurring either within the central galaxy or directly into the ICL component – often contributes \(>20\) per cent of all \(z\approx0\) BCG + ICL stellar particles in both TNG100 and Hydrangea (mean fractions 0.26 and 0.25 respectively), as opposed to typically \(<10\) per cent in both Hz-AGN and The300 GS-7K (mean fractions 0.08 and 0.05 respectively). These in-situ stars are also generally much more centrally concentrated by \(z\approx0\) in TNG100 as compared to Hydrangea. We leave a more thorough investigation of this in-situ star formation to a follow-up paper (though see Section 5.2.3) – for further relevant discussion see e.g. [82], [148], [14], and references therein.

5.2.3 ICL contribution from a cluster’s own satellites↩︎

We consider the lower mean fraction of ICL stellar mass liberated within a cluster’s own halo from satellites joining to-be clusters after \(z\sim3\) in The300 GS-7K as compared to Hz-AGN (0.29 vs 0.62; Figure 4) to partly stem from the comparatively lower resolution of The300 GS-7K, in combination with the minimum stellar particle count of 51 we impose when identifying galaxies (Section 2.2.2). The stellar particles of poorly-resolved galaxies below this threshold are relabelled as “intergalactic”, which when joining clusters are therefore classified as pre-processed ICL. When this minimum stellar particle count for galaxies is relaxed (acknowledging poorly-resolved satellites and so slightly reducing the obtained ICL mass; Section 5.2.1) the mean fraction of \(z=0\) ICL stellar mass liberated from a cluster’s own satellite galaxies in The300 GS-7K rises from 0.29 to 0.38 (and from 0.49 to 0.55 for BCG + ICL stars; the mean ICL and BCG + ICL fractions from tracked satellites in TNG100 and Hydrangea are unaltered by this change). The clusters of the The300 GS-7K sample are also systematically more massive than those of the other samples, and this may also play a role in the typically smaller proportional ICL contribution by a cluster’s own satellites in The300 GS-7K, as the significance of ICL pre-processing is expected to grow with cluster mass ([23]).

Rather than a direct consequence of differing simulation resolution or a discrepancy in mass between the cluster samples, we primarily attribute the lower mean fractions of \(z\approx0\) ICL stellar mass liberated from a cluster’s own satellite galaxies in TNG100 and Hydrangea (0.44 and 0.27 respectively) to a non-negligible ICL contribution from the apparent “in-situ” formation of stars directly into the ICL component in these two simulations. The mean fractions of \(z\approx0\) ICL stellar mass that we link to this channel (i.e. \(z\approx0\) ICL stellar particles first seen within the cluster at \(z<3\) but outside any resolved galaxy) are \(0.14\) and \(0.24\) in TNG100 and Hydrangea respectively (and \(\sim0.01\) in both Hz-AGN and The300 GS-7K) – though we caution that these determined fractions are highly sensitive to the galaxy definition imposed. If this “in-situ” ICL component is ignored, then the mean fractions of the remaining \(z\approx0\) ICL mass contributed by tracked satellites are \(0.51\) and \(0.36\) in TNG100 and Hydrangea respectively. We defer a more comprehensive investigation and inter-simulation comparison of the contribution from this alternate ICL assembly channel to a forthcoming paper (Brown et al. in prep.) – for further relevant discussion see e.g. [159], [148], [15], and references therein.

5.2.4 Liberated fraction of stellar mass↩︎

We find broad agreement between the fitted functions for \(f_\text{lib-ICL}(M_{*})\) from TNG100, The300 GS-7K, and Hydrangea shown in Figure 5. This is despite the different models and methods employed and resolution achieved by each simulation, and also that we do not control for any differences in typical satellite infall time, infall mass ratio, satellite orbital parameters, or typical cluster assembly time between the different simulations.

However, we do note that the fitted function for \(f_\text{lib-ICL}(M_{*})\) from Hz-AGN diverges noticeably from the other simulations, with generally higher mean \(f_\text{lib-ICL}\) values, particularly at low infall masses. This is likely to be a product of the high stellar-to-halo mass ratios ([73]; see also Appendix 8) and low galaxy compactnesses noted in Hz-AGN, with the median (and 16th-84th percentile scatter in) stellar-half mass radii for \(M_{*}\in[1,5]\times10^{10}\,\)M satellites on first cluster infall at \(z\approx0.5\) being \(5.5^{+1.6}_{-0.9}\) kpc in Hz-AGN, as compared to \(2.6^{+2.1}_{-0.6}\), \(3.6^{+0.8}_{-0.7}\), and \(2.3^{+1.9}_{-1.1}\) kpc in TNG100, The300 GS-7K, and Hydrangea respectively. Reduced galaxy compactness and under-massive DM haloes correspond to shallower potentials and so increased vulnerability to tidal stripping (e.g. [160], [161], [155]), from which the higher typical \(f_\text{lib-ICL}\) values for Hz-AGN satellites logically follow.

5.2.5 Infalling satellite mass functions↩︎

Within estimated uncertainty the infalling satellite mass function fits (Figure 6) for the different simulations appear in reasonable agreement, though a slight dearth of low-mass infalling satellites (\(M_{*}\lesssim10^{9.5}\,\)M) is noted for Hz-AGN. We associate this with the weak supernova feedback implementation of Hz-AGN ([73], [75]) and so tentatively link this low-mass satellite deficit with the trace excess of \(M_{*}\sim10^{10.5}\,\)M satellites also seen for Hz-AGN in Figure 6. This slight discrepancy at low-masses could also stem from the other three simulations somewhat over-predicting the abundance of low-mass satellites, with minor excesses of low-mass galaxies previously reported in both TNG and Hydrangea at \(z\lesssim1\) ([82], [162]).

5.2.6 Main progenitors of ICL stars↩︎

Despite the broad agreement between the simulations for the predicted peak contributors of ICL stars, two standout discrepancies are apparent in Figures 7 and 8: Hz-AGN predicts a slightly lower peak ICL contributor mass, and Hydrangea predicts a more significant ICL contribution from \(M_{*}\sim10^{9}-10^{9.5}\,\)M satellites as compared to the other simulations.

The prediction of less massive peak ICL (and BCG + ICL) contributors in Hz-AGN follows directly from the infalling satellite mass function fit for Hz-AGN having its knee at a lower mass than the other simulations (Figure 6; see also Table 2). Although clear trends appeared in the obtained fits for \(f_{\text{lib-ICL}}(M_{*})\) (Figure 5), the average value for \(f_{\text{lib-ICL}}\) never strayed below \(\sim0.05\) at any infall mass for any simulation and also generally did not rise above \(\sim0.5\) (with the fits obtained for \(M^{z=0}_{*,\textrm{Tot}}(M_{*})/M_{*}\) likewise tightly constrained; Appendix 11). This comparatively modest variation cannot significantly influence the shapes of the curves seen in Figure 7, which instead chiefly follow from the competition between rising scarcity with increasing stellar mass, and more massive galaxies each having more stars to potentially contribute. Correspondingly, it can be seen (by comparison between Tables 2 and 3) that both the peak ICL and BCG + ICL contributor infall masses predicted by every simulation are never more than \(\sim0.2\) dex from the value of \(M_{k}\) for the corresponding infalling satellite mass function fit, with the increased efficacy of tidal stripping at lower satellite masses serving only to skew these peaks towards slightly lower masses than would have been obtained by just considering the mass function alone.

The prediction of a more significant ICL contribution from low-mass galaxies (\(M_{*}\lesssim10^{10}\,\)M) in Hydrangea also follows largely from the infalling satellite mass function, which has a particularly steep low-mass slope in Hydrangea (\(\alpha\approx-1.38\); Table 2). A steeper mass function slope at low-masses will result in low-mass galaxies composing a higher fraction of the total stellar mass budget, yielding a prediction for \(f_\text{ICL}(M_{*})\) (equation 2 ) increasingly skewed towards lower-masses, with a heavier low-mass tail [36]⁸. The fit to \(f_\text{lib-ICL}(M_{*})\) for Hydrangea also levels out at a slightly lower value (\(\sim0.05\)) at high masses compared to the other simulations, yielding a slightly lower ICL contribution from more massive satellites, and so further enhancing the relative ICL contribution of less massive galaxies.

5.3 Caveats and limitations↩︎

The predicted ICL contributions of differingly massive satellites depicted in Figures 7 and 8 average over the intrinsic stochasticity of cluster assembly. Although we predict a significant portion of the aggregate ICL of many clusters to be associated with \(M_{*}\gtrsim10^{11}\,\)M progenitor objects, few objects this massive would be anticipated in the assembly history of any individual cluster (Figure 6), with each of these massive objects expected to individually make a substantial ICL mass contribution in absolute terms, even if this is only a small fraction of their infall stellar mass (Figure 5). As such, though we predict the peak contributors of ICL stars to in general be roughly Milky Way mass galaxies, individual clusters may exist in which the dominant ICL progenitor mass deviates towards higher masses, where the content of the ICL will be dominated by the contribution from only a small number of very massive objects. Alternatively, should a cluster’s mass function during assembly be skewed towards lower masses (with a steep low-mass slope), similar to that we find for the Hydrangea sample or as was recently suggested for Hydra I by [163], then for that individual cluster the apex ICL contributor mass may be pushed to lower masses.

For the analysis shown in Figures 7 and 8 we consider only satellite galaxies which fall into the studied clusters. We do not recursively link the pre-assembled diffuse stellar content of groups or smaller clusters that merge into the studied clusters with prior satellites of those groups or clusters, and instead entirely neglect this “pre-processed” ICL component for this analysis. As groups should be assembled from on average less massive progenitors than clusters, we anticipate that if such a recursive analysis were performed the predicted peak ICL (and BCG + ICL) contributor mass would fall slightly. As these “pre-processed” diffuse stars will typically be only weakly bound to the haloes in which they are delivered to clusters (as compared to still-galactic stars), they should be stripped from these halos soon after cluster infall, before the stripping of galactic stars begins ([160]), and so we would expect these “pre-processed” ICL stars to be preferentially deposited at large cluster centric radii. That these “pre-processed” ICL stars should also originate from typically less massive (and so typically more metal-poor; e.g. [164]) galaxies than the more centrally concentrated ICL assembled “in-place” within a cluster from its own satellite galaxies may then contribute to observed radial colour trends in the ICL (e.g. [7]). We directly address the ICL contribution from this “pre-processed” channel and interrogate the origins of these radial ICL trends further in the forthcoming companion paper (Brown et al. in prep.).

6 sec:Conclusions↩︎

In this study we have explored the ICL assembly of low-mass galaxy clusters (\(M_{178c}\sim10^{14}-10^{15}\) M at \(z\approx0\)) using four markedly different hydrodynamical simulations: Hz-AGN, TNG100, The300 GS-7K, and Hydrangea (using samples of \(10-16\) clusters per simulation). We employed a consistent approach for identifying ICL stars, allowing us to probe for significant inter-simulation discrepancies unobstructed by differing ICL extraction methodologies. We have investigated how the typical ICL stellar mass fractions of these clusters compare at \(z\approx0\), how these have evolved since \(z\sim3\), as well as how this evolution is influenced by the adopted methodology for separating the ICL and the BCG, and have also tracked the stars of satellite galaxies joining clusters during this same period in order to quantify the different typical ICL contributions made by galaxies with differing stellar masses at cluster infall. Our main findings can be summarised as follows:

1. The \(z\approx0\) ICL stellar mass fractions of all four simulations are broadly consistent (Figure 1). When uniformly subject to our fiducial ICL definition (based on local stellar density), the ICL fractions of almost every considered cluster are confined to \([0.10,0.19]\) – with median ICL fractions of 0.13, 0.12, 0.18, and 0.17 in Hz-AGN, TNG100, The300 GS-7K, and Hydrangea, respectively. The BCG + ICL stellar mass fractions of each simulation also all span the same range (between \(\sim0.3\) and \(\sim0.7\)), with inter-simulation differences in the average BCG + ICL fraction following directly from differing typical BCG mass fractions, thought to primarily stem from differing typical cluster assembly times and masses between the four samples.
   

2. Whether the average ICL stellar mass fraction remains constant or falls with increasing redshift is determined by the adopted ICL definition, but all four simulations reproduce the same redshift evolution trend when subject to a consistent definition (Figure 3). When considering the main progenitors of \(z\approx0\) clusters, our fiducial ICL definition based on local stellar density relative to the cosmic average matter density, BCG-ICL splits based on spherical aperture cuts that scale with DM halo overdensity radius or BCG + ICL stellar half-mass radius, as well as kinematic BCG-ICL separations all yield no significant redshift evolution in the average ICL stellar mass fraction (at least for \(z\lesssim2\)). Identifying ICL stars using fixed stellar density thresholds or BCG-ICL splits based on spherical apertures of fixed radius (without controlling for evolving central galaxy mass or size) instead yield typical ICL fractions that fall to approximately zero by \(z\sim3\).
   

3. All four simulations agree that a galaxy with less stellar mass on cluster infall can be expected to contribute a greater fraction of its stars to the \(z\approx0\) ICL than a more massive infalling galaxy (Figure 5), when considering both stars liberated by tidal stripping and those ejected during galaxy mergers (and when marginalising over cluster assembly time, satellite orbital parameters, satellite infall time, and satellite infall mass-ratio). However, even the most massive galaxies that join clusters are still predicted to typically eventually contribute \(\gtrsim5\) per cent of their associated stellar mass to the ICL component.
   

4. All four simulations predict the peak contributors of ICL stars in \(M_{178c}\lesssim10^{15}\) M clusters at \(z\approx0\) to be roughly Milky Way mass galaxies (infall stellar mass \(\sim10^{11}\) M; Figure 7), with the peak contributors of BCG + ICL stars only slightly more massive (by \(\sim0.1\) dex) than for the ICL alone. These peak contributor masses are primarily controlled by the infalling satellite mass function during cluster assembly, and are never more than \(\sim0.2\) dex below the characteristic mass of the fitted mass functions (Figure 6).
   

5. All four simulations predict \(\gtrsim50\) (\(\gtrsim80\)) per cent of satellite galaxy sourced ICL stars to be contributed by satellites with infall stellar masses \(>10^{10}\) M (\(>10^{9}\) M; Figure 8). The corresponding satellite infall stellar mass thresholds for the combined BCG + ICL are \(10^{10.5}\) M (for \(>50\) per cent) and \(10^{10}\) M (for \(>80\) per cent).

Our results suggest that the conflicting findings of prior theoretical studies concerning both typical ICL stellar mass fractions and the evolution of these with redshift may entirely stem from the differing ICL extraction methodologies employed by these studies. Subject to a consistent ICL definition, we find no significant discrepancies between the four simulations concerning ICL mass fractions at \(z\approx0\); ICL mass fraction evolution in the progenitors of \(z\approx0\) low-mass clusters back to \(z\sim2\); the relative ICL contributions of differingly massive satellites on cluster infall; or the main progenitor objects of ICL stars.

What remains as yet under-explored is how the considered simulations differ in their predictions for contributions to the ICL via assembly channels beyond just the liberation of stars from cluster satellite galaxies; from the pre-processed ICL accreted along with groups and smaller clusters during cluster assembly, and from stars that form “in-situ” either in the central galaxy or directly into the ICL component. We intend to investigate these alternative ICL assembly channels, as well as the resulting implications for cluster-centric radial trends in stellar age and metallicity, in a forthcoming companion paper (Brown et al. in prep.).

Acknowledgements↩︎

The authors thank the other members of the NottICL Group – and in particular Jesse B. Golden-Marx and Harry Gully – for helpful discussions and comments; and also thank Emanuele Contini for their careful reading of the original manuscript and for their constructive comments which have helped to improve the quality and clarity of the presented work. H. J. Brown thanks Dylan Nelson for their guidance on handling TNG100 data.

H. J. Brown acknowledges support from the UK Science and Technology Facilities Council (STFC) under grant ST/Y509437/1. F. .R. Pearce and N. A. Hatch acknowledge support from the UK STFC under grant ST/X000982/1. Y. M. Bahé acknowledges support from UK Research and Innovation through a Future Leaders Fellowship (grant agreement MR/X035166/1) and financial support from the Swiss National Science Foundation (SNSF) under project “Galaxy evolution in the cosmic web” (200021_213076). J. Butler and N. A. Hatch acknowledge support from the Leverhulme Trust through a Research Leadership Award. W. Cui thanks Comunidad de Madrid for the Atracción de Talento fellowship no. 2020-T1/TIC19882 and Agencia Estatal de Investigación (AEI) for the Consolidación Investigadora Grant CNS2024-154838; he further acknowledges the Project PID2024-156100NB-C21 financed by MICIU/AEI /10.13039/501100011033/FEDER, EU and ERC: HORIZON-TMA-MSCA-SE for supporting the LACEGAL-III (Latin American Chinese European Galaxy Formation Network) project with grant number 101086388 and the science research grants from the China Manned Space Project. A. Knebe is supported by project PID2024-156100NB-C21 financed by MICIU /AEI/10.13039/501100011033 / FEDER, UE, and further thanks Oasis for supersonic.

This work made use of NumPy [165], Matplotlib [166], SciPy [167], h5py [168], pyGAM [169], Scikit-learn [170], and Astropy [171].

The Horizon-AGN simulation was granted access to the HPC resources of CINES under allocations 2013047012, 2014047012 and 2015047012 made by GENCI, and made use of the Infinity cluster, hosted by the Institut d’Astrophysique de Paris. We warmly thank S. Rouberol for running it smoothly. The IllustrisTNG simulations were undertaken with compute time awarded by the Gauss Centre for Supercomputing (GCS) under GCS Large-Scale Projects GCS-ILLU and GCS-DWAR on the GCS share of the supercomputer Hazel Hen at the High Performance Computing Center Stuttgart (HLRS), as well as on the machines of the Max Planck Computing and Data Facility (MPCDF) in Garching, Germany. The authors acknowledge The Red Española de Supercomputación for granting computing time for running the hydrodynamical and DMO simulations of the The Three Hundred galaxy cluster project in the Marenostrum supercomputer at the Barcelona Supercomputing Center and Cibeles Supercomputers through various RES grants. The The Three Hundred HD hydrodynamic simulations (7K and 15K runs) were performed also on the DIaL3 – DiRAC Data Intensive service at the University of Leicester through the RAC15 grant: dp235, and on the Niagara supercomputer at the SciNet HPC Consortium. DIaL3 is managed by the University of Leicester Research Computing Service on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The DiRAC service at Leicester was funded by BEIS (ST/K000373/1), UKRI, STFC capital funding, and STFC operations grants (ST/K0003259/1). DiRAC is part of the UKRI Digital Research Infrastructure. SciNet [172] is funded by Innovation, Science and Economic Development Canada; the Digital Research Alliance of Canada; the Ontario Research Fund: Research Excellence; and the University of Toronto. This work used the DiRAC@Durham facility managed by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital grants ST/K00042X/1, ST/P002293/1, ST/R002371/1 and ST/S002502/1, Durham University and STFC operations grant ST/R000832/1. DiRAC is part of the National e-Infrastructure.

Data Availability↩︎

The data products of the Horizon-AGN simulation are available upon reasonable request through the collaboration’s website: https://www.horizon-simulation.org/. The data products of the TNG100 simulation are publicly available through the IllustrisTNG project’s website: https://tng-project.org/. The data products of The Three Hundred project are available upon reasonable request through the collaboration’s website: https://www.the300-project.org. The data products of the Hydrangea simulations are publicly available through the Leiden Observatory Science Data Repository: https://ftp.strw.leidenuniv.nl/bahe/Hydrangea. Further data products generated from this work are available upon reasonable request from the corresponding author.

7 Recombining fragmented galaxies↩︎

Not all structures identified by Subfind in TNG100 or Cantor in Hydrangea correspond to galaxies, even when considering only those containing stars. As these codes consider any local density peak a seed for potential sub-structure, dense clumps of stars and gas can end up segregated from the bulk of their host galaxy as their own “sub-halo”. As these fragments isolate some of the densest, baryon-dominated regions in massive galaxies, they can not only non-negligibly lower the apparent stellar mass of their host galaxy if sub-structure is neglected (particularly if several are carved away from the same galaxy), but can themselves masquerade as extremely compact “galaxies” in structure catalogues. In TNG100, these fragments can have stellar masses as high as \(\gtrsim10^{9.5}\) M, with sizes of order unity kpc, negligible DM content, and often very high specific star formation rates. These fragments are also long-lived (commonly persisting for several Gyrs) and can multiply by fragmenting further into smaller and smaller objects, such that at times up to \(\sim20\) per cent of the “sub-haloes” with stellar mass between \(10^{8.5}-10^{9.5}\) M in and around some of the TNG100 clusters we consider belong to this class of object, with similarly high abundances reported in Hydrangea previously by [108]. An in-depth investigation into this fragmentation is beyond the scope of this study (for further relevant discussion see e.g. [81], [173], and [174]; see also [108]) – but this behaviour is undesirable for our use case and we seek to correct for this fragmentation when constructing our galaxy catalogues.

In Hydrangea we identify these “anti-hierarchically” formed fragments following the procedure described in [108], in which these objects are labelled “spectres”. Once our procedure for assigning stars to galaxies is nominally complete (prior to neglecting structures with insufficient particle counts; see Section 2.2.2), if the stars of a provisional “galaxy” identified within a sub-halo flagged as a “spectre” are at least partially enclosed by the “galaxy” identified in the spectre’s parent sub-halo, then the “spectre” galaxy is absorbed into its parent (with relevant merger tree links either redirected to the parent or cut as appropriate; see Section 2.3). In the event a spectre has no parent sub-halo, or the spectres’ galaxy stars are not enclosed within the parent galaxy (and similar checks performed if applicable using the parent’s parent sub-halo and so forth also fail) we retain the spectre-flagged object as a prospective galaxy.

A very similar procedure is followed in TNG100, for which the SubhaloFlag field is provided to identify objects thought not of cosmological origin (prototypically clumps formed within galaxies through baryonic processes; [81]). However, we note both a non-negligible number of apparent false positives (i.e. entire infalling satellites flagged as not of cosmological origin despite being otherwise unremarkable) and also a considerable number of compact, DM deficient sub-structures archetypical of the galaxy fragmentation we seek to remedy that are not caught by SubhaloFlag. These false positives can stem from galaxies with tumultuous early histories, such that they are briefly divided into two sub-haloes (one of which briefly drops beneath the DM mass fraction threshold of 0.8 of SubhaloFlag) before recombining. The considerable number of missed fragments we note conceivably stems from our focus on galaxy clusters, which can have very large and distended FoF groups, in combination with the criteria used for SubhaloFlag explicitly omitting objects first detected outside their host FoF group’s virial radius [81]. We consequently do not rely on SubhaloFlag alone to identify these fragments in TNG100 and append our own scheme based on that used in [108]. We flag as fragments (to be handled as described above for the “spectres” of Hydrangea) TNG100 sub-haloes that are first seen in one of the coarsely spaced snapshots we consider (see Section 2.3.1) as a child sub-halo, with a dark-matter fraction less than 0.5, and that are either flagged by SubhaloFlag; are \(\geq75\) per cent composed of particles formerly in their parent sub-halo’s progenitor in the prior coarsely spaced snapshot; or that are sub-structure of a newly-flagged fragment.

8 Cluster samples↩︎

8.1 Cluster properties↩︎

We identify every cluster from the four samples we employ in Tables 4, 5, 6, and 7 (for Hz-AGN, TNG100, The300 GS-7K, and Hydrangea respectively) and present for these a selection of cluster properties. We identify each cluster in Hz-AGN by the AdaptaHOP ID of its central galaxy in snapshot 761; in TNG100 by the ID/index of its FoF group in the structure catalogue for snapshot 99; and in The300 GS-7K and Hydrangea by the ID of the simulation box it is centrally located in.

In Tables 4-7 we present for each cluster the value (at \(z\approx0\)) of \(r_{178c}\), centring on the stellar centre of mass of the identified central galaxy (see Section 2.2.3); the total mass enclosed within this radius at \(z\approx0\) (\(M_{178c}\)); the fractions of this total mass contributed by DM (\(f_\text{DM}\)) and by stellar particles (\(f_{*}\)); as well as the fractions of this stellar mass at \(z\approx0\) contributed by the BCG, satellite galaxies, and the ICL (\(f_{*\text{,\,BCG}}\), \(f_{*\text{,\,Sat}}\), and \(f_{*\text{,\,ICL}}\), respectively; following our fiducial ICL and galaxy identification scheme described in Section 2.2). We also give estimates for the redshifts when each cluster first assembled 50 and 90 per cent of its \(z\approx0\) DM mass (\(z_{50}\) and \(z_{90}\) respectively) as found by fitting monotonic cubic spline functions to the total DM mass at \(r<r_{178c}(z)\) against lookback time for the \(\sim1\) Gyr spaced snapshots we employ in this work (with these fits anchored to the \(z\approx0\) masses) and interpolating. We additionally give stellar mass analogues to the “magnitude gap” (e.g. [30]) defined here as \(M_{1X}=\text{log}_{10}(M_{*,1}/M_{*,X})\) where \(M_{*,1}\) is the BCG stellar mass and \(M_{*,X}\) the stellar mass of the \((X-1)\)th most massive satellite galaxy within \(r_{178c}\). We present both \(M_{12}\) and \(M_{14}\) values for each cluster at \(z\approx0\). To contextualise these “stellar magnitude gap” values, we note the median values of \(M_{12}\) and \(M_{14}\) across all four cluster samples to be \(\sim0.6\) and \(\sim1.0\) respectively. The DM halo concentration values we present (defined here as \(c_\textrm{Halo-178}\equiv r_{178c}/r_s\)) are computed by a similar methodology to that employed by [107]: centring on the density peak of each cluster’s DM halo (as opposed to the BCG centre as elsewhere in our analysis), we fit an NFW profile with scale radius \(r_{s}\) to the spherically averaged DM density in 40 equal log-width radial bins between \(0.05\times r_{178c}\) and \(r_{178c}\). The offsets between the BCG stellar centre of mass and the centre of mass of all DM enclosed by \(r_{178c}\) at \(z\approx0\) are also given (as a percentage of \(r_{178c}\)).

The Hz-AGN cluster sample we employ for this work has significant overlap with that used in . Slight differences between some of the values presented in Table 4 and those values given in table 1 of are due to subtle differences in methodology – chiefly that in potential centres were used for the position of each galaxy, whereas in this work we use stellar centre of mass (slightly shifting the centres used when determining \(r_{178c}\) and \(M_{178c}\)). Some of the values we present for the clusters of the Hydrangea sample in Table 7 slightly differ from those given in table A1 of [107] for similar reasons.

Prior studies have reported possible ICL fraction trends with cluster mass (e.g. [22], [21]), with cluster dynamical state (e.g. [30], [24]), and with the concentration of a cluster’s DM halo (e.g. [175], [15]). Though we investigated potential trends between the ICL stellar mass fraction (under our fiducial ICL definition) and cluster mass, DM halo concentration, as well as various proxy measures of dynamical state (\(z_{50}\), \(z_{90}\), \(M_{12}\), \(M_{14}\), and the offset between the BCG and DM halo centres of mass), we could not identify any significant trends above the intra- and inter-simulation scatter for the restricted mass ranges and modest cluster sample sizes we employ in this study. However, we did note trends between these same cluster properties and the combined BCG + ICL stellar mass fractions – with typically larger BCG + ICL stellar mass fractions seen for clusters that are less massive, that have more concentrated DM haloes, that have smaller BCG-DM centre of mass offsets, that have larger stellar-mass magnitude gap values, and that assembled earlier (consistent with the prior findings of e.g. [24] and [30]).

8.2 Extreme objects↩︎

We highlight cluster 13 in the Hz-AGN sample as a particularly extreme object – in the early stages of two simultaneous major cluster mergers at \(z\approx0\). Clusters 0 and 1 from TNG100 (i.e. the two most massive FoF groups in the simulation box at \(z\approx0\)) are also notable for being highly disturbed, with large populations of massive satellites predominately located to one side of the identified BCG at \(z\approx0\) such that the “central” galaxy in both instances appears noticeably offset from centroid of the local galaxy distribution. Cluster 0 from TNG100 is also noteworthy for having non-negligible late time star formation in its central galaxy – with the median stellar particle age at \(r<0.02\times r_{178c}\) in this cluster being \(<7.5\) Gyr, as compared to \(>10\) Gyr in every other considered TNG100 cluster. This is similar to cluster 310 from the The300 GS-7K sample where the median stellar particle age at \(r<0.02\times r_{178c}\) is \(<4\) Gyr (with 20 per cent of this cluster’s BCG + ICL stars having formed in-situ in either the BCG or ICL). We also highlight TNG100 cluster 3 as unusual, with \(\sim50\) per cent of the \(z\approx0\) BCG + ICL stars in this cluster linked to either BCG or ICL in-situ star formation, significantly higher than any other TNG100 cluster. Among the Hydrangea cluster sample, we highlight CE-4 and CE-15 as extreme objects – with the \(z=0\) central in CE-4 seemingly part of an in-spiralling pair along with a particularly massive satellite galaxy; and with a “satellite” galaxy having just joined the halo of CE-15 in the final \(\sim1\) Gyr before \(z\approx0\) that actually has slightly more stellar mass than what we consider the central galaxy for our analysis (hence the negative \(M_{12}\) value seen for this cluster in Table 7) – with this very slightly more massive object not instead ruled the BCG as it was not found within the sub-halo Cantor ruled the cluster’s “central structure” at \(z\approx0\) (see Section 2.2). We confirm that excluding these extreme objects from our analysis does not considerably alter any of the presented findings.

9 Mass fraction evolution for cluster mass progenitors↩︎

In Sections 3.2 and 3.3, we present how the typical ICL (and BCG + ICL) stellar mass fractions of the four samples of simulated \(z\approx0\) clusters and their \(z>0\) progenitor structures evolve as a function of lookback time. In Figures 9 and 10 (counterparts to Figures 2 and 3 respectively) we present alterative versions of those analyses confined to cluster-mass structures (\(M_{178c}\geq9\times10^{13}\) M) at all times. In Figures 9 and 10 unfilled markers indicate median mass fraction values found for samples of less than 5 objects, for which we also no longer indicate 16th-84th percentile dispersion via the shaded regions.

We note that in Figure 9 and most panels of Figure 10 there is the suggestion that the typical ICL mass fraction of the (small number of) Hydrangea progenitors already of cluster mass begins to rise for lookback times \(\gtrsim8\) Gyr. This chiefly stems from the three structures in the Hydrangea sample that reach cluster-scale mass first (identifiers CE-13, CE-16, and CE-18; see Table 7) having some of the highest ICL mass fractions of the sample at these early lookback times, and so the median ICL fraction is pulled down as additional structures pass the imposed minimum mass threshold. We suspect the high early ICL fractions of these structures linked to the significant amounts of in-situ ICL star formation we note occurring in Hydrangea at \(z>1\), such that an appreciable fraction of the \(z\approx0\) ICL mass of every Hydrangea cluster considered is attributable to this high redshift in-situ star formation. We leave further discussion of this in-situ star formation to the forthcoming companion paper (Brown et al. in prep.).

10 Data and fitting for progenitor galaxies of ICL (and BCG + ICL) stars↩︎

In Figures 5 and 6 we include only curves fitted to the data from each simulation while omitting the original data fitted to for clarity of presentation. We present those data here.

In Figure 11 (counterpart to Figure 5) we present for each simulation individually the mean fraction of associated stellar mass liberated to the ICL (or BCG + ICL) by \(z\approx0\) as a function of satellite infall stellar mass, as found for 0.5 dex wide rolling infall stellar mass bins (bin step 0.15 dex) and to which we fit a cubic spline function, linearly extrapolated beyond the range of masses used for fitting (see Section 4.2.1 for further details). For clarity of presentation, only every other bin is shown in Figure 11. The error bars and shaded regions indicate the dispersion (16th-84th percentile) in mean values and functions fitted among the \(10^{4}\) bootstrap resamples.

In Figure 12 (counterpart to Figure 6) we present for each simulation individually a histogram of satellite infall stellar masses (0.25 dex bin width; counts normalised to mean counts per cluster) to which we fit a [147] function (Equation 1 ). Before fitting, we merge any bins containing fewer than five objects across the entire cluster sample for each simulation with the adjacent lower-mass bin, and when fitting exclude bins including masses that corresponding to poorly resolved galaxies (\(<100\) stellar particles; see Section 4.2.2 for further details). The error bars indicate the dispersion (16th-84th percentile) in bin counts among the \(10^{4}\) bootstrap resamples, and the shaded regions show 16th-84th percentile confidence intervals for the fitted curves (based both on parameter fit uncertainty and the dispersion in fit parameters across the bootstrap resamples).

11 Post-infall satellite star formation↩︎

In Figure 13 we present for each simulation individually the mean ratio between the total stellar mass associated with a satellite galaxy by \(z\approx0\) (see Section 2.3.2) and satellite infall stellar mass as a function of infall mass, considering the same population of galaxies and following the same binning procedure as the analyses depicted in Figures 5 and 11. We fit to these mean mass ratios a cubic spline function in the same manner as our fitted \(f_{\mathrm{lib}}(M_{*})\) functions (see Section 4.2.1 for details), including linear extrapolation beyond the range of masses used for fitting, but with the output constrained to \([1,\infty]\) rather than \([0,1]\) and no longer disregarding low-mass satellites when fitting. For clarity, only every other bin is shown in Figure 13. The error bars and shaded regions indicate the dispersion (16th-84th percentile) in mean values and functions fitted among the \(10^{4}\) bootstrap resamples.

The fitted functions for each of the simulations all have the same fundamental shape – with a peak at intermediate infall masses (\(10^{9.5}-10^{10}\,\)M), falling back to unity at both lower and higher masses. We attribute the differing heights and precise positions of this peak in each panel to the differing sub-grid physics prescriptions and implementations of the different simulations as well as factors we ignore for this analysis (most pertinently differing typical cluster assembly times).

12 Cumulative fractional ICL (and BCG + ICL) contributions direct from simulations↩︎

In Figure 8 we present predictions for the cumulative contribution to the ICL mass sourced from satellite galaxies as a function of satellite infall stellar mass. These predictions are based on fitted functions (as per Equation 2 ) and include extrapolation to unresolved low satellite masses. An equivalent analysis can be obtained directly from the simulations without fitting or extrapolation – which we give in Figure 14, where the cumulative fraction of ICL mass from satellite galaxies (for the aggregate ICL of each simulation’s entire cluster sample) is shown as a function of satellite infall stellar mass (in red). The equivalent predicted curves (based on fitting and extrapolation) are reproduced from Figure 8 for comparison (in blue), normalized relative to the raw data so over-/under-prediction of total contributed ICL mass can be assessed. The shaded regions indicate the dispersion (16th-84th percentile) of equivalent analyses for \(10^{4}\) bootstrap resamples of the infalling satellite population, all also normalized relative to the main analysis of the raw data. The analysis presented in Figure 14 is similar to that presented for Hz-AGN alone in figure 7 of .

Significant broadening of the blue shaded regions can be observed in every panel of Figure 14 for stellar masses beyond \(\sim10^{11}\) M – a consequence of the significant ICL (and BCG + ICL) mass contribution each individual object this massive is expected to make. As such, even small changes to the infalling satellite mass function (encompassed by the shaded regions shown in Figures 6 and 12) can alter the total predicted ICL mass by as much as \(\gtrsim50\) per cent.

References↩︎

---

1. E-mail: Harley.Brown@nottingham.ac.uk↩︎

2. https://www.horizon-simulation.org/↩︎

3. https://www.tng-project.org/↩︎

4. https://the300-project.org/↩︎

5. For this work, \(r_{\Delta c}\) is the radius within which the average DM density is \(\Delta\) times the cosmic critical density (where \(\Delta\) is either 178 or 200), and \(M_{\Delta c}\) is the total mass enclosed by \(r_{\Delta c}\).↩︎

6. We largely adopt the terminology of Subfind and refer to the gravitationally-bound, DM+baryonic structures identified by Subfind, AHF, and Cantor as “sub-haloes”, with one exception: in isolation, we refer to the self-bound DM+baryons of a cluster’s main halo (excised of embedded sub-structure) as the cluster’s “central structure” to avoid implying this structure to be a sub-structure.↩︎

7. We continue to refer to the diffuse stellar component of even the sub-cluster scale progenitor structures as ICL.↩︎

8. Specifically see the accompanying interactive plots they provide, available at https://garrethmartin.github.io/interactive-profiles-ICL/index.html#stripping↩︎