# Probing Intra-Halo Light with Galaxy Stacking in CIBER Images

Yun-Ting Cheng<EMAIL_ADDRESS>Toshiaki Arai  Priyadarshini Bangale James
J. Bock  Tzu-Ching Chang Asantha Cooray  Richard M. Feder Phillip M. Korngut
Dae Hee Lee  Lunjun Liu  Toshio Matsumoto Shuji Matsuura Chi H. Nguyen Kei
Sano Kohji Tsumura Michael Zemcov

###### Abstract

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.

cosmology: observations – diffuse radiation –- large-scale structure of
universe

††software: astropy (Astropy Collaboration et al., 2013), emcee (Foreman-
Mackey et al., 2013), corner (Foreman-Mackey, 2016)

California Institute of Technology, 1200 E. California Boulevard, Pasadena, CA
91125, USA

Institute of Space and Astronautical Science, Japan Aerospace Exploration
Agency, Kanagawa 252-5210, Japan

Center for Detectors, School of Physics and Astronomy, Rochester Institute of
Technology, 1 Lomb Memorial Drive, Rochester, New York 14623, USA

California Institute of Technology, 1200 E. California Boulevard, Pasadena, CA
91125, USA Jet Propulsion Laboratory, California Institute of Technology, 4800
Oak Grove Drive, Pasadena, CA 91109, USA

Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove
Drive, Pasadena, CA 91109, USA California Institute of Technology, 1200 E.
California Boulevard, Pasadena, CA 91125, USA

Department of Physics and Astronomy, University of California, Irvine, CA
92697, U.S.A.

California Institute of Technology, 1200 E. California Boulevard, Pasadena, CA
91125, USA

California Institute of Technology, 1200 E. California Boulevard, Pasadena, CA
91125, USA

Korea Astronomy and Space Science Institute (KASI), Daejeon 305-348, Republic
of Korea

California Institute of Technology, 1200 E. California Boulevard, Pasadena, CA
91125, USA

Institute of Space and Astronautical Science, Japan Aerospace Exploration
Agency, Kanagawa 252-5210, Japan

School of Science and Technology, Kwansei Gakuin University, Sanda, Hyogo
669-1337, Japan

Center for Detectors, School of Physics and Astronomy, Rochester Institute of
Technology, 1 Lomb Memorial Drive, Rochester, New York 14623, USA

Kyushu Institute of Technology, 1-1 Sensui-cho, Tobata, Kitakyushu, Fukuoka
804-8550, Japan

Department of Natural Science, Faculty of Science and Engineering, Tokyo City
University, Tokyo, 158-8557, Japan

Center for Detectors, School of Physics and Astronomy, Rochester Institute of
Technology, 1 Lomb Memorial Drive, Rochester, New York 14623, USA Jet
Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove
Drive, Pasadena, CA 91109, USA

##  1 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 tens or hundreds of kilo-parsecs. The stellar
emission from this material is sometimes referred to as “intra-halo light”
(IHL), or in massive galaxy clusters as “intra-cluster light” (ICL).

The properties of stellar halos across a wide range of mass scales have been
extensively studied using analytical models (e.g., Purcell et al., 2007) and
N-body simulations (e.g., Bullock & Johnston, 2005; Conroy et al., 2007;
Rudick et al., 2009; Cooper et al., 2010, 2013, 2015; Rodriguez-Gomez et al.,
2016; Elias et al., 2018). Several observations have constrained the ICL
content in galaxy clusters (e.g., Lin & Mohr, 2004; Burke et al., 2015;
Gonzalez et al., 2005, 2007, 2005), as well as stellar halos in lower mass
systems by deeply imaging individual galaxies (e.g., Tal et al., 2009;
Martínez-Delgado et al., 2010; Abraham & van Dokkum, 2014; van Dokkum et al.,
2014; Huang et al., 2018) or through stacking (e.g., Zibetti et al., 2005;
D’Souza et al., 2014; Zhang et al., 2019; Wang et al., 2019).

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., 2020). 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., 2020) 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.

In this work, we further constrain the IHL using CIBER broad band 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 $\sim 30,000$ Sloan Digital Sky Survey (SDSS)
photometric galaxies at $z\sim$ 0.2 – 0.5 across five $2\times 2$ deg2 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 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 $\Lambda$CDM cosmology with
$n_{s}=0.97$, $\sigma_{8}=0.82$, $\Omega_{m}=0.26$, $\Omega_{b}=0.049$,
$\Omega_{\Lambda}=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.

Table 1: CIBER Observing Fields

Field Name | R.A. (∘) | Dec. (∘) | Time After Launch (sec) | Number of Frames Used | Integration Time (sec)
---|---|---|---|---|---
Elat10 | 191.50 | 8.25 | 387-436 | 24 | 42.72
Elat30 | 193.94 | 28.00 | 450-500 | 9 | 16.02
BootesB | 218.11 | 33.18 | 513-569 | 29 | 51.62
BootesA | 219.25 | 34.83 | 581-636 | 28 | 49.84
SWIRE(ELAIS-N1) | 241.53 | 54.77 | 655-705 | 25 | 44.50

Note. — We discard the beginning part of Elat30 field integration due to
pointing instability.

##  2 CIBER Experiment

CIBER111https://ciberrocket.github.io/ (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 2009, July
2010, March 2012, and June 2013. The first three CIBER flights were launched
at White Sands Missile Range, New Mexico on a Terrier-Black Brant IX rocket.
These flights reached $\sim$ 330 km apogee with $\sim$ 240 sec 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 sec) 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^{\prime\prime}$ to $2^{\circ}$. For
stacking, CIBER imaging can trace low surface 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\times 1024$
pixel HAWAII-1 HgCdTe detector. The two imagers are identical except for their
$\lambda/\Delta\lambda\sim 2$ filters, which are centered at $1.05$ and $1.79$
$\mu$m222In the first and second CIBER flights, the longer wavelength band is
centered at 1.56 $\mu$m, and thus it is named 1.6 $\mu$m band in previous
CIBER publications (Bock et al., 2013; Zemcov et al., 2014)..

During its fourth flight, CIBER observed eight science fields with $\sim 50$
sec integrations sampled at 1.78 sec 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 result, we only use the last 16 sec of this integration in
our analysis.

##  3 Data Processing

Figure 1:  Images from the SWIRE field in the 1.1 $\mu$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 $\sigma=35^{{}^{\prime\prime}}$
Gaussian kernel to highlight large-scale fluctuations. H: Image in Panel G
after subtracting a fitted 2-D polynomial, also shown smoothed with a
$\sigma=35^{{}^{\prime\prime}}$ 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.

Figure 2:  Same as Fig. 2 in the CIBER 1.8 $\mu$m band.

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

###  3.1 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 non-linear flux dependence, even for short integrations (Bock et
al., 2013). In any pixel that collects more than 5,000 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. 2 and Fig. 2 show the raw images of the
SWIRE field in the CIBER 1.1 and 1.8 $\mu$m bands, respectively.

###  3.2 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 $\sim 0.1$ e- s-1, less than 10 % of the sky
brightness. Panel B of Fig. 2 and Fig. 2 show the dark current maps of CIBER
1.1 and 1.8 $\mu$m bands, respectively.

###  3.3 Pixel Masks

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

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-$\sigma$
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-$\sigma$ 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, time-stream masks, and cosmic ray masks form the instrument mask. In
total, $\sim 10\%$ of pixels are removed by the instrument mask. Panel C of
Fig. 2 and Fig. 2 show the instrument masks in the SWIRE field of 1.1 and 1.8
$\mu$m band, respectively.

####  3.3.2 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_{\rm 1.1}$ and $m_{\rm 1.8}$, from these catalogs, as detailed in Sec. 4.
We mask all point sources brighter than $m_{\rm 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 $\nu I_{\nu}^{\rm th}=1$ nW m-2 sr-1 in the 1.1
$\mu$m band are masked. This choice of threshold value removes $\sim 50\%$ of
pixels in each field. We apply the same masking radius to 1.8 $\mu$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. 2 and Fig. 2 show the SWIRE field source
mask in the CIBER 1.1 and 1.8 $\mu$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-$\sigma$ pixel
clipping mask to identify additional outliers not flagged through the other
methods (e.g., from low-energy cosmic ray events or electronic tranisents).

###  3.4 Flat Fielding

CIBER images have a nonuniform 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 CIBER 1.1 and
1.8 $\mu$m bands are shown in the Panel E of Fig. 2 and Fig. 2, respectively.

###  3.5 Surface Brightness Calibration

Throughout this work, we use nW m-2 sr-1 for the units of surface brightness
($\nu I_{\nu}$). The calibration factor, $C$, that converts photocurrent (e-
s-1) to intensity (nW m-2 sr-1) is derived in the following steps:

  1. 1.

Take the raw images, subtract the dark current template, correct for the flat
field, and apply the instrument and source masks;

  2. 2.

Subtract the mean photocurrent in the unmasked region.

  3. 3.

For each star in the Pan-STARRS catalog, calculate the flux $\nu F_{\nu}$ in
CIBER bands from $m_{\rm 1.1}$ and $m_{\rm 1.8}$.

  4. 4.

Sum the photocurrent in a 5$\times$5 stamp centered on the source
position333We have tested that using 3$\times$3, 5$\times$5, or 7$\times$7
stamp size gives consistent results. Our beam size is approximately twice of
the pixel size, so a 3$\times$3 stamp already has enclosed most of the flux
from a point source..

  5. 5.

Repeat step (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_{\rm 1.1}<16$ for the 1.1
$\mu$m band, and $13.5<m_{\rm 1.1}<17$ for the 1.8 $\mu$m band. These
magnitude ranges are chosen such that the brightest sources that saturate the
detectors (even after non-linear 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. 2 and Fig. 2 show the SWIRE field images masked
by instrument masks at 1.1 and 1.8 $\mu$m, respectively, after flat fielding
and calibration.

###  3.6 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 star light (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. 2 and
Fig. 2 show the mean-subtracted and masked SWIRE images at 1.1 and 1.8 $\mu$m,
respectively. To highlight the large-scale fluctuations, we smooth the images
with a $\sigma=35^{{}^{\prime\prime}}$ Gaussian kernel.

###  3.7 Image Filtering

Although the ZL signal is smooth, a flat-field estimation error may induce a
nonuniform 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 3rd/5th order 2-D polynomial function for the
1.1/1.8 $\mu$m images to filter out any residual large-scale variations (Panel
H of Fig. 2 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).

##  4 External Catalogs

Throughout this work, we used several external source catalogs for (1) masking
bright foreground sources (Sec. 3.3.2); (2) calibration (Sec. 3.5); (3)
modeling the PSF by stacking bright stars in the fields (Sec. 7); and (4)
selecting galaxies for stacking (Sec. 8).

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 $\sim
50^{{}^{\prime\prime}}$ misalignment between the 1.1 and 1.8 $\mu$m images as
they are produced by different telescopes, their astrometry is solved
separately.

###  4.1 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\sim 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_{\rm 1.1}$ and $m_{\rm
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.

###  4.2 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_{\rm 1.1}$ ($m_{\rm 1.8}$)
is derived by linear extrapolation with the 2MASS photometric fluxes in J and
H (H and Ks) bands, respectively. We also use bright stars in 2MASS for
modeling the PSF (see Sec. 7).

###  4.3 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, 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 (Sec. 8.1), as
well as selecting stars for stacking to model the PSF (Sec. 7).

###  4.4 SWIRE Photometric Redshift Catalog

Rowan-Robinson et al. (2008, 2013) performed SED fitting on $\sim 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).

###  4.5 Gaia

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 (Sec. 7.2) and cleaning residual stars
in the galaxy sample selected by SDSS (Sec. 8.1).

###  4.6 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\leqslant 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 $\sim$ 200 clusters from Wen et al.
(2012) over the five CIBER fields.

##  5 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 30723 Mpc3h-3 cubic co-moving box. The dark
matter halos are resolved down to $\sim 3\times 10^{10}M_{\odot}h^{-1}$.

MICECAT is a mock catalog that simulates ideal observations of a 5000 deg2
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 (Carretero et al.,
2015). MICECAT simulates a mass-limited sample complete to $m_{i}\sim 22$ and
$m_{i}\sim 24$ at $z\simeq 0.5$ and $z\simeq 0.9$, respectively (Crocce et
al., 2015). The MICECAT mocks are large enough to permit us to generate up to
$\sim 10^{3}$ independent CIBER field-sized ($2\times 2$ deg2) mock catalogs.
We use modeled magnitudes from MICECAT in Euclid NISP Y and H bands for CIBER
$m_{\rm 1.1}$ and $m_{\rm 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 non-linear
(one-halo) 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. 1.

Select a list of stacking target sources from external catalogs.

  2. 2.

Re-grid each pixel into $N_{\rm sub}\times N_{\rm sub}$ sub-pixels (we use
$N_{\rm sub}=10$ in this work). The intensities of all sub-pixels are assigned
to the same value as the native pixel without interpolation.

  3. 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. 4.

Crop an $N_{\rm size}\times N_{\rm size}$ (at sub-pixel resolution) stamp
centered on the target source. We choose $N_{\rm size}=2401$ in this work,
which corresponds to a $28^{\prime}\times 28^{\prime}$ stamp.

  5. 5.

Repeat steps 3 and 4 for all target sources, average the stamps, and return
the final stacked 2-D image $\Sigma_{\rm stack}(\mathbf{r})$.

The stacked profile $\Sigma_{\rm stack}$ is a convolution of the intrinsic
source profile, $\Sigma_{\rm src}$, the instrument PSF ($PSF_{\rm
instr}$)444Instrument PSF includes all effects from the optics, detector
array, and pointing jitter during the integration., and the pixel function
$PSF_{\rm pix}$:

$\begin{split}\Sigma_{\rm stack}(\mathbf{r})=&\left[\Sigma_{\rm
src}(\mathbf{r})\circledast PSF_{\rm instr}(\mathbf{r})\right]\circledast
PSF_{\rm pix}(\mathbf{r})\\\ =&\Sigma_{\rm src}(\mathbf{r})\circledast
PSF_{\rm stack}(\mathbf{r}),\end{split}$ (1)

where $\mathbf{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_{\rm
stack}(\mathbf{r})\equiv PSF_{\rm instr}(\mathbf{r})\circledast PSF_{\rm
pix}(\mathbf{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 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):

$\displaystyle PSF_{\rm pix}(\mathbf{r})=\left\\{\begin{array}[]{cc}(N_{\rm
sub}-x)(N_{\rm sub}-y)&\\\ \quad\quad\quad\quad\quad\quad\text{if
}|x|,|y|<N_{\rm sub}\\\ 0\quad\quad\text{otherwise}\end{array}\right.$ (5)

As a practical matter, $PSF_{\rm pix}$ can be determined through simulations.
$PSF_{\rm stack}(\mathbf{r})$ can be measured by stacking stars in the field,
where $\Sigma_{\rm src}(\mathbf{r})$ is a delta function, so $\Sigma_{\rm
stack}(\mathbf{r})=PSF_{\rm stack}(\mathbf{r})$. Note that the expression in
the second line of Eq. 1 implies that the intrinsic profile $\Sigma_{\rm
src}(\mathbf{r})$ can be obtained from the stacked profile $\Sigma_{\rm
stack}(\mathbf{r})$ with the knowledge of $PSF_{\rm stack}(\mathbf{r})$,
instead of determining $PSF_{\rm instr}(\mathbf{r})$.

We perform stacking and PSF modeling separately for each field, since
$PSF_{\rm instr}$ is slightly different across the fields due to the varying
pointing performance of the altitude control system during each integration
(c.f. top panel of Fig. 5). After obtaining the 2-D stacked images, we bin
them into 25 logarithmically-spaced 1-D 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.

###  6.2 Covariance Matrix of Stacking Profile

The covariance matrix of the binned 1-D radial stacked profile is calculated
with a jackknife resampling technique. For each stack, we split sources into
$N_{\rm J}=64$ sub-groups based on their spatial coordinates in the image. The
CIBER imager arrays have 1024$\times$1024 pixels, and thus each sub-group
corresponds to sources in a $128\times 128$ pixel sub-region on the array. The
radial profile of the k-th jackknife sample, $\Sigma^{k}_{\rm 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

$\begin{split}&C_{\rm stack}(r_{i},r_{j})=\frac{N_{\rm J}-1}{N_{\rm
J}}\sum_{k=1}^{N_{\rm J}}\left[\Delta\Sigma^{k}_{\rm
stack}(r_{i})\cdot\Delta\Sigma^{k}_{\rm stack}(r_{j})\right]\\\
&\Delta\Sigma^{k}_{\rm stack}(r_{i})\equiv\Sigma^{k}_{\rm
stack}(r_{i})-\Sigma_{\rm stack}(r_{i})\\\ &\Delta\Sigma^{k}_{\rm
stack}(r_{j})\equiv\Sigma^{k}_{\rm stack}(r_{j})-\Sigma_{\rm
stack}(r_{j}),\end{split}$ (6)

where $\Sigma_{\rm 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 ($\ll 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 k-th bootstrap
sample, $\Sigma_{\rm 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

$C_{\rm stack}(r_{i},r_{j})=\frac{1}{N_{\rm B}-1}\sum_{k=1}^{N_{\rm
B}}\left[\Delta\Sigma^{k}_{\rm stack}(r_{i})\cdot\Delta\Sigma^{k}_{\rm
stack}(r_{j})\right].$ (7)

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^{\prime\prime}$ 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.

##  7 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 $PSF_{\rm 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 $PSF_{\rm 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.

###  7.1 Modeling $PSF_{\rm 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
$PSF_{\rm stack}$ robustly on both small and large scales, we construct an
overall star profile from three brightness bins. For the core region
($r<22^{\prime\prime}$), we stack $13<m_{\rm 1.1}<14$ sources in the field;
for intermediate scales, $22^{\prime\prime}<r<40^{\prime\prime}$, we fit a
slope to the stacking profile of $9<m_{J}^{\rm 2MASS}<10$ sources; for outer
radii, we fit another slope to the stacking profile of the brightest
$4<m_{J}^{\rm 2MASS}<9$ sources, and connect the two slopes at
$r=40^{\prime\prime}$ ($m_{J}^{\rm 2MASS}$ 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_{J}^{\rm 2MASS}<9$, $9<m_{J}^{\rm
2MASS}<10$, $13<m_{\rm 1.1}<14$) are statistically consistent. The top panel
of Fig. 3 shows $PSF_{\rm stack}$ from the SWIRE field in the 1.1 $\mu$m band.
The top panel of Fig. 5 shows $PSF_{\rm 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.

###  7.2 Validating $PSF_{\rm 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_{\rm 1.1}<20$), and compare these star stacking profiles with our
$PSF_{\rm stack}$ model.

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

  1. 1.

The source has a parallax measurement $>2\times 10^{-4}$ mas (i.e., distance
$<$ 5 kpc).

  2. 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. 3.

No SDSS galaxies within $0.7^{\prime\prime}$ (sub-pixel grid size) radius
around the source.

  4. 4.

We classify SDSS stars and galaxies using 10 pairs of magnitude differences
between the five 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_{\rm 1.1}<20$,
and compare their stacking profile with the $PSF_{\rm stack}$ model. These
stars span the same brightness range used for galaxy stacking. We down-sample
original 25 radial bins to 15 bins (7 bins for $16<m_{\rm 1.1}<17$ case),
following the same binning used for the galaxy stacking profile (Sec. 8.4).
The results in the 1.1 $\mu$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
$PSF_{\rm stack}$ model. The errors are propagated from the covariance of the
$PSF_{\rm stack}$ model and Gaia star stacks. We also show the $\chi^{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.

Figure 3:  We illustrate the process of constructing and validating the
$PSF_{\rm stack}(r)$ model, in the 1.1 $\mu$m band SWIRE field. Top: star
stacking profile in three different brightness bins (blue, orange, and green),
and the combined $PSF_{\rm 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 $PSF_{\rm
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
$PSF_{\rm stack}(r)$ model. Bottom: comparison of the $PSF_{\rm 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
$PSF_{\rm stack}(r)$ model agree closely with the star stacking profiles, as
shown in Fig. 4. Figure 4:  The difference of the $PSF_{\rm stack}(r)$ model
and the star stacking profiles in all five CIBER fields in the 1.1 $\mu$m
(left) and 1.8 (right) $\mu$m bands ($16<m_{\rm 1.1}<17$ (blue), $17<m_{\rm
1.1}<18$ (orange), $18<m_{\rm 1.1}<19$ (green), and $19<m_{\rm 1.1}<20$
(red)). The $\chi^{2}$ values and their corresponding PTE given in the legend
are consistent with the model. The degrees of freedom for each case is simply
the number of radial bins. Open circles in the top and middle panels represent
negative data points.

###  7.3 Modeling $PSF_{\rm instr}$

Although knowledge of the instrument PSF is not required for reconstructing
the source profile $\Sigma_{\rm src}$ from the stacking profile $\Sigma_{\rm
stack}$, $PSF_{\rm 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. $PSF_{\rm instr}$ is also useful for determining
the masking radius for bright sources (Sec. 3.3.2).

$PSF_{\rm instr}$ is modeled as follows: first, we deconvolve $PSF_{\rm
pix}(\mathbf{r})$ (Eq. 5) from the $PSF_{\rm stack}(\mathbf{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_{\rm instr}$, we fit a $\beta$
model (Cavaliere & Fusco-Femiano, 1978) to the 1-D profile of the deconvolved
image:

$PSF_{\rm
instr}(r)=\left(1+\left(\frac{r}{r_{c}}\right)^{2}\right)^{-3\beta/2}.$ (8)

Though not physically motivated, we find $\beta$ 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_{\rm stack}$.

The bottom panel of Fig. 5 illustrates this procedure in the 1.1 $\mu$m band
of the SWIRE field. The $PSF_{\rm stack}$ model, obtained from star stacks in
three different brightness bins, matches the $\beta$ model of $PSF_{\rm
instr}$ convolved with the pixel function $PSF_{\rm pix}$ (Eq. 5). Our
instrument PSF has comparable size to a pixel (FWHM $\sim 7^{\prime\prime}$).

Figure 5:  Top: $PSF_{\rm stack}$ model for each of the five fields in the 1.1
$\mu$m (solid) and 1.8 $\mu$m (dashed) bands. The variation across fields is
due to the difference in pointing stability. Bottom: demonstration of the
$PSF_{\rm instr}$ reconstruction process. Black data points show the $PSF_{\rm
stack}$ model in the 1.1 $\mu$m band SWIRE field, derived from splicing the
star stacking profile in three different brightness bins (c.f. Fig. 3 top
panel). The blue line is the $PSF_{\rm instr}$ model derived from fitting a
$\beta$ model to $PSF_{\rm stack}$ after deconvolving $PSF_{\rm pix}$ with the
Richardson-Lucy deconvolution algorithm. The orange line shows the convolution
of $PSF_{\rm instr}$ with $PSF_{\rm pix}$ matching the $PSF_{\rm stack}$
model, as a consistency check. Our model for $PSF_{\rm instr}$ is in agreement
with data for $r\lesssim 30^{\prime\prime}$. Our analysis is not susceptible
to the moderate error at larger radii, as $PSF_{\rm instr}$ is only used for
characterizing the clustering signal from nearby galaxies.

##  8 Galaxy Stacking

We stack galaxies within magnitude ranges $16<m_{\rm 1.1}<20$, divided into
several sub-samples spanning $\Delta m_{\rm 1.1}=1$. Our choice of magnitude
bins optimizes the SNR on the stacks, giving sufficient sample sizes for each
source brightness.

###  8.1 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. This criteria prevents 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^{\prime\prime}$, 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_{\odot}$ 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 excludes cluster members from our stacking
sample. Stacking on cluster members introduces extra non-linear 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 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^{14}M_{\odot}$ (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 $\sim 3-5\times$ at 100 arcsec
radius just by excluding galaxies near clusters following our criteria.

###  8.2 Stacking Sub-samples

For the SDSS galaxies within $16<m_{\rm 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 $\mu$m flux in four bins: $16<m_{\rm 1.1}<17$, $17<m_{\rm
1.1}<18$, $18<m_{\rm 1.1}<19$, and $19<m_{\rm 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_{\rm
1.1}<20$ sources to achieve better large-scale sensitivity.

The second set is defined by both the 1.1 $\mu$m apparent magnitude $m_{1.1}$
and the absolute magnitude $M_{1.1}$: $M_{1.1}=m_{1.1}-DM(z)+2.5{\rm
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 $\sim 10^{10.5}$ to
$10^{12}M_{\odot}$, i.e., $\sim L_{*}$ galaxies at this redshift (Muzzin et
al., 2013).

We also infer the host dark matter halo mass associated with our galaxy
samples 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 results are given in Table 2. We derive
the corresponding virial radius, $R_{\rm 200}$ (in co-moving and angular
units), in Table 2. The virial radius is calculated from $R_{\rm
200}=[3M_{h}/(4\pi\cdot 200\rho_{c})]^{1/3}$, where $\rho_{c}$ is the critical
density. We infer the fraction of central galaxies, $f_{\rm cen}$, in each
sub-sample from MICECAT. The distribution of redshift, stellar mass, halo
mass, and virial radius of our sub-samples are summarized in Fig. 6 and Table
2.

Figure 6:  Left: redshift distributions of the 10 galaxy sub-samples used for
stacking. The redshifts are derived from SDSS photometry. Middle-top: stellar
mass distributions for the 5 apparent and absolute magnitude selected bins.
The stellar masses are inferred from Rowan-Robinson et al. (2013) for the
SWIRE field. Middle-bottom: halo mass distributions in 5 apparent and absolute
magnitude selected bins, modeled by applying the same selection criteria to
the MICECAT catalog. Right: distributions of virial radius in co-moving (top)
and observed angular (bottom) units. For visualization purposes, all curves
are normalized by the total number of sources in each sub-sample ($N_{\rm
tot}$). Table 2: Summary of the properties on each stacked galaxy sub-sample
with the $+/-$ values indicating the 68% interval ranges.

Name | Selection Criteria | $N_{\rm gal}$ | $z$ | log $M_{*}$ [$M_{\odot}$] | log $M_{h}$ [M⊙] |  $R_{\rm 200}$ [kpc] |  $R_{\rm 200}$ [arcsec] | $f_{\rm cen}$
---|---|---|---|---|---|---|---|---
mag bin #1 | $16<m_{\rm 1.1}<17$ | 129 | 0.18${}^{+0.04}_{-0.02}$ | 11.6${}^{+0.3}_{-0.3}$ | 12.5${}^{+0.7}_{-0.7}$ | 284${}^{+228}_{-104}$ | 97${}^{+32}_{-19}$ | 0.65
mag bin #2 | $17<m_{\rm 1.1}<18$ | 1173 | 0.21${}^{+0.07}_{-0.04}$ | 11.5${}^{+0.3}_{-0.4}$ | 12.5${}^{+0.7}_{-0.6}$ | 285${}^{+212}_{-105}$ | 66${}^{+38}_{-14}$ | 0.67
mag bin #3 | $18<m_{\rm 1.1}<19$ | 3465 | 0.27${}^{+0.09}_{-0.07}$ | 11.2${}^{+0.4}_{-0.3}$ | 12.4${}^{+0.8}_{-0.6}$ | 254${}^{+202}_{-85}$ | 43${}^{+26}_{-9}$ | 0.62
mag bin #4 | $19<m_{\rm 1.1}<20$ | 31157 | 0.42${}^{+0.17}_{-0.11}$ | 11.1${}^{+0.3}_{-0.5}$ | 12.3${}^{+0.8}_{-0.7}$ | 220${}^{+187}_{-75}$ | 29${}^{+22}_{-6}$ | 0.63
total | $17<m_{\rm 1.1}<20$ | 35795 | 0.40${}^{+0.17}_{-0.14}$ | 11.1${}^{+0.3}_{-0.4}$ | 12.3${}^{+0.8}_{-0.6}$ | 232${}^{+189}_{-80}$ | 40${}^{+31}_{-9}$ | 0.63
high-M/low-z |  $17<m_{\rm 1.1}<18$ $-23<M_{\rm 1.1}<-22$ | 743 | 0.22${}^{+0.04}_{-0.03}$ | 11.6${}^{+0.2}_{-0.4}$ | 12.5${}^{+0.6}_{-0.2}$ | 291${}^{+167}_{-55}$ | 64${}^{+33}_{-17}$ | 0.66
high-M/med-z |  $18<m_{\rm 1.1}<19$ $-23<M_{\rm 1.1}<-22$ | 1274 | 0.34${}^{+0.05}_{-0.05}$ | 11.4${}^{+0.3}_{-0.2}$ | 12.6${}^{+0.5}_{-0.4}$ | 285${}^{+149}_{-58}$ | 42${}^{+22}_{-8}$ | 0.62
high-M/high-z |  $19<m_{\rm 1.1}<20$ $-23<M_{\rm 1.1}<-22$ | 10916 | 0.54${}^{+0.10}_{-0.09}$ | 11.3${}^{+0.3}_{-0.3}$ | 12.6${}^{+0.5}_{-0.3}$ | 270${}^{+139}_{-53}$ | 28${}^{+13}_{-6}$ | 0.66
low-M/low-z |  $18<m_{\rm 1.1}<19$ $-22<M_{\rm 1.1}<-21$ | 1645 | 0.24${}^{+0.05}_{-0.03}$ | 11.1${}^{+0.3}_{-0.2}$ | 12.0${}^{+0.8}_{-0.2}$ | 191${}^{+173}_{-27}$ | 41${}^{+36}_{-7}$ | 0.57
low-M/med-z |  $19<m_{\rm 1.1}<20$ $-22<M_{\rm 1.1}<-21$ | 14730 | 0.38${}^{+0.05}_{-0.05}$ | 11.0${}^{+0.2}_{-0.4}$ | 12.1${}^{+0.8}_{-0.3}$ | 190${}^{+189}_{-32}$ | 27${}^{+27}_{-5}$ | 0.58

Note. — $N_{\rm 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 last four columns are derived by applying the same cuts
to a simulated catalog from MICECAT.

###  8.3 Galaxy Stacking Profile

Figure 7:  The stacked galaxy radial profile from the SWIRE field mag bin #2
in the 1.1 (left) and 1.8 $\mu$m (right) bands. Top: galaxy stacked profile
$\Sigma_{\rm stack}$ (black) and $PSF_{\rm stack}$ model (orange dashed),
scaled to match the innermost radial bin of $\Sigma_{\rm stack}$. The error
bars give the diagonal element of the covariance matrix derived by the
Jackknife method (described in Sec. 6). Middle: the excess profile
($\Sigma_{\rm ex}$, Eq. 9) for the case shown in the top row. The excess is
defined as the difference between the galaxy stacked profile and the $PSF_{\rm
stack}$ model, i.e., the difference of the black data from the orange curve in
the top row. Bottom: the field-averaged excess profile $\Sigma_{\rm 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 $R_{\rm 200}$ values inferred from MICECAT,
respectively. Open circles in all the plots represent negative values.

We calculate 1-D radial profiles from galaxy stacks by averaging pixels in
concentric annuli, as shown in Fig. 7 and Fig. 9. For comparison, we also plot
the expected profile of stacked point sources, $PSF_{\rm 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_{\rm stack}$ profile.

###  8.4 Excess Profile

We define an “excess profile” $\Sigma_{\rm ex}(r)$ as follows:

$\Sigma_{\rm ex}(r)=\Sigma_{\rm stack}(r)-A\cdot PSF_{\rm stack}(r),$ (9)

where the normalization factor $A$ is chosen such that $PSF_{\rm stack}$
matches $\Sigma_{\rm stack}$ at the innermost radial bin $r_{1}$, and thus by
construction, $\Sigma_{\rm ex}(r_{1})=0$, and $A\equiv\Sigma_{\rm
stack}(r_{1})/PSF_{\rm 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

$C_{\rm ex}=\Sigma_{\rm stack}(r_{1})^{2}\left[C_{\rm norm}\left(C_{\rm
stack}\right)+C_{\rm norm}\left(C_{\rm PSF}\right)\right],$ (10)

where $C_{\rm PSF}$ and $C_{\rm stack}$ are the covariance of $PSF_{\rm
stack}$ and $\Sigma_{\rm stack}$, respectively, and

$\begin{split}C_{\rm norm}&\left(C,\frac{\Sigma_{\rm
stack}(r_{i})}{\Sigma_{\rm stack}(r_{1})},\frac{\Sigma_{\rm
stack}(r_{j})}{\Sigma_{\rm stack}(r_{1})}\right)\\\ &=\frac{\Sigma_{\rm
stack}(r_{i})\Sigma_{\rm stack}(r_{j})}{\Sigma_{\rm stack}(r_{1})^{2}}\cdot\\\
&\left[\frac{C(r_{i},r_{j})}{\Sigma_{\rm stack}(r_{i})\Sigma_{\rm
stack}(r_{j})}-\frac{C(r_{i},r_{1})}{\Sigma_{\rm stack}(r_{i})\Sigma_{\rm
stack}(r_{1})}\right.\\\ &\left.-\frac{C(r_{j},r_{1})}{\Sigma_{\rm
stack}(r_{j})\Sigma_{\rm stack}(r_{1})}+\frac{C(r_{1},r_{1})}{\Sigma_{\rm
stack}(r_{1})^{2}}\right]\end{split}$ (11)

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

To fit a model to the measured $\Sigma_{\rm ex}$, we also need the inverse of
$C_{\rm ex}$. However, $C_{\rm 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 regions555Mag
bin # 1 is down-sampled to 7 radial bins as its degree of freedom is limited
by the small number of stacked sources.. After this down-sampling, we derive
the inverse covariance estimator by

$C_{\rm ex}^{-1}=\frac{N_{\rm J}-N_{\rm bin}-2}{N_{\rm J}-1}C_{\rm ex}^{*-1},$
(12)

where $N_{\rm J}=64$, the number of sub-groups used for estimating covariance,
and the number of bins $N_{\rm bin}=15$. $C_{\rm ex}^{*-1}$ is the direct
inverse of the $C_{\rm ex}$ matrix, and the pre-factor in Eq. 12 de-biases the
inverse covariance estimator, as our covariance matrix is derived from our
data (Hartlap et al., 2007)666For mag bin # 1, $N_{\rm bin}=7$, and $N_{\rm
J}=64$ is replaced by $N_{\rm B}=1000$ since we use bootstrap resampling
method in this case..

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 by the covariance (Eq. 10), and thhus 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 in parameter
fitting.

Figure 8: The stacked profile (black data) of each sub-sample stack averaged
over five CIBER fields in the 1.1 $\mu$m (top) and 1.8 $\mu$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 $\Sigma^{\rm gal}$ and the $PSF_{\rm stack}$
model. The purple and brown dashed lines mark the pixel size
(7${}^{{}^{\prime\prime}}$) and $R_{200}$ value inferred from MICECAT. Open
circles represent negative data points.

Figure 9:  The measured (black data) and modeled (red) excess profile
$\Sigma_{\rm ex}$ (black data) of each case shown in Fig. 9. Note the excess
profile is defined by the difference of the stacked profile and $PSF_{\rm
stacked}$ model (orange dotted line). other lines are same as the ones shown
in Fig. 9.

##  9 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), define fitted
profiles (Sec. 9.2), and introduce our model for each component of the stack
(Sec. 9.3). Finally, we describe the model fitting procedure in Sec. 9.4.

###  9.1 Components in Image Space

The raw CIBER image, $I_{\rm raw}$, can be expressed as777For clarification,
$\mathbf{x}$ denotes 2-D coordinate on CIBER images, and $\mathbf{r}$
represents the coordinate that has origin at the source center, which is used
in $PSF_{\rm instr}$ and stacked maps. Since we only consider 1-D radially
averaged profile, $\mathbf{r}$ is replaced by 1-D variable “$r$”.

$\begin{split}I_{\rm raw}(\mathbf{x})=&\left[I_{\rm sig}(\mathbf{x})+I_{\rm
LoS}(\mathbf{x})\right]\circledast PSF_{\rm instr}(\mathbf{r})\cdot
FF(\mathbf{x})\\\ &+I_{\rm DC}(\mathbf{x})+I_{\rm n}(\mathbf{x}),\end{split}$
(13)

where $\mathbf{x}$ is the 2-D pixel coordinate, $FF$ is the flat-field gain,
$I_{\rm DC}$ is the dark current map, and $I_{\rm n}$ is the read noise plus
photon noise. The sky emission is decomposed into $I_{\rm sig}$ and $I_{\rm
LoS}$ terms, where the first term accounts for the signal associated with
stacked galaxies, and $I_{\rm 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^{\prime}_{\rm raw}$:

$I^{\prime}_{\rm raw}(\mathbf{x})=\left[I_{\rm sig}(\mathbf{x})+I_{\rm
LoS}(\mathbf{x})\right]\circledast PSF_{\rm instr}(\mathbf{r})+I^{\prime}_{\rm
n}(\mathbf{x}),$ (14)

where $I^{\prime}_{\rm n}(\mathbf{x})=I_{\rm n}(\mathbf{x})/FF(\mathbf{x})$,
the instrument noise divided by the flat-field response. For simplicity, we
ignore the error in the flat-field estimator in Eq. 14. 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 by the Jackknife method (see Sec. 6.2). We define the
mask $M(\mathbf{x})$ as a binary function set to zero at masked pixels, and
one otherwise. The filtered map is expressed with
$\mathcal{F}\left[I^{\prime}_{\rm raw}(\mathbf{x}),M(\mathbf{x})\right]$,
which is a function of the input map $I^{\prime}_{\rm raw}(\mathbf{x})$ and
mask $M(\mathbf{x})$. As described in Sec. 3.7, we choose $\mathcal{F}$ to be
a 3rd (1.1 $\mu$m)/5th (1.8 $\mu$m) order 2-D polynomial function fitted to
the masked $I^{\prime}_{\rm raw}$ map888Note that the filter map $\mathcal{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(\mathbf{x})$ and
$B(\mathbf{x})$, and a mask $M(\mathbf{x})$,
$\mathcal{F}\left[A(\mathbf{x})+B(\mathbf{x}),M(\mathbf{x})\right]=\mathcal{F}\left[A(\mathbf{x}),M(\mathbf{x})\right]+\mathcal{F}\left[B(\mathbf{x}),M(\mathbf{x})\right]$..
The image used for stacking $I_{\rm map}$ can thus be written as

$\begin{split}I_{\rm map}(\mathbf{x})&=I^{\prime}_{\rm
raw}(\mathbf{x})M(\mathbf{x})-\mathcal{F}\left[I^{\prime}_{\rm
raw}(\mathbf{x}),M(\mathbf{x})\right]M(\mathbf{x})\\\ &=I_{\rm map}^{\rm
sig}(\mathbf{x})+I_{\rm map}^{\rm LoS}(\mathbf{x})+I_{\rm map}^{\rm
n^{\prime}}(\mathbf{x}),\end{split}$ (15)

where

$\begin{split}&I_{\rm map}^{\rm sig}(\mathbf{x})=\left[I_{\rm
sig}(\mathbf{x})\circledast PSF_{\rm instr}(\mathbf{r})\right]M(\mathbf{x})\\\
&-\mathcal{F}\left[I_{\rm sig}(\mathbf{x})\circledast PSF_{\rm
instr}(\mathbf{r}),M(\mathbf{x})\right]M(\mathbf{x}),\end{split}$ (16)

$\begin{split}&I_{\rm map}^{\rm LoS}(\mathbf{x})=\left[I_{\rm
LoS}(\mathbf{x})\circledast PSF_{\rm instr}(\mathbf{r})\right]M(\mathbf{x})\\\
&-\mathcal{F}\left[I_{\rm LoS}(\mathbf{x})\circledast PSF_{\rm
instr}(\mathbf{r}),M(\mathbf{x})\right]M(\mathbf{x}),\end{split}$ (17)

and

$\begin{split}&I_{\rm map}^{\rm n^{\prime}}(\mathbf{x})=I_{\rm
n}^{\prime}(\mathbf{x})M(\mathbf{x})\\\ &-\mathcal{F}\left[I_{\rm
n}^{\prime}(\mathbf{x}),M(\mathbf{x})\right]M(\mathbf{x}).\end{split}$ (18)

###  9.2 Components in the Stack

The stacked profile $\Sigma_{\rm stack}$ can be expressed as the sum of stacks
on the three maps in Eq. 15:

$\Sigma_{\rm stack}(r)=\Sigma_{\rm stack}^{\rm sig}(r)+\Sigma_{\rm stack}^{\rm
LoS}(r)+\Sigma_{\rm stack}^{\rm n}(r).$ (19)

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

We model the stacked galaxy profile as

$\begin{split}\Sigma_{\rm stack}^{\rm sig}(r)=\left[\right.\Sigma_{\rm
stack}^{\rm gal}(r)+\Sigma_{\rm stack}^{\rm 1h}(r)&+\Sigma_{\rm stack}^{\rm
2h}(r)\left.\right]\\\ &-\Sigma_{\rm stack}^{\mathcal{F}}(r),\end{split}$ (20)

where the first three terms are the signal terms, and the last term is the
filtered signal map in Eq. 16. The galaxy profile term, $\Sigma_{\rm
stack}^{\rm gal}$, represents the intrinsic galaxy profile, which includes the
galaxy shape and the extended stellar halo. We decompose the galaxy profile
term, $\Sigma_{\rm stack}^{\rm gal}$, into “core” and “extended” parts:

$\Sigma^{\rm gal}_{\rm stack}(r)=\Sigma_{\rm core}^{\rm gal}(r)+\Sigma_{\rm
ext}^{\rm gal}(r),$ (21)

where the core component is the integrated emission of the $PSF_{\rm stack}$
fitted to the stacking profile, i.e., the $A\cdot PSF_{\rm stack}$ term in Eq.
9, and the extended component is the rest of the galaxy emission:

$\begin{split}\Sigma_{\rm core}^{\rm gal}(r)=&\Sigma_{\rm
stack}(r)-\Sigma_{\rm ex}(r),\\\ \Sigma_{\rm ext}^{\rm gal}(r)=&\Sigma_{\rm
ex}(r)\\\ &-\left[\Sigma_{\rm stack}^{\rm 1h}(r)+\Sigma_{\rm stack}^{\rm
2h}(r)-\Sigma_{\rm stack}^{\mathcal{F}}(r)\right].\end{split}$ (22)

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.
$\Sigma_{\rm stack}^{\rm 1h}$ and $\Sigma_{\rm stack}^{\rm 2h}$ 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 modelled
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 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 $\sim 60$% of stacked
galaxies are central galaxies, and $\sim 40$% are satellite galaxies (See
Table 2).

Figure 10:  Illustration of the components in our model when stacking on a
central (top) or a satellite (bottom) 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.

###  9.3 Modeling the Stacked Galaxy Profile

The stacked galaxy profile $\Sigma_{\rm stack}^{\rm gal}(r)=\Sigma^{\rm
gal}(r)\circledast PSF_{\rm stack}(r)$, is the intrinsic galaxy profile
$\Sigma^{\rm gal}$, including the galaxy shape and the extended stellar halo,
convolved with $PSF_{\rm stack}$. Following Wang et al. (2019), we model
$\Sigma^{\rm gal}$ with a double Sersic function:

$\begin{split}&\Sigma^{\rm gal}(r)=A^{\rm gal}\left(10^{I_{e,1}}{\rm
exp}\left\\{-b_{n1}\left[\left(r/R_{e,1}\right)^{1/n_{1}}-1\right]\right\\}\right.\\\
&+\left.10^{I_{e,2}}{\rm
exp}\left\\{-b_{n2}\left[\left(r/\left(R_{e,1}+R_{e,2}\right)\right)^{1/n_{2}}-1\right]\right\\}\right).\end{split}$
(23)

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 PSF999In 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-samples is given in Table 2..

Our one- and two-halo clustering models, $\Sigma^{\rm 1h}_{\rm stack}$ and
$\Sigma^{\rm 2h}_{\rm stack}$, and the filtered signal
$\Sigma^{\rm\mathcal{F}}_{\rm 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 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 $\Sigma_{\rm stack}^{\rm 1h}$ from MICECAT using
the following steps:

  1. 1.

Select the stacked target in the catalog using the same selection criteria.

  2. 2.

For each target galaxy, generate a source map (using $PSF_{\rm instr}$) for
all galaxies residing in the same halo except for the target galaxy.

  3. 3.

Generate a source mask using the same prescription as our data.

  4. 4.

Stack on the target source position.

  5. 5.

Iterate steps (2)-(4) for all target sources.

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

  1. 1.

Make a CIBER-sized mock image from all the catalog sources with the model
$PSF_{\rm instr}$, and mask it with a source mask generated using the same
masking process applied to the data.

  2. 2.

Fit and subtract a 2-D polynomial map to the image.

  3. 3.

Select the stacked target in the catalog using the same selection criteria as
the real sources.

  4. 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. 5.

Iterate on step (4) to derive a stacked profile of the filtered two-halo
signal.

The resulting stacked profile, $T_{\rm stack}^{\rm 2h-\mathcal{F}}$, is a
model for $\Sigma^{\rm 2h}_{\rm stack}-\Sigma^{\rm\mathcal{F}}_{\rm stack}$,
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_{\rm e,2},A_{\rm 1h},A_{\rm 2h}$) model can be written as

$\begin{split}\Sigma_{\rm stack}&(r,\left\\{R_{\rm e,2},A_{\rm 1h},A_{\rm
2h}\right\\})\\\ &=\Sigma^{\rm gal}(r,\left\\{R_{\rm e,2}\right\\})\circledast
PSF_{\rm stack}(r)\\\ &+A_{\rm 1h}T_{\rm stack}^{\rm 1h}(r)+A_{\rm 2h}T_{\rm
stack}^{\rm 2h-\mathcal{F}}(r).\end{split}$ (24)

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

###  9.4 Model Fitting

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

$\begin{split}\chi^{2}&=\left(\Sigma^{\rm d}_{\rm ex}-\Sigma^{\rm m}_{\rm
ex}\right)^{T}C^{-1}_{\rm ex}\left(\Sigma^{\rm d}_{\rm ex}-\Sigma^{\rm m}_{\rm
ex}\right)\\\ {\rm ln}\mathcal{L}&=-\frac{1}{2}\chi^{2}-\frac{1}{2}{\rm
ln}\left|C_{\rm ex}\right|+{\rm constant},\end{split}$ (25)

where the inverse covariance $C^{-1}_{\rm ex}$ is given by Eq. 12.

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:

${\rm ln}\mathcal{L}=\sum_{i=1}^{N_{\rm field}}{\rm ln}\mathcal{L}_{i}$ (26)

where $N_{\rm 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_{\rm
e,2}$, $A_{\rm 1h}$, and $A_{\rm 2h}$ in the range of [10${}^{-4}R_{\rm 200}$,
$R_{\rm 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$.

##  10 Results

Table 3: Summary of parameter constraints from the joint fit in each case
listed in Table 2. For the cases with less than a 2$\sigma$ detection (95%
confidence interval), we quote the 2$\sigma$ upper bound. For detections, the
$+/-$ values enclose the 68% confidence interval.

| 1.1 $\mu$m | 1.1 $\mu$m | 1.1 $\mu$m | 1.8 $\mu$m | 1.8 $\mu$m | 1.8 $\mu$m
---|---|---|---|---|---|---
Name |  $R_{e,2}$ [arcsec] | $A_{\rm 1h}$ | $A_{\rm 2h}$ |  $R_{e,2}$ [arcsec] | $A_{\rm 1h}$ | $A_{\rm 2h}$
mag bin #1 | $<2.76$ | $<6.06$ | $<48.91$ | $<2.53$ | $<5.72$ | $<58.05$
mag bin #2 | 2.25${}^{+0.14}_{-0.23}$ | $<4.70$ | $<24.22$ | 1.94${}^{+0.12}_{-0.16}$ | $<3.44$ | $<24.76$
mag bin #3 | 1.85${}^{+0.17}_{-0.28}$ | $<4.18$ | $<18.94$ | 1.94${}^{+0.16}_{-0.16}$ | $<2.96$ | $<18.30$
mag bin #4 | 1.85${}^{+0.25}_{-0.21}$ | $<1.16$ | $<6.87$ | 1.63${}^{+0.21}_{-0.14}$ | 0.77${}^{+0.23}_{-0.23}$ | $<6.59$
total | 1.98${}^{+0.17}_{-0.17}$ |  $<1.41$* | $<7.30$ | 1.85${}^{+0.08}_{-0.15}$ | 1.01${}^{+0.24}_{-0.24}$ | $<6.86$
high-M/low-z | 2.30${}^{+0.16}_{-0.29}$ | $<4.76$ | $<25.58$ | 2.17${}^{+0.18}_{-0.18}$ | $<4.2$ | $<33.10$
high-M/med-z | 2.27${}^{+0.37}_{-0.32}$ | $<6.42$ | $<19.53$ | 2.22${}^{+0.19}_{-0.28}$ | 3.37${}^{+1.99}_{-1.17}$ | $<22.76$
high-M/high-z | 1.98${}^{+0.30}_{-0.44}$ | $<1.88$ | $<9.08$ | 1.85${}^{+0.26}_{-0.22}$ | 1.39${}^{+0.43}_{-0.35}$ | $<6.19$
low-M/low-z | 1.98${}^{+0.18}_{-0.30}$ | $<3.18$ | $<16.38$ | 1.89${}^{+0.21}_{-0.17}$ | $<2.77$ | $<17.65$
low-M/med-z | 1.67${}^{+0.29}_{-0.36}$ | $<1.30$ | $<11.30$ | 1.50${}^{+0.21}_{-0.24}$ | $<1.01$ | $<7.58$

Note. — In 1.1 $\mu$m “total” bin, the 68% confidence interval of one-halo
amplitude $A_{\rm 1h}$ is 0.54${}^{+0.42}_{-0.38}$, approximately an $1\sigma$
detection.

Figure 11:  Marginalized parameter constraints from MCMC for each case listed
in 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_{\rm 1h}=1$ and $A_{\rm 2h}=1$,
which are the clustering amplitudes given by MICECAT. The shaded regions show
the total stack over all $17<m_{\rm 1.1}<20$ galaxies.

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 $\chi^{2}$ value between the results
from individual fields and the joint fit using 100 data points for each of the
three parameters (5 fields $\times$ 10 mag bins $\times$ 2 bands). The
resulting $\chi^{2}$ values indicate our fit is internally consistent across
the 5 CIBER fields. In Fig. 9 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 $\Sigma^{\rm gal}$ (Eq. 23)
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.

Figure 12:  The fitted intrinsic galaxy profile $\Sigma^{\rm gal}$ (Eq. 23)
(orange), stacked one-halo (blue) and two-halo (green) profiles in the “total”
magnitude bin averaged over five CIBER fields in the 1.1 $\mu$m (top) and 1.8
$\mu$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.

##  11 Discussion

###  11.1 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_{\rm core}\equiv L_{\rm core}/(L_{\rm core}+L_{\rm ext})$,
where $L_{\rm core}$ and $L_{\rm ext}$ are the total flux in the core and
extension profile (Eq. 22), 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 Perosian 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
radius101010https://www.sdss.org/dr12/algorithms/magnitudes/#mag_petro. With
our galaxy profile, we can calculate the fraction of flux within the same
radius ($f_{\rm petro}$). The results are summarized in Table 5.

We also estimate the missing light fraction with the ‘model magnitude’ given
in SDSS ($f_{\rm 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 fitting111111See
https://www.sdss.org/dr12/algorithms/magnitudes/ for the detailed descriptions
on model magnitude.. 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 $\sim$ 20% of the total galaxy light, a deficit
detected at $\sim 7\sigma$ ($\sim 4\sigma$) 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 ($\sim$ 10 to 20
%). Interestingly, Tal & van Dokkum (2011) probed the radial profile of $z\sim
0.34$ luminous red galaxies (LRGs) in SDSS with a stacking analysis, and they
also found $\sim$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_{*}\sim 10^{11}-10^{12}M_{\odot}$), and model
magnitudes are fitted with a different functional form, we arrive at a similar
fraction of missing flux.

Table 5: 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 $\mu$m | 1.1 $\mu$m | 1.1 $\mu$m | 1.8 $\mu$m | 1.8 $\mu$m | 1.8 $\mu$m
---|---|---|---|---|---|---
Name | $f_{\rm core}$ | $f_{\rm petro}$ | $f_{\rm model}$ | $f_{\rm core}$ | $f_{\rm petro}$ | $f_{\rm model}$
high-M/low-z | 0.79${}^{+0.04}_{-0.02}$ | 0.78${}^{+0.08}_{-0.10}$ | 0.84${}^{+0.11}_{-0.12}$ | 0.81${}^{+0.02}_{-0.02}$ | 0.80${}^{+0.07}_{-0.10}$ | 0.85${}^{+0.10}_{-0.12}$
high-M/med-z | 0.81${}^{+0.04}_{-0.05}$ | 0.74${}^{+0.07}_{-0.13}$ | 0.78${}^{+0.08}_{-0.15}$ | 0.83${}^{+0.04}_{-0.03}$ | 0.75${}^{+0.08}_{-0.11}$ | 0.78${}^{+0.10}_{-0.13}$
high-M/high-z | 0.86${}^{+0.06}_{-0.04}$ | 0.73${}^{+0.07}_{-0.16}$ | 0.77${}^{+0.15}_{-0.19}$ | 0.89${}^{+0.03}_{-0.04}$ | 0.75${}^{+0.07}_{-0.16}$ | 0.79${}^{+0.16}_{-0.18}$
low-M/low-z | 0.84${}^{+0.04}_{-0.02}$ | 0.78${}^{+0.05}_{-0.11}$ | 0.80${}^{+0.10}_{-0.12}$ | 0.85${}^{+0.02}_{-0.03}$ | 0.79${}^{+0.05}_{-0.11}$ | 0.81${}^{+0.10}_{-0.12}$
low-M/med-z | 0.89${}^{+0.04}_{-0.04}$ | 0.78${}^{+0.06}_{-0.16}$ | 0.80${}^{+0.09}_{-0.16}$ | 0.92${}^{+0.03}_{-0.03}$ | 0.80${}^{+0.06}_{-0.15}$ | 0.83${}^{+0.10}_{-0.14}$
total | 0.83${}^{+0.02}_{-0.01}$ | 0.77${}^{+0.03}_{-0.06}$ | 0.80${}^{+0.05}_{-0.06}$ | 0.86${}^{+0.01}_{-0.01}$ | 0.78${}^{+0.03}_{-0.05}$ | 0.81${}^{+0.05}_{-0.06}$

###  11.2 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 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) has studied individual stellar halos out to
100 kpc in more massive galaxies ($10^{11}M_{\odot}\lesssim M_{*}\lesssim
10^{12}M_{\odot}$) at higher redshifts ($z\sim 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\lesssim z\lesssim 0.25$) low-mass galaxies ($9.2\,M_{\odot}<{\rm
log}M_{*}<11.4\,M_{\odot}$) with a stacking analysis on HSC images in r-band.
They stacked galaxies out to $\sim 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_{\rm 90}/R_{\rm 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. 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 $\sim$
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
grey 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 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
$\sim 10$ to $\sim 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.

Figure 13:  Fraction of flux between 10 and 100 kpc from the galaxy profile
derived from CIBER stacking (this work) in the 1.1 (blue) and 1.8 (red) $\mu$m
bands and from HSC stacking. (Wang et al., 2019). The HSC stacking is
performed on low and high concentration populations ($C<2.6$ and $C>2.6$) at
optical wavelengths (r band). The horizontal error bars define the lower and
upper bounds of the stellar mass of each stacking sample. The grey line and
the shaded regions are the median, 16th, and 84th percentile of ex-situ
stellar mass fraction at $z=0$ from Illustris simulations (Rodriguez-Gomez et
al., 2016). The shaded region shows the variance between individual galaxies
in Illustris, whereas for CIBER and HSC, the error bars represent the standard
error on the mean value.

###  11.3 Intra-halo Light Fraction

The fraction of the total emission from a dark matter halo associated with
IHL, $f_{\rm 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 use the extended
galaxy emission beyond 10 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_{\rm IHL}$ can be expressed as

$f_{\rm IHL}=\frac{\sum_{m_{1.1}<20}L(>10\,{\rm
kpc})}{\sum_{m_{1.1}<20}L+\sum_{\rm faint}L},$ (27)

where $\sum_{m_{1.1}<20}L$ is the total light associate with bright galaxies,
and $\sum_{m_{1.1}<20}L(>10\,{\rm kpc})$ is the part of bright galaxy emission
beyond 10 kpc. $\sum_{\rm faint}L$ represents the light from faint galaxies as
well as the unbound stars in the halo, captured in the 1-halo luminosity. Note
that we conservatively assume the 1-halo luminosity arises entirely from
faint, gravitationally bound galaxies. However it is certainly true that some
1-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 $\sum_{\rm faint}L$ can
be described by the total emission in the one-halo term, $L_{\rm 1h}$121212Our
one-halo model also includes the outskirts of bright sources beyond the mask,
but we checked that this component is negligible compared to the faint sources
using the MICECAT simulation.. For the bright sources, we define

$\sum_{m_{1.1}<20}L=L_{\rm gal}\cdot N_{\rm eff},$ (28)

where $L_{\rm 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_{\rm eff}$ accounts for the fact that there are
multiple bright galaxies in the halo, and we infer the average $N_{\rm eff}$
value from MICECAT. For our five stacking sub-samples, we get $N_{\rm
eff}\sim$ 2 to 5. From our fitted galaxy profile, we can also calculate
$L_{\rm gal}(>10\,{\rm kpc})$, and we apply the same $N_{\rm eff}$ to model
the extension from other bright galaxies:

$\sum_{m_{1.1}<20}L(>10\,{\rm kpc})=L_{\rm gal}(>10\,{\rm kpc})\cdot N_{\rm
eff}.$ (29)

This results in

$f_{\rm IHL}=\frac{L_{\rm gal}(>10\,{\rm kpc})/L_{\rm gal}}{1+L_{\rm
1h}/\left(N_{\rm eff}\cdot L_{\rm gal}\right)}.$ (30)

We show our constraints on $f_{\rm IHL}$, as a function of halo mass and
redshift in Fig. 14 and 15, respectively. The halo masses associated with our
galaxies are inferred from the MICECAT simulation, and using the SDSS
photometric redshifts.

Note that the fraction of light beyond 10 kpc (the numerator in Eq. 30) is
shown in Fig. 13, where all five sub-samples have similar value at $\sim
50\%$, and therefore the strong redshift dependency of $f_{\rm IHL}$ in Fig.
15 is from the increase of the one-halo term with redshift. We show the ratio
of one-halo term and the stacked galaxy light in Fig. 16. Note that this
observable quantity tracks the evolution of the one-halo luminosity, but lacks
the $N_{\rm eff}$ term in Eq. 30 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.

We compare our results with $f_{\rm 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 Fig. 15).

Figure 14: The IHL fraction $f_{\rm IHL}$ as a function of halo mass. Blue and
red data points show the constraints from this work in the 1.1 $\mu$m and 1.8
$\mu$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, 2008). The ICL fraction in individual galaxy groups and clusters
from Gonzalez et al. (2005, 2007) and Burke et al. (2015) are shown in black
and grey data points. The two downward arrows give upper limits for the Milky
Way (Carollo et al., 2010) and Andromeda (M31) (Courteau et al., 2011). Figure
15:  $f_{\rm IHL}$ constraints as in Fig. 14, but plotted as a function of
redshift. The masses of the Bullock & Johnston (2005) clusters are
100-1000$\times$ the halo masses associated with our galaxies. Figure 16:  The
ratio of the total one-halo term and stacked galaxy profile term from our
stacking results (blue: 1.1 $\mu$m, red: 1.8 $\mu$m) compared with the MICECAT
simulation (light blue: 1.1 $\mu$m, orange: 1.8 $\mu$m). We observe a somewhat
stronger evolution, causing the fall-off of $f_{\rm IHL}$ with redshift seen
in Fig. 15.

###  11.4 Color of the Galaxy Inner and Outer Regions

We calculate the $m_{\rm 1.1}-m_{\rm 1.8}$ color of the inner and outer region
of the galaxy, defined by the total light inside and outside 10 kpc physical
scale in the fitted galaxy profile. The results are summarized in Table 4.
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. 22. We have no detection of a color difference
between the inner and outer regions in the two CIBER bands. 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_{\rm 1.1}-m_{\rm 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 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_{\rm 1.1}-m_{\rm
1.8}$ extension color. The rest-frame I and J band redshifted to $z\sim 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 4, but we note that Zemcov et al. (2014) reflects the integrated IHL
from all redshifts.

Table 4: Constraints on the color ($m_{\rm 1.1}-m_{\rm 1.8}$) of the galaxy inner and outer components. 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. Name | Inner | Outer
---|---|---
high-M/low-z | 0.43${}^{+0.27}_{-0.25}$ | 0.38${}^{+0.19}_{-0.16}$
high-M/med-z | 0.56${}^{+0.33}_{-0.36}$ | 0.48${}^{+0.10}_{-0.12}$
high-M/high-z | 0.65${}^{+0.40}_{-0.38}$ | 0.62${}^{+0.15}_{-0.11}$
low-M/low-z | 0.40${}^{+0.28}_{-0.26}$ | 0.38${}^{+0.24}_{-0.20}$
low-M/med-z | 0.58${}^{+0.33}_{-0.34}$ | 0.47${}^{+0.28}_{-0.24}$
total | 0.50${}^{+0.14}_{-0.14}$ | 0.49${}^{+0.06}_{-0.07}$
MICECAT | 0.44$\pm$0.07 |
Helgason et al. (2012) | 0.41 |
GISSEL |  | 0.32$\pm$ 0.08
NGC 5907 |  | 1.41$\pm$ 0.61
Zemcov et al. (2014) |  | 0.89${}^{+1.17}_{-1.08}$

###  11.5 One-halo and Two-halo Clustering

The one-halo amplitude is detected in the 1.8 $\mu$m band at the $\sim
4\sigma$ level in the “total” and “high-M/high-z” cases, and at the $\sim
3\sigma$ level in “mag bin #4” and “high-M/med-z” cases. One-halo clustering
is not clearly detected at the 1.1 $\mu$m band since the photocurrent from
sources is lower in this band. The one-halo amplitude $A_{\rm 1h}$ is
consistent with unity to within $\sim 2\sigma$, 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 $\mu$m), we detect extended emission in galaxies. The galaxies being
stacked ($\sim 30,000$ galaxies in total) are split into five sub-samples from
SDSS spanning redshifts $0.2\lesssim z\lesssim 0.5$ and stellar masses
$10^{10.5}M_{\odot}\lesssim M_{*}\lesssim 10^{12}M_{\odot}$, comparable to
$L_{*}$ galaxies at this redshift. We jointly fit a model for the inherent
galaxy light profile and large-scale one- and two-halo clustering.

With the galaxy profile, we estimate that $\sim 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 $\mu$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 $\mu$m band at 4-$\sigma$ significance over the full sample of
galaxies. These results suggest non-linear 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
non-linear clustering in IHL modeling.

The intrinsic galaxy profile fitted from our stacking analysis suggests $\sim
50\%$ of the total galaxy light resides in the outskirts of galaxies at $r>10$
kpc. This result is in agreement with previous HSC measurements at lower
redshifts ($0\lesssim z\lesssim 0.25$) and lower stellar masses
($10^{9.2}M_{\odot}<M_{*}<10^{11.4}M_{\odot}$). The galaxy extension accounts
for significant fraction of luminosity in $L_{*}$ galaxies, but falls off
below $M_{*}\sim 10^{11}M_{\odot}$. We measure an increase in $f_{\rm 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 also suggest that the
galaxy extension is significant for $L_{*}$ galaxies at $0.2\lesssim z\lesssim
0.5$. The flux from the extension, and the one-halo term present in our galaxy
samples indicate that IHL is a significant contribution to EBL and EBL
fluctuation measurements.

We thank the dedicated efforts of the sounding rocket staff at NASA Wallops
Flight Facility and White Sands Missile Range. This work was supported by NASA
APRA research grants NNX07AI54G, NNG05WC18G, NNX07AG43G, NNX07AJ24G,
NNX10AE12G, and 80NSSC20K0595. Initial support was provided by an award to
J.B. from the Jet Propulsion Laboratory’s Director’s Research and Development
Fund. Japanese participation in CIBER was supported by KAKENHI (20.34,
18204018, 19540250, 21340047, 21111004, 24111717, 26800112, and 15H05744) from
Japan Society for the Promotion of Science (JSPS) and the Ministry of
Education, Culture, Sports, Science and Technology (MEXT). Korean
participation in CIBER was supported by the Pioneer Project from Korea
Astronomy and Space Science Institute (KASI). Y.-T. C. acknowledges support by
the Ministry of Education, Taiwan through the Taiwan-Caltech Scholarship. C.
H. N. acknowledges support by NASA Headquarters under the NASA Earth and Space
Science Fellowship Program - Grant 80NSSCK0706. This publication makes use of
data products from the Two Micron All Sky Survey, which is a joint project of
the University of Massachusetts and the Infrared Processing and Analysis
Center/California Institute of Technology, funded by the National Aeronautics
and Space Administration and the National Science Foundation. The Pan-STARRS1
Surveys (PS1) and the PS1 public science archive have been made possible
through contributions by the Institute for Astronomy, the University of
Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its
participating institutes, the Max Planck Institute for Astronomy, Heidelberg
and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns
Hopkins University, Durham University, the University of Edinburgh, the
Queen’s University Belfast, the Harvard-Smithsonian Center for Astrophysics,
the Las Cumbres Observatory Global Telescope Network Incorporated, the
National Central University of Taiwan, the Space Telescope Science Institute,
the National Aeronautics and Space Administration under Grant No. NNX08AR22G
issued through the Planetary Science Division of the NASA Science Mission
Directorate, the National Science Foundation Grant No. AST-1238877, the
University of Maryland, Eotvos Lorand University (ELTE), the Los Alamos
National Laboratory, and the Gordon and Betty Moore Foundation. Funding for
the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan
Foundation, the U.S. Department of Energy Office of Science, and the
Participating Institutions. SDSS-IV acknowledges support and resources from
the Center for High-Performance Computing at the University of Utah. The SDSS
web site is www.sdss.org. SDSS-IV is managed by the Astrophysical Research
Consortium for the Participating Institutions of the SDSS Collaboration
including the Brazilian Participation Group, the Carnegie Institution for
Science, Carnegie Mellon University, the Chilean Participation Group, the
French Participation Group, Harvard-Smithsonian Center for Astrophysics,
Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli
Institute for the Physics and Mathematics of the Universe (IPMU) / University
of Tokyo, the Korean Participation Group, Lawrence Berkeley National
Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-
Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik
(MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE),
National Astronomical Observatories of China, New Mexico State University, New
York University, University of Notre Dame, Observatário Nacional / MCTI, The
Ohio State University, Pennsylvania State University, Shanghai Astronomical
Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma
de México, University of Arizona, University of Colorado Boulder, University
of Oxford, University of Portsmouth, University of Utah, University of
Virginia, University of Washington, University of Wisconsin, Vanderbilt
University, and Yale University. This work has made use of data from the
European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia),
processed by the Gaia Data Processing and Analysis Consortium (DPAC,
https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has
been provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement. This work has made use of
CosmoHub. CosmoHub has been developed by the Port d’Informació Científica
(PIC), maintained through a collaboration of the Institut de Física d’Altes
Energies (IFAE) and the Centro de Investigaciones Energéticas,
Medioambientales y Tecnológicas (CIEMAT) and the Institute of Space Sciences
(CSIC & IEEC), and was partially funded by the ”Plan Estatal de Investigación
Científica y Técnica y de Innovación” program of the Spanish government.

## References

  * Abell (1958) Abell, G. O. 1958, ApJS, 3, 211, doi: 10.1086/190036
  * Abraham & van Dokkum (2014) Abraham, R. G., & van Dokkum, P. G. 2014, PASP, 126, 55, doi: 10.1086/674875
  * Aihara et al. (2011) Aihara, H., Allende Prieto, C., An, D., et al. 2011, ApJS, 193, 29, doi: 10.1088/0067-0049/193/2/29
  * Arnouts et al. (1999) Arnouts, S., Cristiani, S., Moscardini, L., et al. 1999, MNRAS, 310, 540, doi: 10.1046/j.1365-8711.1999.02978.x
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Bernstein (2007) Bernstein, R. A. 2007, ApJ, 666, 663, doi: 10.1086/519824
  * Blanton et al. (2017) Blanton, M. R., Bershady, M. A., Abolfathi, B., et al. 2017, AJ, 154, 28, doi: 10.3847/1538-3881/aa7567
  * Bock et al. (2013) Bock, J., Sullivan, I., Arai, T., et al. 2013, ApJS, 207, 32, doi: 10.1088/0067-0049/207/2/32
  * Bruzual A. & Charlot (1993) Bruzual A., G., & Charlot, S. 1993, ApJ, 405, 538, doi: 10.1086/172385
  * Bullock & Johnston (2005) Bullock, J. S., & Johnston, K. V. 2005, ApJ, 635, 931, doi: 10.1086/497422
  * Burke et al. (2015) Burke, C., Hilton, M., & Collins, C. 2015, MNRAS, 449, 2353, doi: 10.1093/mnras/stv450
  * Carollo et al. (2010) Carollo, D., Beers, T. C., Chiba, M., et al. 2010, ApJ, 712, 692, doi: 10.1088/0004-637X/712/1/692
  * Carretero et al. (2015) Carretero, J., Castander, F. J., Gaztañaga, E., Crocce, M., & Fosalba, P. 2015, MNRAS, 447, 646, doi: 10.1093/mnras/stu2402
  * Cavaliere & Fusco-Femiano (1978) Cavaliere, A., & Fusco-Femiano, R. 1978, A&A, 70, 677
  * Chambers et al. (2016) Chambers, K. C., Magnier, E. A., Metcalfe, N., et al. 2016, arXiv e-prints, arXiv:1612.05560.  https://arxiv.org/abs/1612.05560
  * Conroy et al. (2007) Conroy, C., Wechsler, R. H., & Kravtsov, A. V. 2007, ApJ, 668, 826, doi: 10.1086/521425
  * Cooper et al. (2013) Cooper, A. P., D’Souza, R., Kauffmann, G., et al. 2013, MNRAS, 434, 3348, doi: 10.1093/mnras/stt1245
  * Cooper et al. (2015) Cooper, A. P., Parry, O. H., Lowing, B., Cole, S., & Frenk, C. 2015, MNRAS, 454, 3185, doi: 10.1093/mnras/stv2057
  * Cooper et al. (2010) Cooper, A. P., Cole, S., Frenk, C. S., et al. 2010, MNRAS, 406, 744, doi: 10.1111/j.1365-2966.2010.16740.x
  * Cooray (2016) Cooray, A. 2016, Royal Society Open Science, 3, 150555, doi: 10.1098/rsos.150555
  * Cooray & Sheth (2002) Cooray, A., & Sheth, R. 2002, Phys. Rep., 372, 1, doi: 10.1016/S0370-1573(02)00276-4
  * Cooray et al. (2012) Cooray, A., Smidt, J., de Bernardis, F., et al. 2012, Nature, 490, 514, doi: 10.1038/nature11474
  * Courteau et al. (2011) Courteau, S., Widrow, L. M., McDonald, M., et al. 2011, ApJ, 739, 20, doi: 10.1088/0004-637X/739/1/20
  * Crocce et al. (2015) Crocce, M., Castander, F. J., Gaztañaga, E., Fosalba, P., & Carretero, J. 2015, MNRAS, 453, 1513, doi: 10.1093/mnras/stv1708
  * Domínguez et al. (2011) Domínguez, A., Primack, J. R., Rosario, D. J., et al. 2011, MNRAS, 410, 2556, doi: 10.1111/j.1365-2966.2010.17631.x
  * Driver et al. (2016) Driver, S. P., Andrews, S. K., Davies, L. J., et al. 2016, ApJ, 827, 108, doi: 10.3847/0004-637X/827/2/108
  * D’Souza et al. (2014) D’Souza, R., Kauffman, G., Wang, J., & Vegetti, S. 2014, MNRAS, 443, 1433, doi: 10.1093/mnras/stu1194
  * Elias et al. (2018) Elias, L. M., Sales, L. V., Creasey, P., et al. 2018, MNRAS, 479, 4004, doi: 10.1093/mnras/sty1718
  * Fernandez et al. (2010) Fernandez, E. R., Komatsu, E., Iliev, I. T., & Shapiro, P. R. 2010, ApJ, 710, 1089, doi: 10.1088/0004-637X/710/2/1089
  * Foreman-Mackey (2016) Foreman-Mackey, D. 2016, The Journal of Open Source Software, 1, 24, doi: 10.21105/joss.00024
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * Fosalba et al. (2015a) Fosalba, P., Crocce, M., Gaztañaga, E., & Castand er, F. J. 2015a, MNRAS, 448, 2987, doi: 10.1093/mnras/stv138
  * Fosalba et al. (2015b) Fosalba, P., Gaztañaga, E., Castander, F. J., & Crocce, M. 2015b, MNRAS, 447, 1319, doi: 10.1093/mnras/stu2464
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1, doi: 10.1051/0004-6361/201629272
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1, doi: 10.1051/0004-6361/201833051
  * Gelman & Rubin (1992) Gelman, A., & Rubin, D. B. 1992, Statistical Science, 7, 457, doi: 10.1214/ss/1177011136
  * Gonzalez et al. (2005) Gonzalez, A. H., Zabludoff, A. I., & Zaritsky, D. 2005, ApJ, 618, 195, doi: 10.1086/425896
  * Gonzalez et al. (2007) Gonzalez, A. H., Zaritsky, D., & Zabludoff, A. I. 2007, ApJ, 666, 147, doi: 10.1086/519729
  * Hartlap et al. (2007) Hartlap, J., Simon, P., & Schneider, P. 2007, A&A, 464, 399, doi: 10.1051/0004-6361:20066170
  * Helgason et al. (2012) Helgason, K., Ricotti, M., & Kashlinsky, A. 2012, ApJ, 752, 113, doi: 10.1088/0004-637X/752/2/113
  * Hirata & Choi (2020) Hirata, C. M., & Choi, A. 2020, PASP, 132, 014501, doi: 10.1088/1538-3873/ab44f7
  * Hoffmann et al. (2015) Hoffmann, K., Bel, J., Gaztañaga, E., et al. 2015, MNRAS, 447, 1724, doi: 10.1093/mnras/stu2492
  * Huang et al. (2018) Huang, S., Leauthaud, A., Greene, J. E., et al. 2018, MNRAS, 475, 3348, doi: 10.1093/mnras/stx3200
  * Ilbert et al. (2006) Ilbert, O., Arnouts, S., McCracken, H. J., et al. 2006, A&A, 457, 841, doi: 10.1051/0004-6361:20065138
  * Kashlinsky et al. (2012) Kashlinsky, A., Arendt, R. G., Ashby, M. L. N., et al. 2012, ApJ, 753, 63, doi: 10.1088/0004-637X/753/1/63
  * Kashlinsky et al. (2005) Kashlinsky, A., Arendt, R. G., Mather, J., & Moseley, S. H. 2005, Nature, 438, 45, doi: 10.1038/nature04143
  * Kawara et al. (2017) Kawara, K., Matsuoka, Y., Sano, K., et al. 2017, PASJ, 69, 31, doi: 10.1093/pasj/psx003
  * Keenan et al. (2010) Keenan, R. C., Barger, A. J., Cowie, L. L., & Wang, W. H. 2010, ApJ, 723, 40, doi: 10.1088/0004-637X/723/1/40
  * Kim et al. (2019) Kim, M. G., Matsumoto, T., Lee, H. M., et al. 2019, PASJ, 71, 82, doi: 10.1093/pasj/psz063
  * Korngut et al. (2013) Korngut, P. M., Renbarger, T., Arai, T., et al. 2013, ApJS, 207, 34, doi: 10.1088/0067-0049/207/2/34
  * Koushan et al. (2021) Koushan, S., Driver, S. P., Bellstedt, S., et al. 2021, arXiv e-prints, arXiv:2102.12323.  https://arxiv.org/abs/2102.12323
  * Lang et al. (2010) Lang, D., Hogg, D. W., Mierle, K., Blanton, M., & Roweis, S. 2010, AJ, 139, 1782, doi: 10.1088/0004-6256/139/5/1782
  * Lauer et al. (2020) Lauer, T. R., Postman, M., Weaver, H. A., et al. 2020, arXiv e-prints, arXiv:2011.03052.  https://arxiv.org/abs/2011.03052
  * Levenson et al. (2007) Levenson, L. R., Wright, E. L., & Johnson, B. D. 2007, ApJ, 666, 34, doi: 10.1086/520112
  * Lin & Mohr (2004) Lin, Y.-T., & Mohr, J. J. 2004, ApJ, 617, 879, doi: 10.1086/425412
  * Lucy (1974) Lucy, L. B. 1974, AJ, 79, 745, doi: 10.1086/111605
  * Martínez-Delgado et al. (2010) Martínez-Delgado, D., Gabany, R. J., Crawford, K., et al. 2010, AJ, 140, 962, doi: 10.1088/0004-6256/140/4/962
  * Matsumoto et al. (2015) Matsumoto, T., Kim, M. G., Pyo, J., & Tsumura, K. 2015, ApJ, 807, 57, doi: 10.1088/0004-637X/807/1/57
  * Matsumoto & Tsumura (2019) Matsumoto, T., & Tsumura, K. 2019, PASJ, 71, 88, doi: 10.1093/pasj/psz070
  * Matsumoto et al. (2018) Matsumoto, T., Tsumura, K., Matsuoka, Y., & Pyo, J. 2018, AJ, 156, 86, doi: 10.3847/1538-3881/aad0f0
  * Matsumoto et al. (2011) Matsumoto, T., Seo, H. J., Jeong, W. S., et al. 2011, ApJ, 742, 124, doi: 10.1088/0004-637X/742/2/124
  * Matsuura et al. (2017) Matsuura, S., Arai, T., Bock, J. J., et al. 2017, ApJ, 839, 7, doi: 10.3847/1538-4357/aa6843
  * Mitchell-Wynne et al. (2015) Mitchell-Wynne, K., Cooray, A., Gong, Y., et al. 2015, Nature Communications, 6, 7945, doi: 10.1038/ncomms8945
  * Muzzin et al. (2013) Muzzin, A., Marchesini, D., Stefanon, M., et al. 2013, ApJ, 777, 18, doi: 10.1088/0004-637X/777/1/18
  * Petrosian (1976) Petrosian, V. 1976, ApJ, 210, L53, doi: 10.1086/182301
  * Planck Collaboration et al. (2016) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A13, doi: 10.1051/0004-6361/201525830
  * Purcell et al. (2007) Purcell, C. W., Bullock, J. S., & Zentner, A. R. 2007, ApJ, 666, 20, doi: 10.1086/519787
  * Purcell et al. (2008) —. 2008, MNRAS, 391, 550, doi: 10.1111/j.1365-2966.2008.13938.x
  * Richardson (1972) Richardson, W. H. 1972, Journal of the Optical Society of America (1917-1983), 62, 55
  * Rodriguez-Gomez et al. (2016) Rodriguez-Gomez, V., Pillepich, A., Sales, L. V., et al. 2016, MNRAS, 458, 2371, doi: 10.1093/mnras/stw456
  * Rowan-Robinson et al. (2013) Rowan-Robinson, M., Gonzalez-Solares, E., Vaccari, M., & Marchetti, L. 2013, MNRAS, 428, 1958, doi: 10.1093/mnras/sts163
  * Rowan-Robinson et al. (2008) Rowan-Robinson, M., Babbedge, T., Oliver, S., et al. 2008, MNRAS, 386, 697, doi: 10.1111/j.1365-2966.2008.13109.x
  * Rudick et al. (2009) Rudick, C. S., Mihos, J. C., Frey, L. H., & McBride, C. K. 2009, ApJ, 699, 1518, doi: 10.1088/0004-637X/699/2/1518
  * Rudy et al. (1997) Rudy, R. J., Woodward, C. E., Hodge, T., Fairfield, S. W., & Harker, D. E. 1997, Nature, 387, 159, doi: 10.1038/387159a0
  * Saldana-Lopez et al. (2020) Saldana-Lopez, A., Domínguez, A., Pérez-González, P. G., et al. 2020, arXiv e-prints, arXiv:2012.03035.  https://arxiv.org/abs/2012.03035
  * Sano et al. (2015) Sano, K., Kawara, K., Matsuura, S., et al. 2015, ApJ, 811, 77, doi: 10.1088/0004-637X/811/2/77
  * Sano et al. (2020) Sano, K., Matsuura, S., Yomo, K., & Takahashi, A. 2020, ApJ, 901, 112, doi: 10.3847/1538-4357/abad3d
  * Seo et al. (2015) Seo, H. J., Lee, H. M., Matsumoto, T., et al. 2015, ApJ, 807, 140, doi: 10.1088/0004-637X/807/2/140
  * Skrutskie et al. (2006) Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163, doi: 10.1086/498708
  * Symons et al. (2021) Symons, T., Zemcov, M., Bock, J., et al. 2021, ApJS, 252, 24, doi: 10.3847/1538-4365/abcaa5
  * Tal & van Dokkum (2011) Tal, T., & van Dokkum, P. G. 2011, ApJ, 731, 89, doi: 10.1088/0004-637X/731/2/89
  * Tal et al. (2009) Tal, T., van Dokkum, P. G., Nelan, J., & Bezanson, R. 2009, AJ, 138, 1417, doi: 10.1088/0004-6256/138/5/1417
  * Thompson et al. (2007) Thompson, R. I., Eisenstein, D., Fan, X., Rieke, M., & Kennicutt, R. C. 2007, ApJ, 666, 658, doi: 10.1086/520634
  * Tsumura et al. (2013a) Tsumura, K., Matsumoto, T., Matsuura, S., Sakon, I., & Wada, T. 2013a, PASJ, 65, 121, doi: 10.1093/pasj/65.6.121
  * Tsumura et al. (2013b) Tsumura, K., Arai, T., Battle, J., et al. 2013b, ApJS, 207, 33, doi: 10.1088/0067-0049/207/2/33
  * van Dokkum et al. (2014) van Dokkum, P. G., Abraham, R., & Merritt, A. 2014, ApJ, 782, L24, doi: 10.1088/2041-8205/782/2/L24
  * Wang et al. (2019) Wang, W., Han, J., Sonnenfeld, A., et al. 2019, MNRAS, 487, 1580, doi: 10.1093/mnras/stz1339
  * Wen et al. (2012) Wen, Z. L., Han, J. L., & Liu, F. S. 2012, ApJS, 199, 34, doi: 10.1088/0067-0049/199/2/34
  * Zemcov et al. (2017) Zemcov, M., Immel, P., Nguyen, C., et al. 2017, Nature Communications, 8, 15003, doi: 10.1038/ncomms15003
  * Zemcov et al. (2013) Zemcov, M., Arai, T., Battle, J., et al. 2013, ApJS, 207, 31, doi: 10.1088/0067-0049/207/2/31
  * Zemcov et al. (2014) Zemcov, M., Smidt, J., Arai, T., et al. 2014, Science, 346, 732, doi: 10.1126/science.1258168
  * Zhang et al. (2019) Zhang, Y., Yanny, B., Palmese, A., et al. 2019, ApJ, 874, 165, doi: 10.3847/1538-4357/ab0dfd
  * Zibetti et al. (2005) Zibetti, S., White, S. D. M., Schneider, D. P., & Brinkmann, J. 2005, MNRAS, 358, 949, doi: 10.1111/j.1365-2966.2005.08817.x