Probing Intra-Halo Light with Galaxy Stacking in CIBER Images

The raw imager data provides a time series for each pixel. We fit a slope to the time stream to obtain the photocurrent in each pixel and convert the values from the raw analog-to-digital units (ADU) to e⁻ s⁻¹ using known array gain factors.

The HAWAII-1 detector is linearly responsive to incoming flux over a certain dynamic range. For pixels pointing at bright sources, the detectors saturate and have a nonlinear flux dependence even for short integrations (Bock et al. 2013). In any pixel that collects more than 5000 ADU over the full integration only the first four frames are used in the photocurrent estimate. Hereafter, the term "raw image" refers to the photocurrent map after this linearity correction. Panel A of Figures 1 and 2 show the raw images of the SWIRE field in the CIBER 1.1 and 1.8 μm bands, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the absence of incoming photons, the detectors have a nonzero response, commonly referred to as "dark current," due to thermally produced charge carriers and multiplexer glow. The detector dark current is measured before each flight with the telescopes' cold shutters closed. We obtain a dark-current template for each detector by averaging 11 dark images and then subtracting each template from the corresponding raw images. The dark-current level in CIBER imagers is ∼0.1 e⁻ s⁻¹, less than 10% of the sky brightness. Panel B of Figures 1 and 2 show the dark-current maps of CIBER 1.1 and 1.8 μm bands, respectively.

We mask pixels that meet at least one of the following conditions: (1) a fabrication defect, (2) poor time-stream behavior, (3) abnormal photocurrents compared with other pixels, (4) a cosmic-ray strike, or (5) being on or close to bright point sources on the sky. The pixels satisfying criteria (1)–(4) comprise the "instrument mask" and a "source mask" is composed of pixels with condition (5).

3.3.1. Instrument Mask

3.3.2. Source Mask

CIBER images have a non-uniform response to a constant sky brightness across the detector array, known as the flat-field response. For each CIBER field, the flat-field is estimated by averaging the dark-current-subtracted flight images of the other four sky fields.

A laboratory flat-field measurement was also taken before the flight using a field-filling integrating sphere, a uniform radiance source with a solar spectrum (described in Bock et al. 2013). Ideally, this is a better approach to measure the flat-field since the one derived from stacking flight images contains fluctuations from the other fields that will not average down completely due to the small number of images. However, we found the flat-field from the integrating sphere is not consistent with the flight data on large spatial scales (see Zemcov et al. 2014), and therefore we do not use it in our analysis. The flat-field estimator for the SWIRE field in the CIBER 1.1 and 1.8 μm bands is shown in Panel E of Figures 1 and 2, respectively.

Throughout this work, we use nW m⁻² sr⁻¹ for the units of surface brightness (ν I_ν ). The calibration factor, C, which converts photocurrent (e⁻ s⁻¹) to intensity (nW m⁻² sr⁻¹) is derived in the following steps:

• 1.  
  Take the raw images, subtract the dark-current template, correct for the flat-field, and apply the instrument and source masks;
• 2.  
  Subtract the mean photocurrent in the unmasked region;
• 3.  
  For each star in the Pan-STARRS catalog, calculate the flux νF_ν in CIBER bands from m_1.1 and m_1.8;
• 4.  
  Sum the photocurrent in a 5 × 5 stamp centered on the source position; ¹²
• 5.  
  Repeat steps (3) and (4) for all the selected stars (see below) and take the average value of the flux ratio from (3) and (4) as the calibration factor C.

We select stars in the magnitude range 12.5 < m_1.1 < 16 for the 1.1 μm band, and 13.5 < m_1.1 < 17 for the 1.8 μm band. These magnitude ranges are chosen such that the brightest sources that saturate the detectors (even after nonlinear correction) are excluded. Faint sources are not used because of their low signal-to-noise ratio. We use a different magnitude range for each band as they have different point-source sensitivities. Panel F of Figures 1 and 2 shows the SWIRE field images masked by instrument masks at 1.1 and 1.8 μm, respectively, after flat-fielding and calibration.

The total sky emission is composed of the EBL and various foreground components, including zodiacal light (ZL), diffuse galactic light (DGL), and integrated starlight (ISL) from the Milky Way (Zemcov et al. 2014; Matsuura et al. 2017). ZL is the dominant foreground, approximately an order of magnitude brighter than the EBL (Matsuura et al. 2017). Nevertheless, with its smooth spatial distribution on degree scales, the ZL can be mostly removed by subtracting the mean sky brightness in each field. Panel G of Figures 1 and 2 shows the mean-subtracted and masked SWIRE images at 1.1 and 1.8 μm, respectively. To highlight the large-scale fluctuations, we smooth the images with a σ = 35'' Gaussian kernel.

Although the ZL signal is smooth, a flat-field estimation error may induce a non-uniform ZL residual that cannot be removed by mean subtraction. This residual may dominate over cosmological fluctuations on large scales. Therefore, after removing the mean value in the image, we filter the images by fitting and subtracting a third/fifth-order 2D polynomial function for the 1.1/1.8 μm images to filter out any residual large-scale variations (Panel H of Figures 1 and 2). The filtering will also suppress large-scale cosmological signals and therefore the choice of polynomial order used for filtering is determined by optimizing the trade-off between the reduction of background fluctuations and the large-scale two-halo signal. The effect of filtering on the detected one-halo and galaxy extension terms is small, as our filtering removes fluctuations at much larger scales than these signals and the signal filtering is accounted for in simulations (see Section 9).

We use the Pan-STARRS catalog (Chambers et al. 2016) for masking. Pan-STARRS covers all of the CIBER fields with a depth of m ∼ 20 in the g, r, i, z, y bands. We query the source positions and magnitudes in all five Pan-STARRS bands from their DR1 MeanObject table and derive m_1.1 and m_1.8 with the LePhare SED-fitting software (Arnouts et al. 1999; Ilbert et al. 2006). We use sources that have a y-band measurement and a quality flag (qualityFlag in ObjectThin table) that equals to 8 or 16 for masking.

Some bright stars are not included in the Pan-STARRS catalog, thus we use the 2MASS (Skrutskie et al. 2006) Point Source Catalog (PSC) to get the complete point-source list. For 2MASS sources, m_1.1 (m_1.8) is derived by linear extrapolation with the 2MASS photometric fluxes in the J and H (H and K_s ) bands, respectively. We also use bright stars in 2MASS for modeling the PSF (see Section 7).

We use the Sloan Digital Sky Survey (SDSS) DR13 (Blanton et al. 2017) PhotoObj catalog to get the star/galaxy classification ("type" attribute 6–stars, 3–galaxies) and the galaxy photometric redshift ("Photoz" attribute) for sources in our fields. This information is essential for selecting target galaxies for stacking and inferring their redshift distribution (Section 8.1), as well as selecting stars for stacking to model the PSF (Section 7).

Rowan-Robinson et al. (2008, 2013) performed SED fitting on ∼10⁶ sources in the SWIRE field, based on optical and infrared photometric data from multiple surveys. This provides information on the stellar masses of our stacked galaxies for our analysis (see Section 8.2).

Gaia DR2 (Gaia Collaboration et al. 2016, 2018) provides high-precision astrometry for stars in the Milky Way, which gives high-purity star samples used for both validating the PSF model (Section 7.2) and cleaning residual stars in the galaxy sample selected by SDSS (Section 8.1).

Nearby galaxy clusters along the line of sight introduce extended emission in stacking, so we exclude galaxies that are close to nearby clusters (Section 8.1). We use the cluster catalog from Wen et al. (2012), which compiles 0.05 ≤ z < 0.8 galaxy clusters detected in SDSS-III (Aihara et al. 2011). We also use the Abell cluster samples (Abell 1958) for local galaxy clusters. There are 7 Abell clusters and ∼200 clusters from Wen et al. (2012) over the five CIBER fields.

CIBER imager pixels under-sample the PSF and therefore the surface brightness profile of individual sources is poorly resolved. However, given external source catalogs with high astrometric accuracy, we can stack on a sub-pixel basis and reconstruct the average source profile at scales finer than the native pixel size. This "sub-pixel stacking" technique has been used in previous CIBER imager analyses (Bock et al. 2013; Zemcov et al. 2014) and further investigated recently in the context of optimal photometry (Symons et al. 2021). We summarize the sub-pixel stacking procedure as follows:

• 1.  
  Select a list of stacking target sources from external catalogs;
• 2.  
  Re-grid each pixel into N_sub × N_sub sub-pixels (we use N_sub = 10 in this work). The intensities of all sub-pixels are assigned to the same value as the native pixel without interpolation;
• 3.  
  For each source, unmask pixels associated with its source mask. Pixels masked due to nearby sources or from the instrument mask remain masked;
• 4.  
  Crop an N_size × N_size (at sub-pixel resolution) stamp centered on the target source. We choose N_size = 2401 in this work, which corresponds to a 28' × 28' stamp;
• 5.  
  Repeat steps 3 and 4 for all target sources, average the stamps, and return the final stacked 2D image Σ_stack( r ).

The stacked profile Σ_stack is a convolution of the intrinsic source profile, Σ_src, the instrument PSF (PSF_instr), ¹³ and the pixel function PSF_pix:

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Sigma }}}_{\mathrm{stack}}({\boldsymbol{r}}) & = & \left[{{\rm{\Sigma }}}_{\mathrm{src}}({\boldsymbol{r}}) \circledast {{PSF}}_{\mathrm{instr}}({\boldsymbol{r}})\right] \circledast {{PSF}}_{\mathrm{pix}}({\boldsymbol{r}})\\ & = & {{\rm{\Sigma }}}_{\mathrm{src}}({\boldsymbol{r}}) \circledast {{PSF}}_{\mathrm{stack}}({\boldsymbol{r}}),\end{array}\end{eqnarray} \tag{ 1 }$$

where r = (x, y) is a two-dimensional sub-pixel coordinate system with its origin at the stack center. We define the effective PSF as ${{PSF}}_{\mathrm{stack}}({\boldsymbol{r}})\equiv {{PSF}}_{\mathrm{instr}}({\boldsymbol{r}}) \circledast {{PSF}}_{\mathrm{pix}}({\boldsymbol{r}})$. The pixel function accounts for the fact that sub-pixels retain the value of the original pixels, which is a convolution effect. The pixel function is a matrix with each element proportional to the counts where the sub-pixel and the center sub-pixel that contains the source are within the same native pixel. The position of the center sub-pixel within the native pixel is a uniform probability distribution, and therefore when stacking on a large number of sources, the pixel function converges to the analytic form (Symons et al. 2021):

$$\begin{eqnarray}&&{{PSF}}_{\mathrm{pix}}({\boldsymbol{r}})=\left\{\begin{array}{c}({N}_{\mathrm{sub}}-x)({N}_{\mathrm{sub}}-y)\\ \quad \quad \quad \quad \mathrm{if}\ | x| ,| y| \lt {N}_{\mathrm{sub}}\\ 0\quad \quad \mathrm{otherwise}\end{array}\right.\end{eqnarray} \tag{ 2 }$$

As a practical matter, PSF_pix can be determined through simulations. PSF_stack( r ) can be measured by stacking stars in the field, where Σ_src( r ) is a delta function, so Σ_stack( r ) = PSF_stack( r ). Note that the expression in the second line of Equation (1) implies that the intrinsic profile Σ_src( r ) can be obtained from the stacked profile Σ_stack( r ) with the knowledge of PSF_stack( r ), instead of determining ${{PSF}}_{\mathrm{instr}}({\boldsymbol{r}})$.

We perform stacking and PSF modeling separately for each field, since PSF_instr is slightly different across the fields due to the varying pointing performance of the altitude control system during each integration (see top panel of Figure 5). After obtaining the 2D stacked images, we bin them into 25 logarithmically spaced 1D radial bins. Within each bin, the number of stacked images on each sub-pixel is used for weighting when calculating the average profile in each radial bin. Note that the weight is not the same across sub-pixels since the masks are different for each stacked image.

The covariance matrix of the binned 1D radial stacked profile is calculated with a jackknife resampling technique. For each stack, we split sources into N_J = 64 sub-groups based on their spatial coordinates in the image. The CIBER imager arrays have 1024 × 1024 pixels, thus each sub-group corresponds to sources in a 128 × 128 pixel sub-region on the array. The radial profile of the kth jackknife sample, ${{\rm{\Sigma }}}_{\mathrm{stack}}^{k}$, is obtained from stacking on sources in all the other sub-regions, and then the covariance matrix between radial bin (r_i , r_j ) is given by

$$\begin{eqnarray}\begin{array}{rcl}{C}_{\mathrm{stack}}({r}_{i},{r}_{j}) & = & \displaystyle \frac{{N}_{{\rm{J}}}-1}{{N}_{{\rm{J}}}}\displaystyle \sum _{k=1}^{{N}_{{\rm{J}}}}\left[{\rm{\Delta }}{{\rm{\Sigma }}}_{\mathrm{stack}}^{k}({r}_{i})\cdot {\rm{\Delta }}{{\rm{\Sigma }}}_{\mathrm{stack}}^{k}({r}_{j})\right]\\ {\rm{\Delta }}{{\rm{\Sigma }}}_{\mathrm{stack}}^{k}({r}_{i}) & \equiv & {{\rm{\Sigma }}}_{\mathrm{stack}}^{k}({r}_{i})-{{\rm{\Sigma }}}_{\mathrm{stack}}({r}_{i})\\ {\rm{\Delta }}{{\rm{\Sigma }}}_{\mathrm{stack}}^{k}({r}_{j}) & \equiv & {{\rm{\Sigma }}}_{\mathrm{stack}}^{k}({r}_{j})-{{\rm{\Sigma }}}_{\mathrm{stack}}({r}_{j}),\end{array}\end{eqnarray} \tag{ 3 }$$

where Σ_stack is the average stacked profile of all of the sub-regions.

One of our galaxy stacking samples (mag bin # 1 in Section 8.1) has a small number of sources (≪64 for each field), which makes the covariance estimation from the jackknife method unstable. Therefore we perform bootstrap resampling with N_B = 1000 realizations to calculate the covariance for this case. In this bootstrap, we obtain the radial profile of the kth bootstrap sample, ${{\rm{\Sigma }}}_{\mathrm{stack}}^{k}$, by stacking the same number of sources as the original sample, but the sources are randomly selected from the original sample with replacement. The covariance matrix is then given by

$$\begin{eqnarray}&&{C}_{\mathrm{stack}}({r}_{i},{r}_{j})=\displaystyle \frac{1}{{N}_{{\rm{B}}}-1}\sum _{k=1}^{{N}_{{\rm{B}}}}\left[{\rm{\Delta }}{{\rm{\Sigma }}}_{\mathrm{stack}}^{k}({r}_{i})\cdot {\rm{\Delta }}{{\rm{\Sigma }}}_{\mathrm{stack}}^{k}({r}_{j})\right].\end{eqnarray} \tag{ 4 }$$

In all the other cases, the covariance is derived from jackknife instead of bootstrap resampling since it is numerically expensive to perform a sufficient number of bootstrap realizations given that we have hundreds or thousands of galaxies per field in each stack. We assign galaxies to sub-groups by their spatial positions instead of randomly grouping them to account for large-scale spatial fluctuations.

The first few radial bins within the CIBER 7'' native pixel are highly correlated since all the sub-pixels are assigned to the same value as the native pixel. We also find a high correlation on large angular scales, as the stacking signal is dominated by large-scale spatial variations.

Infrared detectors have a brightness-dependent PSF, the so-called "brighter-fatter effect" (Hirata & Choi 2020). This nonlinearity makes brighter point sources appear broader on the detector array than fainter ones. To model PSF_stack robustly on both small and large scales, we construct an overall star profile from three brightness bins. For the core region (r < 22''), we stack 13 < m_1.1 < 14 sources in the field; for intermediate scales, 22'' < r < 40'', we fit a slope to the stacking profile of $9\lt {m}_{J}^{2\mathrm{MASS}}\lt 10$ sources; for outer radii, we fit another slope to the stacking profile of the brightest $4\lt {m}_{J}^{2\mathrm{MASS}}\lt 9$ sources and connect the two slopes at r = 40'' (${m}_{J}^{2\mathrm{MASS}}$ is the 2MASS J-band Vega magnitude). The choice of magnitude bins and transition radii minimizes the error on all scales. At small radii, using faint stars avoids detector nonlinearity and at large radii, bright stars provide better sensitivity to the extended PSF. For the intermediate scales, we check that the fitted slope from the three star-stacking profiles ($4\lt {m}_{J}^{2\mathrm{MASS}}\lt 9$, $9\lt {m}_{J}^{2\mathrm{MASS}}\lt 10$, 13 < m_1.1 < 14) are statistically consistent. The top panel of Figure 3 shows PSF_stack from the SWIRE field in the 1.1 μm band. The top panel of Figure 5 shows PSF_stack in all five fields in both bands. The slight variation across fields is due to the difference in the pointing stability during each integration, but such motion is common to all sources within an integration.

Zoom In Zoom Out Reset image size

To validate that our PSF model is applicable to the fainter sources of interest, we perform a consistency test by stacking on stars in the Gaia catalog within the same magnitude range as our stacked galaxy samples (16 < m_1.1 < 20) and compare these star-stacking profiles with our PSF_stack model.

To get a clean star sample free of galaxies, we apply the following criteria for selecting stars from Gaia:

• 1.  
  The source has a parallax measurement >2 × 10⁻⁴ mas (i.e., distance <5 kpc);
• 2.  
  No astrometric excess noise is reported in the Gaia catalog (astrometric_excess_noise = 0). Large astrometric excess noise implies the source might be extended rather than a point source;
• 3.  
  No SDSS galaxies within 07 (sub-pixel grid size) radius around the source;
• 4.  
  We classify SDSS stars and galaxies using 10 pairs of magnitude differences between the 5 Pan-STARRS photometric magnitudes (g, r, i, z, and y bands), rejecting sources if they are classified as galaxies by our trained model.

After selecting stars with the above conditions from the Gaia catalog, we stack them in four equally-spaced magnitude bins between 16 < m_1.1 < 20, and compare their stacking profile with the PSF_stack model. These stars span the same brightness range used for galaxy stacking. We down-sample the original 25 radial bins to 15 bins (7 bins for 16 < m_1.1 < 17 case), following the same binning used for the galaxy stacking profile (Section 8.4). The results in the 1.1 μm band SWIRE field are shown on the bottom panel of Figure 3. In Figure 4 we show the difference of Gaia star stacks and the PSF_stack model. The errors are propagated from the covariance of the PSF_stack model and Gaia star stacks. We also show the χ² values and the corresponding probability to exceed (PTE) on all five CIBER fields in both bands. The PSF model shows excellent agreement with the star stacks.

Zoom In Zoom Out Reset image size

Although knowledge of the instrument PSF is not required for reconstructing the source profile Σ_src from the stacking profile Σ_stack, PSF_instr is still needed when we model the clustering signal from a simulated catalog (Section 9), where we make mock galaxy images using the CIBER PSF and pixel gridding. PSF_instr is also useful for determining the masking radius for bright sources (Section 3.3.2).

PSF_instr is modeled as follows: first, we deconvolve PSF_pix(r) (Equation (2)) from the PSF_stack( r ) model with 10 iterations of the Richardson–Lucy deconvolution algorithm (Richardson 1972; Lucy 1974). The deconvolution is unstable at large radii due to noise fluctuations. To get a smooth model for PSF_instr, we fit a β model (Cavaliere & Fusco-Femiano 1978) to the 1D profile of the deconvolved image:

$$\begin{eqnarray}&&{{PSF}}_{\mathrm{instr}}(r)={\left(1+{\left(\displaystyle \frac{r}{{r}_{c}}\right)}^{2}\right)}^{-3\beta /2}.\end{eqnarray} \tag{ 5 }$$

Though not physically motivated, we find the β model is a good empirical description of the extended PSF and requires only two free parameters to achieve acceptable goodness of fit for every PSF_stack.

The bottom panel of Figure 5 illustrates this procedure in the 1.1 μm band of the SWIRE field. The PSF_stack model, obtained from star stacks in three different brightness bins, matches the β model of PSF_instr convolved with the pixel function PSF_pix (Equation (2)). Our instrument PSF has a size comparable to a pixel (FWHM ∼ 7'').

Zoom In Zoom Out Reset image size

The stacking galaxy samples are selected from the SDSS catalog in the CIBER fields. To mitigate systematic effects from confusion, nearby clusters, or misclassified stars in the sample, we reject sources if they meet any of the following criteria:

• 1.  
  Sources are not labeled as galaxies in the SDSS catalog, i.e., the "type" attribute in the SDSS PhotoObj table is not equal to 3;
• 2.  
  Sources are located in the instrument mask;
• 3.  
  Other Pan-STARRS sources exist in the same CIBER pixel;
• 4.  
  The SDSS photometric redshift is less than 0.15. These criteria prevent nearby galaxies from introducing substantial power on large angular scales that would otherwise mimic the clustering signal;
• 5.  
  Sources have nearby Gaia counterparts within 0.7'', i.e., the size of the sub-pixel used in our stacking. These sources are likely to be stars that are misclassified as galaxies in the SDSS catalog;
• 6.  
  Sources are within (1) a 500" radius of any galaxy cluster in Abell (1958) (Section 4.6); or (2) R₂₀₀ of any galaxy cluster with halo mass M_h > 10¹⁴ M_⊙ or redshift z < 0.15 in the SDSS cluster catalog (Wen et al. 2012, Section 4.6). Approximately 10% of the sky area in each field is excluded by this condition.

The last condition mitigates contamination from nearby clusters along the line of sight since they have structures spanning large angular scales, which will produce spurious large-scale extended signals in the stack. Furthermore, as we do not have information on whether a galaxy in SDSS is a member of a large galaxy cluster, the criteria also exclude cluster members from our stacking sample. Stacking on cluster members introduces extra nonlinear one-halo clustering that can overwhelm the linear two-halo clustering signal on large scales.

To quantify the effect of applying this condition, we generate a mock CIBER map from the MICECAT catalog, implementing the same strategies described above to select sources and stacking on the mock maps to measure the one- and two-halo clustering signals (see Section 9 for a detailed description of stacking with MICECAT-generated maps). We tested over a range of halo mass and redshift for selecting clusters and found that excluding sources around clusters with M_h > 10¹⁴ M_⊙ (or redshift z < 0.15) can effectively reduce the one-halo clustering signal on large scales without losing a significant number of sources. For example, for the magnitude range of interest in this work (see Section 8.2), we can reduce the one-halo power by ∼3–5× at 100'' radius just by excluding galaxies near clusters following our criteria.

For the SDSS galaxies within 16 < m_1.1 < 20 that survive all the selection criteria above, we split the sources into two sets. The first set is based on 1.1 μm flux in four bins: 16 < m_1.1 < 17, 17 < m_1.1 < 18, 18 < m_1.1 < 19, and 19 < m_1.1 < 20. Hereafter, these four bins are named "mag bin # 1," "mag bin # 2," "mag bin # 3," and "mag bin # 4," respectively. In addition, we also define a "total stack" with all 17 < m_1.1 < 20 sources to achieve better large-scale sensitivity.

The second set is defined by both the 1.1 μm apparent magnitude m_1.1 and the absolute magnitude M_1.1: ${M}_{1.1}={m}_{1.1}\,-{DM}(z)+2.5{\mathrm{log}}_{10}(1+z)$, where DM is the distance modulus, using SDSS photometric redshifts. The absolute flux serves as a proxy for galaxy size. Galaxies with comparable absolute flux have similar bolometric luminosity, which is correlated with stellar mass, star formation rate, etc. We use these samples to explore the dependence of our results on different galaxy properties. Since the sets approximately correspond to three higher and two lower stellar mass populations, with different redshift distributions, we call them "high-M/low-z," "high-M/med-z," "high-M/high-z," "low-M/low-z," and "low-M/med-z."

In the SWIRE field, we have additional information from a photometric redshift catalog (Rowan-Robinson et al. 2013) based on an SED fit to each galaxy. As the stacked samples from each field are selected with the same criteria, we can assume the galaxy property distributions in the SWIRE field are the same as other fields, and thus infer the stellar mass distribution over all five fields. The log M_* column in Table 2 lists the median and 68% interval stellar mass in the SWIRE field samples from the Rowan-Robinson et al. (2013) catalog. The stellar masses of our samples span from ∼10^10.5 to 10¹² M_⊙, i.e., ∼L_* galaxies at this redshift (Muzzin et al. 2013). In addition, with the stellar mass distribution, we infer the host halo mass of our samples using the mean stellar-to-halo mass relation given by Zu & Mandelbaum (2015), which connects the halo mass to stellar mass with galaxy clustering and lensing measurements. We also derive the corresponding virial radius, R₂₀₀ (in physical and angular units), in Table 2. The virial radius is calculated from ${R}_{200}={[3{M}_{h}/(4\pi \cdot 200{\rho }_{c})]}^{1/3}$, where ρ_c is the critical density.

Table 2. Summary of the Properties on Each Stacked Galaxy Sub-sample with +/- Values Indicating the 68% Interval Ranges

Name	Selection Criteria	Ngal	z	log M* (M⊙)	log Mh (M⊙)	R200 (kpc)	R200 (arcsec)	fcen
mag bin #1	16 < m1.1 < 17	129	${0.18}_{-0.02}^{+0.04}$	${11.6}_{-0.3}^{+0.3}$	${13.8}_{-0.4}^{+0.5}$	${679}_{-181}^{+325}$	${215}_{-57}^{+103}$	0.65
mag bin #2	17 < m1.1 < 18	1173	${0.21}_{-0.04}^{+0.07}$	${11.5}_{-0.4}^{+0.3}$	${13.7}_{-0.6}^{+0.6}$	${584}_{-215}^{+357}$	${163}_{-60}^{+100}$	0.67
mag bin #3	18 < m1.1 < 19	3465	${0.27}_{-0.07}^{+0.09}$	${11.2}_{-0.3}^{+0.4}$	${13.3}_{-0.4}^{+0.5}$	${401}_{-116}^{+178}$	${94}_{-27}^{+42}$	0.62
mag bin #4	19 < m1.1 < 20	31157	${0.42}_{-0.11}^{+0.17}$	${11.1}_{-0.5}^{+0.3}$	${13.0}_{-0.5}^{+0.5}$	${285}_{-86}^{+127}$	${50}_{-15}^{+22}$	0.63
total	17 < m1.1 < 20	35795	${0.40}_{-0.14}^{+0.17}$	${11.1}_{-0.4}^{+0.3}$	${13.1}_{-0.5}^{+0.5}$	${302}_{-93}^{+135}$	${55}_{-17}^{+24}$	0.63
high-M/low-z	17 < m1.1 < 18 − 23 < M1.1 < − 22	743	${0.22}_{-0.03}^{+0.04}$	${11.6}_{-0.4}^{+0.2}$	${13.7}_{-0.5}^{+0.5}$	${608}_{-201}^{+266}$	${168}_{-55}^{+73}$	0.66
high-M/med-z	18 < m1.1 < 19 − 23 < M1.1 < − 22	1274	${0.34}_{-0.05}^{+0.05}$	${11.4}_{-0.2}^{+0.3}$	${13.5}_{-0.3}^{+0.4}$	${447}_{-94}^{+157}$	${89}_{-19}^{+31}$	0.62
high-M/high-z	19 < m1.1 < 20 − 23 < M1.1 < − 22	10916	${0.54}_{-0.09}^{+0.10}$	${11.3}_{-0.3}^{+0.3}$	${13.4}_{-0.4}^{+0.3}$	${325}_{-82}^{+100}$	${50}_{-13}^{+15}$	0.66
low-M/low-z	18 < m1.1 < 19 − 22 < M1.1 < − 21	1645	${0.24}_{-0.03}^{+0.05}$	${11.1}_{-0.2}^{+0.3}$	${13.1}_{-0.3}^{+0.4}$	${359}_{-78}^{+129}$	${90}_{-20}^{+33}$	0.57
low-M/med-z	19 < m1.1 < 20 − 22 < M1.1 < − 21	14730	${0.38}_{-0.05}^{+0.05}$	${11.0}_{-0.4}^{+0.2}$	${12.9}_{-0.4}^{+0.3}$	${275}_{-67}^{+78}$	${51}_{-13}^{+15}$	0.58

Note. N_gal is the total number of galaxies across five CIBER fields in each sub-sample and the redshifts z are derived from SDSS photometry. The quantities on the left side of the double vertical line are derived from a partial set of samples or external catalogs for the sources used in stacks. We infer M_* by matching SWIRE field sources to the catalog from Rowan-Robinson et al. (2013), assuming the same M_* distribution applies to the other four fields. The halo mass and the virial radius are derived with the stellar-to-halo mass relation from Zu & Mandelbaum (2015). The fraction of central galaxies (f_cen) is derived by applying the same cuts to a simulated catalog from MICECAT.

Download table as:  ASCIITypeset image

We note that by selecting galaxies based on absolute or apparent fluxes, our samples will include both central and satellite galaxies. We infer the fraction of central galaxies, f_cen, in each sub-sample from MICECAT by applying the same selection criteria from a MICECAT simulation (i.e., observed magnitude, absolute magnitude and redshift cuts and excluding sources close to nearby clusters). The distribution of redshift, stellar mass, halo mass, virial radius, and f_cen of our sub-samples is summarized in Figure 6 and Table 2.

Zoom In Zoom Out Reset image size

We calculate 1D radial profiles from galaxy stacks by averaging pixels in concentric annuli, as shown in Figures 7 and 8. For comparison, we also plot the expected profile of stacked point sources, PSF_stack, scaled to match the first radial bin of the stacked galaxy profile. In all cases, the galaxy profiles are clearly broader than the PSF_stack profile.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We define an "excess profile" Σ_ex(r) as follows:

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{ex}}(r)={{\rm{\Sigma }}}_{\mathrm{stack}}(r)-A\cdot {{PSF}}_{\mathrm{stack}}(r),\end{eqnarray} \tag{ 6 }$$

where the normalization factor A is chosen such that PSF_stack matches Σ_stack at the innermost radial bin r₁, and thus by construction, Σ_ex(r₁) = 0, and A ≡ Σ_stack(r₁)/PSF_stack(r₁).

Since the excess profile is fixed at r₁, the uncertainties on the galaxy profile and the PSF profile at r₁ have to be accounted for by propagating this error to the other radial bins, and thus the excess profile covariance is given by

$$\begin{eqnarray}&&{C}_{\mathrm{ex}}={{\rm{\Sigma }}}_{\mathrm{stack}}{\left({r}_{1}\right)}^{2}\left[{C}_{\mathrm{norm}}\left({C}_{\mathrm{stack}}\right)+{C}_{\mathrm{norm}}\left({C}_{\mathrm{PSF}}\right)\right],\end{eqnarray} \tag{ 7 }$$

where C_PSF and C_stack are the covariance of PSF_stack and Σ_stack, respectively, and

$$\begin{eqnarray}&&\begin{array}{l}{C}_{\mathrm{norm}}\left(C,\displaystyle \frac{{{\rm{\Sigma }}}_{\mathrm{stack}}({r}_{i})}{{{\rm{\Sigma }}}_{\mathrm{stack}}({r}_{1})},\displaystyle \frac{{{\rm{\Sigma }}}_{\mathrm{stack}}({r}_{j})}{{{\rm{\Sigma }}}_{\mathrm{stack}}({r}_{1})}\right)\\ =\ \displaystyle \frac{{{\rm{\Sigma }}}_{\mathrm{stack}}({r}_{i}){{\rm{\Sigma }}}_{\mathrm{stack}}({r}_{j})}{{{\rm{\Sigma }}}_{\mathrm{stack}}{\left({r}_{1}\right)}^{2}}\cdot \\ \left[\displaystyle \frac{C({r}_{i},{r}_{j})}{{{\rm{\Sigma }}}_{\mathrm{stack}}({r}_{i}){{\rm{\Sigma }}}_{\mathrm{stack}}({r}_{j})}-\displaystyle \frac{C({r}_{i},{r}_{1})}{{{\rm{\Sigma }}}_{\mathrm{stack}}({r}_{i}){{\rm{\Sigma }}}_{\mathrm{stack}}({r}_{1})}\right.\\ \left.-\ \displaystyle \frac{C({r}_{j},{r}_{1})}{{{\rm{\Sigma }}}_{\mathrm{stack}}({r}_{j}){{\rm{\Sigma }}}_{\mathrm{stack}}({r}_{1})}+\displaystyle \frac{C({r}_{1},{r}_{1})}{{{\rm{\Sigma }}}_{\mathrm{stack}}{\left({r}_{1}\right)}^{2}}\right]\end{array}\end{eqnarray} \tag{ 8 }$$

is the covariance for the normalized profile that follows from the product rule for derivatives.

To fit a model to the measured Σ_ex, we also need the inverse of C_ex. However, C_ex is close to singular since our radial bins are highly correlated. Therefore, we reduce the original 25 radial bins to 15 bins by combining highly correlated bins in the inner and outer regions. ¹⁴ After this down-sampling, we derive the inverse covariance estimator by

$$\begin{eqnarray}&&{C}_{\mathrm{ex}}^{-1}=\displaystyle \frac{{N}_{{\rm{J}}}-{N}_{\mathrm{bin}}-2}{{N}_{{\rm{J}}}-1}{C}_{\mathrm{ex}}^{* -1},\end{eqnarray} \tag{ 9 }$$

where N_J = 64, the number of sub-groups used for estimating covariance, and the number of bins N_bin = 15. ${C}_{\mathrm{ex}}^{* -1}$ is the direct inverse of the C_ex matrix, and the pre-factor in Equation (9) de-biases the inverse covariance estimator, as our covariance matrix is derived from our data (Hartlap et al. 2007). ¹⁵

While we have high sensitivity on the small radial bins of both the galaxy stacked profiles and the PSF model, the A value has minimal dependency on the radius chosen for normalization, and the uncertainty of normalization has been accounted for by the covariance (Equation (7)), thus our model parameter inference (Section 9) does not depend on the definition of the excess profile.

We present field-averaged excess profiles in Figure 9. Note that the field-averaged excess profile is only plotted for visualization purposes since the field-to-field PSF variation must be explicitly accounted for in parameter fitting.

Zoom In Zoom Out Reset image size

The raw CIBER image, I_raw, can be expressed as ¹⁶

$$\begin{eqnarray}\begin{array}{rcl}{I}_{\mathrm{raw}}({\boldsymbol{x}}) & = & \left[{I}_{\mathrm{sig}}({\boldsymbol{x}})+{I}_{\mathrm{LoS}}({\boldsymbol{x}})\right] \circledast {{PSF}}_{\mathrm{instr}}({\boldsymbol{r}})\cdot {FF}({\boldsymbol{x}})\\ & & +{I}_{\mathrm{DC}}({\boldsymbol{x}})+{I}_{{\rm{n}}}({\boldsymbol{x}}),\end{array}\end{eqnarray} \tag{ 10 }$$

where x is the 2D pixel coordinate, FF is the flat-field gain, I_DC is the dark-current map, and I_n is the read noise plus photon noise. The sky emission is decomposed into I_sig and I_LoS terms, where the first term accounts for the signal associated with stacked galaxies and I_LoS represents uncorrelated emission from all other sources along the line of sight, including Galactic foregrounds.

After dark-current subtraction and flat-field correction, we retrieve $I{{\prime} }_{\mathrm{raw}}$:

$$\begin{eqnarray}&&I{{\prime} }_{\mathrm{raw}}({\boldsymbol{x}})=\left[{I}_{\mathrm{sig}}({\boldsymbol{x}})+{I}_{\mathrm{LoS}}({\boldsymbol{x}})\right] \circledast {{PSF}}_{\mathrm{instr}}({\boldsymbol{r}})+I{{\prime} }_{{\rm{n}}}({\boldsymbol{x}}),\end{eqnarray} \tag{ 11 }$$

where $I{{\prime} }_{{\rm{n}}}({\boldsymbol{x}})={I}_{{\rm{n}}}({\boldsymbol{x}})/{FF}({\boldsymbol{x}})$, the instrument noise divided by the flat-field response. For simplicity, we ignore the error in the flat-field estimator in Equation (11). In practice, the flat-field estimation uncertainties will not bias the stacking results as they are not correlated with individual stacked sources and the effect on the covariance is accounted for by the jackknife method (see Section 6.2). We define the mask M( x ) as a binary function set to zero for masked pixels and one otherwise. The filtered map is expressed with ${ \mathcal F }\left[I{{\prime} }_{\mathrm{raw}}({\boldsymbol{x}}),M({\boldsymbol{x}})\right]$, which is a function of the input map $I{{\prime} }_{\mathrm{raw}}({\boldsymbol{x}})$ and mask M( x ). As described in Section 3.7, we choose ${ \mathcal F }$ to be a third (1.1 μm)/fifth (1.8 μm)-order 2D polynomial function fitted to the masked $I{{\prime} }_{\mathrm{raw}}$ map. ¹⁷ The image used for stacking I_map can thus be written as

$$\begin{eqnarray}\begin{array}{rcl}{I}_{\mathrm{map}}({\boldsymbol{x}}) & = & I{{\prime} }_{\mathrm{raw}}({\boldsymbol{x}})M({\boldsymbol{x}})-{ \mathcal F }\left[I{{\prime} }_{\mathrm{raw}}({\boldsymbol{x}}),M({\boldsymbol{x}})\right]M({\boldsymbol{x}})\\ & = & {I}_{\mathrm{map}}^{\mathrm{sig}}({\boldsymbol{x}})+{I}_{\mathrm{map}}^{\mathrm{LoS}}({\boldsymbol{x}})+{I}_{\mathrm{map}}^{{\rm{n}}^{\prime} }({\boldsymbol{x}}),\end{array}\end{eqnarray} \tag{ 12 }$$

where

$$\begin{eqnarray}&&\begin{array}{l}{I}_{\mathrm{map}}^{\mathrm{sig}}({\boldsymbol{x}})=\left[{I}_{\mathrm{sig}}({\boldsymbol{x}}) \circledast {{PSF}}_{\mathrm{instr}}({\boldsymbol{r}})\right]M({\boldsymbol{x}})\\ -\ { \mathcal F }\left[{I}_{\mathrm{sig}}({\boldsymbol{x}}) \circledast {{PSF}}_{\mathrm{instr}}({\boldsymbol{r}}),M({\boldsymbol{x}})\right]M({\boldsymbol{x}}),\end{array}\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&\begin{array}{l}{I}_{\mathrm{map}}^{\mathrm{LoS}}({\boldsymbol{x}})=\left[{I}_{\mathrm{LoS}}({\boldsymbol{x}}) \circledast {{PSF}}_{\mathrm{instr}}({\boldsymbol{r}})\right]M({\boldsymbol{x}})\\ -\ { \mathcal F }\left[{I}_{\mathrm{LoS}}({\boldsymbol{x}}) \circledast {{PSF}}_{\mathrm{instr}}({\boldsymbol{r}}),M({\boldsymbol{x}})\right]M({\boldsymbol{x}}),\end{array}\end{eqnarray} \tag{ 14 }$$

and

$$\begin{eqnarray}&&\begin{array}{l}{I}_{\mathrm{map}}^{{\rm{n}}^{\prime} }({\boldsymbol{x}})={I}_{{\rm{n}}}^{\prime} ({\boldsymbol{x}})M({\boldsymbol{x}})\\ -\ { \mathcal F }\left[{I}_{{\rm{n}}}^{\prime} ({\boldsymbol{x}}),M({\boldsymbol{x}})\right]M({\boldsymbol{x}}).\end{array}\end{eqnarray} \tag{ 15 }$$

The stacked profile Σ_stack can be expressed as the sum of stacks on the three maps in Equation (12):

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{stack}}(r)={{\rm{\Sigma }}}_{\mathrm{stack}}^{\mathrm{sig}}(r)+{{\rm{\Sigma }}}_{\mathrm{stack}}^{\mathrm{LoS}}(r)+{{\rm{\Sigma }}}_{\mathrm{stack}}^{{\rm{n}}}(r).\end{eqnarray} \tag{ 16 }$$

The last two terms can be ignored in modeling since they are uncorrelated with the stacked sources, so $\left\langle {{\rm{\Sigma }}}_{\mathrm{stack}}^{\mathrm{LoS}}(r)\right\rangle \,=\left\langle {{\rm{\Sigma }}}_{\mathrm{stack}}^{{\rm{n}}}(r)\right\rangle =0$.

We model the stacked galaxy profile as

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Sigma }}}_{\mathrm{stack}}^{\mathrm{sig}}(r) & = & \left[\right.{{\rm{\Sigma }}}_{\mathrm{stack}}^{\mathrm{gal}}(r)+{{\rm{\Sigma }}}_{\mathrm{stack}}^{1{\rm{h}}}(r)+{{\rm{\Sigma }}}_{\mathrm{stack}}^{2{\rm{h}}}(r)\left.\right]\\ & & -\ {{\rm{\Sigma }}}_{\mathrm{stack}}^{{ \mathcal F }}(r),\end{array}\end{eqnarray} \tag{ 17 }$$

where the first three terms are the signal terms and the last term is the filtered signal map in Equation (13). The galaxy profile term, ${{\rm{\Sigma }}}_{\mathrm{stack}}^{\mathrm{gal}}$, represents the intrinsic galaxy profile, which includes the galaxy shape and the extended stellar halo. We decompose the galaxy profile term, ${{\rm{\Sigma }}}_{\mathrm{stack}}^{\mathrm{gal}}$, into "core" and "extended" parts:

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{stack}}^{\mathrm{gal}}(r)={{\rm{\Sigma }}}_{\mathrm{core}}^{\mathrm{gal}}(r)+{{\rm{\Sigma }}}_{\mathrm{ext}}^{\mathrm{gal}}(r),\end{eqnarray} \tag{ 18 }$$

where the core component is the integrated emission of the PSF_stack fitted to the stacking profile, i.e., the A · PSF_stack term in Equation (6) and the extended component is the rest of the galaxy emission:

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Sigma }}}_{\mathrm{core}}^{\mathrm{gal}}(r) & = & {{\rm{\Sigma }}}_{\mathrm{stack}}(r)-{{\rm{\Sigma }}}_{\mathrm{ex}}(r),\\ {{\rm{\Sigma }}}_{\mathrm{ext}}^{\mathrm{gal}}(r) & = & {{\rm{\Sigma }}}_{\mathrm{ex}}(r)\\ & & -\ \left[{{\rm{\Sigma }}}_{\mathrm{stack}}^{1{\rm{h}}}(r)+{{\rm{\Sigma }}}_{\mathrm{stack}}^{2{\rm{h}}}(r)-{{\rm{\Sigma }}}_{\mathrm{stack}}^{{ \mathcal F }}(r)\right].\end{array}\end{eqnarray} \tag{ 19 }$$

In addition, galaxy clustering will also contribute to the stacked profile, primarily on large scales. We model clustering with the halo model framework (Cooray & Sheth 2002), where large-scale clustering is described by the correlation within (one-halo) and between (two-halo) dark matter halos. ${{\rm{\Sigma }}}_{\mathrm{stack}}^{1{\rm{h}}}$ and ${{\rm{\Sigma }}}_{\mathrm{stack}}^{2{\rm{h}}}$ represent the profile for one- and two-halo clustering, respectively.

In practice, there is no well-defined boundary between the stellar halo of a galaxy and unbound stars in the dark matter halo and the definition of IHL (or ICL) varies in the literature. To some degree, the galaxy extension term and the one-halo term each partially comprise stars not bound to individual galaxies in the halo. Since there are different definitions of IHL (or ICL) and the one-halo term in the literature, here we describe how our modeled components are defined.

In our definition, the galaxy extension describes emission associated with each galaxy, whereas the one-halo term accounts for other galaxies, their extensions, and diffuse stars in the same halo, as illustrated in Figure 10. When we stack on a central galaxy, the galaxy extension term accounts for the extended emission around the stacked galaxy and the one-halo term describes diffuse stars, undetected galaxies, and extension around all the satellite galaxies beyond the masking limit in the same halo. Whereas, when we stack on a satellite galaxy, the galaxy extension term only includes the extended halo around that satellite galaxy and all the other components are described by the one-halo term. In our sample, we estimate that ∼60% of stacked galaxies are central galaxies and ∼40% are satellite galaxies (see Table 2).

Zoom In Zoom Out Reset image size

The stacked galaxy profile ${{\rm{\Sigma }}}_{\mathrm{stack}}^{\mathrm{gal}}(r)={{\rm{\Sigma }}}^{\mathrm{gal}}(r) \circledast {{PSF}}_{\mathrm{stack}}(r)$ is the intrinsic galaxy profile Σ^gal, including the galaxy shape and the extended stellar halo, convolved with PSF_stack. Following Wang et al. (2019), we model Σ^gal with a double Sérsic function:

$$\begin{eqnarray}&&\begin{array}{l}{{\rm{\Sigma }}}^{\mathrm{gal}}(r)={A}^{\mathrm{gal}}\left({10}^{{I}_{e,1}}\exp \left\{-{b}_{n1}\left[{\left(r/{R}_{e,1}\right)}^{1/{n}_{1}}-1\right]\right\}\right.\\ +\ \left.{10}^{{I}_{e,2}}\exp \left\{-{b}_{n2}\left[{\left(r/\left({R}_{e,1}+{R}_{e,2}\right)\right)}^{1/{n}_{2}}-1\right]\right\}\right).\end{array}\end{eqnarray} \tag{ 20 }$$

Wang et al. (2019) performed a stacking analysis on isolated galaxies from Hyper Suprime-Cam (HSC) images and fitted the stacked profile of their high-concentration samples with this model. The first term captures the galaxy shape and the second term models the extended emission. Due to the lack of angular resolution in CIBER data, we are sensitive to the extended profile, and therefore we only vary R_e,2 to fit our stacked profile. We fix all of the other parameters to the best-fit values given by Table 3 of Wang et al. (2019), although when convolved, the total closely follows the PSF. ¹⁸

Our one- and two-halo clustering models, ${{\rm{\Sigma }}}_{\mathrm{stack}}^{1{\rm{h}}}$ and ${{\rm{\Sigma }}}_{\mathrm{stack}}^{2{\rm{h}}}$, and the filtered signal ${{\rm{\Sigma }}}_{\mathrm{stack}}^{{ \mathcal F }}$, are constructed from the MICECAT simulation. MICECAT includes central and satellite galaxies of each halo and each galaxy has a halo ID, enabling us to decouple the one-halo and two-halo contribution in the stacked signal and thus to take into account the complication that we have both central and satellite galaxies in our samples. We model the one-halo term ${{\rm{\Sigma }}}_{\mathrm{stack}}^{1{\rm{h}}}$ from MICECAT using the following steps:

• 1.  
  Select the stacked target in the catalog using the same selection criteria;
• 2.  
  For each target galaxy, generate a source map (using PSF_instr) for all galaxies residing in the same halo except for the target galaxy;
• 3.  
  Generate a source mask using the same prescription as our data;
• 4.  
  Stack on the target source position;
• 5.  
  Iterate steps (2)–(4) for all target sources.

The derived stacked profile provides our template for the one-halo term, ${T}_{\mathrm{stack}}^{1{\rm{h}}}$. The filtered signal term ${{\rm{\Sigma }}}_{\mathrm{stack}}^{{ \mathcal F }}$ accounts for the loss of clustering signal from filtering. ${{\rm{\Sigma }}}_{\mathrm{stack}}^{{ \mathcal F }}$ is the stacked profile on the 2D polynomial filtered map (the second term of Equation (13)), which can be modeled by filtering the simulated map from MICECAT. We model the two-halo term ${{\rm{\Sigma }}}_{\mathrm{stack}}^{2{\rm{h}}}-{{\rm{\Sigma }}}_{\mathrm{stack}}^{{ \mathcal F }}$ after filtering with the following process:

• 1.  
  Make a CIBER-sized mock image from all the catalog sources with the model PSF_instr, and mask it with a source mask generated using the same masking process applied to the data;
• 2.  
  Fit and subtract a 2D polynomial map to the image;
• 3.  
  Select the stacked target in the catalog using the same selection criteria as the real sources;
• 4.  
  Perform stacking with the target source, subtracting all galaxies within the same halo to remove the target galaxy and the one-halo contribution;
• 5.  
  Iterate on step (4) to derive a stacked profile of the filtered two-halo signal.

The resulting stacked profile, ${T}_{\mathrm{stack}}^{2{\rm{h}}-{ \mathcal F }}$, is a model for ${{\rm{\Sigma }}}_{\mathrm{stack}}^{2{\rm{h}}}-{{\rm{\Sigma }}}_{\mathrm{stack}}^{{ \mathcal F }}$, which provides our template for the two-halo term. This process was performed on 400 realizations with CIBER-sized mock images from MICECAT and we take the average stacked profile as the one-halo and filtered two-halo templates. As diffuse stars and faint galaxies below the resolution limit of MICECAT will not be accounted for, we assign free amplitudes to the one-halo and two-halo templates, which are then fit to the observed stacked data. Therefore, our three-parameter (R_e,2, A_1h, A_2h) model can be written as

$$\begin{eqnarray}&&\begin{array}{l}{{\rm{\Sigma }}}_{\mathrm{stack}}(r,\left\{{R}_{{\rm{e}},2},{A}_{1{\rm{h}}},{A}_{2{\rm{h}}}\right\})\\ =\ {{\rm{\Sigma }}}^{\mathrm{gal}}(r,\left\{{R}_{{\rm{e}},2}\right\}) \circledast {{PSF}}_{\mathrm{stack}}(r)\\ \quad +\ {A}_{1{\rm{h}}}{T}_{\mathrm{stack}}^{1{\rm{h}}}(r)+{A}_{2{\rm{h}}}{T}_{\mathrm{stack}}^{2{\rm{h}}-{ \mathcal F }}(r).\end{array}\end{eqnarray} \tag{ 21 }$$

We note that the one- and two-halo profiles already include the PSF convolution in our model.

For each CIBER field and band, we fit the excess profile Equation (6) to a three-parameter model ${{\rm{\Sigma }}}_{\mathrm{ex}}^{{\rm{m}}}(r,\{{R}_{{\rm{e}},2},{A}_{1{\rm{h}}},{A}_{2{\rm{h}}}\})$ (Equation (21)) using a Markov Chain Monte Carlo (MCMC). We assume a Gaussian likelihood, which is given by

$$\begin{eqnarray}\begin{array}{rcl}{\chi }^{2} & = & {\left({{\rm{\Sigma }}}_{\mathrm{ex}}^{{\rm{d}}}-{{\rm{\Sigma }}}_{\mathrm{ex}}^{{\rm{m}}}\right)}^{T}{C}_{\mathrm{ex}}^{-1}\left({{\rm{\Sigma }}}_{\mathrm{ex}}^{{\rm{d}}}-{{\rm{\Sigma }}}_{\mathrm{ex}}^{{\rm{m}}}\right)\\ \mathrm{ln}{ \mathcal L } & = & -\displaystyle \frac{1}{2}{\chi }^{2}-\displaystyle \frac{1}{2}\mathrm{ln}\left|{C}_{\mathrm{ex}}\right|+\mathrm{constant},\end{array}\end{eqnarray} \tag{ 22 }$$

where the inverse covariance ${C}_{\mathrm{ex}}^{-1}$ is given by Equation (9).

We use the fit from individual fields for a consistency check. To provide a best estimate using the combination of all the fields that were observed at once, we also fit to the five CIBER fields using the joint likelihood:

$$\begin{eqnarray}&&\mathrm{ln}{ \mathcal L }=\sum _{i=1}^{{N}_{\mathrm{field}}}\mathrm{ln}{{ \mathcal L }}_{i}\end{eqnarray} \tag{ 23 }$$

where N_field = 5. Note that the PSF model is different for each field, so the information from different fields is combined in the likelihood.

We use the affine-invariant MCMC sampler emcee (Foreman-Mackey et al. 2013) to sample from the posterior distribution. We set flat priors for R_e,2, A_1h, and A_2h in the range of [10⁻⁴ R₂₀₀, R₂₀₀], [0, 50], and [0, 200], respectively. We use an ensemble of 100 walkers taking 1000 steps with 150 burn-in steps. We checked that the chains show good convergence by computing the Gelman-Rubin statistic R (Gelman & Rubin 1992). For all three parameters in all cases, we find R < 1.1.

Given the best-fitting extended galaxy profile, we can calculate the fraction of flux missed in photometric galaxy surveys using a limited aperture. From our model, the fraction of flux within a photometric aperture can be approximated by f_core ≡ L_core/(L_core + L_ext), where L_core and L_ext are the total flux in the core and extension profile (Equation (19)), respectively. In practice, there are various ways to perform photometry. The Petrosian flux (Petrosian 1976) is derived from aperture photometry and thus it is the most straightforward method to compare to our results. The Petrosian flux is defined by the total flux within a multiplicative factor of the Petrosian radius of sources. We obtain the Petrosian radius and Petrosian flux from the SDSS catalog of each stacked galaxy in our sample. In SDSS, the Petrosian flux is calculated by integrating the emission within twice the Petrosian radius. ¹⁹ With our galaxy profile, we can calculate the fraction of flux within the same radius (f_petro). The results are summarized in Table 4.

Table 4. Fraction of Flux in Core Component Compared to Flux Captured in Petrosian and SDSS Model Flux, Assuming the Galaxy Light Profile Follows the Stacking Results in This Work. The Total Row Shows the Weighted Average of the Five Listed Sub-samples

	1.1 μm	1.1 μm	1.1 μm	1.8 μm	1.8 μm	1.8 μm
Name	${f}_{\mathrm{core}}$	fpetro	fmodel	${f}_{\mathrm{core}}$	fpetro	fmodel
high-M/low-z	${0.79}_{-0.02}^{+0.04}$	${0.78}_{-0.10}^{+0.08}$	${0.84}_{-0.12}^{+0.11}$	${0.81}_{-0.02}^{+0.02}$	${0.80}_{-0.10}^{+0.07}$	${0.85}_{-0.12}^{+0.10}$
high-M/med-z	${0.81}_{-0.05}^{+0.04}$	${0.74}_{-0.13}^{+0.07}$	${0.78}_{-0.15}^{+0.08}$	${0.83}_{-0.03}^{+0.04}$	${0.75}_{-0.11}^{+0.08}$	${0.78}_{-0.13}^{+0.10}$
high-M/high-z	${0.86}_{-0.04}^{+0.06}$	${0.73}_{-0.16}^{+0.07}$	${0.77}_{-0.19}^{+0.15}$	${0.89}_{-0.04}^{+0.03}$	${0.75}_{-0.16}^{+0.07}$	${0.79}_{-0.18}^{+0.16}$
low-M/low-z	${0.84}_{-0.02}^{+0.04}$	${0.78}_{-0.11}^{+0.05}$	${0.80}_{-0.12}^{+0.10}$	${0.85}_{-0.03}^{+0.02}$	${0.79}_{-0.11}^{+0.05}$	${0.81}_{-0.12}^{+0.10}$
low-M/med-z	${0.89}_{-0.04}^{+0.04}$	${0.78}_{-0.16}^{+0.06}$	${0.80}_{-0.16}^{+0.09}$	${0.92}_{-0.03}^{+0.03}$	${0.80}_{-0.15}^{+0.06}$	${0.83}_{-0.14}^{+0.10}$
total	${0.83}_{-0.01}^{+0.02}$	${0.77}_{-0.06}^{+0.03}$	${0.80}_{-0.06}^{+0.05}$	${0.86}_{-0.01}^{+0.01}$	${0.78}_{-0.05}^{+0.03}$	${0.81}_{-0.06}^{+0.05}$

Download table as:  ASCIITypeset image

We also estimate the missing light fraction with the "model magnitude" given in SDSS (f_model). Rather than integrating within a certain aperture size, the model magnitude is derived by fitting the galaxy profile with an exponential or de Vaucouleurs functional form, choosing the one with the higher likelihood in the fitting. ²⁰ While it is difficult to apply the same fitting procedure to the sources in CIBER images, we can calculate the ratio between the model flux and the Petrosian flux of each source in the SDSS catalog and thus infer the fraction of missing light in the model flux. We find that both the Petrosian flux, which measures source emission within a limited aperture size, and the model flux derived from fitting a light profile to the small-radii regions of the galaxy, miss ∼20% of the total galaxy light, a deficit detected at ∼7σ (∼4σ) level for Petrosian (model) flux when combing constraints from all five sub-samples. This value is slightly larger than the light fraction in our galaxy extension term (∼10 to 20%). Our results on the missing light fraction in the Petrosian flux are in agreement with previous analytical calculation (Graham et al. 2005). Interestingly, Tal & van Dokkum (2011) probed the radial profile of z ∼ 0.34 luminous red galaxies (LRGs) in SDSS with a stacking analysis and they also found ∼20% of the total light missing at large radii when fitting a Sérsic model to individual galaxies. Although their galaxy samples are at somewhat higher mass (M_* ∼ 10¹¹ − 10¹² M_⊙) and model magnitudes are fitted with a different functional form, we arrive at a similar fraction of missing flux.

The Illustris simulation (Rodriguez-Gomez et al. 2016) traces the dynamics and merger history of stellar particles and estimates the "ex situ" population of stars that formed in other galaxies and were later stripped and accreted into a new galaxy. The shaded region in the left panel of Figure 13 shows the ex situ stellar mass fraction at z = 0 from the Illustris simulation (Rodriguez-Gomez et al. 2016). Although it is difficult to measure the ex situ component in observations, Huang et al. (2018) have studied individual stellar halos out to 100 kpc in more massive galaxies (10¹¹ M_⊙ ≲ M_* ≲ 10¹² M_⊙) at higher redshifts (z ∼ 0.4) in HSC images, finding that the fraction of stellar mass between 10 and 100 kpc is in good agreement with the ex situ fraction constraints from Illustris (Rodriguez-Gomez et al. 2016). In addition, Wang et al. (2019) probe the stellar halo around local (0 ≲ z ≲ 0.25) low-mass galaxies ($9.2\,{M}_{\odot }\lt \mathrm{log}{M}_{* }\lt 11.4\,{M}_{\odot }$) with a stacking analysis on HSC images in the r band. They stacked galaxies out to ∼120 kpc within several stellar mass bins. For each bin, they split the sources into low and high-concentration populations, defined by C < 2.6 and C > 2.6, where C = R₉₀/R₅₀ is the ratio of the radii that contain 90% and 50% of the r-band Petrosian flux.

Zoom In Zoom Out Reset image size

CIBER extends the HSC measurements to higher redshifts and longer wavelength bands. Armed with light profile fits, we can quantify the luminosity fraction in the extended stellar halo around the stacked sources. The left panel of Figure 13 shows the fraction of stellar flux between radii of 10 and 100 kpc, using the fitted galaxy profile from CIBER and HSC (Wang et al. 2019). We observe that ∼50% of the flux originates at galactocentric distances between 10 and 100 kpc. Wang et al. (2019) re-scaled their images to physical units before stacking, whereas in our analysis we stack sources in observed angular units. Therefore, the variations in our measurements are mostly due to the variation of the conversion factor from angular to physical units for each galaxy in our stack. Our constraints are consistent with the HSC results in the highest mass bin.

Both CIBER and HSC are consistent with the ex situ fraction from Illustris at z = 0, but are systematically higher than the median value from Illustris (the gray line in Figure 13). One possible explanation is that the flux between 10 and 100 kpc is not a perfect proxy of the ex situ population for lower mass galaxies. For example, D'Souza et al. (2014) has shown that the transition scale between in situ and ex situ components varies across a wide range from ∼10 to ∼50 kpc, depending on the stellar mass and concentration of the galaxies. Nevertheless, given the limited information in stacking, we use this definition to associate the luminosity from beyond 10 kpc with IHL.

We also show the fraction of flux with a 20 kpc radius cut in the right panel of Figure 13. A radius cut at 20 kpc is a more suitable choice to describe the stripped stellar populations. For example, Milky Way-sized simulations suggest that the infalling stellar debris is recaptured by the galaxy and results in disk thickening at r ≲ 10 kpc (Mazzarini et al. 2020). We note that each galaxy has a different stellar halo profile; interpreting our stacking results requires knowing the stellar halo profile on average over a large sample of galaxies. Given the uncertainty in choosing an average IHL radius, we report our results in both 10 kpc and 20 kpc scales, while showing the full radial range in Figure 14. We find that ∼25% of galaxy fluxes are from outside 20 kpc. The CIBER constraints shown in Figure 13 are summarized in the Appendix.

Zoom In Zoom Out Reset image size

With the galaxy profile from CIBER and HSC, we can estimate the EBL contribution from the extended regions at the redshift of our stacked sources. We model this quantity in the following steps:

• 1.  
  For any given radius cut r_cut, we model the fraction of light beyond r_cut as a function of stellar mass by fitting a line to all CIBER and HSC data points in logarithmic space;
• 2.  
  We estimate total stellar mass density by integrating the stellar mass function from Muzzin et al. (2013; we take their single Schechter function fit with all samples in 0.2≤z < 0.5 bin, approximately the redshift of our sources);
• 3.  
  For each r_cut, we apply the fraction derived in step 1 to the stellar mass function, and integrate to get the stellar mass density from sources outside r_cut;
• 4.  
  Assuming the mass-to-light ratio is the same for all galaxies, the ratio between step 3 and step 2 is our estimate of the EBL fraction from extended sources as a function of r_cut.

The results are shown in Figure 14. We get approximately 30/15% of extended emission in the EBL with r_cut = 10/20 kpc, respectively. Note that these values are close to the fraction in the five individual stellar mass bins from our stacking results. This is expected as our samples are at ∼L* scale, which are the representative population that contains the majority of the total stellar emission of their redshift.

The fraction of the total emission from a dark matter halo associated with IHL, f_IHL, has been investigated with both observation and theoretical modeling (e.g., Lin & Mohr 2004; Gonzalez et al. 2005; Purcell et al. 2007; D'Souza et al. 2014; Burke et al. 2015; Elias et al. 2018). With our stacking results, we can estimate the total halo emission from the sum of the galaxy light and one-halo terms. For the IHL, we consider the extended galaxy emission beyond r_cut = 10/20 kpc of all the bright (m_1.1 < 20) galaxies in the halo, noting that m_1.1 = 20 is also our choice of flux threshold for masking. Therefore, the IHL fraction f_IHL can be expressed as

$$\begin{eqnarray}&&{f}_{\mathrm{IHL}}=\displaystyle \frac{{\sum }_{{m}_{1.1}\lt 20}L(\gt {r}_{\mathrm{cut}})}{{\sum }_{{m}_{1.1}\lt 20}L+{\sum }_{\mathrm{faint}}L},\end{eqnarray} \tag{ 24 }$$

where ${\sum }_{{m}_{1.1}\lt 20}L$ is the total light associated with bright galaxies, and ${\sum }_{{m}_{1.1}\lt 20}L(\gt {r}_{\mathrm{cut}}$) is the part of bright galaxy emission beyond r_cut. ∑_faint L represents the light from faint galaxies as well as the unbound stars in the halo, captured in the one-halo luminosity. Note that we conservatively assume the one-halo luminosity arises entirely from faint, gravitationally bound galaxies. However, it is certainly true that some one-halo light arises from unbound stars as is readily observed in images of massive clusters at low redshift.

From our stacking profile, the faint source emission ∑_faint L can be described by the total emission in the one-halo term, L_1h . ²¹ For the bright sources, we define

$$\begin{eqnarray}&&\sum _{{m}_{1.1}\lt 20}L={L}_{\mathrm{gal}}\cdot {N}_{\mathrm{eff}},\end{eqnarray} \tag{ 25 }$$

where L_gal is the total light in the galaxy profile term from our stacking results, which describes the averaged light of the galaxies within each stacking sample. N_eff accounts for the fact that there are multiple bright galaxies in the halo, and we infer the average N_eff value from MICECAT. For our five stacking sub-samples, we get N_eff ∼ 2–5. From our fitted galaxy profile, we can also calculate L_gal (> r_cut), and we apply the same N_eff to model the extension from other bright galaxies:

$$\begin{eqnarray}&&\sum _{{m}_{1.1}\lt 20}L(\gt {r}_{\mathrm{cut}})={L}_{\mathrm{gal}}(\gt {r}_{\mathrm{cut}})\cdot {N}_{\mathrm{eff}}.\end{eqnarray} \tag{ 26 }$$

This results in

$$\begin{eqnarray}&&{f}_{\mathrm{IHL}}=\displaystyle \frac{{L}_{\mathrm{gal}}(\gt {r}_{\mathrm{cut}})/{L}_{\mathrm{gal}}}{1+{L}_{1{\rm{h}}}/\left({N}_{\mathrm{eff}}\cdot {L}_{\mathrm{gal}}\right)}.\end{eqnarray} \tag{ 27 }$$

We show our constraints on f_IHL as a function of halo mass and redshift in Figures 15 and 16, respectively. The halo masses associated with our galaxies are inferred from the MICECAT simulation and using the SDSS photometric redshifts. The CIBER data points shown in Figures 15 and 16 are summarized in the Appendix.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Note that the fraction of light beyond r_cut (the numerator in Equation (27)) is shown in Figure 13, where the higher redshift bins have slightly higher values. However, in Figure 16, they have lower f_IHL. This is due to the increase of the one-halo term with redshift. We show the ratio of the one-halo term and the stacked galaxy light in Figure 17. Note that this observable quantity tracks the evolution of the one-halo luminosity but lacks the N_eff term in Equation (27) derived from simulations. We compare with the same quantity from the MICECAT simulation, where the one-halo term includes all the unmasked faint galaxies and residual bright source emission outside the mask due to the PSF. We detect a strong redshift evolution of one-halo contribution compared with the MICECAT simulation, which could be attributed to the unbound stars that are not included in MICECAT.

Zoom In Zoom Out Reset image size

We compare our results with f_IHL from previous work, including the Milky Way (Carollo et al. 2010), the Andromeda galaxy (M31; Courteau et al. 2011), the ICL fraction in individual galaxy groups and clusters (Gonzalez et al. 2005, 2007; Burke et al. 2015), and an analytical model (Purcell et al. 2007, 2008). Our results follow a more gradual redshift evolution trend than reported in massive clusters (Burke et al. 2015; see Figure 16).

We calculate the m_1.1−m_1.8 color of the inner and outer region of the galaxy, defined by the total light inside and outside 20 kpc physical scale in the fitted galaxy profile. The results are summarized in Table 5. Note the definition of inner and outer component here is based on the intrinsic profile, which is different from the core/extension separation using the stacked PSF defined in Equation (19). We have no detection of a color difference between the inner and outer regions in the two CIBER bands. We also find similar inner and outer region color with 10 kpc radius cut. Previous measurements in optical bands found that the galaxy outskirts are bluer than their core (e.g., D'Souza et al. 2014; Huang et al. 2018). For comparison, we calculate the m_1.1−m_1.8 color of galaxy cores in MICECAT sources selected from the same criteria, as well as from the empirical galaxy model of Helgason et al. (2012) at z = 0.3, approximately the redshift of our samples. Our inner region color is consistent with these models. To model the extension, we use a collection of elliptical galaxy spectra from the population synthesis package GISSEL (Bruzual & Charlot 1993) redshifted to z = 0.3. We also estimate the extension color using an imaging study on the local spiral galaxy NGC 5907 (Rudy et al. 1997). We use their ratio of I band and J band flux in >1 arcmin regions to approximate the m_1.1−m_1.8 extension color. The rest-frame I and J band redshifted to z ∼ 0.3 (approximately the redshift of our samples) are close to the two CIBER bands. NGC 5907 shows a redder spectrum than our galaxy extension, whereas the elliptical galaxy spectrum template is slightly bluer than our samples. In addition, the IHL constraints from Zemcov et al. (2014) are also given in Table 5, but we note that Zemcov et al. (2014) reflects the integrated IHL from all redshifts.

Table 5. Constraints on the Color (m_1.1 − m_1.8) of the Galaxy Inner and Outer Components

Name	Inner	Outer
high-M/low-z	${0.42}_{-0.17}^{+0.20}$	${0.36}_{-0.31}^{+0.34}$
high-M/med-z	${0.54}_{-0.27}^{+0.25}$	${0.46}_{-0.25}^{+0.24}$
high-M/high-z	${0.65}_{-0.28}^{+0.31}$	${0.61}_{-0.25}^{+0.31}$
low-M/low-z	${0.39}_{-0.18}^{+0.20}$	${0.37}_{-0.37}^{+0.41}$
low-M/med-z	${0.56}_{-0.24}^{+0.23}$	${0.44}_{-0.44}^{+0.50}$
total	${0.49}_{-0.10}^{+0.10}$	${0.47}_{-0.15}^{+0.14}$
MICECAT	0.44 ± 0.07
Helgason et al. (2012)	0.41
GISSEL		0.32 ± 0.08
NGC 5907		1.41 ± 0.61
Zemcov et al. (2014)		${0.89}_{-1.08}^{+1.17}$

Note. The +/− values indicate 68% interval ranges. The total row shows the weighted average of five sub-samples. For comparison, we also show models of core color from MICECAT and an analytical prescription from Helgason et al. (2012) at z = 0.3. For the extension, we compare our results with spectra from a population synthesis code, GISSEL (Bruzual & Charlot 1993), and the outskirts of NGC 5907 redshifted to z = 0.3 (Rudy et al. 1997). The color of EBL fluctuations attributed to redshift-integrated IHL from Zemcov et al. (2014) is also shown.

Download table as:  ASCIITypeset image

The one-halo amplitude is detected in the 1.8 μm band at the ∼4σ level in the "total" and "high-M/high-z" cases, and at the ∼3σ level in "mag bin #4" and "high-M/med-z" cases. One-halo clustering is not clearly detected at the 1.1 μm band since the photocurrent from sources is lower in this band. The one-halo amplitude A_1h is consistent with unity to within ∼2σ, which implies that our one-halo templates built from MICECAT are sufficient to describe the clustering within halos of our stacked samples. However, from our stacking results, it is unclear if this emission actually consists of discrete galaxies as given in the MICECAT simulation. Two-halo clustering is not detected in all cases since the large-scale clustering signal is comparable to the current uncertainties in the measurement.