© The Authors 2024

1. Introduction

Galaxy clusters are the most massive objects in the Universe. Their properties are the result of the complex interplay of multi-scale physical processes and the hierarchical assembly of numerous galaxies. A key component of these objects is the so-called intra-cluster light (ICL) coming from a diffuse distribution of stars. This stellar component was first theorised by Zwicky (1937) and was later confirmed by observations of the nearby Coma cluster (Zwicky 1951). It is likely the by-product of the numerous interactions between galaxies during the assembly and continuous growth of clusters. The study of this diffuse stellar component can therefore give us important information about the assembly history and dynamical state of galaxy clusters and the most massive galaxies in the Universe (e.g. Contini 2021; Montes 2022; Ragusa et al. 2023; Contini et al. 2023).

In recent years, there has been growing interest from the scientific community to study the ICL both from observations (e.g. Furnell et al. 2021; de Oliveira et al. 2022; Werner et al. 2023; Zhang et al. 2023) and theory (e.g. Chun et al. 2022, 2023; Marini et al. 2022; Tang et al. 2023). In observations, this is motivated by the increasing sensitivity of instruments and telescopes, as well as improvements in image-processing techniques. The availability of very deep observations, both from ground-based facilities and from space telescopes has allowed groundbreaking observations such as the detection of ICL at z > 1 using deep infrared imaging data from the Hubble Space Telescope (Joo & Jee 2023), and the first ICL study of a cluster with JWST data (Montes & Trujillo 2022). On the side of simulations, there is also increasing interest due to the availability of more computational resources and the progress in understanding and modelling of the subgrid processes responsible for the formation of stars and their feedback. Cosmological simulations like IllustrisTNG evolve a large enough volume to study the evolution of massive clusters and galaxies (Pillepich et al. 2018a,b), while at the same time resolving the ICL component, too.

In order to study the ICL, the first step is to have some way to define it, which is non-trivial and therefore one of the most important problems in the field¹. Theoretically, the ICL can be defined as stars that are not bound to any particular galaxy but to the potential of the cluster itself. However, in practice, applying this definition in observations is extremely difficult. Common approaches include surface brightness cuts (e.g. Rudick et al. 2011; Feldmeier et al. 2004; Burke et al. 2015) and fitting of functional forms such as Sérsic (e.g. Gonzalez et al. 2005; Janowiecki et al. 2010; Spavone et al. 2018) or Navarro-Frenk-White profiles (Contini & Gu 2020, 2021). For a more in-depth discussion of this matter, exploring the advantages and disadvantages of each definition, the reader is referred to the reviews by Montes (2019, 2022) and Contini (2021).

In simulations, where phase-space information is available, the task of identifying the ICL still poses some issues in separating the different gravitational potentials of galaxies and host clusters (Cañas et al. 2020). Some studies adopt the approach of considering the joint system of the brightest cluster galaxy (BCG) and the ICL surrounding it, rather than trying to separate one component from the other (Gonzalez et al. 2013; Kravtsov et al. 2018; DeMaio et al. 2018; Zhang et al. 2019; DeMaio et al. 2020). Previous efforts to separate BCG+ICL include using the velocity distribution of stellar particles to distinguish between multiple components (Remus et al. 2017), and attempts to distinguish individual particles by estimating their “boundedness” to the BCG (Cui et al. 2014). Pillepich et al. (2018b) argue that, since the physical origin of the BCG and ICL is the same, the best way to separate them is by means of a fixed spatial aperture, so that masses can be estimated and compared unambiguously across observations and simulations. Another physical way to separate the two components is given by the transition radius, that is, the radius at which the diffuse component begins to dominate the stellar mass distribution for BCG+ICL systems in groups and clusters, as recently explored by several authors (e.g. Gonzalez et al. 2021; Chen et al. 2022; Contini et al. 2022; Proctor & Lagos 2024).

One quantity that is widely used to characterise the ICL is its fractional budget, namely the stellar mass (or light) that is contained in this component out of the total stellar mass (or light) of the cluster. This has been the focus of many studies (e.g. Krick & Bernstein 2007; Burke et al. 2015; Contini et al. 2018; Montes & Trujillo 2018), but there is no consensus in the specific value for this fraction; it ranges between less than 10 and more than 50% depending on the study. While this depends significantly on the definition of ICL used (Rudick et al. 2011; Tang et al. 2018; Pillepich et al. 2018b; Kluge et al. 2021; Ellien et al. 2021; Garate-Nuñez et al. 2023), it is a clear conclusion that all massive clusters contain a non-negligible ICL component. Given the ambiguity of the mass fraction, other properties are needed to properly characterise the ICL. Some works include the study of the formation of the ICL, by analysing the age and metallicity of its stellar population (Montes & Trujillo 2018) or its colour gradients and profiles (Morishita et al. 2017; DeMaio et al. 2018), or by studying the importance of mergers in ICL formation (Murante et al. 2007; Contini et al. 2018; Kluge & Bender 2023).

Another consideration that can give insights about the formation of the ICL is the correlation between ICL mass fraction and total cluster mass (noting that these results are also dependent on how the ICL is identified). The observed lack of correlation (e.g. Montes 2022) points to a similar ICL formation efficiency across groups and clusters. There are still some tensions regarding this result, so that this lack of correlation between ICL fraction and cluster total mass cannot be confirmed (see reviews by Contini 2021 or Montes 2022), but recent results by Ragusa et al. (2023) find no significant trend for a wide halo mass range (10^12.5 − 10^15.5 M_⊙). However, for the dynamical state of the cluster, Ragusa et al. (2023) find a significant trend with ICL fraction, that has also been previously found in other studies, both theoretical (Rudick et al. 2006; Contini et al. 2014) and observational (Da Rocha et al. 2008; Montes & Trujillo 2018; Poliakov et al. 2021). More evolved and hence relaxed clusters tend to have higher ICL fractions, so that the dynamical state of a cluster can be a better indicator of ICL fraction. Contini et al. (2023) confirm this finding by showing that the main driver of ICL formation is the halo concentration (which is related to both formation time and dynamical state), in a way that more concentrated halos have higher fractions of ICL.

More recent ICL studies have been devoted to the relation between ICL and the dark matter (DM) component of clusters. By studying their density profiles in the IllustrisTNG simulations, Pillepich et al. (2018b) found that the stellar and DM halos of clusters have a similar shape, which was later confirmed by Montes & Trujillo (2018) in the Hubble Frontier Fields (HFF) clusters. Montes & Trujillo (2019) further observed that there is a tight correlation between the distribution of the diffuse stellar surface density and the surface density of total mass, which effectively means that the ICL can be used to accurately trace the DM distribution. This was theoretically confirmed by Alonso Asensio et al. (2020) using the Cluster-EAGLE simulations (Barnes et al. 2017) and later by Yoo et al. (2022) within a different set of N-body simulations. In contrast, Sampaio-Santos et al. (2021) argue that the question of whether the diffuse light faithfully traces the cluster’s radial matter distribution lacks consistent evidence yet. However, they emphasize that the diffuse light remains an excellent indicator of a cluster’s total mass. These results are particularly interesting since they provide a way of probing the DM in clusters using only deep imaging observations, and hence highlight the importance of further investigating and understanding this diffuse component of clusters.

Driven by these questions about the properties and formation and evolution of the ICL component, which still remain open, we provide a new study of ICL using cosmological simulations of galaxy clusters. In this work, we characterise the intra-cluster light of the clusters in THE THREE HUNDRED data set, that consists of regions of radius 15 h⁻¹ Mpc centred on the 324 most massive objects selected from a cosmological dark matter only simulation of side length 1 h⁻¹ Gpc. These regions have been re-simulated with two different full physics hydrodynamical codes: the smoothed particle hydrodynamics (SPH) code GADGET-X, and the novel meshless hydrodynamic and gravity solver GIZMO (Hopkins et al. 2014; Hopkins 2017). As these simulated regions share the same initial conditions in both data sets, we have the ability to perform code-to-code comparisons on the same clusters. The size of the sample and the mass range included allow for a statistical study of the ICL, analysing also its dependence on different cluster properties. This work just forms the first part in a series of papers. Here we simply start with presenting the data set and a first quantitative analysis of it at redshift z = 0. In a follow-up study (Contreras-Santos, in prep.) we put more focus on the origin of the ICL and direct predictions for observations.

This paper is organised as follows. In Sect. 2 we present THE THREE HUNDRED simulations data set, an overview of the galaxy formation models adopted by the simulations, and a description of how galaxies and dark matter halos are identified. In Sect. 3 we present our method to identify the BCG and ICL components in the simulations, while in Sect. 4 we analyse the joint BCG+ICL component of the simulated clusters and compare their stellar mass with observational estimates. In Sect. 5 we present our resulting ICL mass fraction, and how it depends on cluster total mass and dynamical state. Finally, we study the relation between the ICL and the DM component of clusters in Sect. 6, by comparing their density and velocity dispersion profiles. We summarise and discuss the results of this work in Sect. 7.

2. The Data

2.1. The Three Hundred sample

THE THREE HUNDRED project², first presented in Cui et al. (2018), is a large international collaboration focused on understanding the formation, evolution and properties of massive galaxy clusters using semi-analytic models and hydrodynamical simulations. The data set consists of 324 massive galaxy clusters with masses above ∼8 × 10¹⁴ h⁻¹ M_⊙ drawn from the parent dark matter-only simulation MultiDark Planck 2 (MDPL2; Klypin et al. 2016), which follows the hierarchical assembly of 3840³ dark matter particles in a periodic cosmological volume of 1 h⁻¹ Gpc, with a dark matter particle mass resolution of 1.5 × 10⁹h⁻¹ M_⊙, and a Plummer equivalent softening of 6.5 h⁻¹ kpc. The simulation adopts a ΛCDM cosmology with parameter values Ω_m = 0.307, Ω_b = 0.048, Ω_Λ = 0.693, h = 0.678, σ₈ = 0.823 and n_s = 0.96, in agreement with Planck 2015 results (Planck Collaboration XIII 2016). Dark matter halos in MDPL2 were identified using the ROCKSTAR halo finder (Behroozi et al. 2013), from which the 324 most massive halos were selected to constitute the THE THREE HUNDRED cluster sample.

In order to perform full hydrodynamical simulations of the selected clusters, it was necessary to use multiple levels of refinement across the entire volume of the simulation using the parallel GINNUNGAGAP2 code. This is similar to the so-called zoom technique, but instead of increasing the resolution of a specific region from the parent simulation, the region of interest was kept at the same resolution, and the remainder of the volume was “downgraded” to lower resolution to capture the large scale potential of the box while focusing the computational resources to the region of interest. The “high-resolution” region was delimited by selecting dark matter particles within a sphere of 15 h⁻¹ Mpc centred at the halo of interest at z = 0, which were then traced back to the initial conditions. These were then split into a dark and a baryonic component with the corresponding Ω_m/Ω_b ratio, resulting in a dark matter particle mass of m_p,dm = 1.27 × 10⁹h⁻¹ M_⊙, and an initial gas particle mass of m_p,gas = 2.36 × 10⁸h⁻¹ M_⊙. Each one of the selected clusters, together with the regions around them, was resimulated using GADGET-X and GIZMO-SIMBA, which are described below.

2.2. Hydrodynamical models

2.2.1. Gadget-X

GADGET-X is based on GADGET-3, a tree particle mesh (TreePM) smoothed particle hydrodynamics (SPH) code³, but it includes major improvements over the standard developer version. Hydrodynamic equations are solved using a density-entropy formulation of SPH (Springel & Hernquist 2002) using a Wendland C⁴ kernel (Wendland 1995). Gas cooling in the simulation is implemented following Wiersma et al. (2009), while the adopted star formation model is the same as the one described in Tornatore et al. (2007), and follows the star formation algorithm presented in Springel & Hernquist (2003). Each one of the created star particles represents a simple stellar population (SSP, i.e. all stars have the same metallicity) with the number density following a Chabrier (2003) initial mass function (IMF). For each of these SSPs, stellar lifetimes are considered following Padovani & Matteucci (1993), to account for the appropriate mass-loss and metal release timescale of asymptotic giant branch (AGB) stars, supernovae Ia (SN Ia) and corecollapse supernovae (SN II) individually. These events restore mass and enrich the chemical composition of surrounding gas particles with metal yields of AGB stars are taken from van den Hoek & Groenewegen (1997), for SNe Ia from Thielemann et al. (2003), and for SN II from Woosley & Weaver (1995); being the latter the only contributor to kinetic stellar feedback, implemented as in Springel & Hernquist (2003).

The adopted black hole (BH) growth and active galactic nuclei (AGN) follows Steinborn et al. (2015), which is an extended model of Springel (2005). Black holes in the simulation are implemented as collisionless sink particles that are created at the centre of friends-of-friends (FoF) objects whose mass exceeds 2.5 × 10¹⁰h⁻¹ M_⊙ and do not host a black hole, with a seed mass of M_• = 5 × 10⁶h⁻¹ M_⊙. The black holes mass accretion rate, ṁ_∙, follows a Bondi–Hoyle–Lyttleton (Hoyle & Lyttleton 1939; Bondi & Hoyle 1944; Bondi 1952) scheme, capped at the Eddington limit, multiplied by a factor of α = 100⁴. The feedback energy input into the surrounding is implemented as purely thermal.

This implementation of GADGET-X has been used to study properties of galaxy clusters such as the cool core/non-cool core cluster dichotomy (Rasia et al. 2015); to quantify the abundance and spatial distribution of neutral hydrogen (Villaescusa-Navarro et al. 2016); to study the hot gas pressure profiles, gas clumping and Sunyaev–Zeldovich (SZ) scaling relations (Planelles et al. 2017); and the galaxy cluster X-ray scaling relations (Truong et al. 2018). As part of THE THREE HUNDRED project, GADGET-X simulated clusters have been used to study galaxy and cluster population relations (Cui et al. 2018), the self-similar evolution of galaxy cluster density profiles (Mostoghiu et al. 2019) or the radial and galaxy-halo alignment of dark matter subhalos and satellite galaxies (Knebe et al. 2020) among others.

2.2.2. Gizmo-Simba

The GIZMO-SIMBA runs of THE THREE HUNDRED are performed with the GIZMO code (Hopkins 2015), with the state-of-the-art galaxy formation subgrid models following the SIMBA simulation. SIMBA, first introduced in Davé et al. (2019), is a set of cosmological hydrodynamical simulations run with a modified version of the gravity and hydrodynamics solver GIZMO (Hopkins 2015, 2017), with improved star formation and feedback subgrid models from its predecessor MUFASA (Davé et al. 2016), and including recipes for black hole growth and its corresponding feedback. Although SIMBA on its own is a separate set of simulations, the same galaxy formation model was used to run the 324 selected regions of THE THREE HUNDRED project. We therefore refer to the overall code-model implementation as GIZMO-SIMBA. Below we briefly describe the relevant details of the code and subgrid models.

GIZMO is a massively parallel gravity and hydrodynamics solver first presented in Hopkins (2015) and last updated in Hopkins (2017). The skeleton of the code is largely based upon GADGET-3 (Springel 2005), and while SPH is preserved within the code (amidst some modifications), its main feature is the implementation of so-called Meshless Finite Mass and Volume (MFM and MFV, respectively) methods to solve hydrodynamics. This can be oversimplified as a “mixture” of SPH and moving-mesh approaches that conserves mass, energy, momentum and angular momentum, that is capable of accurately capturing shocks and fluid mixing instabilities, without the need for artificial viscosity.

Photoionisation heating and radiative cooling of gas in SIMBA is computed using the GRACKLE-3.1 library (Smith et al. 2017). Star formation is modelled using an H₂-based approach where the H₂ fraction is computed following Krumholz & Gnedin (2011), accounting as well for variations in numerical resolution, as described in Davé et al. (2016). In its standard implementation, a star-forming gas particle produces a single generation of “stars”, for which its entire mass is converted into a single stellar particle. Stellar feedback accounts for the contribution from AGB stars, and from supernovae type Ia and II. AGB stars and Type Ia SNe produce galactic winds using Davé et al. (2016) two-phase model with a fraction of ejected hot gas to be 30%.

Black holes in the simulations are seeded in FoF galaxies that reach stellar mass M_* ⪆ 10^10.5 M_⊙ and do not already contain a BH particle by converting the star particle closest to the centre of mass into a BH particle. Black hole particles are repositioned to the minimum of the gravitational potential of their host FoF galaxy if it is located within 4 R₀, where R₀ is the size of the black hole kernel used for accretion. Black hole particles are allowed to merge if two of them are located within R₀ and their relative velocity is smaller than three times their mutual escape velocity. SIMBA implements a novel black hole mass accretion model using a two-mode approach for hot and cold gas accretion. In this model, the accretion of hot gas follows the standard Bondi–Hoyle–Lyttleton (Bondi 1952) spherical accretion, while the accretion of cold gas is modelled using the torque-limited approach from Hopkins & Quataert (2011) and its numerical implementation from Anglés-Alcázar et al. (2017).

SIMBA accounts for two different feedback mechanisms associated to black holes: kinetic outflows, which arise from a “radiative” mode at low accretion rates (ṁ_BH) and from “jet” modes at high accretion rates; and feedback from X-ray sources. The addition of feedback from accretion disk X-rays is a novel feature of SIMBA and, although it does not contribute largely to the overall feedback energy release, is key to quenching the most massive objects (Davé et al. 2019). Note that the GIZMO-SIMBA run of THE THREE HUNDRED clusters uses recalibrated parameters for the baryon model due to the different simulation resolutions (see Cui et al. 2022, for details).

2.3. Halo catalogues

We use Amiga Halo Finder (AHF⁵; Gill et al. 2004; Knollmann & Knebe 2009) to identify dark matter halos in the simulations. We use the same halo catalogues as previous studies from THE THREE HUNDRED collaboration for consistency and to take advantage of the estimated properties and derived formation histories of the halos in other studies (e.g. Cui et al. 2018; Mostoghiu et al. 2019). Halos are identified by locating density peaks in the simulation volume in a similar fashion as adaptive mesh refinement (AMR) techniques, in which space is iteratively divided when the total matter density within a cell exceeds a given threshold to constrain the location of overdense regions. Once these density peaks are found, particles of all species (dark matter, stars, gas and black holes) are grouped into gravitationally bound halos. The halos are then spherical regions composed by at least 20 particles with a total matter density of Δ times the redshift-dependent critical density of the Universe, ρ_crit, with a mass of $M_{\text{halo}}=\mathrm{\Delta }\cdot {\rho }_{\text{crit}}(z)\cdot R_{\text{halo}}^3\cdot 4\pi /3$, where R_halo is the corresponding spherical overdensity radius. Two catalogues have been produced with enclosing overdensities of Δ = {200, 500}; physical quantities corresponding to these will be denoted by the subscript 200 and 500, respectively, whenever necessary. Regarding the substructures, AHF defines subhalos as halos that lie within R_halo of a more massive halo, which is called the host halo. The mass of this host halo then includes the masses of all the subhalos contained within it.

3. Identification of BCG and ICL

In general, the BCG is defined to be the brightest galaxy in the cluster, while the ICL is comprised by the stars that are not bound to any galaxy (central or satellite) but to the potential of the cluster itself. However, these definitions are not trivial to implement in observations of galaxy clusters, and so identifying the different components is a challenging task. This way, separating the outer regions of the BCG from the ICL is still an ill-defined problem (Gonzalez et al. 2005; Krick & Bernstein 2007; Jiménez-Teja & Dupke 2016). Moreover, the low surface brightness of the ICL (Mihos et al. 2005; Zibetti et al. 2005; Rudick et al. 2010) and its contamination by foreground and background galaxies make it even more complicated. In simulations, this problem is simplified by the availability of full kinematic information of the baryonic and dark matter particles. However, it is still not trivial to separate the gravitational potential of the host cluster from that of the individual galaxies – and thus the BCG and the ICL. For this reason, different works can follow different approaches, and it is hence important to clarify how this is done. In this section we explain how we define and find the BCG and the ICL in the galaxy clusters of THE THREE HUNDRED simulations.

3.1. Defining the BCG

Since the edge of a BCG is not clearly delimited, we determine it by selecting all the stars inside a sphere of certain radius, using the halo centre as the origin. This way of defining the BCG, based on a fixed physical radius, has often been adopted in simulations. For instance, McCarthy et al. (2010) used a radius of 30 h⁻¹ kpc, whereas in Ragone-Figueroa et al. (2018) three different radii were used too (30 kpc, 50 kpc and 0.1 R₅₀₀). Pillepich et al. (2018b) further consider that this rather arbitrary separation between BCG and ICL is validated by theoretical arguments, since both components are formed by smooth stellar mass accretion and mergers, although in different relative amounts. They use a separation of 30 kpc, and consider also the effect of changing it to 100 kpc. From a more observational approach, Kravtsov et al. (2018) advocate the use of stellar masses defined this way for comparison with observations, and they provide values for nine BCGs using several different radii. Besides, Stott et al. (2010) showed that, using a 50 kpc aperture radius, BCG luminosities can be recovered with less than 5% difference from those obtained in some observational analyses.

Hence, using a fixed aperture allows us to easily identify in the simulation all the stellar particles that are in this central region. As we have seen, different values for this aperture are used throughout the literature. In a previous work with the same simulations, we obtained the BCG for three different fixed apertures: 30, 50 and 70 kpc (Contreras-Santos et al. 2022). In this work, for clarity and simplicity, we will focus on the 50 kpc aperture and, unless specified otherwise, show the results only for this aperture. Please refer to Sect. 3.3 and the Appendix A for results based upon the other apertures.

Despite our choice, there are some caveats of this BCG definition that need to be mentioned. First, we adopt a 3D cut, this is, we consider the stars outside the central sphere on the foreground and background to belong to the ICL, which is not fully reproducible in observations. In this regard, Cui et al. (2022) verified that using the projected 2D cut increases the BCG masses by less than 20%. Additionally, this definition can introduce some bias if applied to a sample with a wider mass range, where the sizes of the central galaxies can change significantly from group to cluster scales. Alternative methods that look for the transition radius from BCG to ICL (e.g. Contini et al. 2022; Chen et al. 2022) can be valuable tools in overcoming these challenges, while offering a better physical understanding of the situation. However, for our study, limited to high masses, a fixed aperture has proven adequate to identify the BCG without introducing any significant bias.

3.2. Defining the ICL

Once the BCG is defined and the particles that assemble it are selected, we can easily find the ICL in our simulations. In order to do this, we select all the stellar particles that are bound to the central cluster, but do neither belong to the BCG nor to any subhalo. This separation of particles can be easily done using the AHF catalogues, which give the IDs of all the particles that belong to each of the halos and subhalos, for the simulations run with both GADGET-X and GIZMO-SIMBA codes. In other words, from all the stellar particles in the halo, we first remove the BCG (whose definition is solely based on a fixed aperture of 50 kpc, not dependent on the halo finder) and then we remove the particles bound to substructures (as indicated by AHF). The remaining particles are the ICL of the cluster. Note that in special situations, such as major mergers, a large secondary halo can bring its own ICL. In this case, as long as the secondary halo is identified as a different object from the main one, all its particles are removed for the counting of the main halo ICL.

In Fig. 1, for one cluster selected as example, we show a 2D projection of all the stellar particles along one axis of our simulations for, from left to right, total stellar mass within R₅₀₀, the ICL and the BCG alone. Top row is for GADGET-X and bottom row for GIZMO-SIMBA. The same colour intensity is used for all the visualizations, with darker tones showing higher densities. In the second column, which corresponds to the ICL, we can see how the substructures in the cluster, that were present in the first column, have been removed, leaving the ICL as a diffuse component. Since we use a 3D aperture to extract the BCG, we still see particles in the innermost 50 kpc when projecting them along one axis, but they are in fact outside the central sphere, either in the front or in the background. We want to note that particles in the outskirts of galaxies might not be identified to belong to the galaxy, and thus still be part of the ICL. Due to projection effects, these would appear in Fig. 1 as small clumps, but would be in fact empty shells, not representing a significant part of the ICL of the cluster. The last column in Fig. 1 shows a zoom-in to the BCG, which, as described before, is selected as all the stellar particles within 50 kpc of the halo centre. The red circle in that panel depicts again a sphere of radius 100 kpc. We note the BCG’s stellar particle distribution is concentrated towards the centre, showing a great reduction towards the outer parts. We are therefore confident that the BCG has properly been excised and separated from the ICL (but see also the discussion in Sect. 3.3).

Fig. 1. For one of the 324 clusters (cluster 6), 2D projection of all the stellar particles along one axis that belong to, from left to right, the total stellar mass within R500, the ICL and the BCG alone. The same colour intensity is used for all the visualizations with darker tones showing higher densities; concentric circles have radii of 100 kpc (inner), 500 kpc (middle), and 1 Mpc (outer). For the rightmost plots only the inner circles (100 kpc radii) are shown.

Finally, comparing the top and bottom rows in Fig. 1, we can see that in general the same structures are generated by both codes, with differences being found only on smaller scales. We can also see that in GADGET-X there are more particles, making the maps look smoother, while by the colour intensity we can see that the mass is very similar in both simulations. This points to different stellar particle resolution between the two simulations, due to their different implementations of star formation. We will discuss this in more detail in the following section.

3.3. Dependence on aperture size

Since our BCG definition involves a “free” parameter, namely the size of the fixed aperture, the resulting ICL also depends on the value used for it. Although we have chosen a value of 50 kpc for this parameter, before continuing with our results we believe it is relevant to show how sensitive the resulting BCG and ICL are to the size of the aperture.

In Fig. 2 we compare the masses of the two components, as computed using different apertures. The upper 2 × 2 subplot compares the BCG mass, the left column (in green) being the results for GADGET-X and the right one (in magenta) for GIZMO-SIMBA. The top row shows a histogram of the ratio M_BCG, 50/M_BCG, 30, where we can see that the distribution is wide, reaching even more than 2, which means that the BCG mass is being doubled when changing the aperture from 30 to 50 kpc. In the second row the ratio shown is M_BCG, 70/M_BCG, 50, and in this case we see that the values are significantly smaller, the median growth being around 30%. This indicates that the difference between the BCGs computed for 70 and 50 kpc apertures is not as relevant as in the first row, which suggests that increasing (moderately) the aperture beyond 70 kpc does not make such a significant difference in the BCG mass.

Fig. 2. Dependence of BCG and ICL mass on aperture size. Top: distribution of the ratio between the BCG mass of the THE THREE HUNDRED simulations obtained for apertures: 30 and 50 kpc (first row), and 50 and 70 kpc (second row). The left column (in green) is for GADGET-X, while the right one (in magenta) is for GIZMO-SIMBA. Bottom: same for ICL mass. Since ICL mass decreases with increased aperture, the ratio is inverted now to keep the values above 1, but the interpretation remains the same.

The lower 2 × 2 subplot in Fig. 2 shows this same comparison but for the ICL mass. Since the BCG mass increases for increased aperture, the ICL mass is reduced. However, to keep the histogram values above 1, we now plot the inverse ratios. The meaning of the plots remains the same, the further the values are from 1 the more the mass is changing from one aperture to the other. We see in this case that the values are much closer to 1 than for the BCG masses; since the ICL extends quite far from the halo centre, the choice of a fixed aperture between 30 and 70 kpc does not affect that much the final mass of this component.

In conclusion, Fig. 2 shows that, while for the BCG the size of the aperture selected can have an important effect on the resulting BCG mass, the ICL mass is not very much dependent on this value. Since the main focus of the paper is on the ICL component of clusters, we believe these differences should not be discussed in the main text of the paper. However, we include an appendix at the end of this work where we show how the main results of the paper, derived in the following sections, are affected by the size of the aperture, for which we use 30, 50 and 70 kpc. From now on in this paper – and unless specified otherwise – we will use a fixed aperture of 50 kpc to define the BCG. For the results for 30 and 70 kpc apertures, we refer the reader to Appendix A. We note that, although both the BCG and ICL masses are affected by the aperture size, by definition the joint component BGC+ICL is independent of it. We will start the next section by studying this component and comparing our findings to those of previous observations.

4. BCG+ICL

In this section we investigate the stellar mass content and distribution of the BCG+ICL, M_{*BCG + ICL}, as a function of cluster halo mass M₅₀₀. We start with comparing against observations for which the BCG+ICL is available to understand similarities and differences intrinsic to the simulations themselves and the implications these may have on other results.

In Fig. 3 we show the stellar mass of the BCG+ICL as a function of halo mass, M₅₀₀. The first row shows the BCG+ICL total stellar mass, the second row shows the total mass inside an aperture of 100 kpc and the bottom row shows the mass content between 10 and 100 kpc. Individual estimates from the simulated clusters are shown as squares, and observational estimates are shown as grey symbols in each panel, as labelled; we note that although multiple groups and clusters may be resolved within the re-simulated regions, these plots only show the main 324 cluster sample. Regarding the total BCG+ICL mass (top row), both simulations are able to reproduce observational estimates from Gonzalez et al. (2013) and Kravtsov et al. (2018) at M₅₀₀ < 10^14.5 M_⊙. However, at larger masses both simulations overpredict observational estimates and best fit line of Kravtsov et al. (2018), with GADGET-X showing a larger spread than GIZMO-SIMBA on the M_{*BCG + ICL} reaching ∼1 dex at M₅₀₀ ≃ 10¹⁵ M_⊙. Note that the total BCG+ICL mass in Kravtsov et al. (2018) corresponds to the mass within R_out, defined as the radius where the error in fitted Sérsic profiles reaches 33% due to background uncertainty. Furthermore, since observations are limited by the sensitivity of the instruments to detect the faint intra-cluster light and are prone to projection effects, our comparison is not fully consistent. Using 3D apertures that include all ICL particles, some excess mass in the simulations can be expected.

Fig. 3. GADGET-X (left) and GIZMO-SIMBA (right) BCG+ICL stellar mass within different spherical regions as a function of total halo mass. First row shows the total BCG+ICL stellar mass, squares denote individual measuremente of each of the 324 main clusters, grey symbols show data from Gonzalez et al. (2013) and Kravtsov et al. (2018), as labelled; the shaded region shows the best-fit from Kravtsov et al. (2018). Second row shows the total BCG+ICL mass within 100 kpc aperture and the third one shows the BCG+ICL stellar mass contained between 10 and 100 kpc apertures. Grey symbols show observational estimations from DeMaio et al. (2020).

The BCG+ICL mass within a 100 kpc aperture (second row in Fig. 3) shows different mass contents for GADGET-X and GIZMO-SIMBA. This difference is only seen in the second row, and not in the third one, which shows the mass between 10 and 100 kpc apertures. This means that the difference between the two simulations is present in the innermost parts of the BCG+ICL ensemble. Compared to observational estimates from DeMaio et al. (2020), for apertures of 100 kpc (second row), GIZMO-SIMBA shows a notable agreement with observations, while GADGET-X shows bigger scatter, with up to ∼2 dex spread. Lastly, regarding the mass content in between 10 and 100 kpc, both simulations show a good agreement with DeMaio et al. (2020) within the halo mass range where observations and simulations overlap. Although the trend is well followed at M₅₀₀ > 10¹⁵ M_⊙, this does not seem to hold for masses below the range contained within the simulations. From this last row we can infer that the observed spread in the 100 kpc aperture estimate comes from the stellar mass contained within the innermost 10 kpc.

We further investigate the BCG+ICL component in Fig. 4, which shows the volumetric-density profile of this component, scaled by cluster radius, R₅₀₀, stacked for the 324 clusters in THE THREE HUNDRED sample. The solid lines show the median density value at a given bin, while the shaded regions indicate the 16th and 84th percentiles, with the different simulations coloured as labelled. For comparison, we also include the density profile of the satellite galaxies (i.e. including only stellar mass) as dotted lines. The inset in Fig. 4 shows only the BCG+ICL component, in logarithmic scale in the X-axis, zooming in for the innermost region of the profile, up to 0.1 R₅₀₀ (showing the location of 10 and 100 kpc as vertical regions for reference).

Fig. 4. Density profile of the BCG+ICL component (solid) and satellite galaxies (dotted) for the 324 main clusters of GADGET-X (green) and GIZMO-SIMBA (magenta) scaled by cluster spherical overdensity radius, R500. Lines shows the median of each bin, and shaded regions delimit the 16th and 84th percentile of each bin. Inset shows only the BCG+ICL component in logarithmic scale in the X-axis, to focus on the inner part of the profile. Here, the vertical shaded regions show for reference the 16th and 84th percentile of the location of 10 and 100 kpc, as labelled.

The profiles show some important features. First, the profiles of both samples display very similar amplitudes and slopes between 0.01 and 0.8 R₅₀₀ with a 1σ scatter in the density of 0.3 dex, suggesting that the stellar distributions are self-similar at least within the mass range of the sample. This not only confirms the agreement in the mass contained in between 10 and 100 kpc between the samples, but also suggests that we should expect similar results for the annulus up to ∼0.8 R₅₀₀, after which the value of the density becomes significantly low, showing the “edge” of the BCG+ICL component. The remarkable agreement of M_{*BCG + ICL} within 10 < R < 100 kpc (see inset) demonstrates that the overall physical processes responsible for the formation of the BCG and the ICL around it are well modelled by both simulations. Regarding the comparison with satellite galaxies, we can see that, in median, satellites start to dominate over the diffuse component at ∼0.35 R₅₀₀, this being slightly closer to the centre for GIZMO-SIMBA. However, the satellites’ profiles show a wide scatter, so that this threshold can range from ∼0.2 to 0.6 R₅₀₀ for the different clusters.

Focusing now on the inset in Fig. 4, we can see that the density profiles of the two simulations’ data sets differ the most within the innermost radii, R < 0.01 R₅₀₀, where GADGET-X clusters show an upturn in the density profile, displaying a slope of ∼ − 2, while GIZMO-SIMBA start to flatten towards the centre, with a slope of –0.9. The density difference between the two data sets’ medians within this region is ≳0.5 dex, showing that GADGET-X produces denser and more massive centres than GIZMO-SIMBA, which explains the difference observed in the 100 kpc aperture measurements of Fig. 3, where GADGET-X displays a large scatter towards higher stellar masses at fixed halo mass. Moreover, the large spread observed for both simulations in M_{*BCG + ICL} within 10 kpc apertures is tightly connected to the numerical resolution of the simulations and their specific implementation of the star formation. On one side, GADGET-X and GIZMO-SIMBA were run with softening lengths of 6.5 h⁻¹ kpc and 5 h⁻¹ kpc, respectively, which means that the kinematics within 10 kpc apertures are very close to the scales that can be considered as well resolved. Some of the differences between the two simulations are likely coming from the different stellar particle resolutions, recalling that the implementation of GADGET-X allows gas elements to produce multiple generations of stars, while for GIZMO-SIMBA the entire mass of a gas particle is converted into a single star. This can make the stars produced by the two simulations to differ significantly in mass, although the initial dark matter and gas mass resolution are the same between both simulations. A direct consequence of distributing the same mass amongst a fewer number of particles is that the system becomes more collisional and therefore the dynamics in the simulation are more prone to be affected by numerical scattering (e.g. Binney & Knebe 2002; Ludlow et al. 2019), as the mass ratio between dark matter particles and star particles is ∼10 for GIZMO-SIMBA and ∼40 for GADGET-X.

5. ICL mass fraction

In the previous sections we focused on explaining how we define both the BCG and ICL of clusters and presenting results for the joint component BCG+ICL. In this section our aim is to characterise the ICL in THE THREE HUNDRED clusters, investigating how much of the total stellar mass of the cluster can be found in this component and how this ratio might depend on different cluster properties.

5.1. Mass fraction of ICL

We first compute the mass fraction of the ICL within the cluster, dividing the stellar mass found in the ICL as defined in Sect. 3 by the total stellar mass within the considered overdensity (M_*, 200 or M_*, 500). In Fig. 5 we show this fraction as a function of the cluster’s halo mass, in the top panel using M₅₀₀ and in the bottom M₂₀₀. The dots show the individual values for all 324 regions, while the lines (solid for GADGET-X and dash-dotted for GIZMO-SIMBA) indicate the median values for the different mass bins. The corresponding shaded regions depict the 16th–84th percentiles, in green for GADGET-X and magenta for GIZMO-SIMBA.

Fig. 5. ICL mass fraction as a function of cluster mass, considering overdensities of 500 (top) and 200 (bottom). The lines show the median values and the shaded regions the 16th–84th percentiles. Results are shown for both GADGET-X (green solid line) and GIZMO-SIMBA (magenta dash-dotted line). Note that the scale is the same for both plots in the Y-axis, but different for the X-axis.

The main feature that stands out in Fig. 5 is that the ICL mass fraction shows no apparent dependence on cluster mass for our considered mass range. The specific value is different between R₂₀₀ and R₅₀₀, with the fraction being higher for the innermost overdensity of R₅₀₀. In this case it has an approximately constant value of f_ICL ∼ 0.36, while for R₂₀₀ the value is around 0.29. Another important result to highlight here is that both hydrodynamical simulations GADGET-X and GIZMO-SIMBA show the same results, with no significant difference between them. This is in agreement with recent studies of ICL at low redshifts, that have also shown no dependency of the ICL mass fraction with cluster mass (Montes 2022; Ragusa et al. 2023), for our considered mass range. Theoretical studies such as Contini et al. (2023) or Proctor & Lagos (2024) further support this scenario, with no discernible correlation for our cluster mass range.

We also show in Fig. 6 the stellar mass fraction in the other two components: the BCG and the rest of the galaxies in the cluster, now focusing only on the overdensity 500. In this case we show the median value for the 324 clusters, together with the 16th–84th percentiles. The green bars are for GADGET-X, while the magenta ones are for GIZMO-SIMBA. We see that the results are similar for the two simulations, with the difference that there is higher scatter in GADGET-X, especially for the BCG, where the median is similar for both but the fraction in GADGET-X can be as high as ∼20%, while it does not reach 15% in GIZMO-SIMBA. This is the same effect we saw in the inner region of the density profiles in Fig. 4, due to the differences in numerical resolution and implementation of star formation between the two codes. Regarding the ICL, we see that the fraction is between 30 and 50% in both cases (this result is the same as Fig. 5, without binning by cluster mass). Finally, for the rest of the galaxies we also see a higher scatter in GADGET-X, with the fraction ranging from 36 to 64%; while for GIZMO-SIMBA the range is between 43 and 67%. For overdensity 200 we are not showing the results here, but the conclusions are very similar, with a smaller ICL fraction (as we could see in Fig. 5), and hence a larger contribution from the satellite galaxies.

Fig. 6. Stellar mass fraction for the BCG, ICL and rest of stars. The bar shows the median value for the 324 clusters, while the error bars indicate 16th–84th percentiles. The fraction is computed dividing the stellar mass content in each of the components by the total stellar mass inside R500, M*, 500.

In general, we can conclude from Fig. 6 that the BCG mass fraction is below 15 and 20% for GIZMO-SIMBA and GADGET-X, respectively. Compared with the ICL fraction, which is between 30 and 50%, this means that, in terms of mass, the BCG+ICL component is dominated by the ICL when considering the whole range up to R₅₀₀ (or R₂₀₀). The whole cluster stellar content, however, is in general dominated by the satellite galaxies, which can account for more than half the stellar mass of the cluster. This is a high percentage when compared to the situation for the total mass of clusters, that is including not only stars but also DM and gas. In this case, the substructures account for a median of only 10% for our THE THREE HUNDRED sample.

How do our results fare against other studies previously presented in the literature, either observationally or theoretically? The ICL mass fraction is a property that has been investigated by several studies, but it shows a great dispersion in the reported values. For 6 clusters at redshift 0.3 < z < 0.6 observed by the Hubble Frontier Fields (HFF) survey, Montes & Trujillo (2018) measured an ICL fraction of ∼7% for R < R₅₀₀, in agreement with previous measurements of the same objects by Morishita et al. (2017). Around redshift 0.2–0.3, Burke et al. (2015) found the ICL fraction to be around 20% in the CLASH survey data, while Krick & Bernstein (2007) estimated the fraction to be 6–20% in 10 nearby clusters. Using DES observations, Zhang et al. (2019) found the contribution of BCG+ICL to be 44 ± 17% of the total cluster stellar luminosity within 1 Mpc, for a sample of 300 clusters at redshift in the range 0.2–0.3. For lower redshifts, Gonzalez et al. (2007, 2013) estimated the BCG+ICL component to make up 33% of the total cluster stellar mass within R₂₀₀, or 25–55% within R₅₀₀. In the Local Universe (z < 0.05), Ragusa et al. (2023) found that the ICL fraction ranges from 5 to 50% using VEGAS images. On the theoretical side, Rudick et al. (2011) used a set of numerical N-body simulations to find that the fraction of the total cluster’s luminosity that is in the ICL ranges from 9 to 36% at z = 0, with a strong dependence on the ICL definition used. The semianalytical models by Contini et al. (2014, 2018) predict an ICL fraction that varies between 20 and 40% depending on the particular implementation adopted. In a detailed study of the IllustrisTNG simulations, Pillepich et al. (2018b) found the BCG+ICL fraction in the total cluster stellar content to be around 50%, for objects with mass above 10¹⁴ M_⊙. Using the Cluster-EAGLE simulations, Alonso Asensio et al. (2020) report a fraction f_ICL = 0.091 ± 0.013 within R₂₀₀. Our result, although in the upper part of the reported range, is still within this range, and thus in overall agreement with previous ICL measurements.

5.2. Influence of dynamical state

In Fig. 5 we studied the ICL mass fraction as a function of cluster mass, and found no significant correlation within the given cluster mass range. We will now investigate the effect of the dynamical state of the clusters. The dynamical state of clusters is a property of special interest, since it can directly affect the cluster mass estimation, and is linked with other properties like halo concentration (Neto et al. 2007). Previous works with THE THREE HUNDRED simulations have been devoted to studying the connection between morphology and dynamical state (De Luca et al. 2021; Capalbo et al. 2021) or its influence in the hydrostatic mass bias (Gianfagna et al. 2023). The dynamical state is also strongly correlated with the mass accretion history of clusters (Mostoghiu et al. 2019; Haggar et al. 2020), and hence its connection with the ICL can provide interesting information regarding the formation channels of the ICL and the cluster itself. Using different sets of simulated galaxy clusters, works like Rudick et al. (2011), Cui et al. (2014) and Contini et al. (2023) found more dynamically evolved clusters to have a higher ICL fraction. The opposite interpretation is given in observational studies by Jiménez-Teja et al. (2018, 2023), claiming that higher fractions are expected in merging systems.

There are different ways to quantify the dynamical state of clusters, depending on the nature and purpose of each study. In simulations, Neto et al. (2007) introduced three different parameters to be used as proxies for relaxation of clusters. These parameters are the centre of mass offset, Δ_r, which is the offset between the positions of the centre of mass of the cluster and the density peak, normalised to the halo radius; the subhalo mass fraction, f_s, which is the fraction of cluster mass contained in subhalos; and the virial ratio η, defined as η = (2T − E_s)/|W|, where T is the total kinetic energy of the cluster, E_s its energy from surface pressure and W its total potential energy (see Cui et al. 2017, for an updated calculation for hydrodynamic simulations). Cui et al. (2018) used these parameters to study the relaxation of the clusters in THE THREE HUNDRED. In a later work, De Luca et al. (2021) thoroughly investigated the dynamical state of these same clusters, comparing theoretical indicators (Δ_r and f_s) with results from morphology classifications. Also with THE THREE HUNDRED clusters, Haggar et al. (2020) introduced the relaxation parameter χ_DS, that combines the three indicated parameters into a single measure of dynamical state, and used it to correlate with the fraction of backsplash galaxies.

Focusing individually on these three theoretical indicators, we present in Fig. 7 the ICL mass fraction of each of our 324 clusters as a function of the indicated parameter. From left to right, the first three panels are for subhalo mass fraction f_s, centre of mass offset Δ_r and virial ratio. The first two parameters have to be minimised for the cluster to be most relaxed, while the virial ratio is 1 for a perfectly relaxed cluster (and hence we plot |1 − η|). The vertical dashed lines indicate the thresholds used in Cui et al. (2018) to distinguish between “relaxed” and “unrelaxed” clusters: 0.1, 0.04 and 0.15, respectively. Clusters with values below the threshold are considered to be relaxed and clusters to the right of the vertical line are labelled as unrelaxed. The value of the Spearman correlation coefficient for each of the plots is indicated in the lower left corner.

Fig. 7. ICL mass fraction within R200 of the cluster as a function of the three different theoretical parameters introduced by Neto et al. (2007) to quantify the dynamical state (from left to right subhalo mass fraction fs, centre of mass offset Δr and virial ratio η), together with the formation time, zform. The first three panels are coloured by zform, while the fourth one is coloured by fs, as indicated in the colourbars. In all panels, the Spearman correlation coefficient, rS, is indicated in a lower corner. The vertical dashed lines indicate the thresholds used in Cui et al. (2018) to distinguish between “relaxed” and “unrelaxed” clusters: 0.1, 0.04 and 0.15, respectively. Clusters to the left of these lines (value lower than the threshold) are relaxed, and those to the right are unrelaxed.

We can see in Fig. 7 that there is a very clear negative correlation for the first two parameters, meaning that the ICL mass fraction is lower for less relaxed clusters. The correlation coefficient is the highest for the subhalo mass fraction. This could be expected given that the definition of the ICL itself is based on removing the subhalo masses (although for the ICL we consider only stellar mass while f_s takes also DM into account). If a cluster has more mass in subhalos, it will also have more mass in satellites. So, if some of the satellites have not been subject to stripping because they are “shielded” by their subhalos, it is likely that the ICL is lower than in a halo with a lower f_s. Given this strong correlation, we fit the points to a straight line and show the results in the first panel of Fig. 7, together with the parameters of the best fit, which shows a slope of −0.36 (note that f_s is in logarithmic scale in this plot). For the virial ratio we can see that the correlation in the plot is not as clear and the coefficient decreases to –0.35. The same trend remains nevertheless, with the clusters farther from virial equilibrium having in general lower ICL mass fractions.

We also include in the last panel of Fig. 7 the correlation with the formation time of the cluster, z_form, which, as defined in Mostoghiu et al. (2019), is computed as the redshift at which M₂₀₀ of the halo is equal to half its value at z = 0, that is M₂₀₀(z_form)/M₂₀₀(z = 0) = 0.5. In this case we see a very clear positive correlation, with r_S ∼ 0.8, indicating that the earlier the cluster was formed, the higher its ICL mass fraction. Clusters that have accreted much of their mass more recently seem to be in turn dominated by satellite galaxies, which can be indicative of events such as major mergers. These results are in agreement with Haggar et al. (2020), where a clear correlation is found between formation time and dynamical state of clusters, and with Contini et al. (2024), who find a connection between the ICL and the formation time of the cluster itself. These correlations can be further seen in the colour coding of the points in Fig. 7. The first three panels are coloured by z_form, while the fourth one is coloured by f_s, with the four plots showing clearly that earlier formed clusters (higher z_form) are in general more relaxed. To avoid repetition, we only show here the plots for GADGET-X, but the same general trends hold for GIZMO-SIMBA. Moreover, although for simplicity it is not explicitly shown here, we also find a positive correlation between the halo concentration (computed by assuming an NFW profile, Navarro et al. 1997) and the ICL fraction, very similarly to the results by Contini et al. (2023). This reinforces our findings, given that clusters formed earlier are more concentrated (e.g. Mostoghiu et al. 2019).

To further explore the physical situation behind Fig. 7 and quantify the dynamical state of clusters in a different way, we now use the information from a previous work about galaxy cluster mergers in THE THREE HUNDRED simulations (Contreras-Santos et al. 2022). In this work, we defined mergers as significant mass increases happening in less than half a dynamical time of the cluster. We apply now this same procedure and consider mergers as increases of at least 25% in the total mass of the cluster⁶, that is (M_f − M_i)/M_i ≥ 0.25, where M_f and M_i are the masses at two different snapshots separated by half a dynamical time. In the top panel of Fig. 8 we plot the ICL fraction against the time since the last merger undergone by the cluster. Although the scatter is high, we find a positive correlation with Spearman coefficient r_S ∼ 0.6, so that clusters that have recently undergone a merger have a lower ICL fraction, in agreement with our previous findings. The lower panel in Fig. 8 shows also the ICL fraction, but now as a function of the number of mergers undergone by the cluster since z ∼ 1.5. The clusters are binned for each value between 0 and 6 mergers, the dots showing the median value and the errors computed as 16th–84th percentiles. The bar plot in the bottom shows the number of clusters within each bin. We see here that, although the scatter is still high and the correlation is weaker, there is also a trend for clusters that have recently experienced more mergers to have a smaller ICL fraction.

Fig. 8. Correlation between ICL mass fraction and merging history of clusters. Top: ICL mass fraction within R200 of the cluster as a function of the time since the most recent merger. The mergers are defined following Contreras-Santos et al. (2022), as mass increases of at least 25% of the initial mass in less than half a dynamical time. Bottom: ICL mass fraction as a function of the number of mergers experienced by the cluster since z ∼ 1.5. The clusters are binned by number of mergers experienced, with the bar plot below indicating the number of clusters in each bin. The dots correspond to the median values, while the errors are 16th–84th percentiles.

Overall, our results point to a scenario where more relaxed clusters have a higher ICL fraction than disturbed or merging clusters, that have formed more recently. Contini et al. (2023, see also Contini et al. 2024 for the extension of these results to z = 2) propose a very similar scenario, where the main driver of ICL formation is the concentration, that also separates dynamically evolved clusters (more concentrated) from younger clusters (less concentrated). Our results are also in agreement with observational findings by Ragusa et al. (2023, see also Da Rocha et al. 2008; Poliakov et al. 2021; Ragusa et al. 2022), who find that clusters with a higher fraction of early type galaxies (more relaxed) show higher ICL fractions. Nevertheless, recent results by Jiménez-Teja et al. (2023, see also Jiménez-Teja et al. 2018; de Oliveira et al. 2022) support the opposite situation, where merging clusters present a higher ICL fraction than more passive and relaxed ones. This highlights that the questions of the origin and formation of the ICL and its connection to the formation of the cluster are still open and remain to be understood. For this paper we prefer to keep our study to z = 0, while we will go into much detail about where the ICL stars are coming from in a future paper already in preparation.

5.3. Selecting the most relaxed and unrelaxed clusters

We have seen that, while the cluster mass has no influence on the ICL mass fraction, this fraction shows a clear correlation with the dynamical state of the cluster, such that more relaxed clusters have a higher ICL fraction. We now want to further quantify this effect. For this we will select the most relaxed and most unrelaxed clusters in the sample, and study how their stellar mass is distributed among the different components. Following Haggar et al. (2020), in order to obtain a continuous, non-binary measure of dynamical state, we combine the three theoretical parameters previously mentioned into a single measure of dynamical state, the so-called “relaxation parameter” of a cluster:

$$\begin{array}{c}\hfill {\chi }_{\mathrm{DS}}=\sqrt{\frac{3}{{\left(\frac{{\mathrm{\Delta }}_r}{0.04}\right)}^2+{\left(\frac{f_s}{0.1}\right)}^2+{\left(\frac{|1-\eta |}{0.15}\right)}^2}}·\end{array}$$(1)

This parameter satisfies that, for “dynamically relaxed” clusters, χ_DS ≳ 1, and the higher χ_DS the more relaxed a cluster is.

Using this parameter we can easily select the most relaxed clusters within our sample, together with the most unrelaxed or disturbed. Given the size of THE THREE HUNDRED sample, that allows us to make subsamples that are still statistically significant, we will select the 50 most extreme clusters in both sides. In Fig. 9 we show, similarly to Fig. 6, the stellar mass fraction in the different components: BCG, ICL and the rest of the galaxies in the clusters. But in this plot it is shown only for the subsamples of the 50 most relaxed and disturbed clusters, that is those with, respectively, highest and lowest χ_DS. The bars indicate the median values, while the error bars show the 16th-84th percentiles. The different colours designate the different codes, green for GADGET-X and magenta for GIZMO-SIMBA. Figure 9 shows very clearly the effect we have previously discussed, with relaxed clusters having a higher ICL fraction. It also highlights how the stellar mass content in satellite galaxies is significantly higher for disturbed clusters. For the BCGs, the error bars show that the difference between relaxed and disturbed clusters is less significant, but we still see a change from ∼12% in relaxed clusters to ∼7% for disturbed ones (almost 50% difference between them). From a different perspective, Fig. 9 shows how the ICL can be used as an indicator of the cluster dynamical state: relaxed clusters can be identified as those with f_ICL > f_rest.

Fig. 9. Stellar mass fraction for the BCG, ICL and rest of stars. The bar shows the median value for the 324 clusters, while the error bars indicate 16th–84th percentiles. The fraction is computed dividing the stellar mass content in each of the components by the total stellar mass inside R500, M*, 500. The straight and dashed bars separate into relaxed and disturbed clusters based on their χDS value at R500. “Relaxed” clusters are the 50 clusters with the highest χDS, while “disturbed” are the 50 clusters with the lowest value of χDS. As before, green bars are for GADGET-X, while magenta ones are for GIZMO-SIMBA results.

In conclusion, and focusing on the ICL, we find that for our whole cluster sample the ICL mass fraction, considering 16th–84th percentiles as the errors, is 0.38 ± 0.10 (0.34 ± 0.09) for GADGET-X (GIZMO-SIMBA). When selecting a subsample of the 50 most relaxed clusters, this fraction increases to $0.{49}_{-0.06}^{+0.05}$ ($0.{45}_{-0.07}^{+0.05}$); while it decreases to $0.{27}_{-0.06}^{+0.08}$ ($0.{24}_{-0.04}^{+0.08}$) for the 50 most disturbed clusters, according to their χ_DS value at z = 0. We note that these fractions are in perfect agreement with those found by Contini et al. (2023) for the most and least concentrated halos in their sample. We explain these differences between the subsamples by their different formation scenarios, with more disturbed clusters being the result of recent merger events, such that galaxies as a whole enter the cluster and remain as bound objects; and relaxed clusters being the result of the disruption of these galaxies and hence slow accretion of their stars into the ICL. In agreement with our results, but with different simulations, Cui et al. (2014) and Rudick et al. (2011) found the ICL fraction to be higher for relaxed clusters, and stated that this can be understood if relaxed clusters are more dynamically evolved than disturbed ones, with the diffuse star particles mainly coming from satellite galaxies undergoing mergers and stripping. This can also be understood within the context of the two-phase formation scenario of galaxy clusters (e.g. Gunn & Gott 1972; Gunn 1977; Ascasibar et al. 2004), that describes an early fast accretion phase building up the central region; followed by a slow accretion phase, where the inner part remains more constant and the mass builds up in the outer regions. This was analysed by Mostoghiu et al. (2019) in THE THREE HUNDRED clusters, concluding that unrelaxed clusters are still in their fast accretion mode, while relaxed clusters have already reached the slower phase. We plan to investigate the origin of the ICL in more detail in a future work, where we will trace the particles back in time. Nevertheless, for the present work we prefer to focus on characterising the ICL at present day, and how it is related with other cluster properties.

6. Relation of ICL to dark matter

In this section we study the relationship between the ICL and the dark matter (DM) component of galaxy clusters. This relationship is particularly interesting since it could provide a way of exploring the DM component by using the ICL, which means probing the DM in clusters using only deep imaging observations. This has been previously investigated using Hubble Space Telescope (HST) data by Montes & Trujillo (2019), finding very promising results pointing to a very similar matter distribution of both components. From the simulations side there has also been a great effort to separate the stellar and DM component of halos, in order to understand their origins and joint evolution. Libeskind et al. (2011) used constrained simulations of the local Universe, while Pillepich et al. (2018b) worked with the cosmological simulation IllustrisTNG. Working with the Cluster-EAGLE simulations, Alonso Asensio et al. (2020) and Deason et al. (2021) also found that the ICL can be used to infer properties of the underlying DM halo, like its shape and boundaries.

6.1. Full cluster sample

In order to explore this relationship, we have computed the volumetric density profile and the velocity dispersion profiles for both the DM and ICL components. The density profile for the ICL is very similar to that in Fig. 4, with the difference that now we are excluding the BCG. For the DM we compute this profile by considering all the DM particles that belong to the central halo. In Fig. 10, the top left panel shows the ratio between these two profiles, ρ_ICL/ρ_DM, stacked for all the 324 clusters. The solid lines are obtained by computing the ratio ρ_ICL/ρ_DM at each R for all the 324 clusters and then computing the median value. The shaded regions show the 16th–84th percentiles, green is for GADGET-X and magenta for GIZMO-SIMBA. We can see that the results are very similar for both simulations, with the ratio showing a reasonably constant slope across different radii, when using logarithmic scale in both axes. The plot in the middle of the upper row shows how this ratio ρ_ICL/ρ_DM can be fitted to a power law of the form y = a ⋅ x^b. The dashed line shows the best fit, while the dotted lines show the 1σ errors. The green lines correspond to the results for GADGET-X, while the ratio and best fit curve for GIZMO-SIMBA are shown in grey, to facilitate the comparison between the two simulations. The plot in the upper right corner is the same as the central one but for GIZMO-SIMBA, with the results for GADGET-X in grey. The legends in the lower part of the plots indicate the parameters of the fit, which are also summarised in Table 1. We want to highlight here the value of the exponent b, which determines the slope of the curves in the upper row of Fig. 10. This parameter takes a value of −1.23  ±  0.04 (−1.13  ±  0.05) for GADGET-X (GIZMO-SIMBA). We note that these results have been obtained considering all the DM particles, that is including also the substructures present in the halo. We have compared the results for both possibilities (excluding or including the substructures) and the conclusions are very similar, with the parameters being compatible within 1σ with each other. For the remaining of this section, the DM mass will refer to the inclusive one, considering all the mass in subhalos as well.

Fig. 10. Relation between ICL and DM radial profiles. Top: ratio between the ICL and DM density profiles. In the left panel, solid lines are the median values, green for GADGET-X and magenta for GIZMO-SIMBA, while the shaded regions are the 16th–84th percentiles. The middle and right panels show the best fit to an exponential curve of the GADGET-X and GIZMO-SIMBA ratios, respectively (the fainter grey lines are for the opposite simulation, for easier comparison between them). The dashed lines depict the best fit (the analytical formula can be seen in the legend), while the dotted lines indicate the 1σ errors. Bottom: same as the top row but for the ratio between ICL and DM velocity dispersion profiles.

Also regarding the values in Table 1, we want to clarify that the parameters – as well as the errors – have been obtained by stacking all the ρ_ICL/ρ_DM profiles for the 324 clusters and then fitting to this one profile. An alternative way of obtaining these values is to fit the profile of each individual cluster and obtain the median best-fit parameters with their percentiles. We prefer the former option since it provides a more robust fit that focuses on the overall trend for our sample, and hence reveals general properties of the relation between ICL and DM. The latter option, that is, fitting for each cluster, considers more the individual variations among clusters and provides cluster-specific values. Although not shown here, we have applied this methodology, and find that the median values of the parameters remain within the 1σ regions, but the errors are considerably larger (3–4 times larger than for the stacked profile for parameter b and 2–3 times for parameter d).

A study by Deason et al. (2021) also explores the connection between ICL and DM in the Cluster-EAGLE simulations but focusing on the “edge” of their distributions, finding that they are both closely related. In that work (cf. their Fig. 4) it can be seen that the logarithmic slope profiles for DM and stars, when stacking all the clusters, show an approximately constant difference between them, which is the same we see in our Fig. 10. Computing the difference between the two logarithmic profiles, the value in Deason et al. (2021) is ∼ − 1, in agreement with the value for the slope b we find (see Table 1).

The bottom row of Fig. 10 focuses on the velocity dispersion profiles of both the DM and ICL components. We compute this profile as the velocity dispersion for all the particles inside each indicated radius (instead of the particles within a thin spherical shell). As for the density profile, in the left panel we show the ratio between the profiles for ICL and DM, σ_v, ICL/σ_v, DM. We see again that the two simulations agree across the majority of radii, and show only slight differences in the innermost part of the cluster, a feature we have already seen and explained in this paper (e.g. Fig. 4). In the two other panels we can see how this ratio can be fitted to a straight line with the form y = c + d ⋅ x (note that the scales are not logarithmic in this case). Again, the values of these free parameters are summarised in Table 1. The value obtained for the slope is 0.02 ± 0.02 (0.09 ± 0.02) for GADGET-X (GIZMO-SIMBA). Both values are significantly low, and even compatible with 0 for GADGET-X, meaning that the ICL and the DM components have velocity dispersion profiles that evolve very similarly.

In a previous study, Montes & Trujillo (2019) used deep observations of galaxy clusters to show that the projected ICL distribution closely follows that of the DM distribution up to ∼140 kpc, thus providing an accurate luminous tracer of DM. This was confirmed by Alonso Asensio et al. (2020) in the Cluster-EAGLE simulations, where they also found that the ratio between the ICL and the total mass density profile has a slope of ∼ − 1, in agreement with what we find for THE THREE HUNDRED simulations (Fig. 10). We still want to note that the two studies mentioned work with projected rather than volumetric quantities, and so the comparison with them should be interpreted with caution. With a different set of simulations, Yoo et al. (2022) developed a new method to measure the similarity between two distributions, further supporting these findings. This result regarding the relation between ICL and DM means that the ICL density profile can be used to infer the DM density profile, which can in turn be integrated to obtain an estimation of the cluster mass. Pillepich et al. (2018b) propose a similar method of inferring halo mass from stellar mass based on the IllustrisTNG simulations.

We further explore the relationship between these two components by introducing the ratio between the ICL and DM velocity dispersion profiles. The velocity dispersion is a very interesting property of clusters, since it can be used as a proxy of cluster mass (see e.g. Allen et al. 2011), and thus for inferring cosmological parameters. Several studies have been devoted to studying this relation between velocity dispersion and cluster mass, and the differences between using the member galaxies or the DM component (see Munari et al. 2013 or Ferragamo et al. 2022 for a study in THE THREE HUNDRED clusters). In turn, the velocity dispersion of the ICL stars can be a useful tool to understand in more detail the kinematics of the cluster, and thus its structure and evolution. Some recent works already present measurements for nearby clusters in the Local Universe (see e.g. Spiniello et al. 2018 for the Fornax cluster or Arnaboldi & Gerhard 2022 for a review). In this work we present the relation between the velocity dispersion of the DM and that of the diffuse stars that follow the gravitational potential of the cluster. We find that the ratio ICL/DM can be described by a straight line with slope close to 0, meaning that they evolve very similarly up to R₅₀₀, so that also one of them can be used to infer the other one (see equation and parameters in Table 1).

6.2. Subsamples of relaxed and disturbed clusters

In the previous Sect. 5.3 we created two subsamples of the 50 most relaxed and disturbed clusters based on their value of the relaxation parameter χ_DS, and we studied how their compositions differed from each other. We found that the ICL fraction is significantly higher for the relaxed clusters, for which this component can even dominate the stellar component of clusters. Now that we have studied the relation of the ICL with the DM component, we also ask ourselves the question of whether this relation depends on the dynamical state of the cluster. In order to check this, we have repeated the previous analysis, this is, the comparison between density and velocity dispersion profiles for ICL and DM components, for these two subsamples of relaxed and disturbed clusters.

We show the results in Fig. 11. The left 2 × 2 panel in this figure is for the density, while the right one is for the velocity dispersion profiles. In each of these 2 × 2 panels, the first column is for the relaxed clusters, GADGET-X on top (in green) and GIZMO-SIMBA on the bottom (in magenta); while the second column is for the disturbed clusters. Similarly to the second and third columns of Fig. 10, in Fig. 11 we show the median value (solid line) of the ratio of the two profiles, ρ_ICL/ρ_DM or σ_v, ICL/σ_v, DM, together with the best fit curve (dashed line) and the 1σ regions for the fits (delimited by dotted lines). In each panel, the grey line corresponds to the opposite subsample for the same simulation, for example in the upper left corner the green lines correspond to “relaxed GADGET-X” clusters, and so the grey lines are the median and best fit for the “disturbed GADGET-X” clusters. This allows for an easier comparison between the relaxed and disturbed subsamples, which is the main goal of this subsection. The equations that describe the fits are indicated in the plot, while the values of the parameters with their 1σ errors are summarised in Table 2.

Fig. 11. For the density (left) and velocity dispersion (right) profiles, ratio between ICL and DM components. Same as Fig. 10 but for the subsamples of the 50 most relaxed and disturbed clusters (see Sect. 5.3). Green is for GADGET-X and magenta for GIZMO-SIMBA simulations. In each panel, the grey line corresponds to the opposite subsample for the same simulation, namely relaxed or disturbed, to allow for an easier comparison between relaxed and disturbed clusters.

Focusing first on the density, we see in Fig. 11 that for relaxed clusters the line is moved up with respect to the disturbed ones, which is a consequence of the ICL fraction being higher for relaxed clusters (see Sect. 5). Looking at the values of the exponents of the fits, parameter b in Table 2, we see that the ratio becomes a bit steeper for the disturbed clusters, although the error is also higher for these clusters and hence the difference is not very significant. Regarding the velocity dispersion, Fig. 11 shows a more relevant difference between relaxed and disturbed clusters. For relaxed clusters the ratio remains almost constant, while the slope is significantly increased for disturbed ones. This difference is more clear for GADGET-X, but is present in both simulations. Parameter d in Table 2 quantifies the slope of this ratio, that for relaxed clusters is compatible with 0 for GADGET-X and a bit higher for GIZMO-SIMBA, but still below the value for the whole cluster sample. This can be explained by the fact that, for relaxed clusters, there are less substructures, less mergers happening and the mass is more concentrated towards the centre. These factors make the velocity dispersion to remain more or less constant beyond a certain point in the profile. For disturbed clusters, substructure and ongoing mergers can cause the velocity dispersion to continue growing with increasing distance from the cluster centre. This effect is more pronounced for the ICL than for the DM, which is less affected by interactions, and hence we see an effect in the ratio σ_v, ICL/σ_v, DM. In the context of our previous conclusions, that relaxed clusters are more dynamically evolved than disturbed ones, this σ_v ratio is constant for relaxed clusters because accreted stars have already had time to relax and trace the underlying gravitational potential. However, for unrelaxed clusters, if the ICL is formed of ejected stars, these have just been stripped from satellites and are not yet fully coupled and in equilibrium with the halo potential.

We have seen in this subsection how the relation between ICL and DM is also affected by the dynamical state of the cluster, affecting both the density and the velocity dispersion profiles, but more importantly the latter. Since the velocity dispersion profile of the ICL can be influenced by many factors, the exact relationship with the dynamical state can be complex and different from cluster to cluster. We have presented here our median results for subsamples of 50 of the most relaxed and disturbed clusters.

7. Conclusions

In this work, we present the first detailed analysis of the intra-cluster light (ICL) within THE THREE HUNDRED simulations project. This consists of a suite of 324 numerically modelled spherical regions centred on the most massive clusters found in a prior DM-only cosmological simulation. As such, it makes a relatively large set of galaxy clusters compared to other state-of-the-art cluster simulations, which typically focus on a handful of objects (for instance, the Cluster-EAGLE simulations, Barnes et al. 2017, include only 30 galaxy clusters). THE THREE HUNDRED clusters constitute a mass-complete sample within the range 14.4 < log(M₅₀₀/M_⊙) < 15.4 at z = 0. The size of the sample allows us to perform a statistical study of the ICL in clusters, as well as to make statistically significant subsamples focusing on different cluster properties. These 324 regions of 15 h⁻¹ Mpc in radius around the cluster centre have been resimulated including full hydrodynamics with two different subgrid physics implementations: GADGET-X, which uses a modified version of the non-public GADGET3 code (Murante et al. 2007; Rasia et al. 2015); and GIZMO-SIMBA, performed with the GIZMO code (Hopkins 2015) and the galaxy formation subgrid models from the SIMBA simulation (Davé et al. 2019). Comparing the results for the two simulations can help to understand the effects of different physical processes and test the robustness of the resulting predictions.

We first define the BCG as the total stellar mass contained within a fixed spherical aperture of 50 kpc (in Sect. 3.3 and Appendix A we study the dependence on the size of this aperture) from the halo centre (located by the AHF halo finder as a peak in the density field). This is a common approach to find the BCG in numerical simulations, with both theoretical and observational arguments supporting it (Pillepich et al. 2018b; Kravtsov et al. 2018). Then, from the rest of the stellar particles, we find the ICL by selecting those that do not belong to any subhalo. In this way we remove the central and the satellite galaxies, leaving us with the particles that are bound to the potential of the cluster itself, but neither the BCG nor any of the satellite galaxies. This diffuse stellar component is generally referred to as ICL (see Fig. 1 for an example representation for one cluster).

As a first check, in Sect. 4 we compare our resulting BCG+ICL mass to observational results by Gonzalez et al. (2013), Kravtsov et al. (2018) and DeMaio et al. (2020) and find our predictions to be in overall agreement with them. The two simulations GADGET-X and GIZMO-SIMBA show highly consistent results down to 0.01 R₅₀₀. Inside this region, the density profile from GADGET-X shows a steeper slope, while GIZMO-SIMBA starts to flatten. Although these scales are very close to the simulations resolution limit, the differences between the two simulations can be attributed to the different stellar formation and AGN implementation (cf. Sect. 2), which eventually leads to higher stellar masses in the very central regions found in the GADGET-X runs.

In Sects. 5 and 6 we solely focus on the ICL of THE THREE HUNDRED clusters. We study its mass fraction and how it depends on halo mass and dynamical state. We also study the relationship between the ICL and the DM components of clusters, verifying whether or not the former could be used to trace the latter. The main results of these sections, which are the main findings of this work, are summarised as follows.

• The ICL mass fraction, computed as M_ICL, Δ/M_*, Δ (with Δ = 200, 500), has a median value of f_ICL, 500 = 0.36 and f_ICL, 200 = 0.29, with 1σ error of ±0.10. This value is independent of halo mass, at least in our considered mass range (Fig. 5).

• The ICL mass fraction shows a clear dependence on the dynamical state of the cluster, which is quantified by the theoretical indicators introduced by Neto et al. (2007): subhalo mass fraction f_s, centre of mass offset Δ_r and – to a lesser extent – with the virial ratio η. With a different approach and based on our previous work (Contreras-Santos et al. 2022), we find a correlation between the ICL fraction and the merging history of the cluster. Clusters that have undergone their last merger more recently present a smaller fraction. There is also a weaker negative correlation between the number of mergers undergone since z ∼ 1.5 and the ICL fraction (see Fig. 8). These results are in agreement with the scenario where relaxed clusters are dynamically more evolved than disturbed ones, and hence star particles in the ICL component mainly come from the satellite galaxies undergoing merging and tidal stripping (Rudick et al. 2011; Cui et al. 2014; Iodice et al. 2020; Ragusa et al. 2021, 2023; Contini et al. 2023).

• The ratio between the volumetric density profiles of the ICL and the DM components follows a power law up to R₅₀₀, with exponent −1.23 ± 0.04 (−1.13 ± 0.05) for GADGET-X (GIZMO-SIMBA). The full equations, together with the 1σ uncertainties of each parameter, can be seen in Table 1. As suggested in previous works (Pillepich et al. 2018b; Alonso Asensio et al. 2020), this relation between the density profiles can be used as a method to infer the DM mass of the halo. Our study, with more than three hundred clusters, provides a robust relation and confirmation of this method.

• The ratio between the velocity dispersion profiles of the ICL and the DM component follows a straight line with slope close to 0, 0.02 ± 0.02 for GADGET-X and 0.09 ± 0.02 for GIZMO-SIMBA (see Table 1 for the full equations). This relation can also be used to infer the DM halo velocity dispersion up to R₅₀₀, which evolves very similarly to that of the ICL.

One limitation of this study that needs to be mentioned is in regard to the selection of the cluster centre, which is in turn used to find the BCG. The halo finder applied by us, namely AHF, locates the halo centre as a peak in the density field, so that we always centre our cluster in the highest density peak, which is in general the BCG. However, some studies have shown that BCGs are not always located at the centre of their clusters, especially when the cluster is unrelaxed (Martel et al. 2014; De Propris et al. 2021). If the BCG is not in the centre, some of its particles might be mislabelled as ICL, and vice versa. Nevertheless, we do not expect these particles to be a majority, especially for the ICL, which is considerably more massive than the BCG, and so we believe our results for the ICL to still be consistent in spite of this.

Finally, we are already planning a future work to study the origin of the ICL in more detail by tracing the particles back in time and investigating where they were formed. Similarly to the work with the IllustrisTNG simulation by Montenegro-Taborda et al. (2023), we will separate the stars into in situ and ex situ, depending on whether they were formed in the ICL itself or somewhere else and then accreted by this diffuse component. This will help us to better understand the relation between the DM component and the ICL in clusters, and in general to gain a better understanding of the process of galaxy cluster formation and evolution.

Acknowledgments

The authors would like to thank Emanuele Contini for carefully reading the manuscript and providing insightful suggestions. This work has been made possible by THE THREE HUNDRED (https://the300-project.org) collaboration. The simulations used in this paper have been performed in the MareNostrum Supercomputer at the Barcelona Supercomputing Center, thanks to CPU time granted by the Red Española de Supercomputación. As part of THE THREE HUNDRED project, this work has received financial support from the European Union’s Horizon 2020 Research and Innovation programme under the Marie Skłodowska-Curie grant agreement number 734374, the LACEGAL project. A.C.S., A.K., I.A.A. and C.D.V. thank the Ministerio de Ciencia e Innovación (MICINN) for financial support under research grants PID2021-122603NB-C21 and PID2021-122603NB-C22. A.K. further thanks Drop Nineteens for “kick the tragedy”. W.C. is supported by the STFC AGP Grant ST/V000594/1 and the Atracción de Talento Contract no. 2020-T1/TIC-19882 granted by the Comunidad de Madrid in Spain. He also thanks the Ministerio de Ciencia e Innovación (Spain) for financial support under Project grant PID2021-122603NB-C21 and ERC: HORIZON-TMA-MSCA-SE for supporting the LACEGAL-III project with grant number 101086388.

References

Appendix A: Dependence on aperture size

In this section we show how the main results shown throughout the paper are affected by the choice of aperture size to separate BCG and ICL. In the previous sections we have used an aperture of 50 kpc from the halo centre to select all the stars that belong to the BCG. Now, we will also use values of 30 and 70 kpc, to see how much the results are influenced by a smaller or larger aperture. Here we like to remind the reader that we apply a 3D aperture, and hence this method is not fully reproducible in observations (see discussion in Sec. 3.1).

A.1. ICL mass fraction

Here we repeat the calculations presented in Sec. 5 for the different fixed apertures considered. In Fig. A.1 we show the dependence of the ICL mass fraction on cluster mass. This is the same as Fig. 5 but including also 30 and 70 kpc and, for simplicity, showing only the median values for each mass bin, without the points themselves and the errors. As before, green depicts GADGET-X results, while magenta is for GIZMO-SIMBA, while lighter to darker colours represent smaller to larger apertures, namely 30, 50 and 70 kpc. Although the values change as expected, that is, a smaller BCG makes a larger ICL, we can see that the same trends remain, with no correlation between the ICL mass fraction and the cluster mass for our considered mass range. We also show only the results for overdensity 500, but note that the same trends remain also for R₂₀₀.

Fig. A.1. Same as Fig. 5, ICL mass fraction as a function of cluster mass, but comparing the values for different aperture sizes to define the BCG, in solid green for GADGET-X and in dash-dotted magenta for GIZMO-SIMBA. Lighter to darker colours indicate increasing size of aperture: 30, 50 and 70 kpc. We only show the values considering overdensity 500.

In Fig. A.2, we show the stellar mass fraction in the two separate components BCG and ICL. This is the fraction of the total stellar mass that belongs to these components, computed as f_BCG, 500 = M_BCG/M_*, 500 or the same for f_ICL. Unlike in Fig. 6 (see Section 5), here we do not show the values for the satellite galaxies, that is stars that do not belong to the BCG neither the ICL, because this value is independent of the aperture size used to define the BCG, which does not affect the joint component BCG+ICL. This figure is similar to Fig. A.1, but in Fig. A.2 we focus on the median values rather than the dependece on cluster halo mass, and we also present results for the BCG mass fraction. We can see that the median value of f_BCG, which was around 10 per cent for 50 kpc, decreases to ∼7 per cent for 30 kpc and increases to ∼13 per cent for the largest aperture of 70 kpc. For the ICL, since the absolute values are larger, the relative changes from one aperture to the other are smaller, with the same trends as for the BCG. We do not show here the fractions when selecting only the most relaxed and disturbed clusters (see Fig. 9), but we note that the trends remain the same regardless of the aperture used.

Fig. A.2. Stellar mass fraction for the BCG and ICL comparing the aperture size used to define the BCG. This is computed dividing the stellar mass content in each of the components by the total stellar mass inside R500, M*, 500. The dots show the median value for the 324 clusters, while the error bars indicate 16th-84th percentiles. Green circles are for GADGET-X, magenta diamonds for GIZMO-SIMBA.

A.2. Relation of ICL to dark matter

We also repeat now the results of Section 6, where we showed the relation between ICL and DM in the clusters, by comparing their mass density and velocity dispersion profiles (see Fig. 10). For the density profile, changing the aperture does not affect the profile except by how far in the inner part the profile gets, the shape remains always the same. This is because of the way the density profile is constructed, and thus there is no reason to repeat it for different sizes of the aperture. However, for the velocity dispersion, the aperture size can change the results, and hence we are depicting this change in Fig. A.3. In this plot, the left column shows the median values of the ratio of the two profiles, σ_v, ICL/σ_v, DM, with lighter to darker colours being for 30, 50 and 70 kpc respectively. The shaded regions are the 16th-84th percentiles. The upper row, in green, is for GADGET-X, while the lower one, in magenta, is for GIZMO-SIMBA. We can see that this ratio does not change very significantly from one aperture to another. The main effect is that, for larger apertures, the ratio is increased, which means that the ICL profile is higher (the DM profile is not affected by the aperture size). This can be explained by the fact that particles in the innermost region, selected only for smaller apertures, have velocities that are more similar within each other, such that the velocity dispersion profile decreases towards the cluster centre.

Fig. A.3. Dependence of the relation between ICL and DM profiles on the aperture size used to define the BCG. First column: ratio between the ICL and DM velocity dispersion profiles, comparing the results for the three different apertures considered to define the BCG and thus the ICL: 30, 50 and 70 kpc. In green for GADGET-X and in magenta for GIZMO-SIMBA. Second column: fit of the previous plots to a straight line, with the equations of these lines indicated in the legend.

The right column in Fig. A.3 shows the fit of these ratios to a straight line with the form y = c + d ⋅ x. The values of these parameters can be seen in the legend of the plot, as well as summarised together with their 1 σ errors in Table A.1. As it could be expected from the left column, the lines are very similar between each other, showing only very slight changes. Moreover, in Table A.1 we can see that, for 30, 50 and 70 kpc, the values of the c and d parameters are within the 1 σ intervals of each other, so that they are all compatible. Hence, we can say that choosing a different size of the fixed aperture that encloses the BCG does not make a significant difference regarding the relation between ICL and DM velocity dispersion profiles.

All Tables

All Figures

Fig. 1. For one of the 324 clusters (cluster 6), 2D projection of all the stellar particles along one axis that belong to, from left to right, the total stellar mass within R500, the ICL and the BCG alone. The same colour intensity is used for all the visualizations with darker tones showing higher densities; concentric circles have radii of 100 kpc (inner), 500 kpc (middle), and 1 Mpc (outer). For the rightmost plots only the inner circles (100 kpc radii) are shown.

In the text

Fig. 2. Dependence of BCG and ICL mass on aperture size. Top: distribution of the ratio between the BCG mass of the THE THREE HUNDRED simulations obtained for apertures: 30 and 50 kpc (first row), and 50 and 70 kpc (second row). The left column (in green) is for GADGET-X, while the right one (in magenta) is for GIZMO-SIMBA. Bottom: same for ICL mass. Since ICL mass decreases with increased aperture, the ratio is inverted now to keep the values above 1, but the interpretation remains the same.

In the text

Fig. 3. GADGET-X (left) and GIZMO-SIMBA (right) BCG+ICL stellar mass within different spherical regions as a function of total halo mass. First row shows the total BCG+ICL stellar mass, squares denote individual measuremente of each of the 324 main clusters, grey symbols show data from Gonzalez et al. (2013) and Kravtsov et al. (2018), as labelled; the shaded region shows the best-fit from Kravtsov et al. (2018). Second row shows the total BCG+ICL mass within 100 kpc aperture and the third one shows the BCG+ICL stellar mass contained between 10 and 100 kpc apertures. Grey symbols show observational estimations from DeMaio et al. (2020).

In the text

Fig. 4. Density profile of the BCG+ICL component (solid) and satellite galaxies (dotted) for the 324 main clusters of GADGET-X (green) and GIZMO-SIMBA (magenta) scaled by cluster spherical overdensity radius, R500. Lines shows the median of each bin, and shaded regions delimit the 16th and 84th percentile of each bin. Inset shows only the BCG+ICL component in logarithmic scale in the X-axis, to focus on the inner part of the profile. Here, the vertical shaded regions show for reference the 16th and 84th percentile of the location of 10 and 100 kpc, as labelled.

In the text

	Fig. 5. ICL mass fraction as a function of cluster mass, considering overdensities of 500 (top) and 200 (bottom). The lines show the median values and the shaded regions the 16th–84th percentiles. Results are shown for both GADGET-X (green solid line) and GIZMO-SIMBA (magenta dash-dotted line). Note that the scale is the same for both plots in the Y-axis, but different for the X-axis.
In the text

	Fig. 6. Stellar mass fraction for the BCG, ICL and rest of stars. The bar shows the median value for the 324 clusters, while the error bars indicate 16th–84th percentiles. The fraction is computed dividing the stellar mass content in each of the components by the total stellar mass inside R500, M*, 500.
In the text

Fig. 7. ICL mass fraction within R200 of the cluster as a function of the three different theoretical parameters introduced by Neto et al. (2007) to quantify the dynamical state (from left to right subhalo mass fraction fs, centre of mass offset Δr and virial ratio η), together with the formation time, zform. The first three panels are coloured by zform, while the fourth one is coloured by fs, as indicated in the colourbars. In all panels, the Spearman correlation coefficient, rS, is indicated in a lower corner. The vertical dashed lines indicate the thresholds used in Cui et al. (2018) to distinguish between “relaxed” and “unrelaxed” clusters: 0.1, 0.04 and 0.15, respectively. Clusters to the left of these lines (value lower than the threshold) are relaxed, and those to the right are unrelaxed.

In the text

Fig. 8. Correlation between ICL mass fraction and merging history of clusters. Top: ICL mass fraction within R200 of the cluster as a function of the time since the most recent merger. The mergers are defined following Contreras-Santos et al. (2022), as mass increases of at least 25% of the initial mass in less than half a dynamical time. Bottom: ICL mass fraction as a function of the number of mergers experienced by the cluster since z ∼ 1.5. The clusters are binned by number of mergers experienced, with the bar plot below indicating the number of clusters in each bin. The dots correspond to the median values, while the errors are 16th–84th percentiles.

In the text

Fig. 9. Stellar mass fraction for the BCG, ICL and rest of stars. The bar shows the median value for the 324 clusters, while the error bars indicate 16th–84th percentiles. The fraction is computed dividing the stellar mass content in each of the components by the total stellar mass inside R500, M*, 500. The straight and dashed bars separate into relaxed and disturbed clusters based on their χDS value at R500. “Relaxed” clusters are the 50 clusters with the highest χDS, while “disturbed” are the 50 clusters with the lowest value of χDS. As before, green bars are for GADGET-X, while magenta ones are for GIZMO-SIMBA results.

In the text

Fig. 10. Relation between ICL and DM radial profiles. Top: ratio between the ICL and DM density profiles. In the left panel, solid lines are the median values, green for GADGET-X and magenta for GIZMO-SIMBA, while the shaded regions are the 16th–84th percentiles. The middle and right panels show the best fit to an exponential curve of the GADGET-X and GIZMO-SIMBA ratios, respectively (the fainter grey lines are for the opposite simulation, for easier comparison between them). The dashed lines depict the best fit (the analytical formula can be seen in the legend), while the dotted lines indicate the 1σ errors. Bottom: same as the top row but for the ratio between ICL and DM velocity dispersion profiles.

In the text

Fig. 11. For the density (left) and velocity dispersion (right) profiles, ratio between ICL and DM components. Same as Fig. 10 but for the subsamples of the 50 most relaxed and disturbed clusters (see Sect. 5.3). Green is for GADGET-X and magenta for GIZMO-SIMBA simulations. In each panel, the grey line corresponds to the opposite subsample for the same simulation, namely relaxed or disturbed, to allow for an easier comparison between relaxed and disturbed clusters.

In the text

	Fig. A.1. Same as Fig. 5, ICL mass fraction as a function of cluster mass, but comparing the values for different aperture sizes to define the BCG, in solid green for GADGET-X and in dash-dotted magenta for GIZMO-SIMBA. Lighter to darker colours indicate increasing size of aperture: 30, 50 and 70 kpc. We only show the values considering overdensity 500.
In the text

Fig. A.2. Stellar mass fraction for the BCG and ICL comparing the aperture size used to define the BCG. This is computed dividing the stellar mass content in each of the components by the total stellar mass inside R500, M*, 500. The dots show the median value for the 324 clusters, while the error bars indicate 16th-84th percentiles. Green circles are for GADGET-X, magenta diamonds for GIZMO-SIMBA.

In the text

Fig. A.3. Dependence of the relation between ICL and DM profiles on the aperture size used to define the BCG. First column: ratio between the ICL and DM velocity dispersion profiles, comparing the results for the three different apertures considered to define the BCG and thus the ICL: 30, 50 and 70 kpc. In green for GADGET-X and in magenta for GIZMO-SIMBA. Second column: fit of the previous plots to a straight line, with the equations of these lines indicated in the legend.

In the text