© ESO, 2014

1. Introduction

Since its first discovery by Zwicky (1951) and to the most recent works (Guennou et al. 2012; Burke et al. 2012; Adami et al. 2012), the intracluster light (ICL) has gained increasing interest because it can help us understand both the assembly history of galaxy clusters and its contribution to the baryonic budget. The ICL consists of stars that are bound to the cluster potential after being stripped off member galaxies as they interacted and merged with either the brightest cluster galaxy (BCG) or the other member galaxies (Murante et al. 2004; Sommer-Larsen et al. 2005; Monaco et al. 2006; Murante et al. 2007; Conroy et al. 2007; Puchwein et al. 2010; Rudick et al. 2011; Cui et al. 2014; Contini et al. 2014). The ICL signature can be seen in the surface brightness (SB) profile of the BCG as an excess of light with respect to the typical r^1/4 law (de Vaucouleurs 1953). Gonzalez et al. (2005) show that a double r^1/4 model provides a better fit to the BCG+ICL SB profile and that the ICL has a more concentrated profile than does the total cluster light (see also Zibetti et al. 2005).

The origin of the ICL strictly connects it to the evolutionary history of the clusters, thus, we can recall the assembly history of the clusters by studying the ICL properties. The ICL colors can provide information on the timescales involved in ICL formation and on its progenitors when compared to BCG colors. Some works found that ICL colors are consistent with those of the BCG (e.g., Zibetti et al. 2005; Krick & Bernstein 2007; Pierini et al. 2008; Rudick et al. 2010), suggesting that the ICL has been triggered by ongoing interactions among cluster members and the BCG. The merging cluster in the sample of Pierini et al. (2008) and some compact groups (Da Rocha & Mendes de Oliveira 2005) represent an exception showing bluer colors for the ICL, hinting at either in-situ star formation or blue dwarf disruption after interaction.

Usually the ICL is found to be strongly aligned with the position angle (PA) of the BCG (Gonzalez et al. 2005; Zibetti et al. 2005), but there are cases of misalignment and/or prominent features or plumes (Mihos et al. 2005; Krick & Bernstein 2007). Studying the connections between the ICL spatial distribution and the presence of cluster substructures can shed light on the origin of the ICL and its connection to the assembly history of the cluster. ICL plume-like structures that bridge the BCG and other galaxies, arcs and tidal streams of ICL have been found by many works (e.g., Gregg & West 1998; Calcáneo-Roldán et al. 2000; Feldmeier et al. 2004; Krick et al. 2006; Da Rocha et al. 2008). According to simulations these features trace recent interactions and merger events between galaxies and clusters and they are supposed to last only ~1.5 times their dynamical timescale because of disruption by cluster tidal field (Rudick et al. 2009). Adami et al. (2005) and Krick & Bernstein (2007) also found an association between ICL sources and infalling groups of galaxies, and they used it to infer the dynamical evolution of the clusters.

Beside characterizing the ICL properties and the specific evolution of a single cluster, the ICL can be put in a much more comprehensive context by determining its contribution to the total stellar cluster mass and, as a consequence, to the baryon fraction. Observational studies show fractions of ICL ranging from a few percent of the total light up to half of it (Feldmeier et al. 2004; Da Rocha & Mendes de Oliveira 2005; Zibetti et al. 2005; Krick & Bernstein 2007; Gonzalez et al. 2007; Da Rocha et al. 2008; Guennou et al. 2012; Burke et al. 2012; Adami et al. 2012), depending on enclosing radius and cluster mass. On top of this, there is no common definition of ICL among both observational works and simulations. Ideally, the ICL consists of the residual light after having subtracted the contribution of all galaxies, including the BCG. However, both choosing the separation between the BCG and the ICL and determining the best fit model of member galaxies is a difficult task. As a consequence some studies prefer to focus on a BCG+ICL map and mask other members (Gonzalez et al. 2005, 2007), while other authors chose to mask all galaxies down to different arbitrary surface brightness levels (Zibetti et al. 2005; Krick & Bernstein 2007; Burke et al. 2012), and finally Da Rocha & Mendes de Oliveira (2005) and Guennou et al. (2012) remove all the galaxy contribution via a wavelet technique. Different ICL detection methods can suffer from different systematics or bias thus providing discordant ICL fractions as shown for simulations (Cui et al. 2014). This variety of ICL definitions can explain part of the lack of a general consensus on the effective role played by the ICL in the cluster baryon budget.

Moreover the fraction of ICL can correlate with cluster properties, such as mass, projected distance, and redshift, depending on the dominant process and epoch at which they occur (see Krick & Bernstein 2007, for a comprehensive description of the origin of these correlations). Guennou et al. (2012) find only a weak correlation between the ICL content and the cluster velocity dispersion/mass, and there is no variation in the amount of ICL between z = 0.4 and z = 0.8. The absence or mildness of these trends is also confirmed at lower redshifts, i.e., z< 0.3, (Zibetti et al. 2005; Krick & Bernstein 2007). These findings are inconsistent with most of the previous results from both cosmological and analytical simulations, which generally agree with an increasing ICL fraction as cluster mass grows (Murante et al. 2004; Lin & Mohr 2004; Purcell et al. 2007; Watson et al. 2012). However, recent simulations suggest a much weaker dependence of the ICL fraction on cluster mass (Murante et al. 2007; Dolag et al. 2010; Puchwein et al. 2010; Martel et al. 2012; Cui et al. 2014).

Apparently ICL is a promising and complementary way to understand the mechanisms occurring in galaxy cluster and their constituents, however there are two main disadvantages. First, the ICL features typically have extremely faint surface brightnesses of ~1% of the brightness of the night sky, making their study extremely difficult. Second, the surface brightness dimming increases with redshift as (1 + z)⁴. As a consequence, detecting the ICL is very difficult and there are only a few detections at z> 0.3 (Jee 2010; Guennou et al. 2012; Burke et al. 2012; Adami et al. 2012; Giallongo et al. 2014).

In this paper we present our ICL detection and measurement method and we obtained the results from optical images of MACS1206.2-0847 (hereafter MACS1206), one cluster in the Cluster Lensing And Supernova survey with Hubble (CLASH) sample (Postman et al. 2012). Overall, this cluster is one of the most massive (M₂₀₀ = 1.41 × 10¹⁵M_⊙) among the CLASH sample and it is located at a medium redshift, z ~ 0.44, with plenty of ancillary information, so it is a suitable case for testing the performances of our ICL detection method. The CLASH survey comprises 25 massive clusters of galaxies in the redshift range 0.2 ≲ z ≲ 0.9. Among these, 14 have been selected for spectroscopic follow-up at the VLT. At completion, both photometric and dynamical properties of each cluster will be available, allowing the study of ICL and its connection to cluster properties over a wide redshift range. Using deep multiband images from Subaru, we studied the colors and the morphology of the ICL in MACS1206, as well as its connection to cluster substructures and its contribution to the total baryon budget. We then compare these results with those we obtain applying different ICL detection methods, in order to explore the advantages and disadvantages of each method and to reveal possible systematics in each method.

In Sect. 2 we show the dataset we used and the details of the reduction, and in Sect. 3 we explain our ICL detection and measurement method. Section 4 describes our results in terms of both ICL properties and its contribution to the total cluster light and mass. We discuss our results in Sect. 5, and in Sect. 6 we draw our conclusions and future prospects.

Throughout this paper we use H₀ = 70 km s^-1 Mpc^-1, Ω_M = 0.3, and Ω_Λ = 0.7, which gives 5.685 /′′at z = 0.44, the distance of MACS1206.

2. Data

CLASH is one of the three multicycle treasury programs of HST targeting 25 relaxed galaxy clusters with mass range 5−30 × 10¹⁴M_⊙  and redshift range 0.2 ≲ z ≲ 0.9 and providing images for each cluster in 16 passbands using WFC3/UVIS, WFC3/IR, and ACS/WFC (see Postman et al. 2012, for a detailed description of the survey). MACS1206 is part of the CLASH sample, and it has also been selected for the CLASH-VLT follow-up proposal (Rosati et al., in prep.) and for Subaru imaging for the weak lensing program (Umetsu et al. 2012). We chose this cluster as the test case for our analysis because it is the first cluster for which VLT data reduction is completed, thus we have a wealth of both photometric and spectroscopic information. In this section we describe the dataset at our disposal and the reduction techniques.

2.1. Photometry

We analyzed deep BVRcIcz images obtained with the Suprime-Cam mounted at Subaru telescope and available in the Subaru archive, SMOKA¹. A full description of the observations can be found in Umetsu et al. (2012), while for a detailed explanation of data reduction we refer the reader to Nonino et al. (2009). Here we only provide a brief description. The typical seeing in the final sky-subtracted images varies from 0.58′′ in the z band up to 1.01′′ in the B band with exposure times ranging between 1.6 ks and 3.6 ks with a pixel scale of 0.2 ′′ pixel^-1. The limiting magnitudes are m_B = 26.5, m_V = 26.5, m_Rc = 26.2, m_Ic = 26.0, and m_z = 25.0 mag for a 3σ limiting detection within a 2′′ diameter aperture, see Table 1 for a summary of our photometric dataset.

Sky subtraction and diffuse low-level light-pattern removal are crucial because part of the ICL can be removed in these steps of the data reduction. As described in Nonino et al. (2009), we carefully determine the background by a back-and-forth process. First, we detect sources in a preliminary stacked image, the area covered by each source is enlarged by 20%, and the corresponding segmentation map is used to flag the same pixels in each original image. Flagged pixels in each individual image are replaced by a random value normally distributed with mean and standard deviation obtained by a ~30″ × 30″ box surrounding each pixel, excluding flagged pixel values. Finally, each resulting image is wavelet-transformed, and the background of each image corresponds to the lowest order plane of the wavelet transformation. To ensure that this process does not affect our estimation of the ICL, we use our BCG+ICL map of MACS1206, see Sects. 3.1 and 4, as a control map. Only 0.37% of the BCG+ICL map pixels having a value higher than 3 × σ_sky fall out of the enlarged segmentation map, where σ_sky refers to the σ of the residuals after sky subtraction as estimated in an area free of any source contamination.

None of these pixels is recognized as a source by SExtractor; i.e., these few pixels are randomly distributed, and they most probably represent fluctuations. If we restrict this analysis to a 3′× 2′ area surrounding the BCG, then the percentage of outlier pixels decreases to 0.09%. Thus, the enlarged mask used in the background subtraction process ensures that no pixels associated to the ICL has been oversubtracted. As a consequence, background subtraction does not affect our ICL estimation, and we consider σ_sky as our limit to detect the ICL. As a further check, we applied the SBlimit method (see Sect. 4.3) to the F625W HST stacked image, i.e., the closest HST filter to the Rc Subaru band, and we cross-correlate it with the corresponding Rc band image. This way we can check whether the spatial distribution of the ICL down to different SB levels is the same in both images. According to the cross-correlation analysis, the optimal x,y shift to match the two images is zero for all the SB levels. Given that the HST image has been reduced in an independent way, i.e., using a different background subtraction process, this ensures that we did not remove any real low surface brightness sources during the data reduction.

The stellar point spread functions (PSFs) were measured from a combination of unsaturated stars with S/N ≥ 50 and ellipticity ≤0.1. Here ellipticity is defined as (1 − a)/(1 + a), where a is the source aspect ratio; i.e., an ellipticity of 0.1 corresponds roughly to an aspect ratio of ~0.8. The point sources are detected and modeled using SExtractor and PsfeX softwares (Bertin & Arnouts 1996; Bertin 2011), and their PSF model is derived solely from the robust combination of their resampled input vignettes. In the following analysis, this PSF model is convolved with the best fit model of each galaxy obtained as described in Sect. 3.1.

The B and Rc broad-band filters nicely probe the spectral region across the 4000 Å  break at the cluster redshift, thus the (B − Rc) color is a good indicator of the galaxy average star formation (SF) history, and it can constrain the characteristics of the bulk of its stellar population. We use this color to derive information on the ICL properties.

We obtained magnitudes in each band and the relative colors for all detected sources. These data were used to derive photometric redshifts, z_photo, using a method based on neural networks: the multi layer perceptron with quasi newton algorithm (MLPQNA) (Brescia et al. 2013). This method was calibrated on a subsample of objects with spectroscopic redshifts, and it was applied to the whole dataset with available and reliable BVRcIcz band magnitudes down to m_Rc = 25.0 (see Biviano et al. 2013; Mercurio et al., in prep., for a detailed description on the z_photo estimation). The validation process with spectroscopically measured redshifts makes the estimated z_photo insensitive to photometric systematic errors and more robust than methods based on spectral energy distribution (SED) fitting because the neural network method do not depend on synthesis models or on photometric zero point accuracy. Tests of the MLPQNA based on a combination of parameters from different surveys estimate an excellent accuracy of Δ_zphoto = 0.004 × (1 + z_spec) (Cavuoti et al. 2012; Brescia et al. 2013).

2.2. Spectroscopy

Though our work is based on the imaging data described in the previous section, we also take advantage of the information from the spectroscopic dataset of CLASH/VLT to interpret our results. Here we only give the basic description of this dataset and refer the reader to Rosati et al. (in prep.) and references therein for the details. The CLASH/VLT program is the VLT/VIMOS follow-up of 12/25 CLASH clusters, and it comprises a total of 98 pointings that were obtained in the spectral range of 3700–97 000 Å  using the medium-resolution (MR) and low-resolution (LR) grisms, yielding spectral resolutions of 580 and 180, respectively.

In the case of MACS1206 12 masks (4 MR, 8 LR) were observed for a total exposure time of 10.7 h. Additional spectra were obtained at VLT/FORS2, Magellan telescope, and from literature/archival data (Lamareille et al. 2006; Jones et al. 2004; Ebeling et al. 2009). The final dataset contains 2749 objects with reliable redshift estimates, z_spec with an average error of 75, and 153 km s^-1 for spectra in MR and LR modes.

We measure the main spectral features in the observed spectral range, i.e., Dn(4000), Hδ, [OII], OIII, and Hα. Joining this information to the (B − Rc) color allows us to classify each source according to its stellar population (see Mercurio et al. 2004). In particular, two classes of galaxies will be relevant for discussing our results (see Sect. 5):

• 1.

  Passive galaxies: sources with Dn(4000) > 1.45 and EW(Hδ) < 3.0 Å;

• 2.

  Red Hδ: sources with Dn(4000) > 1.45 and EW(Hδ) > 3.0 Å.

2.3. Cluster membership

We need to distinguish between cluster members and fore- and background sources both when using the ICL detection method for MACS1206 (see Sect. 4) and when determining the cluster total light (see Sect. 4.2). Photometric information is complementary to the spectroscopic one, thus allowing a cluster member association that is complete down to m_Rc = 25.

The cluster membership for each object is assigned according to its spectroscopic redshift, when available, or to its photometric redshift combined with a color–color cut. We refer the reader to Biviano et al. (2013) for a detailed description of membership assignment, but here we summarize the main steps. Briefly, spectroscopic members with 18 ≤ m_R ≤ 23 were defined according to the peak+gap (P+G) method of Fadda et al. (1996). Photometric members were selected among all the sources having a photometric redshift in the range 0.34 ≤ z_spec/photo ≤ 0.54 and satisfying one of the following color–color cut in the (B − V) and (Rc − Ic) diagram:

3. ICL detection

As already mentioned, the ICL consists of the residual light after having removed all the light contribution of galaxies. Ideally, this can be obtained by subtracting each galaxy best-fitting model, choosing among many different light profiles, e.g., de Vaucouleurs, Sérsic (Sérsic 1963, 1968), Exponential disk, and any combination of them. Unfortunately it is not always possible to perfectly fit the galaxies, such that the final residuals are not artifacts due to a bad subtraction. As a consequence, most works favor masking galaxies down to an arbitrary surface brightness level or subtract a direct image via wavelet transformation. In our approach we both subtract the best fit model and mask whenever the fit is not satisfying.

3.1. Method

We developed an automated method based on the software GALAPAGOS (Barden et al. 2012), which makes extensive use of the code GALFIT (Peng et al. 2010). GALAPAGOS detects sources in the target image using SExtractor, estimates sky background, creates postage stamp images for all detected sources, prepares object masks, and finally performs Sérsic fitting with GALFIT. We refer the reader to Barden et al. (2012) for more details, but here we focus only on those steps that are of key importance for our goal. The source detection was performed with a double pass of SExtractor, one for the bright sources and the second for the faintest ones, and then the code recognizes whether to discard or to keep a faint source, depending on its position with respect to the nearest bright source. This minimizes the number of missing or mistaken faint sources.

We set the startup parameter file in order to extract faint source with at least S/N ≥ 1σ_sky. We removed the sky background estimation step since we worked with sky-subtracted images; however, if this step is included, the sky is generally estimated as 0.000 ± 0.001. This support the goodness of our global sky subtraction. The most important step of this code is the postage stamp creation: in this step GALAPAGOS centers the image section on the source of prime interest and optimizes the area to also include the neighboring galaxies. This enables GALFIT to simultaneously fit all sources that contribute to the total light in each section, thus providing a better fit of each contributing source and removing light coming from the outer envelopes of close companions. This cleans the final residual image and ideally provides the light contribution coming only from ICL. It is worth noticing that GALAPAGOS forces GALFIT to fit a single Sérsic model to each source. The initial guess for the Sérsic model parameters corresponds to the SExtractor estimates of x_image, y_image, mag_best, f(flux_radius), and theta_image. In many cases a single Sérsic model is a good approximation, but sometimes it can represent a poor fit, as described in the following. As a last step, the code creates the final output catalog containing both SExtractor and GALFIT information for each source.

At this point we developed an IDL code, GALtoICL, able to go the other way around: from single postage stamps to a final global residual image that we call the BCG+ICL map, where global refers to the original full science image as in contrast to individual postage stamps. From now onward we will use the acronyms GBIma, GBRes, and GBFit when referring to the full original science, residual, and best-fit model images. The code is composed of four main steps:

• 1.

  creation of a GALFIT parameter files fora 1000 × 1000 pixels section of the GBIma;

• 2.

  creation of the GBFit image and the GBRes image;

• 3.

  extraction of those sources with a high percentage of high residuals and manual intervention;

• 4.

  creation of the final BCG+ICL map.

At first all sources are listed according to their χ² and their best fit model parameters are stored. Then a number of GALFIT set-up files containing at most 50 sources each are created until they account for all sources filling the 1000 × 1000 section, i.e., ~1150 × 1150 at MACS1206 redshift. The choice of 50 sources to be modeled in a 1000 × 1000 pixel section corresponds to the best compromise of N_gals and area that GALFIT is able to deal with due to memory issues. All parameters of each source profile are kept fixed since they correspond to their best fit model, and we run GALFIT in model mode, i.e., no fitting, only model image creation based on input parameters. To check whether our conversion from (x,y) postage coordinates to (X,Y) GBIma coordinates is well determined we made some tests allowing (X,Y) to vary within ±2 pixels to account for possible errors in centering the sources. We do not find the need for any (X,Y) marginal correction and thus we rely on our coordinate transformation.

Then, all models in each 1000 × 1000 pixels section are put together to obtain the final GBFit which is then subtracted to the original GBIma to obtain the GBRes image. Bright stars are excluded from the GBFit because they might show strong residuals in case of saturation and they need specific masking. The code allows you to interactively check the GBFit, and the GBRes using DS9, to update the GBFit if necessary, and to run again GALFIT. This is the only step at which manual intervention is possible. The reason for it is well explained in Fig. 1 where we show two examples of GALFIT performances on postage stamps: from left to right we show the original image, its best fit model, and the residuals. Top panels refer to a a clean fit case, while bottom panels show a case with a high percentage of high residuals. Most of the times we get large residuals because a single Sérsic model is not enough to properly describe the galaxies so more components are needed.

To identify the sources with bad fitting residuals in an automated way, we compared the distribution of pixels values in a region of pure sky with that of the residual image. Figure 2 shows these distributions. Those pixels deviating more than 1, 2, 3, 4, and 5σ_sky are flagged, and through SExtractor segmentations maps are connected to the source they belong to. At this point one can choose either to simply mask them or to perform manual fitting, to up-date the model and to re-run GALFIT to create a better GBFit and GBRes images. As a final step, the code allows adding ad-hoc masks to those automatically created to fix bad pixels, i.e., bright saturated stars and spikes. The code is meant to provide BCG+ICL maps; i.e., it does not create the best fit model of the BCG, though one can also choose to obtain only “ICL” maps, i.e., also subtracting the BCG best fit model.

Once the final model achieved the code produces:

• 1.

  final GBFit image;

• 2.

  final GBRes image;

• 3.

  IDs list of deviant sources;

• 4.

  mask images;

• 5.

  final BCG+ICL map with the deviating pixels masked at 1, 2, 3, 4, and 5σ_sky levels.

The whole process, GALAPAGOS+GALtoICL, can be iterated twice in order to identify the bright and well deblended sources at first and then to also model those very faint sources, especially the faint/small satellites of the BCG. To do this, one can choose the “ICL” maps mode and feed again GALAPAGOS with them.

The parameters of the GBFit can be used as a benchmark for other observed bands by running each 50-source GALFIT set-up file in optimize mode, i.e., allowing X,Y,R_e, and Mag to change within a certain range.

Fig. 1 GALFIT residuals examples. From left to right: original image, best fit model, and residuals. Top panels refer to a clean fit case, while bottom panels show a case with a high percentage of high residuals.

Fig. 2 Comparison of pixel-value distribution in the residual image (red dashed line) with that of an empty area (black solid line), i.e., free of source contamination, to identify deviant sources. See text for details.

3.2. Detection efficiency

Before applying our detection method to the real images, we tested its efficiency in detecting faint diffuse-light sources. We generated fake faint sources with different surface brightnesses and randomly introduced them into our real Rc-band images. We also want to determine our ability to deblend and identify these faint sources from close bright companions, thus a low percentage of these fake sources are forced to lie close to a bright one. We then ran our code on these real+simulated images.

The artificial faint sources are modeled as de Vaucouleurs profiles with total magnitude ranging from 21.5 to 24.5 and effective radius varying from 20 to 60 pixels, i.e., ~25–70 at z = 0.44, the cluster MACS1206 redshift. The choice of these parameters translates into surface brightness values ranging between 28 and 32 mag/arcsec² within a 2′′ diameter aperture (28 and 30 mag/arcsec² for the blended sources). In the Local Universe the ICL is usually detected in the V-band, as the light surviving a surface brightness level cut-off, typically μ_V = 26.5,27.5 mag/arcsec² (Feldmeier et al. 2004; Mihos et al. 2005; Krick & Bernstein 2007). To compare our results with these studies, we transformed these V-band SB levels into the corresponding ones at z = 0.44 in the Rc-band; i.e., we added the surface brightness cosmological dimming 2.5·log (1 + z)⁴ and applied the k-correction for different bands. The latter term was determined by running the GALAXEV code on stellar population synthesis models (Bruzual & Charlot 2003) for a solar metallicity with formation redshift zf = 3, a Chabrier initial mass function (IMF) (Chabrier 2003), and accounting for the stellar population evolution. Metallicity and formation redshift values are chosen according to the similarity between typical ICL colors and those of the BCGs (Zibetti et al. 2005; Krick & Bernstein 2007; Pierini et al. 2008; Rudick et al. 2010). The resulting SB levels are μ_Rc(z = 0.44) = 28.87,29.87 mag/arcsec² respectively, while our 1σ_sky level corresponds to μ_1σsky = 30.9 mag/arcsec², thus our Rc-band images are deep enough to detect typical diffuse light sources redshifted to the considered cluster distance.

In Fig. 3 we show our results in terms of SExtractor detection efficiency as a function of the Rc-band surface brightness for both the complete sample of artificial faint sources, i.e., both the randomly positioned ones and those lying close to bright companions, and only the well deblended sources. We set up the SExtractor parameter such that a minimum significant area of five pixels for a 1.5σ detection threshold is requested. As a reference, we show the 1σ_sky surface brightness and the surface brightness limits μ_V(z = 0) = 26.5; 27.5 transformed into the corresponding Rc-band value at z = 0.44.

We note that the detection efficiency for the deblended sample is 100% at SB values well far beyond the lowest μ_V(z = 0) SB level; moreover, the detection efficiency at sky level is almost 50%. If we consider only the range of SB for which we also have blended sources, then the detection efficiency is still more than 70%. These tests ensure that the combination of these deep Subaru images and our detection method is good enough to allow diffuse light source detections for our test-case cluster MACS1206.

The efficency in recovering the initial parameters, such as R_e, Sérsic index, PA, and ellipticity, should also be tested. We used our sample of artificial sources to estimate our ability to recover the original parameter value as a function of the surface brightness as measured within a 2′′ diameter aperture. We split our sample in two subsets, μ_Rc,2″ap ≤ 26.5 and 26.5 <μ_Rc,2″ap< 30.5, to highlight the presence of trends with the SB, if any. Table 2 summarizes our results in terms of the median, low, and high quartiles of the distribution of either the difference or the ratio between the retrieved and the original parameters for each subsample. We do not find any strong trend in the median value as a function of SB, while the errors on the median value tend to increase as we move from high to low surface brightness sources.

This result is in good agreement with Barden et al. (2012) where they used a larger sample of simulated data set-up, i.e., ~10³ more galaxies, in order to achieve enough statistical significance and to test the recoverability with GALAPAGOS of source parameters and its dependence on neighbors. Barden et al. (2012) show that GALAPAGOS has optimal performances for bright galaxies, i.e., μ_input ≤ 22.5, while its efficiency decreases at faint magnitudes, i.e., μ_input> 22.5, and high Sérsic indices, i.e., 2.5 <n< 8.0 (see the left panel of their Fig. 14). Generally speaking, there is no systematic trend or bias for the mean recovered parameter value, while the accuracy gets worse from bright to faint sources. As far as the influence of neighboring galaxies is concerned, Barden et al. (2012) show that GALAPAGOS results do not depend on either the magnitude of or the distance from the next neighbor (see their Fig. 16). Given the agreement on parameters’ retrival tests, we did not repeat this test and rely on their conclusions.

Both the absence of systematic trends and the satisying accuracy level ensure that the recovered GBFit will not be significantly affected by our parameters retrival ability.

Fig. 3 SExtractor detection efficiency as a function of Rc-band surface brightness magnitude. Black diamonds refer to the complete sample of fake faint sources and blue triangles to deblended sources, see text for details. The dotted line corresponds to the sky surface brightness and dashed and dot-dashed lines to the surface brightness limits μV(z = 0) = 26.5; 27.5 after accounting for both surface brightness dimming and k-correction to transform them into Rc-band limits, see text for details.

Fig. 4 Rc band image of the MACS1206 core (left panel), its GBFit image (central panel), and the final BCG+ICL map masked down to 1σsky level.

4. Results: MACS1206 the test case

Our test case cluster, MACS1206, is located at RA = 12^h06^m1228, Dec = –08°48′02.̋4 (J2000), and z = 0.44, and it was originally part of the Most Massive Galaxy Clusters survey (MACS Ebeling et al. 2001). It was codified with morphological class 2, i.e., good optical/X-ray alignment and concentric contours (Ebeling et al. 2010), and this relaxed appearance made it a good target for the CLASH survey. Umetsu et al. (2012) show that there is only a small offset, i.e., 1′′, between the DM peak of mass and the location of the BCG, which also coincides with the X-ray peak emission (Ebeling et al. 2009). The excellent agreement between the mass profile of MACS1206 as derived by the kinematical analysis Biviano et al. (2013) and the lensing analysis Umetsu et al. (2012) is a further indication that this cluster is dynamically relaxed. The relaxed status of the cluster is also confirmed by the absence of a significant level of substructures as found by Lemze et al. (2013).

We notice that despite this general relaxed condition, MACS1206 displays an elongated large-scale structure (LSS) along the NW-SE direction, (Umetsu et al. 2012). This preferred direction is well aligned with the position angle (PA) of the BCG and is also traced by a few infalling groups, as revealed by the dynamical analysis of Girardi et al. (in prep.). The cluster has a velocity dispersion σ_vel = 1087 km s^-1 as estimated by the dynamical analysis of Biviano et al. (2013), from which we also infer a virial mass M₂₀₀ = 1.41 × 10¹⁵M_⊙ that is in good agreement with the results from weak/strong lensing (Umetsu et al. 2012), and it corresponds to R₂₀₀ = 1.98 Mpc.

We ran GALtoICL in the iterated mode on the Rc band image of MACS1206 and used the GBFit as the benchmark model to be adapted for the B-band. After obtaining the first tentative GBFit, we allowed an interactive check and manual intervention in case of large residuals. Specifically, for each galaxy showing a high level of residuals, we proceeded this way: we checked its z_spec, if available and consistent with cluster membership; we performed a detailed manual fit and updated the GBFit; while whenever there was not spectroscopic information we masked at different σ_sky levels. When improving the model by manual fitting, we generally added a second component to the single Sérsic model. Close-enough initial guesses for each component parameter are important for obtaining a reliable fit, so we took advantage of the SExtractor+PsfEx softwares combination that allows spheroid+disk decomposition for each extracted source. The estimated MAG_SPHEROID/DISK, SPHEROID/ DISK_REFF_IMAGE, SPHEROID/DISK_ASPECT_IMAGE, and SPHEROID_SERSICN values are then used as a first guess for GALFIT. Tests on simulated galaxies show that manual intervention reduces by 1.5–2.0 times the number of masked pixels, while providing similar improvement for the residuals in the outermost area of the source segmentation map, i.e., where the signal starts to blur into sky and small differences in the residuals become important for low SB sources.

In Fig. 4 we show the Rc-band image of the MACS1206 core (left panel), its GBFit image (central panel), and the final BCG+ICL map masked down to 1σ_sky level. The galaxy contribution to the light has been removed efficiently, and only 4.8% of the pixels needed to be masked down to 1σ_sky level (only 1.4% when choosing 5σ_sky level).

In the following we report the results we obtained using the masking down to 3σ_sky for the Subaru data that corresponds to μ_Rc ~ 29.3 mag/arcsec^-2 at z = 0.44.

Fig. 5 ICL properties: SB profile and residuals to the best fit (top panels), the ellipticity (central panel), and the PA (bottom panel) as a function of the distance from the center. The dotted and dashed lines in the top panel refer to the SB at 1σsky level and to the best fit model for a de Vaucouleurs profile, respectively. Red squares correspond to those points for which the isophotal analysis did not converge, and the dot-dashed line indicates the psf FWHM limit.

4.1. ICL properties

We performed the classical isophotal analysis of the BCG+ICL using the IRAF² task ellipse. We kept the center position fixed and let the ellipticity and PA vary. Figure 5 shows the SB profile and residuals to the de Vaucouleurs best fit (top panels), the ellipticity (central panel), and the PA (bottom panel) as a function of the distance from the center. We performed a fit of the SB profile with the typical de Vaucouleurs profile. Looking at the residuals, it is clear that the r^1/4 law is a poor representation of the data, and in the outer region of the BCG, R ≥ 40, there is an excess of light with respect to the fit. This excess of light increases as we move farther away from the center and it is the signature of the ICL. At this distance the ellipticity has increased up to ϵ ~ 0.55, while the PA basically has a constant value of PA ~ − 74° (degrees measured counterclokwise from N direction). We notice that these points are located in the regime where the SB reaches the sky level. A close inspection of the BCG+ICL maps reveals an asymmetric elongation of the ICL in the SE direction, thus we suppose that in the SE direction we might also be to detect the ICL at these distances.

Fig. 6 Left panel: zoom of the Rc-band BCG+ICL map of MACS1206 smoothed with a Gaussian kernel of 3 × 3 pixels. We overlaid the slits along the SE (blue) and NW (red) direction from which we extract the SB profiles. The green cross and circle correspond to the location of the BCG and the second brightest galaxy respectively. Right panel: (B − Rc) color map of BCG+ICL. Slits are overlaid as in the left panel.

Fig. 7 Top left panel: SB profile of the Rc-band BCG+ICL map along the SE (blue) and NW (red) directions. Points from the extra slit along the SE direction are highlighted with a yellow circle, while the σsky level is shown by the dotted line. The generic Sérsic (gs), generic Sérsic with high index (gshn), double de Vaucouleurs (dd), de Vaucouleurs plus generic Sérsic (ds), and de Vaucouleurs plus generic Sérsic with low index (dsln) best fit models are shown by the dot-dot-dot-dashed, long-dashed, dot-dashed, short-dashed, and solid lines, respectively. Top right panel: zoomed version of the SB profile in the radial range 20 ≤ R ≤ 100 to highlight the asymmetric radial distribution of the SB profile. Bottom Left panel: fit residuals along each direction for the single component profiles, i.e., the generic Sérsic (circles) and the generic Sérsic with high index (triangles). Bottom right panel: fit residuals along each direction for the double component profiles, i.e., the double de Vaucouleurs (upside-down triangles), de Vaucouleurs plus generic Sérsic (stars), and de Vaucouleurs plus generic Sérsic with low index (squares).

To verify the presence of an asymmetric light distribution, we extracted the SB profile from two slits along the PA: one in the SE direction and the other in the NW direction. In the left panel of Fig. 6, we show a smoothed version of the BCG+ICL map for the Rc band with the slits overlaied. We located two slits along the SE direction: the main one coinciding with the BCG major axis and an extra slit following the ICL elongation toward the second brightest galaxy. We extracted the SB profile from each slit and show it in the top left panel of Fig. 7. Points are color coded according to the slit they belong to. To separate the two slits along the SE direction, we highlight those points obtained from the SE extra slit. The SB profiles along each direction show similar behavior within r ~ 60, while at greater distances the SB profile in the SE direction is systematically above the one in the NW direction. Moreover, at r ≥ 100, the SB profile in the NW direction blurs into the sky regime, while there is still a signal in the SE direction. We also detect a signal from the extra slit even if it is at sky level.

Both SB profiles show excess with respect to the single de Vaucouleurs best fit model, so we tried different models to describe the light profiles: 1) a generic Sérsic profile that either constrains or not the allowed range for the Sérsic index (Oemler 1976; Carter 1977; Schombert 1986; Stott et al. 2011); 2) a double de Vaucouleurs model (Gonzalez et al. 2005); and 3) a composite de Vaucouleurs plus generic Sérsic profile with either free n or within a constrained range of allowed values (Seigar et al. 2007). The generic Sérsic best fit profile (n = 3.16) gives even worse results than the single de Vaucouleurs one, especially in the outer region where the ICL contribution becomes important. The double de Vaucouleurs profile improves the fit even though there is still an excess of light that cannot be fit in the outer region. This light excess can be better appreciated in the zoomed version of the SB profile in the righthand panel of Fig. 7. Color code and line types are the same as in the lefthand panel, but we show only the SB profile at . In contrast, both the composite de Vaucouleurs plus generic Sérsic profiles and the single generic Sérsic profile with 4 <n< 8 manage to also fit the light excess at large distances. The de Vaucouleurs plus generic Sérsic with high index profiles provides the best . The bottom panels show the residulas of single component fitted profiles (left) and composite fitted profiles (right).

In Table 3 we list the best fit parameters for each profile. We notice that both the PA and the ellipticity, ϵ = 1 − q, show a narrow range of values among all the adopted profiles: −70° ≲ PA ≲ − 80° and 0.59 ≲ ϵ ≲ 0.49. This also suggests that in the case of a two-component profile, the BCG and the ICL show good alignment irrespective of the model choice in agreement with the findings of Gonzalez et al. (2005), Zibetti et al. (2005). For a single component fit, the effective radius ranges between ~20 and ~80 , while when we adopt a composite profile, the component associated with the BCG has , whereas the ICL one is less concentrated and has larger effective radius: 37 ≲ r_e,ICL ≲ 175.

As mentioned above, we chose to use the Rc band GBFit as the benchmark model to be adapted for the B-band, this enabled us to create a color BCG+ICL map. We degraded the Rc-band image to the same PSF as that of the B band, i.e., the one with the worst seeing. To transform the PSF of the Rc-band we estimated the kernel function K(r) such that: PSF_{Rc band} (r)∗K(r) = PSF_{B band}(r), where the symbol ∗ denotes a convolution and only unsaturated stars were used. Sky uncertainties are very challenging in creating color maps, in particular at very low SB they can significantly affect the final color even if they are very small, i.e., at μ_V = 28.5 mag/arcsec^-2 an offset of 1σ_sky transforms into an uncertainty of ~0.2 mags in the (B − Rc) color, while at 2 mag brighter the uncertainty is only 0.02. For this reason we rely only on those pixels with μ_V≤29.5.

Fig. 8 (B − Rc) color profile of the BCG+ICL. Points are color coded as in Fig. 7. As a reference we overplot the mean (B − Rc) color of cluster member galaxies within R = 300 (dotted line) and within R500 (dashed line). The shaded area corresponds to the standard deviation of satellite colors within R = 300.

In the righthand panel of Fig. 6, we show the (B − Rc) color map for the BCG+ICL, the color bar shows exactly the color value that ranges from 2.3 in the very core of the BCG, down to 1.5 at distances greater than 50 . As a reference we overlaid the same slits we used in the SB profile analysis. At first glance the map shows a color gradient from redder to bluer colors when moving from the core of the BCG toward the outer regions which are ICL dominated. We quantified this trend extracting the mean color along the slits and in Fig. 8 we show the mean color as a function of the distance from the BCG center in bins of 5 , points are color coded as in the previous plots. The errors correspond to the standard deviation of colors in each bin, as expected in the outer regions the large spread in colors shows the difficulty to retrieve reliable colors at very shallow SB. There is a bluening trend from the BCG center toward outer regions such that the ICL colors tend to be much more like to those of the outer envelope of the BCG rather than its central region. This is consistent with previous results (e.g., Zibetti et al. 2005; Rudick et al. 2010). However the BCG+ICL is reliably detected only out to r = 50 in the B band, i.e., 2σ detection, thus the bluening trend is milder if we consider only the safe detection region. A linear fit to the color profile out to r = 50 returns a slope of −0.16 ± 0.12 in Δ (B − Rc)/Δlog (r), which is compatible with zero gradient or very weak negative gradient. As a reference, we overplot the mean (B − Rc) color of cluster member galaxies within R = 300 and within R₅₀₀. The shaded area corresponds to the standard deviation of satellite colors within R = 300, which is approximately the same for satellites within R₅₀₀. We note that BCG+ICL colors within the safe detection region, i.e., r ~ 50, are in good agreement with those of the satellite galaxies residing in the core of the cluster.

The color profile in the two directions is in good agreement within the error bars, but we note that the innermost point, , in the SE direction tends to be bluer than the corresponding one along the NW direction, though within 1σ. This bluening is confirmed by the presence of [OII] emission in the BCG spectrum obtained by our team with FORS2 as part of the program 090.A-0152(A) see Grillo et al. (in prep.). This [OII] emission line has already been noted by Ebeling et al. (2009), and it was interpreted as evidence in favor of MACS1206 being a CC cluster. However, a careful inspection of HST data reveals the presence of both a compact source and an inner core spiral arm at ~1′′, i.e., ~6 , which are completely blended to the BCG center in the Subaru data thanks to their pixel scale. Both these features are embraced in the spectrum aperture and may be responsible for the [OII] emission. The lefthand panel of Fig. 9 shows the HST F140W image of the BCG center (see Postman et al. 2012; Koekemoer et al. 2011, for the description of HST image observation and data reduction), a small source in the SE direction close to the BCG center is highlighted. In the righthand panel we show the same region but for the F475W filter, whose transmission curve brackets the [OII] emission redshifted at the cluster redshift. In this bluer filter, the blue compact source is visible and separated from the BCG center. This filter also highlights the presence of a sort of spiral arm in the very center of the BCG extending only in the SE direction. Given that this structure is present only in one direction, it is more probable that it is a residual of stripped material.

Our data may suggest that this [OII] emission can be associated to the blue compact source and/or peculiar features blended with the BCG core emission, but we cannot exclude the presence of a moderate/weak CC. Whether MACS1206 is a CC or not is far beyond the purpose of this paper, thus we refer the reader to Appendix A for a brief discussion of this point. For the sake of completeness, we should mention the possibility of the blue compact source being a fore- or background source, while the spiral arm seems connected to the BCG center.

Fig. 9 Left panel: F140W image of the BCG center, the green cross point is located at the BCG center, and the presence of a small source in the SE direction is highlighted by a red arrow. Right panel: same as above but for the F475W filter.

4.2. ICL contribution to the total mass budget

We determined the BCG+ICL fraction as a function of the cluster-centric radius. We extracted the total flux within a set of circular apertures from both the BCG+ICL map and the total members map. To create the total members map, we needed to assign membership to each source in the field of view, and we rely on the cluster membership as described in Sect. 2.3. We masked all the light contribution from fore- and back-ground galaxies down to 1σ_sky, while bright stars were identified using the CLASS_STAR parameter of SExtractor, i.e., CLASS_STAR > 0.98, and we created an ad hoc mask to ensure spike coverage.

In the lefthand panel of Fig. 10, we show the BCG+ICL contribution to the total cluster light within each circular aperture of radius R. Error bars are estimated as in Djorgovski & King (1984): we divide each aperture into eight sections and estimate the total flux in each sector. The error bars represent the rms of total flux in each sector thus taking the possible lumpiness of light distribution in each aperture into account.

We note that at 100 the BCG+ICL contributes more than 50%, while at R ~ 350 it drops down to ~20% of the light within that circular aperture. This BCG+ICL percentage is also confirmed by the analysis of the dark matter profile decomposition performed by Grillo et al. (in prep.) at a similar radial distance.

Fig. 10 Left panel: BCG+ICL contribution to the total cluster light within each circular apertures of radius R as derived from the residual map obtained using the GALtoICL code. Right panel: ICL fraction as a function of the cluster-centric distance for different surface brightness levels and different ICL measurement methods. Empty symbols refer to the SB limit method and filled ones to the GALtoICL code. Circles, triangles, and squares correspond to μRc = 26.5, 27.5, and 28.5 mag/arcsec2 surface brightness levels, respectively. The dotted line at R ~ 1300 kpc indicates R500.

In our approach we extract BCG+ICL maps because it is not trivial to distinguish between the two components so we decided to avoid any a priori separation. However, we can quantify the ICL contribution by combining the de Vaucouleurs + Sérsic profile parameters that best fit the SB profile of the BCG+ICL (see Sect. 4.1) and a proper M/L conversion. Our (B − Rc) color analysis shows that the ICL color tends to be similar to that of the BCG outer envelope; i.e., it can be treated as a red/passive source. To derive the M/L conversion for the ICL, we then determine the best fit of the relation between the stellar masses of red cluster member galaxies; i.e., 2.0 ≤ (B − Rc) ≤ 2.5 and the total Rc magnitude of their best fit model we obtained with GALAPAGOS: (4)where we use stellar masses by Annunziatella et al. (in prep.), i.e., obtained by SED fitting using the MAGPHYS software (da Cunha et al. 2008), based on the 2007 version of the BC03 models (Bruzual & Charlot 2003; Bruzual 2007) with Chabrier IMF (Chabrier 2003) and assuming a set of exponentially declining star formation histories and random bursts superimposed to them. Applying this relation to the total Rc magnitude of the de Vaucouleurs plus generic Sérsic best fit model, we obtain M_ICL = (9.9 ± 3.8) × 10¹¹M_⊙ and M_BCG = (4.0 ± 2.1) × 10¹¹M_⊙.

By summing all the galaxy stellar masses of cluster members down to log (M/M_⊙) = 9.5, i.e., the stellar mass completeness limit corresponding to 23 mag in Rc band (Annunziatella et al., in prep., see text for details), out to R₅₀₀ and that of the BCG as obtained using the above calibration, we obtain the total stellar mass of the cluster, M_{∗ , 500} = (1.7 ± 0.7) × 10¹³M_⊙. Error bars on M_{∗ , 500} were obtained by summing in quadrature the typical galaxy stellar mass error and errors from the standard bootstrap technique. The critical radius R₅₀₀ was determined using the NFW profile for M₂₀₀ = (1.4 ± 0.2) × 10¹⁵M_⊙ and c₂₀₀ = 5.8 ± 1.1 as obtained by the lensing analysis of Umetsu et al. (2012), and we get R₅₀₀ = 1.3 Mpc, which means M₅₀₀ = 1.0 × 10¹⁵M_⊙. The ICL contains 5.9 ± 1.8%,of the stars within R₅₀₀, while the BCG+ICL contribution to M_{∗ , 500} is 8.2 ± 2.5%. As a further check, we estimated the light contained in the de Vaucouleurs + Sérsic best fit model, i.e., in the BCG+ICL components, out to R₅₀₀, and we summed the light of each member galaxy out to R₅₀₀, rather modeling them, to obtain the total cluster light out to R₅₀₀. The corresponding BCG+ICL and ICL fractions are 6.3 ± 0.6% and 4.3 ± 0.2%, respectively. These values are in good agreement with those obtained converting the BCG+ICL total magnitudes into stellar masses within the error bars.

The corresponding contribution of stars, f_∗ to the total mass of the cluster, taking the ICL contribution also into account, is then (M_{∗ , 500} + M_ICL)/M₅₀₀ = 0.0177 ± 0.006. We should also recall that the total galaxy stellar mass within R₅₀₀ is affected by projection effects that tend to increase its value. If we consider a spherical cluster having MACS1206 values for M₂₀₀ and c₂₀₀ and extending out to 3 × R₂₀₀, then the 2D projected mass within R₅₀₀ is 1.56 × M₅₀₀. Considering this projection effect, then M_{∗ , 500, deproj} = 1.18 × 10¹³M_⊙ corresponding to f_{∗ , deproj} = 0.0116 ± 0.006, where we have excluded the BCG from the correction since it lies in the center of the cluster.

4.3. Comparison with the surface brightness method

We now compare these results with those obtained using a different definition of the ICL. We determine the ICL fraction by applying the same approach of many works in the literature (Krick & Bernstein 2007; Burke et al. 2012, and references therein): choosing an arbitrary SB cut-off level below which pixels are masked and counting all the light above this level as the ICL. This ICL definition is a very naive way to separate galaxy light and ICL, but it is the most suitable definition from the operational point of view and for comparison purpose. Moreover, we will be able to explore the advantages and disadvantages of each method and to reveal possible systematics.

We produced ICL maps using SExtractor segmentation maps: we set up the THRESH_TYPE parameter to absolute mode and we chose three different SB cut-off thresholds: 26.5, 27.5, and 28.5 mag/arcsec². This way the sources are extracted only down to each SB level, and the segmentation maps correspond to the galaxy light to be masked. In the ICL maps, those pixels associated to either a source counterpart in the segmentation maps, or stars, fore-, and background galaxies, or sky areas were masked. All the remaining pixels are considered as ICL.

In Fig. 11 we show the Rc-band ICL map down to 26.5, 27.5, and 28.5 mag/arcsec² and the total cluster light map from top left to bottom right. These images show the same asymmetric light distribution along the SE-NW direction in the proximity of the BCG as we found with the GALtoICL code. These images only have a display purpose, to quantify the ICL fraction we sum-up all the flux contained in circular apertures out to R₅₀₀ for each image in Fig. 11.

Fig. 11 Images show only cluster members’ light below given surface brightness levels, which is considered as ICL. The surface brightness limits correspond to μRc = 26.5, 27.5, 28.5 mag/arcsec2 and total cluster light from top left to bottom right. The black circle in the bottom right panel corresponds to R500

In the right panel of Fig. 10 we show the ICL contribution to the total light for each SB level: 26.5, 27.5, and 28.5 mag/arcsec². The fraction of ICL shows a common trend among all SB levels: it has a steep increase from the core out to R ~ 100 kpc where it reaches its maximum, then it shows a plateau. Given that the BCG+ICL fraction as obtained with the GALtoICL code accounts for more than 50% of the light at R ~ 100 kpc and then drops quite rapidly. Then the plateau trend at larger radii can only be justified as a light contribution from the other member galaxies. As a further confirmation, we masked the BCG+ICL map with a circle centered on the BCG and a radius corresponding to the typical distance at which the BCG SB profile reaches 26.5, 27.5, and 28.5 mag/arcsec², i.e., R ~ 15, 30, and 50 kpc. We then extracted the light in the same aperture as before and determined its contribution to the total light. This is shown in the righthand panel of Fig. 10, where different symbols correspond to different SB masking levels. We notice that at large radii, i.e., ~300 , the ICL contribution drops to 10–15% depending on the adopted SB limit. This suggests that most of the ICL is concentrated in the close surroundings of the BCG, while at greater distances the ICL constribution is not significant.

By comparing the ICL fraction as obtained from the GALtoICL code and the SBlimit method we note that even at small distances, i.e., at R ~ 50 kpc, there is a significant difference between them. Moreover, the general trend toward increasing the ICL fraction out to R ~ 80–100 kpc is still present, but then at larger radii, the ICL fraction drops down to a low percentage instead of showing an almost constant value. This reinforces the idea that the SB limit method can be contaminated by the light coming from the outer regions of cluster member galaxies. Despite this, the SB limit method is still the easiest way to compare the observational results to the expected values from the simulations or to other observational studies. Thus we applied the SB limit method for μ_Rc = 28.87,29.87, i.e., the μ_V(z = 0) = 26.5,27.5 SB levels transformed into Rc-band at z = 0.44 (see Sect. 3.2). The corresponding ICL fraction at R₅₀₀ are 12.5 ± 0.6% and 4.7 ± 0.4%, respectively.

These ICL fraction are based on our deepest and best ICL detection filter, the Rc band, but we have multiband imaging of this cluster, so we decided to determine the first ICL SED to measure ICL stellar mass fraction. We used the Rc-band masks for each SB levels as reference masks on the others bands, i.e., B,V, and Ic (having adapted masks to differences in seeing conditions among different bands), we then mask stars, fore-, and background in each band according to their detections down to 1σ_sky level. Finally, we extracted the light that survived to the masking and that is associated to the member galaxies according to SExtractor segmentation maps within R₅₀₀ in each band.

In the top panel of Fig. 12, we show the SED of the total cluster and that of the ICL for different SB limits: μ_Rc = 26.5, 27.5, 28.5 mag/arcsec². We performed a fit to these SEDs using the software MAGPHYS and in the top panel of Fig. 12 we show the SED best-fitting models for the cluster and ICL. In the bottom panel of Fig. 12 we plot the residuals between the observed fluxes in each band and those obtained using the SEDs best-fitting models for each SB level. The ICL mass fraction obtained from the SED fits range between 20% and 55% depending on the chosen SB level and qualitatively in agreement with the SBlimit values. We did not repeat the same exercise for the μ_Rc(z = 0.44) = 28.87,29.87 because the corresponding B and V band masks already cover all the galaxies; i.e., at these SB levels we reach the sky regime.

Fig. 12 Top panel: SED of the total cluster light within R500 (red empty circles) and that of the ICL within R500 for different SB limits: μRc(z = 0.44) = 26.5, 27.5, 28.5 mag/arcsec2 (violet filled circles, blue filled triangles, and cyan filled squares, respectively). Bottom panel: residuals between the observed fluxes in each band and those obtained using the SEDs best-fitting models for each SB level.

5. Discussion

We developed an automated method to create BCG+ICL maps and measured a diffuse intracluster component in MACS1206. We confirm previous findings on general ICL properties: 1) a composite profile fits the data best (Gonzalez et al. 2005; Zibetti et al. 2005; Seigar et al. 2007), though we find that a de Vaucouleurs plus Sérsic profile provides a better fit than a double de Vaucouleurs one; 2) BCG and ICL position angles agreee within a few degrees (Gonzalez et al. 2005; Zibetti et al. 2005) and both agree with the global cluster elongation and its filament (Umetsu et al. 2012; Girardi et al., in prep.); and 3) ICL colors agree with those of the outer envelope of the BCG (Zibetti et al. 2005; Krick et al. 2006; Pierini et al. 2008; Rudick et al. 2010).

Disentangling the BCG component from the ICL is one of hardest tasks when studying the diffuse light, and for this reason we preferred to create BCG+ICL maps. However, to quantify the ICL properties and its contribution to the total cluster light, we separate it from the BCG. We tried different profiles, either single or composite ones by combining the de Vaucouleurs and the Sérsic profiles. Ellipticies and PA show a narrow range of values both in case of a single and composite profiles, while the effective radius show a wider range depending on the adopted profile. In the case of a single component fit, the effective radius ranges between ~20 and ~80 , while when we adopt a composite profile, the component associated with the BCG has 15 ≲ r_e,BCG ≲ 32, whereas the ICL one is less concentrated and it has a larger effective radius: 37 ≲ r_e,ICL ≲ 175.

Ascaso et al. (2011) analyzed a sample of BCGs at a similar redshift, and they fitted them with both a single de Vaucouleurs and a generic Sérsic profile. They find ⟨ r_e,deVauc ⟩ = 19 ± 10 and ⟨ r_e,Sers ⟩ = 23 ± 15. Their mean effective radii are in good agreement with our results if we consider that MACS1206 has a higher X-ray luminosity than that of the Ascaso et al. (2011) sample, i.e., L_{X,0.1 − 2.4 keV} = 24.3 × 10⁴⁴ erg s^-1, and that larger BCGs are located in more massive clusters. Similarly, Stott et al. (2011) find ⟨ r_e,deVauc ⟩ = 27 ± 2 and ⟨ r_e,Sers ⟩ = 57 ± 16 at higher redshift, i.e., z ~ 1. Concerning the effective radius of the outer component for the double de Vaucouleurs fit, we find a small radius when compared to Gonzalez et al. (2005). Their mean effective radii of the ICL component is ~160 kpc though 20% of their sample have r_e,ICL< 50 kpc, thus small ICL effective radii are not ruled out. We should also consider that our double de Vaucouleurs profile is not able to properly fit the outer component (see residuals in bottom panel of Fig. 7), so it might be that we are also underestimating R_e. In contrast, the effective radius of the outer component for the de Vaucouleurs + Sérsic profiles has a higher value, ~140 kpc.

The most peculiar feature of the ICL in MACS1206 is its asymmetric radial distribution: there is an excess of ICL in the SE direction. Peculiar streams of ICL are supposed to last only ~1.5 times their dynamical timescale in the cluster according to simulations (Rudick et al. 2009) because of disruption by cluster tidal field. More generally, the streams found in the cluster core live only τ_{ICL survival} ≤ 1 Gyr owing to the strong tidal fields they are subject to. Thus the galaxy/ies from which this material has been stripped away should have interacted with the BCG no later than a Gyr ago. Moreover, the ICL enhancement along the SE direction extends out to the second brightest galaxy, which is classified as an Hδred galaxy, i.e., poststarburst galaxies (PSBs). The spectral properties of PSB galaxies can only be reproduced by either models of galaxies in a quiescent phase soon after a starburst (τ_PSB ≤ 1.5 Gyr) or by models where a regular star formation has been halted in an abrupt way (Poggianti et al. 1999). Recently, Pracy et al. (2013) have shown that Hδ equivalent width radial profiles in local PSBs can be reproduced by merger simulation at even shorter ages after the peak of the starburst: 0.2–0.75 Gyr. The ICL survival timescale and that of PSBs are in good agreement, thus the ICL stream along the SE direction can be interpreted as the stars stripped from the second brigthest galaxy which has crossed the cluster, sunk to the center, and interacted with the BCG. We note that the second brightest galaxies is aligned with the ICL extra slit PA along the SE direction (see Fig. 6). The dynamical analysis of MACS1206 has highlighted the presence of a preferred direction that is traced by both the passive and Hδred galaxies with PA_{Hδ/ Passive} ~ 110° (measured counter-clockwise from north) (Girardi et al., in prep.). Matching our BCG/ICL PA estimates, we find 101° ≤ PA_BCG/ICL ≤ 109° which is similar to this preferred direction, thus suggesting a further connection between the ICL and the infalling direction of the PSBs population. This scenario is also supported by the presence of an elongated large scale structure (LSS) around the cluster whose major axis runs along the NW-SE direction, 15° ≤ PA_LSS ≤ 30° measured N of W (Umetsu et al. 2012). Matching our PA estimates to the same reference system as Umetsu et al. (2012), we find 11° ≤ PA_BCG/ICL ≤ 19°, depending on the assumed BCG+ICL best fit profile. Thus both the BCG and the ICL are oriented along the same axis as that of the LSS, and this also holds when comparing the ellipticity of the LSS and of the BCG+ICL. As a consequence the BCG of MACS1206 should have experienced a strong interaction that dates back to at least τ_{past merger} ≤ 1.5 Gyr ago. This interaction might also involve the second brightest galaxy and may have occurred along the preferential NW-SE direction.

Both observation and simulations suggest that short-lived major mergers can produce a significant fraction of the ICL (Burke et al. 2012; Burke & Collins 2013; Murante et al. 2007; Laporte et al. 2013; Contini et al. 2014). If we consider the extreme case of the second brightest galaxy merging into the BCG of MACS1206, we can determine the dynamical friction timescale and compare it with the light travel time to z = 0. If the former is shorter than the latter, then we can roughly estimate the second brightest galaxy contribution to the ICL at the end of the merging process. The dynamical friction timescale for a galaxy of mass M_gal at a given initial radius R_in that spirals into the center of the cluster potential well on a circular orbit with velocity V_c is given by Eq. (5) (Binney & Tremaine 1987): (5)where ln(Λ) is the Coulomb logarithm, . In the cluster core, the impact parameter, b_max, is roughly 100 kpc; the typical circular velocity is km s^-1, where we used the velocity dispersion obtained by Biviano et al. (2013); and the second brightest galaxy has M_gal ~ M_{gal, ∗}/f_baryon,gal ~ 10^11.5/ 0.05 ~ 6.3 × 10¹², where we used the galaxy stellar mass obtained by Annunziatella et al. (in prep.) and the typical baryon fraction of early-type galaxies (Hoekstra et al. 2005; Jiang & Kochanek 2007). Thus, ln(Λ) ~ 2.2 and τ_df ~ 2.7 Gyr, given the projected radial distance between the second brightest galaxy and the BCG, R_in ~ 300. Nath (2008) find similar dynamical timescale values for a massive galaxy (M_gal = 3 × 10¹²M_⊙) embedded in a rich cluster (M_cl = 10¹⁵M_⊙) at a similar initial radius. Equation (5) is based on strong approximation, i.e., circular orbit and point-like object. Boylan-Kolchin et al. (2008) take the effect of an extended object with different orbital parameters on the τ_df estimate into account and find that standard approximation tend to shorten the dynamical friction timescale. They also provide a fitting formula to determine the merging timescale due to dynamical friction as a function of both the satellite-to-host halo mass ratio and the satellite orbital properties (see their Eq. (5)). If we consider the host halo as mainly composed of the BCG+ICL, M_host = 10^12.1, and we assume the same baryon fraction as for the second brightest galaxy, then our mass ratio is M_sat/M_host = 10^11.5/ 10^12.1 ~ 0.25. Allowing the initial circularity and the initial orbital energy parameter to vary in the same validity range as Boylan-Kolchin et al. (2008), i.e., 0.33–1.0 and 0.65–1.0, respectively, we obtain 1.0 ≲ τ_merge,df ≲ 6.0 Gyr with a ⟨ τ_merge,df ⟩ ~ 2.6. The light travel time to z = 0 is ~4.6 Gyr, thus there is enough time for the second brightest galaxy to merge into the BCG, if this is the case.

The fraction of ICL coming from galaxies that merged with the BCG ranges between 5% to 30% for the most massive clusters depending on the simulation set-up (Murante et al. 2007; Puchwein et al. 2010; Laporte et al. 2013; Contini et al. 2014). If the second brightest galaxy is going to merge with the BCG, then it will release 1.6 − 9.5 × 10¹⁰M_⊙ to the ICL by z = 0. This corresponds to ~1–10% of the ICL at z = 0.44, and this increase is well within the error bars. A similar consideration can be made in terms of f_ICL, which would become ~5.9–6.4%.

We quantified the mass contribution of the BCG+ICL to the stellar cluster mass within the critical radius R₅₀₀ as ~8%, and this value is in good agreement with the general trend of decreasing BCG+ICL mass (light) fraction with increasing cluster mass (Lin & Mohr 2004; Gonzalez et al. 2007, 2013, G13 hereafter). For comparison purposse in the bottom lefthand panel of Fig. 13, we show BCG+ICL fraction of light (mass) within R₅₀₀ as a function of cluster mass for both MACS1206 and the Gonzalez et al. (2013) cluster sample. Gonzalez et al. (2013) provides BCG+ICL luminosity fractions, while we estimate the mass BCG+ICL fraction. According to Cui et al. (2014), luminosity-weighted and mass-weighted ICL fractions are in good agreement, especially at the high cluster mass end of their sample; i.e., the ratio of luminosity to mass fractions at M₅₀₀ ~ 10¹⁵M_⊙ is consistent with one when AGN feedback is taken into account. We also show the predicted cluster mass M₅₀₀ lower limit for the CLASH sample according to the M − T_X best fit relation of Mahdavi et al. (2013) and to the CLASH cluster selection T_X ≥ 5 keV. We note that the expected cluster mass range covered by the CLASH sample will fill the lack of observational data at the high-mass end, thus allowing this kind of study on a wider cluster mass range and with a well constrained total cluster mass estimate. On top of this, the CLASH/VLT sample will also span a wider range in cosmic time, so we will be able to study the BCG+ICL contribution to the cluster stellar mass disentangling between halo mass and redshift dependences, if any. In the bottom righthand panel of Fig. 13 we show the BCG+ICL fraction as a function of redshift, and the G13 sample is color-coded according to their M₅₀₀ as separated into three bins: M₅₀₀ ≤ 2 × 10¹⁴M_⊙, 2 × 10¹⁴ ≤ M₅₀₀ ≤ 3 × 10¹⁴M_⊙, and M₅₀₀ ≥ 4 × 10¹⁴M_⊙.

Fig. 13 Top panel: stellar baryon fraction as a function of M500 for both MACS1206 and the cluster sample of Gonzalez et al. (2013, G13 hereafter). The (orange square) red triangle refers to the (de-)projected f∗ for MACS1206, while (upside-down gray triangles), open circles refer to the (de-)projected G13 sample. The (green dashed) blue solid line corresponds to the (de-)projected best fit relation from G13, while the dot-dashed line indicates the predicted cluster mass M500 lower limit for the CLASH sample (see text for details). Bottom left panel: BCG+ICL fraction of light/mass within R500 as function of cluster mass for both MACS1206 and the G13 cluster sample, symbols and lines as above. Bottom right panel: BCG+ICL fraction of light/mass within R500 as function of cluster redshift. G13 sample is color-coded according to their M500. Blue, green, and red circles correspond to M500 ≤ 2 × 1014M⊙, 2 × 1014 ≤ M500 ≤ 3 × 1014M⊙, and M500 ≥ 4 × 1014M⊙, respectively.

We notice that the ICL stellar mass (light) of MACS1206 represents ~72 (70)% that of the BCG+ICL assuming our best fit model parameters and the adopted mass to light conversion. Though using a different composite profile, we obtain similar results to Gonzalez et al. (2005) with a high percentage of the light residing in the outer component, the one associated to the ICL. As a consequence, the ICL contribution on small scales is very important, though on larger scales it becomes less significant. This is clearly shown in the righthand panel of Fig. 10 once we adopt a SB threshold on our BCG+ICL maps; on the contrary, applying the same SB limit to the original image shows a plateau of the ICL fraction at large radii. This highlights the systematic error in the ICL contribution estimate depending on the adopted method: light from the outer envelopes of member galaxies can significantly affect the ICL fraction when using the SB limit method. This effect is stronger at lower SB limits, but even at the higher SB limit, the estimated ICL fraction is twice what is obtained with the GALtoICL method. Once again we stress the importance of removing all the light from galaxy members that can affect the real ICL contribution. Unfortunately, the SBlimit method is the best way to compare results among observational works and simulations.

We find good agreement between our ICL fractions at Rc-band SB levels corresponding to μ_V(z  =  0) ≥26.5 mag/arcsec² and those expected from simulations. For a cluster with the same M₅₀₀ as MACS1206, Cui et al. (2014) estimates ICL fraction at R₅₀₀ of 10–20% and 5–10% for μ_V(z = 0) = 26.5 and 27.5, respectively, depending on the adopted simulation, i.e., with either gas cooling, star forming, and supernova feedback or including AGN feedback, thus showing good agreement with our results. Rudick et al. (2011) has simulated clusters with a narrower mass range; but if we consider their most massive cluster B65, M₂₀₀ = 6.5 × 10¹⁴M_⊙, the ICL fraction for μ_V(z = 0) = 26.5 is nearly 12% within 1.5×R₂₀₀ (see left panel of their Fig. 3). Given that they claim only a smaller increase in the ICL fraction within R₅₀₀, these values are in good agreement with our results.

A direct comparison with observational works is less trivial owing to different ICL enclosing radius or lack of cluster total mass information. For instance, Feldmeier et al. (2004) find a ICL fraction of ~10 (2)% above μ_V(z = 0) = 26.5(27.5) mag/arcsec² for a set of clusters located at z ~ 0.17. These values are in good agreement with our ICL fraction of ~12 (4)% at Rc-band SB levels corresponding to μ_V(z = 0) ≥ 26.5(27.5) mag/arcsec², thus suggesting a lack of evolution in the ICL fraction with cosmic time. This result agrees with the absence of strong variation in the amount of ICL between z = 0 and z = 0.8 reported by Guennou et al. (2012) and other authors (Krick & Bernstein 2007). However, we should recall that this comparison is done regardless of the cluster total mass and/or ICL enclosing radius. In contrast, we should mention that most of the simulation studies report a significant increase in the ICL with time. Irrespective of the formation redshift of the ICL, simulations show that roughly 60–80% of the ICL present at z = 0 is built up at z< 1 (Murante et al. 2007; Rudick et al. 2011; Contini et al. 2014). Both the simulation and observation suggest that part of the ICL origins from tidal disruption of intermediate-mass galaxies as they interact with the BCG or the other most massive galaxies in the cluster (Willman et al. 2004; Murante et al. 2007; Coccato et al. 2011; Martel et al. 2012; Giallongo et al. 2014). This scenario is supported by the analysis of environmental dependence of the galaxy mass function of MACS1206 see Annunziatella et al. (in prep.).

We estimate the total star contribution to the baryon fraction and both our projected and de-projected f_∗ are in good agreeement with the results of the recent analysis of Gonzalez et al. (2013) where they also considered the effects of projection. More generally, our values agree with previous studies and the general trend of low f_∗ for the most massive clusters (Andreon 2010; Zhang et al. 2011; Laganá et al. 2011; Lin et al. 2012; Gonzalez et al. 2013). In the top panel of Fig. 13 we show f_∗ as a function of M₅₀₀ for both MACS1206 and the cluster sample of Gonzalez et al. (2013). As a reference, we show the (de-)projected best fit relation from G13 and the predicted cluster mass M₅₀₀ lower limit for the CLASH sample as in the bottom lefthand panel. We note that our estimate of f_∗ is in excellent agreement with the expectation from the best fit relation of G13. Once again we stress that at completion, CLASH/VLT will enlarge the baseline of the f_∗ − M₅₀₀ relation with the advantage of having well constrained cluster total mass.

Adding the gas fraction f_g = 0.144 ± 0.025 as estimated by Ettori et al. (2009) to the stellar component, we obtain the total baryon fraction f_b = 0.156 ± 0.026, to be compared with f_b = 0.167 (0.154), as expected from WMAP7 (Planck) results (Planck Collaboration XVI 2014; Komatsu et al. 2011). The comparison with Planck results is less straightforward owing to different cosmological parameters that have a strong impact, as shown by Gonzalez et al. (2013). Our total baryon fraction is 7% below the expected value but well within 1σ. Generally speaking, this result agrees with the trend toward increasing (decreasing) f_gas (f_∗) with cluster total mass, thus supporting the idea of a less efficient star formation at the high end of the cluster mass function (Andreon 2010; Zhang et al. 2011; Laganá et al. 2011; Lin et al. 2012; Gonzalez et al. 2013, and references therein).

6. Summary and conclusions

We have developed an automated method of extracting BCG+ICL light maps in a refined way: GALtoICL. By applying this technique to MACS1206:

• 1.

  We highlighted the presence of an extra component, i.e., theICL, when studying the SB profile of the BCG.This component appears to be asymmetric in radial distribution,and we interpret it as evidence of a past merger. We linked the ICLproperties to those of the cluster substructures, and in this way wereconstructed the most recent cluster assembly history.

• 2.

  We estimated the BCG+ICL mass fraction and the (de-) projected f_∗ of MACS1206 to be in good agreement with recent results that suggest a lowering in star formation efficency at higher cluster masses.

• 3.

  We estimated the sole ICL contribution with two different methods, GALtoICL and the SBlimit methods, and compared their results. The SBlimit method provides ICL fractions that are systematically larger than those obtained with the GALtoICL method due to the light contamination of member galaxies other than the BCG. The GALtoICL method removes this contamination by fitting simultaneously galaxies, thus providing safe ICL detection, and it also highlights the presence of features such as plumes in the ICL. As a disadvantage, the GALtoICL method is much more time-consuming than simpler methods, such as the SB limit definition, and it can only be applied to a small field of view.

• 4.

  Based on the SBlimit method, we obtained the first tentative ICL SED. The ICL mass fraction that we obtained by the SED fitting is in qualitativly good agreement with those simply obtained by fluxes in the single reference broadband filter Rc.

The high-quality dataset, the new refined ICL detection method, and the comparison of different ICL detection methods are the most striking novelties of this work. Deep multiband photometry allowed us to securely detect the ICL at a relatively high redshift, z = 0.44, while the spectroscopic information allowed us to select cluster members, determine their masses down to log (M/M_⊙) = 9.5, and thus obtain an accurate estimate of the cluster stellar mass, BCG+ICL stellar mass, and f_∗. The wide spectroscopic dataset also permitted the ICL properties to be associated to the dynamical analysis of MACS1206 and thus reconstruct its assembly history. While a single data point cannot give statistical relevance to our results or allow to draw strong conclusions, at completion the CLASH/VLT survey will provide a high-quality dataset over a wide redshift range, thus enabling us to constrain both the role of the ICL in the baryon budget and the f_∗ − M₅₀₀ relation.

This work also highlighted the importance of a common definition of ICL to allow comparison among both observational and numerical works. Simple ICL definitions, such as the SBlimit method, might be easier to compare, but they do not retrieve the real ICL properties because of contamination effects.

Acknowledgments

We thank the anonymous referee for constructive comments that helped us to improve the manuscript. V.P. is grateful to Monaco, P., Murante, G., and De Grandi, S. for useful discussion and comments. V.P. acknowledges the grant “Cofinanziamento di Ateneo 2010” and financial support from PRIN-INAF2010 and MIUR PRIN2010-2011 (J91J12000450001). W.C. acknowledges a fellowship from the European Commission’s Framework Programme 7, through the Marie Curie Initial Training Network CosmoComp (PITN-GA-2009-238356), supports from ARC DP130100117 and from the Survey Simulation Pipeline (SSimPL; http://ssimpl-universe.tk/). A.F. acknowledges the support by INAF through VIPERS grants PRIN 2008 and PRIN 2010. The work of LAM was carried out at Jet Propulsion Laboratory, California Institute of Technology, under a contract with NASA. Support for A.Z. is provided by NASA through Hubble Fellowship grant #HST-HF-51334.01-A awarded by STScI.

References

Appendix A: BCG [OII] emission line

Our team obtained a medium-resolution spectrum of the BCG with FORS2 as part of the program 090.A-0152(A) see Grillo et al. (in prep.). We measured the [OII] equivalent width (EW) from an aperture of ~1.5′′, i.e., 9 diameter, around the peak emission of the BCG flux-calibrated spectrum: EW_OII = −4.9 ± 3.2 Å. This corresponds to L_[OII] = 7.4 ± 4.8 × 10⁴⁰ erg s^-1, having multiplied the EW by the flux density of the best-fitting SED at 3727 Å. The level of our [OII] emission line detection is very low, in contrast to what is expected for a strong-to-moderate cool core (CC) (Crawford et al. 1995) and in agreement with normal BCG showing no/low [OII] emission (Samuele et al. 2011). This [OII] emission line has already been noted by Ebeling et al. (2009), and it was interpreted as evidence in favor of MACS1206 being a CC cluster. Ebeling et al. (2009) also note that the [OII] emission was at a much lower level than typically observed in large CC clusters, thus flagging MACS1206 as a moderate CC cluster. Using a different parameter, Baldi et al. (2012) also classify MACS1206 as a CC cluster even if the temperature profile is approximately constant around kT ~ 10 keV. This kind of temperature is very high as compared to typical CC central temperatures, i.e., 3–4 keV (Finoguenov et al. 2001), and it also has too low a central metallicity, i.e., 0.25 (Cavagnolo et al. 2009, see also the ACCEPT web site³) with respect to typical CC. Cavagnolo et al. (2009) also estimated the central cooling time, τ₀ ~ 1 Gyr, and the central entropy, K₀ ~ 70 keV cm², of MACS1206. These values are borderline between the absence of CC and the presence of a weak CC according to the multiparameter analysis of Hudson et al. (2010).

All Tables

All Figures

	Fig. 1 GALFIT residuals examples. From left to right: original image, best fit model, and residuals. Top panels refer to a clean fit case, while bottom panels show a case with a high percentage of high residuals.
In the text

	Fig. 2 Comparison of pixel-value distribution in the residual image (red dashed line) with that of an empty area (black solid line), i.e., free of source contamination, to identify deviant sources. See text for details.
In the text

Fig. 3 SExtractor detection efficiency as a function of Rc-band surface brightness magnitude. Black diamonds refer to the complete sample of fake faint sources and blue triangles to deblended sources, see text for details. The dotted line corresponds to the sky surface brightness and dashed and dot-dashed lines to the surface brightness limits μV(z = 0) = 26.5; 27.5 after accounting for both surface brightness dimming and k-correction to transform them into Rc-band limits, see text for details.

In the text

	Fig. 4 Rc band image of the MACS1206 core (left panel), its GBFit image (central panel), and the final BCG+ICL map masked down to 1σsky level.
In the text

Fig. 5 ICL properties: SB profile and residuals to the best fit (top panels), the ellipticity (central panel), and the PA (bottom panel) as a function of the distance from the center. The dotted and dashed lines in the top panel refer to the SB at 1σsky level and to the best fit model for a de Vaucouleurs profile, respectively. Red squares correspond to those points for which the isophotal analysis did not converge, and the dot-dashed line indicates the psf FWHM limit.

In the text

Fig. 6 Left panel: zoom of the Rc-band BCG+ICL map of MACS1206 smoothed with a Gaussian kernel of 3 × 3 pixels. We overlaid the slits along the SE (blue) and NW (red) direction from which we extract the SB profiles. The green cross and circle correspond to the location of the BCG and the second brightest galaxy respectively. Right panel: (B − Rc) color map of BCG+ICL. Slits are overlaid as in the left panel.

In the text

Fig. 7 Top left panel: SB profile of the Rc-band BCG+ICL map along the SE (blue) and NW (red) directions. Points from the extra slit along the SE direction are highlighted with a yellow circle, while the σsky level is shown by the dotted line. The generic Sérsic (gs), generic Sérsic with high index (gshn), double de Vaucouleurs (dd), de Vaucouleurs plus generic Sérsic (ds), and de Vaucouleurs plus generic Sérsic with low index (dsln) best fit models are shown by the dot-dot-dot-dashed, long-dashed, dot-dashed, short-dashed, and solid lines, respectively. Top right panel: zoomed version of the SB profile in the radial range 20 ≤ R ≤ 100 to highlight the asymmetric radial distribution of the SB profile. Bottom Left panel: fit residuals along each direction for the single component profiles, i.e., the generic Sérsic (circles) and the generic Sérsic with high index (triangles). Bottom right panel: fit residuals along each direction for the double component profiles, i.e., the double de Vaucouleurs (upside-down triangles), de Vaucouleurs plus generic Sérsic (stars), and de Vaucouleurs plus generic Sérsic with low index (squares).

In the text

	Fig. 8 (B − Rc) color profile of the BCG+ICL. Points are color coded as in Fig. 7. As a reference we overplot the mean (B − Rc) color of cluster member galaxies within R = 300 (dotted line) and within R500 (dashed line). The shaded area corresponds to the standard deviation of satellite colors within R = 300.
In the text

	Fig. 9 Left panel: F140W image of the BCG center, the green cross point is located at the BCG center, and the presence of a small source in the SE direction is highlighted by a red arrow. Right panel: same as above but for the F475W filter.
In the text

Fig. 10 Left panel: BCG+ICL contribution to the total cluster light within each circular apertures of radius R as derived from the residual map obtained using the GALtoICL code. Right panel: ICL fraction as a function of the cluster-centric distance for different surface brightness levels and different ICL measurement methods. Empty symbols refer to the SB limit method and filled ones to the GALtoICL code. Circles, triangles, and squares correspond to μRc = 26.5, 27.5, and 28.5 mag/arcsec2 surface brightness levels, respectively. The dotted line at R ~ 1300 kpc indicates R500.

In the text

	Fig. 11 Images show only cluster members’ light below given surface brightness levels, which is considered as ICL. The surface brightness limits correspond to μRc = 26.5, 27.5, 28.5 mag/arcsec2 and total cluster light from top left to bottom right. The black circle in the bottom right panel corresponds to R500
In the text

Fig. 12 Top panel: SED of the total cluster light within R500 (red empty circles) and that of the ICL within R500 for different SB limits: μRc(z = 0.44) = 26.5, 27.5, 28.5 mag/arcsec2 (violet filled circles, blue filled triangles, and cyan filled squares, respectively). Bottom panel: residuals between the observed fluxes in each band and those obtained using the SEDs best-fitting models for each SB level.

In the text

Fig. 13 Top panel: stellar baryon fraction as a function of M500 for both MACS1206 and the cluster sample of Gonzalez et al. (2013, G13 hereafter). The (orange square) red triangle refers to the (de-)projected f∗ for MACS1206, while (upside-down gray triangles), open circles refer to the (de-)projected G13 sample. The (green dashed) blue solid line corresponds to the (de-)projected best fit relation from G13, while the dot-dashed line indicates the predicted cluster mass M500 lower limit for the CLASH sample (see text for details). Bottom left panel: BCG+ICL fraction of light/mass within R500 as function of cluster mass for both MACS1206 and the G13 cluster sample, symbols and lines as above. Bottom right panel: BCG+ICL fraction of light/mass within R500 as function of cluster redshift. G13 sample is color-coded according to their M500. Blue, green, and red circles correspond to M500 ≤ 2 × 1014M⊙, 2 × 1014 ≤ M500 ≤ 3 × 1014M⊙, and M500 ≥ 4 × 1014M⊙, respectively.

In the text