Assembly of the Intracluster Light in the Horizon-AGN Simulation

The diffuse stellar component of galaxy clusters made up of intergalactic stars is termed the intracluster light (ICL). Though 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 $\sim1\times 10^{14}$ M$_{\odot}$ and $\sim7\times 10^{14}$ M$_{\odot}$ at $z\approx0$) 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\approx0$ 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.


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 theorised and then discovered by Zwicky (1937Zwicky ( , 1951)), this low surface-brightness cluster component has been traced out to ≳ 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 et al. 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 et al. 2006, Mihos et al. 2017), with this final formation channel commonly called ICL pre-processing.Though once considered potentially significant formation channels, total disruption of dwarf galaxies (e.g.Purcell et al. 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 et al. 2012, Melnick et al. 2012, DeMaio et al. 2018, and Gullieuszik et al. 2020; ★ E-mail: Harley.Brown@nottingham.ac.uk 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  ∼ 0, having potentially superseded the overall contribution from mergers sometime around  ∼ 1 (see Joo &Jee 2023, Jiménez-Teja et al. 2024 andreferences 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  ≈ 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, Contini et al. 2024a, and Chun et al. 2024).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 and so should be distributed similarly within the cluster.It has thus been proposed that the ICL might be used as a visible dark matter tracer (e.g.Montes & Trujillo 2019, Alonso Asensio et al. 2020, and 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 contextualised, necessitating investigations into the predictions of numerical simulations concerning the ICL.
Though 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  ∼ 0 (though not unanimously), whilst 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 onto an existing N-body dark matter simulation, whilst hydrodynamical simulations instead numerically solve equations for gravity, thermodynamics, and hydrodynamics, for dark matter, gas, and stars simultaneously.Though 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  = 0 should be galaxies with stellar masses ≳ 10 10.5 M ⊙ (Contini et al. 2014(Contini et al. , 2019) ) and that the ICL can be heavily influenced by a small number of massive progenitors -between ∼ 1 and ∼ 10 per cent the mass of the BCG (Cooper et al. 2015).Harris et al. (2017) re-simulated a Fornax-like cluster (halo mass ∼ 4 × 10 13 M ⊙ ) from a dark matter only N-body simulation by replacing haloes infalling after  = 1.65 with full galaxy models, and found > 60 per cent of the  = 0 ICL of this cluster to have been sourced from just two massive objects (stellar masses ∼ 5 × 10 10 M ⊙ ).By employing a similar "galaxy replacement technique" Chun et al. (2023) studied 6 clusters (virial masses of order 10 14 M ⊙ ) and found the ICL in all but the most relaxed cluster to be dominated by stars from galaxies with infall stellar masses ≳ 10 10 M ⊙ .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 ( 200 [M ⊙ ]) < 14.8) to be galaxies with stellar masses between 10 10 M ⊙ and 10 11 M ⊙ on cluster infall, though also noted that the ICL of particularly massive or un-relaxed clusters could have significant contributions from galaxies with infall stellar masses > 10 11 M ⊙ as well.Ahvazi et al. (2024b) conducted a study on 39 groups and clusters (virial masses between 5 × 10 12 M ⊙ and 2 × 10 14 M ⊙ ) 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 M ⊙ and 10 11 M ⊙ , 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 ∼ 0.1 and ∼ 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 M ⊙ and 10 10 M ⊙ .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 ⊙ .For example, Montes et al. (2021) suggested -based on considerations of colour profiles -that the ICL of A85 ( ≈ 0.05,  200 ≈ 1.7 × 10 15 M ⊙ ) was built up mainly by the stripping of satellites with stellar masses of order 10 10 M ⊙ .Similarly, Montes & Trujillo (2014, 2018) found (also based on colours) that the ICL of several massive ( 200 > 10 14.3 M ⊙ ) clusters at 0.3 <  < 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×10 13 ≲  500 / M ⊙ ≲ 9×10 14 ) they considered at 0.29 <  < 0.89 should have originated from galaxies with stellar mass > 10 10.4 M ⊙ .
Conversely, Morishita et al. (2017) found the ICL of six Hubble Frontier Field clusters (0.3 <  < 0.6,  500 ≳ 10 15 M ⊙ ) to be dominated by moderately old stars (∼ 1 to 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 ∼ 10 9.5 M ⊙ since  ∼ 1, instead of more massive galaxies.Likewise, Gu et al. (2020) studied portions of the ICL in the Coma cluster ( ≈ 0.024,  200 ≈ 5.1 × 10 14 M ⊙ : 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 ( ≈ 0.004,  200 ≈ 1.4 × 10 14 M ⊙ : Urban et al. 2011) -and suggested that these stars may have primarily originated from the accretion of low-mass galaxies (stellar masses less than ∼ 3 × 10 9 M ⊙ ).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  = 0 in 11 clusters with dark matter halo masses between ∼ 1 × 10 14 M ⊙ and ∼ 7 × 10 14 M ⊙ at  ≈ 0. By tracking stars associated with galaxies that fall into these clusters during this time, we are able to link much of the  ≈ 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 9 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 organised as follows: In Section 2 we present the details of the Horizon-AGN simulation (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 (2.2).In Section 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 (3.1), then on the contributions of progenitor objects with differing stellar masses on cluster infall to the ICL (3.2).We conclude by summarising the main results of our investigation in Section 4. Throughout this paper, we assume a ΛCDM cosmology with ℎ ≡  0 /100 km s −1 Mpc −1 = 0.704.Unless stated otherwise, distances are given in proper rather than co-moving units.

Horizon-AGN
The full details of the Horizon-AGN simulation can be found in Dubois et al. (2014Dubois et al. ( , 2016)), and are reproduced only in brief here.
The Horizon-AGN simulation is a cosmological-volume hydrodynamical simulation which employs 1024 3 dark matter (DM) particles (each with mass 8 × 10 7 M ⊙ ) in a cubic volume of side-length 100 ℎ −1 Mpc (co-moving) with periodic boundary conditions.A standard ΛCDM cosmological model (Peebles & Ratra 2003) is used with total matter density Ω m = 0.272, baryon density Ω b = 0.045, dark energy density Ω Λ = 0.728, dark matter power spectrum amplitude  8 = 0.81, Hubble constant  0 = 70.4km s −1 Mpc −1 , and spectral index  s = 0.967, compatible with the Seven-Year Wilkinson Microwave Anisotropy Probe cosmology (Komatsu et al. 2011), and also that of the Planck Collaboration et al. ( 2014) within ∼ 10 per cent relative variation.The simulation utilises Ramses (Teyssier 2002) -an adaptive mesh refinement Eulerian hydrodynamics code -with an initially uniform 1024 3 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 Δ = 1 kpc after seven levels of refinement.
Gas heating from a uniform UV background begins after  = 10, following Haardt & Madau (1996), and gas is allowed to cool to 10 4 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 onto 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 −3 .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 6 M ⊙ .
The Horizon-AGN simulation was not calibrated to the local Universe, except for choosing blackhole feedback parameters.Despite this, Kaviraj et al. (2017) reported that the simulation is in generally good agreement with observations (using observational data from 0 <  < 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 10.5 M ⊙ at all epochs is noted compared to observations, though agreement can be achieved by considering observational uncertainties due to cosmic variance.

AdaptaHOP and Treemaker
For identifying structures of both DM particles (i.e.haloes) and star particles (i.e.galaxies), the AdaptaHOP structure finder (Aubert et al. 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 6 M ⊙ , this gives a minimum (detectable) galaxy stellar mass of ∼ 10 8 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× 178 of the halo centre (where  178 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 maximises the fraction of particles shared by the would-be linked structures.Between  ≈ 3 and  ≈ 0 the mergers trees we use in this study were built using a time resolution of ∼ 0.05 Gyr.
For this study the galaxy considered to be the BCG of each cluster is selected at  ≈ 0 ( = 0.0556) as the most massive galaxy within 0.1× 178 of the cluster halo centre.Prior to  ≈ 0 the galaxy classified as the BCG of the cluster in each snapshot is the main progenitor of the  ≈ 0 BCG.We place the border of each cluster in each snapshot at a distance away from the BCG equal to  178 of the cluster halo (see Section 2.2.2) and classify all galaxies within this border besides the BCG as satellite galaxies.

Tagging and tracking the stars of infalling galaxies
For this investigation we select 11 clusters from the Horizon-AGN simulation, with  ≈ 0 ( = 0.0556) DM halo masses (within  178 i.e.  178 ) in the range ∼ 1 × 10 14 to ∼ 7 × 10 14 M ⊙ 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 ∼ 1 Gyr, down to  ≈ 0 ( = 0.0556) from as far back as when a main progenitor of the cluster BCG is first identified (which is always before  ∼ 1.5 and almost always before  ∼ 2) or  ∼ 3 -whichever is later.As later described in Section 3.1, we never find more than ∼ 0.1 per cent of the  ≈ 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 andContini et al. 2024a) which typically found the vast majority of the  = 0 ICL mass of similarly massive clusters to have still been contained in galaxies at  ∼ 1.
We place the border of each cluster at  178 and consider galaxies to have become satellites of a cluster once they pass  178 .Rather than use the values provided by AdaptaHOP directly, for each cluster we instead fit a monotonic cubic spline to the  178 values of each cluster from AdaptaHOP as a function of look-back time (using snapshots with ∼ 1 Gyr coarse spacing), allowing us to account for and smooth out any temporary and nonphysical fluctuations in cluster size (such as those caused by cluster mergers).
To identify the progenitor objects of the  ≈ 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  ≈ 0. Infalling galaxies are identified as galaxies at  >  178 in any coarse snapshot whose main descendant in the next coarse snapshot is at  <  178 .This infaller galaxy identification is done chronologically, so galaxies which have previously fallen into the cluster but whose orbit apocentre still lies outside  178 can be identified (and repeatcounting 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 ∼ 4 per cent of all  ≈ 0 ICL stars across all 11 clusters (and ∼ 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 ∼ 4000 infalling galaxies tracked, we modify the merger-trees of ∼ 15 per cent after initial infall -most commonly to patch breaks while galaxies are partway through merging with the BCG.

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  ≈ 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 Figure 1; the left and centre panels show the projected positions for an arbitrary line-of-sight at  ≈ 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), whilst 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 Figure 1 is a velocity dispersion plot (speed relative to the BCG plotted against distance from BCG centre, normalised by the 3D velocity dispersion of all cluster star particles and  178 respectively) for the same ICL star particles as shown in the centre panel of Figure 1 and employing the same scheme for colours.Comparison between the centre and 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  ≈ 0 cluster.The DM and stellar masses are those within  178 at  ≈ 0 and include substructure, and  BCG ,  Sat , and  ICL are the fractions of this stellar mass in the BCG, satellite galaxies, and ICL respectively and are also depicted in Figure 2. The ICL sub-fractions given in the last three columns are also depicted in Figure 3 right panels shows correspondence between features in position and velocity space, as expected for ICL stars recently sourced from the same infalling galaxy.

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 16 th and 84 th percentiles from these repeat analyses.

Stellar mass budget and ICL pre-processing
Figure 2 depicts how the  ≈ 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  ≈ 0 ICL (BCG + ICL) stellar mass fraction of the clusters is 13 (45) per cent, the 16 th percentile is 11 (40) per cent, and the 84 th 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 andKluge et al. 2021) -and agree fairly well with the findings of other theoretical studies employing similar ICL definitions, such as that by Rudick et al. (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 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.≳ 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 ∼ 10 14 M ⊙ (e.g.Pillepich et al. 2018, Brough et al. 2024, and 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, and 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.
Figure 3 depicts the fractions of  ≈ 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 of ICL stars that were previously part of any cluster galaxy are also included.
The median fraction of  ≈ 0 ICL stars formerly in satellite galaxies (and never part of the BCG) is 44 per cent (16 th percentile 43 per cent and 84 th percentile 52 per cent), and the median fraction that were formerly part of the BCG is 36 per cent (16 th percentile 26 per cent and 84 th percentile 41 per cent).That never less than 40 per cent of the  ≈ 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  ≈ 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  200 ≥ 10 14 M ⊙ clusters in the TNG300 simulation (Nelson et al. 2021), they found the mean fraction of ICL stars (at  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 IGL of infalling groups) are shown in grey.Right: Scatter plot of the speeds, , of intergalactic stars (relative to the BCG and scaled by the 3D velocity dispersion of the cluster stars,   ) against distance from BCG centre, , (scaled by  178 ) for the same cluster and using the same colours as the middle panel. = 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.Fractions of  ≈ 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 preprocessed ICL) also shown.Dotted lines indicate 16 th and 84 th percentiles, and solid lines the maximum, minimum, and median values.

Satellite
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-Taborda et al. (2023) or Chun et al. (2024).
Using SAMs Contini et al. (2024a) obtained a median  = 0 ICL contribution from mergers of only ∼ 10 per cent for clusters with halo masses between ∼ 10 14 M ⊙ and ∼ 10 14.5 M ⊙ .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  ≈ 0 ICL stars identified by Adap-taHOP as previously being part of a cluster galaxy actually entered that cluster as part of a galaxy.For example, ∼ 20 per cent of the  ≈ 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 Figure 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  ≈ 0 ICL stars that have never been classed as part of a cluster galaxy is 16 per cent (16 th percentile 12 per cent and 84 th 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 M ⊙ and ∼ 10 14.5 M ⊙ .
The pre-processed fractions shown in Figure 3 feature one significant outlier: ID 13, the only cluster with a pre-processed fraction ≳ 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  ≈ 0 this cluster is uniquely part-way through two simultaneous major cluster mergers).It is also of note that the cluster with the lowest  ≈ 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 Figure 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, andBahé et al [in preparation]).We investigated these minor channels and find that at most ∼ 0.1 per cent of the stacked  ≈ 0 ICL of the 11 clusters could be from in-situ star formation 1 , and never find more than ∼ 0.1 per cent of the  ≈ 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  ≈ 0 ICL of all 11 clusters to be from unresolved galaxies with stellar masses < 10 9 M ⊙ .

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  ≈ 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,  lib , of that progenitor.Figure 4 shows the mean  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 ⊙ ) ≈ 8.2 and log( * /M ⊙ ) ≈ 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 fifteenth bin is shown in Figure 4.The equivalent analysis for the combined BCG + ICL system is also shown in faint.
As we expect any features in the  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  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  178 in consecutive coarse snapshots and never returning before  ≈ 0 -as well as those that only cross  178 for the first time at  ≈ 0. We exclude these galaxies as we consider them distinct from the population which fall into the clusters ≳ 1 Gyr prior to  ≈ 0 and are allowed ample opportunity to be processed by the cluster.By repeating our analysis with these galaxies included we verify the 1 We estimate the in-situ fraction as the fraction of  ≈ 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  ≈ 0 cluster (i.e outside  178 ), and have never (in any coarse snapshot) been classed as part of a stellar structure by AdaptaHOP.Any to-be  ≈ 0 ICL stars potentially formed outside of a galaxy but beyond  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.best-fitting function remains nearly unaffected, bar a small (order 0.01) approximately uniform shift to lower values of  lib .
From Figure 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 Figure 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  ≈ 0 ICL stars contributed by more massive progenitor galaxies having previously (between progenitor infall and  ≈ 0) been part of the BCG of their cluster (not shown in Figure 4): ∼ 25 per cent of all  ≈ 0 ICL stars (from all 11 clusters) linked to progenitors with infall stellar masses between 10 8.5 M ⊙ and 10 9.5 M ⊙ 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 M ⊙ and 10 11.5 M ⊙ .
The dominant role of violent mergers as the mechanism that liber- ates 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  ≈ 0 typically having significantly higher  lib values compared with those massive infallers that survive as satellites (not shown in Figure 4).The median (mean) value for  lib (for the ICL alone) for all 10.5 < log 10 ( * /M ⊙ ) < 11.5 infallers from any of the 11 clusters that did not merge with their BCG before  ≈ 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  lib for the ICL alone in Figure 4. We caution that the shape of the fit for  lib against infall stellar mass seen in Figure 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  lib fit for the ICL system alone would move upwards to more closely resemble that seen in Figure 4 for the combined BCG + ICL system.Though 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.

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,  * , of galaxies which fell into any of the 11 clusters after we began to monitor ICL assembly is shown in Figure 5.A fit to these data is also shown, which adopts the form of a Schechter (1976) function, where  is the low-mass-end slope of the function, Φ  the normalisation, and   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 ⊙ .The resulting influence of resolution effects can be noted in Figure 5 for infall stellar masses ≲ 10 9 M ⊙ .To mitigate the impact of these resolution effects on the resulting fitted function we therefore only use data from  * > 10 9 M ⊙ when fitting the Schechter function.
For consistency with the analysis described in Section 3.2.1,galaxies just falling into a cluster at  ≈ 0 or those that only "skimmed" a cluster rather than infalling are excluded from the histogram shown in Figure 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.(2)

ICL contribution as a function of infall stellar mass
The result is shown in Figure 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 Figure 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  ≈ 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 Figure 6.The median (mean) ratio between the total associated stellar mass at  ≈ 0 and the infall stellar mass of a progenitor for all ∼ 4000 infaller galaxies across all 11 clusters was ∼ 1.01 (∼ 1.18).
It can be seen in Figure 6 that the dominant contributors of stars to the ICL are galaxies with infall stellar masses between ∼ 10 10.5 M ⊙ and ∼ 10 11 M ⊙ .This is in good agreement with the general conclusion of most prior theoretical studies (e.g.Contini et al. 2014Contini et al. , 2019;;Chun et al. 2023Chun et al. , 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 Figure 6 -at log 10 ( * /M ⊙ ) = 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 Figure 5 (log 10 (  /M ⊙ ) = 11.13 +0.08 −0.07 ).Additionally, we note that the position of the peak in Figure 6 varies very little 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 ⊙ ) = 11.2 +0.1 −0.3 ), despite the two fits used for  lib as a function of infall mass being significantly different (as can be seen in Figure 4).This is a consequence of both these fits for  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 Figure 5) instead varying by several orders of magnitude over this same infall mass range.As such, the specific fit used for  lib as a function of infall mass has only a minor effect on the shape of the curve seen in Figure 6, which is instead governed almost entirely by the infall stellar mass function.
Though Figure 6 shows that the expected contribution of infalling objects with  * ≳ 10 11 M ⊙ 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 Figure 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  ≈ 0 ICL stars that were formerly part of the BCG (shown in Figure 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 Figures 4 and 6 ignores pre-processed ICL.If this analysis were repeated with the 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 Figures 4 and 6 to be substantially altered, with this extra attributed ICL mass enhancing the fractional contribution of infalling galaxies with  * ≳ 10 11 M ⊙ (as well as the  lib values determined for these galaxies), potentially shifting the location of the peak seen in Figure 6 to a slightly higher mass.
The top panel of Figure 7 shows the (normalised) cumulative contribution to the stacked, not pre-processed ICL of all 11 clusters from progenitor objects with infall stellar masses,  * , less than : 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 Figure 6.Dashed lines are once again used to indicate extrapolation beyond the range of infall stellar masses used for fitting.The corrected curve is normalised to unity, and the uncorrected curve is normalised relative to the corrected curve (adjusted by a constant factor equal to the mean ratio between total associated stellar mass at  ≈ 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 Figure 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 Figure 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 ⊙ , and ∼ 90 cent is predicted to come from progenitors with stellar masses ≳ 10 9 M ⊙ .Even after compensating for resolution effects using our extrapolated fitting function approach, that galaxies with  * < 10 9 M ⊙ contribute only ≲ 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  * ∼ 10 9 M ⊙ , further improvements in resolution will yield no significant changes in the bulk properties of the ICL and only minor changes (≲ 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 ∼ 10 7 M ⊙ 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 ≳ 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 ∼ 10 9 M ⊙ and that ≳ 50 per cent typically came from progenitors with stellar masses greater than ∼ 10 10 M ⊙ (though with significant dispersion).
It is worth noting that the shaded regions in the Figure 7 (indicating estimated uncertainties based on bootstrapping) only significantly broaden for stellar masses beyond ∼ 10 11 M ⊙ , 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 Figure 7 is somewhat coincidental: the Schechter function fitted in Figure 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 Figure 7 falls slightly below the uncorrected curve at the high mass end.
We emphasize that the analysis represented in Figure 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 Figure 7, particularly for very high mass pro-genitor objects (scarce few of which should fall into such a 10 14 M ⊙ 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.

CONCLUSIONS
Using the Horizon-AGN simulation (Dubois et al. 2014), we have studied the assembly of the intracluster light (ICL).We investigated ICL assembly in 11 simulated clusters with  ≈ 0 dark matter halo masses between ∼ 1 × 10 14 M ⊙ and ∼ 7 × 10 14 M ⊙ , tracking the stars of galaxies that fell into these clusters over the past ∼ 10 Gyr in order to quantify the differing contributions made to the  ≈ 0 ICL by galaxies with differing stellar masses on cluster infall.Our main findings can be summarised as follows: (i) On average, 44 per cent of the  ≈ 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 (Figure 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  ≈ 0 ICL stars that were pre-processed and never part of a cluster galaxy was 16 per cent (Figure 3), though in the least relaxed and most massive cluster ≳ 40 per cent of  ≈ 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 (≳ 60 per cent; Figure 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 ⊙ (Figure 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 ∼ 10 9 M ⊙ and ∼ 10 12 M ⊙ , 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 ≳ 10 9 M ⊙ , even after compensating for resolution effects (Figure 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 a ∼ 10 14 M ⊙ 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.
Though 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 Figure 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, and Golden-Marx et al. 2023).We intend to investigate this further in a future work.

Figure 1 .
Figure 1.Left: Projected positions of all the stars in the vicinity of an example  ≈ 0 Horizon-AGN galaxy cluster (ID 19 in Table1) 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 IGL of infalling groups) are shown in grey.Right: Scatter plot of the speeds, , of intergalactic stars (relative to the BCG and scaled by the 3D velocity dispersion of the cluster stars,   ) against distance from BCG centre, , (scaled by  178 ) for the same cluster and using the same colours as the middle panel.

Figure 2 .
Figure 2. The  ≈ 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 16 th and 84 th percentiles, and solid lines the maximum, minimum, and median values.
2023) we do not exclude ICL stars stripped from satellite galaxies that later merge with the BCG or otherwise do not survive until  = 0. Chun et al. (2024) reported a median value of ∼ 50 per cent for the fraction of the  = 0 ICL stars of 14 < log( 200 /M ⊙ ) < 14

Figure 3 .
Figure3.Fractions of  ≈ 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 preprocessed ICL) also shown.Dotted lines indicate 16 th and 84 th percentiles, and solid lines the maximum, minimum, and median values.

Figure 4 .
Figure 4. Mean fraction of liberated stars,  lib , as a function of progenitor galaxy infall stellar mass,  * , 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 addedSavitzky & 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.

Figure 5 .
Figure 5. Stellar masses at infall,  * , for galaxies that fall into any of the 11 clusters after  ∼ 2. A Schechter function (Equation 1) is fitted using only bins with  * > 10 9  ⊙ ; bins used for fitting are shown in orange and bins ignored for fitting are shown in red.The best-fitting parameters are presented in the upper-right and the corresponding fit line is shown in brown.The error bars and shaded regions around the fit line indicate estimated uncertainties based on bootstrapping.

Figure 6 .
Figure 6.Fractional contribution of galaxies with different stellar masses,  * , 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.

Figure 7 .
Figure7.Top: Normalised cumulative contribution to the stacked (not preprocessed) ICL of the 11 clusters by galaxies with stellar mass,  * , on cluster infall less than  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. .