Probing Intra-Halo Light with Galaxy Stacking in CIBER Images

We study the stellar halos of $0.2\lesssim z \lesssim 0.5$ galaxies with stellar masses spanning $M_*\sim 10^{10.5}$ to $10^{12}M_\odot$ (approximately $L_*$ galaxies at this redshift) using imaging data from the Cosmic Infrared Background Experiment (CIBER). A previous CIBER fluctuation analysis suggested that intra-halo light (IHL) contributes a significant portion of the near-infrared extragalactic background light (EBL), the integrated emission from all sources throughout cosmic history. In this work, we carry out a stacking analysis with a sample of $\sim$30,000 Sloan Digital Sky Survey (SDSS) photometric galaxies from CIBER images in two near-infrared bands (1.1 and 1.8 $\mu$m) to directly probe the IHL associated with these galaxies. We stack galaxies in five sub-samples split by brightness, and detect an extended galaxy profile, beyond the instrument point spread function (PSF), derived by stacking stars. We jointly fit a model for the inherent galaxy light profile, plus large-scale one- and two-halo clustering to measure the extended galaxy IHL. We detect non-linear one-halo clustering in the 1.8 $\mu$m band, at a level consistent with numerical simulations. Our results on the galaxy profile suggest that $\sim 50\%$ of the total galaxy light budget in our galaxy sample resides in the outskirts of the galaxies at $r>10$ kpc. We describe this extended emission as IHL and and are able to study how this fraction evolves with cosmic time. These results are new in the near-infrared wavelength at the $L_*$ mass scale, and suggest that IHL has a significant contribution to the integrated galactic light, and to the amplitude of large-scale background fluctuations.


INTRODUCTION
In the standard cosmological paradigm, galaxies grow hierarchically through merger and accretion. Galaxies accreting onto more massive systems become disrupted and stars stripped away from their parent galaxies become redistributed in the merged dark matter halo. This results in extended stellar halos that are known to span ycheng3@caltech.edu tens or hundreds of kilo-parsecs. The stellar emission from this material is sometimes referred to as "intrahalo light" (IHL), or in massive galaxy clusters as "intracluster light" (ICL).
An independent way to study the aggregate emission from diffuse sources like IHL is through measurements of the extragalactic background light (EBL), which encodes the integrated emission from all sources across cosmic history (Cooray 2016). Absolute optical and near-infrared EBL photometry has proven challenging as measurements must tightly control systematic errors and carefully model and subtract local foregrounds (e.g., Kawara et al. 2017;Zemcov et al. 2017;Matsuura et al. 2017;Matsumoto et al. 2018;Lauer et al. 2021). Several authors (Bernstein 2007;Levenson et al. 2007;Tsumura et al. 2013a;Matsumoto et al. 2015;Sano et al. 2015;Zemcov et al. 2017;Matsuura et al. 2017;Sano et al. 2020;Lauer et al. 2021) have reported potential detections above the integrated galaxy light (IGL) derived from galaxy counts (Keenan et al. 2010;Domínguez et al. 2011;Helgason et al. 2012;Driver et al. 2016;Saldana-Lopez et al. 2020;Koushan et al. 2021), which may indicate the existence of extragalactic emission missed in source-counting surveys.
Additionally, EBL fluctuation analyses have also consistently reported excess fluctuations over those expected from the IGL (Kashlinsky et al. 2005; Thompson et al. 2007;Matsumoto et al. 2011;Kashlinsky et al. 2012;Cooray et al. 2012;Zemcov et al. 2014;Mitchell-Wynne et al. 2015;Seo et al. 2015;Kim et al. 2019;Matsumoto & Tsumura 2019). One explanation is emission from the epoch of reionization (Kashlinsky et al. 2005;Matsumoto et al. 2011;Kashlinsky et al. 2012;Mitchell-Wynne et al. 2015), while other studies suggest IHL contributes most of the excess fluctuations (Cooray et al. 2012). In particular, Zemcov et al. (2014) interpret imaging data from the Cosmic Infrared Background Experiment (CIBER) as arising from an IHL intensity comparable to the IGL at near-infrared wavelengths. This result would imply that stars diffusely scattered in dark matter halos may account for a non-negligible fraction of the near-IR cosmic radiation budget. The absorption spectra from blazars constrain the EBL column density along the line of sight (e.g., Aharonian et al. 2006Aharonian et al. , 2007MAGIC Collaboration et al. 2008;Abdo et al. 2010;Ackermann et al. 2012;H. E. S. S. Collaboration et al. 2017;Ackermann et al. 2018;Abeysekara et al. 2019;Acciari et al. 2019;Abdalla et al. 2020). While IHL is generally produced at low redshifts, improving the uncertainties in its redshift history helps place IHL in the context of these constraints.
In this work, we further constrain the IHL using CIBER broadband imaging. Rather than studying EBL intensity fluctuations as in Zemcov et al. (2014), we perform a stacking analysis to directly probe the stellar halos around galaxies. We stack a sample of ∼ 30, 000 Sloan Digital Sky Survey (SDSS) photometric galaxies at z ∼ 0.2-0.5 across five 2 × 2 deg 2 fields. Our samples span a range of stellar masses at approximately L * scales at this redshift (Muzzin et al. 2013). Although we only study stellar halos around a subset of galaxies rather than the aggregate population as probed by fluctuations, stacking provides a direct path to probe the IHL associated with this sample. Stacking complements fluctuation measurements by probing the relationship between individual galaxies and their stellar halos. Stacking also allows us to investigate how stellar halos depend on host-galaxy properties, e.g, stellar mass, redshift, etc. A complementary fluctuation analysis of these same data is currently in progress.
This paper is organized as follows. First, we introduce CIBER in Sec. 2 and the data processing in Sec. 3. Sec. 4 and Sec. 5 describe the external data sets used in this work, including observed and simulated source catalogs. Sec. 6 details the stacking procedure and Sec. 7 describes the point spread function (PSF) model. The stacking results are presented in Sec. 8. Sec. 9 introduces the theoretical model we use to fit the data and the parameter fitting procedure. The results on model parameter constraints are given in Sec. 10, and further discussion is presented in Sec. 11. Sec. 12 summarizes the paper. Throughout this work, we assume a flat ΛCDM cosmology with n s = 0.97, σ 8 = 0.82, Ω m = 0.26, Ω b = 0.049, Ω Λ = 0.69, and h = 0.68, consistent with the measurement from Planck (Planck Collaboration et al. 2016). All fluxes are quoted in the AB magnitude system.

CIBER EXPERIMENT
CIBER 1 (Zemcov et al. 2013) is a rocket-borne instrument designed to characterize the near-infrared EBL. CIBER consists of four instruments: two wide-field imagers (Bock et al. 2013), a narrow-band spectrometer (Korngut et al. 2013), and a low-resolution spectrometer (Tsumura et al. 2013b). CIBER has flown four times in February, 2010July, 2012March, and 2013 June. The first three CIBER flights were launched at White Sands Missile Range, New Mexico on a Terrier- Note-We discard the beginning part of the Elat30 field integration due to pointing instability.
Black Brant IX rocket. These flights reached ∼ 330 km apogee with ∼ 240 s of exposure time, and the payload was recovered for future flights. The fourth flight was a non-recovery flight launched 3:05 UTC 2013 June 6 from Wallops Flight Facility, Virginia on a four-stage Black Brant XII rocket. The payload reached 550 km altitude, much higher than the two-stage rocket used in the previous three flights. This gives more exposure time (335 s) for observing more science fields with long integrations to achieve better sensitivity and systematics control. This work presents the first science results from the CIBER fourth flight imager data. The data from previous flights have been studied with a fluctuation analysis, published in Zemcov et al. (2014). With a large field of view and low sky background above the atmosphere, CIBER imaging provides fidelity on angular scales from 7 to 2 • . For stacking, CIBER imaging can trace lowsurface-brightness emission on degree angular scales providing a unique dataset compared with ground-based or small field-of-view space-borne studies. Each CIBER imager uses a 1024 × 1024 pixel HAWAII-1 HgCdTe detector. The two imagers are identical except for their λ/∆λ ∼ 2 filters, which are centered at 1.05 and 1.79 µm 2 .
During its fourth flight, CIBER observed 8 science fields with ∼ 50 s integrations sampled at 1.78 s intervals. We discard the first three fields in this analysis due to contamination from airglow that produces a strong non-uniform emission across the images that requires aggressive filtering which also significantly reduces our signal (Zemcov et al. 2014). Table 1 summarizes the sky coordinates and the integration time of the five science fields used in this work. In the beginning of the Elat30 integration, the rocket's pointing was not stable which has the effect of smearing the PSF on the sky. As a re-2 In the first and second CIBER flights, the longer wavelength band is centered at 1.56 µm, and thus it is named the 1.6 µm band in previous CIBER publications (Bock et al. 2013;Zemcov et al. 2014).
sult, we only use the last 16 s of this integration in our analysis.

DATA PROCESSING
In this section, we describe the data reduction from the raw flight data to the final images used for stacking.

Raw Time Stream to 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 −1 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 Fig. 1 and Fig. 2 show the raw images of the SWIRE field in the CIBER 1.1 and 1.8 µm bands, respectively.

Dark Current
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 −1 , less than 10 % of the sky brightness. Panel B of Fig. 1 and Fig. 2 show the dark-current maps of CIBER 1.1 and 1.8 µm bands, respectively.  Figure 1. Images from the SWIRE field in the 1.1 µm band. A: the raw image of the photoccurent map. B: dark-current template constructed from dark images before the flight. C: instrument mask encoding the pixels with fabrication defects, unusual photocurrents, and cosmic-ray contamination. D: source mask for bright stars and galaxies in the 2MASS and Pan-STARRS catalogs. E: flat-field estimator from averaging the other four sky fields. F: raw image after dark-current subtraction, flat-field correction, and calibration. G: image in Panel F after (constant) background removal and masking. This image is smoothed with a σ = 35 Gaussian kernel to highlight large-scale fluctuations. H: image in Panel G after subtracting a fitted 2D polynomial, also shown smoothed with a σ = 35 Gaussian kernel. Compared to Panel G, we see that the large-scale background fluctuations have been reduced after filtering. This is the final product of the data reduction pipeline. We mask pixels that meet at least one of the following conditions: (1) a fabrication defect, (2) poor timestream 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).

Instrument Mask
Pixels with fabrication defects and significant multiplexer glow are mostly distributed near the edges or corners of each quadrant on the detector arrays. They exhibit pathologies in their photocurrent response and can be found by comparison to the population of normal pixels. We perform a 3σ clipping on stacked dark images (the same dataset used for a dark-current template in Sec. 3.2) to identify these pixels.
During integration, some cosmic-ray events or electronic transients leave a step feature in the time stream. We use a 100σ clip on each time stream to pick out pixels that show these abrupt changes during an integration. Sometimes cosmic-ray events also leave a comet-like structure on the array and these regions are also masked. The union of the pathological pixel, timestream masks, and cosmic-ray masks form the instrument mask. In total, ∼ 10% of pixels are removed by the instrument mask. Panel C of Fig. 1 and Fig. 2 show the instrument masks in the SWIRE field of the 1.1 and 1.8 µm bands, respectively.

Source Mask
To remove bright foreground stars and galaxies in our fields, we use position and brightness information from the Pan-STARRS and 2MASS catalogs (see Sec. 4 for details). We further derive source magnitudes in the two CIBER bands, m 1.1 and m 1.8 , from these catalogs, as detailed in Sec. 4. We mask all point sources brighter than m 1.1 = 20, choosing a masking radius for each source derived as follows. With the modeled instrument PSF (Sec. 7.3), the masking radius is chosen such that for each source, pixels with intensity brighter than νI th ν = 1 nW m −2 sr −1 in the 1.1 µm band are masked. This choice of threshold value removes ∼ 50% of pixels in each field. We apply the same masking radius to 1.8 µm band sources. The same masking function is also applied to simulations to account for residual emission from bright sources outside the masks and the unmasked faint populations. Panel D of Fig. 1 and Fig. 2 shows the SWIRE field source mask in the CIBER 1.1 and 1.8 µm bands, respectively.
The final mask we apply to the data is the union of the instrument mask and source mask. After applying these masks, we apply a final 3σ pixel clipping mask to identify additional outliers not flagged through the other methods (e.g., from low-energy cosmic-ray events or electronic tranisents).

Flat-Fielding
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-currentsubtracted 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 Fig. 1 and Fig. 2, respectively.

Surface Brightness Calibration
Throughout this work, we use nW m −2 sr −1 for the units of surface brightness (νI ν ). The calibration factor, C, which converts photocurrent (e − s −1 ) to intensity (nW m −2 sr −1 ) 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 3 ; 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 Fig. 1 and Fig. 2 shows the SWIRE field images masked by instrument masks at 1.1 and 1.8 µm, respectively, after flat-fielding and calibration.

Background Removal
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 Fig. 1 and Fig. 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.

Image Filtering
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 Fig. 1 and Fig. 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 a much larger scales than these signals and the signal filtering is accounted for in simulations (see Sec. 9).
To match the catalog sources to our data, we fit the astrometry coordinates of our images with the online software nova.astrometry.net (Lang et al. 2010). For each image, we solve for the astrometry in four quadrants separately to mitigate the effect of image distortion. Since there is a fixed ∼ 50 misalignment between the 1.1 and 1.8 µm images as they are produced by different telescopes, their astrometry is solved separately.

Pan-STARRS
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 SEDfitting 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.

2MASS
Some bright stars are not included in the Pan-STARRS catalog, and 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 Sec. 7).

SDSS
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, 3galaxies) 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 (Sec. 8.1), as well as selecting stars for stacking to model the PSF (Sec. 7). -Robinson et al. (2008-Robinson et al. ( , 2013 performed SED fitting on ∼ 10 6 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 Sec. 8.2).

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

Nearby Cluster Catalog
Nearby galaxy clusters along the line of sight introduce extended emission in stacking, so we exclude galaxies that are close to nearby clusters (Sec. 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.

SIMULATION CATALOG-MICECAT
In addition to the observed source catalogs, we make use of the MICECAT simulated galaxy catalog (Fosalba et al. 2015a,b;Hoffmann et al. 2015) to estimate the signal from galaxy clustering. MICECAT is a product of the N -body cosmological simulation MICE Grand Challenge run (MICE-GC), which has 70 billion dark matter particles in a 3072 3 Mpc 3 h −3 cubic co-moving box. The dark matter halos are resolved down to ∼ 3 × 10 10 M h −1 .
MICECAT is a mock catalog that simulates ideal observations of a 5000 deg 2 light cone covering 0 < z < 1.4. MICECAT builds on MICE-GC by combining a halo occupation distribution (HOD) with subhalo abundance matching (SHAM) to calibrate to observed luminosity functions and clustering . MICECAT simulates a mass-limited sample complete to m i ∼ 22 and m i ∼ 24 at z 0.5 and z 0.9, respectively . The MICECAT mocks are large enough to permit us to generate up to ∼ 10 3 independent CIBER field-sized (2 × 2 deg 2 ) mock catalogs. We use modeled magnitudes from MICECAT in Euclid NISP Y and H bands for CIBER m 1.1 and m 1.8 , respectively, since the NISP filters are similar to the CIBER imager bands.
MICECAT simulates both central and satellite galaxies generated with its HOD+SHAM model, which allows us to model the linear (two-halo) and nonlinear (onehalo) clustering in the stacking signal separately. We use the radial shapes derived from MICECAT stacking to fit the one-halo and two-halo amplitudes in our stacking data. Details on modeling galaxy clustering in the stacking signals are further described in Sec. 9. 6. STACKING 6.1. Sub-pixel Stacking 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 (P SF instr ) 4 , and the pixel function P SF pix : 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 P SF stack (r) ≡ P SF instr (r) P SF pix (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 subpixel 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): As a practical matter, P SF pix can be determined through simulations. P SF stack (r) can be measured by stacking stars in the field, where Σ src (r) is a delta function, so Σ stack (r) = P SF stack (r). Note that the expression in the second line of Eq. 1 implies that the intrinsic profile Σ src (r) can be obtained from the stacked profile Σ stack (r) with the knowledge of P SF stack (r), instead of determining P SF instr (r). We perform stacking and PSF modeling separately for each field, since P SF 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 Fig. 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.

Covariance Matrix of Stacking Profile
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, and 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, Σ k stack , 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 where Σ stack is the average stacked profile of all of the sub-regions. One of our galaxy stacking samples (mag bin # 1 in Sec. 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, Σ k stack , 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 (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.

PSF MODELING
An accurate model for the PSF is essential for quantifying the galaxy extension from stacking images. As stars are point sources on the sky, we measure the PSF of each field by stacking stars in the same CIBER field. The radial profile of star stacks gives P SF stack (Eq. 1), which accounts for all effects that distribute the light from a point source to the stacked profile, including spreading by the instrument optical system and detectors, pointing instability during integration, astrometry uncertainties, and the pixel function P SF pix . Since we use bright stars in the CIBER fields to model the PSF, the uncertainty on the PSF is subdominant to our galaxy stacked profiles.

Modeling P SF stack
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 P SF 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 < m 2MASS J < 10 sources; for outer radii, we fit another slope to the stacking profile of the brightest 4 < m 2MASS J < 9 sources and connect the two slopes at r = 40 (m 2MASS J 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 < m 2MASS J < 9, 9 < m 2MASS J < 10, 13 < m 1.1 < 14) are statistically consistent. The top panel of Fig. 3 shows P SF stack from the SWIRE field in the 1.1 µm band. The top panel of Fig. 5 shows P SF 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.

Validating P SF stack
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 P SF 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 −4 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 0.7 (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 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 P SF 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 (Sec. 8.4). The results in the 1.1 µm band SWIRE field are shown on the bottom panel of Fig. 3. In Fig. 4 we show the difference of Gaia star stacks and the P SF stack model. The errors are propagated from the covariance of the P SF stack model and Gaia star stacks. We also show the χ 2 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.

Modeling P SF instr
Although knowledge of the instrument PSF is not required for reconstructing the source profile Σ src from the stacking profile Σ stack , P SF instr is still needed when we model the clustering signal from a simulated catalog (Sec. 9), where we make mock galaxy images using the CIBER PSF and pixel gridding. P SF instr is also useful for determining the masking radius for bright sources (Sec. 3.3.2).   We illustrate the process of constructing and validating the P SF stack (r) model, in the 1.1 µm band SWIRE field. Top: star-stacking profile in three different brightness bins (blue, orange, and green), and the combined P SF stack (r) model (black dashed curve) derived from splicing these three stacking profiles together at the radii marked by the black vertical dashed lines. The black data points show the binned P SF stack (r) and the error bars propagated from their original star stacks. The filled data points and the three colored solid curves are the data used in the P SF stack (r) model. Bottom: comparison of the P SF stack (r) model with the stacking profiles from fainter stars selected from Gaia. The four chosen brightness bins match the ones used in galaxy stacking. The P SF stack (r) model agrees closely with the star-stacking profiles as shown in Fig. 4. P SF instr is modeled as follows: first, we deconvolve P SF pix (r) (Eq. 2) from the P SF 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 P SF instr , we fit a β model (Cavaliere & Fusco-Femiano 1978)  the deconvolved image: 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 P SF stack . The bottom panel of Fig. 5 illustrates this procedure in the 1.1 µm band of the SWIRE field. The P SF stack model, obtained from star stacks in three different brightness bins, matches the β model of P SF instr convolved with the pixel function P SF pix (Eq. 2). Our instrument PSF has a size comparable to a pixel (FWHM ∼ 7 ).

GALAXY STACKING
We stack galaxies within magnitude ranges 16 < m 1.1 < 20, divided into several sub-samples spanning ∆m 1.1 = 1. Our choice of magnitude bins optimizes the The variation across fields is due to the difference in pointing stability. Bottom: demonstration of the P SFinstr reconstruction process. Black data points show the P SF stack model in the 1.1 µm band SWIRE field, derived from splicing the starstacking profile in three different brightness bins (c.f. Fig. 3 top panel). The blue line is the P SFinstr model derived from fitting a β model to P SF stack after deconvolving P SFpix with the Richardson-Lucy deconvolution algorithm. The orange line shows the convolution of P SFinstr with P SFpix matching the P SF stack model, as a consistency check. Our model for P SFinstr is in agreement with data for r 30 . Our analysis is not susceptible to the moderate error at larger radii, as P SFinstr is only used for characterizing the clustering signal from nearby galaxies.
SNR on the stacks, giving sufficient sample sizes for each source brightness.

Source Selection Criteria
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: • 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; • Sources are located in the instrument mask; • Other Pan-STARRS sources exist in the same CIBER pixel; • 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; • 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; • Sources are within (1) a 500" radius of any galaxy cluster in Abell (1958) (Sec. 4.6); or (2) R 200 of any galaxy cluster with halo mass M h > 10 14 M or redshift z < 0.15 in the SDSS cluster catalog (Wen et al. 2012, Sec. 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 Sec. 9 for a detailed description of stacking with MICECATgenerated 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 14 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 Sec. 8.2), we can reduce the one-halo power by ∼ 3 − 5× at 100 radius just by excluding galaxies near clusters following our criteria.
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.5log 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/highz," "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 12 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 200 (in physical and angular units), in Table 2. The virial radius is calculated from 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 subsamples is summarized in Fig. 6 and Table 2.

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

Excess Profile
We define an "excess profile" Σ ex (r) as follows: where the normalization factor A is chosen such that P SF stack matches Σ stack at the innermost radial bin r 1 , and thus by construction, Σ ex (r 1 ) = 0, and A ≡ Σ stack (r 1 )/P SF stack (r 1 ).
Since the excess profile is fixed at r 1 , the uncertainties on the galaxy profile and the PSF profile at r 1 have to be accounted for by propagating this error to the other radial bins, and thus the excess profile covariance is given by (7) where C PSF and C stack are the covariance of P SF stack and Σ stack , respectively, and 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 5 . After this down-sampling, we derive the inverse covariance estimator by where N J = 64, the number of sub-groups used for estimating covariance, and the number of bins N bin = 15. C * −1 ex is the direct inverse of the C ex matrix, and the 5 Mag bin # 1 is down-sampled to 7 radial bins as its degree of freedom is limited by the small number of stacked sources. pre-factor in Eq. 9 de-biases the inverse covariance estimator, as our covariance matrix is derived from our data (Hartlap et al. 2007) 6 . 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 (Eq. 7), thus our model parameter inference (Sec. 9) does not depend on the definition of the excess profile.
We present field-averaged excess profiles in Fig. 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.

MODELING THE GALAXY PROFILES
We model the galaxy profile with three components as follows. We start by decomposing the stacked profile in image space (Sec. 9.1), defining fitted profiles (Sec. 9.2), and introducing our model for each component of the stack (Sec. 9.3). Finally, we describe the model fitting procedure in Sec. 9.4. 6 For mag bin # 1, N bin = 7, and N J = 64 is replaced by N B = 1000 since we use bootstrap resampling method in this case.

Components in Image Space
The raw CIBER image, I raw , can be expressed as 7 where x is the 2D pixel coordinate, F F 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 retrive I raw : (11) where I n (x) = I n (x)/F F (x), the instrument noise divided by the flat-field response. For simplicity, we ignore the error in the flat-field estimator in Eq. 11. In practice, the flat-field estimation uncertainties will not bias

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 (fcen) is derived by applying the same cuts to a simulated catalog from MICECAT.
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 Sec. 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 F [I raw (x), M (x)], which is a function of the input map I raw (x) and mask M (x). As described in Sec. 3.7, we choose F to be a third (1.1 µm)/fifth (1.8 µm)-order 2D polynomial function fitted to the masked I raw map 8 . The image used for stacking I map can thus be written as where 8 Note that the filter map F can be decomposed into the sum of three filter maps because the polynomial fitting is a linear operation, i.e., given two maps A(x) and B(x), and a mask M (x),

Components in the Stack
The stacked profile Σ stack can be expressed as the sum of stacks on the three maps in Eq. 12: Σ stack (r) = Σ sig stack (r) + Σ LoS stack (r) + Σ n stack (r).
The last two terms can be ignored in modeling since they are uncorrelated with the stacked sources, so Σ LoS stack (r) = Σ n stack (r) = 0. We model the stacked galaxy profile as where the first three terms are the signal terms and the last term is the filtered signal map in Eq. 13. The galaxy profile term, Σ gal stack , represents the intrinsic galaxy profile, which includes the galaxy shape and the extended stellar halo. We decompose the galaxy profile term, Σ gal stack , into "core" and "extended" parts: Σ gal stack (r) = Σ gal core (r) + Σ gal ext (r), Top: galaxy stacked profile Σ stack (black) and P SF stack model (orange dashed), scaled to match the innermost radial bin of Σ stack . The error bars give the diagonal element of the covariance matrix derived by the jackknife method (described in Sec. 3). Middle: the excess profile (Σex, Eq. 6) for the case shown in the top row. The excess is defined as the difference between the galaxy stacked profile and the P SF stack model, i.e., the difference of the black data from the orange curve in the top row. Bottom: the field-averaged excess profile Σex for mag bin #2, derived from the weighted average of the excess profile in the five individual fields. The improved sensitivity from combining fields can be seen compared to the middle row. The purple and brown dashed lines mark the pixel size and the median R200 values inferred from MICECAT, respectively. Open circles in all the plots represent negative values.
where the core component is the integrated emission of the P SF stack fitted to the stacking profile, i.e., the A · P SF stack term in Eq. 6, and the extended component is the rest of the galaxy emission: 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 (twohalo) dark matter halos. Σ 1h stack and Σ 2h stack 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 Fig. 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).

Modeling the Stacked Galaxy Profile
The stacked galaxy profile Σ gal stack (r) = Σ gal (r) P SF stack (r), is the intrinsic galaxy profile Σ gal , including the galaxy shape and the extended stellar halo, convolved with P SF stack . Following Wang et al. (2019), we model Σ gal with a double Sersic function: 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   Figure 8. Stacked profile (black data) of each sub-sample stack averaged over 5 CIBER fields in the 1.1 µm (top) and 1.8 µm (bottom) bands. Red lines and shaded regions indicate the median and 68% confidence interval of the joint fit constrained through MCMC, respectively. The blue, green, and orange solid lines show the best-fit model of the stacked one-halo, two-halo, and galaxy profile term from MCMC. The orange dashed and dotted lines show the best-fit intrinsic galaxy profile Σ gal and the P SF stack model. The purple and brown dashed lines mark the pixel size (7 ) and R200 value inferred from MICECAT. Open circles represent negative data points.   Figure 9. Measured (black data) and modeled (red) excess profile Σex (black data) of each case shown in Fig. 8. Note the excess profile is defined by the difference of the stacked profile and P SF stacked model (orange dotted line). Other lines are same as the ones shown in Fig. 8. values given by Table 3 of Wang et al. (2019), although when convolved, the total closely follows the PSF 9 . Our one-and two-halo clustering models, Σ 1h stack and Σ 2h stack , and the filtered signal Σ F stack , 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 onehalo 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 Σ 1h stack from MICECAT using the following steps: 1. Select the stacked target in the catalog using the same selection criteria; 9 In Wang et al. (2019), the values of R e,1 and R e,2 are reported in terms of x e,1 = R e,1 /R 200 and x e,2 = R e,2 /R 200 . R 200 is the projected virial radius of the host dark matter halo in angular units and its value for each sub-sample is given in Table 2. 2. For each target galaxy, generate a source map (using P SF 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 1h stack . The filtered signal term Σ F stack accounts for the loss of clustering signal from filtering. Σ F stack is the stacked profile on the 2D polynomial filtered map (the second term of Eq. 13), which can be modeled by filtering the simulated map from MICECAT. We model the two-halo term Σ 2h stack −Σ F stack after filtering with the following process: 1. Make a CIBER-sized mock image from all the catalog sources with the model P SF instr , and mask galaxy. The dark regions show the galaxy extensions associated with each galaxy, and the light blue and green regions show diffuse stars in the halos that are not tightly bound to any galaxy. The white parts with black dashed boundaries show the masked regions. The smaller galaxies without masks are fainter than the masking cutoff. The magenta stars and the orange regions show the stacked galaxy and its extension. The blue regions represent the one-halo term, and the green regions show the two-halo term contributed by emission from other halos. When stacking on a central galaxy, the one-halo term includes the satellite galaxy extensions beyond the masking radius as well as faint satellite galaxies and their stellar halos. When stacking on a satellite galaxy, the one-halo term includes the extensions of both the central and the satellite galaxies beyond their masks as well as the fainter satellite galaxies.
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 2h−F stack , is a model for Σ 2h stack − Σ F stack , which provides our template for the twohalo 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 threeparameter (R e,2 , A 1h , A 2h ) model can be written as We note that the one-and two-halo profiles already include the PSF convolution in our model.

Model Fitting
For each CIBER field and band, we fit the excess profile Eq. 6, to a three-parameter model Σ m ex (r, {R e,2 , A 1h , A 2h }) (Eq. 21) using a Markov Chain Monte Carlo (MCMC). We assume a Gaussian likelihood, which is given by where the inverse covariance C −1 ex is given by Eq. 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: 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 −4 R 200 , R 200 ], [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.

RESULTS
We show the MCMC results in Fig. 11 and Table 3, for all cases listed in Table 2. As a sanity check, we calculate the χ 2 value between the results from individual fields and the joint fit using 100 data points for each of the three parameters (5 fields × 10 mag bins × 2 bands). The resulting χ 2 values indicate our fit is internally consistent across the 5 CIBER fields. In Fig. 8 and Fig. 9, we show the stacked and excess profile data averaged over five fields, respectively, along with the marginalized one-halo, two-halo, and galaxy profile model from the joint fit. Fig. 12 shows the fitted intrinsic galaxy profile Σ gal (Eq. 20) and the one-and two-halo terms in the "total" magnitude bin, also averaged over five fields. The field-averaged profiles are only shown for visualization purposes; when we fit the data with MCMC, the information is combined in the likelihood function rather than in data space.

Missing Light in Galaxy Photometry
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 (Eq. 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 10 . With our galaxy profile, we can calculate the fraction of flux within the same radius (f petro ). The results are summarized in Table 4.
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 11 . 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 Sersic model to individual galaxies. Although their galaxy samples are at somewhat higher mass (M * ∼ 10 11 − 10 12 M ) and model magnitudes are fitted with a different functional form, we arrive at a similar fraction of missing flux.

Extended Stellar Halo
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 Fig. 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 11 M M * 10 12 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) Note-For the cases with less than a 2σ detection (95% confidence interval), we quote the 2σ upper bound. For detections, the +/− values enclose the 68% confidence interval. In 1.1 µm "total" bin, the 68% confidence interval of one-halo amplitude A 1h is 0.54 +0.42 −0.38 , approximately an 1σ detection. low-mass galaxies (9.2 M < logM * < 11.4 M ) 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 90 /R 50 is the ratio of the radii that contain 90% and 50% of the r -band Petrosian flux. 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 Fig. 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 Fig. 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 ex-  Table 2. The data points and error bars are the median and 68% confidence intervals from MCMC. Black data points show the joint fit from all five fields, with colored points for the individual fields. The gray horizontal lines in the middle and bottom panels mark A 1h = 1 and A 2h = 1, which are the clustering amplitudes given by MICECAT. The shaded regions show the total stack over all 17 < m1.1 < 20 galaxies.
ample, 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 Fig. 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 Fig. 14. We find that ∼ 25% of galaxy fluxes are from outside 20 kpc. The CIBER constraints shown in Fig. 13 are summarized in Table 6.

EBL from Extended Stellar Halos
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.  Figure 12. Fitted intrinsic galaxy profile Σ gal (Eq. 20) (orange), stacked one-halo (blue) and two-halo (green) profiles in the "total" magnitude bin averaged over five CIBER fields in the 1.1 µm (top) and 1.8 µm (bottom) bands. We convert the angular scale to physical units (kpc) using the median conversion factor inferred from MICECAT (Table 2). Solid lines and shaded regions indicate the median and 68% confidence interval of the joint fit constrained through MCMC, respectively.
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 Fig. 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.

Intra-halo Light Fraction
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 where m1.1<20 L is the total light associated with bright galaxies, and m1.1<20 L(> r 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 12 . For the bright sources, we define m1.1<20 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: This results in We show our constraints on f IHL , as a function of halo mass and redshift in Fig. 15 and 16, respectively. The Fraction of EBL intensity from galaxy extension as a function of rcut. This is estimated with the light profile fits from CIBER (this work) and HSC (Wang et al. 2019), and the stellar mass function from Muzzin et al. (2013).
halo masses associated with our galaxies are inferred from the MICECAT simulation and using the SDSS photometric redshifts. The CIBER data points shown in Fig. 15 and 16 are summarized in Table 7.
Note that the fraction of light beyond r cut (the numerator in Eq. 27) is shown in Fig. 13, where the higher redshift bins have slightly higher values. However, in Fig. 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 Fig. 17. Note that this observable quantity tracks the evolution of the one-halo luminosity but lacks the N eff term in Eq. 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 MICE-CAT simulation, which could be attributed to the unbound stars that are not included in MICECAT.
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(Gonzalez et al. , 2007Burke et al. 2015), and an analytical model (Purcell et al. 2007(Purcell et al. , 2008. Our results follow a more gradual redshift evolution trend than reported in massive clusters (Burke et al. 2015, see Fig. 16).

Color of the Galaxy Inner and Outer Regions
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 Eq. 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 f IHL (20 < r < 100 kpc) Figure 15. The IHL fraction fIHL as a function of halo mass. The IHL is defined by the light beyond a radius rcut around the galaxy. Here we consider three different rcut values: 10 kpc (left) and 20 kpc (right). Blue and red data points show the constraints from this work in the 1.1 µm and 1.8 µm bands, respectively. Dark and light green shaded regions denote the 68% and 95% variations among galaxies from an analytical model at z = 0 (Purcell et al. 2007(Purcell et al. , 2008. The ICL fraction in individual galaxy groups and clusters from Gonzalez et al. (2005Gonzalez et al. ( , 2007 and Burke et al. (2015) are shown in black and gray data points. The two downward arrows give upper limits for the Milky Way (Carollo et al. 2010) and Andromeda (M31) (Courteau et al. 2011). f IHL (20 < r < 100 kpc) Figure 16. fIHL constraints as in Fig. 15, but plotted as a function of redshift. The masses of the Burke et al. (2015) clusters are 100-1000× the halo masses associated with our galaxies.
A. & 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.
11.6. One-halo and Two-halo Clustering 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.
12. CONCLUSIONS By stacking galaxies from CIBER imaging data in two near-infrared bands (1.1 and 1.8 µm), we detect ex-  Figure 17. Ratio of the total one-halo term and stacked galaxy profile term from our stacking results (blue: 1.1 µm, red: 1.8 µm) compared with the MICECAT simulation (light blue: 1.1 µm, orange: 1.8 µm). We observe a somewhat stronger evolution, causing the fall-off of fIHL with redshift seen in Fig. 16. 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 A. & 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.
tended emission in galaxies. The galaxies being stacked (∼ 30, 000 galaxies in total) are split into five subsamples from SDSS spanning redshifts 0.2 z 0.5 and stellar masses 10 10.5 M M * 10 12 M , comparable to L * galaxies at this redshift. We jointly fit a model for the inherent galaxy light profile and large-scale oneand two-halo clustering.
With the galaxy profile, we estimate that ∼ 20% of total light is missing in galaxy photometry due to the use of limited apertures, in agreement with previous estimates from the literature. We do not detect a 1.1-1.8 µm color difference in the inner and outer region of our galaxy samples.
While we do not detect two-halo clustering, we detect one-halo clustering in the 1.8 µm band at 4σ significance over the full sample of galaxies. These results suggest nonlinear clustering could have a significant impact on modeling the IHL, but is not accounted for in previous fluctuation analysis by Zemcov et al. (2014). An IHL fluctuation model with one-halo clustering (e.g., Fernandez et al. 2010) is needed to fully account for the nonlinear clustering in IHL modeling.
The intrinsic galaxy profile fitted from our stacking analysis suggests ∼ 50%/25% of the total galaxy light resides at r > 10/20 kpc, respectively. This result is in agreement with previous HSC measurements at lower redshifts (0 z 0.25) and lower stellar masses (10 9.2 M < M * < 10 11.4 M ). The galaxy extension accounts for significant fraction of luminosity in L * galaxies, but falls off below M * ∼ 10 11 M . We extrapolate the fraction of extended galaxy light we measure to all galaxy mass scales and assuming a Schechter luminosity function, we find ∼ 30%/15% of the total galaxy light are from r > 10/20 kpc, respectively. We measure a moderate increase in f IHL with cosmic time, which we attribute to the decrease in one-halo contribution within the dark matter halo of our stacked samples. The previous fluctuation study using CIBER data (Zemcov et al. 2014) found that the IHL has comparable intensity to the IGL in the near-infrared EBL. While our study cannot constrain the whole IHL contribution to the EBL since we only study galaxies from a certain range of redshift and masses, our results suggest that ∼ L * galaxy at 0.2 z 0.5 have an extended light profile which contributes appreciable IHL to their host halos. As ∼ L * galaxies are the representative population, which contain most of the IGL emission, the flux from the extension, and the one-halo term present in our galaxy samples, both need to be properly accounted for in future EBL photometry and fluctuation measurements.  (Lang et al. 2010), LePHARE (Arnouts et al. 1999;Ilbert et al. 2006) APPENDIX A. EXTENSION AND IHL FRACTION Table 6 summarizes the fraction of light beyond 10 and 20 kpc, assuming our fitted light profile. These are the data presented in Fig. 13. Table 7 summarize the f IHL values with r cut= 10 and 20 kpc, assuming our fitted light profile and the onehalo contribution from the MICECAT. These are the data presented in Fig. 15 and 16. Note-These are the values shown in Fig. 13. The total row shows the weighted average of the five listed sub-samples.