# Cosmic Near-Infrared Background Tomography with SPHEREx Using Galaxy Cross-
Correlations

Yun-Ting Cheng<EMAIL_ADDRESS>Tzu-Ching Chang

###### Abstract

The extragalactic background light (EBL) consists of integrated light from all
sources of emission throughout the history of the Universe. At near-infrared
wavelengths, the EBL is dominated by stellar emission across cosmic time;
however, the spectral and redshift information of the emitting sources is
entangled and cannot be directly measured by absolute photometry or
fluctuation measurements. Cross-correlating near-infrared maps with tracers of
known redshift enables EBL redshift tomography, as EBL emission will only
correlate with external tracers from the same redshift. Here we forecast the
sensitivity of probing the EBL spectral energy distribution as a function of
redshift by cross-correlating the upcoming near-infrared spectro-imaging
survey, SPHEREx, with several current and future galaxy redshift surveys.
Using a model galaxy luminosity function, we estimate the cross-power spectrum
clustering amplitude on large scales, and forecast that the near-infrared EBL
spectrum can be detected tomographically out to $z\sim 6$. We also predict a
high significance measurement ($\sim 10^{2}-10^{4}\sigma$) of the small-scale
cross-power spectrum out to $z\sim 10$. The amplitudes of the large-scale
cross-power spectra can constrain the cosmic evolution of the stellar
synthesis process through both continuum and the line emission, while on the
non-linear and Poisson-noise scales, the high sensitivity measurements can
probe the mean spectra associate with the tracer population across redshift.
111ⓒ 2021. All rights reserved.

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

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) is the aggregate light from all
sources of emission across cosmic time. EBL measurements have been made from
gamma-rays to radio (see e.g., Cooray, 2016; Mattila & Väisänen, 2019, for
recent reviews). At near-infrared wavelengths, the EBL is mostly produced by
redshifted ultraviolet and optical stellar emission, and thus carries
essential information on the history of stellar synthesis processes in our
universe. However, observations also suggest other sources of the near-
infrared EBL, including diffuse light in the dark matter halos from stripped
stars, sometimes referred to as intra-halo light (IHL, Cooray et al., 2012;
Zemcov et al., 2014; Cheng et al., 2021), and early formation of stars and
galaxies from the epoch of reionization and cosmic dawn (Kashlinsky et al.,
2005; Matsumoto et al., 2011; Kashlinsky et al., 2012; Mitchell-Wynne et al.,
2015). Possible near-infrared EBL emission from 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 candidate (e.g.,
Gong et al., 2016; Caputo et al., 2020) have also been proposed.

Extensive studies have attempted to probe the near-infrared EBL through
different methods. Galaxy counts measures emission from resolved galaxies and
the observations are extrapolated to estimate contributions from faint
galaxies below the detection limit (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). This sets a lower bound to the integrated galaxy light
(IGL) component of EBL, as EBL from diffuse emission or sources not associated
with galaxies are not included. Direct measurements using absolute photometry
capture all emission in the EBL (Bernstein, 2007; Levenson et al., 2007;
Tsumura et al., 2013; Matsumoto et al., 2015; Sano et al., 2015; Zemcov et
al., 2017; Matsuura et al., 2017; Sano et al., 2020; Lauer et al., 2020),
however, absolute photometry is challenging since systematic errors and
foreground contamination have to be tightly controlled. Fluctuation analysis
is an alternative approach that is less susceptible to foreground
contamination (e.g. zodiacal light, Galactic cirrus), since foregrounds and
the EBL signal have distinct spatial and spectral correlation features
(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); nevertheless inferring the absolute intensity of the EBL from
fluctuation measurements depends on model assumptions. Another limitation of
absolute photometry and fluctuation measurements is the inferred EBL redshift
resolution. As the measured intensity is a projection of all emission along
the line of sight, one way to infer the redshift dependency of the EBL is to
measure the opacity of gamma-ray photons from individual blazars. Near-
infrared photons along the line of sight will interact with gamma-ray photons
by pair production, and thus the redshift composition of EBL can be
constrained by the absorption features in blazar spectra (e.g., Aharonian et
al., 2006, 2007; MAGIC Collaboration et al., 2008; Abdo et al., 2010;
Ackermann et al., 2012; H. E. S. S. Collaboration et al., 2017; Ackermann et
al., 2018; Abeysekara et al., 2019; Acciari et al., 2019; Abdalla et al.,
2020). This method enables a redshift tomography of the EBL by observing
blazars from different distances. However, these estimates have low spectral
resolution and depend on the assumption of intrinsic blazar spectra.

As the EBL intensity contains a mixture of emitting sources from all
redshifts, and different sources have distinct spectral features, it will be
insightful to decompose the EBL signal by its redshift and spectral
dependencies. Cross-correlation provides a path to perform EBL tomography
while retaining spectral information. By cross-correlating with external
sources (e.g., a galaxy catalog) that trace the underlying large-scale
structure at a certain redshift, one can extract the IGL and other components
associated with the large-scale structure, such as the IHL, on large scales.
The cross-correlation technique has previously been proposed and applied to
derive redshift distribution of photometric catalogs (Newman, 2008; Ménard et
al., 2013; McQuinn & White, 2013), perform redshift tomography of broadband
imaging surveys (Schmidt et al., 2015; Chiang & Ménard, 2019; Chiang et al.,
2019), probe the cosmic infrared background (Serra et al., 2014) and the
thermal or kinematic Sunyaev-Zel’dovich signal from the cosmic microwave
background (e.g., Hill et al., 2016; Chiang et al., 2020). Cross correlation
has also been proposed to probe the direct collapse blackholes (Cappelluti et
al., 2013; Mitchell-Wynne et al., 2016) and the axion decay (Creque-
Sarbinowski & Kamionkowski, 2018) in the optical and near-infrared EBL. In
line intensity mapping (see Kovetz et al., 2017, for a review), where the
emission from a certain spectral line is used to trace the three-dimensional
large-scale structure, cross-correlation is also a powerful technique to
measure line emission from 21 cm (e.g., Chang et al., 2010; Masui et al.,
2013), CO (e.g., Pullen et al., 2013), [C ii] (e.g., Pullen et al., 2018; Yang
et al., 2019), and Ly$\alpha$ (e.g., Croft et al., 2016, 2018) with much
reduced impact from foreground and systematics contaminations.

SPHEREx (Doré et al., 2014, 2018a) is an approved NASA MIDEX mission that will
carry out an all-sky spectro-imaging survey at near-infrared wavelengths (see
Sec. 5 for details). SPHEREx will achieve unparalleled spectral resolution,
sensitivity, and sky coverage for spectrally mapping the near-infrared sky,
which provides an exceptional dataset for studying EBL through fluctuation
analysis.

Next generation spectroscopic and photometric galaxy surveys will come online
in the next few years, including DESI (DESI Collaboration et al., 2016), the
Rubin Observatory LSST (LSST Science Collaboration et al., 2009), Euclid
(Amendola et al., 2018), and the Roman Space Telescope (Spergel et al., 2015).
These surveys will provide galaxy catalogs in unprecedented depth and synergy
with the SPHEREx spectro-imaging dataset. In this work, we forecast
constraints on EBL enabled by cross-correlations of SPHEREx and upcoming
galaxy surveys. SPHEREx spectro-images encode the spatial and spectral
information of the EBL, whereas galaxies trace the three-dimensional cosmic
structures with redshift information. With cross-correlation, we perform
redshift tomography of the EBL spectrum, which constrains the redshift
composition and spectral response of the integrated background emission, as
well as its spatial scale dependence. These measurements primarily consist of
the IGL, and will provide valuable insight on the cosmic star formation and
stellar mass history. Furthermore, information on the IHL or other possible
EBL emitting sources that trace the large-scale structure (e.g. the decaying
dark matter candidates) can also be revealed by the cross-correlation
measurements.

A recent study by Scott et al. (2021) investigated the ultraviolet to optical
EBL constraints from cross correlation, using images from the future
ultraviolet survey Cosmological Advanced Survey Telescope for Optical-
ultraviolet Research (CASTOR; Cote et al., 2019) for lower redshift, and
SPHEREx for the epoch of reionization. Their study focuses on the UV broadband
emission and the Ly$\alpha$ line. In this work, we conduct a comprehensive
forecast of constraining the EBL spectra from SPHEREx with all accessible
wavelength and redshift ranges.

This paper is organized as follows. We first introduce the redshift tomography
formalism of cross-power spectrum estimation in Sec. 2. We model components of
the EBL signal in Sec. 3. Sec. 4 details the evaluation metrics and the
fiducial cases we choose to present the results. Sec. 5 introduces SPHEREx,
and Sec. 6 describes the galaxy surveys considered in this work. The results
are presented in Sec. 7. Discussion and further science interpretations are
given in Sec. 8 and Sec. 9, respectively. Sec. 10 concludes 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.

##  2 Formalism

We detail the analytical expression for a cross power spectrum, and its error
estimation, between a continuous density field (such as a diffuse intensity
map) and a discrete density tracer sample (such as a galaxy catalog). An
intensity field $I$ and a galaxy tracer field $g$ can each be expanded in
terms of spherical harmonics $Y_{\ell m}$ as

$\delta I\left(\hat{n}\right)=\sum_{\ell,m}a^{I}_{\ell m}Y_{\ell
m}\left(\hat{n}\right),$ (1)

and

$\delta g\left(\hat{n}\right)=\sum_{\ell,m}a^{g}_{\ell m}Y_{\ell
m}\left(\hat{n}\right).$ (2)

The angular cross power spectrum of these two fields, $C_{\ell}^{\rm x}$,
binned in a width of $\Delta\ell$, is given by

$C_{\ell}^{\rm x}=\left\langle\frac{1}{\left(2\ell+1\right)f_{\rm
sky}}\sum_{m=-\ell}^{\ell}\left(a^{I}_{\ell m}\right)^{*}a^{g}_{\ell
m}\right\rangle_{\ell\in\Delta\ell},$ (3)

with a variance of

$\left(\delta C_{\ell}^{\rm
x}\right)^{2}=\frac{1}{N_{\ell}}\left[(C_{\ell}^{\rm
x})^{2}+C_{\ell}^{I}C_{\ell}^{g}\right],$ (4)

where $N_{\ell}=\Delta\ell\left(2\ell+1\right)f_{\rm sky}$ is the number of
$\ell$ modes in the band power spectrum at $\ell$ with a width of
$\Delta\ell$, and $f_{\rm sky}$ is the fraction of sky area used in the cross
correlation. $C_{\ell}^{g}$, $C_{\ell}^{I}$, and $C_{\ell}^{\rm x}$ are the
galaxy auto power spectrum, intensity field auto power spectrum, and their
cross power spectrum, respectively.

Eq. 4 is the Gaussian variance that assumes all multipole modes are
independent. This is a good approximation on large scales, but on small scales
in the Poisson regime, non-Gaussian contribution to the covariance needs to be
accounted for. Therefore, we also consider the trispectrum contribution to the
variance. The cross-power spectrum covariance is thus

${\rm Cov}\left(C_{\ell}^{\rm x},C_{\ell^{\prime}}^{\rm
x}\right)=\delta^{K}_{\ell\ell^{\prime}}\left(\delta C_{\ell}^{\rm
x}\right)^{2}+\frac{\mathcal{T}^{gIgI}(\ell,\ell^{\prime})}{\Omega_{\rm
sur}},$ (5)

where $\delta^{K}$ is the Kronecker delta, $\Omega_{\rm sur}=4\pi f_{\rm sky}$
is the survey area in solid angle, and
$\mathcal{T}^{gIgI}(\ell,\ell^{\prime})$ is the trispectrum of the galaxy and
intensity field (Schaan et al., 2018).

###  2.1 Galaxy Sample Auto Power Spectrum $C_{\ell}^{g}$

We describe a galaxy auto power spectrum as comprised of clustering and (shot)
Poisson-noise terms:

$C_{\ell}^{g}=C_{\ell,\rm clus}^{g}+C_{\ell,\rm P}^{g}.$ (6)

Using the Limber approximation, the clustering term can be expressed as

$C_{\ell,\rm clus}^{g}=\int
dz\frac{H(z)}{c}\frac{f^{g}(z)f^{g}(z)}{\chi^{2}(z)}b_{g}^{2}(z)P_{\rm
m}\left(k=\frac{\ell+\frac{1}{2}}{\chi(z)},z\right),$ (7)

where $H$, $\chi$, $c$, $f^{g}(z)$ are the Hubble parameter, co-moving
distance, speed of light, and the galaxy selection function, respectively.
$P_{m}$ is the matter power spectrum. In this work, we focus on the linear
clustering regime and restrict our analysis to large scales at $k<0.2$ h/Mpc.
We assume $b_{g}$ is a scale-independent linear bias. In practice, with an
assumed cosmological model, $b_{g}$ can be measured from the amplitude of the
galaxy auto spectrum. We consider the limit where galaxies are from a narrow
redshift bin ($\Delta z^{g}\ll 1$) centered at $z^{g}$, and use the following
approximation:

$f^{g}(z)=\begin{cases}\frac{1}{\Delta z^{g}}\quad\quad\text{if}\quad
z_{g}-\frac{\Delta z^{g}}{2}<z<z_{g}+\frac{\Delta z^{g}}{2}\\\
0\quad\quad\quad\text{Otherwise},\end{cases}$ (8)

and

$C_{\ell,\rm clus}^{g}=\frac{1}{\Delta
z^{g}}\frac{H(z^{g})}{c\chi^{2}(z)}b_{g}^{2}(z)P_{\rm
m}(k=\frac{\ell+\frac{1}{2}}{\chi(z^{g})},z^{g}).$ (9)

The galaxy Poisson noise power is the reciprocal of the surface number density
$dN_{g}/d\Omega$,

$C_{\ell,\rm P}^{g}=\left(\frac{dN_{g}}{d\Omega}\right)^{-1}=\left(\Delta
z^{g}\frac{dN_{g}}{dzd\Omega}\right)^{-1},$ (10)

where the second equality uses the same assumption of a narrow redshift bin.

###  2.2 Intensity Field Auto Power Spectrum $C_{\ell}^{I}$

The angular auto power spectrum of an intensity field includes contributions
from clustering, Poisson noise, and instrument noise,

$C_{\ell}^{I}=C_{\ell,\rm clus}^{I}+C_{\ell,\rm P}^{I}+C_{\ell,n}^{I}.$ (11)

Using the Limber approximation, the clustering power spectrum of the intensity
map at frequency $\nu$ is

$\begin{split}C_{\ell,\rm clus}^{I}&=\int
dz\frac{H(z)}{c}\frac{f^{I}(z)f^{I}(z)}{\chi^{2}(z)}b_{I}^{2}(z,\nu)\\\
&\cdot\left[\frac{d(\nu I_{\nu})(z,R(\nu))}{dz}\right]^{2}P_{\rm
m}\left(k=\frac{\ell+\frac{1}{2}}{\chi(z)},z\right),\end{split}$ (12)

where $b_{I}$ is the large-scale bias, and $f^{I}$ the redshift selection
function. We take $f^{I}(z)=1$ for simplicity as the intensity field is
comprised from emission over a large redshift range. $\nu I_{\nu}(z,R(\nu))$
is the intensity of emitting sources at redshift $z$ observed by a filter with
frequency response function $R(\nu)$, i.e.,

$\frac{d(\nu I_{\nu})}{dz}(z,R(\nu))=\frac{1}{\int d\nu R(\nu)}\int d\nu
R(\nu)\frac{d(\nu I_{\nu})}{dz}(z,\nu),$ (13)

where $\nu I_{\nu}(z,\nu)$ (in units of nW m-2 sr-1) is the intrinsic spectral
energy distribution of the EBL from sources at redshift $z$ observed at
frequency $\nu$. In this work, we consider filters with a narrow spectral
width ($\Delta\nu/\nu\ll 1$), and use the approximation:

$\frac{d(\nu I_{\nu})}{dz}(z,R(\nu))=\frac{d(\nu I_{\nu})}{dz}(z,\nu).$ (14)

The EBL is the aggregate intensity from all sources across cosmic time,
defined as

$\nu I_{\nu}(\nu)=\int dz\frac{d(\nu I_{\nu})}{dz}(z,\nu).$ (15)

$d(\nu I_{\nu})/dz$ is the redshift-dependent EBL intensity in the same unit
of as $\nu I_{\nu}$. In the optical to near-infrared wavelengths, the emitting
sources contributing to the measured intensity include galaxies, quasars,
stars, from all luminosity ranges and from all redshifts (e.g., Cooray, 2016;
Mattila & Väisänen, 2019).

We model the integrated galaxy light (IGL) as the main emitting source of the
EBL, which can be considered as a lower bound to $\nu I_{\nu}(z,\nu)$. We
model the continuum and spectral line emission of the IGL spectrum separately,

$\begin{split}\left.\frac{d(\nu I_{\nu})}{dz}\right|_{\rm
IGL}(z,\nu)&=\left.\frac{d(\nu I_{\nu})}{dz}\right|_{\rm cont}(z,\nu)\\\
&+\sum_{\rm line}\left.\frac{d(\nu I_{\nu})}{dz}\right|_{\rm
line}(z,\nu).\end{split}$ (16)

The continuum can be expressed in terms of a volume averaged galaxy
emissivity, consisting of the galaxy luminosity function $\Phi(m,z,\nu_{\rm
rf})=dN/dV/dm$, where $N$ is the galaxy number count, $V$ is the co-moving
volume, $m$ is the AB magnitude, and $\nu_{\rm rf}$ is the rest-frame
frequency. The spectrum of the continuum intensity can be written as

$\begin{split}\left.\frac{d(\nu I_{\nu})}{dz}\right|_{\rm
cont}(z,\nu)&=\int_{m_{\rm th}(z,\nu)}^{\infty}dm\Phi(m,z,\nu(1+z))\\\
&\cdot\nu F_{\nu}(m)\frac{d\chi}{dz}(z)D_{A}^{2}(z),\end{split}$ (17)

where $F_{\nu}(m)=3631\cdot 10^{m/2.5}$ Jy is the specific flux density at
magnitude $m$, $d\chi/dz=c/H(z)$, and $D_{A}$ is the co-moving angular
diameter distance, which equals the co-moving distance in a flat universe. The
lower bound of the magnitude integration, $m_{\rm th}(z,\nu)$, is the masking
magnitude threshold, below which sources are masked in the diffuse intensity
map. Masking reduces the foreground contributions to the cross power spectrum
Poisson noise, as discussed in Sec. 8.2. The models for $\Phi$, $b_{I}$, and
the choice of $m_{\rm th}$ used in this work are detailed in Sec. 3.

For spectral lines, we build the emission model based on the halo model
formalism and the mass $M$ and line luminosity $L_{\rm line}$ relation:

$\begin{split}\left.\frac{d(\nu I_{\nu})}{dz}\right|_{\rm
line}&(z,\nu)=\int^{M_{\rm th}(z)}_{M_{\rm min}}dM\frac{dn}{dM}(M,z)\frac{\nu
L_{\rm line}(M,z)}{4\pi D_{L}^{2}(z)}\\\
&\cdot\frac{d\chi}{dz}(z)D_{A}^{2}(z)\delta^{D}(\nu-\nu_{\rm rf}^{\rm
line}/(1+z)),\end{split}$ (18)

where $dn/dM$ is the halo mass function (Sheth et al., 2001), $D_{L}$ is the
luminosity distance, $\delta^{D}$ is the Dirac delta function, and $\nu_{\rm
rf}^{\rm line}$ is the rest frame frequency of the spectral line. Here we
treat the intrinsic spectral line width as a delta function, as the intrinsic
line width in a galaxy is unresolved in SPHEREx’s low-resolution spectral
bands. The lower limit of the integration, $M_{\rm min}$, is the minimum halo
mass hosting sources with line emission, and we set $M_{\rm
min}=10^{8}M_{\odot}$/h. The upper limit, $M_{\rm th}(z)$, is determined by
the line masking threshold and detailed in Sec. 8.2.

The line intensity measured at frequency band $\nu$ with width $\Delta\nu$ is

$\left.\frac{d(\nu I_{\nu})}{dz}\right|_{\rm
line}(z,\nu,\Delta\nu)=\frac{1}{\Delta z^{\rm line}}\left.\nu
I_{\nu}\right|_{\rm line}(z^{\rm line}),$ (19)

where

$\begin{split}\left.\nu I_{\nu}\right|_{\rm line}(z)&=\int^{M_{\rm
th}(z)}_{M_{\rm min}}dM\frac{dn}{dM}(M,z)\\\ &\cdot\frac{\nu L_{\rm
line}(M,z)}{4\pi D_{L}^{2}(z)}\frac{d\chi}{d\nu}(z)D_{A}^{2}(z),\end{split}$
(20)

$z^{\rm line}=\nu_{\rm rf}^{\rm line}/\nu-1$, $\Delta z^{\rm
line}=\Delta\nu\left(\nu_{\rm rf}^{\rm line}/\nu^{2}\right)$, and
$d\chi/d\nu=c(1+z)^{2}/(\nu_{\rm rf}^{\rm line}H(z))$. Thus,

$\begin{split}\left.\frac{d(\nu I_{\nu})}{dz}\right|_{\rm
line}(z,\nu,\Delta\nu)&=\frac{c}{4\pi(1+z)^{2}H(z)}\frac{\nu}{\Delta\nu}\\\
&\cdot\int^{M_{\rm th}(z)}_{M_{\rm min}}dM\frac{dn}{dM}(M,z)L_{\rm
line}(M,z).\\\ \end{split}$ (21)

Note the line intensity is inversely proportional to the observed spectral
width $\Delta\nu$.

The intensity Poisson noise comprises galaxy continuum and line emission, and
foreground emission. Here we consider Galactic stars as foreground sources
that contribute to the fluctuation power spectrum, and thus

$C_{\ell,\rm P}^{I}=\left.C_{\ell,\rm P}^{I}\right|_{\rm cont}+\sum_{\rm
line}\left.C_{\ell,\rm P}^{I}\right|_{\rm line}+\left.C_{\ell,\rm
P}^{I}\right|_{\rm star}.$ (22)

In the limit of narrow spectral channel width, the Poisson-noise components
are given by

$\begin{split}\left.C_{\ell,\rm P}^{I}\right|_{\rm cont}=\int dz&\int_{m_{\rm
th}(z,\nu)}^{\infty}dm\Phi(m,z,\nu(1+z))\\\ &\cdot\left[\nu
F_{\nu}(m)\right]^{2}\frac{d\chi}{dz}(z)D_{A}^{2}(z),\end{split}$ (23)

$\begin{split}&\left.C_{\ell,\rm P}^{I}\right|_{\rm line}=\sum_{\rm
line}\frac{1}{\Delta z^{\rm line}}\frac{H(z^{\rm line})}{c\chi^{2}(z^{\rm
line})}\int^{M_{\rm th}(z)}_{M_{\rm min}}dM\\\ &\cdot\frac{dn}{dM}(M,z^{\rm
line})\left[\frac{\nu L_{\rm line}(M,z^{\rm line})}{4\pi D_{L}^{2}(z^{\rm
line})}\frac{d\chi}{d\nu}(z^{\rm line})D_{A}^{2}(z^{\rm
line})\right]^{2},\end{split}$ (24)

and

$\left.C_{\ell,\rm P}^{I}\right|_{\rm star}=\int_{m_{\rm
th}(z,\nu)}^{\infty}dm\frac{dN_{\rm star}}{dmd\Omega}(m,\nu)\left[\nu
F_{\nu}(m)\right]^{2},$ (25)

where $dN_{\rm star}/dm/d\Omega$ is the star count per magnitude per solid
angle.

Here we ignore the correlation between continuum and and line emissions, and
model them as independent source populations that contribute to the Poisson
noise, while in reality the continuum and line emissions can come from the
same or nearby sources that are spatially correlated. The impact of this
correlation is negligible in narrow band measurements ($\Delta\nu/\nu\ll 1$)
considered in this work, as the line emission in a given frequency channel
only comes from a narrow redshift range ($\Delta z^{\rm line}\ll 1$).

We assume instrument noise to be Gaussian with a variance
$\sigma_{n}^{2}(\nu)$, so the noise power spectrum is given by

$C_{\ell,n}^{\rm I}=\sigma_{n}^{2}(\nu)\Omega_{\rm pix},$ (26)

where $\Omega_{\rm pix}$ is the pixel size. Note that $\sigma_{n}^{2}(\nu)$
has the same units as $\nu I_{\nu}$.

###  2.3 Galaxy-Intensity Cross Power Spectrum $C_{\ell}^{x}$

The cross power spectrum also contains clustering and Poisson noise
components:

$C_{\ell}^{\rm x}=C_{\ell,\rm clus}^{\rm x}+C_{\ell,\rm P}^{\rm x}.$ (27)

The clustering term in the Limber approximation has the expression:

$\begin{split}C_{\ell,\rm clus}^{\rm x}&=\int
dz\frac{H(z)}{c}\frac{f^{g}(z)f^{I}(z)}{\chi^{2}(z)}r_{\rm
x}(z)b_{g}(z)b_{I}(z)\\\ &\cdot\frac{d(\nu I_{\nu})}{dz}(z,\nu)P_{\rm
m}(k=\frac{\ell+\frac{1}{2}}{\chi(z)},z),\end{split}$ (28)

where $r_{\rm x}$ is the cross correlation coefficient between galaxy sample
and intensity field. We set $r_{\rm x}=1$ throughout this work as we consider
linear scales only.

We separate the clustering term into continuum and line emission,

$C_{\ell,\rm clus}^{\rm x}=\left.C_{\ell,\rm clus}^{\rm x}\right|_{\rm
cont}+\sum_{\rm line}\left.C_{\ell,\rm clus}^{\rm x}\right|_{\rm line}.$ (29)

In the limit of a narrow observed spectral channel ($\Delta\nu/\nu\ll 1$) and
redshift width ($\Delta z^{g}\ll 1$), we get

$\begin{split}\left.C_{\ell,\rm clus}^{\rm x}\right|_{\rm
cont}&=\frac{H(z^{g})}{c\chi^{2}(z^{g})}b_{g}(z^{g})b_{I,\rm cont}(z^{g})\\\
&\cdot\left.\frac{d(\nu I_{\nu})}{dz}\right|_{\rm cont}(z^{g},\nu)P_{\rm
m}\left(k=\frac{\ell+\frac{1}{2}}{\chi(z^{g})},z^{g}\right),\end{split}$ (30)

$\begin{split}\left.C_{\ell,\rm clus}^{\rm x}\right|_{\rm
line}&=\frac{H(z^{g})}{c\chi^{2}(z^{g})}\frac{\Delta z^{\rm g,line}}{\Delta
z^{g}\Delta z^{\rm line}}b_{g}(z^{g})b_{I,\rm line}(z^{g})\\\ &\cdot\left.\nu
I_{\nu}\right|_{\rm line}(z^{g})P_{\rm
m}\left(k=\frac{\ell+\frac{1}{2}}{\chi(z^{g})},z^{g}\right),\end{split}$ (31)

where $\Delta z^{\rm g,line}$ is the redshift range where the galaxy tracer
redshift $z^{g}$ overlaps with the line redshift $z^{\rm line}$.

The Poisson cross-power spectrum is proportional to the total intensity of the
tracer galaxies; in real space, this is nearly equivalent to stacking the
intensity at the position of tracer galaxies, or the conditional correlation
function. The Poisson noise cross power spectrum can similarly be written as
the sum of continuum and line contributions:

$\begin{split}C_{\ell,\rm P}^{\rm
x}&=\left(\frac{dN_{g}}{d\Omega}\right)^{-1}\Delta z^{g}\left.\frac{d(\nu
I_{\nu})}{dz}\right|_{g}(z^{g},\nu)\\\ &=\left.C_{\ell,\rm P}^{\rm
x}\right|_{\rm cont}+\left.C_{\ell,\rm P}^{\rm x}\right|_{\rm
line}.\end{split}$ (32)

We use “$g$” to represent the intensity of the tracer (galaxy) sample. The
continuum and line Poisson noise components are given by

$\left.C_{\ell,\rm P}^{\rm x}\right|_{\rm
cont}=\left(\frac{dN_{g}}{d\Omega}\right)^{-1}\Delta z^{g}\left.\frac{d(\nu
I_{\nu})}{dz}\right|_{{\rm cont},g}(z^{g},\nu),$ (33)

and

$\left.C_{\ell,\rm P}^{\rm x}\right|_{\rm
line}=\left(\frac{dN_{g}}{d\Omega}\right)^{-1}\frac{\Delta z^{g,\rm
line}}{\Delta z^{\rm line}}\left.\nu I_{\nu}\right|_{{\rm
line},g}(z^{g},\nu),$ (34)

where the last term in both equations is the integration of Eq. 17 and 20 over
the galaxy sample, respectively. Here we assume galaxy tracers are the
brightest continuum and line emitting sources, and thus the tracer population
in our model corresponds to all the galaxies above a limiting magnitude
$m_{g}(z,\nu)$ and its corresponding halo mass $M_{g}(z)$. For a given masking
threshold $m_{\rm th}(z,\nu)$ (corresponding to halo mass threshold $M_{\rm
th}(z)$), if the tracer catalog is deeper than the masking threshold, i.e.,
$m_{g}>m_{\rm th}$ ($M_{g}<M_{\rm th}$), we get

$\begin{split}\left.\frac{d(\nu I_{\nu})}{dz}\right|_{{\rm
cont},g}(z,\nu)=\int_{m_{\rm
th}(z,\nu)}^{m_{g}(z,\nu)}&dm\Phi(m,z,\nu(1+z))\\\ \cdot\nu
F_{\nu}(m)&\frac{d\chi}{dz}(z)D_{A}^{2}(z),\end{split}$ (35)

and

$\begin{split}\left.\nu I_{\nu}\right|_{{\rm line},g}(z,\nu)=\int^{M_{\rm
th}(z)}_{M_{g}(z)}dM&\frac{dn}{dM}(M,z)\\\ \cdot\frac{\nu L_{\rm
line}(M,z)}{4\pi D_{L}^{2}(z)}&\frac{d\chi}{d\nu}(z)D_{A}^{2}(z).\end{split}$
(36)

On the other hand, if $m_{g}\leqslant m_{\rm th}$ ($M_{g}\geqslant M_{\rm
th}$), all tracer galaxies are masked, and thus

$\left.\frac{d(\nu I_{\nu})}{dz}\right|_{{\rm cont},g}(z,\nu)=\left.\nu
I_{\nu}\right|_{{\rm line},g}(z,\nu)=0,$ (37)

and the cross Poisson power vanishes.

###  2.4 Galaxy-Intensity Cross Trispectrum $\mathcal{T}^{gIgI}$

We model the non-Gaussian contribution in the Poisson regime using the
trispectrum between the galaxy and intensity fields, given by (Schaan et al.,
2018)

$\mathcal{T}^{gIgI}(\ell,\ell^{\prime})=\left[\left.C_{\ell,\rm
P}^{I}\right|_{{\rm cont},g}+\sum_{\rm line}\left.C_{\ell,\rm
P}^{I}\right|_{{\rm line},g}\right]\cdot\left(C_{\ell,\rm P}^{g}\right)^{2},$
(38)

where $C_{\ell,\rm P}^{g}$ is described by Eq. 10, and $\left.C_{\ell,\rm
P}^{I}\right|_{\rm cont,g}$ and $\left.C_{\ell,\rm P}^{I}\right|_{\rm line,g}$
are the integration of Eq. (23) and Eq. (24) over the galaxy sample,
respectively. In the Poisson regime, the
$\mathcal{T}^{gIgI}(\ell,\ell^{\prime})$ function is a constant for all values
of $\ell$ and $\ell^{\prime}$ pairs, and specifies the off-diagonal terms of
the covariance matrix (Eq. 5) that describes the correlation between different
multipole modes. If $m_{g}>m_{\rm th}$ ($M_{g}<M_{\rm th}$),

$\begin{split}\left.C_{\ell,\rm P}^{I}\right|_{{\rm cont},g}=\int
dz\int_{m_{\rm th}(z,\nu)}^{m_{g}(z,\nu)}&dm\Phi(m,z,\nu(1+z))\\\
&\cdot\left[\nu
F_{\nu}(m)\right]^{2}\frac{d\chi}{dz}(z)D_{A}^{2}(z)\end{split}$ (39)

and

$\begin{split}\left.C_{\ell,\rm P}^{I}\right|_{{\rm line},g}&=\sum_{\rm
line}\frac{1}{\Delta z^{\rm line}}\frac{H(z^{\rm line})}{c\chi^{2}(z^{\rm
line})}\\\ &\cdot\int^{M_{\rm th}(z)}_{M_{g}(z)}dM\cdot\frac{dn}{dM}(M,z^{\rm
line})\\\ &\cdot\left[\frac{\nu L_{\rm line}(M,z^{\rm line})}{4\pi
D_{L}^{2}(z^{\rm line})}\frac{d\chi}{d\nu}(z^{\rm line})D_{A}^{2}(z^{\rm
line})\right]^{2}.\end{split}$ (40)

Otherwise, if $m_{g}\geqslant m_{\rm th}$ ($M_{g}\leqslant M_{\rm th}$), we
get

$\left.C_{\ell,\rm P}^{I}\right|_{{\rm cont},g}=\left.C_{\ell,\rm
P}^{I}\right|_{{\rm line},g}=\mathcal{T}^{gIgI}=0.$ (41)

##  3 Emission Model

We assume the integrated galaxy light (IGL) constitutes the majority of the
near-IR emission, and model the spectral energy distribution as a continuum
and spectral line components. In addition, we consider a model for quasars to
study constraints on quasar emission in the Poisson regime. Note we do not add
the quasar model to the IGL model, since IGL, as described below, already
includes the quasar contribution. We also consider Galactic stars as a
foreground emission.

###  3.1 Galaxy Spectral Energy Distribution

####  3.1.1 Continuum Emission

Our model for the IGL continuum emission is based on a luminosity function
prescription from Helgason et al. (2012), where the luminosity function,
$\Phi(m,z,\nu_{\rm rf})$, has a Schechter functional form (Schechter, 1976)
that depends on magnitude $m$, redshift $z$, and rest-frame frequency
$\nu_{\rm rf}$. The parameters are calibrated to observations in several
ultraviolet to mid-infrared wavelength bands, and we interpolate the Schechter
parameters to all wavelengths from their model values.

Our model for the continuum intensity bias, $b_{I,\rm cont}(z,\nu)$, is also
from Helgason et al. (2012), who uses a halo occupation distribution (HOD)
framework (Zheng et al., 2005) and assumes no frequency dependence:

$b_{I,\rm cont}(z)=\frac{\int dM\frac{dn}{dM}(M,z)b(M,z)\left\langle N_{\rm
gal}\right\rangle}{\int dM\frac{dn}{dM}(M,z)\left\langle N_{\rm
gal}\right\rangle},$ (42)

where $\left\langle N_{\rm gal}\right\rangle$ is the total halo occupation
number. Following Helgason et al. (2012), we use the HOD parameters given by
the SDSS measurement in Zehavi et al. (2011).

####  3.1.2 Line Emission

In this work, we consider five prominent optical spectral lines: H$\alpha$
(656.3 nm), [O iii] (500.7 nm), H$\beta$ (486.1 nm), [O ii] (372.7 nm), and
Ly$\alpha$ (121.6 nm). These lines are of particular interest as they are
visible in SPHEREx bands across a range of redshifts.

Our line emission model is built on an empirical halo mass and line luminosity
relation $L_{\rm line}(M,z)$. We adopt and summarize the prescription from
Gong et al. (2017) below. First, we assume a simple linear relation between
star formation rate (SFR) and halo mass, and calibrate the scaling factor
using the star formation rate density constraints from Hopkins & Beacom
(2006). We then further assume a linear relation between the line luminosity
and SFR, using SFR–$L_{\rm line}$ relations from Kennicutt (1998) and Ly et
al. (2007) for H$\alpha$, [O ii] and [O iii] emission:

$\displaystyle\frac{\rm SFR}{M_{\odot}/\rm{yr}}$ $\displaystyle=7.9\times
10^{-42}\frac{L_{H\alpha}}{\rm{erg/s}},$ (43) $\displaystyle\frac{\rm
SFR}{M_{\odot}/\rm{yr}}$ $\displaystyle=1.4\times
10^{-41}\frac{L_{[O\,II]}}{\rm{erg/s}},$ (44) $\displaystyle\frac{\rm
SFR}{M_{\odot}/\rm{yr}}$ $\displaystyle=7.6\times
10^{-42}\frac{L_{[O\,III]}}{\rm{erg/s}}.$ (45)

For H$\beta$, we assume a fixed line ratio of $H\beta/H\alpha=0.35$
(Osterbrock & Ferland, 2006). As Ly$\alpha$ emission is not discussed in Gong
et al. (2017), we use a model with a constant halo mass to luminosity
M–$L_{{\rm Ly}\alpha}$ ratio at halo masses of $M=10^{8}$–$10^{15}$ Mpc/h, and
calibrate the scaling factor to match the mean Ly$\alpha$ intensity of the
analytical model from Pullen et al. (2014, Eq. 51).

The line emission bias factor, $b_{I,\rm line}$, is given by the the halo
bias, $b(M,z)$, weighted by line luminosity:

$b_{I,\rm line}(z)=\frac{\int_{M_{\rm min}}^{M_{\rm
th}(z)}dM\frac{dn}{dM}(M,z)b(M,z)L_{\rm line}(M,z)}{\int_{M_{\rm min}}^{M_{\rm
th}(z)}dM\frac{dn}{dM}(M,z)L_{\rm line}(M,z)}.$ (46)

We use the same lower and upper integration limits ($M_{\rm min}$ and $M_{\rm
th}(z)$) as Eq. 18.

###  3.2 Quasar Spectrum

Our quasar emission model assumes a quasar luminosity function and an averaged
quasar spectrum. We adopt the fiducial model of the quasar luminosity function
in Conroy & White (2013), and use the quasar spectrum template compiled in
Mas-Ribas & Hennawi (2018) that splices together a composite spectra from
Vanden Berk et al. (2001), Richards et al. (2006), and Lusso et al. (2015).

###  3.3 Foreground Stars

Galactic stars contribute to the Poisson noise (Eq. 25) in the auto spectrum
of the intensity field. We model star counts, $dN_{\rm star}/dm/d\Omega$ in
Eq. 25, using the stellar population synthesis code,
Trilegal222http://stev.oapd.inaf.it/cgi-bin/trilegal (Girardi et al., 2002,
2005). We use their model in the UBVRIJHKLMN photometric system at the north
galactic pole ($\ell=0^{\circ}$, $b=90^{\circ}$).

##  4 Cross-Correlation Forecast

###  4.1 Evaluation Metrics

We focus on constraining the redshift-dependent IGL spectrum, $d(\nu
I_{\nu})/dz(\nu,z)$, from cross-correlating galaxies and the intensity field.
The clustering and Poisson-noise amplitudes in the cross power spectrum are
proportional to $b_{I}(z)d(\nu I_{\nu})/dz(z,\nu)$ and $\left.d(\nu
I_{\nu})/dz\right|_{\rm g}(z,\nu)$, respectively (Eq. 28, 32). Therefore we
forecast the signal-to-noise ratio (SNR) on these two quantities obtained from
the cross power spectrum. Note $\left.d(\nu I_{\nu})/dz\right|_{\rm g}(z,\nu)$
represents the IGL spectrum from the specific galaxies (tracers) used in the
cross-correlation measurements, and does not equal the cosmic mean IGL
spectrum, $d(\nu I_{\nu})/dz(z,\nu)$.

As the matter power spectrum is linear on large scales, here defined as
$k\lesssim 0.2$ h/Mpc, we evaluate constraints on $b_{I}(z)d(\nu
I_{\nu})/dz(z,\nu)$ from the amplitude of all $\ell_{\rm clus}^{\rm
min}<\ell<\ell_{\rm clus}^{\rm max}$ modes, where $\ell_{\rm clus}^{\rm
max}=0.2{\rm[h/Mpc]}\chi(z)-1/2$, and the minimum accessible mode $\ell_{\rm
clus}^{\rm min}=51$, which corresponds to SPHEREx’s field of view of 3.5
degree on the smaller side, as we conservatively assumes modes larger than the
field size will be susceptible to the zodiacal light foreground filtering
process and other large-scale systematcis and become inaccessible. The SNR on
the clustering amplitude is

${\rm SNR}_{\rm clus}(b_{I}(z)\frac{d(\nu
I_{\nu})}{dz}(z,\nu))=\sqrt{\sum_{\ell\in\left[\ell_{\rm clus}^{\rm
min},\ell_{\rm clus}^{\rm max}\right]}\left[\frac{C^{\rm x}_{\ell,\rm
clus}}{\delta C_{\ell}^{\rm x}}\right]^{2}}.$ (47)

In the Poisson regime, defined as $k\gtrsim 1$ h/Mpc, we use $\ell$ modes in
the range of $\ell_{\rm P}^{\rm min}<\ell<\ell_{\rm P}^{\rm max}$ to constrain
$\left.\frac{d(\nu I_{\nu})}{dz}\right|_{\rm g}(z,\nu)$, where $\ell_{\rm
P}^{\rm min}=1{\rm[h/Mpc]}\chi(z)-1/2$, and $\ell_{\rm P}^{\rm max}=10^{5}$,
which is the highest multipole mode available given SPHEREx pixel size of
$6.2^{\prime\prime}$. Incorporating the effects from both Gaussian and non-
Guassian terms, the SNR on the Poisson amplitude is given by

$\begin{split}&{\rm SNR}_{\rm P}(\left.\frac{d(\nu
I_{\nu})}{dz}\right|_{g}(z,\nu))\\\ &=\frac{C^{\rm x}_{\ell,\rm
P}}{\sqrt{\frac{1}{N_{\ell>\ell_{\rm min}}}\sum_{\ell\in\left[\ell_{\rm
P}^{\rm min},\ell_{\rm P}^{\rm max}\right]}\left(\delta C_{\ell}^{\rm
x}\right)^{2}+\frac{\mathcal{T}^{gIgI}}{\Omega_{\rm sur}}}},\end{split}$ (48)

where $\Omega_{\rm sur}$ is the survey area solid angle used for cross
correlation.

Since $b_{I}(z)d(\nu I_{\nu})/dz(z,\nu)$ and $\left.d(\nu
I_{\nu})/dz\right|_{g}(z,\nu)$ are the main quantities of interest, we
forecast below their constraints as a function of redshift and wavelength from
upcoming surveys.

###  4.2 Masking

In practice, pixels that contain bright stars and bright galaxies will have to
be masked in data processing to reduce sample variance noise in the cross
power spectrum, and the resulting power spectrum amplitude will depend on the
depth of the masking limit. To model the effect of masking on the power
spectrum, we assume a photometric source catalog with sufficient depth
provides the location and fluxes of bright stars and galaxies to be masked.
The masking magnitude threshold, $m_{\rm th}^{\rm ph}$, is set by the limiting
magnitude of a photometric band at frequency $\nu^{\rm ph}$ (or the
corresponding wavelength $\lambda^{\rm ph}$).

For simplicity, we mask all stars and galaxies brighter than our masking
threshold (see Sec. 8.2 for discussion). Since only the flux density
information is needed, we can use a photometric source catalog for this
purpose, which is usually deeper than a spectroscopic catalog. Given $m_{\rm
th}^{\rm ph}$ and $\nu^{\rm ph}$, we adopt the abundance matching concept
(Kravtsov et al., 2004; Vale & Ostriker, 2004; Conroy et al., 2006; Behroozi
et al., 2010; Moster et al., 2010) to infer the masking magnitude threshold
$m_{\rm th}(\nu,z)$ for continuum emission at a given frequency $\nu$; this
assumes that the brightest $N$ sources in a given frequency bin are also the
brightest ones in another frequency bin. We therefore first calculate the
number density of galaxies being masked at $\nu^{\rm ph}$:

$N=\int^{m_{\rm th}^{\rm ph}}_{-\infty}dm\Phi(m,z,\nu^{\rm ph}(1+z)),$ (49)

and then at each frequency solve for the $m_{\rm th}(\nu,z)$ value which gives

$N=\int^{m_{\rm th}(\nu,z)}_{-\infty}dm\Phi(m,z,\nu(1+z)).$ (50)

Analogously, for line emission modeled with a line luminosity and halo mass
relation, given the halo mass function we can derive its maximum halo mass
$M_{\rm th}(z)$ by solving for

$N=\int^{\infty}_{M_{\rm th}(z)}dM\frac{dn}{dM}(M,z).$ (51)

Throughout this work, we set the fiducial masking threshold to $m_{\rm
th}^{\rm ph}=20$ at frequency $\nu^{\rm ph}=3\times 10^{14}$ Hz (1 $\mu$m).
Our model suggests that $\sim 3\%$ of the SPHEREx pixel ($6.2^{\prime\prime}$)
contains stars or galaxies above this limit, and the intensity Poisson noise
from unmasked stars is approximately an order of magnitude below the galaxy
Poisson noise. For simplicity we ignore the sensitivity loss from mode mixing
due to masking, and discuss the implications in Sec. 8.2.

###  4.3 Fiducial Case

We consider cross power spectra between SPHEREx intensity map and several
current and future galaxy redshift surveys. For simplicity, we present results
with a few chosen parameters as our fiducial case.

Fig. 1 shows the redshift-dependent IGL spectra, $d(\nu I_{\nu})/dz$, from our
model (Sec. 3) at fiducial redshifts defined below. The spectra include
continuum and spectral line emission, and are presented with SPHEREx spectral
binning (Sec. 5) of 96 channels over $0.75-5\mu$m333SPHEREx recently modified
its design to have 102 frequency channels covering the same wavelength range.
As this has negligible impact on our results, we use their 96 channel
configuration with published sensitivity forecasts for this work.. The IGL
spectral amplitude builds up with cosmic time. Since SPHEREx has higher
spectral resolution at $\lambda>3.8$ $\mu$m (see Sec. 5), the line intensity
at longer wavelength is stronger compared to the continuum as the line signal
is inversely proportional to the spectral width (Eq. 19). The model we used
for IGL continuum (Helgason et al., 2012) is only calibrated to observations
at $z\lesssim 5$ and rest frame wavelengths $\lambda>0.15$ $\mu$m. The higher
redshifts and shorter wavelengths models are thus derived from extrapolation.

At low redshifts ($z<3$), cross correlations between SPHEREx and spectroscopic
redshift surveys with high redshift accuracy ($\sigma_{z}/(1+z)\sim
10^{-4}-10^{-3}$) allow selecting galaxies in thin redshift slices. We choose
a redshift bin width of $\Delta z^{g}/(1+z^{g})=0.03$, and calculate the
cross-correlation at five fiducial redshifts, $z^{g}=[0.25,0.5,1,2,3]$.
Between a redshift of $0$ and $3$ there are $\sim 46$ such measurements. We
choose a fiducial masking threshold $m_{\rm th}^{\rm ph}=20$ at $\lambda^{\rm
ph}=1$ $\mu$m for all these cases, except for at $z=0.25$ where we use $m_{\rm
th}^{\rm ph}=18$. This threshold is chosen to optimize the trade off between
losing signal and reducing foreground emission (see Sec. 8.2 for more
discussion). SPHEREx will achieve an all-sky point source sensitivity of
$m_{\rm AB}\sim 20$ at $\lambda/\Delta\lambda\sim 4$, and therefore our
fiducial masking depth is feasible with the SPHEREx source catalog.

At higher redshifts ($z>3$), spectroscopic catalogs are either unavailable or
lacking sufficient number density and/or sky coverage, and therefore we
consider cross-correlating SPHEREx with photometric redshift catalogs, which
have a redshift accuracy of $\sigma_{z}/(1+z)\sim 0.02-0.05$. We choose eight
fiducial redshifts for the calculation: $z^{g}=[3,4,5,6,7,8,9,10]$. At
$3\leqslant z\leqslant 6$, we use a redshift bin width of $\Delta
z^{g}/(1+z^{g})=0.1$, which results in $\sim 6$ such measurements. At $z>6$,
the redshift bins are determined by the expected redshift resolution from the
Lyman-break dropout technique (see Sec. 6.2 for details). We also set the
masking threshold to $m_{\rm th}^{\rm ph}=20$ at $\lambda^{\rm ph}=1$ $\mu$m.

Figure 1:  Redshift-dependent IGL spectra from our model (Sec. 3) at fiducial
redshifts, using SPHEREx spectral binning. The spectra contain continuum and
line emission. The spectral lines from left to right are Ly$\alpha$, [O ii],
H$\beta$, [O iii], and H$\alpha$, respectively, as labeled in the $z=6$ case
(in green). At $z<7$, H$\beta$ and [O iii] lines are in adjacent spectral
channel and thus appear as one blended line feature; whereas at $z\geq 7$, the
two lines redshifted into SPHEREx’s high spectral resolution bands (see Fig.
2), which split the lines by several channels, and thus two peaks are
distinguishable in the plot. Note here we present our model on infinitesimally
thin redshift slices, and ignore the intrinsic linewidth and broadening
effects of all spectral lines, so the linewidth shown only reflects the impact
of SPHEREx spectral resolution.

##  5 SPHEREx Intensity Mapping

SPHEREx is a NASA MIDEX mission scheduled to launch in 2024 (Doré et al.,
2014, 2018a)444http://spherex.caltech.edu. SPHEREx will carry out the first
all-sky near-infrared spectro-imaging survey from 0.75 to 5 $\mu$m. It has a
pixel size of $6.2^{\prime\prime}$, and a wavelength-dependent spectral
resolution: a spectral resolving power of $R=41$ at wavelengths between 0.75
and 2.42 $\mu$m, (48 spectral channels), $R=35$ between 2.42 and 3.82 $\mu$m
(16 spectral channels), $R=110$ between 3.82 and 4.42 $\mu$m (16 spectral
channels), and $R=130$ between 4.42 and 5.00 $\mu$m (16 spectral channels),
The top panel of Fig. 2 shows the SPHEREx spectral resolution as a function of
wavelength (spectral channel), and the bottom panel is the expected noise rms
per spectral channel per pixel for all-sky and the two $\sim 100$ deg2 deep
surveys555Data downloaded from: https://github.com/SPHEREx/Public-
products/blob/master/Surface_Brightness_v28_base_cbe.txt, which we use to
calculate the noise contribution to cross-power spectrum (Eq. 26). Here we
adopt SPHEREx’s optimistic sensitivity model, noting that the pessimistic
forecast is $\sim 0.5$ mag lower. For SPHEREx all-sky survey, we only use 75%
of the the full sky coverage to avoid the Galactic contamination.

Figure 2: Top: SPHEREx spectral resolution of each wavelength channel. The
vertical dashed lines mark the six wavelength bands, each contains 16
channels. Bottom: SPHEREx surface brightness sensitivity per spectral channel
for all sky (blue) and deep field (orange) in a $6.2^{\prime\prime}$ sky
pixel.

##  6 Redshift Surveys

We consider cross-correlating intensity maps from SPHEREx and several current
and future galaxy surveys. Table 1 summarizes the survey parameters of all
galaxy surveys. The model for source number density as a function of redshift
is shown in Fig. 3.

Table 1: Summary of the sky area and redshift coverage of the surveys used in this work. survey name |  $A_{\rm sur}$ [deg2] | $f_{\rm sky}$ | redshift range
---|---|---|---
Intensity Maps |  |  |
SPHEREx | 30,940 | 0.75 | N/A
SPHEREx deep | 2$\times$100 | $4.85\times 10^{-3}$ | N/A
Spectroscopic Surveys |  |  |
BOSS CMASS | 10,000 | 0.242 | $0<z<1$
eBOSS LRG | 7,500 | 0.182 | $0<z<1.2$
eBOSS ELG | 620 | 0.0150 | $0<z<1.5$
eBOSS QSO | 7,500 | 0.182 | $0<z<3.5$
DESI BGS | 14,000 | 0.339 | $0<z<0.5$
DESI LRG | 14,000 | 0.339 | $0.6<z<1.2$
DESI ELG | 14,000 | 0.339 | $0.6<z<1.7$
DESI QSO | 14,000 | 0.339 | $0.6<z<4.2$
Euclid ELG | 15,000 | 0.364 | $0.65<z<2.05$
Euclid deep ELG | 10 | $2.42\times 10^{-4}$ | $0.65<z<2.05$
Roman ELG | 2,200 | 0.053 | $0.5<z<3$
SPHEREx spec. | 30,940 | 0.75 | $0<z<4.6$
Photometric Surveys |  |  |
Rubin phot. | 18,000 | 0.436 | $0<z<4$
Rubin phot. |  |  |
(gold sample) | 18,000 | 0.436 | $0<z<3$
Rubin Lyman break | 18,000 | 0.436 | $2.4<z<3.7$
| 18,000 | 0.436 | $3.7<z<4.8$
| 18,000 | 0.436 | $4.8<z<5.9$
| 18,000 | 0.436 | $5.9<z<6.7$
Roman phot. | 2,200 | 0.053 | $0<z<6$
Roman Lyman Break | 2,200 | 0.053 | $5.5<z<6.5$
| 2,200 | 0.053 | $6.5<z<7.5$
| 2,200 | 0.053 | $7.5<z<8.5$
| 2,200 | 0.053 | $8.5<z<9.5$
| 2,200 | 0.053 | $9.5<z<10.5$
Figure 3:  Model of number density of the spectroscopic and photometric galaxy
catalogs from each survey considered in this work. The values are compiled
from the literature or predicted with our model, detailed in Sec. 6. We show
the number density per co-moving volume on the top panel, and the number
density per redshift and solid angle on the bottom panel. The vertical dashed
lines mark the fiducial redshifts used in this work.

###  6.1 Spectroscopic Galaxy Surveys

We consider the following current and planned future spectroscopic redshift
surveys: SDSS BOSS/eBOSS, DESI, Euclid, the High Latitude Survey with the
Roman Space Telescope. Each survey has one or multiple types of density
tracers. Below we describe the parameters we use for each tracer catalog,
including the number density, sky coverage, and the galaxy bias factors.

####  6.1.1 SDSS BOSS/eBOSS

SDSS BOSS (Dawson et al., 2013) and eBOSS (Dawson et al., 2016) are
spectroscopic surveys targeting several types of galaxy tracers. In this work,
we use four spectroscopic samples from the SDSS: the BOSS CMASS survey, eBOSS
LRG (luminous red galaxies), eBOSS ELG (emission line galaxies), and eBOSS QSO
(quasars) surveys.

For BOSS CMASS, eBOSS LRG, and eBOSS QSO, we use the predicted number density
from Dawson et al. (2016) Table 1 ($z_{\rm conf}>1$ for LRG, and QSO_CORE New
+ Known for QSO); for eBOSS ELG, the number density is given in Table 4 of
Raichoor et al. (2017). We assume a tophat functional form of $dN/dV(z)$, the
number density per co-moving volume, between the redshift bins quoted in their
tables. Our model for galaxy bias, $b_{g}$, also follows the prescription of
Dawson et al. (2016) and Raichoor et al. (2017): $b_{g}({\rm
LRG})=1.7/\sigma_{8}(0)/\sigma_{8}(z)$, $b_{g}({\rm
ELG})=1.0/\sigma_{8}(0)/\sigma_{8}(z)$, $b_{g}({\rm QSO})=0.53+0.29(1+z)^{2}$,
and we set $b_{g}({\rm CMASS})=2$ according to Nuza et al. (2013).

####  6.1.2 DESI

The Dark Energy Spectroscopic Instrument (DESI) (DESI Collaboration et al.,
2016) is an ongoing spectroscopic survey that observes four types of tracers
across 14,000 deg2: bright galaxy samples (BGS), luminous red galaxies (LRG),
star-forming emission line galaxies (ELG), and quasars (QSO). The predicted
number density for each tracer catalog is given in Table 2.3, 2.5, and 2.7 of
(DESI Collaboration et al., 2016), and we assume a tophat functional form of
$dN/dV(z)$ between the redshift bins quoted in their tables. For LRG, ELG, and
QSO, we use the same bias model as the eBOSS: $b_{g}({\rm
LRG})=1.7/\sigma_{8}(0)/\sigma_{8}(z)$, $b_{g}({\rm
ELG})=1.0/\sigma_{8}(0)/\sigma_{8}(z)$, $b_{g}({\rm QSO})=0.53+0.29(1+z)^{2}$,
and we set $b_{g}({\rm BGS})=2$, same as the BOSS CMASS case.

####  6.1.3 Euclid

The Euclid spectroscopic survey (Laureijs et al., 2011) will map ELGs
(H$\alpha$ emitters) over 15,000 deg2 between redshifts $0.9\lesssim z\lesssim
2$. The expected number density as a function of redshift is given by Table 3
of Amendola et al. (2018) (we use the $n_{2}$ reference case), and we assume a
tophat $dN/dV(z)$ between the redshift bins quoted in their table. For galaxy
bias, we use the model from Merson et al. (2019): $b_{g}(z)=0.7z+0.7$.

In addition, Euclid will also have deep fields that overlap with the SPHEREx
deep field at the north ecliptic pole. We therefore consider a case of cross-
correlating the Euclid deep field catalog with SPHEREx deep field images. For
Euclid deep field sources, we adopt the same bias and redshift coverage as the
wide-field ELGs, and assume a constant number density of $n_{g}=4\times
10^{-3}$ (h/Mpc)3 over this redshift range.

####  6.1.4 Roman Space Telescope

The Roman Space Telescope (Spergel et al., 2015) plans a spectroscopic survey
of ELGs using the H$\alpha$ and [O iii] lines as part of the High Latitude
Survey (HLS; Doré et al., 2018b). Compared to Euclid, the Roman Space
Telescope-HLS covers a smaller sky area (2200 deg2) but with better point
source sensitivity. We adopt the ELG galaxy number density predictions from
Zhai et al. (2019). We use the H$\alpha$ number density (Table 1) for $z<2$,
the [O iii] number density (Table 2) for $2<z<3$, and assume a tophat
$dN/dV(z)$ between the redshift bins quoted in their tables. For both spectral
lines we use their “dust model fit at high redshifts” and $1\times 10^{-16}$
erg s-1 cm-2 flux limit case.

####  6.1.5 SPHEREx Spectroscopic Catalog

In addition to diffuse imaging, SPHEREx will also produce a spectroscopic
galaxy catalog (Doré et al., 2014, 2018a). We therefore also consider the case
of cross-correlating the SPHEREx galaxy catalog with its intensity maps. For
the SPHEREx all-sky survey, we use the predicted number density and galaxy
bias of the SPHEREx source catalog at $z<4.6$ from their public products
666https://github.com/SPHEREx/Public-
products/blob/master/galaxy_density_v28_base_cbe.txt. Specifically, we use the
forecast with galaxy photometric redshift error of
$0.01<\sigma_{z}/(1+z)<0.03$, and consider a tophat functional form of
$dN/dV(z)$ between the quoted redshift bins. The SPHEREx deep field with
higher sensitivity will produce a deeper source catalog, and we leave its
application to future work.

###  6.2 Photometric Galaxy Surveys

We consider cross-correlating SPHEREx and two upcoming photometric surveys:
the Rubin Observatory Legacy Survey of Space and Time (LSST) and the Roman
Space Telescope High Latitude Survey (HLS). Despite the larger redshift
uncertainties in photometric samples, they have the advantage of reaching
higher redshifts and fainter sources.

####  6.2.1 Rubin Observatory LSST

The Rubin Observatory LSST will carry out an 18,000 deg2 photometric galaxy
survey at $0<z<4$ in optical bands (LSST Science Collaboration et al., 2009).
The expected 5-$\sigma$ point source depth after a ten-year observation is a
magnitude of 27.5 at r band777https://www.lsst.org/scientists/keynumbers . A
subset of galaxies from $0<z<3$ with higher redshift accuracy
($\sigma_{z}/(1+z)<0.05$) constitute the “gold sample,” which has a depth of
25.3 in i band (LSST Science Collaboration et al., 2009). We predict the
number density of these two samples by applying these magnitude thresholds to
our galaxy luminosity function model from Helgason et al. (2012).

Photometric redshift estimates are not available for $z>4$ sources from the
Rubin Observatory LSST; nevertheless, redshift information of high-z sources
can be determined by the Lyman-break drop-out technique. We thus estimate
Lyman-break galaxy samples in a few broad redshift bands (listed in Table 1)
defined by the LSST filter boundaries and the redshifted Ly$\alpha$ wavelength
(rest-frame $121.6$ nm). For these Lyman-break selected galaxies, we estimate
their number density by applying the same magnitude threshold (27.5 in r band)
to our luminosity function model. We use Eq. 42 to estimate galaxy bias in
both the photometric and Lyman-break samples.

####  6.2.2 The Roman Space Telescope

The planned photometric galaxy survey as part of the Roman Space Telescope HLS
(Doré et al., 2018b) covers 2200 deg2 with a 5-$\sigma$ depth of $\sim 26.5$
in Y, J, and H band (Spergel et al., 2015). We apply these magnitude
thresholds to our IGL luminosity function model to predict the number density
of the HLS photometric galaxies from $z=0$ to $z=6$. AT $z>6$, we use the
expected Lyman-break selected galaxy number density from the Roman Space
Telescope Science Investigation Team EXPO (S. Finkelstein, private
communication)888See slides from the presentation given by S. Finkelstein in
the conference Astronomy in the 2020s: Synergies with WFIRST, available at the
time of writing at
https://www.stsci.edu/~dlaw/WFIRST2020s/slides/finkelstein.pdf. These number
density values are consistent with the numerical simulation prediction from
Waters et al. (2016) at $z\geqslant 8$. We also use Eq. 42 to estimate galaxy
biases in both photometric and high-redshift samples.

##  7 Results

###  7.1 An Example Cross Power Spectrum

Figure 4:  Examples of cross power spectra (black) of SPHEREx 3 $\mu$m image
and SPHEREx spectroscopic galaxy catalog at $z=0.25$, $1$, and $3$. Blue and
orange lines are the clustering and Poisson-noise components, respectively. We
constrain the redshift evolution of the combination of bias and redshift-
dependent intensity of the EBL and of tracer galaxies, $b_{I}(z)d(\nu
I_{\nu})/dz(z,\nu)$ and $\left.d(\nu I_{\nu})/dz\right|_{g}(z,\nu)$,
respectively, from the cross-power spectrum amplitudes in the clustering (blue
band) and Poisson regimes (orange band), respectively; their SNR values are
given in the legend. See Sec. 4.1 for details.

In Fig. 4, we show a set of example cross power spectra of a all-sky SPHEREx
intensity map at 3 $\mu$m and SPHEREx spectroscopic catalog at $z=0.25$, $1$,
and $3$, respectively. The cross power spectra are given by Eq. 27 with error
bars estimated by Eq. 4. In the Poisson-noise regime, we expect high
significance measurements of the cross power spectra at all three redshifts
due to the large number of modes on small scales, while in the clustering
regime, we expect significant detection at the lower two redshifts only. For
the $z=3$ case, the low tracer number density gives a high galaxy Poisson
noise level that reduces sensitivity on the cross power spectrum (see Sec. 8.1
for discussions).

We aim to extract two quantities related to the redshift-dependent IGL
spectra, $b_{I}(z)d(\nu I_{\nu})/dz(z,\nu)$ and $\left.d(\nu
I_{\nu})/dz\right|_{g}(z,\nu)$, from the large (clustering) and small
(Poisson-noise) scale amplitudes of the cross power spectrum. Their respective
SNRs are denoted as ${\rm SNR}_{\rm clus}$ and ${\rm SNR}_{\rm P}$,
respectively. Note that for ${\rm SNR}_{\rm P}$, we also include the non-
Gaussian effect from the trispectrum term (Eq. 48).

###  7.2 Sensitivity to the Clustering Term

Figure 5:  Constraints on the amplitude of cross power spectra of SPHEREx and
galaxy surveys on large scales at $z\leqslant 3$. Top: $b_{I}(z)d(\nu
I_{\nu})/dz(z,\nu)$ from our model using spectroscopic (black) and photometric
(grey, shifted down by $\times 0.5$ for presentation purposes) galaxy samples.
The photometric surveys have wider redshift bins and thus lower spectral
resolution. Middle: SNR forecast on $b_{I}(z)d(\nu I_{\nu})/dz(z,\nu)$ (${\rm
SNR}_{\rm clus}$ in Eq. 47) per SPHEREx spectral channel from cross-
correlating SPHEREx images with different spectroscopic surveys. For the
surveys that have multiple tracer catalogs, we plot the one that gives the
highest SNR values (BOSS CMASS and DESI BGS for $z=0.25$; BOSS CMASS for
$z=0.5$; eBOSS QSO and DESI ELG for $z=1$; eBOSS QSO and DESI QSO for $z=2$
and $3$). Solid lines are the cases using SPHEREx all-sky survey, and dashed
lines are the cases with two SPHEREx deep fields, which have lower instrument
noise but smaller sky coverage. Bottom: same as the middle rows but with
photometric surveys.

The top row of Fig. 5 shows our model for $b_{I}(z)d(\nu I_{\nu})/dz(z,\nu)$
as a function of the 96 SPHEREx spectral bins and at the five fiducial
$z\leqslant 3$ redshifts and masking thresholds. This quantity can be derived
from the amplitude of the cross power spectra on large scales. We consider
both spectroscopic and photometric galaxy surveys, where we use different
redshift bin widths ($\Delta z^{g}/(1+z^{g})=$ 0.03 and 0.1 for spectroscopic
and photometric, respectively) to account for the respective redshift
accuracies. For cross-correlation with photometric galaxies, the SPHEREx
spectral line features are smoothed out due to the wide redshift bins.

The middle and bottom rows of Fig. 5 show the SNR forecasts on $b_{I}(z)d(\nu
I_{\nu})/dz(z,\nu)$, i.e., ${\rm SNR}_{\rm clus}$ in Eq. 47, per SPHEREx
spectral channel by cross-correlating SPHEREx maps with different
spectroscopic (middle) and photometric (bottom) galaxy catalogs. For surveys
that have multiple tracer catalogs, we plot the ones that give the highest SNR
values. With spectroscopic galaxies, we expect to achieve a ${\rm SNR}_{\rm
clus}>5$ out to $z\sim 2$. However, at $z=3$, due to the low number density of
tracers, ${\rm SNR}_{\rm clus}\sim 1$ to 3. Photometric galaxies, on the other
hand, have a much higher source density, and can reach a ${\rm SNR}_{\rm
clus}$ of $\sim 10(5)$ with Rubin observatory LSST (Roman Space Telescope HLS)
at $z=3$. Dashed lines indicate results from SPHEREx deep fields, which have
lower instrument noise but smaller sky coverage. According to our model, the
SNR in the deep fields are much lower than the all-sky cases, due to their
small $f_{\rm sky}$ that results in a large sample variance noise
contribution. Note that with our chosen redshift binning there are $\sim 37$
measurements at $z<2$ and $\sim 9$ measurements at $2<z<3$.

Figure 6:  Constraints on the amplitude of cross power spectra of SPHEREx and
galaxy surveys on large scales at $4<z<10$. Top: $b_{I}(z)d(\nu
I_{\nu})/dz(z,\nu)$ from our model using photometric surveys. We cut off the
spectrum below the Lyman-$\alpha$ (121.6 nm) for $z>6$ since those high-energy
photons will be absorbed by neutral gas in the intergalactic medium. Bottom:
SNR forecast on $b_{I}(z)d(\nu I_{\nu})/dz(z,\nu)$ (${\rm SNR}_{\rm clus}$ in
Eq. 47) per SPHEREx spectral channel from cross-correlating SPHEREx images
with different photometric surveys. The Lyman break selected galaxy bins from
Rubin Observatory LSST and the Roman Space Telescope are also shown in the
panels close to their central redshifts. Solid lines are the cases using
SPHEREx all-sky survey, and dashed lines are the cases with two SPHEREx deep
fields.

Fig. 6 shows the same forecast as Fig. 5 at high redshifts ($4<z<10$) with
photometric surveys. Note with our chosen redshift binning there are $\sim 6$
measurements available between $3<z<6$. For $z>6$ where the universe has not
yet fully ionized, we cut off the spectrum below the Lyman-$\alpha$ wavelength
(121.6 nm rest frame) since those high-energy photons will be absorbed by the
neutral gas in the intergalactic medium.

The Lyman break selected galaxy bins from Rubin Observatory LSST and Roman
Space Telescope are also shown in the panels close to their central redshift.
Note that the Lyman break galaxy samples have a much wider redshift bin, and
therefore they have higher sensitivity but lower spectral resolution on the
inferred IGL spectrum.

Our estimate suggests that the IGL spectrum can be robustly measured over the
full $0.75-5\mu$m wavelength range across redshift, to a significance level of
a few $\sigma$ per SPHEREx spectral channel out to $z=6$ by cross-correlating
SPHEREx all-sky and upcoming photometric galaxy surveys.

###  7.3 Sensitivity to the Poisson-noise Term

Figure 7:  Constraints on the cross power spectrum Poisson noise amplitude at
$z\leqslant 3$. Top: the IGL spectrum from all galaxies, $d(\nu I_{\nu})/dz$
(black), and from the subset of the galaxies used for cross-correlations,
$\left.d\nu I_{\nu}/dz\right|_{g}$ (colored), which is proportional to the
Poisson noise amplitude. Solid lines give the cases using the SPHEREx all-sky
survey, and dashed lines are the cases with the two SPHEREx deep fields, which
have lower instrument noise but smaller sky coverage. Bottom: SNR forecast on
$\left.d(\nu I_{\nu})/dz\right|_{g}$ (${\rm SNR}_{\rm P}$ in Eq. 48) per
SPHEREx spectral channel from cross-correlating SPHEREx images with different
spectroscopic surveys. For surveys that have multiple tracer catalogs, we use
the same tracer as in Fig. 5, which gives the highest SNR values. We use the
fiducial masking depths as described in Sec. 4.3. In the CMASS case at
$z=0.25$ and $z=0.5$, our chosen masking threshold is deeper than the galaxy
catalog ($m_{\rm th}>m_{g}$), and thus the cross Poisson power vanishes (see
Eq. 37). Therefore we only have clustering and no Poisson-noise measurement
for these two cases.

Fig. 7 shows our forecast on the cross power spectrum Poisson noise amplitude
of SPHEREx and galaxy surveys at $z\leqslant 3$. The top row compares the
redshift-depedent IGL spectrum from all sources, $d(\nu I_{\nu})/dz$, and from
the subset of sources in the galaxy samples used for cross-correlation,
$\left.d(\nu I_{\nu})/dz\right|_{g}$, which is proportional to the Poisson
noise amplitude. Naturally, deeper catalogs capture a larger fraction of the
total IGL in the cross Poisson noise amplitude. The bottom row shows our SNR
forecast on $\left.d(\nu I_{\nu})/dz\right|_{g}$, i.e. ${\rm SNR}_{\rm P}$, in
Eq. 48. In most cases, the cross-correlation results in a highly significant
${\rm SNR}_{\rm P}$ ($\sim 10^{2}-10^{4}$) because of the large number of
modes available on small scales. Dashed lines indicate the cross-correlation
of the two SPHEREx deep field images with a Euclid deep field. Compared to
all-sky, the deep fields capture a higher fraction of the IGL, as shown in the
top row, but the SNR is lower due to the small sky coverage.

Figure 8:  Constraints on the cross power spectrum Poisson noise amplitude for
$4<z<10$. Top: the IGL spectrum from all galaxies, $d(\nu I_{\nu})/dz$
(black), and from the subset of the galaxies used for cross-correlations,
$\left.d\nu I_{\nu}/dz\right|_{g}$ (colored), which is proportional to the
Poisson noise amplitude. Solid lines give the cases using the SPHEREx all-sky
survey, and dashed lines are the cases with the two SPHEREx deep fields, which
have lower instrument noise but smaller sky coverage. The Lyman break selected
galaxy bins from Rubin Observatory LSST and Roman Space Telescope are also
shown in the panels close to their central redshifts. We cut off the spectrum
below the Lyman-$\alpha$ wavelength (121.6 nm) for $z>6$ since those high
energy photons will be absorbed by neutral gas in the intergalactic medium.
Bottom: SNR forecast on $\left.d(\nu I_{\nu})/dz\right|_{g}$ (${\rm SNR}_{\rm
sh}$ in Eq. 48) per SPHEREx spectral channel from cross-correlations. The kink
in SNR above $\lambda\sim 3.8$ $\mu$m is due to the change in SPHEREx
resolving power and thus instrument noise per channel (see Fig. 2).

Fig. 8 shows the same forecast as Fig. 7 at high redshifts ($4<z<10$) with
photometric surveys and the Lyman-break selected samples. We observe a slight
kink in SNR at $\lambda\sim 3.8$ $\mu$m due to the change in SPHEREx spectral
resolution and thus instrument noise level (see Fig. 2). In practice, one can
bin the measurements at $\lambda>3.8$ $\mu$m spectrally to increase the
sensitivity. The same feature is not evident at the low redshift cases as
shown in Fig. 5, because instrumental noise only becomes the dominant source
of uncertainty in $\delta C_{\ell}^{x}$ (Eq. 4) at higher redshifts. Our
forecast suggests that high significance measurements with SNR $\sim 100$, can
be achieved out to $z=6$ with upcoming photometric galaxy surveys.

Note that $\left.d(\nu I_{\nu})/dz\right|_{g}$ is the mean spectrum of sources
in the galaxy catalog, which is not the same as the mean IGL spectrum from the
averaged emission from all sources. Nevertheless, for simplicity, we assume
any subset of galaxies has the same mean spectrum as the IGL model used in
this calculation. We also highlight potentially different measurements from
the cross Poisson noise term with an example of quasars in Sec. 9.1.

##  8 Discussion

###  8.1 Error Estimation

Figure 9:  Cross power spectra (black solid lines) and their error components
(colored lines) from cross-correlating SPHEREx all-sky map at the 3-$\mu$m
channel with SPHEREx spectroscopic catalogs at $z=0.25$ (left), $1$ (middle),
and $3$ (right), using a masking threshold of $m_{\rm th}^{\rm ph}=18$. Grey
solid and dashed lines are the clustering and Poisson-noise terms in the cross
power spectrum, respectively. The thick black dashed lines show the total
Gaussian noise power spectrum variance (i.e., the sum of all colored lines) as
a function of $\ell$ bin, and the horizontal bars on the signal cross power
spectra (short black solid lines) mark the $\ell$ binning. The thick blue and
orange horizontal bars are the total Gaussian variance on the cross-power
spectrum clustering and Poisson-noise amplitudes. They are computed by combing
the multipole modes within the $\ell$ range of $[\ell^{\rm min}_{\rm
clus},\ell^{\rm max}_{\rm clus}]$ (blue band) and $[\ell^{\rm min}_{\rm
P},\ell^{\rm max}_{\rm P}]$ (orange band), respectively. The amplitude of the
non-Gaussian variance from trispectrum is shown by the green dashed horizontal
bar (see Eq. 47 and 48). The two lower redshift cases are background-limited
where the cross spectrum errors are dominated by $C_{\ell,\rm
clus}^{I}C_{\ell,\rm clus}^{g}/N_{\ell}$, whereas at $z=3$, the sensitivity is
limited by the galaxy Poisson noise due to its low tracer number density.

The variance of a cross power spectrum is the sum of several contributing
factors in Eq. 4, including instrument noise, galaxy Poisson (shot) noise, and
sample variance. Fig. 9 compares three cross power spectra and their different
error components from cross-correlating SPHEREx all-sky maps at the 3-$\mu$m
channel, with SPHEREx spectroscopic catalogs at $z=0.25$, $1$, and $3$, using
a masking threshold of $m_{\rm th}^{\rm ph}=18$. The Gaussian error on the
cross power spectrum (Eq. 4) comprises a quadratic sum of individual power
spectrum terms:

$\begin{split}\left(\delta C_{\ell}^{\rm
x}\right)^{2}=\frac{1}{N_{\ell}}&\left[(C_{\ell}^{\rm
x})^{2}+C_{\ell}^{I}C_{\ell}^{g}\right]\\\
=\frac{1}{N_{\ell}}&\left[\right.(C_{\ell}^{\rm x})^{2}+\left(C_{\ell,\rm
clus}^{I}+C_{\ell,\rm P}^{I}+C_{\ell,n}^{I}\right)\\\ &\cdot\left(C_{\ell,\rm
clus}^{g}+C_{\ell,\rm P}^{g}\right)\left.\right]\\\
=\frac{1}{N_{\ell}}&\left[\right.(C_{\ell}^{\rm x})^{2}+C_{\ell,\rm
clus}^{I}C_{\ell,\rm clus}^{g}+C_{\ell,\rm clus}^{I}C_{\ell,\rm P}^{g}\\\
&+C_{\ell,\rm P}^{I}C_{\ell,\rm clus}^{g}+C_{\ell,\rm P}^{I}C_{\ell,\rm
P}^{g}\\\ &+C_{\ell,n}^{I}C_{\ell,\rm clus}^{g}+C_{\ell,n}^{I}C_{\ell,\rm
P}^{g}\left.\right].\end{split}$ (52)

First, we see that the cross-spectrum amplitude decreases with redshift
following the redshift dependence of the IGL spectrum (Fig. 1) and the
underlying matter density fluctuations. For the two lower redshift cases, the
dominant noise term is $C_{\ell,\rm clus}^{I}C_{\ell,\rm clus}^{g}/N_{\ell}$.
While $C_{\ell,\rm clus}^{g}$ is determined by the underlying cosmology, the
only way to further bring down the error on cross-power spectrum is by
reducing the amplitude of $C_{\ell,\rm clus}^{I}$, which is dominated by
aggregate emission along the line of sight. Therefore, these measurements are
in what we refer to as background-limited regime, where background here refers
to all emissions not from the redshift of interest (and thus include
foreground, too). We can apply a deeper mask to suppress the background
(foreground). However, this will also reduce the signal, and therefore there
exists a masking threshold that optimizes the trade-off between background
(foreground) and signal removal. We detail this comparison in Sec. 8.2.

At $z=3$, the dominant noise term becomes $C_{\ell,\rm clus}^{I}C_{\ell,\rm
P}^{g}/N_{\ell}$ due to the increased galaxy Poisson noise power spectrum
amplitude, $C_{\ell,\rm P}^{g}$, in available surveys. This means that a
deeper catalog with a higher tracer number density can effectively improve the
sensitivity.

Figure 10:  Cross power spectra (black solid lines) and their error components
(colored lines) from cross-correlating SPHEREx all-sky map at the 3-$\mu$m
channel with Euclid ELG (left) and DESI QSO (middle), and SPHEREx deep-field
map at 3-$\mu$m with the Euclid deep field ELG catalog (right). All three
cases are at $z=1$ with a masking threshold of $m_{\rm th}^{\rm ph}=18$. The
curves are the same as in Fig. 9. The Euclid ELG (left) and DESI QSO (middle)
are comparable in survey area size, but the former has lower tracer number
density and thus higher Poisson noise that suppresses the sensitivity. The
Euclid deep field case (right) has a higher ELG number density than the wide
field case (left), however it has a lower SNR due to the larger sample
variance from its smaller sky coverage.

Fig. 10 compares three different galaxy surveys at $z=1$ with the same masking
threshold of $m_{\rm th}^{\rm ph}=18$ for the SPHEREx 3-$\mu$m channel. The
first two, Euclid ELG and DESI QSO surveys, have a similar survey area,
$f_{\rm sky}$, and thus similar $N_{\ell}$’s, but Euclid ELG has a higher
source number density. As a result, the two cases have similar cross spectrum
amplitudes in the clustering regime, but DESI QSO has a higher Poisson noise
amplitude, and the noise terms that contain $C_{\ell,\rm P}^{g}$ (all dashed
lines) are therefore higher in the DESI QSO case.

We can also compare the Euclid wide and deep field cross-correlations with the
SPHEREx all-sky and deep fields, respectively. As shown on the left and right
panels of Fig. 10, we see that the terms associated with galaxy Poisson noise
(all dashed lines) are lower in the deep fields because of the higher galaxy
number density. However, the deep fields cannot access the lower $\ell$ modes,
and the overall noise is higher in the deep field due to the larger sample
variance (small $N_{\ell}$) from its smaller sky coverage.

Figure 11:  Cross power spectra (black solid lines) and their error components
(colored lines) from cross-correlating SPHEREx all-sky maps at the 3-$\mu$m
channel with DESI BGS at $z=0.25$ using a masking threshold of $m_{\rm
th}^{\rm ph}=14$ (left), $18$ (middle), and $22$ (right), respectively. The
curves are the same as in Fig. 9. The SNR is improved by raising the masking
threshold from $m_{\rm th}^{\rm ph}=14$ to $m_{\rm th}^{\rm ph}=18$, since the
line-of-sight contamination (the $C_{\ell,\rm clus}^{\rm I}$ and $C_{\ell,\rm
P}^{\rm I}$ terms) decreases, while signal reduction is negligible (black
solid lines). However, if we further increase the masking threshold to $m_{\rm
th}^{\rm ph}=22$, the SNR becomes lower due to the significant loss of signal
by masking.

Fig. 11 demonstrates the impact of masking depth, using the cross-correlation
of SPHEREx all-sky maps at 3 $\mu$m and DESI BGS galaxies at $z=0.25$ with
three different masking depths, $m_{\rm th}^{\rm ph}=14,18,22$. Comparing
$m_{\rm th}^{\rm ph}=14$ and $m_{\rm th}^{\rm ph}=18$, the signals are almost
the same (black solid line) but the total Gaussian error (black dashed line)
is lower with a deeper mask; this suggests that by masking sources to
$14<m_{\rm th}^{\rm ph}<18$, a significant portion of the foreground can be
removed, while the signal reduction is negligible which leads to an
improvement of SNR. However, if we further increase the masking depth to
$m_{\rm th}^{\rm ph}=22$, the signal is drastically reduced, which means a
significant fraction of IGL at this redshift is from galaxies with $18<m_{\rm
th}^{\rm ph}<22$, and we cannot get a higher SNR by this deeper masking depth.

###  8.2 Masking Depth

Figure 12:  Top: the IGL spectrum $d(\nu I_{\nu})/dz$ in the SPHEREx 3-$\mu$m
channel as a function of masking depth $m_{\rm th}^{\rm ph}$ at fiducial
masking wavelength $\lambda^{\rm ph}=1$ $\mu$m. Bottom: ${\rm SNR}_{\rm clus}$
as a function of $m_{\rm th}^{\rm ph}$ with different galaxy surveys shown in
Fig. 5.

We consider masking all sources from an external catalog that are brighter
than a flux threshold, and therefore the choice of masking depth is a trade-
off between signal loss and foreground reduction. A deeper masking depth tends
to remove more foreground sources contributing to the intensity maps, but at
the same time, bright galaxies in the target redshift range of interest can
also be masked. Fig. 12 demonstrates how the IGL amplitude and the error
depends on masking depth. Here we consider masking sources brighter than
magnitude threshold $m_{\rm th}^{\rm ph}$ at the fiducial masking wavelength
$\lambda^{\rm ph}=1$ $\mu$m. At lower redshift, the IGL amplitude strongly
depends on the masking depths $m_{\rm th}^{\rm ph}$. At higher redshift, the
signals are almost unchanged with varying masking depth, as the sources that
dominate the IGL at these redshifts are below the masking thresholds. In other
words, at low redshift, there is a masking depth that optimizes the tradeoff
between signal loss and foreground reduction. At higher redshift, deeper masks
perform better since they reduce more foreground emission while having
negligible impact on the signal.

Here we only consider masking all stars and galaxies brighter than a magnitude
threshold. In practice, if the external catalog provides redshift information
of individual sources, a redshift-dependent magnitude threshold might also
help to minimize foreground contamination from lower-redshift clustering
signals. We leave this investigation to future work.

We ignore the effect of pixel loss due to masking in our forecast. In reality,
masking results in highly non-uniform images, which not only reduces the
number of multipole modes in the data, but also introduces mode mixing, the
coupling of signal on different scales that needs to be corrected for in
analysis. The reduction in power spectral sensitivity due to masking as a
function of multipole depends on the source luminosity function and clustering
amplitude, the masking size around each source, and the mode-coupling
correction method. A more detailed study requires realistic simulated images,
which is beyond the scope of this work.

###  8.3 Caveats of the Abundance Matching Model

In our model, we use the abundance matching technique to relate the underlying
halo mass distribution to the luminosity of galaxy tracer continuum and line
emissions. In other words, if we have N galaxy tracers, we assume they are in
the most massive N dark matter halos, and are also the strongest emitters of
continuum and spectral lines in all frequencies. In reality, the relation has
a large scatter. It has been shown that red and blue galaxies have different
spatial clustering amplitudes and trace the matter density field differently
(Zehavi et al., 2011). Each galaxy has a distinct color and spectral line
strength depending on its stellar composition, metallicity, dust attenuation,
and other intersetllar medium properties. Moreover, the tracer galaxies are
typically selected by certain spectral features, making them unrepresentative
of the average EBL emitters. For example, the ELG samples will bias towards
star-forming galaxies which have strong emission lines. The de-correlation
between matter and continuum and line emissions will not only introduce
uncertainties to our model, but also reduce the cross correlation signal on
small scales. As discussed in Chiang et al. (2019), the stochasticity in the
relation between mass, tracers, and EBL emission could be a possible
explanation to their low EBL intensity bias values. Quantifying these de-
correlations requires a coherent model for mass and light, which we leave to
future work.

###  8.4 Intensity-bias Degeneracy

We derive constraints on $b_{I}(z)d(\nu I_{\nu})/dz(z,\nu)$ from the cross
power spectrum amplitude in the clustering regime, but we do not further
consider decoupling these two terms. In practice, we can apply the methods
used in Chiang et al. (2019) to jointly fit the redshift and frequency
dependencies of bias and intensity, using all the cross spectrum measurements.
In Chiang et al. (2019), they calibrate the total intensity to the intensity
of resolved sources observed in the local universe, as the sources are
sufficiently faint to account for most of the intensity from all emitters. In
our case, we can cross-correlate SPHEREx maps with a deep catalog at low
redshift, and similarly assume $\left.d(\nu I_{\nu})/dz\right|_{g}\approx
d(\nu I_{\nu})/dz$ to break the bias-intensity degeneracy. According to our
model, SPHEREx spectroscopic catalog at $z=0.1$ contains $\sim 95\%$ of the
total IGL intensity, and thus this approach is feasible for our
interpretation, assuming little redshift evolution of the EBL composition.

##  9 Science Interpretation

###  9.1 Poisson Noise Constraints on the Tracer Spectrum

In the Poisson-noise sensitivity forecast presented in Sec. 7.3, we assume the
tracer population has the same averaged spectrum as our IGL model spectrum.
This is not necessarily true if we select a certain type of source for cross
correlation (e.g., ELG, QSO), where the Poisson noise signal will be the
averaged spectrum of the tracer sources. Therefore we can extract potentially
distinct information from the Poisson-noise spectrum. As a demonstration, we
consider a case of cross-correlating SPHEREx with quasars from DESI at $z=2$,
using a quasar spectrum and luminosity function model described in Sec. 3.2.
The constraints on the averaged quasar spectrum, $\left.d(\nu
I_{\nu})/dz\right|_{Q}$, from the cross Poisson spectrum is shown in Fig. 13.
Using the small-scale cross-correlation information of SPHEREx with the quasar
sample, we can obtain a high-SNR measurement on the average quasar spectrum at
near-infrared wavelengths across redshift.

Given the high SNR value from the Poisson-noise measurements in all cases
considered in this work, it is also feasible to gain more information from
Poisson-noise by splitting the tracers into sub-samples by different
properties. For example, we can select galaxies by their stellar mass, star
formation rate, morphology, and constrain their average spectrum individually.

Figure 13:  Average quasar’s spectrum constraints (per SPHEREx spectral bin)
from cross-correlating SPHEREx with DESI quasar samples at $z=2$. The blue
line is the modeled quasar spectrum template from Mas-Ribas & Hennawi (2018).

###  9.2 Science from Non-Linear Clustering

In this work, we only consider information from the linear and Poisson noise
scales of the cross power spectrum. Non-linear clustering scales also contain
crucial information on galaxy evolution as well as the underlying cosmology.
For example, emission from satellite galaxies and the contribution of IHL can
be constrained by the amplitude and shape of the non-linear power spectrum
(e.g., Cooray et al., 2012; Zemcov et al., 2014; Cheng et al., 2021). As we
have demonstrated that a high-SNR measurement can be achieved by cross-
correlating SPHEREx with current or upcoming galaxy surveys, it is promising
to extract more constraints from non-linear scales.

###  9.3 Science from Near-Infrared EBL Tomography

The rest-frame near-infrared light from galaxies are dominated by low mass
stars, and thus the near-infrared photometry has been used to measure the
stellar mass function across redshifts (Madau & Dickinson, 2014). Most of the
current constraints on stellar mass function are based on deep photometric
samples (Muzzin et al., 2013; Wright et al., 2018; Thorne et al., 2021). While
individual galaxies can be resolved down to the $\sim 10^{-2}-10^{-3}M*$ mass
limit at lower redshift (e.g. Muzzin et al., 2013), at $z\gtrsim 3$,
photometric surveys are only complete down to $\sim M*$ scales, and the faint
end estimation are extrapolated from the bright end with a functional fit
(Thorne et al., 2021). SPHEREx EBL measurement, $d\nu I_{\nu}/dz(\nu,z)$,
provides complimentary information to the current galaxy detection-based
measurement. Although individual galaxy spectra cannot be extracted from the
diffused maps, SPHEREx EBL measurements contain emission from all sources with
high spectral and redshift resolutions, and cover a large sky area. They can
be used to calibrate the stellar mass function model built from the resolved
galaxies, and probe emission from diffuse origins such as the IHL. At higher
redshift, SPHEREx covers the rest-frame optical and ultraviolet spectra and
helps constrain the star formation history (Madau & Dickinson, 2014; Doré et
al., 2014, 2018a).

###  9.4 Science from Spectral Line Intensity Mapping

While quoting the SNR, we do not separate the significance of spectral lines
from continuum emission in our results. Our forecast suggests that spectral
lines can be detected with high significance at lower redshift ($z\lesssim 3$)
in both the clustering and Poisson-noise regimes. Line emissions contain
valuable information on the galaxy properties. For example, the H$\alpha$,
Ly$\alpha$, and [O ii] line flux can be used to infer the star formation rate
(Brinchmann et al., 2004; Sobral et al., 2013; Gunawardhana et al., 2013;
Khostovan et al., 2015); the Balmer decrement, the ratio of H$\alpha$ and
H$\beta$ line flux, is a crucial indicator of dust extinction (Domínguez et
al., 2013). The SPHEREx cross-correlation measurements of line emissions in
the intensity mapping regime have the potential to address these galaxy
properties across redshifts. Moreover, with the high spectral resolution and
sensitivity, SPEHREx offers an unique probe of the three-dimensional large-
scale structure with line intensity mapping (e.g. Kovetz et al., 2017) for
cosmological constraints in the high-redshift Universe.

###  9.5 Cosmological Constraints

So far we have focused on, EBL spectrum constraints derived from cross-
correlation, which only depend on the large-scale and Poisson noise amplitude
of the cross power spectrum. In general, abundant cosmological information can
be inferred from the power spectrum; for example, baryon acoustic oscillation
signals can be extracted from the cross power spectrum, which can then be used
to constrain the growth and geometry of the universe (Eisenstein et al.,
2005).

##  10 Conclusion

We forecast the sensitivity of near-infrared EBL tomography using spectro-
imaging from the upcoming SPHEREx mission through cross-correlation with
several current and future galaxy surveys, spanning redshifts of $0.25<z<10$.
We consider IGL as the only emitting source that constitutes the near-infrared
EBL, and build a model for IGL continuum and spectral line emissions. We then
model the cross power spectrum of SPHEREx images and galaxy surveys as a
function of SPHEREx spectral channel and galaxy redshift.

From the amplitudes of the cross power spectrum on linear scales, we infer the
redshift-dependent EBL spectrum, multiplied by a bias factor, $b_{I}(z)d(\nu
I_{\nu})/dz(z,\nu)$. Our forecast suggests that this quantity can be
constrained to a significance of at least a few $\sigma$ out to $z=6$, using
the SPHEREx all-sky survey in cross-correlation with a combination of
spectroscopic (at $z\lesssim 3$) and photometric galaxy surveys (at
$0.25<z<10$). At $3\lesssim z\lesssim 10$, photometric galaxy catalogs from
the upcoming Rubin Observatory and the Roman Space Telescope can constrain
$b_{I}(z)d(\nu I_{\nu})/dz(z,\nu)$, albeit at a lower redshift and spectral
resolution.

The Poisson noise level of the cross power spectrum represents the average
intensity of sources used in cross correlation. Our forecast shows that the
Poisson noise level can be extracted at high significance ($\gtrsim 10^{2}$)
out to $z\sim 10$, suggesting the SPHEREx all-sky and deep-field data can
provide high sensitivity measurements of the average near-infrared spectrum of
any selected population of tracer sources.

In summary, a high-sensitivity tomographic measurement of the EBL spectrum can
be achieved by cross-correlating SPHEREx with current and future galaxy
spectroscopic and photometric surveys. Making use of the clustering and
Poisson noise of the cross power spectra, we expect further cosmological and
astrophysical information can be extracted from this rich dataset.

We are grateful to Jamie Bock, Yi-Kuan Chiang, Emmanuel Schann, Olivier Doré,
Richard Feder, Lluis Mas-Ribas, Daniel Masters and Rogier Windhorst for
helpful discussions and comments on the draft. We are in debt to Lluis Mas-
Ribas for providing the quasar spectrum template, Emmanuel Schann for
suggesting the trispectrum correction, and Steven Finkelstein for providing
the Roman Space Telescope High Latitude Survey galaxy number count estimates.
Y.-T.C. acknowledges support by the Ministry of Education, Taiwan through the
Taiwan-Caltech Scholarship. T.-C.C. acknowledges support from the JPL
Strategic R&TD awards. Part of this work was done at Jet Propulsion
Laboratory, California Institute of Technology, under a contract with the
National Aeronautics and Space Administration.

## References

  * Abdalla et al. (2020) Abdalla, H., Adam, R., Aharonian, F., et al. 2020, MNRAS, 494, 5590, doi: 10.1093/mnras/staa999
  * Abdo et al. (2010) Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, ApJ, 723, 1082, doi: 10.1088/0004-637X/723/2/1082
  * Abeysekara et al. (2019) Abeysekara, A. U., Archer, A., Benbow, W., et al. 2019, ApJ, 885, 150, doi: 10.3847/1538-4357/ab4817
  * Acciari et al. (2019) Acciari, V. A., Ansoldi, S., Antonelli, L. A., et al. 2019, MNRAS, 486, 4233, doi: 10.1093/mnras/stz943
  * Ackermann et al. (2012) Ackermann, M., Ajello, M., Allafort, A., et al. 2012, Science, 338, 1190, doi: 10.1126/science.1227160
  * Ackermann et al. (2018) Ackermann, M., Ajello, M., Baldini, L., et al. 2018, ApJS, 237, 32, doi: 10.3847/1538-4365/aacdf7
  * Aharonian et al. (2006) Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2006, Nature, 440, 1018, doi: 10.1038/nature04680
  * Aharonian et al. (2007) Aharonian, F., Akhperjanian, A. G., Barres de Almeida, U., et al. 2007, A&A, 475, L9, doi: 10.1051/0004-6361:20078462
  * 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
  * Behroozi et al. (2010) Behroozi, P. S., Conroy, C., & Wechsler, R. H. 2010, ApJ, 717, 379, doi: 10.1088/0004-637X/717/1/379
  * Bernstein (2007) Bernstein, R. A. 2007, ApJ, 666, 663, doi: 10.1086/519824
  * Brinchmann et al. (2004) Brinchmann, J., Charlot, S., White, S. D. M., et al. 2004, MNRAS, 351, 1151, doi: 10.1111/j.1365-2966.2004.07881.x
  * Cappelluti et al. (2013) Cappelluti, N., Kashlinsky, A., Arendt, R. G., et al. 2013, ApJ, 769, 68, doi: 10.1088/0004-637X/769/1/68
  * Caputo et al. (2020) Caputo, A., Vittino, A., Fornengo, N., Regis, M., & Taoso, M. 2020, arXiv e-prints, arXiv:2012.09179.  https://arxiv.org/abs/2012.09179
  * Chang et al. (2010) Chang, T.-C., Pen, U.-L., Bandura, K., & Peterson, J. B. 2010, Nature, 466, 463, doi: 10.1038/nature09187
  * Cheng et al. (2021) Cheng, Y.-T., Arai, T., Bangale, P., et al. 2021, arXiv e-prints, arXiv:2103.03882.  https://arxiv.org/abs/2103.03882
  * Chiang et al. (2020) Chiang, Y.-K., Makiya, R., Ménard, B., & Komatsu, E. 2020, ApJ, 902, 56, doi: 10.3847/1538-4357/abb403
  * Chiang & Ménard (2019) Chiang, Y.-K., & Ménard, B. 2019, ApJ, 870, 120, doi: 10.3847/1538-4357/aaf4f6
  * Chiang et al. (2019) Chiang, Y.-K., Ménard, B., & Schiminovich, D. 2019, ApJ, 877, 150, doi: 10.3847/1538-4357/ab1b35
  * Conroy et al. (2006) Conroy, C., Wechsler, R. H., & Kravtsov, A. V. 2006, ApJ, 647, 201, doi: 10.1086/503602
  * Conroy & White (2013) Conroy, C., & White, M. 2013, ApJ, 762, 70, doi: 10.1088/0004-637X/762/2/70
  * Cooray (2016) Cooray, A. 2016, Royal Society Open Science, 3, 150555, doi: 10.1098/rsos.150555
  * Cooray et al. (2012) Cooray, A., Smidt, J., de Bernardis, F., et al. 2012, Nature, 490, 514, doi: 10.1038/nature11474
  * Cote et al. (2019) Cote, P., Abraham, B., Balogh, M., et al. 2019, in Canadian Long Range Plan for Astronomy and Astrophysics White Papers, Vol. 2020, 18, doi: 10.5281/zenodo.3758463
  * Creque-Sarbinowski & Kamionkowski (2018) Creque-Sarbinowski, C., & Kamionkowski, M. 2018, Phys. Rev. D, 98, 063524, doi: 10.1103/PhysRevD.98.063524
  * Croft et al. (2018) Croft, R. A. C., Miralda-Escudé, J., Zheng, Z., Blomqvist, M., & Pieri, M. 2018, MNRAS, 481, 1320, doi: 10.1093/mnras/sty2302
  * Croft et al. (2016) Croft, R. A. C., Miralda-Escudé, J., Zheng, Z., et al. 2016, MNRAS, 457, 3541, doi: 10.1093/mnras/stw204
  * Dawson et al. (2013) Dawson, K. S., Schlegel, D. J., Ahn, C. P., et al. 2013, AJ, 145, 10, doi: 10.1088/0004-6256/145/1/10
  * Dawson et al. (2016) Dawson, K. S., Kneib, J.-P., Percival, W. J., et al. 2016, AJ, 151, 44, doi: 10.3847/0004-6256/151/2/44
  * DESI Collaboration et al. (2016) DESI Collaboration, Aghamousa, A., Aguilar, J., et al. 2016, arXiv e-prints, arXiv:1611.00036.  https://arxiv.org/abs/1611.00036
  * 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
  * Domínguez et al. (2013) Domínguez, A., Siana, B., Henry, A. L., et al. 2013, ApJ, 763, 145, doi: 10.1088/0004-637X/763/2/145
  * Doré et al. (2014) Doré, O., Bock, J., Ashby, M., et al. 2014, ArXiv e-prints.  https://arxiv.org/abs/1412.4872
  * Doré et al. (2018a) Doré, O., Werner, M. W., Ashby, M. L. N., et al. 2018a, arXiv e-prints, arXiv:1805.05489.  https://arxiv.org/abs/1805.05489
  * Doré et al. (2018b) Doré, O., Hirata, C., Wang, Y., et al. 2018b, arXiv e-prints, arXiv:1804.03628.  https://arxiv.org/abs/1804.03628
  * 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
  * Eisenstein et al. (2005) Eisenstein, D. J., Zehavi, I., Hogg, D. W., et al. 2005, ApJ, 633, 560, doi: 10.1086/466512
  * Girardi et al. (2002) Girardi, L., Bertelli, G., Bressan, A., et al. 2002, A&A, 391, 195, doi: 10.1051/0004-6361:20020612
  * Girardi et al. (2005) Girardi, L., Groenewegen, M. A. T., Hatziminaoglou, E., & da Costa, L. 2005, A&A, 436, 895, doi: 10.1051/0004-6361:20042352
  * 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
  * Gong et al. (2017) Gong, Y., Cooray, A., Silva, M. B., et al. 2017, ApJ, 835, 273, doi: 10.3847/1538-4357/835/2/273
  * Gunawardhana et al. (2013) Gunawardhana, M. L. P., Hopkins, A. M., Bland-Hawthorn, J., et al. 2013, MNRAS, 433, 2764, doi: 10.1093/mnras/stt890
  * H. E. S. S. Collaboration et al. (2017) H. E. S. S. Collaboration, Abdalla, H., Abramowski, A., et al. 2017, A&A, 606, A59, doi: 10.1051/0004-6361/201731200
  * Helgason et al. (2012) Helgason, K., Ricotti, M., & Kashlinsky, A. 2012, ApJ, 752, 113, doi: 10.1088/0004-637X/752/2/113
  * Hill et al. (2016) Hill, J. C., Ferraro, S., Battaglia, N., Liu, J., & Spergel, D. N. 2016, Phys. Rev. Lett., 117, 051301, doi: 10.1103/PhysRevLett.117.051301
  * Hopkins & Beacom (2006) Hopkins, A. M., & Beacom, J. F. 2006, ApJ, 651, 142, doi: 10.1086/506610
  * 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
  * Kennicutt (1998) Kennicutt, Robert C., J. 1998, Annual Review of Astronomy and Astrophysics, 36, 189, doi: 10.1146/annurev.astro.36.1.189
  * Khostovan et al. (2015) Khostovan, A. A., Sobral, D., Mobasher, B., et al. 2015, MNRAS, 452, 3948, doi: 10.1093/mnras/stv1474
  * 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, arXiv e-prints, arXiv:2102.12323.  https://arxiv.org/abs/2102.12323
  * Kovetz et al. (2017) Kovetz, E. D., Viero, M. P., Lidz, A., et al. 2017, ArXiv e-prints.  https://arxiv.org/abs/1709.09066
  * Kravtsov et al. (2004) Kravtsov, A. V., Berlind, A. A., Wechsler, R. H., et al. 2004, ApJ, 609, 35, doi: 10.1086/420959
  * 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
  * Laureijs et al. (2011) Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, arXiv e-prints, arXiv:1110.3193.  https://arxiv.org/abs/1110.3193
  * Levenson et al. (2007) Levenson, L. R., Wright, E. L., & Johnson, B. D. 2007, ApJ, 666, 34, doi: 10.1086/520112
  * LSST Science Collaboration et al. (2009) LSST Science Collaboration, Abell, P. A., Allison, J., et al. 2009, arXiv e-prints, arXiv:0912.0201.  https://arxiv.org/abs/0912.0201
  * Lusso et al. (2015) Lusso, E., Worseck, G., Hennawi, J. F., et al. 2015, MNRAS, 449, 4204, doi: 10.1093/mnras/stv516
  * Ly et al. (2007) Ly, C., Malkan, M. A., Kashikawa, N., et al. 2007, ApJ, 657, 738, doi: 10.1086/510828
  * Madau & Dickinson (2014) Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415, doi: 10.1146/annurev-astro-081811-125615
  * MAGIC Collaboration et al. (2008) MAGIC Collaboration, Albert, J., Aliu, E., et al. 2008, Science, 320, 1752, doi: 10.1126/science.1157087
  * Mas-Ribas & Hennawi (2018) Mas-Ribas, L., & Hennawi, J. F. 2018, AJ, 156, 66, doi: 10.3847/1538-3881/aace5f
  * Masui et al. (2013) Masui, K. W., Switzer, E. R., Banavar, N., et al. 2013, ApJ, 763, L20, doi: 10.1088/2041-8205/763/1/L20
  * 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
  * McQuinn & White (2013) McQuinn, M., & White, M. 2013, MNRAS, 433, 2857, doi: 10.1093/mnras/stt914
  * Ménard et al. (2013) Ménard, B., Scranton, R., Schmidt, S., et al. 2013, arXiv e-prints, arXiv:1303.4722.  https://arxiv.org/abs/1303.4722
  * Merson et al. (2019) Merson, A., Smith, A., Benson, A., Wang, Y., & Baugh, C. 2019, MNRAS, 486, 5737, doi: 10.1093/mnras/stz1204
  * Mitchell-Wynne et al. (2016) Mitchell-Wynne, K., Cooray, A., Xue, Y., et al. 2016, ApJ, 832, 104, doi: 10.3847/0004-637X/832/2/104
  * Mitchell-Wynne et al. (2015) Mitchell-Wynne, K., Cooray, A., Gong, Y., et al. 2015, Nature Communications, 6, 7945, doi: 10.1038/ncomms8945
  * Moster et al. (2010) Moster, B. P., Somerville, R. S., Maulbetsch, C., et al. 2010, ApJ, 710, 903, doi: 10.1088/0004-637X/710/2/903
  * Muzzin et al. (2013) Muzzin, A., Marchesini, D., Stefanon, M., et al. 2013, ApJ, 777, 18, doi: 10.1088/0004-637X/777/1/18
  * Newman (2008) Newman, J. A. 2008, ApJ, 684, 88, doi: 10.1086/589982
  * Nuza et al. (2013) Nuza, S. E., Sánchez, A. G., Prada, F., et al. 2013, MNRAS, 432, 743, doi: 10.1093/mnras/stt513
  * Osterbrock & Ferland (2006) Osterbrock, D. E., & Ferland, G. J. 2006, Astrophysics of gaseous nebulae and active galactic nuclei
  * 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
  * Pullen et al. (2013) Pullen, A. R., Chang, T.-C., Doré, O., & Lidz, A. 2013, ApJ, 768, 15, doi: 10.1088/0004-637X/768/1/15
  * Pullen et al. (2014) Pullen, A. R., Doré, O., & Bock, J. 2014, ApJ, 786, 111, doi: 10.1088/0004-637X/786/2/111
  * Pullen et al. (2018) Pullen, A. R., Serra, P., Chang, T.-C., Doré, O., & Ho, S. 2018, MNRAS, 478, 1911, doi: 10.1093/mnras/sty1243
  * Raichoor et al. (2017) Raichoor, A., Comparat, J., Delubac, T., et al. 2017, MNRAS, 471, 3955, doi: 10.1093/mnras/stx1790
  * Richards et al. (2006) Richards, G. T., Lacy, M., Storrie-Lombardi, L. J., et al. 2006, ApJS, 166, 470, doi: 10.1086/506525
  * 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
  * Schaan et al. (2018) Schaan, E., Ferraro, S., & Spergel, D. N. 2018, Phys. Rev. D, 97, 123539, doi: 10.1103/PhysRevD.97.123539
  * Schechter (1976) Schechter, P. 1976, ApJ, 203, 297, doi: 10.1086/154079
  * Schmidt et al. (2015) Schmidt, S. J., Ménard, B., Scranton, R., et al. 2015, MNRAS, 446, 2696, doi: 10.1093/mnras/stu2275
  * Scott et al. (2021) Scott, B., Upton Sanderbeck, P., & Bird, S. 2021, arXiv e-prints, arXiv:2104.00017.  https://arxiv.org/abs/2104.00017
  * 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
  * Serra et al. (2014) Serra, P., Lagache, G., Doré, O., Pullen, A., & White, M. 2014, A&A, 570, A98, doi: 10.1051/0004-6361/201423958
  * Sheth et al. (2001) Sheth, R. K., Mo, H. J., & Tormen, G. 2001, MNRAS, 323, 1, doi: 10.1046/j.1365-8711.2001.04006.x
  * Sobral et al. (2013) Sobral, D., Smail, I., Best, P. N., et al. 2013, MNRAS, 428, 1128, doi: 10.1093/mnras/sts096
  * 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
  * Thompson et al. (2007) Thompson, R. I., Eisenstein, D., Fan, X., Rieke, M., & Kennicutt, R. C. 2007, ApJ, 666, 658, doi: 10.1086/520634
  * Thorne et al. (2021) Thorne, J. E., Robotham, A. S. G., Davies, L. J. M., et al. 2021, MNRAS, doi: 10.1093/mnras/stab1294
  * 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
  * Vale & Ostriker (2004) Vale, A., & Ostriker, J. P. 2004, MNRAS, 353, 189, doi: 10.1111/j.1365-2966.2004.08059.x
  * Vanden Berk et al. (2001) Vanden Berk, D. E., Richards, G. T., Bauer, A., et al. 2001, AJ, 122, 549, doi: 10.1086/321167
  * Waters et al. (2016) Waters, D., Di Matteo, T., Feng, Y., Wilkins, S. M., & Croft, R. A. C. 2016, MNRAS, 463, 3520, doi: 10.1093/mnras/stw2000
  * 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
  * Wright et al. (2018) Wright, A. H., Driver, S. P., & Robotham, A. S. G. 2018, MNRAS, 480, 3491, doi: 10.1093/mnras/sty2136
  * Yang et al. (2019) Yang, S., Pullen, A. R., & Switzer, E. R. 2019, MNRAS, 489, L53, doi: 10.1093/mnrasl/slz126
  * 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
  * Zehavi et al. (2011) Zehavi, I., Zheng, Z., Weinberg, D. H., et al. 2011, ApJ, 736, 59, doi: 10.1088/0004-637X/736/1/59
  * 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
  * Zhai et al. (2019) Zhai, Z., Benson, A., Wang, Y., Yepes, G., & Chuang, C.-H. 2019, MNRAS, 490, 3667, doi: 10.1093/mnras/stz2844
  * Zheng et al. (2005) Zheng, Z., Berlind, A. A., Weinberg, D. H., et al. 2005, ApJ, 633, 791, doi: 10.1086/466510