# Near-infrared Extragalactic Background Light Fluctuations on Nonlinear
Scales

Yun-Ting Cheng<EMAIL_ADDRESS>James J. Bock

###### Abstract

Several fluctuation studies on the near-infrared extragalactic background
light (EBL) find an excess power at tens of arcminute scales ($\ell\sim
10^{3}$). Emission from the intra-halo light (IHL) has been proposed as a
possible explanation for the excess signal. In this work, we investigate the
emission from the integrated galaxy light (IGL) and IHL in the power spectrum
of EBL fluctuations using the simulated galaxy catalog MICECAT. We find that
at $\ell\sim 10^{3}$, the one-halo clustering from satellite galaxies has
comparable power to the two-halo term in the IGL power spectrum. In some
previous EBL analyses, the IGL model assumed a small one-halo clustering
signal, which may result in overestimating the IHL contribution to the EBL. We
also investigate the dependence of the IGL$+$IHL power spectrum on the IHL
distribution as a function of redshift and halo mass, and the spatial profile
within the halo. Our forecast suggests that the upcoming SPHEREx deep field
survey can distinguish different IHL models considered in this work with high
significance. Finally, we quantify the bias in the power spectrum from the
correlation of the mask and the signal, which has not been accounted for in
previous analyses.

cosmology: Diffuse radiation –- Near infrared astronomy – Large-scale
structure of universe – Galaxy evolution – Cosmic background radiation

††software:  astropy (Astropy Collaboration et al., 2013, 2018)

California Institute of Technology, 1200 E. California Boulevard, Pasadena, CA
91125, 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

##  1 Introduction

The extragalactic background light (EBL) arises from the aggregate emission
from all sources across cosmic time, and has been measured across the
electromagnetic spectrum from gamma rays to the radio. The EBL at different
wavelengths originates from different astrophysical sources, and reflects a
wide range of physical processes in the universe (see, e.g., Cooray, 2016;
Mattila & Väisänen, 2019, for recent reviews). At near-infrared wavelengths,
the integrated galaxy light (IGL) is thought to be primarily produced by
redshifted ultraviolet and optical stellar emission from $z\lesssim 2$.
However, several EBL absolute photometry measurements (Tsumura et al., 2013;
Matsumoto et al., 2015; Sano et al., 2015; Matsuura et al., 2017; Zemcov et
al., 2017; Sano et al., 2020; Lauer et al., 2021; Carleton et al., 2022) have
reported excess emission above the IGL (Keenan et al., 2010; Domínguez et al.,
2011; Helgason et al., 2012; Driver et al., 2016; Koushan et al., 2021;
Saldana-Lopez et al., 2021), suggesting potential near-infrared EBL
contributions from faint and diffuse populations. EBL fluctuations as a
function of wavelength and angular scale probe the EBL sources that trace the
underlying matter density field. Near-infrared fluctuation analyses have also
reported excess fluctuations over the expected IGL signal at tens of arcminute
scales (Kashlinsky et al., 2005; Thompson et al., 2007; Matsumoto et al.,
2011; Kashlinsky et al., 2012; Cooray et al., 2012; Zemcov et al., 2014; Seo
et al., 2015; Mitchell-Wynne et al., 2015; Kim et al., 2019; Matsumoto &
Tsumura, 2019).

As the scales of excess fluctuations correspond to galaxy nonlinear clustering
at low redshifts ($z\lesssim 0.5$), some studies suggest that the excess
signal is mostly contributed by stripped stars in dark matter halos, referred
to as intra-halo light (IHL; Cooray et al., 2012; Zemcov et al., 2014; Cheng
et al., 2021), while others explain it with the emission from first stars and
galaxies from the epoch of reionization (EoR; Santos et al., 2002; Salvaterra
& Ferrara, 2003; Kashlinsky et al., 2005; Matsumoto et al., 2011; Kashlinsky
et al., 2012; Mitchell-Wynne et al., 2015; Sun et al., 2021). More exotic
candidates responsible for the excess EBL have also been proposed, such as
Population III stars and accretion disks around their stellar mass black holes
from $z>7$ (Windhorst et al., 2018), direct collapse black holes from the dark
ages (Yue et al., 2013, 2014) and the decay of axion-like particles as dark
matter candidates (e.g., Gong et al., 2016; Kalashev et al., 2019; Caputo et
al., 2021).

In previous fluctuation analyses (e.g. Kashlinsky et al., 2012; Cooray et al.,
2012; Zemcov et al., 2014), the nonlinear clustering from satellite galaxies
is based on the analytical halo occupation distribution (HOD) prescription
from Helgason et al. (2012, hereafter H12). However, this model adopts the
shape of the IGL power spectrum from galaxy number counts, and assumes the
galaxy clustering is independent of luminosity, which results in an
underestimate of the nonlinear clustering power. Fernandez et al. (2010)
demonstrate that nonlinear bias from halo emission at the EoR has a
significant effect on the EBL power spectrum, which cannot be ignored.
Although the galaxy bias is much smaller at lower redshifts than during the
EoR, it is still important to investigate the effects of nonlinear clustering
from satellite galaxies and the IHL. Therefore, a more realistic modeling at
nonlinear clustering scales is essential for interpreting the EBL fluctuation
data, particularly for upcoming experiments that will observe EBL fluctuations
with high sensitivity and wide spectral coverage such as SPHEREx (Doré et al.,
2014, 2018), Euclid (Amendola et al., 2018) and the Roman Space Telescope
(Spergel et al., 2015). For example, SPHEREx deep fields will map the EBL in
200 deg2 across 0.75–5 $\mu$m in 102 spectral channels, promising highly
accurate EBL spatial and spectral measurements that will improve constraints
on different components in the EBL signal.

In this work, we use the simulated galaxy catalog MICECAT (Sec. 2.1) to model
the EBL power spectrum from IGL and IHL at low-redshift ($z\leqslant 1.4$),
which captures the bulk of the near-infrared IGL intensity and fluctuations.
MICECAT contains information on individual halos and their central and
satellite galaxies, quantifying the nonlinear clustering from satellite
galaxies, which is ideal for studying the effects on the power spectrum with
different IHL models. We investigate the EBL power spectrum with IHL models
spanning a wide range of redshift, halo mass, and spatial profile. This allows
us to understand the dependence of EBL fluctuations on IHL model parameters
and constraints on IHL models with future observations.

This paper is organized as follows. Sec. 2 introduces MICECAT and our
application of the catalog to our calculations. Sec. 3 presents the MICECAT
IGL model and its power spectrum. In Sec. 4, we consider several different IHL
models, and analyze their power spectra as a function of the parameters of the
model. We summarize our results and discuss contexture outlook and future
measurements in Sec. 5. 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
Planck cosmology (Planck Collaboration et al., 2016). All fluxes are quoted in
the AB magnitude system.

##  2 Simulations

Our EBL model builds on an $N$-body simulation catalog – MICECAT (detailed in
Sec. 2.1). As it models the three-dimensional distribution of dark matter
halos and galaxies, MICECAT contains IGL clustering information on both large
(linear) and small (nonlinear) spatial scales. The simulation catalog provides
a flexible basis to implement different IGL and IHL models, as well as the
effects from observation or data processing such as the instrument point-
spread function (PSF) and source masking.

###  2.1 MICECAT

We use the MICECAT (v2.0) simulated galaxy catalog (Fosalba et al., 2015a, b;
Hoffmann et al., 2015)111For details on MICECAT v2.0, see the documentation:
https://www.dropbox.com/s/0ffa8e7463n8h1q/README_MICECAT_v2.0_for_new_CosmoHub.pdf?dl=0
to model the EBL fluctuations from the IGL and IHL. MICECAT is a product of
the $N$-body cosmological simulation Marenostrum Institut de Ciéncies de
l’Espai 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 simulates ideal observations of one octant of the sky ($\sim$ 5000
deg2; $0^{\circ}<{\rm R.A.}<90^{\circ}$, $0^{\circ}<{\rm decl.}<90^{\circ}$),
and covers the redshift range of $0<z<1.4$ without repeating the simulation
box. MICECAT builds on MICE-GC by combining an HOD with subhalo abundance
matching to calibrate to observed luminosity functions and clustering
(Carretero et al., 2015). MICECAT simulates a mass-limited sample, with deeper
coverage in two sub-regions (decl.$>30^{\circ}$; decl.$<30^{\circ}$ and
$30^{\circ}<{\rm decl.}<60^{\circ}$), which is complete to $m_{i}\sim 24$ at
all redshifts ($0<z<1.4$). In this work, we only use their simulation in these
two deep subregions.

For each galaxy, MICECAT provides its magnitude in several optical and near-
infrared bands (SDSS, DES, Euclid, etc), which allows us to model the IGL
wavelength dependence. All galaxies in MICECAT have a halo ID and a
central/satellite galaxy flag, so we can identify the central galaxy and all
satellite galaxies in each halo and use them to model the linear (two-halo)
and nonlinear (one-halo) clustering. Yang et al. (2009) develop an algorithm
to group observed galaxies in the Sloan Digital Sky Survey (SDSS) into halos
and identify central and satellite galaxies for each halo. As a validation, we
compare the conditional stellar mass functions for both central and satellite
galaxies in MICECAT and the group catalog from Yang et al. (2009), and find
that the two catalogs are consistent across halo masses of $10^{12}\lesssim
M\lesssim 10^{14}\,M_{\odot}$.

###  2.2 Modeling the EBL with MICECAT

In this work, we only model the EBL in a single observing band, and leave
wavelength dependence studies to future work. We select the Euclid NISP Y band
from MICECAT, which has an effective wavelength of $\lambda_{\rm eff}=1.07$
$\mu$m. All simulations throughout this work are run on $2\times 2$ deg2
fields selected from the two deep subregions in MICECAT. For each field we use
$1024\times 1024$ pixels, with a $7{{}^{\prime\prime}}\times
7{{}^{\prime\prime}}$ pixel size, matching CIBER imager data (Bock et al.,
2013; Zemcov et al., 2014; Cheng et al., 2021). Our choice of field and pixel
size corresponds to multipoles $10^{2}\lesssim\ell\lesssim 10^{5}$, suitable
for modeling the nonlinear clustering signal that dominates at tens of
arcminute scales ($\ell\sim 10^{3}$). We use 35 such fields (covering 140 deg2
in total) to obtain statistics with a negligible sample variance at the scales
of interest. When forecasting sensitivity, we also consider a 200 deg2 survey,
the total area of the SPHEREx deep fields near the north and south ecliptic
poles.

We find an excess number of low-redshift faint galaxies above a Schechter
luminosity function (Schechter, 1976) in MICECAT. It is unclear whether these
faint sources are consistent with observations or are artifacts from the
simulation. H12 shows that galaxy number counts are consistent with the
Schechter function fit to the faint end ($m_{\rm AB}\sim 26$ at $\sim 1$
$\mu$m), while in Beare et al. (2019) (Fig. 11), the rest-frame $K$ band
luminosity function at $z\sim 3$ from observation shows a faint-end excess
above the Schechter function fit. These sources may be excess and would
potentially bias the IGL clustering signal, especially on nonlinear scales as
many of them are satellite galaxies in dark matter halos. Furthermore, in
reality, faint galaxies are more likely to experience tidal disruptions during
the course of their evolution history, which will then be described by the IHL
component in our model. Therefore, to make a conservative estimate of the IGL,
we discard all $M_{r}>-18$ MICECAT galaxies, where $M_{r}$ is the SDSS r-band
absolute magnitude given by MICECAT. Fig. 1 shows the MICECAT luminosity
function at $z=0.2$ in the Euclid NISP Y band, and a model of IGL luminosity
function from H12, fitted to a Schechter function from observed galaxy counts.
We also determined that in the same band, the $M_{r}>-18$ population only
contributes $\sim 10\%$ of the total IGL in MICECAT at $z\lesssim 0.5$; there
are almost no $M_{r}>-18$ sources at $z\gtrsim 1$ (Fig. 2). In addition, we
also checked that discarding $M_{r}>-18$ sources only reduces $\sim 10\%$ of
the IGL power spectrum at $\ell\sim 10^{3}$. Note that galaxy tidal disruption
is a continuous process with no clear boundary between small satellite
galaxies and IHL, so we incorporate all $M_{r}>-18$ galaxies into our model of
the IHL component.

Figure 1:  MICECAT luminosity function at $z=0.2$ in the Euclid NISP Y band
(black), and the $M_{r}<-18$ population (red) used in this work. Error bars
give the sample variance of the 140 deg2 of MICECAT data we used in this
study. The blue line shows a model of luminosity function at the same
wavelength and redshift from H12, fitting a Schechter function to observed
galaxy counts. The MICECAT luminosity function on the faint end shows an
excess above the model, and thus we discard all $M_{r}>-18$ sources in our
analysis. Figure 2:  Fraction of the MICECAT IGL from $M_{r}<-18$ galaxies as
a function of redshift. The galaxies fainter than $M_{r}=-18$ only constitute
$\sim 10\%$ of the IGL at $z\lesssim 0.5$, and this fraction decreases to zero
at $z\gtrsim 1$. Figure 3:  Top: fraction of IGL intensity from $z<1.4$
galaxies in H12 model as a function of masking threshold $m_{\rm th}$ at $1$
(blue), $3$ (green), and $5$ (red) $\mu$m. The masking threshold is defined at
1 $\mu$m for all three wavelengths, and the corresponding masking depth at
other wavelengths is derived by the abundance matching technique (Cheng &
Chang, 2022). Bottom: fraction of IGL power spectrum from linear clustering
and Poisson fluctuations at $\ell=10^{3}$ from $z<1.4$ galaxies. We use the
power spectrum model from Cheng & Chang (2022) built from the H12 IGL model.

We use the H12 IGL model to evaluate the EBL signal from higher redshift
galaxies that are not included in MICECAT ($z>1.4$). The top panel of Fig. 3
shows the fraction of IGL intensity from $z<1.4$ galaxies at three near-
infrared wavelengths ($1$, $3$, and $5$ $\mu$m) as a function of flux
threshold in AB magnitude $m_{\rm th}$. We consider this fraction with a range
of $m_{\rm th}$ at 1 $\mu$m, and infer the masking depth at other wavelengths
with the abundance matching recipe from Cheng & Chang (2022). We also compare
the IGL two-halo and Poisson power spectrum at $\ell=10^{3}$ from $z<1.4$
galaxies with the same power spectrum from the total IGL, shown in the bottom
panel of222In Fig. 3, the ratio at $3$ $\mu$m drops faster than the other two
wavelengths is because in the H12 model, the $L^{*}$, characteristic
luminosity in the Schechter function, has a maximum at rest-frame $\sim 1.5$
$\mu$m, which makes the emission around this wavelength more sensitive to
masking. At $z<1.4$, $3$ $\mu$m band emission is dominated by sources around
this wavelength, and thus its intensity reduces more by masking than 1 and 5
$\mu$m. Fig. 3. We use the power spectrum model from Cheng & Chang (2022),
which is also based on the H12 IGL model. We set $m_{\rm th}=20$ as our
fiducial case (see Sec. 3.2), and thus we estimate that at $\sim 1$ $\mu$m,
the MICECAT catalog captures $\sim 70\%$ and $>95\%$ of the total IGL
intensity and power spectrum, respectively. With the cases considered in this
work, the higher redshift sources have negligible clustering power at the
scales of interest. However, at longer wavelengths or with a deeper masking
depth, the IGL will be weighted toward higher redshifts, so MICECAT will not
be sufficient to model their signal.

##  3 IGL Fluctuations

We first quantify IGL fluctuations from MICECAT galaxies as point sources, and
discuss the effect of PSF broadening in Sec. 4.4. We decompose the linear and
nonlinear clustering signals in MICECAT using a halo model framework (Cooray &
Sheth, 2002, Sec. 3.1), and consider masking bright galaxies to increase
sensitivity on clustering fluctuations (Sec. 3.2).

###  3.1 Halo Model Framework

The power spectrum of the MICECAT map consists of clustering and Poisson
terms:

$C_{\ell}=C_{\ell,{\rm clus}}+C_{\ell,{\rm P}}.$ (1)

The clustering power spectrum captures projected large-scale structure
fluctuations. Using the Limber approximation,

$C_{\ell,{\rm clus}}=\int dz\frac{H(z)}{c\,\chi^{2}(z)}P_{\rm clus}(k,z),$ (2)

where $H$, $\chi$, $c$, and $P_{\rm clus}$ are the Hubble parameter, co-moving
distance, speed of light, and 3D clustering power spectrum, respectively. The
Poisson term arises from the discreteness of galaxies:

$C_{\ell,{\rm P}}=\frac{1}{\Omega_{\rm sur}}\sum_{i}\left(\nu
F_{\nu,i}\right)^{2},$ (3)

where $\Omega_{\rm sur}$ is the survey angular size, $F_{\nu,i}$ is the
specific flux of the $i$th galaxy, and we sum over all galaxies in the map.

In the halo model framework, $P_{\rm clus}$ can be decomposed into two-halo
(linear) and one-halo (nonlinear) terms, $P_{\rm 2h}$ and $P_{\rm 1h}$, which
describe the correlations between different halos and within the same halo,
respectively:

$P_{\rm clus}(k,z)=P_{\rm 1h}(k,z)+P_{\rm 2h}(k,z),$ (4)

and we can also express the angular power spectrum in terms of these two
components:

$C_{\ell,{\rm clus}}=C_{\ell,{\rm 1h}}+C_{\ell,{\rm 2h}},$ (5)

where

$C_{\ell,{\rm 1h}}=\int dz\frac{H(z)}{c\,\chi^{2}(z)}P_{\rm 1h}(k,z),$ (6)

and

$C_{\ell,{\rm 2h}}=\int dz\frac{H(z)}{c\,\chi^{2}(z)}P_{\rm 2h}(k,z).$ (7)

The one- and two-halo 3D power spectrum can be expressed as333We drop the
redshift symbol $z$ hereafter for clarity.

$P_{\rm 1h}(k)=\frac{1}{\bar{n}_{g}^{2}}\int
dM\frac{dn}{dM}I^{2}_{h}(M)u^{2}(k|M),$ (8)

and

$P_{\rm 2h}(k)=P_{\rm lin}(k)\left[\frac{1}{\bar{n}_{g}}\int
dM\frac{dn}{dM}I_{h}(M)b_{h}(M)u(k|M)\right]^{2},$ (9)

where $\bar{n}_{g}$ is the mean number density of galaxies, $\frac{dn}{dM}$ is
the halo mass function, $I_{h}(M)$ and $b_{h}(M)$ are the mean intensity and
bias from all sources within a mass $M$ halo, respectively, $P_{\rm lin}$ is
the linear matter power spectrum, and $u(k|M)$ is the Fourier transform of the
normalized intensity radial profile of mass $M$, which approaches $1$ on large
scales (small $k$), and rolls off at the scale of the halo. We directly
evaluate $C_{\ell,{\rm 1h}}$ and $C_{\ell,{\rm 2h}}$ from the MICECAT catalog
instead of performing the integration in Eq. 8 and 9, described below.

As each galaxy in MICECAT has a halo ID, we can calculate the one- and two-
halo correlation functions from all galaxy pairs within the same halo and
across different halos, respectively. The power spectra are the Fourier
transform of the two correlation functions. Calculating the correlation
function is computationally expensive, as it requires $O(N^{2})$ operations
for $N$ galaxies. Therefore, we use an approximate estimator to derive one-
and two-halo power spectra directly from 2D maps as follows. First, we make an
intensity map of all the central and satellite galaxies in MICECAT (hereafter,
referred to as the full map). The power spectrum of the full map gives the
total clustering signal $C_{\ell,{\rm clus}}$ that contains both one- and two-
halo terms. Next, we make another map by moving all satellite galaxies to the
same location as their central galaxies (hereafter, referred to as the central
map). We define the power spectrum of the central map as our two-halo power
spectrum estimator, $\hat{C}_{\ell,{\rm 2h}}$. The one-halo power spectrum
estimator is defined as $\hat{C}_{\ell,{\rm 1h}}=C_{\ell,{\rm
clus}}-\hat{C}_{\ell,{\rm 2h}}$.

Our central map does not contain one-halo correlations, and its halo profile
is a Dirac delta function as all fluxes are from the halo center, and thus
$u(k|M)=1$ for all $k$. Therefore, $\hat{C}_{\ell,{\rm 1h}}$ and
$\hat{C}_{\ell,{\rm 2h}}$ is given by the integration of $\hat{P}_{\rm 1h}$
and $\hat{P}_{\rm 2h}$ with Eq. 6 and 7, respectively, where

$\hat{P}_{\rm 2h}(k)=P_{\rm lin}(k)\left[\frac{1}{\bar{n}_{g}}\int
dM\frac{dn}{dM}I_{h}(M)b_{h}(M)\right]^{2},$ (10)

and

$\hat{P}_{\rm 1h}(k)=P_{\rm clus}(k)-\hat{P}_{\rm 2h}(k).$ (11)

We calculate $C_{\ell,{\rm clus}}$ and $\hat{C}_{\ell,{\rm 2h}}$ from the
power spectra of the full and central-only map, respectively. From Eq. 9 and
10, we note that on large scales where $u(k|M)\sim 1$, $\hat{P}_{\rm
2h}(k)\approx P_{\rm 2h}(k)$, our estimator converges to the exact expression.
On small scales, $\hat{P}_{\rm 2h}(k)>P_{\rm 2h}(k)$, and thus $\hat{P}_{\rm
1h}(k)<P_{\rm 1h}(k)$, so our one/two-halo power spectrum estimator is a
lower/upper bound of the true power spectrum, respectively.

The scales at which $\hat{C}_{\ell,{\rm 2h}}$ deviates from $C_{\ell,{\rm
2h}}$ are the angular size of halos. We estimate this scale by stacking halos
in MICECAT. With our fiducial masking threshold $m_{\rm th}=20$ (Sec. 3.2),
the stacked halo profile from satellite galaxies has an FWHM $\sim
15{{}^{\prime\prime}}$ ($\ell\sim 4\times 10^{4}$, black dashed line in Fig.
4). Therefore, at the scales of interest (tens of arcminute scales), the error
from our one- and two-halo power spectrum approximation assuming $u(k|M)=1$ is
negligible. As an additional check, we randomly select $10\%$ of the galaxies
in a $2\times 2$ deg2 MICECAT field, and compute the two-point correlation
function of the two-halo term in the full map and the central map. The results
of the two maps are consistent at $\gtrsim 100$ arcsec ($\ell\lesssim 6500$),
and thus the error from our approximated two-halo term is negligible on the
scales of interest. We estimate the Poisson term from the high-$\ell$ power
spectrum, where the clustering signal is negligible, and we subtract Poisson
fluctuations from the total power spectrum to get the clustering terms.

###  3.2 Masking

We consider masking bright galaxies before measuring the power spectrum. With
simulations, we can perform a “perfect masking” by removing bright source
fluxes from the IGL map. In reality, masking is usually implemented by
removing pixels near bright sources. Pixel masking leaves residual flux
outside the masking region that depends on the PSF, and will also cause mode
coupling in the power spectrum. We set the fiducial masking scheme to perfect
masking with an AB magnitude threshold $m_{\rm th}=20$ in the Euclid NISP Y
band. All the results shown in this paper are using this fiducial masking
scheme if not otherwise specified. The effect of the PSF and pixel masking
will be discussed in Sec. 4.4.

###  3.3 IGL Power Spectrum

Figure 4: The IGL power spectrum from MICECAT with no masking (top) and with a
“perfect” masking threshold $m_{\rm th}=18$ (middle) and $20$ (bottom). Gray
data points give the power spectrum of the total IGL. Brown data points show
the two-halo power spectrum estimator from the central map. Error bars in all
cases give the projected uncertainty in a 200 deg2 observation. The green and
orange dashed lines give Poisson fluctuations in the full and central maps,
respectively, which are fitted at high $\ell$. Green and orange solid lines
give the clustering power spectrum of the full and central maps, respectively,
defined by the difference between the total power spectrum and the Poisson
term. The blue line gives the one-halo clustering estimator, which is the
difference between the total clustering signal (green line) and the two-halo
clustering estimator (orange line). The black dashed line is the angular scale
of the average halo profile from stacking, which also corresponds to the
turnover $\ell$ mode for the one-halo power spectrum. Our approximation of
one- and two-halo power spectra ($\hat{C}_{\ell,{\rm 1h}}$ and
$\hat{C}_{\ell,{\rm 2h}}$) is valid for scales larger than this dashed line.

Fig. 4 presents the components of the IGL power spectrum from MICECAT with no
masking, $m_{\rm th}=20$ masking (fiducial case), and $m_{\rm th}=18$ masking
(similar to the masking depth adopted in Zemcov et al. (2014)). The one- and
two-halo components are derived from the full and central maps (see Sec. 3.1).

With deeper masking, the clustering signal is dominated by smaller halos, and
thus the crossover of one- and two-halo power spectrum occurs at larger $\ell$
modes. We note that the crossover scale also depends on the redshift and halo
mass distribution of $f_{\rm IHL}$. If there is more IHL emission from more
massive or lower redshift halos, the average angular size of IHL will be
larger and thus the crossover will be at smaller $\ell$ modes. Fig. 5 compares
the MICECAT power spectrum in redshift slices $0.4<z<0.6$ and $1.2<z<1.4$ at
$m_{\rm th}=20$, where we observe the crossover scale of one- and two-halo
terms is larger at lower redshift.

Figure 5:  The one-(blue) and two-halo (orange) clustering power spectrum from
MICECAT with sources at $0.4<z<0.6$ (top) and $1.2<z<1.4$ (bottom) with
masking threshold $m_{\rm th}=20$. The crossover angular scale of the two
terms (black dashed lines) is smaller at higher redshifts.

In both cases with masking in Fig. 4, we see that one- and two-halo terms have
comparable power at $\sim 10{{}^{\prime}}$ scales ($\ell\sim 10^{3}$). This
indicates that when modeling EBL fluctuations at this scale, the non-linear
clustering from IGL needs to be properly handled. In some previous fluctuation
analyses (e.g., Kashlinsky et al., 2012; Cooray et al., 2012; Zemcov et al.,
2014), the IGL model is adopted from H12, which uses the HOD model fitted to
the SDSS galaxy number counts. This model assumes galaxy clustering is
independent of galaxy luminosity, which results in a fixed shape in the IGL
(one- and two-halo) power spectrum, where different wavelengths and masking
limits only affect the overall amplitude. Furthermore, for the power spectrum
of the IGL intensity field studied here, the contribution from each galaxy is
weighted by its flux, whereas in the case of SDSS galaxy clustering, all
galaxies have equal weights in the power spectrum. These simplifications in
H12 lead to a small one-halo term that equals the two-halo power at $\sim
4^{\prime}$ ($\ell\sim 2500$). For example, as shown in Fig. 4, at the masking
threshold used in the analysis of Zemcov et al. (2014, $m_{\rm th}=18$), our
model predicts a larger contribution from the one-halo term (the one- and two-
halo power spectrum crossover at $\ell\sim 800$), implying that the nonlinear
clustering from satellite galaxies in this previous analysis is
underestimated, and thus their excess power has to be explained by other
sources such as the IHL.

##  4 Intra-Halo Light

Next, we add IHL to our model. We first define the IHL fraction, $f_{\rm
IHL}$, as the fractional added luminosity from IHL to the halo. In the near
infrared, we consider only IGL and IHL contributions to the total halo
emission, so

$f_{\rm IHL}=\frac{F_{\nu,{\rm IHL}}}{F_{\nu,{\rm IGL}}+F_{\nu,{\rm IHL}}}.$
(12)

In reality, $f_{\rm IHL}$ will depend on redshift, halo mass, and the time
evolution of the halo. Different models give different IHL properties in each
halo, resulting in a different scale-dependent contribution to the one- and
two-halo IGL+IHL power spectra. For simplicity, we set our fiducial IHL model
to the case where all halos have the same $f_{\rm IHL}$ regardless of their
redshift and mass, and assume the IHL spatial distribution follows a
Navarro–Frenk–White (NFW) profile (Navarro et al., 1996). In the following
subsections, we apply different redshift, mass, and spatial profile dependence
to the fiducial model, and discuss their impact on the power spectrum.

We first investigate the dependence of the clustering power spectrum on
$f_{\rm IHL}$ values in the fiducial case (Fig. 6). On scales larger than the
typical halo size ($\ell\sim 5\times 10^{3}$), both IGL (one-halo $+$ two-
halo) and IHL trace the matter density field, so the IGL$+$IHL power spectra
converge to the pure two-halo IGL power spectrum (as the one-halo term is
suppressed on large scales) scaled by the square of the total intensity
$\frac{F_{\nu,{\rm IGL}}+F_{\nu,{\rm IHL}}}{F_{\nu,{\rm
IGL}}}=(\frac{1}{1-f_{\rm IHL}})^{2}$. Both satellite galaxies (one-halo term)
and the IHL produce excess power at the halo scale. With our fiducial IHL
model, the IHL gives a slightly higher excess signal compared to the satellite
galaxies (see bottom panel of Fig. 6), where their difference in the power
spectrum is $\sim 10\%$ if $f_{\rm IHL}=0.5$, and less than $5\%$ with two
other lower $f_{\rm IHL}$ cases ($f_{\rm IHL}=0.1,0.2$) at the halo scale. Our
model indicates that in a 200 deg2 survey limited by sample variance, this
difference in the power spectrum can be detected with high significance in all
three $f_{\rm IHL}$ values considered here.

Figure 6:  Top: clustering power spectrum of IGL$+$IHL (one-halo $+$ two-halo)
with IHL fractions $f_{\rm IHL}=0$ (black), $0.1$ (blue), $0.2$ (orange), and
$0.5$ (green). Here, we show the fiducial case where all halos have the same
$f_{\rm IHL}$, and the IHL follows an NFW profile. Error bars denote sample
variance for a 200 deg2 survey. Middle: ratio of the IGL$+$IHL power spectrum
to the two-halo IGL power spectrum, normalized by $\left(\frac{1}{1-f_{\rm
IHL}}\right)^{2}$. We only show the low-$\ell$ modes where there is sufficient
sensitivity on the two-halo power spectrum. The one-halo power in the IHL
model closely follows the one-halo power in the MICECAT simulation with small
departures at $\ell\sim 5000$ indicated below. Bottom: ratio of the IGL$+$IHL
and IGL (one-halo $+$ two-halo) power spectra of the same three $f_{\rm IHL}$
cases. The normalization factor $(\frac{1}{1-f_{\rm IHL}})^{2}$ accounts for
the change in the two-halo amplitude with $f_{\rm IHL}$.

Zemcov et al. (2014) measured fluctuations with the CIBER instrument,
reporting excess fluctuation power over the IGL on angular scales
$500<\ell<2000$. However, as seen in the middle panel of Fig. 6, both MICECAT
and the IHL model predict a significant one-halo contribution to the total
power spectrum, increasing the total power by factors between 1.2 and 2.2 on
these scales. Accounting for the one-halo contribution thus reduces the
inferred IHL luminosity.

We tried fitting the CIBER 1.1 $\mu$m power spectrum measurements to our
fiducial IGL$+$IHL model, but our constant $f_{\rm IHL}$ model does not
accurately fit the shape of the measured power spectrum. The overall amplitude
is set by points at $\ell>10^{4}$, which have high statistical weight, and the
fixed shape then underpredicts the measured power at $\ell<3000$. This
suggests that our fiducial model assuming a constant $f_{\rm IHL}$ for all
halos is not a good description to the data, and thus a more sophisticated
model that includes redshift and halo mass dependence or other free parameters
is essential for building a realistic IHL model for fluctuation measurements,
which we leave to future work.

###  4.1 Redshift Dependent $f_{\rm IHL}$

W apply a redshift dependence to our $f_{\rm IHL}$ model while using the same
fixed mass dependence and the NFW profile as the fiducial model. Therefore, in
our model, the majority of the IHL emission is associated with the
characteristic sizes of halos ($\sim 10^{12}\,M_{\odot}$ at $z=0$) that
dominate the IGL emission. Although it is difficult to observe IHL in typical-
sized halos individually due to their low surface brightness, some model-
dependent constraints have been placed from EBL fluctuation measurements
(Cooray et al., 2012; Zemcov et al., 2014). Fig. 7 shows the $f_{\rm IHL}$–$z$
relation from Cooray et al. (2012), where their model is fitted to 3.6 and 4.5
$\mu$m fluctuation maps from Spitzer. In addition, the same model in two near-
infrared bands in CIBER (1.1 and 1.6 $\mu$m) is also shown in Fig. 7. Zemcov
et al. (2014) study the near-infrared fluctuations with CIBER, and constrain
the IHL using a model similar that of Cooray et al. (2012). For comparison,
Fig. 7 also shows the stellar mass fraction of the intra-cluster light (ICL)
from Contini et al. (2018) (their “STANDARD” model). We see that in Cooray et
al. (2012) model, the $f_{\rm IHL}$ increases with redshift for $z<1.5$, and
the redshift dependence is stronger for longer wavelengths; whereas Contini et
al. (2018) suggest that for massive galaxy clusters, $f_{\rm ICL}$ decreases
with redshift. While there are very limited observational constraints on the
IHL in a typical-sized halo, and the dependence of $f_{\rm IHL}$ on redshift
and halo mass is highly degenerate in EBL fluctuation measurements, more
sensitive observations on individual halos or numerical simulations are
required to explain the $f_{\rm IHL}$ color dependence and the opposite
redshift trend in the IHL and the IGL. As a demonstration in this work, we
adopt the $f_{\rm ICL}$ from Contini et al. (2018) for all halos at a given
redshift as our model for the $f_{\rm IHL}$–$z$ relation.

Cross correlating the 2D EBL images with galaxies at known redshift (Cheng &
Chang, 2022) enables redshift tomography, since EBL emission only correlates
with external tracers from the same redshift. In Fig. 7 we also show a
forecast of cross correlation with the SPHEREx galaxy catalog, with a
projected number density given by their public
products444https://github.com/SPHEREx/Public-
products/blob/master/galaxy_density_v28_base_cbe.txt, using photometric
redshift errors of $0.01<\sigma_{z}/(1+z)<0.03$. To implement cross
correlation with MICECAT, we consider redshift slices with width $\Delta
z=0.03\times(1+z)$, and generate a tracer galaxy catalog by choosing the same
number of brightest galaxies as the SPHEREx number density in the Euclid $Y$
band. The cross power spectrum between the IGL$+$IHL intensity map and the
tracer catalog at $z=0.2$, $0.5$, and $1$ shown in the bottom panel of Fig. 7
readily detects the redshift dependence of $f_{\rm IHL}$ as well as the subtle
effect of the NFW profile in the IHL model compared with one-halo clustering.

Figure 7:  Top: our model of the $f_{\rm IHL}$–$z$ relation, adopting the
empirical ICL fraction ($f_{\rm ICL}$) from Contini et al. (2018) (black). The
IHL model from Cooray et al. (2012) in four near-infrared bands is also shown
for comparison (colored lines). The three colored dots correspond to the
redshift of the power spectra shown in the bottom panel. Bottom: ratio of
galaxy-IGL$+$IHL and galaxy-IGL cross power spectra at $z=0.2$ (blue), $0.5$
(green), and $1$ (red). Error bars denote sample variance in a 200 deg2
survey. The dashed lines are the ratio of IGL$+$IHL and IGL flux at each
redshift, the expected ratio $\frac{1}{1-f_{\rm IHL}}$ on large scales.

###  4.2 Mass Dependent $f_{\rm IHL}$

Next, we investigate the effect of a halo mass dependence on $f_{\rm IHL}$.
Fig. 8 shows the three $f_{\rm IHL}$–$M$ relations considered in this work.
The first model is the fiducial case where $f_{\rm IHL}$ is independent of
halo mass. In the second and the third model, we use a power law with slope
$\alpha=2.5$ and $0.1$ (i.e., $f_{\rm IHL}(M)\propto M^{2.5}$ and $f_{\rm
IHL}(M)\propto M^{0.1}$), respectively, for $M<10^{12.5}\,M_{\odot}$, and
extending the value at $10^{12.5}\,M_{\odot}$ to all $M>10^{12.5}\,M_{\odot}$
halos. The $\alpha=2.5$ model approximates the $f_{\rm IHL}(M)$ function from
an analytical model at $z=0$ (Purcell et al., 2007, 2008), whereas
$\alpha=0.1$ is the best-fit power-law slope found in Zemcov et al. (2014),
who fit the $f_{\rm IHL}$–$M$ relation from the near-infrared fluctuations in
CIBER images. To compare the three cases, we normalize the $f_{\rm IHL}(M)$
function to the same mean IHL fraction, $\bar{f}_{\rm IHL}=0.2$, where
$\bar{f}_{\rm IHL}$ is defined as the fraction of total EBL flux from the IHL,
and we only consider IGL and IHL as the source of total EBL emission.

Figure 8:  Three models of the $f_{\rm IHL}$–$M$ relations considered in this
work. The constant model (blue) is our fiducial case where all the halos have
the same IHL fraction. The other two models are a power law with slope
$\alpha=2.5$ (orange) and $0.1$ (green), respectively, for
$M<10^{12.5}\,M_{\odot}$, and extending the value at $10^{12.5}\,M_{\odot}$
(black dashed line) to all $M>10^{12.5}\,M_{\odot}$ halos. The three curves
are normalized such that the mean IHL fraction $\bar{f}_{\rm IHL}=0.2$, where
$\bar{f}_{\rm IHL}$ is defined as the fraction of total EBL flux from the IHL.

We calculate the ratio of the galaxy-IGL$+$IHL and the galaxy-IGL cross power
spectra for the three $f_{\rm IHL}(M)$ models at $z=0.2$, $0.5$, and $1$, as
shown in Fig. 9. We consider the same redshift slices and the same galaxy
catalog as Sec. 4.1. The ratio is scaled by $\frac{1}{1-\bar{f}_{\rm IHL}}$,
which is the expected scaling of two-halo clustering on large scales. Our
model suggests that with the same $\bar{f}_{\rm IHL}$ value, the mass
dependence of $f_{\rm IHL}$ has a significant impact on fluctuations at
cluster angular scales. In the $\alpha=2.5$ case, the majority of the IHL
emission is from massive halos, which results in large fluctuation power on
one-halo scales compared to the other two cases.

Figure 9:  Ratio of the galaxy-IGL$+$IHL and the galaxy-IGL cross power
spectrum of the three $f_{\rm IHL}(M)$ models shown in Fig. 8 at $z=0.2$
(left), $0.5$ (middle), and $1$ (right). Error bars are the sample variance in
a 200 deg2 survey. The values are scaled by the mean ratio of the IGL$+$IHL
and the IGL flux ($\frac{1}{1-\bar{f}_{\rm IHL}}$), which is the expected
ratio on large scales. Here, we use $\bar{f}_{\rm IHL}=0.2$.

###  4.3 Power Spectrum with Different IHL Profiles

We compare the power spectrum from different IHL profiles while using constant
$f_{\rm IHL}$–$z$ and $f_{\rm IHL}$–$M$ relations in our fiducial model. We
consider an isothermal halo profile, which has a more extended light
distribution compared to the fiducial NFW profile shown in Fig. 10. In Fig. 11
we show the ratio of the IGL$+$IHL and the IGL auto power spectra for the two
models. In a 200 deg2 survey, our forecast shows that the two IHL profile
models can be distinguished with high significance. In reality, IHL is formed
from tidal disruption interactions, and thus the average distribution given in
an intensity mapping observation could be more complicated than the two simple
models considered here. While Fig. 11 is based on auto-spectra, measurements
based on galaxy redshift cross spectra will also be useful in studying IHL
profiles.

Figure 10:  An example of a $5^{\prime}\times 5^{\prime}$ image gridded in our
$7{{}^{\prime\prime}}\times 7{{}^{\prime\prime}}$-sized pixel. We choose a
region centered at a galaxy cluster in MICECAT. Left: IGL intensity field.
Middle: IHL$+$IGL field with fiducial IHL model that follows an NFW profile.
Left: IHL$+$IGL field with fiducial IHL model but using isothermal halo
profile. The isothermal model has a more extended halo profile than the NFW
model. Figure 11:  Ratio of IGL$+$IHL and IHL power spectra using NFW (blue)
and isothermal (orange) IHL profiles with $f_{\rm IHL}=0.1$ (top), $0.2$
(middle), and $0.5$ (bottom). We assume $f_{\rm IHL}$ is independent of
redshift and mass as in our fiducial model. Error bars are the sample variance
in a 200 deg2 survey. We normalize the ratio by the factor $(\frac{1}{1-f_{\rm
IHL}})^{2}$ (the ratio of the square of total intensity), which is the
expected value on large scales.

###  4.4 Effects of the PSF and Masking

In practice, masking to remove pixels contaminated by sources is always
imperfect due to the broadening effect of the PSF. The PSF and pixel masks
result in two systematic effects: (1) residual fluxes outside the mask
(hereafter, referred to as PSF halo) introduce extra intensity and
correlations at the scale of mask size, and (2) pixel masks not only make the
observed image highly irregular but are also correlated with unmasked sources
(see, e.g., Rotti et al., 2021). Hereafter, this effect is referred to as
mask-signal correlation. Both factors introduce mode coupling to the power
spectrum. Therefore, the power spectrum estimation will be biased if these
effects are not properly taken into account in the mode-coupling correction
procedure.

To estimate these systematic effects, we incorporate a Gaussian PSF into the
MICECAT simulated images and then mask pixels. We implement the mode-coupling
correction algorithm from Zemcov et al. (2014), where the mode-coupling matrix
is derived by simulating the response from pure-tone Gaussian random fields.
However, in these simulations, the mask and the signal are not correlated.
Furthermore, residual flux from imperfect masking (PSF halos) can introduce a
bias. We investigate the resulting bias with these simulations. We compare the
power spectrum of four different cases:

  1. 1.

Fiducial case (case 1): no PSF ($\sigma_{\rm PSF}=0$); perfect masking (i.e.,
the masking is performed by removing bright sources). In this case, there is
no systematic error from either the mask-signal correlation or PSF halos.

  2. 2.

Pixel masking-only case (case 2): no PSF ($\sigma_{\rm PSF}=0$); pixel
masking. In this case, there is systematic effect from the mask-signal
correlation but no PSF halos.

  3. 3.

Pixel masking and $\sigma_{\rm PSF}=7{{}^{\prime\prime}}$ case (case 3):
Gaussian PSF with $\sigma_{\rm PSF}=7{{}^{\prime\prime}}$ (size of one pixel);
pixel masking. In this case, there are systematic errors from both mask-signal
correlation and PSF halos.

  4. 4.

Pixel masking and $\sigma_{\rm PSF}=70{{}^{\prime\prime}}$ case (case 4): same
as case 3 but with $\sigma_{\rm PSF}=70{{}^{\prime\prime}}$, which exaggerates
both effects.

The four cases are illustrated in Fig. 12. In case 2–4, we use the same mask
and the same mode-coupling matrix. The masks are circles centered at sources
with a magnitude-dependent radius, and we mask all galaxies brighter than 18
AB magnitude in the Euclid $Y$ band with the radius–magnitude relation from
Zemcov et al. (2014).

Figure 12: Illustration of the four masking cases in Sec. 4.4. The top and
bottom rows show a source on an image before and after masking, respectively.
In cases 1 and 2, there is no PSF so the source is a point source on the
image; in cases 3 and 4, we consider a smaller and a larger Gaussian PSF. Case
1 implements perfect masking that completely removes the flux from the source,
whereas for the other three cases, we use pixel masking and the white circle
denotes the masking region. Cases with a PSF leave low-level residual halos
that extend beyond the mask.

Fig. 13 compares the power spectrum of the four cases described above with
different $f_{\rm IHL}$ values. The power spectrum of the pixel masking-only
case (case 2) is systematically lower than the fiducial case (case 1) below
one-halo scales ($\ell\sim 10^{3}$) due to the unaccounted mask-signal
correlation. The two cases with a non-negligible PSF (cases 3 and 4) have
higher power than the pixel masking-only case at all scales due to the PSF
halos. On large scales ($\ell\lesssim 10^{3}$), the power spectrum is
proportional to the square of total intensity, and therefore the $\sigma_{\rm
PSF}=70{{}^{\prime\prime}}$ case has the highest power because it has more
residual flux outside the mask. On small scales, the power spectrum of the
$\sigma_{\rm PSF}=70{{}^{\prime\prime}}$ case is lower than the $\sigma_{\rm
PSF}=7{{}^{\prime\prime}}$ case due to the uncorrected PSF roll-off. For all
four example scenarios considered here, our simulation suggests that the mask-
signal correlation and PSF halos only cause a $\sim 10\%$ bias on the power
spectrum. However, we note that the bias from these systematic effects will in
the general case depend on the PSF, masking radius, masking depth, mode-
coupling correction algorithm, and the underlying clustering of the signal.

Figure 13:  Top: power spectra of the four cases described in Sec. 4.4 with
four $f_{\rm IHL}$ values. Black: fiducial case (case 1); blue: pixel masking-
only case (case 2); orange: pixel masking and $\sigma_{\rm
PSF}=7{{}^{\prime\prime}}$ case (case 3); green: pixel masking and
$\sigma_{\rm PSF}=70{{}^{\prime\prime}}$ case (case 4). Error bars are the
sample variance in a 200 deg2 survey. Bottom: fractional error with respect to
the fiducial case.

##  5 Conclusion

In this work, we modeled the near-infrared EBL from the IGL and IHL using the
numerical simulation MICECAT, focusing on fluctuation signals at non-linear
clustering scales. The MICECAT galaxy catalog provides the galaxy–halo
relation for quantifying IGL non-linear clustering from satellite galaxies. We
also modeled the signal from the IHL, and investigated the power spectrum
dependence on different IHL models. Finally, we considered the effects of the
PSF and pixel masking, and showed the systematic bias on the power spectrum if
these effects are not properly taken into account. In summary:

  1. 1.

The one-halo clustering from satellite galaxies and the two-halo clustering
have comparable power at tens of arcminute scales ($\ell\sim 10^{3}$).
Previous studies have found excess fluctuations at these scales in near- and
mid-infrared, and explain the excess signal with emission from the EoR
(Kashlinsky et al., 2012) or the low-redshift IHL (Cooray et al., 2012; Zemcov
et al., 2014). However, we find their IGL models underpredict the nonlinear
clustering from satellite galaxies, which can account for some of the excess
power.

  2. 2.

We investigated the IHL power spectrum with different redshift, halo mass, and
IHL profile dependence. Our forecast shows that with a 200 deg2 survey, such
as the SPHEREx deep field, the models considered in this work can be
distinguished with high significance, assuming that sample variance is the
dominant source of uncertainty.

  3. 3.

We simulated the systematic bias on the power spectrum from the mask-signal
correlation and residual bright source fluxes in the presence of PSF and pixel
masking. As a demonstration, we apply the same masking and mode-coupling
correction procedure as Zemcov et al. (2014) to the MICECAT simulated map. Our
results show that the systematic bias from these effects is modest in our
example case. In practice, this bias depends on the specifics of the mask,
PSF, and the mode-coupling correction algorithm, as well as the underlying
clustering of the signal.

While we only explored some simple IHL models in this work, our results show
that upcoming data from SPHEREx promises to probe detailed properties of the
IHL from EBL fluctuation measurements. Here, we detail a few future directions
toward more sophisticated signal models and analysis tools to further improve
IGL and IHL constraints.

First, it is crucial to develop a detailed model of the IGL and IHL that
describes their dependence on redshift, halo mass, and spatial profile. In
this work, we explore this dependence separately by only varying one factor in
each test case. In reality, the IHL properties as a function of redshift, halo
mass, and spatial profile are likely to be a more complicated relation rather
than a separable function. In addition, since the IHL arises from stripped
stars in halos, its properties and evolution history is tightly related to the
IGL. Therefore, it is essential to model the IGL and IHL in a coherent manner.
Hydrodynamic simulations that simultaneously capture large-scale cosmological
clustering and the small-scale physics of the IHL and satellite galaxies will
provide useful models for future EBL measurements.

Furthermore, our current model based on MICECAT is limited to $z<1.4$, which
is not sufficient to model the signal at longer wavelengths ($\sim 3-5$
$\mu$m) where $z\sim 2$ galaxies dominate the IGL and IHL emission.

Modeling the wavelength dependence of the IGL and IHL and other EBL components
is also an essential building block for analyzing data from future
experiments. For example, SPHEREx will carry out a low-resolution spectroscopy
of the EBL from $0.75-5$ $\mu$m. The large number of available band auto and
cross correlations provide a new means to study both the EBL redshift history
and source spectra, so it is essential to have a detailed model of the EBL
component wavelength and redshift dependence.

In addition, to prepare for upcoming EBL observations, it is also important to
explore optimal statistics to extract the desired information from the data.
For example, in Sec. 4.1, we discuss that cross correlation will be a powerful
tool to study the redshift dependence of the EBL. Besides redshift tomography,
we can also cross-correlate EBL maps with tracers selected by a certain
brightness or galaxy type to infer the IGL and IHL associated with those
tracers.

We are grateful to Mike Zemcov for providing analysis on refitting our model
to CIBER data. We thank Asantha Cooray, Tzu-Ching Chang, and the CIBER and
SPHEREx collaborations for helpful discussions. We thank Mike Zemcov, Shuji
Matsuura, and Richard Feder for valuable comments that improve the manuscript.
This work was supported by NASA grants NNX16AJ69G and 80NSSC20K0595. 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

  * Amendola et al. (2018) Amendola, L., Appleby, S., Avgoustidis, A., et al. 2018, Living Reviews in Relativity, 21, 2, doi: 10.1007/s41114-017-0010-3
  * 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
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Beare et al. (2019) Beare, R., Brown, M. J. I., Pimbblet, K., & Taylor, E. N. 2019, ApJ, 873, 78, doi: 10.3847/1538-4357/ab041a
  * Bock et al. (2013) Bock, J., Sullivan, I., Arai, T., et al. 2013, ApJS, 207, 32, doi: 10.1088/0067-0049/207/2/32
  * Caputo et al. (2021) Caputo, A., Vittino, A., Fornengo, N., Regis, M., & Taoso, M. 2021, J. Cosmology Astropart. Phys, 2021, 046, doi: 10.1088/1475-7516/2021/05/046
  * Carleton et al. (2022) Carleton, T., Windhorst, R. A., O’Brien, R., et al. 2022, AJ, 164, 170, doi: 10.3847/1538-3881/ac8d02
  * 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
  * Cheng & Chang (2022) Cheng, Y.-T., & Chang, T.-C. 2022, ApJ, 925, 136, doi: 10.3847/1538-4357/ac3aee
  * Cheng et al. (2021) Cheng, Y.-T., Arai, T., Bangale, P., et al. 2021, ApJ, 919, 69, doi: 10.3847/1538-4357/ac0f5b
  * Contini et al. (2018) Contini, E., Yi, S. K., & Kang, X. 2018, MNRAS, 479, 932, doi: 10.1093/mnras/sty1518
  * 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
  * 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
  * Doré et al. (2014) Doré, O., Bock, J., Ashby, M., et al. 2014, arXiv e-prints, arXiv:1412.4872.  https://arxiv.org/abs/1412.4872
  * Doré et al. (2018) Doré, O., Werner, M. W., Ashby, M. L. N., et al. 2018, arXiv e-prints, arXiv:1805.05489.  https://arxiv.org/abs/1805.05489
  * 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
  * 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
  * 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
  * Gong et al. (2016) Gong, Y., Cooray, A., Mitchell-Wynne, K., et al. 2016, ApJ, 825, 104, doi: 10.3847/0004-637X/825/2/104
  * Helgason et al. (2012) Helgason, K., Ricotti, M., & Kashlinsky, A. 2012, ApJ, 752, 113, doi: 10.1088/0004-637X/752/2/113
  * Hoffmann et al. (2015) Hoffmann, K., Bel, J., Gaztañaga, E., et al. 2015, MNRAS, 447, 1724, doi: 10.1093/mnras/stu2492
  * Kalashev et al. (2019) Kalashev, O. E., Kusenko, A., & Vitagliano, E. 2019, Phys. Rev. D, 99, 023002, doi: 10.1103/PhysRevD.99.023002
  * 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
  * 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
  * Koushan et al. (2021) Koushan, S., Driver, S. P., Bellstedt, S., et al. 2021, MNRAS, 503, 2033, doi: 10.1093/mnras/stab540
  * Lauer et al. (2021) Lauer, T. R., Postman, M., Weaver, H. A., et al. 2021, ApJ, 906, 77, doi: 10.3847/1538-4357/abc881
  * 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. (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
  * Mattila & Väisänen (2019) Mattila, K., & Väisänen, P. 2019, Contemporary Physics, 60, 23, doi: 10.1080/00107514.2019.1586130
  * Mitchell-Wynne et al. (2015) Mitchell-Wynne, K., Cooray, A., Gong, Y., et al. 2015, Nature Communications, 6, 7945, doi: 10.1038/ncomms8945
  * Navarro et al. (1996) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563, doi: 10.1086/177173
  * 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
  * Rotti et al. (2021) Rotti, A., Bolliet, B., Chluba, J., & Remazeilles, M. 2021, MNRAS, 503, 5310, doi: 10.1093/mnras/stab469
  * Saldana-Lopez et al. (2021) Saldana-Lopez, A., Domínguez, A., Pérez-González, P. G., et al. 2021, MNRAS, 507, 5144, doi: 10.1093/mnras/stab2393
  * Salvaterra & Ferrara (2003) Salvaterra, R., & Ferrara, A. 2003, MNRAS, 339, 973, doi: 10.1046/j.1365-8711.2003.06244.x
  * 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
  * Santos et al. (2002) Santos, M. R., Bromm, V., & Kamionkowski, M. 2002, MNRAS, 336, 1082, doi: 10.1046/j.1365-8711.2002.05895.x
  * Schechter (1976) Schechter, P. 1976, ApJ, 203, 297, doi: 10.1086/154079
  * 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
  * Spergel et al. (2015) Spergel, D., Gehrels, N., Baltay, C., et al. 2015, arXiv e-prints, arXiv:1503.03757.  https://arxiv.org/abs/1503.03757
  * Sun et al. (2021) Sun, G., Mirocha, J., Mebane, R. H., & Furlanetto, S. R. 2021, MNRAS, 508, 1954, doi: 10.1093/mnras/stab2697
  * 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. (2013) Tsumura, K., Matsumoto, T., Matsuura, S., Sakon, I., & Wada, T. 2013, PASJ, 65, 121, doi: 10.1093/pasj/65.6.121
  * Windhorst et al. (2018) Windhorst, R. A., Timmes, F. X., Wyithe, J. S. B., et al. 2018, ApJS, 234, 41, doi: 10.3847/1538-4365/aaa760
  * Yang et al. (2009) Yang, X., Mo, H. J., & van den Bosch, F. C. 2009, ApJ, 695, 900, doi: 10.1088/0004-637X/695/2/900
  * Yue et al. (2013) Yue, B., Ferrara, A., Salvaterra, R., Xu, Y., & Chen, X. 2013, MNRAS, 433, 1556, doi: 10.1093/mnras/stt826
  * Yue et al. (2014) —. 2014, MNRAS, 440, 1263, doi: 10.1093/mnras/stu351
  * Zemcov et al. (2017) Zemcov, M., Immel, P., Nguyen, C., et al. 2017, Nature Communications, 8, 15003, doi: 10.1038/ncomms15003
  * Zemcov et al. (2014) Zemcov, M., Smidt, J., Arai, T., et al. 2014, Science, 346, 732, doi: 10.1126/science.1258168