Assembly of the intracluster light in the
HORIZON-AGN simulation

Garreth Martin, Frazer R. Pearce, Nina A. Hatch, Yannick M. Bahé, Yohan
Dubois

► To cite this version:

Garreth Martin, Frazer R. Pearce, Nina A. Hatch, Yannick M. Bahé, Yohan Dubois. Assembly of the intracluster light in the HORIZON-AGN simulation. Monthly Notices of the Royal Astronomical Society, 2024, 534, pp.431-443. 10.1093/mnras/stae2084. insu-04851504

HAL Id: insu-04851504

https://insu.hal.science/insu-04851504v1

Submitted on 21 Dec 2024

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
ROYAL ASTRONOMICAL SOCIETY

MNRAS 534, 431–443 (2024)

https://doi.org/10.1093/mnras/stae2084

Advance Access publication 2024 September 06

# Assembly of the intracluster light in the HORIZON-AGN simulation

Harley J. Brown<sup>1, ∗</sup> Garreth Martin<sup>1, †</sup> Frazer R. Pearce,<sup>1</sup> Nina A. Hatch<sup>1, ‡</sup> Yannick M. Bahé<sup>1, §</sup> and Yohan Dubois<sup>3</sup>

Accepted 2024 September 3. Received 2024 August 12; in original form 2024 June 20

## ABSTRACT

The diffuse stellar component of galaxy clusters made up of intergalactic stars is termed the intracluster light (ICL). Although there is a developing understanding of the mechanisms by which the ICL is formed, no strong consensus has yet been reached on which objects the stars of the ICL are primarily sourced from. We investigate the assembly of the ICL starting approximately 10 Gyr before $z=0$ in 11 galaxy clusters (halo masses between ~ $1 \times 10^{14}$ and ~ $7 \times 10^{14}$ $M_{\odot}$ at $z \approx 0$) in the HORIZON-AGN simulation. By tracking the stars of galaxies that fall into these clusters past cluster infall, we are able to link almost all of the $z \approx 0$ ICL back to progenitor objects. Satellite stripping, mergers, and pre-processing are all found to make significant contributions to the ICL, but any contribution from *in situ* star formation directly into the ICL appears negligible. Even after compensating for resolution effects, we find that approximately 90 per cent of the stacked ICL of the 11 clusters that is not pre-processed should come from galaxies infalling with stellar masses above $10^9 M_{\odot}$, with roughly half coming from infalling galaxies with stellar masses within half a dex of $10^{11} M_{\odot}$. The fact that the ICL appears largely sourced from such massive objects suggests that the ICL assembly of any individual cluster may be principally stochastic.

**Key words:** methods: numerical – galaxies: clusters: general – galaxies: evolution.

## 1 INTRODUCTION

The diffuse ensemble of stars that permeates the intergalactic space within a galaxy cluster is known as the intracluster light, or ICL. Said to have been first theorized and then discovered by Zwicky (1937, 1951), this low surface brightness cluster component has been traced out to $\gtrsim 1$ Mpc from the cluster centre in stacked observational studies (Zhang et al. 2019), and is expected to contribute significantly to the stellar mass budget of any cluster – between ~ 5 and ~ 50 per cent (varying between clusters and with the exact definition used for the ICL; see Mihos 2015; Contini 2021; or Montes 2022).

It is currently thought that the ICL is formed through three primary channels: from the stellar detritus of massive galaxy mergers (e.g. Murante et al. 2007; Contini, Yi & Kang 2018), from stars stripped from satellite galaxies within the cluster by gravitational interactions (e.g. Willman et al. 2004; Rudick et al. 2009), and also from the accreted intragroup light (IGL) of galaxy groups that previously merged with the cluster (e.g. Rudick, Mihos & McBride 2006; Mihos et al. 2017), with this final formation channel commonly called ICL pre-processing. Although once considered potentially significant formation channels, total disruption of dwarf galaxies (e.g. Purcell, Bullock & Zentner 2007) and *in situ* formation of stars into the ICL (e.g. Puchwein et al. 2010) are now thought to typically only minimally contribute to the ICL (Martel, Barai & Brito 2012;

Melnick et al. 2012; DeMaio et al. 2018; Gullieuszik et al. 2020; but also see Ahvazi et al. 2024a and Bahé et al. [in preparation] as well).

Which or whether any of the three primary formation channels for the ICL is dominant is not yet fully understood, though there is a developing consensus in support of the two-phase formation scenario (Kluge et al. 2020; Golden-Marx et al. 2023) in which stripping in general contributes the most by $z \sim 0$, having potentially superseded the overall contribution from mergers sometime around $z \sim 1$ (see Joo & Jee 2023; Jiménez-Teja et al. 2024, and references therein). Numerous recent studies – both theoretical (e.g. Tang et al. 2023; Contini et al. 2024a) and observational (e.g. Montes & Trujillo 2018; Yoo et al. 2021) – support the idea of stripping being the generally dominant formation channel at $z \approx 0$, especially for lower mass haloes (Contini et al. 2024b), whereas with more massive haloes pre-processing is expected to play an increasingly significant role (Ragusa et al. 2023; Chun et al. 2024; Contini et al. 2024a). For any specific individual cluster, however, the dominant source as well as other properties of the ICL will be intimately linked with the dynamical state of that specific cluster (Chun et al. 2023).

It has been shown that ICL can provide remarkable insight into the accretion history of a cluster, with the ICL serving as a fossil record for cluster dynamics (Montes 2022). Furthermore, as the constituent stars of the ICL are collision-less, with their motions governed primarily by the overall cluster potential rather than that of individual galaxies, it is expected that these stars will behave much like dark matter (DM) and so should be distributed similarly within the cluster. It has thus been proposed that the ICL might be used as a visible DM tracer (e.g. Montes & Trujillo 2019; Alonso Asensio

<footnote n="1"><sup>1</sup>School of Physics & Astronomy, University of Nottingham, University Park, Nottingham NG7 2RD, UK</footnote>

<footnote n="2"><sup>2</sup>Laboratory of Astrophysics, Institute of Physics, Ecole Polytechnique Fédérale de Lausanne (EPFL), Observatoire de Sauverny, CH-1290 Versoix, Switzerland</footnote>

<footnote n="3"><sup>3</sup>Institut d'Astrophysique de Paris, CNRS, 98 bis blvd Arago, F-75014 Paris, France</footnote>

<footnote>* E-mail: Harley.Brown@nottingham.ac.uk</footnote>
et al. 2020; Yoo et al. 2024) and also that it might reveal key features of the cluster potential such as the splash-back radius (e.g. Deason et al. 2020; Gonzalez et al. 2021).

Advances in recent years, both in technique and instrumentation, have pushed observational surface brightness limits to sufficiently low levels to enable revolutionary new ICL studies (see Mihos 2019 for a recent brief review) such as the first ICL study using JWST data by Montes & Trujillo (2022), and that by Kluge et al. (2024) using Euclid Early Release Observations to study the ICL of the Perseus cluster. However, these observational studies require theoretical counterparts so that their results might be interpreted and contextualized, necessitating investigations into the predictions of numerical simulations concerning the ICL.

Although there is a developing consensus on the dominant formation channel of the ICL, a consensus has not yet been reached on the objects from which the stars of the ICL are primarily sourced. Theoretical studies generally agree on more massive galaxies as the more likely source for the majority of ICL stars by $z \sim 0$ (though not unanimously), while observational studies remain more strongly divided.

Numerical simulations of galaxy clusters can be classified into a number of distinct types which adopt different methods (see Somerville & Davé 2015 for a review), including semi-analytical models (SAMs) and numerical hydrodynamic simulations. SAMs couple a simplified, analytical description of baryonic physics post hoc on to an existing N-body DM simulation, while hydrodynamical simulations instead numerically solve equations for gravity, thermodynamics, and hydrodynamics, for DM, gas, and stars simultaneously. Although relatively computationally inexpensive, baryons are not explicitly modelled in SAMs and so the physical processes responsible for the formation of the ICL are included only through approximate, simplified models. Hydrodynamical models (see Vogelsberger et al. 2020 for a review) are considerably more computationally expensive than SAMs, but allow direct and self-consistent modelling of the physical processes that produce the ICL. However, limited computational resources force these simulations to compromise between total simulation volume and the resolution at which the physics can be explicitly simulated. Realistically modelling ICL formation requires simulation volumes large enough to provide cosmological context but this comes at the cost of a decreased resolution, which reduces the accuracy of the simulation at the small scale and may lead to non-physical numerical resolution effects.

Theoretical studies using SAMs have previously found that the main contributors of stars to the ICL by $z=0$ should be galaxies with stellar masses $\gtrsim 10^{10.5} M_\odot$ (Contini et al. 2014; Contini, Yi & Kang 2019) and that the ICL can be heavily influenced by a small number of massive progenitors – between $\sim 1$ and $\sim 10$ per cent the mass of the brightest cluster galaxy (BCG) (Cooper et al. 2015). Harris et al. (2017) re-simulated a Fornax-like cluster (halo mass $\sim 4 \times 10^{13} M_\odot$) from a DM only N-body simulation by replacing haloes infalling after $z=1.65$ with full galaxy models, and found $> 60$ per cent of the $z=0$ ICL of this cluster to have been sourced from just two massive objects (stellar masses $\sim 5 \times 10^{10} M_\odot$). By employing a similar ‘galaxy replacement technique’, Chun et al. (2023) studied six clusters (virial masses of order $10^{14} M_\odot$) and found the ICL in all but the most relaxed cluster to be dominated by stars from galaxies with infall stellar masses $\gtrsim 10^{10} M_\odot$. Continuing this work, Chun et al. (2024) reported the typical progenitors of the IGL/ICL in 84 simulated groups and clusters ($13.6 < \log_{10}(M_{200}[M_\odot]) < 14.8$) to be galaxies with stellar masses between $10^{10}$ and $10^{11} M_\odot$ on cluster infall, though also noted that the ICL of particularly massive

or unrelaxed clusters could have significant contributions from galaxies with infall stellar masses $> 10^{11} M_\odot$ as well. Ahvazi et al. (2024b) conducted a study on 39 groups and clusters (virial masses between $5 \times 10^{12}$ and $2 \times 10^{14} M_\odot$) in the TNG50 cosmological (magneto)hydrodynamical simulation (Nelson et al. 2021) and found that half the IGL/ICL across all the considered systems was brought in by galaxies with stellar masses between $10^{10}$ and $10^{11} M_\odot$, with the ICL of some systems almost entirely originating from objects as or more massive than the Milky Way.

Contrarily, when Tang et al. (2023) analysed mock images of massive clusters at redshifts between $\sim 0.1$ and $\sim 1$ from the TNG100 simulation (Nelson et al. 2021) they found the ICL to lie closest in the age–metallicity plane to satellites with stellar masses between $10^8$ and $10^{10} M_\odot$. They consequently suggested that these intermediate-mass galaxies were the main source of the ICL at these redshifts, as opposed to more massive galaxies. When considering the tension between their results and those of other theoretical studies, however, they also note that the methods they employed to extract galaxy stellar masses in these mock observations (from Tang et al. 2021) could yield lower masses than would be found by traditional substructure extraction algorithms.

Many observational studies support the picture generally presented by simulations – that the stars of the ICL are primarily sourced from galaxies with stellar masses equal to or greater than approximately $10^{10} M_\odot$. For example, Montes et al. (2021) suggested – based on considerations of colour profiles – that the ICL of A85 ($z \approx 0.05$, $M_{200} \approx 1.7 \times 10^{15} M_\odot$) was built up mainly by the stripping of satellites with stellar masses of order $10^{10} M_\odot$. Similarly, Montes & Trujillo (2014, 2018) found (also based on colours) that the ICL of several massive ($M_{200} > 10^{14.3} M_\odot$) clusters at $0.3 < z < 0.6$ appeared to primarily come from Milky Way-like objects, and DeMaio et al. (2018) found (on the basis of colour matching) that 75 per cent of the IGL/ICL in the groups and clusters ($3 \times 10^{13} \lesssim M_{500}/M_\odot \lesssim 9 \times 10^{14}$) they considered at $0.29 < z < 0.89$ should have originated from galaxies with stellar mass $> 10^{10.4} M_\odot$.

Conversely, Morishita et al. (2017) found the ICL of six Hubble Frontier Field clusters ($0.3 < z < 0.6$, $M_{500} \gtrsim 10^{15} M_\odot$) to be dominated by moderately old stars ($\sim 1–3$ Gyr), and found the colours of these stars to be more consistent with these having been stripped mainly from cluster galaxies with stellar masses less than $\sim 10^{9.5} M_\odot$ since $z \sim 1$, instead of more massive galaxies. Likewise, Gu et al. (2020) studied portions of the ICL in the Coma cluster ($z \approx 0.024$, $M_{200} \approx 5.1 \times 10^{14} M_\odot$; Gavazzi et al. 2009), unearthing a very old and metal-poor stellar population – similar to that found by Williams et al. (2007) for intracluster stars in the Virgo cluster ($z \approx 0.004$, $M_{200} \approx 1.4 \times 10^{14} M_\odot$; Urban et al. 2011) – and suggested that these stars may have primarily originated from the accretion of low-mass galaxies (stellar masses less than $\sim 3 \times 10^9 M_\odot$). However, the authors of both these studies note that their analyses consider only typical galaxy colours and metallicities, and invoke the presence of strong colour gradients and lower outskirt metalliticites in more massive galaxies – and that the stars in these galaxies at large galactocentric radii should be more easily stripped and so make up a disproportionate fraction of the stars these galaxies contribute to the ICL – as a potential explanation for the tension between their results and the emerging consensus from theoretical studies.

The tension that remains between predictions for the main progenitors of the ICL – both between theoretical and observational studies, and between different studies within each of these categories – merits further investigation. A small number of prior studies have already investigated the primary progenitors of ICL stars as predicted by hydrodynamical simulations (e.g. Ahvazi et al. 2024b). However,
their conclusions are likely, at least in part, dependent on the specific
physical models used by each different simulation code and there is
therefore worth in investigating this area with additional simulations.

In this paper, we present a new study of the origin and formation
of the ICL using the HORIZON-AGN simulation (Dubois et al. 2014).
We investigate the assembly of the ICL starting approximately
10 Gyr before z = 0 in 11 clusters with DM halo masses between
~1 × 10¹⁴ and ~7 × 10¹⁴ M☉ at z ≈ 0. By tracking stars associ-
ated with galaxies that fall into these clusters during this time, we
are able to link much of the z ≈ 0 ICL with progenitor galaxies
and so investigate the relative contributions to the ICL of galaxies
with differing stellar masses preceding cluster infall. Additionally,
by extrapolating our results for intermediate and high stellar mass
progenitor objects (≈ 10⁹ M☉) down to lower masses to compensate
for resolution effects, we estimate the range of galaxy stellar masses
responsible for the bulk of the ICL in a typical cluster.

This paper is organized as follows: In Section 2, we present the
details of the HORIZON-AGN simulation (Section 2.1) as well as
those of the ADAPTAHOP structure finder and TREEMAKER merger
tree builder, and then describe how we employ these for our
investigation into the assembly of the ICL (Section 2.2). In Sec-
tion 3, we describe and discuss the results of our investigation
– first on the stellar mass budgets of the 11 clusters and what
fraction of ICL stars could formerly be found in galaxies within
the clusters (Section 3.1), then on the contributions of progenitor
objects with differing stellar masses on cluster infall to the ICL
(Section 3.2). We conclude by summarizing the main results of
our investigation in Section 4. Throughout this paper, we assume a
$\Lambda$CDM cosmology with $h = H_0/100\text{kms}^{-1}\text{Mpc}^{-1} = 0.704$. Unless
stated otherwise, distances are given in proper rather than co-moving
units.

2 METHODS

## 2.1 HORIZON-AGN

The full details of the HORIZON-AGN simulation can be found in
Dubois et al. (2014, 2016), and are reproduced only in brief here.

The HORIZON-AGN simulation is a cosmological-volume hydro-
dynamical simulation which employs 1024³ DM particles (each with
mass 8 × 10⁷ M☉) in a cubic volume of side-length 100 h⁻¹ Mpc
(co-moving) with periodic boundary conditions. A standard $\Lambda$CDM
cosmological model (Peebles & Ratra 2003) is used with total matter
density $\Omega_m = 0.272$, baryon density $\Omega_b = 0.045$, dark energy den-
sity $\Omega_\Lambda = 0.728$, DM power spectrum amplitude $\sigma_8 = 0.81$, Hubble
constant $H_0 = 70.4 \text{ km s}^{-1} \text{Mpc}^{-1}$, and spectral index $n_s = 0.967$,
compatible with the *Seven-Year Wilkinson Microwave Anisotropy
Probe* cosmology (Komatsu et al. 2011), and also that of the Planck
Collaboration XVI (2014) within ~10 per cent relative variation.
The simulation utilizes RAMSES (Teyssier 2002) – an adaptive mesh
refinement Eulerian hydrodynamics code – with an initially uniform
1024³ cell gas grid, 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 $\Delta x = 1$ kpc after seven levels of refinement.

Gas heating from a uniform ultraviolet background begins after
z = 10, following Haardt & Madau (1996), and gas is allowed to
cool to 10⁴ K via H and He collisions with a contribution from
metals using a Sutherland & Dopita (1993) model. Feedback from
active galactic nuclei is a key process for regulating the stellar
mass content of massive galaxies (Dubois et al. 2016), and active
galactic nuclei are modelled in HORIZON-AGN by the accretion

of gas on to supermassive black holes following a Bondi–Hoyle–
Lyttleton accretion rate (Bondi 1952) capped at the Eddington rate,
with switching between jet (‘radio’) and heating (‘quasar’) feedback
modes according to accretion rate (Dubois et al. 2012). Star formation
is modelled using a Schmidt law, with 2 per cent star formation
efficiency (Kennicutt 1998), and is only allowed in regions with gas
number density exceeding 0.1 H cm⁻³. Feedback is included from
Type Ia and Type II supernovae, as well as stellar winds, through
mass, energy, and metal release. The star particles (i.e. the resolution
elements of the ICL) have an initial mass resolution of approximately
2 × 10⁶ M☉.

The HORIZON-AGN simulation was not calibrated to the local
Universe, except for choosing black hole feedback parameters.
Despite this, Kaviraj et al. (2017) reported that the simulation
is in generally good agreement with observations (using observa-
tional data from 0 < z < 6) concerning e.g. predicted evolution
of luminosity functions, stellar mass functions, the star formation
main sequence, and the cosmic star formation history. A slight
overabundance of galaxies with stellar masses ≲ 10¹⁰.⁵ M☉ at all
epochs is noted compared to observations, though agreement can be
achieved by considering observational uncertainties due to cosmic
variance.

## 2.2 Tracking galaxies

2.2.1 ADAPTAHOP and TREEMAKER

For identifying structures of both DM particles (i.e. haloes) and
star particles (i.e. galaxies), the ADAPTAHOP structure finder (Aubert,
Pichon & Colombi 2004) is employed (and specifically the version
updated by Tweed et al. 2009 for building merger trees), the details
of which are reproduced only in brief here.

Particles are grouped into structures on the basis of local particle
density, calculated using a nearest neighbours approach with 20
neighbour particles. No unbinding procedure is employed when
creating these particle groups. Stellar structures and substructures
are assembled from the groups created by linking star particles with
densities > 178× the total matter density (of the entire simulated
volume) according to their closest local density maxima. Saddle
points in the density field between these groups are then used to link
these groups together into stellar structures (i.e. galaxies), which are
then subdivided hierarchically on the basis of increasing density to
identify substructure. The distance used for force softening is ~ 2 kpc
(hence substructures smaller than this are considered irrelevant),
and stellar structures with ≤ 50 star particles are neglected. As the
star particles have a mass of ~ 2 × 10⁶ M☉, this gives a minimum
(detectable) galaxy stellar mass of ~ 10⁸ M☉. We address the
potential contribution to the ICL from galaxies less massive than
this detection threshold in Section 3.1.

Essentially the same procedure is used to identify DM structures
(and substructures), though the density threshold used is instead 80×
the total matter density, and the minimum membership threshold for
a structure to not be neglected is raised to 100 particles. Galaxies
and DM haloes are identified independently and only afterwards are
galaxies linked to host haloes. The main galaxy of a halo is defined
to be the most massive galaxy within 0.1 × r₁₇₈ of the halo centre
(where r₁₇₈ is the radius from the halo centre within which the average
DM density is 178× the critical density).

The TREEMAKER algorithm (originally developed by Hatton et al.
2003; see also Tweed et al. 2009) is used to generate merger trees
for the galaxies identified by ADAPTAHOP. In brief, each structure in
each simulation snapshot is (if possible) linked to a main descendent
structure in the next snapshot and a main progenitor in the previous snapshot using a merit function which maximizes the fraction of particles shared by the would-be linked structures. Between $z \approx 3$ and $z \approx 0$ the mergers trees we use in this study were built using a time resolution of $\sim 0.05$ Gyr.

For this study, the galaxy considered to be the BCG of each cluster is selected at $z \approx 0$ ($z = 0.0556$) as the most massive galaxy within $0.1 \times r_{178}$ of the cluster halo centre. Prior to $z \approx 0$, the galaxy classified as the BCG of the cluster in each snapshot is the main progenitor of the $z \approx 0$ BCG. We place the border of each cluster in each snapshot at a distance away from the BCG equal to $r_{178}$ of the cluster halo (see Section 2.2.2) and classify all galaxies within this border besides the BCG as satellite galaxies.

### 2.2.2 Tagging and tracking the stars of infalling galaxies

For this investigation, we select 11 clusters from the HORIZON-AGN simulation, with $z \approx 0$ ($z = 0.0556$) DM halo masses (within $r_{178}$ i.e. $M_{178}$) in the range $\sim 1 \times 10^{14}$ to $\sim 7 \times 10^{14}$ M$_\odot$ and investigate their ICL assembly. We restrict our analysis to cubes, side-length 8 Mpc, centred on the BCG of each cluster. For each cluster, we perform our analysis using a coarse time resolution of $\sim 1$ Gyr, down to $z \approx 0$ ($z = 0.0556$) from as far back as when a main progenitor of the cluster BCG is first identified (which is always before $z \sim 1.5$ and almost always before $z \sim 2$) or $z \sim 3$ – whichever is later. As later described in Section 3.1, we never find more than $\sim 0.1$ per cent of the $z \approx 0$ ICL stars of any of the 11 clusters to have assembled before we begin to monitor ICL assembly. This is in good agreement with previous theoretical studies (e.g. Willman et al. 2004; Contini et al. 2024a) which typically found the vast majority of the $z = 0$ ICL mass of similarly massive clusters to have still been contained in galaxies at $z \sim 1$.

We place the border of each cluster at $r_{178}$ and consider galaxies to have become satellites of a cluster once they pass $r_{178}$. Rather than use the values provided by ADAPTAHOP directly, for each cluster we instead fit a monotonic cubic spline to the $r_{178}$ values of each cluster from ADAPTAHOP as a function of look-back time (using snapshots with $\sim 1$ Gyr coarse spacing), allowing us to account for and smooth out any temporary and non-physical fluctuations in cluster size (such as those caused by cluster mergers).

To identify the progenitor objects of the $z \approx 0$ ICL, we track the star particles of galaxies falling into the 11 clusters starting from the coarse snapshot immediately proceeding galaxy infall down to $z \approx 0$. Infalling galaxies are identified as galaxies at $r > r_{178}$ in any coarse snapshot whose main descendant in the next coarse snapshot is at $r < r_{178}$. This infaller galaxy identification is done chronologically, so galaxies which have previously fallen into the cluster but whose orbit apocentre still lies outside $r_{178}$ can be identified (and repeat-counting avoided) by confirming that each candidate infaller galaxy does not have a previously identified infaller as a main progenitor. When a galaxy infalling for the first time is identified, all of its stars in the coarse snapshot immediately proceeding infall are tagged as associated with that progenitor, and the galaxy then followed through all remaining coarse snapshots using the merger trees.

Any new star particles born in the descendent of an infaller are also tagged as associated with that infaller, provided that the infaller galaxy is still part of the main progenitor branch of the descendant galaxy those star particles are born in. This restriction is imposed to avoid double-counting (so these new star particles are only linked to a single progenitor rather than potentially several progenitors that merged after infall). We find these ‘newborn’ star particles – those

born during or after galaxy infall – to make up $\sim 4$ per cent of all $z \approx 0$ ICL stars across all 11 clusters (and $\sim 6$ per cent of all not pre-processed ICL stars), so include them in our analysis for completeness, attributing them as contributed to the ICL by the progenitor object they are tagged to and equivalent to the stars of the same galaxy that formed long before cluster infall. However, by repeating our analysis while ignoring these newborn stars we confirm that whether or not they are included does not meaningfully alter our conclusions.

We use the merger trees generated by TREEMAKER to follow galaxies past cluster infall, but the extreme environment of a cluster makes errors in linking progenitor and descendent objects somewhat common, leading to some missing or incorrect links. Whenever we identify a break in the merger tree of a tracked galaxy or find that no star particles are shared between a tracked galaxy and its supposed descendent in the next coarse snapshot, and another galaxy can also be found in the next coarse snapshot which contains $> 50$ per cent of the tracked galaxy’s star particles, we overrule TREEMAKER and classify the latter galaxy as the tracked galaxy’s proper descendent. Of the $\sim 4000$ infalling galaxies tracked, we modify the merger trees of $\sim 15$ per cent after initial infall – most commonly to patch breaks while galaxies are partway through merging with the BCG.

### 2.2.3 Identifying the ICL

We identify the ICL using the output from ADAPTAHOP. Specifically, any star particles in the clusters not classed by ADAPTAHOP as part of a stellar structure are considered to be part of the ICL (and so the definition for the ICL employed in this study is one based on instantaneous three-dimensional density). The ICL fractions (i.e. fraction of cluster stellar mass attributed to the ICL rather than cluster galaxies) according to this definition of the ICL for the 11 clusters at $z \approx 0$ are shown in Table 1, along with a selection of other cluster properties.

It bears acknowledging that different ICL identification methodologies (e.g. different structure finders) will disagree on the exact stellar content of the ICL, so our choice of structure finder will have some impact on our results. However, prior studies have shown broadly good agreement across a wide range of structure finder codes (for halo properties such as virial mass, bulk velocity, rotation velocity, and the presence and location of substructure: see Knebe et al. 2011, Onions et al. 2012 and references therein), and Brough et al. (2024) demonstrated that ADAPTAHOP yields ICL properties in good agreement with those recovered using other structure finder codes, such as SUBFIND (Dolag et al. 2009) which includes an unbinding procedure. As such, we do not anticipate our choice of structure finder will significantly influence our results. We address the related issue of splitting the BCG and the ICL in Section 3.1.

An example of one of the 11 clusters (ID 19 in Table 1) is shown in Fig. 1; the left and centre panels show the projected positions for an arbitrary line of sight at $z \approx 0$ of all stars and only ICL stars, respectively. The stars shown in the left panel are coloured according to whether they were classified by ADAPTAHOP as being part of the BCG (red), a satellite galaxy (green), or the ICL (blue), while those shown in the central panel are coloured according to associated progenitor. ICL star particles unable to be associated with a progenitor object (such as those which did not enter the cluster as part of a galaxy, i.e. pre-processed ICL) are shown in black; otherwise the specific colours used do not relate to any physical properties of the stars or their progenitor. The right panel in Fig. 1 is
<table><caption>**Table 1.** Various properties of the 11 simulated clusters. The assembly redshifts are those when a main progenitor of the cluster halo first has a total DM mass equal to or exceeding 50 per cent of that of the $z \approx 0$ cluster. The DM and stellar masses are those within $r_{178}$ at $z \approx 0$ and include substructure, and $f_{BCG}$, $f_{Sat}$, and $f_{ICL}$ are the fractions of this stellar mass in the BCG, satellite galaxies, and ICL, respectively, and are also depicted in Fig. 2. The ICL sub-fractions given in the last three columns are also depicted in Fig. 3.</caption><thead><tr><th>Cluster ID</th><th>DM halo mass (10<sup>14</sup> M<sub>☉</sub>)</th><th>Stellar mass (10<sup>12</sup> M<sub>☉</sub>)</th><th>Assembly redshift</th><th>f<sub>BCG</sub></th><th>f<sub>Sat</sub></th><th>f<sub>ICL</sub></th><th>ICL fraction directly from satellites</th><th>ICL fraction previously in BCG</th><th>ICL fraction from pre-processing</th></tr></thead><tbody><tr><td>1</td><td>1.44</td><td>4.74</td><td>0.415</td><td>0.304</td><td>0.566</td><td>0.130</td><td>0.523</td><td>0.339</td><td>0.138</td></tr><tr><td>9</td><td>0.973</td><td>2.84</td><td>1.24</td><td>0.506</td><td>0.352</td><td>0.142</td><td>0.442</td><td>0.476</td><td>0.0818</td></tr><tr><td>13</td><td>6.86</td><td>19.6</td><td>0.122</td><td>0.132</td><td>0.740</td><td>0.129</td><td>0.446</td><td>0.115</td><td>0.438</td></tr><tr><td>19</td><td>5.32</td><td>14.4</td><td>0.594</td><td>0.338</td><td>0.495</td><td>0.168</td><td>0.439</td><td>0.394</td><td>0.167</td></tr><tr><td>46</td><td>1.17</td><td>3.54</td><td>0.566</td><td>0.336</td><td>0.549</td><td>0.115</td><td>0.511</td><td>0.368</td><td>0.121</td></tr><tr><td>48</td><td>1.03</td><td>3.28</td><td>0.745</td><td>0.353</td><td>0.518</td><td>0.129</td><td>0.435</td><td>0.423</td><td>0.142</td></tr><tr><td>49</td><td>3.09</td><td>8.66</td><td>1.07</td><td>0.346</td><td>0.514</td><td>0.139</td><td>0.522</td><td>0.359</td><td>0.119</td></tr><tr><td>71</td><td>1.34</td><td>4.78</td><td>0.369</td><td>0.316</td><td>0.572</td><td>0.112</td><td>0.437</td><td>0.298</td><td>0.264</td></tr><tr><td>132</td><td>1.20</td><td>4.08</td><td>0.270</td><td>0.310</td><td>0.588</td><td>0.101</td><td>0.409</td><td>0.283</td><td>0.308</td></tr><tr><td>174</td><td>2.09</td><td>6.17</td><td>0.233</td><td>0.248</td><td>0.623</td><td>0.128</td><td>0.584</td><td>0.190</td><td>0.226</td></tr><tr><td>183</td><td>1.08</td><td>3.64</td><td>0.332</td><td>0.348</td><td>0.540</td><td>0.112</td><td>0.428</td><td>0.399</td><td>0.173</td></tr></tbody></table>

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Left: Projected positions of all the stars in the vicinity of an example $z \approx 0$ HORIZON-AGN galaxy cluster (ID 19 in Table 1) using a projection depth of 8 Mpc (centred on the BCG) and coloured according to whether a star is considered part of the BCG (red), a satellite galaxy (green), part of the ICL (blue), or is outside of the cluster (grey). Middle: Projected positions of only intergalactic stars for the same cluster as shown on the left. ICL stars associated with the same progenitor galaxy are shown in the same colour; those not associated with a progenitor infaller (e.g. pre-processed ICL) are shown in black. Intergalactic stars outside of the cluster (e.g. the ICL of infalling groups) are shown in grey. Right: Scatter plot of the speeds, $v$, of intergalactic stars (relative to the BCG and scaled by the 3D velocity dispersion of the cluster stars, $\sigma_v$) against distance from BCG centre, $r$, (scaled by $r_{178}$) for the same cluster and using the same colours as the middle panel.</figcaption></figure>

a velocity dispersion plot (speed relative to the BCG plotted against distance from BCG centre, normalized by the 3D velocity dispersion of all cluster star particles and $r_{178}$, respectively) for the same ICL star particles as shown in the centre panel of Fig. 1 and employing the same scheme for colours. Comparison between the centre and right panels shows correspondence between features in position and velocity space, as expected for ICL stars recently sourced from the same infalling galaxy.

### 2.2.4 Uncertainty estimation

Throughout the remainder of this paper, all indicated uncertainties are estimated using bootstrapping. Specifically, we produce 5000 bootstrap resamples of our population of ~4000 tracked infaller galaxies and their relative contributions to the ICL (sampled with replacement to produce infaller galaxy resamples of equal size to the original population), and repeat our analysis on each of these resamples. The indicated uncertainties are found using the 16th and 84th percentiles from these repeat analyses.

## 3 RESULTS

### 3.1 Stellar mass budget and ICL pre-processing

Fig. 2 depicts how the $z \approx 0$ stellar mass budgets of each of the 11 clusters are split between satellite galaxies, the BCG, and the ICL. The combined BCG + ICL fractions are also included. The median $z \approx 0$ ICL (BCG + ICL) stellar mass fraction of the clusters is 13 (45) per cent, the 16th percentile is 11 (40) per cent, and the 84th percentile is 14 (49) per cent.

These ICL fractions fall within the wide range of values reported in the literature – between ~5 and ~50 per cent (Contini 2021; see also Tang et al. 2018; Kluge et al. 2021) – and agree fairly well with the findings of other theoretical studies employing similar ICL definitions, such as that by Rudick, Mihos & McBride (2011) who found ICL fractions between ~9 and ~15 per cent while employing an ICL definition based on instantaneous density in their investigation of a selection of simulated clusters with masses of order ~$10^{14}$ M$_{☉}$. These ICL fractions also broadly agree with many of those found for similarly massive clusters by studies which employed different ICL
<figure><img src="image_3.png" /><figcaption>**Figure 2.** The $z \approx 0$ stellar mass budgets of the 11 simulated clusters, divided between satellite galaxies, the BCG, and the ICL. Combined BCG + ICL stellar mass fractions also shown. Dotted lines indicate 16th and 84th percentiles, and solid lines the maximum, minimum, and median values.</figcaption></figure>

definitions – both theoretical (e.g. ~9–14 per cent; Ahvazi et al. 2024b) and observational (e.g. ~7–15 per cent; Mihos et al. 2017) – though some studies instead report much higher ICL fractions (e.g. $\ge$ 30 per cent, Pillepich et al. 2018; 20–40 per cent, Furnell et al. 2021).

A major driver of the large variation in ICL fractions found by both theoretical and observational ICL studies is the wide variety of different methods employed to measure the ICL (Brough et al. 2024). These different ICL measures differ most significantly in how they separate the ICL from the BCG, with large variations seen in the ICL fraction depending on where the border between the ICL and BCG is drawn. This perspective is supported by the broad agreement seen in BCG + ICL stellar fractions between studies – generally around ~50 per cent for clusters with masses of order $\sim 10^{14}$ M$_\odot$ (e.g. Pillepich et al. 2018; Brough et al. 2024; Contreras-Santos et al. 2024, in addition to this study) – even when measured ICL fractions differ significantly. Consequently, rather than implement an uncertain border between them, some previous studies have opted to simply not separate the BCG from the ICL and consider the combined system only (e.g. Pillepich et al. 2018; Zhang et al. 2019; Chun et al. 2023). In similar acknowledgement of this definition issue, we repeat the remainder of our analysis while likewise not separating the BCG from the ICL, and where relevant present this alternative analysis alongside our primary analysis for comparison (and always refer to this combined system as BCG + ICL). This definition issue also has significant implications for any study of the radial dependence of ICL sourced from different formation channels, so we defer investigating such a dependence to future work.

Fig. 3 depicts the fractions of $z \approx 0$ ICL stars in the 11 clusters that were previously part of a satellite galaxy (but never part of the BCG), previously part of the BCG, or that have never been part of a cluster galaxy (i.e. were pre-processed). The combined fractions

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Fractions of $z \approx 0$ ICL stars formerly part of satellite galaxies, previously part of the BCG, or that were never part of a cluster galaxy or an infalling galaxy (i.e. pre-processed ICL) in the 11 clusters. Combined fractions from either the BCG or directly from satellite galaxies (i.e. not pre-processed ICL) also shown. Dotted lines indicate 16th and 84th percentiles, and solid lines the maximum, minimum, and median values.</figcaption></figure>

of ICL stars that were previously part of any cluster galaxy are also included.

The median fraction of $z \approx 0$ ICL stars formerly in satellite galaxies (and never part of the BCG) is 44 per cent [16th percentile (43 per cent) and 84th percentile (52 per cent)], and the median fraction that were formerly part of the BCG is 36 per cent [16th percentile (26 per cent) and 84th percentile (41 per cent)]. That never less than 40 per cent of the $z \approx 0$ ICL stars in any of the 11 clusters were liberated directly from satellite galaxies within the cluster (and had never been part of the BCG) indicates stripping must play a significant role in ICL assembly. Likewise, that on average more than a third of $z \approx 0$ ICL stars were formerly part of the BCG of their cluster implies violent BCG mergers can play a similarly significant role. We revisit the importance of BCG mergers for ICL assembly in Section 3.2.1.

When Montenegro-Taborda et al. (2023) investigated the ICL assembly of $M_{200} \ge 10^{14}$ M$_\odot$ clusters in the TNG300 simulation (Nelson et al. 2021), they found the mean fraction of ICL stars (at $z=0$) that had been stripped from surviving satellites to be 38.1 per cent. The value we find for the median fraction of ICL stars directly stripped from satellite galaxies (~44 per cent) is somewhat larger than this, as would be expected since (unlike Montenegro-Taborda et al. 2023) we do not exclude ICL stars stripped from satellite galaxies that later merge with the BCG or otherwise do not survive until $z=0$. Chun et al. (2024) reported a median value of ~50 per cent for the fraction of the $z=0$ ICL stars of $14 < \log(M_{200}/M_\odot) < 14.4$ clusters that could be linked to the merger tree of the cluster BCG. As this value combines together the ICL contributions from BCG mergers and from the tidal stripping of galaxies prior to merging with the BCG, it is unsurprising that this value is larger than the median fraction of the ICL we find to have previously been part of the BCG (~36 per cent). We thus do not
consider our results in tension with those of Montenegro-Tabora et al. (2023) or Chun et al. (2024).

Using SAMs, Contini et al. (2024a) obtained a median $z = 0$ ICL contribution from mergers of only ~10 per cent for clusters with halo masses between ~$10^{14}$ and ~$10^{14.5}$ M$_{\odot}$. This low fraction does not appear compatible with the median fraction of ICL stars we find to have formerly been part of the BCG (~36 per cent), as mergers should be the primary mechanism for liberating these stars. We speculate that this disagreement may stem from the simple prescription employed for mergers by the utilized SAM, which does not allow BCG stars to be added to the ICL (instead just moving 20 per cent of the stellar mass of a merging satellite into the ICL component when a BCG 'merger' occurs). The aforementioned issue of studies which employ differing methodologies for separating the BCG and the ICL determining significantly differing properties for the ICL may also play some part in this apparent incompatibility.

We note that not all of the $z \approx 0$ ICL stars identified by ADAPTAHOP as previously being part of a cluster galaxy actually entered that cluster as part of a galaxy. For example, ~20 per cent of the $z \approx 0$ ICL stars across all 11 clusters classed by ADAPTAHOP as formerly being in satellite galaxies could not be linked to a progenitor infaller galaxy (as they neither entered the cluster as part of that infaller nor were they formed in its descendent). These stars would often appear distinct from the rest of their supposed galaxy in phase space, and were frequently only classified as part of a cluster galaxy for a single coarse snapshot, having originally entered the cluster while not classed by ADAPTAHOP as part of a galaxy (i.e. as pre-processed ICL). We suspect the transient galaxy membership of these stars to be a result of ADAPTAHOP not employing any kind of unbinding procedure. The median fraction (across the 11 clusters) of ICL stars classed as formerly being part of satellite galaxies, the BCG, or either that could not be linked to a progenitor infaller galaxy were 19, 27, and 23 per cent, respectively. As these stars could not be associated with a progenitor infaller, they are not considered (along with all other pre-processed ICL) in our later analysis of the differing contributions of infalling galaxies of varying masses to the ICL.

As a further consequence of ADAPTAHOP not employing an unbinding procedure, we note that any ICL stars on radial orbits that are coincidentally caught within the BCG at the time of a coarse snapshot would be classed by ADAPTAHOP as part of the BCG, and hence be added to the category referred to as 'previously in BCG' in Fig. 3. As a result, we caution that this category is not an ideal proxy for the ICL contribution from violent mergers with the BCG and should instead be interpreted as an upper limit for the contribution from this channel.

Among the 11 clusters, the median fraction of $z \approx 0$ ICL stars that have never been classed as part of a cluster galaxy is 16 per cent [16th percentile (12 per cent) and 84th percentile (28 per cent)]. This is in middling agreement with the pre-processed ICL fractions found (using SAMs) by Contini et al. (2024a), who found a mean pre-processed fraction of ~25 per cent (with significant scatter) for clusters with halo masses between ~$10^{14}$ and ~$10^{14.5}$ M$_{\odot}$.

The pre-processed fractions shown in Fig. 3 feature one significant outlier: ID 13, the only cluster with a pre-processed fraction $\gtrsim 30$ per cent. As shown in Table 1, this is also the most massive cluster (both in DM and in stellar mass), that with the greatest fraction of its stellar mass contained in satellites, as well as that with the lowest assembly redshift (at $z \approx 0$ this cluster is uniquely part-way through two simultaneous major cluster mergers). It is also of note that the cluster with the lowest $z \approx 0$ pre-processed fraction – ID 9 – is the least massive (both in DM and in stellar mass), that with the lowest fraction of its stellar mass contained in satellites, as well as that with

the highest assembly redshift. These findings appear compatible with those of Chun et al. (2023), who previously noted a trend towards higher BCG + ICL stellar mass fractions in more relaxed clusters.

In addition to pre-processing, we note that the category labelled as such in Fig. 3 also includes minor contributions from sub-threshold galaxies not detectable by the structure finder, ICL already present at high redshifts, and possibly from star formation directly into the ICL (Puchwein et al. 2010; Ahvazi et al. 2024a; and Bahé et al., in preparation). We investigated these minor channels and find that at most ~0.1 per cent of the stacked $z \approx 0$ ICL of the 11 clusters could be from *in situ* star formation,<fnref n="1" /> and never find more than ~0.1 per cent of the $z \approx 0$ ICL of any cluster to have assembled before we begin monitoring (typically ~0.01 per cent). Based on the analysis described in Section 3.2.2, we estimate less than 5 per cent of the stacked $z \approx 0$ ICL of all 11 clusters to be from unresolved galaxies with stellar masses < $10^9$ M$_{\odot}$.

## 3.2 Contribution of infalling galaxies to the ICL

### 3.2.1 Fraction of stars liberated from infalling galaxies

We refer to the fraction of star particles tagged to the same progenitor infaller galaxy that become part of the $z \approx 0$ ICL – regardless of the mechanism that adds those stars to the ICL and considering both stars that were part of the galaxy prior to infall as well as those that later formed in a descendant of the infaller (see Section 2.2.2) – as the liberated fraction, $f_{\text{lib}}$, of that progenitor. Fig. 4 shows the mean $f_{\text{lib}}$ values of tracked infaller galaxies from all 11 clusters in rolling infall stellar mass bins (0.5 dex width) together with their best-fitting cubic spline fit. We use a total of 500 bins between $\log(M_*/M_{\odot}) \approx 8.2$ and $\log(M_*/M_{\odot}) \approx 12.1$ (the range of infall stellar masses seen in the 11 clusters) though bins containing fewer than 30 galaxies from the original sample are always ignored for fitting. Outside the infalling galaxy stellar mass range used for fitting, the fit line is linearly extrapolated. The shaded regions indicate the dispersion of fit lines generated in the same fashion for each of the bootstrap resamples. For clarity, only every 15th bin is shown in Fig. 4. The equivalent analysis for the combined BCG + ICL system is also shown in faint.

As we expect any features in the $f_{\text{lib}}$ fractions on scales smaller than a few tenths of a dex in mass to be the result of noise rather than physically meaningful features, we apply a Savitzky & Golay (1964) filter to the mean $f_{\text{lib}}$ values using a 0.1 dex wide stellar mass window before fitting. Additionally, for this analysis we ignore progenitor galaxies which only ‘skimmed’ a cluster – entering and then exiting $r_{178}$ in consecutive coarse snapshots and never returning before $z \approx 0$ – as well as those that only cross $r_{178}$ for the first time at $z \approx 0$. We exclude these galaxies as we consider them distinct from the population which fall into the clusters $\gtrsim 1$ Gyr prior to $z \approx 0$ and are allowed ample opportunity to be processed by the cluster. By repeating our analysis with these galaxies included, we verify the

<footnote n="1">¹We estimate the *in situ* fraction as the fraction of $z \approx 0$ ICL stars that were first seen with an age less than the elapsed time since the previous coarse snapshot, have never (in any coarse snapshot) been seen outside their $z \approx 0$ cluster (i.e. outside $r_{178}$), and have never (in any coarse snapshot) been classed as part of a stellar structure by ADAPTAHOP. Any to-be $z \approx 0$ ICL stars potentially formed outside of a galaxy but beyond $r_{178}$ – such as in a group that was later accreted by the cluster – are considered pre-processed ICL instead. We regard this estimate as an upper limit as it may include contamination from other minor channels, such as stars that are both born and liberated into the ICL between coarse snapshots.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Mean fraction of liberated stars, $f_{\text{lib}}$, as a function of progenitor galaxy infall stellar mass, $M_*$, for galaxies that fell into any of the 11 clusters. Calculated using 0.5 dex wide rolling bins. A fit to the bins containing more than 30 galaxies with an added Savitzky & Golay (1964) filter is shown. The error bars and shaded regions indicate estimated uncertainties based on bootstrapping. The faint line and symbols show the equivalent fractions for the combined BCG + ICL system.</figcaption></figure>

best-fitting function remains nearly unaffected, bar a small (order 0.01) approximately uniform shift to lower values of $f_{\text{lib}}$.

From Fig. 4 it can be seen that a less massive infalling galaxy is expected to contribute a higher fraction of its stars to the ICL. Despite the most massive galaxies being those that contribute the smallest fraction of their infall stellar mass to just the ICL alone, these galaxies are expected to make major contributions to the combined BCG + ICL system – both in terms of absolute mass and relative to their infall stellar mass.

As they experience stronger dynamical friction, more massive infalling galaxies should fall into the centre of a cluster faster (Contini 2021). Less massive galaxies, on the other hand, will have shallower potential wells and so should be more easily stripped of their stars by gravitational interactions within the cluster (Read et al. 2006). These two phenomena can explain the main trends seen in Fig. 4: very massive infalling galaxies contribute only a small fraction of their stars to the ICL as these galaxies are less easily stripped of their stars during infall, and are expected to spiral into the cluster centre (to merge with the BCG) relatively quickly. Conversely, the least massive infalling galaxies should be the most easily stripped and are expected to typically travel into the cluster centre comparatively slowly, providing ample opportunity for stripping to siphon a large fraction of their stellar mass into the ICL. This perspective is supported by a notably higher fraction of the $z \approx 0$ ICL stars contributed by more massive progenitor galaxies having previously (between progenitor infall and $z \approx 0$) been part of the BCG of their cluster (not shown in Fig. 4): ~25 per cent of all $z \approx 0$ ICL stars (from all 11 clusters) linked to progenitors with infall stellar masses between $10^{8.5}$ and $10^{9.5}$ M$_\odot$ were previously part of the BCG of their cluster, rising to ~37 per cent for stars from progenitors with infall stellar masses between $10^{10.5}$ and $10^{11.5}$ M$_\odot$.

The dominant role of violent mergers as the mechanism that liberates stars from more massive galaxies into the ICL (rather than

gradual stripping) is additionally supported by massive infalling galaxies that merge with the BCG before $z \approx 0$ typically having significantly higher $f_{\text{lib}}$ values compared with those massive infallers that survive as satellites (not shown in Fig. 4). The median (mean) value for $f_{\text{lib}}$ (for the ICL alone) for all $10.5 < \log_{10}(M_*/M_\odot) < 11.5$ infallers from any of the 11 clusters that did not merge with their BCG before $z \approx 0$ is ~0.04 (~0.06), as opposed to ~0.10 (~0.13) for progenitors that did merge with the BCG. More massive infallers quickly merging with the BCG and so violent mergers becoming an increasingly important channel for adding their stars to the ICL may also explain the flattening that can be seen at the high-mass end of the fit for $f_{\text{lib}}$ for the ICL alone in Fig. 4.

We caution that the shape of the fit for $f_{\text{lib}}$ against infall stellar mass seen in Fig. 4 for the ICL alone is suspected to be highly sensitive to the specific ICL definition employed. As described in Section 3.1, different ICL definitions differ most significantly on the border between the BCG and ICL. If an alternative ICL definition were employed that moved the border between the two towards smaller cluster-centric radii, we expect the $f_{\text{lib}}$ fit for the ICL system alone would move upwards to more closely resemble that seen in Fig. 4 for the combined BCG + ICL system. Although exploring this effect of ICL definition further is beyond the scope of this study, we intend to investigate this in a future work using a broader range of simulations.

### 3.2.2 Infalling galaxy mass function

Along with determining the typical contribution to the ICL made by a single galaxy infalling with a particular stellar mass, quantifying the overall expected contribution to the ICL from the entire population of infalling galaxies of a given stellar mass also requires assessing how often galaxies that massive join clusters. A histogram of the stellar masses on infall, $M_*$, of galaxies which fell into any of the 11 clusters after we began to monitor ICL assembly is shown in Fig. 5. A fit to these data is also shown, which adopts the form of a Schechter (1976) function,

$$ \Phi(M_*) \cdot dM_* = \Phi_n \left( \frac{M_*}{M_k} \right)^\alpha \exp(-M_*/M_k) \cdot dM_* \quad (1) $$

where $\alpha$ is the low-mass-end slope of the function, $\Phi_n$ the normalization, and $M_k$ corresponds to the ‘knee’ of the function (i.e. the infall stellar mass when the function exhibits a rapid change in slope). Due to the star particle mass in HORIZON-AGN and the minimum membership threshold of the ADAPTAHOP structure finder, the minimum stellar mass of a detectable galaxy is ~$10^8$ M$_\odot$. The resulting influence of resolution effects can be noted in Fig. 5 for infall stellar masses ~$10^9$ M$_\odot$. To mitigate the impact of these resolution effects on the resulting fitted function, we therefore only use data from $M_* > 10^9$ M$_\odot$ when fitting the Schechter function.

For consistency with the analysis described in Section 3.2.1, galaxies just falling into a cluster at $z \approx 0$ or those that only ‘skimmed’ a cluster rather than infalling are excluded from the histogram shown in Fig. 5. We have verified that the best-fitting parameters are nearly unchanged when these galaxies are included, so that our conclusions are unaffected by this choice.

### 3.2.3 ICL contribution as a function of infall stellar mass

By combining our infall stellar mass function fit ($\Phi(M_*)$; from Fig. 5) with our fit for the liberated fraction as a function of infall stellar mass ($f_{\text{lib}}(M_*)$; from Fig. 4), and scaling by infall stellar
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Stellar masses at infall, $M_*$, for galaxies that fall into any of the 11 clusters after $z \sim 2$. A Schechter function (equation 1) is fitted using only bins with $M_* > 10^9 M_{\odot}$; bins used for fitting are shown in lighter orange and bins ignored for fitting are shown in darker red. The best-fitting parameters are presented in the upper right with the corresponding fit line plotted over the bins. The error bars and shaded regions around the fit line indicate estimated uncertainties based on bootstrapping.</figcaption></figure>

mass ($M_*$), we produce a fit for the fraction of the stacked, not pre-processed ICL of the 11 clusters associated with progenitors of a particular infall stellar mass, i.e.

$$f(M_*) = \frac{M_* f_{\text{lib}}(M_*) \Phi(M_*)}{\int M_* f_{\text{lib}}(M_*) \Phi(M_*) \, dM_*} \quad (2)$$

The result is shown in Fig. 6 (normalized to per dex in infall stellar mass), with dotted lines used to indicate extrapolation beyond the fit range of the constituent functions (which allows us to compensate for the resolution limit of the simulation). The equivalent analysis for the combined BCG + ICL system is also shown in Fig. 6 with a fainter blue line.

We find no clear trend (and considerable scatter) between the infall stellar mass of each galaxy and the fraction of the stars associated with that infaller that were formed during or after cluster infall. We therefore assume a constant ratio between the total stellar mass associated with a progenitor at $z \approx 0$ and its infall stellar mass, and consequently neglect post-infall star formation in equation (2) as this constant ratio has no impact on the shape of the curve seen in Fig. 6. The median (mean) ratio between the total associated stellar mass at $z \approx 0$ and the infall stellar mass of a progenitor for all $\sim 4000$ infaller galaxies across all 11 clusters was $\sim 1.01$ ($\sim 1.18$).

It can be seen in Fig. 6 that the dominant contributors of stars to the ICL are galaxies with infall stellar masses between $\sim 10^{10.5}$ and $\sim 10^{11} M_{\odot}$. This is in good agreement with the general conclusion of most prior theoretical studies (e.g. Contini et al. 2014, 2019; Chun et al. 2023, 2024; Ahvazi et al. 2024b) – that approximately Milky Way mass galaxies are the main progenitors of the ICL.

We highlight that the location of the peak in Fig. 6 – at $\log_{10}(M_*/M_{\odot}) = 10.94^{+0.13}_{-0.07}$ (uncertainty estimated by bootstrapping) – is very close to the location of the ‘knee’ of the Schechter function fitted in Fig. 5 ($\log_{10}(M_k/M_{\odot}) = 11.13^{+0.08}_{-0.07}$). Additionally, we note that the position of the peak in Fig. 6 varies very little

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Fractional contribution of galaxies with different stellar masses, $M_*$, on cluster infall to the build-up of the (not pre-processed) ICL, stacked across all 11 clusters. Dotted lines indicate extrapolation beyond the range of masses used for fitting. The shaded regions indicate estimated uncertainties based on bootstrapping. The faint line shows the equivalent fractional contribution of infalling galaxies to the combined BCG + ICL system.</figcaption></figure>

between the main version of our analysis and the alternative version in which the BCG and ICL are not separated (peak at $\log_{10}(M_*/M_{\odot}) = 11.2^{+0.1}_{-0.3}$), despite the two fits used for $f_{\text{lib}}$ as a function of infall mass being significantly different (as can be seen in Fig. 4). This is a consequence of both these fits for $f_{\text{lib}}$ only varying by a factor of order unity over a nearly four orders of magnitude change in infall stellar mass, with the fitted infall stellar mass function (shown in Fig. 5) instead varying by several orders of magnitude over this same infall mass range. As such, the specific fit used for $f_{\text{lib}}$ as a function of infall mass has only a minor effect on the shape of the curve seen in Fig. 6, which is instead governed almost entirely by the infall stellar mass function.

Although Fig. 6 shows that the expected contribution of infalling objects with $M_* \gtrsim 10^{11} M_{\odot}$ to the ICL rapidly declines with any further increase in mass, this is entirely due to how rarely we expect such supremely massive objects to fall into clusters (as can be noted from Fig. 5); much of the mass of such a massive galaxy would need to be assembled through mergers, and so can only be assembled at the heart of a rich group or cluster, hence such objects only join clusters very rarely (as part of major cluster mergers). However, if such a massive object does join a cluster we anticipate that it would produce a significant contribution to the ICL of that cluster, even before considering the pre-processed ICL that should accompany it. This view agrees with the findings of Cooper et al. (2015) and Harris et al. (2017): that the ICL of a cluster can be heavily influenced by only a small number of massive progenitors, being largely built up stochastically as these rare massive objects infrequently join the cluster. This stochasticity may also explain the large dispersion found in the fractions of $z \approx 0$ ICL stars that were formerly part of the BCG (shown in Fig. 3), as violent mergers between the BCG and the most massive infallers are expected to be the primary mechanism for liberating BCG stars into the ICL.

It warrants restating that the analysis depicted in Figs 4 and 6 ignores pre-processed ICL. If this analysis were repeated with the
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Top: Normalized cumulative contribution to the stacked (not pre-processed) ICL of the 11 clusters by galaxies with stellar mass, $M_*$, on cluster infall less than M after correcting for resolution effects (shown in blue). Dotted lines are used to indicate extrapolation beyond the range of masses used for fitting. Shown in red (for comparison) is a cumulative sum plot for the stacked (not pre-processed) ICL taken directly from the 11 simulated clusters (with no corrections made for resolution effects), normalized relative to the corrected curve. The shaded regions indicate estimated uncertainties based on bootstrapping. Bottom: Same as the top panel but for the combined BCG + ICL system.</figcaption></figure>

pre-processed ICL from accreted groups/clusters instead attributed to the central galaxies of those groups/clusters, we would expect the shapes of the fits seen in Figs 4 and 6 to be substantially altered, with this extra attributed ICL mass enhancing the fractional contribution of infalling galaxies with $M_* \gtrsim 10^{11} M_\odot$ (as well as the $f_{\text{lib}}$ values determined for these galaxies), potentially shifting the location of the peak seen in Fig. 6 to a slightly higher mass.

The top panel of Fig. 7 shows the (normalized) cumulative contribution to the stacked, not pre-processed ICL of all 11 clusters from progenitor objects with infall stellar masses, $M_*$, less than M: the red curve is derived directly from the simulation (without correcting for resolution effects), and the blue curve is derived from the fitted function depicted in Fig. 6. Dashed lines are once again used

to indicate extrapolation beyond the range of infall stellar masses used for fitting. The corrected curve is normalized to unity, and the uncorrected curve is normalized relative to the corrected curve (adjusted by a constant factor equal to the mean ratio between total associated stellar mass at $z \approx 0$ and progenitor infall stellar mass to account for star formation during and after cluster infall), so that the total ICL mass encapsulated by both curves can be directly compared. The bottom panel of Fig. 7 shows the same as the top panel but for the combined BCG + ICL system instead.

It can be seen in the top panel of Fig. 7 that ~50 per cent of the (not pre-processed) ICL is predicted to come from galaxies with infall stellar masses within half a dex of $10^{11} M_\odot$, and ~90 per cent is predicted to come from progenitors with stellar masses $\gtrsim 10^9 M_\odot$. Even after compensating for resolution effects using our extrapolated fitting function approach, that galaxies with $M_* < 10^9 M_\odot$ contribute only $\lesssim 10$ per cent of the total ICL mass suggests that the bulk properties of the ICL should be insensitive to contributions from these low-mass galaxies. This also tentatively implies that once a simulation achieves sufficient resolutions to resolve galaxies with $M_* \sim 10^9 M_\odot$, further improvements in resolution will yield no significant changes in the bulk properties of the ICL and only minor changes ($\lesssim 10$ per cent) in the total ICL masses found.

When Puchwein et al. (2010) previously investigated the effects of altering the mass resolution of a hydrodynamical simulation on the ICL masses identified in simulated clusters, beyond a star-particle mass of $\sim 10^7 M_\odot$ the ICL masses of their simulated clusters appear to have converged — a suggestion which supports our findings here. Likewise, our findings agree well with those of Ahvazi et al. (2024b), who found that $\gtrsim 90$ per cent of the ICL in the groups and clusters from the higher resolution TNG50 simulation that they investigated typically came from progenitor objects with stellar masses greater than $\sim 10^9 M_\odot$ and that $\gtrsim 50$ per cent typically came from progenitors with stellar masses greater than $\sim 10^{10} M_\odot$ (though with significant dispersion).

It is worth noting that the shaded regions in the Fig. 7 (indicating estimated uncertainties based on bootstrapping) only significantly broaden for stellar masses beyond $\sim 10^{11} M_\odot$, indicating that the substitution of even only a very small number of progenitor objects this massive for less massive ones (or vice versa) can appreciably alter the total ICL mass of all 11 clusters combined – emphasizing the significant ICL contributions objects this massive can make when they do join clusters. We do note, however, that the seeming strong agreement between the corrected and uncorrected curves at the high mass end in the top panel of Fig. 7 is somewhat coincidental: the Schechter function fitted in Fig. 5 slightly exaggerates the scarcity of the most massive objects, with this coincidentally diminishing the expected overall contribution to the ICL from high mass infallers by approximately as much as unresolved low-mass galaxies are expected to contribute. The mass function fit slightly exaggerating the scarcity of high mass infallers is also why the corrected curve in the bottom panel of Fig. 7 falls slightly below the uncorrected curve at the high-mass end.

We emphasize that the analysis represented in Fig. 7 is for the combined ICL/BCG + ICL of all 11 clusters stacked together. We warn that if this analysis were repeated on any individual cluster the stochastic nature of cluster assembly could cause the shape of the resulting curve to significantly diverge beyond what is encompassed by the shaded regions in Fig. 7, particularly for very high mass progenitor objects (scarce few of which should fall into such a $10^{14} M_\odot$ cluster and so the entire contribution of the high infall mass regime may come from only one or two objects). Equivalent curves for each of the 11 clusters individually are presented in Appendix A.
# 4 CONCLUSIONS

Using the HORIZON-AGN simulation (Dubois et al. 2014), we have studied the assembly of the ICL. We investigated ICL assembly in 11 simulated clusters with $z \approx 0$ DM halo masses between $\sim 1 \times 10^{14}$ and $\sim 7 \times 10^{14}$ M$_{\odot}$, tracking the stars of galaxies that fell into these clusters over the past $\sim 10$ Gyr in order to quantify the differing contributions made to the $z \approx 0$ ICL by galaxies with differing stellar masses on cluster infall. Our main findings can be summarized as follows:

(i) On average, 44 per cent of the $z \approx 0$ ICL stars in the 11 studied clusters were previously part of a satellite galaxy within the cluster halo (but never part of the BCG); another 36 per cent on average were formerly part of the BCG (Fig. 3). Such significant fractions in both instances indicate both satellite stripping and violent BCG mergers play significant roles in ICL assembly.

(ii) In situ formation of stars directly into the ICL appears insignificant in the HORIZON-AGN simulation; we discern an upper limit of order 0.1 per cent on the fraction of the stacked ICL of the 11 clusters formed by this channel. The remainder of the ICL is virtually all pre-processed. The average fraction of $z \approx 0$ ICL stars that were pre-processed and never part of a cluster galaxy was 16 per cent (Fig. 3), though in the least relaxed and most massive cluster $\gtrsim 40$ per cent of $z \approx 0$ ICL stars fell into this category. Inversely, the lowest pre-processed fraction was seen in the most relaxed and least massive cluster.

(iii) A galaxy with less stellar mass on cluster infall can be expected to contribute a greater percentage of its stars to the ICL than a more massive infalling galaxy. The most massive galaxies that join clusters do, however, typically contribute a very high fraction of their infall stellar mass to the combined BCG + ICL system ($\gtrsim 60$ per cent; Fig. 4).

(iv) Roughly half of the stacked, not pre-processed ICL of the 11 clusters came from progenitor objects with infall stellar masses within half a dex of $10^{11}$ M$_{\odot}$ (Fig. 7). This is the case even after compensating for resolution effects. As the mean fraction of stars liberated from infalling galaxies only varied by a factor of order unity between infall masses of $\sim 10^9$ and $\sim 10^{12}$ M$_{\odot}$, the location of this peak is largely determined by the infalling galaxy stellar mass function instead. Any object as or more massive than the Milky Way that joins a cluster is expected to make a sizeable contribution to the ICL of that cluster and – as such massive objects joining clusters should be infrequent events – this suggests ICL assembly may be chiefly stochastic.

(v) 90 per cent of the bulk ICL of the 11 clusters which was not pre-processed could be attributed to progenitors with infall stellar masses $\gtrsim 10^9$ M$_{\odot}$, even after compensating for resolution effects (Fig. 7). As the ICL appears virtually complete even without any contribution from less massive galaxies (at least within the validity of our stellar mass function extrapolation), we expect the bulk properties of the ICL to be insensitive to the low-mass galaxy population.

In summary, we have shown that the main progenitors of the ICL should be massive galaxies, and more specifically those with stellar masses on cluster infall close to that of the Milky Way. Despite less massive galaxies being more numerous and adding a larger fraction of their stellar mass on cluster infall to the ICL, the overall contribution of this population to the ICL of an $\sim 10^{14}$ M$_{\odot}$ halo mass cluster appears rather small. Conversely, though very massive galaxies join clusters only rarely, they make major contributions to the mass of the ICL when they do so, despite this conferred ICL mass only being a small fraction of their infall stellar mass.

Although we do not expect the low-mass galaxy population to meaningfully influence the bulk properties of the ICL, the potential remains for this population to have a significant impact on the radial dependence of ICL properties. We expect much of the mass of more massive infalling galaxies to be deposited close to the centre of the cluster (as implied in Fig. 4), potentially allowing the ICL in the outskirts of a cluster to be dominated by stars from less massive progenitors – an idea supported observationally by prior studies noting metallicity and colour gradients in the ICL (e.g. DeMaio et al. 2018; Gu et al. 2020; Golden-Marx et al. 2023). We intend to investigate this further in a future work.

# ACKNOWLEDGEMENTS

We thank the anonymous referee for their careful reading of our original manuscript and for their constructive comments which have helped to improve the quality and clarity of the presented work. We also thank Jesse B. Golden-Marx for helpful discussions and comments. This work made use of NUMPY (Harris et al. 2020), MAT-PLOTLIB (Hunter 2007), SCIPY (Virtanen et al. 2020), H5PY (Collette et al. 2023), ASTROPY (The Astropy Collaboration 2022), and PYGAM (Servén et al. 2018). HJB gratefully acknowledges support from the UK Science and Technology Facilities Council (STFC) under grant ST/Y509437/1. GM, FRP, and NAH gratefully acknowledge support from the UK STFC under grant ST/X000982/1. NAH gratefully acknowledges support from the Leverhulme Trust through a Research Leadership Award. YMB gratefully acknowledges financial support from the Swiss National Science Foundation (SNSF) under funding reference 200021_213076. This work was granted access to the HPC resources of CINES under allocations 2013047012, 2014047012, and 2015047012 made by GENCI. This work has made use of the Infinity cluster, hosted by the Institut d'Astrophysique de Paris. We warmly thank S. Rouberol for running it smoothly.

# DATA AVAILABILITY

The raw 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 generated from this work are available upon reasonable request from the corresponding author.

# REFERENCES

* Ahvazi N., Sales L. V., Navarro J. F., Benson A., Boselli A., D'Souza R., 2024a, preprint (arXiv:2403.04839)

* Ahvazi N., Sales L. V., Doppel J. E., Benson A., D'Souza R., Rodriguez-Gomez V., 2024b, MNRAS, 529, 4666

* Alonso Asensio I., Dalla Vecchia C., Bahé Y. M., Barnes D. J., Kay S. T., 2020, MNRAS, 494, 1859

* Aubert D., Pichon C., Colombi S., 2004, MNRAS, 352, 376

* Bondi H., 1952, MNRAS, 112, 195

* Brough S. et al., 2024, MNRAS, 528, 771

* Chun K., Shin J., Smith R., Ko J., Yoo J., 2023, ApJ, 943, 148

* Chun K., Shin J., Ko J., Smith R., Yoo J., 2024, ApJ, 969, 142

* Collette A. et al., 2023, h5py/h5py: 3.8.0-aarch64-wheels, available at: https://doi.org/10.5281/zenodo.7568214

* Contini E., 2021, Galaxies, 9, 60

* Contini E., De Lucia G., Villalobos A., Borgani S., 2014, MNRAS, 437, 3787

* Contini E., Yi S. K., Kang X., 2018, MNRAS, 479, 932

* Contini E., Yi S. K., Kang X., 2019, ApJ, 871, 24
Contini E., Rhee J., Han S., Jeon S., Yi S. K., 2024a, AJ, 167, 7

Contini E., Han S., Jeon S., Rhee J., Yi S. K., 2024b, ApJ, 962, L10

Contreras-Santos A. et al., 2024, A&A, 683, A59

Cooper A. P., Gao L., Guo Q., Frenk C. S., Jenkins A., Springel V., White S. D. M., 2015, MNRAS, 451, 2703

Deason A. J. et al., 2020, MNRAS, 500, 4181

DeMaio T., Gonzalez A. H., Zabludoff A., Zaritsky D., Connor T., Donahue M., Mulchaey J. S., 2018, MNRAS, 474, 3009

Dolag K., Borgani S., Murante G., Springel V., 2009, MNRAS, 399, 497

Dubois Y., Devriendt J., Slyz A., Teyssier R., 2012, MNRAS, 420, 2662

Dubois Y. et al., 2014, MNRAS, 444, 1453

Dubois Y., Peirani S., Pichon C., Devriendt J., Gavazzi R., Welker C., Volonteri M., 2016, MNRAS, 463, 3948

Furnell K. E. et al., 2021, MNRAS, 502, 2419

Gavazzi R., Adami C., Durret F., Cuillandre J.-C., Ilbert O., Mazure A., Pelló R., Ulmer M. P., 2009, A&A, 498, L33

Golden-Marx J. B. et al., 2023, MNRAS, 521, 478

Gonzalez A. H., George T., Connor T., Deason A., Donahue M., Montes M., Zabludoff A. I., Zaritsky D., 2021, MNRAS, 507, 963

Gu M. et al., 2020, ApJ, 894, 32

Gullieuszik M. et al., 2020, ApJ, 899, 13

Haardt F., Madau P., 1996, ApJ, 461, 20

Harris K. A. et al., 2017, MNRAS, 467, 4501

Harris C. R. et al., 2020, Nature, 585, 357

Hatton S., Devriendt J. E. G., Ninin S., Bouchet F. R., Guiderdoni B., Vibert D., 2003, MNRAS, 343, 75

Hunter J. D., 2007, Comput. Sci. Eng., 9, 90

Jiménez-Teja Y., Dupke R. A., Lopes P. A. A., Dimauro P., 2024, ApJ, 960, L7

Joo H., Jee M. J., 2023, Nature, 613, 37

Kaviraj S. et al., 2017, MNRAS, 467, 4739

Kennicutt R. C. Jr, 1998, ApJ, 498, 541

Kluge M. et al., 2020, ApJS, 247, 43

Kluge M., Bender R., Riffeser A., Goessl C., Hopp U., Schmidt M., Ries C., 2021, ApJS, 252, 27

Kluge M. et al., 2024, preprint (arXiv:2405.13503)

Knebe A. et al., 2011, MNRAS, 415, 2293

Komatsu E. et al., 2011, ApJS, 192, 18

Martel H., Barai P., Brito W., 2012, ApJ, 757, 48

Melnick J., Giraud E., Toledo I., Selman F., Quintana H., 2012, MNRAS, 427, 850

Mihos J. C., 2015, Proc. IAU Symp. 11, Intragroup and Intracluster Light. Kluwer, Dordrecht, p. 27

Mihos J. C., 2019, preprint (arXiv:1909.09456)

Mihos J. C., Harding P., Feldmeier J. J., Rudick C., Janowiecki S., Morrison H., Slater C., Watkins A., 2017, ApJ, 834, 16

Montenegro-Taborda D., Rodriguez-Gomez V., Pillepich A., Avila-Reese V., Sales L. V., Rodríguez-Puebla A., Hernquist L., 2023, MNRAS, 521, 800

Montes M., 2022, Nat. Astron., 6, 308

Montes M., Trujillo I., 2014, ApJ, 794, 137

Montes M., Trujillo I., 2018, MNRAS, 474, 917

Montes M., Trujillo I., 2019, MNRAS, 482, 2838

Montes M., Trujillo I., 2022, ApJ, 940, L51

Montes M., Brough S., Owers M. S., Santucci G., 2021, ApJ, 910, 45

Morishita T., Abramson L. E., Treu T., Schmidt K. B., Vulcani B., Wang X., 2017, ApJ, 846, 139

Murante G., Giovalli M., Gerhard O., Arnaboldi M., Borgani S., Dolag K., 2007, MNRAS, 377, 2

Nelson D. et al., 2021, Computational Astrophysics and Cosmology, 6, 29

Onions J. et al., 2012, MNRAS, 423, 1200

Peebles P. J. E., Ratra B., 2003, Rev. Mod. Phys., 75, 559

Pillepich A. et al., 2018, MNRAS, 475, 648

Planck Collaboration XVI, 2014, A&A, 571, A16

Puchwein E., Springel V., Sijacki D., Dolag K., 2010, MNRAS, 406, 936

Purcell C. W., Bullock J. S., Zentner A. R., 2007, ApJ, 666, 20

Ragusa R. et al., 2023, A&A, 670, L20

Read J. I., Wilkinson M. I., Evans N. W., Gilmore G., Kleyna J. T., 2006, MNRAS, 366, 429

Rudick C. S., Mihos J. C., McBride C., 2006, ApJ, 648, 936

Rudick C. S., Mihos J. C., Frey L. H., McBride C. K., 2009, ApJ, 699, 1518

Rudick C. S., Mihos J. C., McBride C. K., 2011, ApJ, 732, 48

Savitzky A., Golay M. J. E., 1964, Anal. Chem., 36, 1627

Schechter P., 1976, ApJ, 203, 297

Servén d., Brummitt C., Abedi H., hlink, 2018, dswah/pyGAM: v0.8.0, available at: https://doi.org/10.5281/zenodo.1476122

Somerville R. S., Davé R., 2015, ARA&A, 53, 51

Sutherland R. S., Dopita M. A., 1993, ApJS, 88, 253

Tang L., Lin W., Cui W., Kang X., Wang Y., Contini E., Yu Y., 2018, ApJ, 859, 85

Tang L., Lin W., Wang Y., Napolitano N. R., 2021, MNRAS, 508, 3321

Tang L., Lin W., Wang Y., Li J., Lan Y., 2023, ApJ, 959, 104

Teyssier R., 2002, A&A, 385, 337

The Astropy Collaboration, 2022, ApJ, 935, 167

Tweed D., Devriendt J., Blaizot J., Colombi S., Slyz A., 2009, A&A, 506, 647

Urban O., Werner N., Simionescu A., Allen S. W., Böhringer H., 2011, MNRAS, 414, 2101

Virtanen P. et al., 2020, Nat. Methods, 17, 261

Vogelsberger M., Marinacci F., Torrey P., Puchwein E., 2020, Nat. Rev. Phys., 2, 42

Williams B. F. et al., 2007, ApJ, 656, 756

Willman B., Governato F., Wadsley J., Quinn T., 2004, MNRAS, 355, 159

Yoo J., Ko J., Kim J.-W., Kim H., 2021, MNRAS, 508, 2634

Yoo J. et al., 2024, ApJ, 965, 145

Zhang Y. et al., 2019, ApJ, 874, 165

Zwicky F., 1937, ApJ, 86, 217

Zwicky F., 1951, PASP, 63, 61

## APPENDIX A: INDIVIDUAL CLUSTER CUMULATIVE ICL CONTRIBUTION CURVES

Fig. 7 shows the cumulative contribution to the ICL/BCG + ICL from increasingly massive progenitors for the stacked ICL/BCG + ICL of all 11 clusters combined. The uncertainties included in that plot (estimated based on bootstrapping and indicated by the shaded areas) are thus not thought to be representative of the anticipated cluster-to-cluster scatter should an equivalent analysis be performed on individual clusters. We present this equivalent analysis for each of the 11 clusters individually in Fig. A1, showing for each cluster the fraction of the (not pre-processed) $z \approx 0$ ICL mass (top panel) or BCG + ICL mass (bottom panel) linked to progenitors with infall stellar mass, $M_*$, less than $M$ (taken directly from the simulation with no corrections made for resolution effects). The grey shaded regions are bounded by the 16th and 84th percentiles of the individual cluster curves. To facilitate comparison, the corresponding (uncorrected) curves for the ICL/BCG + ICL of all 11 clusters stacked from Fig. 7 (normalized to unity) are also shown. The ID numbers referenced are those presented in Table 1.

In the top panel of Fig. A1, the individual cluster curves for IDs 9 and 19 are worthy of particular note. Only ~30 per cent of the total $z \approx 0$ ICL mass of ID 9 is from progenitors with $M_* \gtrsim 10^{10.5} M_\odot$, and so the ICL of this specific cluster is primarily sourced from lower mass progenitors. Conversely, ~50 per cent of the total $z \approx 0$ ICL mass of ID 19 is from progenitors with $M_* \gtrsim 10^{11}$ hence the ICL of this particular cluster is dominated by contributions from infalling galaxies more massive than the Milky Way. These two clusters serve as prime examples of how the main progenitors of ICL stars in individual clusters can diverge from what
<figure><img src="image_9.png" /><figcaption>**Figure A1.** Top: Normalized cumulative contribution to the (not pre-processed) $z \approx 0$ ICL of each of the 11 clusters by galaxies with stellar mass, $M_*$, on cluster infall less than $M$ (with no corrections made for resolution effects). The grey shaded region is bounded by the 16th and 84th percentiles of the individual cluster curves. The corresponding curve for the ICL of all 11 clusters stacked is also shown (as the thicker red line). Bottom: Same as the top panel but for the combined BCG + ICL system.</figcaption></figure>

is expected from considering an ensemble of clusters due to the stochastic nature of cluster assembly. Several instances can also be seen in the top panel of Fig. A1 of $\gtrsim 15$ per cent of the (not pre-processed) ICL mass of an individual cluster being linked to a single massive ($M_* \gtrsim 10^{11} M_\odot$) progenitor. This supports the perspective that – although they generally join clusters too infrequently to be the dominate overall contributors of ICL stars – individual massive galaxies can potentially make major ICL contributions when they do (rarely) join clusters.

This paper has been typeset from a TeX/LATEX file prepared by the author.