The Transition Region between Brightest Cluster Galaxies and Intra-Cluster Light in Galaxy Groups and Clusters

As mentioned in Section 1, the BCG+ICL light distribution is, observationally, derived by profile fitting with the use of functional forms, usually double/triple Sérsic (Sérsic 1968) or composite profiles. In order to be consistent with the growths of bulge and disk in the semi-analytic model, in Contini & Gu (2021) we modelled them by assuming a Jaffe profile (Jaffe, 1983) for the bulge, and an exponential one for the disk. It is worth stressing that BCGs are normally ellipticals, but in some cases they have time to grow a new disk with the absence of major mergers and a strong AGN feedback. In our sample, most of them can be classified as early-type galaxies with negligible disks, but as stated above, to be consistent with the growths of bulge and disk, we consider both components in the analysis.

The Jaffe profile reads as follows:

where $M_{J}$ is the total mass of the bulge and $r_{J}$ the scale length. By integrating within a sphere of radius $r$, we obtain the mass

The exponential profile for the disk, instead, has the following form:

where $M_{D,tot}$ is the total stellar mass in the disk, and $R_{D}$ the scale radius of the disk.

The distribution of the ICL is modelled by assuming a modified version of the NFW profile (Navarro et al., 1997):

where $\rho_{0}$ is the characteristic density of the halo at the time of its collapse and $R_{s}$ is its scale radius. The concentration of a DM halo is defined as

where $R_{200}$ is the virial radius of the halo. The modification with respect to an NFW profile described in Equation 4 is the introduction of a free parameter in Equation 5. We define the ICL concentration as

For further details about the profiles, and in particular about the free parameter $\gamma$, we refer the reader to Contini & Gu (2020, 2021). We just stress here that, being the following analysis focused at the present time, the value of $\gamma$ used is 3, which means that the ICL is more concentrated than the DM halo in which it resides. As discussed in Contini & Gu (2020), this is in good agreement with recent theoretical works (e.g. Harris et al. 2017; Contini et al. 2018; Pillepich et al. 2018) and observational studies (e.g. Montes & Trujillo 2019). Moreover, we consider stellar haloes surrounding the BCG as being part of the ICL. This caveat is important in the context of comparisons with other studies, both observational and theoretical, which might consider stellar haloes to be either part of the BCG or even a third component.

In the following analysis we use two radii: (a) a scale radius of BCG+ICL system, $R_{BCG+ICL,X\%}$, which is defined as the radius that contains $X\%$ of the total BCG+ICL stellar mass, and the transition radius, $R_{transition}$, defined as the radius where the differential mass in ICL is 90% of the total (BCG+ICL). Both radii are derived by using the bulge, disk and ICL profiles described above.

In this paper we take advantage of two sets of N-body simulations, from which we extracted the mergers trees that are the input of our semi-analytic model. The two sets are different in many ways: one comprises 27 high-resolution zoom-in cluster simulations, and the other one is a cosmological simulation with a box size of around 170 Mpc/h. Below the important details of both simulations.

DIANOGA. This set of simulations has been described in detail in Contini et al. 2012 (but further information can be found in Bassini et al. 2020), while here we just mention the most important features. DIANOGA comprises 27 high-resolution numerical simulations of regions around galaxy clusters, carried out assuming the following cosmological parameters (see Bonafede et al. 2011; Fabjan et al. 2011): $\Omega_{m}=0.24$ for the matter density parameter, $\Omega_{\rm bar}=0.04$ for the contribution of baryons, $H_{0}=72\,{\rm km\,s^{-1}Mpc^{-1}}$ for the present-day Hubble parameter, $n_{s}=0.96$ for the primordial spectral index, and $\sigma_{8}=0.8$ for the normalization of the power spectrum. The latter is expressed as the r.m.s. fluctuation level at $z=0$, within a top-hat sphere of $8\,h^{-1}$Mpc radius. For all simulations the mass of each DM particle in the high resolution region is $10^{8}\,M_{\odot}/h$, and a Plummer-equivalent softening length is fixed to $\epsilon=2.3$ kpc/h in physical units at $z<2$, and in comoving units at higher redshift. The data are stored in 93 output times, spanning a redshift range between 60 and 0. Main haloes have been identified using a standard friends-of-friends (FOF) algorithm, while the algorithm SUBFIND (Springel et al., 2001) has been used to decompose each FOF group into a set of disjoint substructures, identified as locally overdense regions in the density field of the background halo. Similarly to previous studies, only substructures that contain at least 20 bound particles are considered to be genuine. At $z=0$, our sample of halos counts 361 groups and clusters with virial mass larger than $10^{13}\,M_{\odot}/h$ and up to more than $10^{15}\,M_{\odot}/h$.

NJ167. NJ167 is a DM-only cosmological simulation of ($166.67{\rm\,Mpc/h})^{3}$ box, carried out with $1024^{3}$ DM particles from redshift $z=127$ to the present time, by using an improved version of the TreePM code P-GADGET-3 (Springel, 2005). It can be considered a medium/high-resolution simulation since the mass of each particle is $3.7\cdot 10^{8}M_{\odot}/h$ and the softening length $\epsilon=4$ kpc/h. The data have been stored in 128 discrete snapshots and the following Planck 2018 (Planck Collaboration et al., 2020) cosmology has been used: $\Omega_{m}=0.31$ for the total matter density, $\Omega_{\Lambda}=0.69$ for the cosmological constant, $n_{s}=0.97$ for the primordial spectral index, $\sigma_{8}=0.81$ for the power spectrum normalization, and $h=0.68$ for the normalized Hubble parameter. Also for this simulation we used SUBFIND to isolate subhaloes that could retain at least 20 bound particles. Our final sample of haloes at $z=0$ counts 1450 groups and clusters with virial mass larger than $10^{13}\,M_{\odot}/h$ and up to almost $10^{15}\,M_{\odot}/h$.

The goal of this section is to link an observable radius, $R_{BCG+ICL}$ scale, with the transition radius, i.e. the radius at which the ICL dominates the total distribution (BCG+ICL), by studying their relation with other two observables, the total BCG+ICL and halo masses.

In Figure 1 we plot the two radii as a function of the total BCG+ICL mass, for different percentages used in the definition of the BCG+ICL scale radius and separately for the two samples used (different columns). It appears clear, at first glance, that there is not much difference between the two sets of simulations, except for some scatters. What really appears interesting is that, while the transition radius (blue lines) is independent of the BCG+ICL mass, the scale radius (black lines), independently of the percentage chosen to define it, shows a clear trend by increasing from low to high mass BCG+ICL systems. Moreover, both radii show a considerable scatter at all stellar masses. Another important feature is that increasing the percentage when defining $R_{BCG+ICL}$ leads to sample $R_{transition}$ at decreasing BCG+ICL stellar masses. Indeed, the BCG+ICL scale radius (small percentage) in the top panels samples the transition radius of intermediate-high mass systems ($12.0<\log M_{BCG+ICL}<12.6$), while high percentages, such as 65% used in the bottom panels, sample the transition radius of low mass systems ($11.0<\log M_{BCG+ICL}<11.5$).

Similar results and trends are found in Figure 2, which shows $R_{BCG+ICL}$ and $R_{transition}$ as a function of halo mass. Again, $R_{transition}$ is independent of the halo mass, and the trend of the percentage used in the definition of $R_{BCG+ICL,X\%}$ with mass still holds, in the sense that small percentages sample the transition radius of BCG+ICL systems in high mass haloes, while larger percentages sample those in low mass haloes. Given the fact that the transition radius is independent of both BCG+ICL and halo masses, we have calculated the typical ranges of $R_{transition}$ for different percentages used in its definition, and report them in Table 1 separately for DIANOGA and NJ167 simulations. Regardless the percentage adopted, the ranges are all considerably wide, and the mean value clearly increases with increasing percentage. When the mass of the BCG and ICL are the same, i.e. 50% in the definition of $R_{transition}$, the region of transition can be very close to the halo centre. Indeed, the mean value for DIANOGA (NJ167) is $19\pm 8$ ($26\pm 16$), which means that in some cases the transition radius can be as small as around 10 kpc. On the other hand, in the extreme case where we consider as transition the distance at which the ICL accounts for 90% of the total (differential) BCG+ICL mass (the case shown in the plots), the transition radius can be as large as 100 kpc, with an average value of around 60 kpc. These results are consistent with what we found in Contini & Gu (2021) and with several observational measurements (Zibetti et al., 2005; Seigar et al., 2007; Iodice et al., 2016; Montes et al., 2021; Gonzalez et al., 2021). We will discuss this point in more detail in Section 4.

The aim of the analysis in this section is to find an answer (if possible) to the following questions: assuming we measure the total BCG+ICL stellar mass, is there a way to separate the ICL from the BCG in a statistical manner? If yes, what would the typical error be as a function of the total BCG+ICL and/or halo masses? The answer to the first question can be found if and only if the ICL mass correlates well with the total BCG+ICL mass, while the answer to the second one strongly depends on how significant the correlation would be. Given the large samples of haloes we can count on, and their wide ranges in terms of mass, we can address these points in a statistical way.

Figure 3 shows the relation between the stellar mass in ICL versus the total BCG+ICL mass, for our two samples, Dianoga (black stars) and NJ167 (red squares). Overplotted the linear fits (black and red lines), done separately for the two sets of data and the details of which are reported in Table 2. The two quantities appear to be in a significant linear relation, with an important scatter for objects with $\log M_{BCG+ICL}\lesssim 12.0$ ¹¹1This scatter is mainly due to the higher number of groups (which host less massive BCG+ICL systems) with respect to clusters (which host more massive BCG+ICL systems) and to the dispersion in the relation between the halo mass and its concentration, with the latter being an important parameter of the ICL profile., but with small scatters in the fit (see the values in Table 2). The two fits are somewhat different, probably due to the different number of haloes and the intrinsic different cosmology of the two simulations, but the slopes of the fits are consistent to each other within $3\sigma$ (and so are the intercepts).

In order to answer to the questions above, we now want to compare the true value of the ICL mass with that extrapolated from the fit, by knowing just the total BCG+ICL mass. This is done in Figure 4, where we use the relations reported in Table 2 to infer the values of ICL_fit to plot against the true values, ICL_true. Clearly, the data have to stay as close as possible to the blue line in the plot, which represents a one-to-one relation. Indeed, systems with a true value of ICL mass larger than $\log M_{ICL}\sim 11.6-11.7$ tend to populate very close to the blue line, while for smaller ICL masses the scatter is quite large and with a trend of an overestimation of the true values (the value from the fit larger than the true one).

To quantify how much the inferred mass of the ICL from the fit can be reliable (or not), and the trend of the scatter with both BCG+ICL and halo masses, in Figure 5 we plot the increase/decrease, in percentage, between the value of the ICL mass inferred from the fit and the true ICL mass. The left panel shows this information for both samples (Dianoga with black stars and NJ167 with red squares) as a function of BCG+ICL mass, while in the right panel the information is as a function of halo mass. In both panels, the blue line represents a zero increase along BCG+ICL or halo masses. As expected from the previous figure, at low BCG+ICL or halo masses, the scatter is wide (see the comment of Figure 3) and the value of the ICL inferred from the fit can be three (and even more) times higher than the true value. However, the important features in this figure is shown by the right sides of both panels. Indeed, for $\log M_{BCG+ICL}>12$ and for halo masses on cluster scales ($\log M_{Halo}>14$), the scatter is limited between $\pm$30%. We will come back on this point in the next section with a full discussion of the implications of these results.

The definition of a transition region from the BCG to the ICL is not as trivial as it could seem. Indeed, the two components overlap in a region that can be quite extended depending on the profiles they follow. The observational methods used to isolate the ICL have all pros and cons, as argued in Section 1, which might bring to underestimate/overestimate the amount of stellar mass in ICL. A simple cut in surface brigthness, for instance, could severely underestimate the ICL in the innermost regions that are BCG dominated and similarly, a profile fitting might give the same problem simply because it cannot be precise in the region where the two components overlap.

Our model provides the mass of the BCG and that of the ICL separately, as well as their distributions, which means that we have all the information required to quantify the transition radius of a given BCG+ICL system. We have studied the relations between the transition radius, $R_{transition}$, with BCG+ICL and halo masses, assuming that $R_{transition}$ is the distance from the centre where the ICL starts to account for 90% of the total (locally) BCG+ICL mass. We also tried several other percentages and in all cases (see Table 1), the transition radius does not depend on BCG+ICL and halo masses. In the case shown in figures 1 and 2, where the percentage is 90%, the transition radius has an average value of around 60 kpc and a scatter of about $\pm$40 kpc (the precise values depend on the set of simulations), which means that, on average, it can be as small as 20 kpc and as large as 100 kpc. These values are in good agreement with previous results (discussed below in this section), but given the consistency and the large scatter around the mean value, they also point out the difficulty of finding a proxy that can constrain the region of transition for a given BCG+ICL system.

The transition radius, as mentioned above, can be rougly estimated also from observational studies via profile fitting, similarly to what our model does. A way that can help in measuring $R_{transition}$ is the use of a proxy that is directly observable, such as a radius that contains a given percentage of the total BCG+ICL mass. The advantage of this approach consists on the fact that there would not be any need to separate the two components. If this radius, which we called BCG+ICL scale radius, $R_{BCG+ICL}$, depends on the BCG+ICL or halo masses, once overplotted with $R_{transition}$, the region of the overlap is an indication of the transition radius for a given BCG+ICL/halo mass range obtained directly without separating the two components. Moreover, given the fact that the scale radius can be defined by using different percentages, if changing them provides more overlap with the transition radius, we can probe different BCG+ICL/halo mass ranges.

This was exactly our phisolophy in the analysis done in Section 3.1, where we plot $R_{BCG+ICL}$ by using different percentages together with the transition radius. In our analysis we used 20%, 50% and 65% to define the scale radius $R_{BCG+ICL}$, and the general trend we find is that an increasing percentage traces the transition radius of decreasing BCG+ICL/halo masses. The scale radius is found to increase with BCG+ICL/halo mass, independently of the percentage chosen. This trend has a quite simple explanation: first of all, the mass in a small halo is distributed in a smaller region with respect to larger haloes and, secondly, less massive BCG+ICL systems are likely to be hosted by group size haloes, which are also more concentrated than cluster size haloes. A higher concentration of the halo, according to our ICL profile (see Equation 4), reflects in a more concentrated ICL, which in turns means a smaller scale radius.

Quantitatively speaking, the scale radius $R_{BCG+ICL}$ defined by using different percentages cannot be used as a tracer of the transition radius $R_{transition}$ for different BCG+ICL/halo masses, given the large scatters in the relations. Nevertheless, our findings suggest that the region of transition between the BCG and the ICL does not depend on the total mass or halo mass, being universal independently of the system considered. This is an important result, considering the fact that an average of 60 kpc can be considered a very small shell within which the BCG dominates the distribution. That region is even smaller if we define the ICL to dominate the distribution by considering smaller percentages. For instance, if we define $R_{transition}$ as the radius where the local contribution of the ICL is just higher than 50% of the total BCG+ICL mass, $R_{transition}$ becomes, on average, around 20-25 kpc with a scatter of about 10-15 kpc (see Table 2).

Unfortunately, for the reasons discussed aboved, $R_{transition}$ cannot be taken blindly as the distance after which most of the ICL is distributed for all the range of BCG+ICL/halo masses considered, because the region from the centre to the transition radius contains, proportionally, different percentages of ICL from group-size to cluster-size haloes. Indeed, as we can infer from figures 1 and 2, while in less massive systems (either in terms of BCG+ICL or halo masses) $R_{transition}$ is much larger than, for example $R_{BCG+ICL,20-50\%}$, in more massive systems it is much smaller. This means that the formers have much more ICL concentrated in the innermost regions (as we already discussed) than the latters. Put in other words, while on cluster scales the transition radius can be considered as the radius from which most of the total ICL is distributed, on group scales it cannot be, because on these scales a larger fraction (with respect to the total) of ICL distributes in the same region where the mass of the BCG does.

Our average values of the transition radius are in good agreement with several observational results. Among the first measurements, it is surely worth mentioning $R_{transition}\sim 80$ kpc by Zibetti et al. (2005). The authors, by stacking the profiles of around 700 galaxy clusters taken from the SDSS, found a flattening of the profile at that radius, and interpreted it as the transition to an ICL dominated region. Such a value was later confirmed by other works (Gonzalez et al. 2005; Seigar et al. 2007; Iodice et al. 2016 and others) which showed a similar break of the BCG+ICL profile in the range 60-80 kpc. More recently, Zhang et al. (2019) found similar numbers with a sample of 300 galaxy clusters in the redshift range $0.2<z<0.3$. These authors analyzed the BCG+ICL radial profile imposing a triple Sérsic model and concluded that the transition from the BCG to the ICL is just outside 100 kpc, arguing that this can be interpreted as a potential cut in order to separate the two components. Their conclusion is in good agreement with the discussion above, where we suggest the transition radius as a cut to separate BCG and ICL on cluster scales. Close numbers have been found also by Montes et al. (2021) and Gonzalez et al. (2021), $\sim 75$ kpc and $\sim 70$ kpc, respectively.

Even more recently, Chen et al. (2021) analyzed a sample of around 3000 clusters at $0.2<z<0.3$ from the SDSS and decomposed the stellar surface mass profile of BCG+ICL in three components, assuming a de Vaucouleurs inner profile (BCG dominated region), a outer profile that follows the DM distribution (ICL dominated), and a transitional component that distributes mostly between 70 kpc and 200 kpc. This transitional component, argued by the authors to be necessary in order to fully explain the distribution of the ICL, accounts for around 17% of the total diffuse mass on scales between 50 kpc and 300 kpc, and peaks at 100 kpc with 25%. They interpreted this region as the transition between the BCG dominated part of the distribution, where the ICL forms from processes that are connected to the growth of the BCG itself, such as mergers and stripping in the central regions, and the ICL dominated part of the distribution, where the ICL forms via processes that are not linked to the growth of the BCG, such as stripping of infalling satellites and pre-processing within infalling groups, that are rather connected to the richness of the cluster. The extent of what the authors called sphere of influence (see their paper for more details) can shed some light on the processes responsible for the ICL formation.

In Section 4.1 we have seen that the region of transition between a BCG and an ICL dominated parts of the stellar mass distribution of the two components together can be identified by defining a transition radius, which is found to be independent of both BCG+ICL and halo masses. However, despite there is no relation between such a radius and either of the two masses, the amount of ICL (with respect to the total BCG+ICL mass) in the region of overlap of the two components can be very different from halo to halo, due to the mass-concentration relation. This means that, even though the transition radius gives an idea of where the ICL starts to dominate the distribution, it cannot be directly taken as a physical separation between the two components in the whole halo (or BCG+ICL) mass range investigated.

This conclusion leads to the discussion of the analysis done in Section 3.2 where, by means of the information provided by our semi-analytic model and our BCG/ICL profile modelling, we analyzed the relation between the ICL mass alone and the total BCG+ICL mass, with the goal to provide a way to separate the two components. Clearly, given the fact that the information comes from a semi-analytic model, there are caveats that are worth being discussed, and we will do it below in this section.

From Figure 3 and Table 2 we can see that the $M_{ICL}-M_{BCG+ICL}$ relation is very tight, scattered at low $M_{BCG+ICL}$ due to the increased number of objects at those scales compared to more massive systems. Two questions arise naturally: given that the $M_{ICL}-M_{BCG+ICL}$ relation is tight, is it possible to extract the ICL mass from it? And if so, what is the typical error that one would make by doing it? We know the total BCG+ICL mass, which is also an observable, and we extract the ICL mass from the relation found in Figure 3. How does the ICL mass extracted from that relation compare with the true ICL mass? Put in other words, we can test the reliability of the fits by comparing the ICL mass extracted from them and the true one. If the two ICL masses are comparable within a reasonable scatter, the fits are reliable and the procedure can be used to separate observed BCG and ICL masses.

We have addressed this point in Figure 4, aiming to find an answer to the questions above. We made use of the fits reported in Table 2 to derive the ICL mass extracted from the $M_{ICL}-M_{BCG+ICL}$ relation and compared it with the true ICL mass. We found that for systems having a true ICL mass larger than $\log M_{ICL}\sim 11.6-11.7$, the fits provide a reasonable estimate of it, while below that value the scatter is quite large and the fits tend to overestimate the true value of the ICL mass. These low values of ICL masses likely belong to systems with low BCG+ICL mass, which mostly reside in group size haloes.

In order to quantify the level of accuracy of the fits in extracting the ICL mass from the total, as a function of BCG+ICL and halo masses, in Figure 5 we showed the difference in percentage between the ICL mass extracted from the fits and the true one, both as a function of BCG+ICL and halo masses. The result is quite clear: in systems with low BCG+ICL mass (or halo mass) the scatter is very large, and such a method can underpredict the true value of the ICL mass by 50% or, even worse, it can overpredict the true value by even a factor of three in the worst cases. Clearly, for these objects the fits do not provide a reasonable estimate of the ICL mass, but, at higher BCG+ICL or halo masses, the picture changes. For BCG+ICL systems with $\log M_{BCG+ICL}>12$, or for halo masses on cluster scales, $\log M_{Halo}>14$, the difference in percentage between the extracted and true values is confined between $\pm$30%. Assuming that the scatters in these plots are not mainly driven by the number statistics, $\pm$30% can be considered a reasonable error for a tentative way to separate the BCG from the ICL, in terms of mass estimate.

Typical errors on the observed estimates of the ICL mass/luminosity in the literature can be higher than the statistical error quoted above. As mentioned a few times in this paper, the most common observational methods to isolate the ICL from all the other contributions, BCG included, are all affected by uncertainties that are more or less significative depending on the method used and the quality of the data. Our $M_{ICL}-M_{BCG+ICL}$ relation, at least on cluster scale, offers the opportunity to have an indicative idea of the amount of ICL in objects of a given BCG+ICL or halo masses. It is worth noting that the $M_{ICL}-M_{BCG+ICL}$ relation (and so the fits) is independent of the particular profiles used to describe the mass distribution of BCG and ICL, since it accounts for the total ICL mass within the virial radius $R_{200}$. This means that the prediction of our model can be observationally tested, no matter the profile used to describe the mass distribution of the two components. Nevertheless, there are caveats that must be discussed.

In fact, the validity of the results presented in this work depends on the accuracy of our semi-analytic model in predicting the right amounts of the BCG and ICL masses. Starting from the latter, it has been shown in several works, theoretical (e.g., Contini et al. 2014; Tang et al. 2018; Henden et al. 2020; Cañas et al. 2020; Alonso Asensio et al. 2020) but most importantly observational (e.g., Montes & Trujillo 2014; Edwards et al. 2016; Harris et al. 2017; Alamo-Martínez & Blakeslee 2017; Montes & Trujillo 2018; DeMaio et al. 2018; Edwards et al. 2020; DeMaio et al. 2020; Spavone et al. 2020; Raj et al. 2020; Ragusa et al. 2021; Furnell et al. 2021; Chen et al. 2021; Yoo et al. 2021), that our semi-analytic model provides a good description of the formation and evolution of the ICL. Not just its mass at different redshifts, or the fraction as a function of halo mass, but several other properties that have been confirmed by observational studies. Our model is able to predict when the bulk of the ICL starts to form, its colors and metallicity, and the main mechanisms responsible for its formation and evolution. For what concerns the mass of the BCG, this is a long-standing problem. Semi-analytics models, in general, have always had in the past the problem of overpredicting the low and high mass end of the mass/luminosity function, i.e. including the most massive galaxies such as the BCGs. As time passed by, the problem has been largely alleviated by the inclusion of feedback, and in regards of the discussion here, by the inclusion of the ICL formation itself (see, e.g., Contini et al. 2014, 2017 and references therein). In the light of these reasons, we believe that the predictions in this work are reliable and directly comparable with observational measurements.