Dark Energy Survey Year 1 results: Detection of Intra-cluster Light at Redshift $\sim$ 0.25

TL;DR

This paper reports the first detection and analysis of intra-cluster light in galaxy clusters at redshift ~0.25 using DES data, revealing its properties, distribution, and relation to cluster mass.

Contribution

It introduces new methods for masking and stacking to measure faint ICL, and demonstrates its self-similarity and correlation with cluster mass distribution.

Findings

ICL extends to at least 1 Mpc from cluster centers.

ICL accounts for about 44% of total stellar luminosity within 1 Mpc.

ICL brightness correlates with cluster mass distribution.

Abstract

Using data collected by the Dark Energy Survey (DES), we report the detection of intracluster light (ICL) with $\sim300$ galaxy clusters in the redshift range of 0.2-0.3. We design methods to mask detected galaxies and stars in the images and stack the cluster light profiles, while accounting for several systematic effects (sky subtraction, instrumental point-spread function, cluster selection effects and residual light in the ICL raw detection from background and cluster galaxies). The methods allow us to acquire high signal-to-noise measurements of the ICL and central galaxies (CGs), which we separate with radial cuts. The ICL appears as faint and diffuse light extending to at least 1 Mpc from the cluster center, reaching a surface brightness level of 30 mag arcsec$^{-2}$. The ICL and the cluster CG contribute to $44\%\pm17$\% of the total cluster stellar luminosity within 1 Mpc. The ICL color is overall consistent with that of the cluster red sequence galaxies, but displays the trend of becoming bluer with increasing radius. The ICL demonstrates an interesting self-similarity feature -- for clusters in different richness ranges, their ICL radial profiles are similar after scaling with cluster $R_\mathrm{200m}$, and the ICL brightness appears to be a good tracer of the cluster radial mass distribution. These analyses are based on the DES redMaPPer cluster sample identified in the first year of observations.

Tables

Table 1. Table 1 Fitted parameters of the ICL r − limit-from 𝑟 r- band radial profile to Sersic models, I ​ ( r ) = I ​ e × exp ​ ( − b n ​ ( r / R ​ e ) 1 / n ) 𝐼 𝑟 𝐼 𝑒 exp subscript 𝑏 𝑛 superscript 𝑟 𝑅 𝑒 1 𝑛 I(r)=Ie\times\mathrm{exp}(-b_{n}(r/Re)^{1/n}) .

	$I{e}^{\ast }$	$n$	$Re$
	(flux/arcsec2)
First Sersic Component	$9830\pm 162$	$1.34\pm 0.023$	$9.13\pm 0.24$ kpc
Second Sersic Component	$8846\pm 1046$	$3.07\pm 0.08$	$52.1\pm 2.2$ kpc
Third Sersic Component	$9.1\pm 3.3$	$2.1\pm 0.4$	$2.6\pm 0.7$ Mpc

Equations

$$f(r)=2\pi {\int }_0^r{r}^{\prime }l(r^{\prime })\mathrm{d}{r}^{\prime }.$$

$$\nabla (g-r)\equiv \frac{\mathrm{d}(g-r)}{\mathrm{d}\mathrm{Log}(r)}=-0.152\pm 0.027,$$

$$\begin{array}{rl}\xi (r|M,c)& =\max ({\xi }_{NFW}(r|M,c),{\xi }_{2h}(r|M)),\\ \Sigma (r|M,c,z)& ={\int }_{-\infty }^{\infty }d\chi {\rho }_m(z)\xi (\sqrt{r^2+{\chi }^2}|M,c).\end{array}$$

Full text

1. Introduction

The central galaxy (CG) of a galaxy cluster tends to be surrounded by an extended light envelope (Zwicky 1951, 1952; Matthews et al. 1964; Morgan & Lesh 1965). Studies indicate that this light envelope extends to hundreds of kiloparsecs and sometimes encloses several galaxies, especially if the cluster is experiencing a merging process (see reviews in Mihos 2004; Lauer et al. 2014; Mihos 2015). Given its diffuse nature (Schombert 1986, 1987, 1988) and the fact that it may enclose multiple galaxies, it seems more reasonable to consider this envelope as a genuine component of galaxy clusters. The diffuse light envelope is thus frequently referred to as the intracluster light (ICL).

Despite the conceptual difference between CGs and the ICL, it is almost impossible to observationally separate them, because the outskirts of CGs naturally blend into the ICL. Outside the outskirts of CGs, the ICL is extremely faint, which poses significant challenges for separating it from CGs, and for further characterizing its distribution and properties. In this paper, we do not attempt to dissect the ICL and CGs. Instead, we follow a similar convention to Pillepich et al. (2018), wherein we simultaneously characterize CGs and the ICL. We use ”ICL” to qualitatively refer to the diffuse light outside a few kpcs of the CG centers, and when making quantitative assessments, ”ICL” and ”CG” are quantitatively distinguished by radial cuts.

Studies, especially of the halo occupation distribution, indicate that the ICL potentially makes up a significant fraction of the cluster stellar light. Through modeling the growth histories of cluster member galaxies or cluster CGs, a few analyses have found that much of the stellar mass accreted by a galaxy cluster is missing from the luminosity or stellar mass sum of cluster galaxies or CGs, and the ICL has been considered the key to this missing mass (Lin & Mohr 2004; Monaco et al. 2006; Burke et al. 2012; Contini et al. 2014; Conroy et al. 2007). According to the extrapolation in Behroozi et al. (2013), the ICL may contain $\sim$ 2.5-5 times the stellar mass of the CG for a $1\times 10^{14}M_{\odot}$ cluster. Although a quantitative conclusion about ICL stellar content is yet to be reached in simulations and observational studies, many estimate that the ICL and CGs consist of up to $\sim$ 10 % to 50 % of the total cluster stellar content (e.g., Zibetti et al. 2005; Gonzalez et al. 2007; Burke et al. 2015; Montes & Trujillo 2018).

The origin and distribution of the ICL have been explored theoretically or through simulation studies. Most of these explain the origin of the ICL as stars that are originally formed inside galaxies and became dispersed into the intracluster space during galaxy interactions (see reviews in Contini et al. 2014; DeMaio et al. 2018). Tidal stripping (Gallagher & Ostriker 1972), galaxy disruption (Guo et al. 2011) and relaxation of galaxy mergers (Murante et al. 2007) are possible astrophysical processes that produce the ICL in this way. Surprisingly, in-situ star formation also appears to be a viable contributor to the ICL (Puchwein et al. 2010; Tonnesen & Bryan 2012). The ICL likely has formed through multiple channels, but the relative contributions from the channels are still being explored. Depending on the main formation mechanism, simulated ICL exhibits different color and spatial distributions (Contini et al. 2018b), and the total amount of the ICL stellar mass varies (Contini et al. 2018a), which provide clues for testing ICL formation hypotheses with observations.

Observational studies of the ICL have been performed with both Hubble and ground-based telescopes. These include targeted imaging observations of individual (e.g., Krick et al. 2006; Toledo et al. 2011; Giallongo et al. 2014; Montes & Trujillo 2014; Presotto et al. 2014) or a sample of clusters (e.g., Krick & Bernstein 2007; Morishita et al. 2017; DeMaio et al. 2018; Montes & Trujillo 2018), or statistical studies of a few hundreds of clusters based on wide field survey data (Zibetti et al. 2005). Many works focused on nearby (e.g., Mihos et al. 2017) or intermediate-redshift clusters, but high-redshift ICL observations to $z\sim 1.0$ (Burke et al. 2012, 2015; Ko & Jee 2018) have been reported. Advances in integral-field spectroscopy also bring forward spectroscopic studies of the ICL stellar population (Adami et al. 2016; Edwards et al. 2016; Gu et al. 2018; Johnston et al. 2018). The ICL has been found to contain an old stellar population and extends to several hundreds of parsecs from cluster centers, although cluster-to-cluster variation (e.g., Krick & Bernstein 2007), cluster dynamic state (Jiménez-Teja et al. 2018), redshift evolution in addition to methodological differences and sometimes differences in ICL definitions (Cui et al. 2014; Jiménez-Teja & Dupke 2016; Tang et al. 2018) may cause different conclusions.

Due to its low surface brightness and diffuseness, studying the ICL presents significant challenges. Simulation ICL studies generally require high resolution (Contini et al. 2014) as well as proper baryonic physics such as active galactic nucleus (AGN) feedback (Martizzi et al. 2016; Pillepich et al. 2018) to trace the dispersion of stellar particles. Observational studies using the Hubble Space Telescope (HST) tend to be pointed observations and benefit from the high resolution and low sky background level of the images but are limited to small cluster samples. ICL studies with ground-based telescopes benefit from the richer observing resources yet face more challenges in data processing. Studies of individual galaxy clusters require long exposure time, often dozens of hours even with large-aperture telescopes. Analyses of the data then require a meticulous account of systematic effects such as flat-fielding and sky background subtraction, and the telescope and imaging camera point spread function (PSF) is a confusing effect for interpreting ICL detection results at large cluster radius (Bernstein 2007; Sandin 2014; Tang et al. 2018).

Studying the statistical distribution of the ICL in an ensemble of clusters is one approach toward easing systematic effects from ground-based observations. A particularly exciting method is to combine the images of dozens to several hundreds of galaxy clusters by aligning the centers of the CGs and measuring the ICL in these combined images, also known as “stacking”. With this approach, Zibetti et al. (2005) analyzed $\sim 600$ clusters with Sloan Digital Sky Survey (SDSS) data and detected the ICL at a surface brightness level of 32 mag/arcsec2. Krick et al. (2006) adopted a similar ICL analysis method, and Tal & van Dokkum (2011) and Wang et al. (2018) applied the stacking idea to analyze the light envelopes of luminous red galaxies and the CGs of galaxy groups. Notably, the stacking method is less affected by, but not completely immune to, systematic effects. Issues such as sky background subtraction and the PSF still plague the method and need to be carefully evaluated.

In this paper, we use a “stacking”, i.e., an averaging approach in which we analyze hundreds of clusters for ICL detection and characterization. We use observational data sets from the Dark Energy Survey (DES), which is in the process of collecting imaging data of an unprecedented combination of large area (5000 deg2 of the sky) and depth. DES expects to find 30,000 galaxy clusters using their population of member galaxies. Although this paper focuses on ICL detection in clusters of redshift $0.2-0.3$ identified in the first year of DES observations (DES Y1), the final DES data set will eventually allow us to apply our methods to a higher redshift range and possibly characterize the redshift evolution of ICL.

We develop a “stacking” approach in that we study the ICL by stacking the light profiles derived from individual cluster images, rather than stacking the images of the clusters as in previous works (Zibetti et al. 2005; Krick et al. 2006; Tal & van Dokkum 2011). Our method achieves results equivalent to the latter method upon our initial method testing, but it is computationally more efficient and versatile. The data sets and methods used in this paper are described in Sections 2 and 3 respectively. Because PSF is a possible effect that influences ICL interpretation, we derive the DES Dark Energy Camera PSFs in Section 4. Section 5 and Section 6 present our results on ICL detection and ICL properties. In Section 7, we discuss how systematic effects, such as sky background subtraction and the PSF, influence the interpretation of the results. Section 8 discusses and summarizes the analyses. Throughout the paper, we assume a flat $\Lambda$CDM cosmology with $h=0.7$.

2. Data

2.1. The redMaPPer Cluster Catalogs

We study the ICL in a galaxy cluster sample identified with the redMaPPer algorithm (Rykoff et al. 2014). redMaPPer searches for galaxy clusters by considering galaxy colors, luminosities, and spatial distributions. A notable feature of the algorithm is that it looks for cluster galaxies that match the colors of the red sequence galaxies in spectroscopically confirmed galaxy clusters. The result is a highly complete and pure cluster sample with precise photometric redshifts and a cluster richness estimator that scales well with cluster mass (Rykoff et al. 2012; Rozo & Rykoff 2014; Rozo et al. 2015a, b; Melchior et al. 2017; McClintock et al. 2019). A detailed description of the algorithm can be found in Rykoff et al. (2016). The final products of the redMaPPer algorithm includes a cluster catalog, a cluster member galaxy catalog, and a random point catalog that records the sky coverage of the cluster-finding algorithm. Both the cluster catalog and the random point catalog are used in this paper.

Specifically, we select a subset of $0.2<z<0.3$ redMaPPer clusters discovered in DES Year 1 coadd observations (McClintock et al. 2019). Because of a slightly better performance in identifying cluster CGs, we use a DES internal version of the redMaPPer sample based on DESDM coadd catalogs, rather than the nominal redMaPPer sample based on multiobject fitting photometry (Drlica-Wagner et al. 2018). We eliminate galaxy clusters that have flags of bright stars, bright galaxies, adjacency to the Large Magellanic Cloud, etc. (bad region mask $>$ 2 in Drlica-Wagner et al. 2018) within 526″(2000 DES coadd image pixels) of the redMaPPer centers. In particular, these eliminations reduce undesirable imaging features such as streaks of saturated stars or improper sky estimation around nearby galaxies, which would lower the overall signal-to-noise ratio (S/N) of ICL detection. A DES object depth cut is also applied, with details described in Section  2.3, which requires object detection to be highly complete in the cluster regions. These criteria remove $\sim$ 2/3 of the redMaPPer cluster sample. In total, we analyze $\sim$ 280 galaxy clusters with a mean mass of $2.50\times 10^{14}\mathrm{M_{\odot}}$, the redshift and richness distribution of which are shown in Figure 1.

We use the redMaPPer random point catalog, which is a sample of random locations that traces the sky coverage of the redMaPPer cluster sample and matches the sky coordinates of the redMaPPer cluster search area. The main purpose of using it is to estimate light profiles of random fields in the redMaPPer footprint. As described in Section 3, the light profiles of the random fields are subtracted from the cluster light profiles to eliminate sky-subtraction residual. For this purpose, we draw a sample of 2000 random points, which is $\sim 7$ times the size of the cluster catalog to avoid significantly increasing profile uncertainties. We also test the ICL detection procedures with the random points as described in Section 3.3. For that purpose, we draw a second set of random points of the same size as the cluster sample.

2.2. DES Finalcut Images

This analysis makes use of fully processed single-exposure images, i.e., Finalcut images, from the official DES image-processing pipeline (Morganson et al. 2018). The corresponding pipeline handles bias and flat corrections, astrometry and photometry calibrations, as well as the masking of bad pixels and artifacts including cosmic rays and satellite trails. Because of the significantly improved sky-subtraction procedure applied to DES images collected in the first three years (Y3), we make use of the Finalcut images from the DES Y3 data processing campaign, yet limited to the region of Y1 redMaPPer identifications .

During this Y3 data processing campaign, the imaging sky background is estimated and subtracted from the Finalcut images using a Principal Component Analysis (PCA) method developed in Bernstein et al. (2017). The method evaluates the sky level across the whole focal plane image of 3 deg2. The full exposure sky background is described as the linear combination of a small set (three or four) of fixed patterns, and the coefficient for each is determined independently for each exposure. This PCA method is designed to capture the possible realizations of signals produced by the uniform zodiac and atmospheric background with gentle gradients across the image, but it does not attempt to remove light coming from individual objects. Because the sky background is evaluated over the 3 deg2 focal plane, it is unlikely to be severely influenced by bright stars, ghosts, dust reflection, and ICL, as long as they do not affect a significant fraction of the focal plane. We further investigate the effect of sky subtraction on ICL detection in Section 7.

2.3. The Object Catalog

To enable the detection of ICL, we mask regions associated with objects already detected by the DES coaddition pipeline.

The object-detection scheme from the DES coadd image is described in Drlica-Wagner et al. (2018); Abbott et al. (2018); Morganson et al. (2018). In short, the SExtractor software (Bertin & Arnouts 1996) is employed to detect objects in the DES coadd images, yielding coordinates as well as photometry measurements for the detected objects. Of particular interest to this analysis, SExtractor evaluates a Kron magnitude for each object, which is the flux enclosed within 2.5 Kron radii (Kron 1980), denoted as MAG_AUTO. The coadd catalog generated from DES Y1 data reaches a 10 $\sigma$ magnitude limit of $\sim$ 22.5 mag in the $i-$band. Although not relevant to this paper, multiepoch, multiobject fitting photometry is also acquired for the detected objects using the ngmix code (Sheldon 2014; Jarvis et al. 2016).

We select a subset of clusters around which object detection is highly complete and hence the objects can be correctly masked. The requirement is that the $i$-band completeness limit (see definition in Zhang et al. 2017) within 526″(2000 DES coadd image pixels) of the redMaPPer CGs reaches 22.4 mag. Most of the redMaPPer clusters identified in Y1 data satisfy this requirement.

3. Methods

The centerpiece of this analysis is the measurement of ICL profiles. The derivation of ICL light profiles is based on DES images and consists of two steps, image processing and profile estimation. In this section, we provide the details of the method components and present a test to demonstrate the effectiveness of the methods in determining an accurate ICL light profile. The general flow of the methods is illustrated in Figure 2.

3.1. Image Processing

The ICL analysis starts with processing DES single-exposure images. For each of the clusters, we download DES Finalcut images and the corresponding weight maps overlapping a circular region of 0°.15 radius centered on the redMaPPer central galaxy. These images are projected and averaged with weights to create one coadd image for each cluster, using the SWARP software (Bertin et al. 2002). Specifically, the images are projected with a TAN projection system (see a review in Calabretta & Greisen 2002), and the SWARP sky background function is turned off. These coadded images for each of the clusters are the input to computing the profiles of ICL. Notably, we do not use the coadded images from the official DES coaddition pipeline, which has incorporated a sky-subtraction process by the SExtractor software and thus complicates ICL detection.

To improve the detection of faint diffuse light, we mask detected objects above a magnitude limit. These objects are selected from the DES coadd catalog (Section 2.3) to be (1) above the 5 $\sigma$ magnitude limit in the $i$-band, i.e., $\texttt{MAGERR\_AUTO\_I}<0.218$, and (2) i.e., brighter than 22.4 mag in the $i-$band, i.e., $\texttt{MAG\_AUTO\_I}<22.4$. According to the study of the relation between DES catalog depth and completeness in Zhang et al. (2017), objects of $i-$band $\texttt{MAG\_AUTO}<22.4$ mag and $magerr\_auto<0.218$ (flux measurement at 5$\sigma$) reach a $\sim$ 99.8% completeness. We expect the depth and completeness relation in DES Y1 data to be similar to those in Zhang et al. (2017) because of consistent DES imaging/coadding strategies over the years.

For each of the objects, we mask the region enclosed within its 2.5 Kron radius (Kron 1980) ellipse. The major and minor axes of the ellipse as well as its inclination are computed by SExtractor in the DES coadd catalog. Analytical studies indicate that more than 90% of the galaxy luminosity is enclosed within this elliptical aperture (Graham & Driver 2005). To ensure that bright stars are properly masked, we set a minimum masking radius along the major axis of the ellipse to be $48-2\times mag\_auto\_i$ pixels (1 pixel = 0.263″) according to a visual check of the relation between Kron magnitudes and apertures.

The CG identified by the redMaPPer algorithm is excluded in the masking procedure.

3.2. Light profile estimation

We measure the light profiles of each of the clusters from the coadded images. The light flux profiles are computed as the averaged pixel values in radial annuli centered on the redMaPPer CGs, excluding the masked regions. The profiles are then averaged, i.e. ”stacked”, for an ensemble of clusters to detect ICL.

The cluster selection criteria eliminate undesirable features near clusters because cluster centers are selected to be 0°.15 away from bright stars or nearby galaxies, but the extended profiles of these undesirable features may still contribute to ICL detection, which is stronger at the cluster outskirt, causing a nonflat radial trend (see the top panels of Figures 3 and 4). To offset this trend, we subtract the light profiles of redMaPPer random points selected with the same criteria in terms of sky coverage (Section 2.1).

We assign the clusters to 40 sky regions111This number is high enough to estimate standard deviations, while keeping each region sufficiently large to include noise caused by large-scale structures, e.g., superclusters, filaments, etc. using a kmeans clustering algorithm (Steinhaus 1956; MacQueen 1967) implemented for celestial coordinates222https://github.com/esheldon/kmeans_radec. The 2000 redMaPPer random points are also assigned to the 40 regions. We measure the light profiles of each of the random points and then average them in each of the 40 kmeans sky regions. For each of the cluster light profiles, we subtract from it the averaged random light profile of its corresponding kmeans region. To further reduce cluster-to-cluster variations of the background level, for each cluster, we derive the average flux level value at the edge of the image of each cluster which is 0°.145 from the cluster centers, and subtract it from the overall cluster light radial profile. Thus by design, for each cluster, its radial flux profile reaches 0 at 0.145 deg, which is 1.73 to 2.34 Mpc from the cluster centers depending on the cluster redshift, which is outside the $R_{200c}$ of a typical redMaPPer cluster (McClintock et al. 2019).

The above region associations are also used to estimate light profile uncertainties. After the subtraction procedures, for each of the 40 aforementioned kmeans regions, we average the cluster light profiles in the region to acquire a ”stacked” measurement. The profiles of the 40 kmeans regions are then combined using the jackknife technique (Lovász et al. 1986) to estimate the mean and uncertainties of the final stacked profile.

In Figure  3, we show the final stacked profile of the clusters before and after the subtraction procedures, estimated with the jackknife technique. The subtraction steps eliminate sky residual and reduce correlated noise in the final stacked profile, which is more evident in the random test shown in the next section. If not explicitly stated, uncertainties presented in this paper are all estimated using the jackknife approach. Note that, by construction, they thus do not represent the cluster-to-cluster variation of the profiles.

3.3. Test with Random Points

We first apply the methods to a set of redMaPPer random points of the same size as the cluster sample. These are selected from the overall set of random points in the same way as our cluster sample from the redMaPPer catalog. The set of random points goes through the image- and profile-processing steps described in earlier sections. Figure 4 shows the stacked light profiles of this random point sample, before and after the profile-subtraction procedures (Section 3.2). The stacked light profiles before the subtraction procedures show a nonflat trend from the center of the image to the outskirt. After image inspections and noticing that random points/clusters are at a distance away from bright stars or nearby galaxies by design, we interpret the trend as the sky footprint selection effect of the redMaPPer algorithm. To our relief, the stacked light profile of random points is consistent with 0 after the subtraction procedures. Using the jackknife method, we estimate the $1\sigma$ surface brightness uncertainties to be $\sim$31.5 mag arcsec*-2* at $r=1~{}\mathrm{arcmin}$ and $\sim$33 mag arcsec*-2* at $r=7~{}\mathrm{arcmin}$, which means that our methods are sensitive enough to detect light above these surface brightness levels.

Notably, the noise fluctuations of the random points appear to be highly correlated at different radii, as shown in Figure 5. We notice that any image imperfections, such as scattered light, uneven sky-subtraction residual in different exposures, bright stars, etc., are likely to affect a significant portion of a single-exposure image (see the Appendix for some examples), causing correlated noise in the stacked profiles. The correlation is especially prominent at small radii, due to the smaller area. However, because the correlated fluctuations appear to be within the 3$\sigma$ uncertainty range in the random test (Figure 4 bottom panel), we do not consider it further in this paper.

To test the accuracy of the ICL measurement using our method, we inject a simulated ICL+CG distribution into the images of the random points. The ICL+CG radial distribution is simulated as a combination of three Sersic profiles with parameters described in Section 5.1, and injected into the processed and masked images of the random points, which then go through the light profile estimation procedures described in Section 3.2. Poisson photon noise from the ICL+CG model is not included in the injection to create an ideal situation of low-noise ICL+CG observation.

The result of this test is presented in Figure 6, showing the recovered radial profile of the simulated ICL+CG in comparison and the injected model. Because we assume that the ICL flux reduces to 0 at $\sim$ 0°.145, the ICL+CG model shown in the figure has also been subtracted by its model values at $\sim$ 0°.145. Our method has recovered the injected ICL+CG model with excellent accuracy out to $\sim$ 7 arcmin from the injection center, down to a surface brightness level of 30 mag/arcsec2.

4. DECam Point-Spread Function (PSF)

Literature studies show that the PSF of imaging instruments mounted on ground-based telescopes have low surface brightness wings that extend tens of arcminutes (Moffat 1969; King 1971; Racine 1996; Bernstein 2007). Detection of ICL is not to be confused with this feature. In this section, we measure and model the DECam (Flaugher et al. 2015) PSF in the $g-$ and $r-$ bands as a tension test of our methods, and also to determine whether or not scattered light from an extended PSF could have any significant effect on our detection of ICL.

To measure the central component of the PSF, we select a set of stars with magnitudes $19<r<19.3$, as measured by the DES MAG_PSF quantity. These point-source objects are separated from galaxies using “the modest classifier” from the DES Y1 coadd catalog (Drlica-Wagner et al. 2018). To further improve the purity of the point-source selection, we apply a color cut $\texttt{MAG\_PSF\_R}-\texttt{MAG\_PSF\_I}>1.5$ as recommended in Bechtol et al. (2015). This set of stellar objects is below the saturation limit of the DES single exposures in $g$ and $r$, and is thus useful for quantifying the central component of the PSF inside 10″. We derive the stacked light profiles of these stellar objects as well as their uncertainties using methods outlined in Section 3.

To study the farthest reaches of the DECam PSF, we use a set of very bright stellar objects selected from the Gaia 333https://www.gaia-eso.eu database (Gaia Collaboration et al. 2016). The first set contains 57 bright stellar objects of $\sim 8$ magn in the Gaia $G-$band. Although these objects are saturated in DES images in the central 3″, they give high S/N when studying the outskirts of the PSF beyond 20″. Finally, to derive the PSF radial profile between 10″and 20″, we select a set of moderately bright stellar objects which contains 100 stars of $\sim$14 mag in the Gaia $G-$band. These objects are also saturated in their cores in DES images, but still have sufficient S/N at the radii where the unsaturated stars have no signal left above the noise level. To compute the light profiles of these bright objects, we adopt the same image-processing procedure as described in Section 3, and then a similar profile estimation procedure excluding the subtraction of random point profiles and the requirement of no bright stars within 526$\arcsec$. The light profile uncertainties of the bright star are estimated with a bootstrap approach by randomly sampling the stars.

Figure 7 shows the radial profiles of the PSF in the $g-$ and $r-$ bands for a composite of faint, intermediate, and bright stars. The PSF is well characterized as a Gaussian of FWHM = 1.15″(FWHM = 1.043 ″for $r$, and 1.18 ″for $g$) at $r<1.1$″, and by a three-component broken power law, $f\approx r^{-4.13},~{}\>1.1\arcsec<r<7.5\arcsec,~{}f\approx r^{-2.46},~{}7.5\arcsec<r<60\arcsec$, and $f\approx r^{-1.16},~{}60\arcsec<r<450\arcsec$. The effective seeing was slightly different for the $g-$ and $r-$ bands, yielding a slight difference in profiles in the inner arcsecond. A Moffat profile fits the inner profile to somewhat larger radii (1.5″) versus a Gaussian (1.1″). The red shaded band gives a typical profile range for a 14th magnitude $r$-band star at intermediate radii, and the blue shaded band gives the typical spread in $g$-band scattered light at larger radii. The slopes are identical in the $g$ and $r$-bands within the uncertainties.

With these derived PSF profiles, for a point source, we estimate that about 80% of the light is contained within a radius of 1″, 98.3% within 5″, 99.6% within 50″, and 99.9% within 400″. Note that for the typical clusters measured here, 1 arcsec = 3.943 kpc (at redshift 0.25). In section 7.2, we further discuss the PSF effect on the interpretation of ICL results.

5. ICL+CG Profile

In figure 8, we show the derived light profiles of the individual clusters as well as their stacked light profile. These are the raw ICL+CG profile measurements in this analysis. Although on small scales within 0.5$\arcmin$ of the cluster centers the individual clusters show significant signals of ICL+CG, stacking enables the light profile measurement of ICL at distances beyond $\sim 0.5\arcmin$ out to $\sim 6\arcmin$. In the rest of this section, we consider the residual light contribution from cluster galaxies and derive results such as ICL+CG color and ICL+CG integrated luminosity. Although this paper is aimed at analyzing the properties of ICL, we do not attempt to quantitatively dissect ICL and CGs in the analyses. We qualitatively refer to the diffuse light outside a few kpcs of the CG center as the ICL, and apply a ICL/CG radius cut when making quantitative evaluations.

5.1. “Pure” ICL+CG profile

With the procedures described in Section 3, we derive the ICL+CG light profiles of each cluster. Unfortunately, not all of the light in the raw profile is associated with the ICL or CG. The light from the cluster foreground and background structures has been accounted for with the random point subtraction process in Section 3. Faint, unmasked cluster galaxies and the cluster galaxy light leaking outside the masking apertures constitute a cluster light residual in the raw ICL+CG measurements. We estimate the residual from the light profile of cluster galaxies.

We first derive the stacked light profile of non-CG cluster galaxies (cluster satellite galaxies). The methods described in Section 3, including the background residual and random light profile subtractions, are applied to the same set of clusters, but masking only objects brighter than 18 mag. Subtracting the previously derived light profiles with masking to 22.4 mag gives the radial light profile of non-CG cluster galaxies between 18 and 22.4 mag.

For each of the clusters, we assume the luminosity of non-CG cluster galaxies (cluster satellite galaxies) to follow a single Schechter function (Schechter 1976), and derive (1) the luminosity fraction between the 18 and 22.4 mag range and (2) the luminosity fraction beyond 22.4 mag. For the Schechter function parameters, we assume a faint-end slope of -1.0 and the characteristic magnitude of 1.29 $L_{\odot}/h^{2}$ at redshift 0.25 (Hansen et al. 2009). We further assume that the light profile of cluster satellite galaxies in different luminosity ranges have the same radial distributions (same shape of the radial light profiles). With the above luminosity fraction estimates and the measurement of non-CG cluster galaxy light profile between 18 and 22.4 mag, we derive the light profiles of cluster satellite galaxies and the light profile of satellite galaxies fainter than 22.4 mag. The cluster satellite galaxies below 22.4 magnitude makes up $\sim$ 4% of the total cluster satellite galaxy light in the redshift range of 0.2-0.3.

Our masking apertures (2.5 Kron radii) may not fully enclose the light of cluster galaxies due to their extended radial profiles. Graham & Driver (2005) estimated that the unmasked component makes up $0.1-9.6\%$ of the Sérsic (Sérsic 1963) galaxy profiles. We aggressively assume the masking residual to be 9.6% of the cluster galaxy light. From the raw ICL profile of each cluster, we subtract the light of faint, unmasked cluster satellite galaxies below 22.4 magnitude and the residual light of the masked galaxies assuming that they both have the same radial distribution as the general cluster satellite galaxy population.

We correct the derived “pure” ICL+CG light profiles for wavelength shifting (K-correction) and passive redshift evolution (e-correction) to an observer redshift of 0.25 using a stellar population template444http://www.baryons.org/ezgal/ (Bruzual & Charlot 2003; Mancone & Gonzalez 2012) with a single starburst of metallicity $Z=0.008$ at $z=3.0$. These corrected profiles are interpolated onto a physical radius grid and then averaged to compute a stacked “pure” ICL+CG profile, using the jackknife method.

The radial profiles of ICL+CG are shown in Figure 9. The raw measurement is dominated by pure ICL+CG light to $\sim 100$ kpc, beyond which the cluster galaxy residual becomes relevant. The subtraction of the cluster galaxy residual introduces significant noise to the pure ICL+CG detection. Nevertheless, the stacked ICL+CG profile is measured with high S/N to 1 Mpc, reaching 30 mag arcsec*-2*.

Similar to the finding in Kravtsov et al. (2018), the “pure” ICL+CG radial profile is best described with a combination of three Sersic models (Figure 10): a dominant core of Sersic index $1.34\pm 0.023$ and radius $9.13\pm 0.24$ kpc within the inner 10 kpc, a bulge between 30 and 100 kpc with Sersic index $3.1\pm 0.08$ and radius $52.1\pm 2.2$ kpc, and a diffuse component dominant outside 200 kpc with index $2.1\pm 0.4$ and radius $2.6\pm 0.7$ Mpc. The fitted Sersic profile parameters are listed in Table 1.

5.2. ICL+CG Integrated Luminosity

A significant fraction of the cluster stellar luminosity may be contained in ICL+CG, which is a debated topic in the literature. With the ICL+CG surface brightness profile derived in the previous sections, we examine the integrated luminosity profile of ICL+CG at different radii as

[TABLE]

$l(r^{\prime})$ is the ICL+CG surface luminosity profile from the previous subsection, corrected to an observer frame of redshift 0.25. We compute the integrated ICL+CG luminosity for both the “raw” and “pure” ICL+CG radial profiles, the results of which are shown in Figure 11.

These integrated ICL+CG luminosity profiles are measured with high S/N out to 1 Mpc. Within $\sim$ 30 kpc, the integrated luminosity increases steeply with radius because of the contribution from the cores of the CGs. Outside $\sim 30$ kpc, the increase slows down, indicating that the light profile has transitioned into a diffuse component.

The integrated luminosity within CGs and ICL is an important quantity in understanding CG and ICL formation. For example, Zhang et al. (2016) speculated that 30 to 60% of the stars that merged into the CGs need to be deposited outside 32 kpc at redshift 0 in order to explain the redshift evolution of the CGs. We find that the ICL integrated luminosity between 32 and 200 kpc makes up $\sim$ 42% of the total ICL+CG luminosity. The ICL integrated luminosity between 32 kpc and 1 Mpc is even more significant, $\sim$ 3 times as luminous as the luminosity enclosed within 32 kpc. Interestingly, the integrated luminosity of the ICL does not seem to converge within the radius range of 1 Mpc. For a more thorough consensus of ICL total luminosity, future studies may have to investigate an even larger radius range.

We also compute the integrated luminosity of the non-CG cluster galaxies (cluster satellite galaxies), using the light profile measurements of these galaxies in section 5.1. The derivation of this satellite galaxy luminosity is susceptible to modeling uncertainties of the galaxy luminosity function and the galaxy-masking apertures. Nevertheless, the results are shown in Figure 12, in comparison with the integrated luminosities of ICL+CG.

Within 200 kpc, the integrated luminosity of ICL+CG surpasses the total luminosity of cluster satellite galaxies (non-CG cluster galaxies). Outside 200 kpc, the satellite galaxy luminosity becomes paramount to that of ICL+CG. We further estimate the ratio of ICL+CG in the total cluster stellar luminosity, which is computed as the sum of ICL+CG and non-CG cluster galaxies. We refer to this quantity as the ICL+CG fraction, as in Burke et al. (2015). The ICL+CG fraction gradually drops when we enlarge the radius of the luminosity integration, reaching $44\pm 17$% at 1 Mpc.

5.3. ICL+CG color

Using our measurements of the ICL+CG profiles in both the $g-$ and $r-$bands, we further derive the ICL+CG color profile (Figure 13). The $g-$band ICL+CG measurement has a greater level of noise compared to the $r-$band measurement (Figure 13 upper panel), because of less accurate estimation of cluster galaxy residual. Because the color measurements require high signal to noise flux measurements in both the $g$ and $r$ bands, our color derivation is only robust to 90 kpc from the center. Outside 90 kpc, the color measurements are dominated by noise.

We compare the ICL+CG colors to those of the cluster CGs and the satellite red sequence galaxies selected by redMaPPer. These galaxy colors are derived with the MODEL_MAG magnitudes in the DES database (Drlica-Wagner et al. 2018), K-corrected to $z=0.25$ as is done for the ICL+CG profiles. Overall, the core component of the ICL+CG profile appears to be consistent with the CG measurement in the DES database, but the ICL color beyond $\sim 20$ kpc is more consistent with the average color of the satellite red sequence galaxies. The ICL+CG color also displays a radial trend of becoming bluer at a larger radius. We find that the color trend in the central 10 kpc of the ICL+CG combination is likely caused by the PSF differences in the $g-$ and $r-$bands (see section 7.2 for detailed discussions). Between 10 and 90 kpc, a $\chi^{2}$ minimization gives a significant radial gradient,

[TABLE]

indicating a robust color radial trend.

6. Cluster Mass Dependence

6.1. More ICL in Richer Clusters

The redMaPPer algorithm adopts richness, denoted by $\lambda$, as a mass proxy, which is defined as the number of cluster red sequence galaxies above a luminosity threshold. Using this richness quantity as a cluster mass proxy, we examine the variation of ICL surface brightness with mass.

The $\sim$ 300 clusters studied in the paper are further divided into three subsamples, $20\leq\lambda<30$, $30\leq\lambda<50$, and $\lambda\geq 50$, representing three cluster mass ranges with mean masses of $10^{14.09}$, $10^{14.30}$, and $10^{14.61}$ $\mathrm{M_{\odot}}/h$ (Melchior et al. 2017). We derive the stacked ICL+CG profiles of these three subsamples as shown in Figure 14. The ICL is more luminous in richer and hence more massive clusters. The core components within $\sim$ 10 kpc of the ICL+CG profile appear similar among the three cluster richness subsamples, in agreement with the inside-out growth scenario (e.g., van Dokkum et al. 2010; van der Burg et al. 2015) wherein the cores of the CGs form early and the accreted CG stellar content is deposited onto the outskirt at a later time. The ”inside-out” trend is further demonstrated in the integrated luminosity profiles (Figure 15). Richer clusters are more abundant in the total ICL+CG luminosity, but the differences start outside 20 kpc and rapidly enlarge at larger radii.

6.2. ICL+CG and Cluster Mass Distribution

Interestingly, after scaling the ICL+CG profiles by the radii of clusters, the radial profiles of ICL+CG in different richness ranges appear similar (Figure 14). In this exercise, instead of stacking together ICL+CG radial profiles in terms of physical distances, we stack them in radial bins of $r/R_{200m}$, with $R_{200m}$ derived from the mass-richness relation in Melchior et al. (2017). This is the first evidence of ICL distribution being “self-similar”, as in the self-similarity of cluster mass or gas distributions (e.g., see a review in Kravtsov & Borgani 2012) wherein their radial profiles are indistinguishable after scaling with a characteristic radius.

We further examine the resemblance between the ICL and cluster mass distributions, by comparing their radial profiles. The state-of-the-art constraint of cluster mass distribution (including dark and baryonic matter) comes from weak lensing studies, yet cluster surface mass profile is not a direct observable from weak lensing555Weak lensing directly measures excess matter density enclosed within a radius, rather than the surface mass density at that radius., and its reconstruction can be noisy. In this comparison, we use a cluster surface mass density model from McClintock et al. (2019) 666https://tmcclintock.github.io/code/cluster_toolkit, which is a combination of a Navarro-Frenk-White (NFW) model (Navarro et al. 1997) and a two-halo matter correlation model, projected onto the plane of the sky:

[TABLE]

In the above equations, $\xi_{NFW}(r|M,c)$ and $\xi_{2h}(r|M)$ represent the cluster matter correlation functions in the one-halo and two-halo regimes of a cluster with mass $M$ and concentration $c$. $\rho_{m}(z)$ represents the universe mean matter density at redshift $z$. $\Sigma(r|M,c,z)$ is the surface mass density at distance $r$ on the plane of the sky, integrated along the line-of-sight distance $\chi$. We compute $\Sigma(r|M,c,z)$ at the mean masses of the three cluster richness subsamples, $10^{14.09}$, $10^{14.30}$ and $10^{14.61}$ $\mathrm{M_{\odot}}/h$, assuming the concentration to be 5 and the redshift to be 0.25.

The ratio between the ICL+CG luminosity and the cluster surface mass density is shown in Figure 16. This luminosity-to-mass ratio is at its highest within 100 kpc. As found of the cluster matter distribution in multiwavelength observations (Newman et al. 2013b, a), the core of galaxy clusters is dominated by the stellar mass of the CGs, which would explain the high ICL to cluster mass ratio near the center.

The ratio drops with enlarging radius. Between 100 kpc and 1 Mpc, the ICL luminosity-to-cluster-mass ratio no longer displays a noticeable radial trend, which is $10^{-2.85\pm 0.09}$ at $r=215$ kpc and $10^{-2.84\pm 0.4}$ at $r=985$ kpc, in the medium richness bin. Because the cluster surface mass density (McClintock et al. 2019) and the ICL surface brightness (Section 5.1) both drop by a factor of $\sim 10$ in this radial range, the flatness of the ICL+CG luminosity-to-mass ratio indicates that the ICL luminosity distribution closely follows the cluster mass distribution. Moreover, throughout the 100 kpc to 1 Mpc radial range, the ICL+CG luminosity to cluster mass ratios are similar among the three cluster richness subsamples, a further evidence that the ICL+CG luminosity profile is self-similar and traces the overall cluster mass distribution.

The 100 kpc to 1 Mpc radial range is also coincident with where the ICL+CG profile transitions into an extended diffuse component (Figure 10). Thus, this transition potentially marks a physical separation of ICL and CG.

7. Systematics

7.1. Sky Subtraction

Sky subtraction near bright galaxies or stars is a known difficulty in processing imaging data. It is common for the processing software or pipelines to overestimate the sky level around bright objects because of light contamination from the objects, and hence oversubtract the sky background. This kind of effect is particularly detrimental to the detection of low surface brightness light near bright galaxies, and in our case, the detection of faint diffuse light around the CGs.

In this paper, the sky levels are estimated and subtracted during three processes:

Sky level is estimated for each exposure image and subtracted before the single-exposure images are coadded. 2. 2.

In section 3.2, residual sky background is estimated and subtracted again by sampling the light profiles of random points that cover the footprint of the cluster sample to account for cluster selection effects (i.e. we only look for clusters in sky regions that are not in immediate proximity to bright stars or nearby galaxies). 3. 3.

In section 3.2, residual background light is again estimated as the flux level at a distance of $\sim$ 7 arcmin (about $\sim$1.6 Mpc at redshift 0.25) from the cluster center and then subtracted to reduce the cluster-to-cluster variations of the ICL detection.

In the above steps, (2) and (3) are by definition either far away or not correlated with cluster locations, hence we do not consider them to be contaminated by the light of the CGs. As previously mentioned in Section 2.2, the sky level estimation in step (1) is derived through modeling the variation of light over the 3 deg2 DECam focal plane and thus should be relatively insensitive to the presence of individual bright galaxies or stars.

Nevertheless, we investigate whether or not the sky level in step (1) may have been overestimated near the cluster center. In this exercise, we look into the sky background maps estimated for the single-exposure images covering the redMaPPer CGs, and directly compare the sky values estimated at the location of the cluster center and at locations away from the center.

For each of the $\sim$ 300 clusters analyzed in the paper, we search for single-exposure images that cover at least the central $2.06\times 2.06$ arcmin2 region of the cluster777The size of each DES CCD image is approximately $8.77\times 17.5$ arcmin2, and varies slightly depending on the CCD location on the focal plane.. Because of the survey nature of DES, the central region of each cluster may not fully appear within one exposure image, and we would like to avoid comparing sky levels in different exposures as they may vary significantly. For those single-exposure images that fully overlap with the cluster’s central $2.06\times 2.06$ arcmin2 regions, we compute the sky flux difference between the cluster center location and locations at distances away from the centers, with the differences rescaled to a magnitude zero point of 30 for each exposure.

Figure 17 shows the sky value differences between the cluster center and locations that are at distances of 0.14, 0.45, and 1.36 arcmin, corresponding to physical distance separations of 33, 106, and 320 kpc at $z=0.25$. If the sky level is indeed overestimated at the cluster center, we would expect the differences to be shifted to the positive side, and the shifts may become larger at larger distances. However, neither of these trends appear to be present in this test. The mean of the differences is at a flux level of $\sim 0.1$ arcsec*-2* between the cluster centers and locations 1.36 arcmin away from the centers, significantly under the ICL flux limit in this work (30 mag $\mathrm{arcsec}^{-2}$ at $\sim$ 1 Mpc corresponds to a flux of $\sim 1$ arcsec*-2* with a magnitude zero point of 30).

7.2. PSF Effect

A frequently debated topic in the study of galaxy low surface brightness light envelope is the effect of PSFs. As noted in Moffat (1969); King (1971); Racine (1996); Bernstein (2007), the telescope and instrument point spread functions almost always have an extended component that radially extends several arcseconds or arcminutes. This is also shown to be true for the DECam images on the Blanco telescope in Section 4. In previous studies about whether or not faint light envelope light exists around galaxies (D’Souza et al. 2014; Duc et al. 2015) including the Milky Way and other spiral (Sackett et al. 1994; Zibetti et al. 2004; Zibetti & Ferguson 2004; van Dokkum et al. 2014; Zheng et al. 2015) or elliptical galaxies (Tal & van Dokkum 2011), de Jong (2008); Sandin (2014, 2015) showed that the extended component of the telescope and camera PSF may partially or largely account for the observed galaxy light halo and its color gradient.

In this section, we consider whether or not the extended component of DECam/Blanco PSFs can explain the ICL measurements in this paper. The measurement of the extended wings of the PSF and ICL profiles do not have sufficient S/N to allow the smooth deconvolution of PSF from ICL. To evaluate the effect of PSFs, we creates a simulated ICL profile, convolve it with a PSF, and examine the difference in the ICL profiles before and after the convolution (Figure 18). The unconvolved ICL profile is modeled with the three Sersic components as described in Section  5.1, and we convolve it to the extended PSF model derived in Section 4. We find that the PSF convolution only significantly alters the shape of the ICL within the inner 10 kpc, and is negligible outside 100 kpc at redshift 0.25.

To estimate the effect of PSFs on the ICL color measurements, we compute the color shifts when convolving the ICL profile model to the $g-$ and $r-$band PSF models from Section  4. Driven by the wider PSF FWHM in the $g-$ band compared to that of the $r-$band, PSF convolution produces an artificial trend of bluer color at a large radius within the inner 10 kpc, which informs us not to analyze the ICL color gradient in this radial region. The color shift caused by PSF convolution is negligible outside 10 kpc.

7.3. Masking of Cluster Galaxies

As discussed in Section 5.1, the raw ICL measurements contain light from faint unmasked cluster galaxies, as well as residual light outside the masks of the already masked galaxies. These contributions are subtracted from the raw ICL measurements to derive pure ICL profiles, by assuming (1) a luminosity function of cluster galaxies and that (2) 9.6% of the galaxies’ light leaks to the outside of the masks.

As discussed in Section 5.1, we assume a single Schechter luminosity function to estimate the light residual from faint unmasked cluster galaxies. As with the parameters of the Schechter function, given recent measurement uncertainties of the cluster galaxy luminosity function (e.g., Lan et al. 2016; Ricci et al. 2018), we find that varying the faint-end slope parameter by $\pm 0.2$ or the characteristic magnitude parameter by $\pm 0.5$ mag changes the estimated fraction of residual light by $\sim 1\%$.

Notably, we have ignored possible deviations from a single Schechter function caused by faint dwarf cluster galaxies. Lan et al. (2016) estimated that the faint component of cluster galaxies has a normalization factor about 8% of the brighter ones, with a characteristic magnitude $\sim 3.5$ mag fainter and a faint-end slope of $\sim-2.0$. Such a faint dwarf galaxy component makes up less than 1% of the total light of the bright component, and hence ignored in our analysis.

The second contamination to the ICL measurement, the residual light outside the masks of the already masked galaxies, is assumed to be 9.6% of the total light of those galaxies, which is likely an overestimation according to Graham & Driver (2005). We do not see a significant reduction in the raw ICL+CG brightness when adjusting the masking apertures from 2.5 Kron radii to 3.5 Kron radii. Thus, we expect the ”pure ICL”+CG profile computed in this paper to be an underestimation, and the truth ”pure” ICL+CG profile to be between the ”pure” and ”raw” measurements.

Finally, our cluster galaxy residual subtraction process does not consider the existence of the recently discovered ultra-diffuse galaxies (see van Dokkum et al. 2015, and the referencing papers). Neither have these galaxies been given consideration in the cluster galaxy luminosity function models used in the paper. Therefore, our ”pure” ICL+CG profile contains light from these ultradiffuse galaxies. Distinguishing them will require a better understanding of these lesser-known objects in the future.

8. Summary and Discussion

This paper develops methods to measure low surface brightness ICL centered on the cluster CG, which are tested with random points and simulated ICL+CG profiles, and also applied to the study of the DECam PSFs in the $g-$ and $r-$bands for averaged DES Y3 observations. We estimate the light residual in raw measurements of ICL+CG profiles from cluster satellite galaxies and consider systematic effects from the sky background subtraction procedure, from the extended wings of the DECam PSFs and from assumptions about the luminosity function of cluster satellite galaxies.

Following the high S/N detection of ICL+CG radial profiles, we study a variety of ICL properties – radial distribution, integrated luminosity, color profile, and the connection between ICL distribution and the cluster mass distribution. The ICL extends out to 1 Mpc from the cluster centers at a surface brightness level of 30 mag arcsec*-2*. The ICL outside 30 kpc of the CGs, is $\sim 3$ times as luminous as the ICL+CG central component within 30 kpc. The ICL also displays an interesting self-similarity feature wherein their radial profiles appear similar after scaling by the cluster $R_{200m}$, and the ICL radial distribution closely follows the overall cluster mass distribution.

8.1. Comparison to Previous Observational Studies

A very similar work to this one is described in Zibetti et al. (2005) (hereafter Z05) wherein $g$, $r$, and $i$ images from the Sloan Digital Sky Survey (SDSS) of $\sim 600$ galaxy clusters at redshift 0.2-0.3 were stacked, and ICL+CG profiles out to 700 kpc from the CG centers were derived. Our analyses shared some similar ICL findings to Z05 –both works have found the ICL to have a shallower profile than the CGs in the core, bluer ICL $g-r$ color at a greater distance, and more luminous ICL in richer clusters, yet the core component of the CGs have a similar luminosity regardless of the cluster richness.

Not all results in the two works agree, however. Notably, the ICL profile derived with our method is much brighter than that in Z05, as shown in Figure 19. While Z05 estimated the ICL surface brightness to be 32 mag arcsec*-2* in the $r-$band at around 700 kpc, our results show that the ICL surface brightness is above 32 mag/arcsec*-2* even at 1 Mpc. Several methodological differences may have contributed to this discrepancy. Most importantly, Z05 chose to evaluate the sky background level as the average flux in a 100 kpc annulus centered at 1 Mpc from the cluster center, which offsets the ICL flux level to 0 in that radial range. Our nominal method offsets the ICL flux level to 0 at $\sim$ 1.8 Mpc, and if we offset our measured ICL flux to 0 in the same radial range as Z05, our ICL surface brightness drops to $\sim$ 31.5 mag arcsec*-2* at 700 kpc, close to the estimation in Z05. The rest of the differences may be due to the applications of K- and e-corrections (wavelength shifting with redshift and galaxy passive evolution), as well as cluster selection methods. Overall, our analysis is based on a much improved cluster sample with better techniques to quantify sky background and PSFs and account for sky selection effects associated with the cluster sample.

Many other studies have investigated ICL properties to a great distance, albeit with fewer clusters (e.g., Krick et al. 2006; Gonzalez et al. 2007; Krick & Bernstein 2007; Rudick et al. 2010; Guennou et al. 2012; Giallongo et al. 2014; Montes & Trujillo 2014; Burke et al. 2015). With deep images from Hubble or ground-based telescopes, these works have measured the ICL to surface brightness levels beyond 25 mag arcsec*-2* and out to at least 100 kpc from the CG centers, sometimes close to a surface brightness limit of 30 mag arcsec*-2*.

A frequently investigated ICL property in these works is the fractional contribution of the ICL to the total cluster stellar content. Around redshift 0.2-0.3, there is a great dispersion in the reported values. Burke et al. (2015) found the ICL fraction to be around 20% at redshift 0.25 in the CLASH survey data, while Krick & Bernstein (2007) estimated the fraction to be 6%-20% in 10 nearby clusters, or around 10% at redshift 0.25. Going to a slightly lower redshift range of 0.16-0.20, Feldmeier et al. (2004) estimated the ICL fraction to be around 10% with four Abell clusters. Note that these three works have defined the ICL as the diffuse light fainter than a surface brightness limit, usually around 25-26 mag arcsec*-2* in $B-$ or $V-$bands, which roughly correspond to $>$ 20 kpc in our analysis. On the higher side, Gonzalez et al. (2007, 2013) estimated the ICL+CG to make up 33% of the total cluster stellar mass within $r_{200}$ in the redshift range of 0.03-0.13, or 20%-50% in the redshift range of 0.05-0.24 within $r_{500}$. Toledo et al. (2011) found the ICL+CG fraction to be $>$40% within 500 kpc for one cluster at redshift 0.29, and Morishita et al. (2017) found the fraction within 300 kpc to be between 15% and 60%. The reported ICL fraction varies with the radial range of the comparison. For one cluster of redshift 0.44, Presotto et al. (2014) estimated the ICL+CG fraction to be $>$ 50% within 100 $\mathrm{kpc/h}$, or 20 % within 350 $\mathrm{kpc/h}$, or 8.2 %within $R_{500}$. Our measurement of the ICL+CG luminosity fraction to be $44\pm 17$% at radius 1 Mpc is on the high end compared to the previous observations, although still within the range of those reports.

Another frequently visited ICL property is the color, which informs us about its stellar population. Most imaging studies have found ICL+CG colors to be comparable to an old stellar population with a radial trend of becoming bluer at a greater distance (e.g., Montes & Trujillo 2014; Presotto et al. 2014; Mihos et al. 2017; Huang et al. 2018a; DeMaio et al. 2018), which is further supported by spectroscopic studies of the ICL+CG stellar composition (Adami et al. 2016; Edwards et al. 2016; Johnston et al. 2018). Our results of the ICL color being consistent with the cluster red sequence galaxies and becoming bluer with distance agree with those previous conclusions (e.g., Z05, as shown in Figure 13). The color gradient of ICL in particular indicates that the ICL is likely produced from disruptions of dwarf galaxies or tidal stripping of cluster member galaxies. The former mechanism is expected to create an ICL color gradient as the disrupted dwarfs at different radii have different masses and colors. The later mechanism, the tidal stripping of galaxies, becomes stronger closer to the cluster center and will be able to strip the redder, inner parts of the cluster galaxies. Major mergers are not likely to be the main ICL formation mechanism, as it would give rise to a more uniform ICL color profile (e.g. Eigenthaler & Zeilinger 2013; Contini et al. 2018b).

8.2. Comparison to Simulations

A few simulation studies have attempted to include the ICL (Rudick et al. 2006; Murante et al. 2007; Purcell et al. 2007; Guo et al. 2011; Rudick et al. 2011) but balancing between the properties of the general cluster galaxy population and the properties of CG and ICL is a challenge (Puchwein et al. 2010).

Using semianalytical methods, Contini et al. (2014, 2018a) investigated different pathways to ICL formation, such as stripping and merging relaxation, and found that regardless of these mechanisms, the ICL has experienced rapid redshift evolution in the late times. At redshift 0, the ICL+CG may make up 30%-50% of the total cluster stellar content, or close to such a quantity at redshift 0.25.

Recent developments in hydrodynamics simulations are also producing encouraging ICL results. In high-resolution dark matter and baryon resimulations of massive cluster-sized halos incorporating AGN feedback, Martizzi et al. (2014) were able to identify diffuse ICL envelopes, extending to a few hundreds of kpcs around the CGs. Between surface brightness limits of 25 and 27 mag arcsec*-2* in $V-$band, ICL makes up 20%-60% of the total light of the CG and ICL combination. The surface brightness limits applied in Martizzi et al. (2014) roughly correspond to a 30-80 kpc radial range. Our estimation of the fraction of the ICL in the ICL+CG combination in such a radial range is $28\%$ (Figure 11), which agrees with the reports in Martizzi et al. (2014).

Pillepich et al. (2018) carried out a detailed study of the ICL with the IllustrisTNG hydrodynamics simulation suites. In particular, Pillepich et al. (2018) investigated the ICL radial distribution, and found the ICL stellar mass profile to be as shallow as that of the dark matter, with a similar power-law radial dependence. The ICL stellar mass outside 30 kpc, or outside 100 kpc of the cluster center, scales with cluster mass with a power-law index close to 1. Our observation that ICL traces the overall cluster mass distribution supports such a finding. As with the amount of the ICL, for a $3\times 10^{14}M_{\odot}$ cluster, Pillepich et al. (2018) found that the ICL stellar mass outside 30 kpc of the CG center is about three time as massive as the inner 30 kpc of the CG+ICL combination ($3\times 10^{12}M_{\odot}$ versus $1\times 10^{12}M_{\odot}$), and the ICL outside 100 kpc is approximately equal to the CG+ICL stellar mass within 100 kpc ($2\times 10^{12}M_{\odot}$). The ICL+CG fraction in the total cluster stellar content is 50 %, which is on the high end of the previous observational studies. However, both the Pillepich et al. (2018) and the Contini et al. (2014, 2018a) estimations are in excellent agreement with the results from our analysis.

Modeling-wise, Behroozi et al. (2013, 2018) also included an ICL component, called intra-halo light (IHL), in their efforts to provide a self-consistent halo model to quantify the average galaxy star formation history from various galaxy observables. For a halo of mass $1\times 10^{14}\mathrm{M_{\odot}}$, Behroozi et al. (2018) also found that ICL/IHL has a rather late formation time, increasing rapidly by $\sim 3$ times between redshift 0 and 1. Near redshift 0, the ICL/IHL contributes $\sim 30$ % of the total halo stellar mass.

8.3. ICL Self-similarity

A rather interesting result from this analysis is that the ICL radial profile appears self-similar – it scales with cluster $R_{200m}$ inferred from cluster richness. In addition, starting from 100 kpc and out to 1 Mpc, the ICL profile appears to trace a theoretical cluster mass profile model. This is one of the first direct evidence establishing a connection between the ICL and the cluster mass radial distribution. Prior to this work, Montes & Trujillo (2019) examined the surface brightness contours of the ICL and the weak lensing mass maps of six galaxy clusters studied by the Hubble frontier program and found (1) visual similarity between the two and moreover, (2) compatible contour shapes with a quantitative shape estimator, in various radial ranges. Montes & Trujillo (2019) noted the potential of using the ICL observation to trace the cluster mass distribution. The results from our analysis provide additional direct evidence to such a conclusion.

Prior to the analyses in this paper and those presented in Montes & Trujillo (2019), hints of a connection between the ICL and cluster mass distributions can be found in a few observational and simulation studies that examine the radial profile of the ICL or the scaling between the ICL stellar mass and cluster mass. For example, Zibetti et al. (2005) found that the ICL profile is reasonably approximated by an NFW model (Navarro et al. 1997). A few other observational works have noted a stronger correlation between cluster mass and ICL stellar mass or luminosity outside the CG cores (e.g., DeMaio et al. 2018; Huang et al. 2018b). While the accuracy of the ICL tracing cluster mass distribution awaits further investigations, this similarity may provide another channel to explore the origin of the ICL, if different ICL formation mechanisms, e.g., tidal stripping or dwarf galaxy disruption, have different radial or cluster mass dependence. Furthermore, as noted in Montes & Trujillo (2019), this phenomenon, if confirmed, provides another extraordinary opportunity to observationally constrain the cluster mass distribution.

Another possible explanation of this phenomenon is that the ICL simply traces the luminosity distribution of cluster satellite galaxies, which then traces the cluster mass distribution. We have examined the radial light profile of the cluster satellite galaxies derived in Section  5.1 and compared it to the theoretical cluster mass profile model, but cannot determine if the ICL is a better or worse tracer of cluster mass than the satellite galaxies, given the large uncertainties of those measurements. It would be worthwhile to revisit this topic with a larger cluster sample in the future.

8.4. Outlook

With the successful detection ofthe ICL in this paper, an interesting topic to explore next is its redshift evolution, which will provide more clues about its origin. The depth and wavelength coverage of the DES produce volume-limited cluster samples up to redshift $\sim 0.65$ in Y1 data and up to redshift $\sim 0.8$ with the first three years of data (Y3). The increased volume of the redMaPPer cluster sample at a higher redshift and hence the increased sample size partially offset the dimming effect of the ICL with distance and may allow for ICL detection up to redshift $\sim 0.65$, which will provide preliminary answers to whether or not the ICL has experienced much redshift evolution.

The method developed in this paper can also be applied to studying low surface brightness light around other types of galaxies in DES data. We note that we have employed this method to analyze the faint light halo around luminous red galaxies (Rozo et al. 2016) in an upcoming study (Leung in preparation).

Finally, based on the ICL surface brightness results presented in this work, the ICL may become a significant systematic effect for cluster weak lensing studies with future cosmic surveys such as the LSST, because of ICL contamination in the flux measurements of weak lensing source galaxies. For this aspect of the ICL study, we refer the readers to Gruen et al. (2018).

Some of the ancillary data products shown in the Figures of this paper are available on the DES data release page https://des.ncsa.illinois.edu/releases/other/paper-data. Readers interested in comparing to these results are encouraged to check out the release page or contact the corresponding author for additional information. We thank the anonymous referee for the constructive comments and Stefano Zibetti for kindly sharing the data products in his paper for comparison with this work. We also thank Rebecca Bernstein, Gary Bernstein, Peter Behroozi, Risa Weschler, Tom McClintock and Eduardo Rozo for insightful discussions. We use Monte Carlo Markov Chain sampling in this analysis, which is performed with the PYMC (Patil et al. 2010) and emcee (Foreman-Mackey et al. 2013) Python packages. Funding for the DES Projects has been provided by the U.S. Department of Energy, the U.S. National Science Foundation, the Ministry of Science and Education of Spain, the Science and Technology Facilities Council of the United Kingdom, the Higher Education Funding Council for England, the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign, the Kavli Institute of Cosmological Physics at the University of Chicago, the Center for Cosmology and Astro-Particle Physics at the Ohio State University, the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University, Financiadora de Estudos e Projetos, Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento Científico e Tecnológico and the Ministério da Ciência, Tecnologia e Inovação, the Deutsche Forschungsgemeinschaft and the Collaborating Institutions in the Dark Energy Survey. The Collaborating Institutions are Argonne National Laboratory, the University of California at Santa Cruz, the University of Cambridge, Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas-Madrid, the University of Chicago, University College London, the DES-Brazil Consortium, the University of Edinburgh, the Eidgenössische Technische Hochschule (ETH) Zürich, Fermi National Accelerator Laboratory, the University of Illinois at Urbana-Champaign, the Institut de Ciències de l’Espai (IEEC/CSIC), the Institut de Física d’Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig-Maximilians Universität München and the associated Excellence Cluster Universe, the University of Michigan, the National Optical Astronomy Observatory, the University of Nottingham, The Ohio State University, the University of Pennsylvania, the University of Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the University of Sussex, Texas A&M University, and the OzDES Membership Consortium. Based in part on observations at Cerro Tololo Inter-American Observatory, National Optical Astronomy Observatory, which is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation. The DES data management system is supported by the National Science Foundation under Grant Numbers AST-1138766 and AST-1536171. The DES participants from Spanish institutions are partially supported by MINECO under grants AYA2015-71825, ESP2015-66861, FPA2015-68048, SEV-2016-0588, SEV-2016-0597, and MDM-2015-0509, some of which include ERDF funds from the European Union. IFAE is partially funded by the CERCA program of the Generalitat de Catalunya. Research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Program (FP7/2007-2013) including ERC grant agreements 240672, 291329, and 306478. We acknowledge support from the Australian Research Council Centre of Excellence for All-sky Astrophysics (CAASTRO), through project number CE110001020, and the Brazilian Instituto Nacional de Ciência e Tecnologia (INCT) e-Universe (CNPq grant 465376/2014-2). This manuscript has been authored by Fermi Research Alliance, LLC under Contract No. DE-AC02-07CH11359 with the U.S. Department of Energy, Office of Science, Office of High Energy Physics. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes.

Appendix A Examples of Noise Sources in ICL Detection

While examining the masked images of the clusters, we identified a few cases that increase the noise level of ICL detection. Three of them are noted below as examples. These noted contaminations are not structures associated with the redshift 0.2-0.3 clusters, therefore we expect their net effects to be canceled out in the process of ”stacking” and background subtraction. However, these cases can contribute correlated noise on a scale of a few pixels, or produce the noise spikes seen in the ICL figures.

• •

Bright stars or nearby galaxies (e.g., galaxies below redshift 0.05). The unmasked light residual in the outskirts of these objects can be bright enough to cause significant contamination to the faint ICL detection on a scale of tens of arcseconds. Examples of such cases can be seen in Figure 20. Note that we have implemented a random point subtraction process to account for the fact that the cluster sample selection method avoids looking for clusters in close proximity to these objects.

• •

Coadding discontinuity in the images. This can be caused by different background residual levels in the processed single-exposure images, or CCD readout differences in the same exposure. Because DES does not purposely place clusters at certain locations in the camera field of view, this effect should be canceled out in the ”stacking” process. An example of such a case can be seen in Figure 20.

• •

Transient objects (e.g., astroids). We noticed these objects as they may appear to be particularly bright in one exposure but disappear in another. These objects are often not cataloged in the DES coadd object database, as the object-detection process is only based on the coadded images of DES $r-$, $i-$ and $z-$band observations. For the same reason, these contamination cases are more common in the $g-$band rather than in the $r-$band. Again, their contamination can be canceled out by ”stacking”. An example of such a case is shown in Figure 21.