# Probing the Diffuse Lyman-alpha Emission on Cosmological Scales: Ly$\alpha$
Emission Intensity Mapping Using the Complete SDSS-IV eBOSS Survey

Xiaojing Lin<EMAIL_ADDRESS><EMAIL_ADDRESS>Department of Astronomy, Tsinghua University, Beijing 100084, China  Zheng
Zheng<EMAIL_ADDRESS>Department of Physics and Astronomy,
University of Utah, UT 84112, USA  Zheng Cai<EMAIL_ADDRESS>Department of Astronomy, Tsinghua University, Beijing 100084, China

(Received xx, xx; Revised xx, xx; Accepted xx)

###### Abstract

Based on the Sloan Digital Sky Survey Data Release 16, we have detected the
large-scale structure of Ly$\alpha$ emission in the Universe at redshifts
$z=2$–3.5 by cross-correlating quasar positions and Ly$\alpha$ emission
imprinted in the residual spectra of luminous red galaxies. We apply an
analytical model to fit the corresponding Ly$\alpha$ surface brightness
profile and multipoles of the redshift-space quasar-Ly$\alpha$ emission cross-
correlation function. The model suggests an average cosmic Ly$\alpha$
luminosity density of $6.6_{-3.1}^{+3.3}\times 10^{40}{\rm
erg\,s^{-1}cMpc^{-3}}$, a $\sim 2\sigma$ detection with a median value about
8–9 times those estimated from deep narrowband surveys of Ly$\alpha$ emitters
at similar redshifts. Although the low signal-to-noise ratio prevents us from
a significant detection of the Ly$\alpha$ forest-Ly$\alpha$ emission cross-
correlation, the measurement is consistent with the prediction of our best-fit
model from quasar-Ly$\alpha$ emission cross-correlation within current
uncertainties. We rule out the scenario that these Ly$\alpha$ photons mainly
originate from quasars. We find that Ly$\alpha$ emission from star-forming
galaxies, including contributions from that concentrated around the galaxy
centers and that in the diffuse Ly$\alpha$ emitting halos, is able to explain
the bulk of the the Ly$\alpha$ luminosity density inferred from our
measurements. Ongoing and future surveys can further improve the measurements
and advance our understanding of the cosmic Ly$\alpha$ emission field.

intensity mapping – diffuse Lyman alpha emission – large scale structure

††journal: ApJS

##  1 Introduction

The filamentary structure of the cosmic web, which links galaxies to the
intergalactic medium (IGM), is predicted to be a rich reservoir of nearly
pristine gas (e.g., Fumagalli et al., 2011; Giavalisco et al., 2011).
Reprocessed radiation from quasars or the ultraviolet (UV) background will
ionize hydrogen atoms in the circumgalactic medium (CGM) and IGM (e.g.,
Gallego et al., 2021; Borisova et al., 2016; Lujan Niemeyer et al., 2022), and
the recombination of the ionized hydrogen will produce fluorescent Ly$\alpha$
emission, especially in the high-redshift Universe (Cantalupo et al., 2008; Li
et al., 2021). Extended Ly$\alpha$ emission is expected due to the high cross-
section of Ly$\alpha$ photons for resonant scatterings by neutral hydrogen
(Zheng et al., 2011a).

Direct imaging of the IGM Ly$\alpha$ emission is challenging because of its
low surface brightness (SB)(Cantalupo et al., 2005). One solution is to search
around local ionized sources, such as luminous quasars, which reside at the
densest regions of the cosmic web. The diffuse gas emission can be enhanced by
orders of magnitude, leading to the discovery of the enormous Ly$\alpha$
nebulae (ELANe) (Cantalupo et al., 2014; Hennawi et al., 2015; Cai et al.,
2017, 2018; Arrigoni Battaia et al., 2018). These extrema of Ly$\alpha$
nebulosities have Ly$\alpha$ surface brightness $\geq$ $10^{-17}{\rm
erg\,s^{-1}cm^{-2}arcsec^{-2}}$ and Ly$\alpha$ luminosity $\geq$ $10^{44}{\rm
erg\,s^{-1}}$, with Ly$\alpha$ size greater than 200 kpc. Recently, the
progress in wide-field integral field spectrographs extends the detectability
of CGM/IGM with low surface brightness. The most advanced facilities, such as
the Keck Cosmic Web Imager (KCWI, Morrissey et al., 2018) and the MultiUnit
Spectroscopic Explorer (MUSE, Bacon et al., 2010), can reach a surface
brightness of a few $\times 10^{-19}{\rm erg\,s^{-1}cm^{-2}arcsec^{-2}}$,
making it possible for observational probes into emission from the CGM/IGM in
the vicinity of bright sources (KCWI: Borisova et al. 2016; Arrigoni Battaia
et al. 2016; Cai et al. 2019, etc.; MUSE: Wisotzki et al. 2018; Bacon et al.
2021; Kusakabe et al. 2022, etc.). Large amounts of individual Ly$\alpha$
halos around strong Ly$\alpha$ emitters (LAEs), have been detected thanks to
these state-of-the-art instruments (e.g., Wisotzki et al., 2016; Leclercq et
al., 2017). Moreover, a recent discovery unveiled that star-forming galaxies
generally have Ly$\alpha$ halos by investigating Ly$\alpha$ emission around
UV-selected galaxies (Kusakabe et al., 2022).

On scales up to several Mpc from the central bright sources, no direct
observational evidence for diffuse gas emissions is found so far. The
predicted Ly$\alpha$ surface brightness at $z\geq 3$ stimulated by the diffuse
ionizing background, is on the order of $10^{-20}$ erg s-1 cm-2 arcsec-2
(Gould & Weinberg, 1996; Cantalupo et al., 2005; Kollmeier et al., 2010;
Witstok et al., 2019). Currently, this goes far beyond the capability of the
most advanced instruments on individual detections. The technique of line
intensity mapping (Kovetz et al., 2017) is expected to exceed current
observational limits, by mapping large scale structures with integrated
emission from spectral lines originating from galaxies and the diffuse IGM,
but without resolving discrete objects. Its application on 21-cm H I emission
has revealed a promising prospect for observing the low-density cosmic web
(Masui et al., 2013; Anderson et al., 2018; Tramonte et al., 2019; Tramonte &
Ma, 2020).

Ly$\alpha$ line can also be used for intensity mapping. Ly$\alpha$ intensity
mapping experiments can provide viable complimentary approaches to testing
many theoretical predictions on the diffuse emission from IGM filaments (Silva
et al., 2013, 2016; Heneka et al., 2017; Elias et al., 2020), bringing new
insights into the evolution of the Universe independent of cosmological
hydrodynamic simulations (Gallego et al., 2018, 2021; Croft et al., 2016,
2018). Croft et al. (2016) measured the large-scale structure of Ly$\alpha$
emission by the cross-correlation between Ly$\alpha$ surface brightness,
extracted from the spectra of luminous red galaxies (LRGs), and quasars in
Sloan Digital Sky Survey (SDSS)/Baryon Oscillation Spectroscopic Survey
(BOSS). If the Ly$\alpha$ emission originates from star formation in faint
Ly$\alpha$ emitting galaxies, the star formation rate density (SFRD) inferred
from the measurement would be $\sim$30 times higher than those from narrow-
band (NB) LAE surveys but comparable to dust corrected UV estimates, if nearly
all the Ly$\alpha$ photons from these galaxies escape without dust absorption
(Croft et al., 2016). They updated their measurements in Croft et al. (2018)
using SDSS Data Release 12 (DR12). After careful examination for possible
contaminations and systematics, the corrected cross-correlation is $\sim$50
percent lower than the DR10 result of Croft et al. (2016). They also performed
the cross-correlation of Ly$\alpha$ emission and Ly$\alpha$ forest as
complementary evidence, which presents no signal, and claimed that quasars
would dominate the Ly$\alpha$ surface brightness within 15 $h^{-1}$ cMpc.

Inspired by the cross-correlation technique in Croft et al. (2016, 2018), we
measure the Ly$\alpha$ surface brightness on scales of several Mpc from
quasars using the most up-to-date LRG spectra and quasar catalog in SDSS DR16,
much larger samples than those in Croft et al. (2018). In Section 2 we
introduce the data samples used in this work. We compute the quasar-Ly$\alpha$
emission cross-correlation and obtain the projected surface brightness profile
in Section 3. In Section 4, Ly$\alpha$ forest-Ly$\alpha$ emission cross-
correlation is carried out as a complementary measurement. In Section 5 we
perform simple analysis on our results and investigate possible Ly$\alpha$
sources for our detected signals. Our methods to remove potential
contamination are presented in Appendix A.

Throughout this paper, we adopt a spatially flat $\Lambda$ cold dark matter
($\Lambda$CDM) cosmological model according to the Planck 2018 results (Planck
Collaboration et al., 2020), with $H_{0}=100h\,{\rm km\,s^{-1}Mpc^{-1}}$ with
$h=0.674$, $\Omega_{m}$ = 0.315, $\Omega_{b}h^{2}$ = 0.0224, and
$\Omega_{c}h^{2}$ = 0.120. We use pMpc (physical Mpc) or pkpc (physical kpc)
to denote physical distances and cMpc (comoving Mpc) to denote comoving Mpc.

##  2 Data Samples

The data used in this study are selected from the final eBOSS data in the SDSS
Data Release 16 (DR16; Ahumada et al. 2020), the fourth data release of the
fourth phase of the Sloan Digital Sky Survey (SDSS-IV), which contains SDSS
observations through August 2018. As the largest volume survey of the Universe
up to date, the eBOSS survey is designed to study the expansion and structure
growth history of the Universe and constrain the nature of dark energy by
spectroscopic observation of galaxies and quasars. The spectrograph for SDSS-
IV eBOSS covers a wavelength range of 3650-10,400Å, with a resolution of
$\lambda/\Delta\lambda\sim$ 1500 at 3800Å and $\sim$2500 at 9000Å. There are
1000 fibers per 7-square-degree plate, and each fiber has a diameter of 120
$\mu$m, i.e., 2 arcsec in angle. There are two spectrographs, with each
collecting data from 500 fibers, roughly 450 dedicated to science targets and
50 for flux calibration and sky-background subtraction. The eBOSS data from
SDSS DR16 also include spectra obtained using the SDSS-I/II spectrographs
covering 3800Å to 9100Å.

In this work, we correlate the residual flux in the galaxy spectra (after
subtracting bestfit galaxy spectral templates) with quasars and Ly$\alpha$
forest to extract information of high-redshift Ly$\alpha$ emission imprinted
in the galaxy fiber spectra. We describe the quasar catalog, the LRG spectra,
and the Ly$\alpha$ forest samples used in this work.

###  2.1 Quasar Catalog

The SDSS DR16 quasar catalog (DR16Q; Lyke et al. 2020a), the largest selection
of spectroscopically confirmed quasars to date, contains 750,414 quasars in
total, including 225,082 new quasars observed for the first time. DR16Q
includes different redshift estimates generated by different methods, such as
the SDSS spectroscopic pipeline, visual inspection and principal component
analysis (PCA). It also provides a “primary” redshift for each quasar, which
is selected from, most preferably, the visual inspection redshift, or,
alternatively, the SDSS automated pipeline redshift. In this work we adopt the
“primary” redshift and apply a redshift restriction of $2.0\leq z<3.5$. This
redshift cut is also adopted in Croft et al. (2016, 2018) due to the
spectrograph cutoff for low-redshift Ly$\alpha$ emission and the limited
number of observed quasars at higher redshifts. Further, we exclude quasars
with redshift estimates of “catastrophic failures”, if their PCA-based
redshift estimates have a velocity difference of $|\Delta v|>3000$ km s-1 from
the “primary” redshift. We end up with 255,570 quasars in total, with a median
redshift of 2.40.

###  2.2 LRG Spectra

For one of the main projects of the SDSS surveys, a large sample of LRGs have
been observed spectroscopically to detect the baryon acoustic oscillations
(BAO) feature. BOSS has been conducted during 2009–2014, producing two
principal galaxy samples, LOWZ and CMASS (Reid et al., 2015). The BOSS LOWZ
galaxy sample targets the 343,160 low-redshift galaxy population spanning
redshifts $0.15<z<0.43$, to extend the SDSS-I/II Cut I LRG sample (Eisenstein
et al., 2001) by selecting galaxies of dimmer luminosity. The BOSS CMASS
galaxy sample targets 862,735 higher-redshift ($0.43<z<0.75$) galaxies. It
used similar color-magnitude cuts to those utilized by the Cut-II LRGs from
SDSS-I/II and the LRGs in 2SLAQ (Cannon et al., 2006), but with the galaxy
selection towards the bluer and fainter galaxies. Operated over 2014–2019, the
eBOSS LRG sample (Ahumada et al., 2020) extends the high-redshift tail of the
BOSS galaxies, with 298,762 LRGs covering a redshift range of $0.6<z<1.0$.

We select 1,389,712 LRG spectra from the combination of BOSS LOWZ sample, BOSS
CMASS sample and eBOSS LRG sample. These LRG spectra have been wavelength-
calibrated, sky-subtracted, flux-calibrated, and are the co-added ones of at
least three individual exposures, with a uniform logarithmic wavelength grid
spacing of $\Delta\log_{10}\lambda=10^{-4}$ (about 69 ${\rm km\,s^{-1}}$ per
pixel). Each spectrum has an inverse variance per pixel to estimate the
uncertainty, which incorporates photon noise, CCD read noise, and sky-
subtraction error. Bad pixels are flagged by pixel mask information, and we
use AND_MASK provided by SDSS to rule out bad pixels in all exposures.

Each LRG spectrum has a best-fitting model spectrum by performing a rest-frame
PCA using four eigenspectra as the basis (Bolton et al., 2012). A set of trial
redshifts are explored by shifting the galaxy eigenbasis and modelling their
minimum-chi-squared linear combination. A quadratic polynomial is added to fit
some low-order calibration uncertainties, such as the Galactic extinction,
intrinsic extinction and residual spectrophotometric calibration errors. For
each fiber, any objects along the corresponding line of sight that falls
within the fiber aperture can have their emission imprinted in the spectrum.
For example, the LRG fiber may capture the signal of diffuse Ly$\alpha$
emission originated from high-redshift galaxies and intergalactic medium, and
this is the signal we intend to extract in this work.

In the following analysis we only use the pixels from 3647Å to 5470Å in the
observed frame, corresponding to Ly$\alpha$ emission in the redshift range
$2.0<z<3.5$.

###  2.3 Ly$\alpha$ Forest

The Ly$\alpha$ forest
samples111https://data.sdss.org/sas/dr16/eboss/lya/Delta_LYA/ used in this
work are selected from the “Ly$\alpha$ regions”, $\lambda_{\rm
RF}\in[1040,1200]$Å, of 210,005 BOSS/eBOSS quasar spectra ranging from $z=2.1$
to $z=4$ (du Mas des Bourboux et al., 2020), where $\lambda_{\rm RF}$
represents the wavelength in quasar’s rest frame. Broad absorption line
quasars (BAL QSOs), bad observations, and spectra whose Ly$\alpha$ regions
have less than 50 pixels are all excluded. Then every three original pipeline
spectral pixels ($\Delta\log_{10}\lambda\sim 10^{-4}$) are rebinned
($\Delta\log_{10}\lambda\sim 3\times 10^{-4}$) for the purpose of measuring
Ly$\alpha$ correlations.

For each spectral region the flux-transmission fields are estimated by the
ratio of the observed flux, $f_{q}$, to the mean expected flux, $\langle
F_{q}\rangle$ (du Mas des Bourboux et al., 2020):

$\delta_{f}(\lambda)=\frac{f_{q}}{\langle F_{q}\rangle}-1.$ (1)

The pipeline deals with Ly$\alpha$ forest with identified damped Ly$\alpha$
systems (DLAs) cautiously. Pixels where a DLA reduces the transmission by more
than 20% are masked, and the absorption in the wings is corrected using a
Voigt profile following the procedure described in Noterdaeme et al. (2012).
Besides, we also mask $\pm 50$Å regions around the DLA positions predicted by
Ho et al. (2021), to ensure that DLA contamination is removed. The number of
the remaining Ly$\alpha$ forest pixels is $\sim 3.4\times 10^{7}$, with a
median redshift of 2.41.

##  3 Quasar-Ly$\alpha$ emission Cross-correlation

As the SDSS fiber would capture signals from high-redshift background sources,
the LRG residual spectra, with the bestfit galaxy model spectra subtracted,
may have Ly$\alpha$ emission from the high-redshift galaxies and IGM/CGM
superposed. However, the signals are overwhelmed by noises in most cases.
Cross-correlating the residual spectrum pixels with quasar positions is
equivalent to stacking the Ly$\alpha$ signal in the quasar neighborhood.
Suppressing the noise, the cross-correlation technique makes it possible to
exceed current observation limits and detect diffuse Ly$\alpha$ emission with
dimmer luminosities (Croft et al., 2016, 2018).

In this section we perform and analyze the quasar-Ly$\alpha$ cross-correlation
using the quasar catalog and LRG spectra mentioned in Section 2. In Section
3.1 we describe the detailed measurement of the two-dimensional cross-
correlation as a function of the separations along and perpendicular to the
line-of-sight direction. We measure the corresponding projected surface
brightness profile in Section 3.2 and multipoles of the redshift-space two-
point correlation function in Section 3.3.

###  3.1 Cross-correlation Transverse and Parallel to the Line of Sight

Firstly we split the LRGs into 885 subsamples based on their angular
positions, identified by the HEALPix (Górski et al., 2005) number with
Nside=16, which makes it convenient to search for neighboring quasars within a
limited sky region. After obtaining a quasar-LRG spectrum pixel pair with an
angular separation of $\theta$, we can compute the line-of-sight separation
$r_{\|}$ and transverse separation $r_{\perp}$ between these two objects:

$\displaystyle r_{\|}$ $\displaystyle=$ $\displaystyle\left[D_{\rm C}(z_{\rm
Ly\alpha})-D_{\rm C}(z_{\rm q})\right]\cos\frac{\theta}{2},$ (2)
$\displaystyle r_{\perp}$ $\displaystyle=$ $\displaystyle\left[D_{\rm
M}(z_{\rm Ly\alpha})+D_{\rm M}(z_{\rm q}\right]\sin{\frac{\theta}{2}},$ (3)

where $D_{\rm C}$ is the line-of-sight comoving distance as a function of
redshift $z$, $D_{\rm M}$ is the transverse comoving distance as a function of
redshift $z$, $z_{\rm q}$ is the quasar redshift and $z_{\rm Ly\alpha}$ is the
redshift of Ly$\alpha$ emission converted from the wavelength of the LRG
spectrum pixel, i.e., $z_{\rm Ly\alpha}=\lambda/\lambda_{\rm Ly\alpha}-1$ with
$\lambda_{\rm Ly\alpha}=1215.67$Å.

Following Croft et al. (2016), we estimate the quasar-Ly$\alpha$ emission
surface brightness cross-correlation, $\xi_{q\alpha}(r_{\perp},r_{\|})$, by
summing over all quasar-LRG spectrum pixel pairs separated by $r_{\|}$ along
the line-of-sight direction and $r_{\perp}$ along the transverse direction
within a certain bin:

$\xi_{q\alpha}(r_{\|},r_{\perp})=\frac{1}{\sum_{i=1}^{N(\vec{r})}w_{ri}}\sum_{i=1}^{N(\vec{r})}w_{ri}\Delta_{\mu,ri},$
(4)

where $N(\vec{r})$ is the number of LRG spectrum pixels within the separation
bin centered at the position $\vec{r}=(r_{\perp},r_{\|})$ and
$\Delta_{\mu,ri}=\mu_{ri}-\langle\mu(z)\rangle$ denotes the fluctuation of
Ly$\alpha$ surface brightness for the $i$-th pixel in this bin. Here,
$\mu_{ri}$ is the residual surface brightness calculated by subtracting the
bestfit galaxy model spectra from the observed LRG spectra and dividing the
residuals by the angular area of SDSS fiber, and $\langle\mu(z)\rangle$ is the
average residual surface brightness at each redshift (Figure 1), obtained by
stacking the surface brightness of all residual LRG spectra in the observed
frame. The spectral interval $\Delta\log_{10}\lambda=10^{-4}$ (about 69 ${\rm
km\,s^{-1}}$ per pixel) in the SDSS spectra is kept when we compute
$\langle\mu(z)\rangle$. The pixel weight $\omega_{ri}$ is the inverse variance
of the flux, $1/\sigma_{ri}^{2}$, for valid pixels and zero for masked pixels.
To avoid the stray light contamination from quasars on the CCD, similar to
Croft et al. (2016), we exclude any LRG spectrum once it is observed within 5
fibers or less away from a quasar fiber (i.e., $\Delta_{\rm fiber}\leq 5$), as
discussed in Appendix A.1. A more detailed analysis of the potential
contamination in our measurement and the correction to possible systematics
are discussed in Appendix A.

Figure 1: The average residual surface brightness $\langle\mu(z)\rangle$,
obtained by averaging all the individual residual spectra in the observed
frame after subtraction of the bestfit galaxy model sepctra from the LRG
spectra. The gray regions, centered at 3934Å and 3969Å spanning 30Å,
respectively, and 4358Å spanning 40Å, are masked for zero-redshift Ca H&K
lines and the strong Mercury G line from streetlamps.

Note that the average residual surface brightness shown in Figure 1 differs
from that in Croft et al. (2016), mainly due to improved algorithms in flux-
calibration and extraction for
DR16222https://www.sdss.org/dr16/spectro/pipeline/#ChangesforDR16. Regardless,
the strong features at the zero redshift calcium H and K lines and Mercury G
line remain the largest excursions. This difference has little impact on our
following analysis, since the residual continuum contributes only to the
statistical noise in the measurement, not the signal. The feature around 4050Å
might be due to one of the sky emission lines, Hg I 4047Å. We decide not to
mask it, as we also do not specially deal with regions where other sky lines
may reside. Since it is not the flux of $\langle\mu(z)\rangle$ itself but the
fluctuation level $\Delta_{\mu}$ relative to it matters (see Equation 4), this
feature, shared by all individual residual spectra, will not affect our cross-
correlation measurements.

We show the quasar-Ly$\alpha$ emission cross-correlation on a linear scale in
Figure 2. The contours are somewhat stretched along the $r_{\|}$ direction for
$r_{\perp}$ below a few $h^{-1}$Mpc. Croft et al. (2016) quantified the
redshift-space anisotropies by assuming a linear $\Lambda$CDM correlation
function shape distorted by a peculiar velocity model, which includes standard
linear infall for large-scale flows and a small-scale random velocity
dispersion. In fact the elongation in the $r_{\|}$ direction can be caused by
a combination of multiple factors, including the intrinsic velocity dispersion
of quasars in their host halos, the intrinsic velocity dispersion of the
sources of Ly$\alpha$ emission, and quasar redshift uncertainties. The
uncertainty in quasar redshifts primarily comes from systematics offsets
between measured redshifts adopting different indicators, which can sometimes
become large due to the complexity of physical processes related to broad
emission lines. That makes it difficult to precisely and accurately distangle
a systemic redshift. For example, the variation of quasar redshift offsets
between Z_PCA (redshift estimated by PCA) and Z_MgII (redshift indicated by
MgII emission lines) in DR16 can reach over $\pm$ 500 ${\mathrm{km~{}s^{-1}}}$
(Pâris et al., 2018; Lyke et al., 2020b; Brodzeller & Dawson, 2022), which
corresponds to $\sim\pm 4.7h^{-1}$cMpc.

Figure 2: The quasar-Ly$\alpha$ emission cross-correlation as a function of
$r_{\perp}$ and $r_{\|}$. To reduce noise in the image, the data is smoothed
with a 2D Gaussian kernel with a standard deviation of 4 $h^{-1}$cMpc.
Potential light contamination is removed by pixel veto. For display, the
pattern is mirrored along $r_{\perp}=0$.

###  3.2 Projected Ly$\alpha$ Emission Surface Brightness

In this subsection, we measure the projected Ly$\alpha$ SB profile in a psedo-
narrow band by collapsing the 2D cross-correlation along the line-of-sight
direction. There are previous studies of Ly$\alpha$ surface brightness profile
around quasars. In order to compare to our derived profile, we first summarize
those observations.

Cai et al. (2019) studied quasar circumgalactic Ly$\alpha$ emission using KCWI
observations of 16 ultraluminous Type I QSOs at $z=2.1-2.3$. They integrated
over a fixed velocity range of $\pm$1000 km s-1 around the centroid of
Ly$\alpha$ nebular emission to calculate the SB. The median Ly$\alpha$ SB
profile in their work can be described by the following power-law profile
centered at the QSO at projected radius $r_{\perp}$ of 15-70 pkpc, which we
denote as ${\rm SB_{C}}$:

$\displaystyle\mathrm{SB}_{\mathrm{C}}(z\approx 2.3)$ $\displaystyle=3.7\times
10^{-17}\times(r_{\perp}/10\ \mathrm{pkpc})^{-1.8}$ (5)
$\displaystyle\mathrm{erg}\
\mathrm{s}^{-1}\mathrm{~{}cm}^{-2}\operatorname{arcsec}^{-2}.$

Borisova et al. (2016) found large Ly$\alpha$ nebulae on the spatial extent of
$>$100 pkpc from a MUSE snapshot survey on 17 radio quiet QSOs at $z>3.1$.
Twelve of them are selected specifically for their study from the catalog of
Véron-Cetty & Véron (2010), as the brightest radio-quiet quasars known in the
redshift range of $z=3.0-3.3$, and the other five at $z=3.6-4.0$ are selected
originally for studying absorption line systems in quasar spectra. They fixed
the width of their pseudo-NB images to the maximum spectral width of the
Ly$\alpha$ nebulae, with a median of 43Å. The median of their integrated SB
profiles, denoted as ${\rm SB}_{\rm B}$ here, can be described as:

$\displaystyle\mathrm{SB}_{\mathrm{B}}(z\approx 3.1)$ $\displaystyle=3.2\times
10^{-17}\times(r_{\perp}/10\ \mathrm{pkpc})^{-1.8}$ (6)
$\displaystyle\operatorname{erg}\mathrm{s}^{-1}\mathrm{~{}cm}^{-2}\operatorname{arcsec}^{-2}.$

Besides, Croft et al. (2018) used a power-law,

$\displaystyle{\rm SB_{Croft}}(z\approx 2.55)$ $\displaystyle=3.5\times
10^{-19}\times\left(r_{\perp}/{\rm cMpc}\right)^{-1.5}$ (7) $\displaystyle{\rm
erg}\ {\rm s}^{-1}{\rm cm}^{-2}{\rm arcsec}^{-2},$

to follow the broad trend seen in the data.

If we make a simple correction for cosmological surface brightness dimming to
$z=2.40$, the median redshift of our quasar sample, by scaling with a factor
of $(1+z)^{4}$, the above SB profiles become

$\displaystyle\mathrm{SB}_{\mathrm{C}}(z\approx 2.40)=$ $\displaystyle
3.3\times 10^{-17}\times(r_{\perp}/10\ \mathrm{pkpc})^{-1.8}$ (8)
$\displaystyle\mathrm{erg}\
\mathrm{s}^{-1}\mathrm{~{}cm}^{-2}\operatorname{arcsec}^{-2},$
$\displaystyle\mathrm{SB}_{\mathrm{B}}(z\approx 2.40)=$ $\displaystyle
6.8\times 10^{-17}\times(r_{\perp}/10\ \mathrm{pkpc})^{-1.8}$
$\displaystyle\operatorname{erg}\mathrm{s}^{-1}\mathrm{~{}cm}^{-2}\operatorname{arcsec}^{-2},$

and

$\displaystyle{\rm SB_{Croft}}(z\approx 2.40)$ $\displaystyle=4.16\times
10^{-19}\times\left(r_{\perp}/{\rm cMpc}\right)^{-1.5}$ (9) $\displaystyle{\rm
erg}\ {\rm s}^{-1}{\rm cm}^{-2}{\rm arcsec}^{-2}.$

To properly compare our measured SB with these previous work, we firstly
collapse the 2D cross-correlation measurement in Section 3.1 along $r_{\|}$ to
obtain the SB as a function of $r_{\perp}$. We integrate the cross-correlation
over a fixed line-of-sight window of $\pm 1000$ km s-1, corresponding to a
window spanning $\pm 4$Å around $\lambda_{{\rm Ly}\alpha}\approx 1216$Å in the
$z=2.40$ quasar rest frame, or a window of $\pm 9.37h^{-1}{\rm cMpc}$ around
the quasar.

We use the jackknife method to compute the standard deviation of the obtained
SB, by drawing a jackknife sample set from the 885 LRG subsamples and perform
a cross-correlation with the quasar sample. The covariance matrix $C_{ij}$ can
be written as:

$\displaystyle C_{ij}(r_{\perp,i},r_{\perp,j})=\frac{n-1}{n}\times$ (10)
$\displaystyle\sum_{k=1}^{n}\left[{\rm
SB}_{k}\left(r_{\perp,i}\right)-\overline{\rm
SB}\left(r_{\perp,i}\right)\right]\left[{\rm
SB}_{k}\left(r_{\perp,j}\right)-\overline{\rm
SB}\left(r_{\perp,j}\right)\right],$

where ${\rm SB}_{k}\left(r_{\perp,i}\right)$ is the surface brightness in bin
$i$ centered at the transverse separation $r_{\perp,i}$ for the jackknife
sample $k$, $\overline{\rm SB}(r_{\perp,i})$ denotes the surface brightness
measured from the full LRG data set, and the number of jackknife samples, $n$,
is 885.

As shown in Figure 3, we have a detection of the SB profile at projected
radius $r_{\perp}$ ranging from $\sim$0.1$h^{-1}$cMpc to
$\sim$100$h^{-1}$cMpc.

The SB profile within $r_{\perp}\leq 0.5\ h^{-1}{\rm cMpc}$ appears to be
consistent with the observations of the QSO nebulae on smaller scales in Cai
et al. (2019) and Borisova et al. (2016), and on scales of $1h^{-1}{\rm
cMpc}\leq r_{\perp}\leq 10h^{-1}{\rm cMpc}$ our profile broadly agrees with
the power-law fit in Croft et al. (2018).

Figure 3:  Projected Ly$\alpha$ surface brightness profile (red points) around
quasars obtained from our cross-correlation measurement. For comparison, the
power-law fit (Equation 9) from the intensity mapping result in Croft et al.
(2018) is shown as the black dashed line. The SB profiles from observations of
Ly$\alpha$ emission around quasars on smaller scales are shown as the green
shaded region (representing the range of 25th and 75th percentiles in Cai et
al. 2019) and purple points (Borisova et al., 2016), with the green and purple
dashed lines denoting the power-law fit and extrapolation (Equation 8). In the
bottom panel, the measured SB is shown in linear scale.

###  3.3 Multipoles of the Redshift-Space Two-Point Correlation Function

In addition to measure the quasar-Ly$\alpha$ emission cross-correlation
function (a.k.a. two-point correlation function; 2PCF) in bins of $r_{\perp}$
and $r_{\|}$, to better describe its shape, we further measure the cross-
correlation in bins of $s$ and $\mu$, where $s$ is the separation between
quasars and Ly$\alpha$ pixels, i.e., $s=\sqrt{r_{\perp}^{2}+r_{\|}^{2}}$, and
$\mu$ is cosine of the angle between $\vec{s}$ and the line-of-sight
direction, $\mu=r_{\|}/s$.

The redshift-space 2PFC $\xi(s,\mu)$ can be expanded into multipoles, with the
multipole moment $\xi_{\ell}$ calculated by (Hamilton, 1992):

$\xi_{l}(s)=\frac{2l+1}{2}\int_{-1}^{1}\xi(s,\mu)\mathcal{L}_{\ell}(\mu)\mathrm{d}\mu,\\\
$ (11)

where $\mathcal{L}_{\ell}$ is the $\ell$-th order Legendre polynomial. In the
linear regime (Kaiser, 1987), there are three non-zero components of the
redshift-space 2PCF: the monopole $\xi_{0}$, the quadrupole $\xi_{2}$ and the
hexadecapole $\xi_{4}$,

$\xi(s,\mu)=\sum_{\ell=0,2,4}\xi_{\ell}(s)\mathcal{L}_{\ell}(\mu).$ (12)

At small transverse separations, however, the redshift-space 2PCF is affected
by the small-scale non-linear effect, such as the Finger-of-God (FoG) effect,
and also the quasar redshift uncertainty in our cases. To reduce the small-
scale contamination, we follow McCarthy et al. (2019) to adopt the truncated
forms of the multipoles by limiting the calculation to large transverse
separations ($r_{\perp}>r_{\perp,{\rm cut}}$),

$\hat{\xi}_{\ell}=\frac{2\ell+1}{2}\int_{-\mu_{\max}}^{\mu_{\max}}\xi(s,\mu)\mathcal{L}_{\ell}(\mu)d\mu,$
(13)

where $\mu_{\max}=\sqrt{1-(r_{\perp,\mathrm{cut}}/s)^{2}}$. The transformation
between $\mathbf{\xi}=\left(\xi_{0},\xi_{2},\xi_{4}\right)^{T}$ and
$\mathbf{\hat{\xi}}=\left(\hat{\xi}_{0},\hat{\xi}_{2},\hat{\xi}_{4}\right)^{T}$
can be described using a $3\times 3$ matrix $\mathbf{R}$:

$\mathbf{\hat{\xi}}=\mathbf{R}\mathbf{\xi},$ (14)

where

$R_{\ell
k}=\frac{2\ell+1}{2}\int_{-\mu_{\max}}^{\mu_{\max}}\mathcal{L}_{\ell}(\mu)\mathcal{L}_{k}(\mu)d\mu\quad{\rm
for~{}}\ell,k=0,2,4.$ (15)

In our measurement we set $r_{\perp,{\rm cut}}=4h^{-1}{\rm cMpc}$ to ensure
that bulk of small-scale contamination is excluded. The multipole measurements
will be presented along with the modeling results.

###  3.4 Modeling the Quasar-Ly$\alpha$ Emission Cross-correlation

In Croft et al. (2016), the amplitude of the measured quasar-Ly$\alpha$
emission cross-correlation, if modelled by relating Ly$\alpha$ emission to
star-forming galaxies, would imply a value of Ly$\alpha$ emissivity to be
comparable to that inferred from the cosmic SFRD without dust correction,
appearing too high compared with the predictions from the Ly$\alpha$
luminosity functions (LF) of Ly$\alpha$ emitting galaxies. In Croft et al.
(2018), with the correction to the systematic effect from quasar clustering
and the complementary measurement of Ly$\alpha$ forest-Ly$\alpha$ emission
cross-correlation, the detected Ly$\alpha$ emission is found to be explained
by Ly$\alpha$ emission associated with quasars based on populating a large
hydrodynamic cosmological simulation. In this subsection we will revisit both
scenarios by constructing a simple analytic model to describe the measured
Ly$\alpha$ intensity, and argue that the observed Ly$\alpha$ emission cannot
be only contributed by quasars. The simple model can also be applied to
Ly$\alpha$ forest-Ly$\alpha$ emission cross-correlation, and our corresponding
prediction and detailed analysis are presented in Section 4.

We assume that the Ly$\alpha$ emission from sources clustered with quasars
contribute the bulk of the detected signals on large scales, while on small
scales the Ly$\alpha$ photons associate with the central quasar count.
Supposing that $\langle\mu_{\alpha}\rangle$ is the mean surface brightness of
Ly$\alpha$ emission, $b_{q}$ and $b_{\alpha}$ are the linear bias factors of
quasars and Ly$\alpha$ sources, respectively, in the linear regime the non-
vanishing multipoles of the redshift-space quasar-Ly$\alpha$ emission cross-
correlation are given by

$\displaystyle{\xi}_{0}(s)=b_{q}b_{\alpha}\langle\mu_{\alpha}\rangle
f_{\beta,0}{\xi}_{mm}(r),$ (16)
$\displaystyle{\xi}_{2}(s)=b_{q}b_{\alpha}\langle\mu_{\alpha}\rangle
f_{\beta,2}\left[{\xi}_{mm}(r)-\bar{\xi}_{mm}(r)\right],$
$\displaystyle{\xi}_{4}(s)=b_{q}b_{\alpha}\langle\mu_{\alpha}\rangle
f_{\beta,4}\left[\xi_{mm}(r)+\frac{5}{2}\bar{\xi}_{mm}(r)-\frac{7}{2}\bar{\bar{\xi}}_{mm}(r)\right],$

where (e.g., Percival & White, 2009)

$\displaystyle f_{\beta,0}$
$\displaystyle=1+\frac{1}{3}\left(\beta_{q}+{\beta}_{\alpha}\right)+\frac{1}{5}\beta_{q}{\beta}_{\alpha},$
(17) $\displaystyle f_{\beta,2}$
$\displaystyle=\frac{2}{3}\left(\beta_{q}+{\beta}_{\alpha}\right)+\frac{4}{7}\beta_{q}{\beta}_{\alpha},$
$\displaystyle f_{\beta,4}$
$\displaystyle=\frac{8}{35}\beta_{q}{\beta}_{\alpha},$

and (e.g., Hawkins et al., 2003)

$\displaystyle\bar{\xi}(r)=\frac{3}{r^{3}}\int_{0}^{r}\xi\left(r^{\prime}\right)r^{\prime
2}\mathrm{~{}d}r^{\prime},$ (18)
$\displaystyle\bar{\bar{\xi}}(r)=\frac{5}{r^{5}}\int_{0}^{r}\xi\left(r^{\prime}\right)r^{\prime
4}\mathrm{~{}d}r^{\prime}.$

Note that $r$ is for the distance in the real space and $s$ denotes the
distance in the redshift space and in the above expressions $r=s$. Then the
model for the truncated two-point correlation function $\hat{\xi}$ can be
obtained according to Equations (14) and (15).

The redshift-space distortion parameter $\beta_{q}$ for quasars depicts the
redshift-space anisotropy caused by peculiar velocity,
$\beta_{q}=\Omega_{m}^{0.55}(z=2.4)/b_{q}$. We fix $b_{q}=3.64$ according to
Font-Ribera et al. (2013). The redshift-space distortion parameter
$\beta_{\alpha}$ for Ly$\alpha$ emission is similarly defined. We set
$b_{\alpha}=b_{q}$ for the case that the main contributors to Ly$\alpha$
emission are clustered quasars and $b_{\alpha}=3$ for the case that Ly$\alpha$
emission is dominated by contributions from star-forming galaxies. A value of
3 appears to be a good estimate of the luminosity-weighted bias $b_{\alpha}$
for star-forming galaxies. Following Croft et al. (2016), we find that
$b_{\alpha}$ is within $\sim$5% of 3 with different low halo mass cuts and
different prescriptions of the stellar mass-halo mass relation at $z\sim 2.4$
(e.g., Moster et al., 2010, 2013; Behroozi et al., 2019). In both scenarios we
leave $\beta_{\alpha}$ and $\langle\mu_{\alpha}\rangle$ as free parameters to
be fitted. We note that $\beta_{\alpha}$ can potentially include additional
effects other than the Kaiser effect, such as the Ly$\alpha$ radiative
transfer on clustering (Zheng et al., 2011a).

We also model the Ly$\alpha$ SB profile. As discussed in Section 3.2, previous
observations indicate that the small-scale SB profile can be well described by
a power law with an index of $-1.8$. We therefore decompose the full SB
profile into two components: the one-halo term ${\rm SB_{1h}}$ dominated by
Ly$\alpha$ emission associated with the central quasars and the two halo-term
${\rm SB_{2h}}$ by the clustered Ly$\alpha$ sources,

$\displaystyle{\rm SB_{1h}}$ $\displaystyle={\rm
SB_{0}}\left(\frac{r_{\perp}}{1h^{-1}{\rm cMpc}}\right)^{-1.8},$ (19)
$\displaystyle\mathrm{SB}_{2\mathrm{h}}$ $\displaystyle=\frac{\rho_{\rm
Ly\alpha}}{4\pi(1+z)^{2}}\int_{\pi_{\rm min}}^{\pi_{\rm
max}}\xi(r_{\perp},r_{\|})dr_{\|}$ $\displaystyle=\frac{\rho_{\rm
Ly\alpha}}{4\pi(1+z)^{2}}b_{q}b_{\alpha}\left(f_{\beta,0}w_{p,0}+f_{\beta,2}w_{p,2}+f_{\beta,4}w_{p,4}\right).$

Here $\xi$ is the linear correlation function between quasars and Ly$\alpha$
emission sources (quasars or star-forming galaxies) in redshift space,
$\rho_{\rm
Ly\alpha}=4\pi\langle\mu_{\alpha}\rangle[H(z)/c]\lambda_{\alpha}(1+z)^{2}$ is
the comoving Ly$\alpha$ luminosity density (Croft et al., 2016), $\pi_{\max}$
and $\pi_{\min}$ correspond to $\pm$9.37$h^{-1}$cMpc, the width of the pseudo-
narrow band used in § 3.2. The projected cross-correlation function is put in
the form of the projected multipoles, which are calculated as

$\displaystyle w_{p,0}(r_{\perp})$
$\displaystyle=\int_{\pi_{\min}}^{\pi_{\max}}\xi_{mm}(r)\mathcal{L}_{0}(\mu)dr_{\|},$
(20) $\displaystyle w_{p,2}(r_{\perp})$
$\displaystyle=\int_{\pi_{\min}}^{\pi_{\max}}\left[{\xi}_{mm}(r)-\bar{\xi}_{mm}(r)\right]\mathcal{L}_{2}(\mu)dr_{\|},$
$\displaystyle w_{p,4}(r_{\perp})$
$\displaystyle=\int_{\pi_{\min}}^{\pi_{\max}}\left[\xi_{mm}(r)+\frac{5}{2}\bar{\xi}_{mm}(r)-\frac{7}{2}\bar{\bar{\xi}}_{mm}(r)\right]\mathcal{L}_{4}(\mu)dr_{\|}.$

with $r=\sqrt{r_{\perp}^{2}+r_{\|}^{2}}$ and $\mu=r_{\|}/r$.

With three free parameters (${\rm SB_{0}}$, $\beta_{\alpha}$, and
$\langle\mu_{\alpha}\rangle$), we perform a joint fit to the three (large-
scale) multipoles and the projected SB profile, assuming that the Ly$\alpha$
sources in the model are mainly quasars and star-forming galaxies,
respectively, as discussed in Section 3.4.1 and Section 3.4.2.

####  3.4.1 Star-forming Galaxies as Ly$\alpha$ Sources

In the case that Ly$\alpha$ emission is dominated by the contribution from
galaxies, we fix $b_{\alpha}=3$. The best-fit results for the multipoles and
the SB profile are shown in Figure 4 and Figure 5. Given the uncertainties in
the measurements, the model provides a reasonable fit and shows a broad
agreement with the trend in the data. The middle panel of Figure 6 shows a
reconstructed 2D image of the redshift-space linear cross-correlation function
from the best-fit model. If it is subtracted from the measurement (left
panel), the residual (right panel) is dominated by the small-scale clustering
that we do not model.

The constraints on the three parameters are presented in Figure 7. The
parameter representing the amplitude of the one-halo term is loosely
contrained, ${\rm SB_{0}}=3.49_{-2.02}^{+2.27}\times 10^{-20}{\rm
erg\,s^{-1}cm^{-2}arcsec^{-2}}$. The parameter $\langle\mu_{\alpha}\rangle$,
proportional to the comoving Ly$\alpha$ emissivity or luminosity density, is
constrained at the 2$\sigma$ level,
$\langle\mu_{\alpha}\rangle=1.13^{+0.57}_{-0.53}\times 10^{-21}{\rm erg\
s^{-1}cm^{-2}\textup{\AA}^{-1}arcsec^{-2}}$. The redshift-space distortion
parameter has a high probability density of being negative but with a tail
toward positive values, $\beta_{\alpha}=0.07^{+1.65}_{-0.73}$. Given its
uncertainty, the value is consistent with that from the Kaiser effect,
$\Omega_{m}(z=2.4)^{0.55}/b_{\alpha}\simeq 0.32$, and we are not able to tell
whether there is any other effect (e.g., caused by radiative transfer; Zheng
et al. 2011a).

We note that fitting the clustering measurements leads to an anti-correlation
between $\langle\mu_{\alpha}\rangle$ and $\beta_{\alpha}$ (Eq. 16 and Eq. 17;
Fig. 7). If $\beta_{\alpha}$ is restricted to the formal value of $\sim$0.32
from the Kaiser effect, the constraints on $\langle\mu_{\alpha}\rangle$ become
$1.09_{-0.24}^{+0.25}\times 10^{-21}{\rm erg\
s^{-1}cm^{-2}\textup{\AA}^{-1}arcsec^{-2}}$, a nearly 4$\sigma$ detection. If
we set the upper limit of $\beta_{\alpha}$ to be 0.32 to allow room for
radiative transfer effect (e.g., Zheng et al., 2011a), the constraints change
to $\langle\mu_{\alpha}\rangle=1.44_{-0.38}^{+0.45}\times 10^{-21}{\rm erg\
s^{-1}cm^{-2}\textup{\AA}^{-1}arcsec^{-2}}$. In the following discussions, to
be conservative, we take the $\langle\mu_{\alpha}\rangle$ constraints without
these restrictions.

The constrained $\langle\mu_{\alpha}\rangle$ corresponds to a comoving
Ly$\alpha$ luminosity density of $\rho_{\rm Ly\alpha}=6.6_{-3.1}^{+3.3}\times
10^{40}{\rm erg\,s^{-1}cMpc^{-3}}$. This value is about 3.6 times lower than
that in Croft et al. (2016) or $\sim 2.2$ times lower than that in Croft et
al. (2018). With the lower amplitude, the fractional uncertainty is larger.
The comparison is shown in Figure 8. We also show the Ly$\alpha$ luminosity
densities at different redshifts calculated by integrating the Ly$\alpha$ LFs
of LAEs down to low luminosity. For example, LFs in Ouchi et al. (2008, 2010)
are integrated down to $L_{\rm Ly\alpha}=0$ with the best-fit Schechter
parameters for $z=3.1,3.7,5,7$ and $6.6$; that in Drake et al. (2017a) down to
$\log[L_{\rm Ly\alpha}/({\rm erg\,s^{-1}})]=41.0$; that in Sobral et al.
(2018) down to $1.75\times 10^{41}{\rm erg\,s^{-1}}$. These quoted Ly$\alpha$
luminosity densities are inferred without separating the contribution of
potential AGNs except at the very luminous end (see Wold et al. 2017 for a
two-component fit). The luminous end is usually excluded in the parametrized
fits to the Ly$\alpha$ LFs, but they do not contribute much to the total
Ly$\alpha$ luminosity density due to their rather low number density. The
quoted LAE Ly$\alpha$ luminosity densities in Figure 8 should have included
the potential contribution of relatively faint AGNs (with AGNs detected in
X-ray and radio contributing at a level of a few percent; Sobral et al. 2018).
At $z\sim 2.4$, our inferred Ly$\alpha$ luminosity density is about one order
of magnitude higher than that inferred from the LAE LF, although they can be
consistent within the uncertainty.

We further show the H$\alpha$-converted Ly$\alpha$ luminosity density as done
in Wold et al. (2017), which is obtained by scaling the H$\alpha$ luminosity
density measured in the HiZELS survey (Sobral et al., 2013) with an escape
fraction of 5% and a correction about 10%(15%) for AGN contribution at $z<1$
($z>1$). The cosmic Ly$\alpha$ luminosity density measured by Chiang et al.
(2019) through broad-band intensity mapping is also shown, which probes the
total background including low surface brightness emission by spatially cross-
correlating photons in far-UV and near-UV bands with spectroscopic objects.
They claimed that their derived cosmic Ly$\alpha$ luminosity density is
consistent with cosmic star formation with an effective escape fraction of 10%
assuming that all of the Ly$\alpha$ photons originate from star formation.
Combined our measurement with the results of Chiang et al. (2019), it appears
that the cosmic Ly$\alpha$ luminosity density grows with redshift over
$0\lesssim z\lesssim 2.5$, and more data points at different redshifts are
expected to confirm this trend.

If we assume that all the Ly$\alpha$ emission originates from star formation,
we can convert our inferred Ly$\alpha$ luminosity density to a SFRD, by using
a simple conversion (Kennicutt, 1998),

$\rho_{\rm SFR}/(M_{\odot}\mathrm{yr}^{-1}\mathrm{cMpc}^{-3})=\frac{\rho_{\rm
Ly\alpha}/({\rm erg~{}s^{-1}~{}cMpc^{-3}})}{1.1\times 10^{42}({\rm
erg\,s^{-1}})/(M_{\odot}{\rm yr}^{-1})}.$ (21)

This gives $\rho_{\rm SFR}=0.06\pm
0.03M_{\odot}\mathrm{yr}^{-1}\mathrm{cMpc}^{-3}$, higher than that from
integrating LAE LFs, as shown in Figure 9. The value is on the low end of the
cosmic star formation rate density based on UV and infrared observations
(e.g., Robertson et al., 2015). However, we emphasize that the
Ly$\alpha$-converted $\rho_{\rm SFR}$ in this case should be treated as a
lower limit for estimates of the intrinsic star formation, since no correction
is applied to account for dust extinction and Ly$\alpha$ escape fraction. The
comparison in Figure 9 is simply to highlight the high amplitude of Ly$\alpha$
emission inferred from the quasar-Ly$\alpha$ emission cross-correlation.

Figure 4: Modified monopole, quadrupole and hexadecapole of the quasar-
Ly$\alpha$ emission cross-correlation (see Equation 13) and their fitting
results based on the galaxy-dominated model (see Section 3.4.1). The points
represent our measurements with jackknife error bars. The solid curves denote
modelled modified multipoles with parameters randomly drawn from their
posterior probability distributions, among which the thickest ones correspond
to the best-fits. The modified multipoles remove any information within
$r<r_{\perp,{\rm cut}}$, i.e., the gray-shaded regions, to avoid small-scale
contamination.  Figure 5: Ly$\alpha$ SB profile. The data points are from
integrating the measured quasar-Ly$\alpha$ emission cross-correlation function
along the line of sight, and the solid curve is the best-fit SB profile for
the galaxy-dominated model depicted in Section 3.4.1. The dashed lines denote
the best-fit one-halo and two-halo term, respectively, and the shaded region
represents the $\pm 1\sigma$ range. In the bottom panel, linear scale is used
in the $y$-axis.  Figure 6: The measured (left panel), best-fit (middle panel)
and residual (right panel) quasar-Ly$\alpha$ emission cross-correlation as a
function of $r_{\|}$ and $r_{\perp}$. The model fit is only to large-scale
signals by using the modified multipoles (Equation 13). The best-fit pattern
shown here is reconstructed from the corresponding multipoles (Equation 16)
with the best-fit parameters. The residual is obtained by subtracting the
best-fit model from the measurement, with elongated distortion along the
$r_{\|}$ direction on small scales, the small-scale anisotropy not included in
our model. All the three images are smoothed using a 2D Gaussian kernel with a
standard deviation of 4$h^{-1}$cMpc.  Figure 7: The probability distribution
of parameters $\langle\mu_{\alpha}\rangle$, $\beta_{\alpha}$ and ${\rm
SB_{0}}$ as a result of the joint fit to the modified multipoles and SB
profile of quasar–Ly$\alpha$ emission cross-correlation, with an assumption
that star-forming galaxies dominate the large-scale Ly$\alpha$ emission and
thus $b_{\alpha}=3$. $\langle\mu_{\alpha}\rangle_{-21}$ is
$\langle\mu_{\alpha}\rangle$ in units of $10^{-21}{\rm
erg\,s^{-1}cm^{-2}\textup{\AA}^{-1}arcsec^{-2}}$, and ${\rm SB}_{-20}$ is
${\rm SB_{0}}$ in units of $10^{-20}{\rm erg\,s^{-1}cm^{-2}arcsec^{-2}}$. The
dashed lines in the histograms denote 16th, 50th and 84th percentiles of the
marginalized distributions. Figure 8: Ly$\alpha$ luminosity density $\rho_{\rm
Ly\alpha}$. The red star shows the value inferred from our quasar-Ly$\alpha$
emission measurement, assuming that the detected Ly$\alpha$ emission is due to
star-forming galaxies with a typical luminosity-weighted bias of
$b_{\alpha}=3$. As a comparison, we also show the values with previous
intensity mapping measurements (Croft et al., 2016, 2018; Chiang et al., 2019)
and those from integrating the Ly$\alpha$ LFs of LAEs (Ouchi et al., 2008,
2010; Drake et al., 2017a; Sobral et al., 2018; Hu et al., 2019; de La
Vieuville et al., 2019) and scaling H$\alpha$ luminosities with an escape
fraction of 5% (Sobral et al., 2013; Wold et al., 2017).  Figure 9:  Same as
Fig. 8 but the Ly$\alpha$ luminosity density $\rho_{\rm Ly\alpha}$ is
converted to star formation rate density $\rho_{\rm SFR}$ under the assumption
that Ly$\alpha$ emission is purely caused by star formation. Given the effect
of dust extinction and Ly$\alpha$ escape fraction, these Ly$\alpha$-converted
$\rho_{\rm SFR}$ values should be considered as lower limits of the intrinsic
star formation. The orange shaded region represents the parameterized model
for the evolving star formation rate density in Robertson et al. (2015), based
on infrared and ultraviolet observations.

####  3.4.2 Quasars as Ly$\alpha$ Sources

In the case that Ly$\alpha$ emission is dominated by the contribution from
quasars, we make a simple assumption that the quasars involved are almost the
same, with a typical Ly$\alpha$ luminosity $L_{q,\alpha}$ and a comoving
number density $n_{q}$, so that $\rho_{\rm Ly\alpha}=L_{\alpha,q}n_{q}$.

We calculate $n_{q}$ by integrating the luminosity evolution and density
evolution (LEDE) model (Ross et al., 2013) of the optical quasar luminosity
function (QLF), fitted using data from SDSS-III DR9 and allowing luminosity
and density to evolve independently. The QLF gives the number density of
quasars per unit magnitude, and its integration over the magnitude range from
$M_{i}[z=2]=-30$ to $M_{i}[z=2]=-18$ yields $n_{q}\approx 1.34\times
10^{-4}h^{3}{\rm Mpc}^{-3}$

With the analytical model in this quasar-dominant scenario, we jointly fit
both the measured cross-correlation multipoles and the SB profile, where
$b_{\alpha}$ is fixed to be $b_{q}$ and $\rho_{\rm Ly\alpha}$ is interpreted
to be $L_{\alpha,q}n_{q}$, leaving $L_{q,\alpha}$, $\beta_{\alpha}$ and ${\rm
SB_{0}}$ as free parameters. Our joint fitting result, presented in Figure 10,
indicates that the required quasar Ly$\alpha$ luminosity under the above
assumption should be $\log[L_{q,\alpha}/({\rm
erg\,s^{-1}})]=45.12^{+0.18}_{-0.27}$. The best-fit value is even brighter
than some ultraluminous quasars usually targeted to search for enormous
nebulae (e.g., $\sim 10^{43}$ – $\lesssim 10^{45}{\rm erg\,s^{-1}}$ in Cai et
al. 2018). Such a high Ly$\alpha$ luminosity per quasar makes the quasar-
dominated model unlikely to work.

Figure 10: The probability distribution of $L_{q,\alpha}$ and $\beta_{\alpha}$
from the joint fit to the multipoles and SB profiles from the measured quasar-
Ly$\alpha$ emission cross-correlation. The mean quasar Ly$\alpha$ luminosity
$L_{q,\alpha}$ is in units of ${\rm erg\,s^{-1}}$. Note that there are
actually three parameters, $L_{q,\alpha}$, $\beta_{\alpha}$ and ${\rm
SB_{0}}$, in the model, but here we focus on the constraints on $L_{q,\alpha}$
and $\beta_{\alpha}$. See the text for detail.

Our modeling result appears to be inconsistent with the quasar-dominated model
in Croft et al. (2018). In their model, the Ly$\alpha$ SB profile on scales
above $\sim 1h^{-1}{\rm Mpc}$ is well reproduced (see their Fig.10).
Ly$\alpha$ emission in their model is presented as Ly$\alpha$ SB as a function
of gas density and distance from the quasar, while the total Ly$\alpha$
luminosity per quasar is not given. The luminosity, however, can be dominated
on scales $\lesssim 1h^{-1}{\rm Mpc}$, which is not shown in their figure.
Fortunately, panel (b) in their Figure 8 (“Model Q”) enables an estimation of
the mean quasar Ly$\alpha$ luminosity (R. Croft, private communication). With
a mean Ly$\alpha$ SB $\langle\mu_{\alpha}\rangle=7.0\times 10^{-22}{\rm
erg\,s^{-1}cm^{-2}\textup{\AA}^{-1}arcsec^{-2}}$ (their section 5.1) from the
slice with thickness of $40h^{-1}{\rm Mpc}$ (corresponding to observed spread
of $\sim 29$Å in Ly$\alpha$ emission) and side length $400h^{-1}{\rm Mpc}$
($\sim 2.04\times 10^{4}{\rm arcsec}$ at $z\sim 2.5$), we obtain the total
Ly$\alpha$ luminosity in the slice to be $\sim 4.0\times 10^{47}{\rm
erg\,s^{-1}}$. As there are about 100 quasars in the slice, the average
Ly$\alpha$ luminosity in “Model Q” of Croft et al. (2018) is $\sim 4.0\times
10^{45}{\rm erg\,s^{-1}}$, which agrees well with our result here.

In conclusion, the modeling results from our analytical models rule out the
quasar-dominated scenario. For the galaxy-dominated scenario, however, both
our measurement and that in Croft et al. (2018) imply that the detected
Ly$\alpha$ signals cannot be explained simply by emission from currently
observed LAEs. There must be additional Ly$\alpha$ emitting sources other than
theses LAEs. We will explore the possibilities in Section 5 after presenting
the Ly$\alpha$ forest-Ly$\alpha$ emission cross-correlation results in Section
4.

##  4 Ly$\alpha$ forest-Ly$\alpha$ emission cross-correlation

Ly$\alpha$ forest, as a probe of the cosmic density field, can be used as an
alternative tracer, more space-filling than quasars, to detect diffuse
Ly$\alpha$ emission on cosmological scales. The Ly$\alpha$ forest-Ly$\alpha$
emission cross-correlation can provide additional information for
understanding the origin of the Ly$\alpha$ emission.

Following Croft et al. (2018), we measure the Ly$\alpha$ forest-Ly$\alpha$
emission cross-correlation in a way similar to quasar-Ly$\alpha$ emission
cross-correlation:

$\xi_{f\alpha}(r,\mu)=\frac{1}{\sum_{i=1}^{N(\vec{r})}w_{ri,\alpha}w_{ri,f}}\sum_{i=1}^{N(\vec{r})}w_{ri,\alpha}w_{ri,f}\Delta_{\mu,ri}\delta_{f,{ri}},$
(22)

where $N(\vec{r})$ is the number of Ly$\alpha$ forest-Ly$\alpha$ emission
pixel pairs within the bin centered at the separation $\vec{r}=(r,\mu)$.
$\Delta_{\mu,ri}$ is the fluctuation of Ly$\alpha$ emission SB (from the
residual LRG spectra) for the $i$-th pixel pair in this bin, and
$\delta_{f,ri}$ is the flux-transmission field of Ly$\alpha$ forest in the
quasar spectra. The weights $w_{ri,\alpha}$ of Ly$\alpha$ emission pixels are
the same as in Equation (4), and the weights for Ly$\alpha$ forest pixels
$w_{ri,f}=1/\sigma_{ri,f}^{2}$, where $\sigma_{ri,f}^{2}$ is the pixel
variance due to instrumental noise and large scale structure, with the latter
accounting for the intrinsic variance of the flux-transmission field.

Likewise, we can decompose the 2D Ly$\alpha$ forest-Ly$\alpha$ emission cross-
correlation into the monopole, quadrupole and hexadecapole moments. To avoid
spurious correlation induced by same-half-plate pixel pairs, we only use pixel
pairs residing on different half-plates and reject signals within $|r_{\|,{\rm
cut}}|=4$cMpc, as discussed in Appendix A.2. Similar to what we do with the
quasar-Ly$\alpha$ emission cross-correlation, we define the modified
multipoles of the Ly$\alpha$ forest-Ly$\alpha$ emission cross-correlation to
be

$\displaystyle\hat{\xi}_{f\alpha,\ell}(s)=\frac{2\ell+1}{2}$
$\displaystyle\left(\int_{-1}^{-\mu_{\min}}\xi_{f\alpha}(s,\mu)\mathcal{L}_{\ell}(\mu)d\mu\right.$
(23) $\displaystyle+$
$\displaystyle\left.\int_{\mu_{\min}}^{1}\xi_{f\alpha}(s,\mu)\mathcal{L}_{\ell}(\mu)d\mu\right),$

where $\mu_{\min}=|r_{\|,{\rm cut}}|/s$. Like Equation (15), the original and
modified multipoles are connected through
$\hat{\xi}_{f\alpha}=\mathbf{R^{\prime}}\xi_{f\alpha}$, where the element of
the transformation matrix $\mathbf{R^{\prime}}$ takes the form of

$\displaystyle R^{\prime}_{\ell k}=\frac{2\ell+1}{2}$
$\displaystyle\left(\int_{-1}^{-\mu_{\min}}\mathcal{L}_{\ell}(\mu)\mathcal{L}_{k}(\mu)d\mu\right.$
(24) $\displaystyle+$
$\displaystyle\left.\int^{1}_{\mu_{\min}}\mathcal{L}_{\ell}(\mu)\mathcal{L}_{k}(\mu)d\mu\right),$

with $\ell,k=0$, 2, and 4.

The analytical model for the Ly$\alpha$ forest-Ly$\alpha$ emission cross-
correlation is similar to the one for the quasar-Ly$\alpha$ emission cross-
correlation, and we only need to replace $b_{q}$ and $\beta_{q}$ in equations
(16) and (17) with $b_{f}$ and $\beta_{f}$, respectively. Here $b_{f}$ is the
Ly$\alpha$ forest transmission bias, evolving with redshift as
$b_{f}(z)=b_{f}(z_{\mathrm{ref}})[(1+z)/(1+z_{\mathrm{ref}})]^{\gamma_{\alpha}}$
with $\gamma_{\alpha}=2.9$, and $\beta_{f}$ is the redshift distortion
parameter for Ly$\alpha$ forest, $\beta_{f}=fb_{\eta}/b_{f}$, where $f$ is the
linear growth rate of structure and $b_{\eta}$ is the velocity bias of
Ly$\alpha$ forest (e.g., Seljak, 2012; Blomqvist et al., 2019). We fix
$b_{\eta}=-0.225$ and $\beta_{f}=1.95$ at a reference redshift of $z_{\rm
ref}=2.34$ according to the quasar-Ly$\alpha$ forest cross-correlation result
in du Mas des Bourboux et al. (2020), yielding $b_{f}=-0.119$ at $z=2.41$.

Figure 11: Modified monopole, quadrupole and hexadecapole of the Ly$\alpha$
forest-Ly$\alpha$ emission cross-correlation as a function of the Ly$\alpha$
forest-Ly$\alpha$ emission pixel pair separation. The data points are the
measurements, and the solid curves are the predictions using parameters in
literature to describe the Ly$\alpha$ forest and parameters
$\langle\mu_{\alpha}\rangle$ and $\beta_{\alpha}$ derived from fits to the
quasar-Ly$\alpha$ emission cross-correlation under the galaxy-dominated
scenario. The various solid curves are the predicted modified multipoles from
randomly drawing $\langle\mu_{\alpha}\rangle$ and $\beta_{\alpha}$ from their
posterior probability distribution, with the thickest ones from adopting the
best-fit parameters.

Given the small transmission bias $b_{f}$ of Ly$\alpha$ forest, the expected
Ly$\alpha$ forest-Ly$\alpha$ emission cross-correlation level at $\sim
10h^{-1}$cMpc is $\sim 5\%$ of the quasar-Ly$\alpha$ emission cross-
correlation. The subsequent low signal-to-noise ratio would lead to weak
parameter constraints from fitting the Ly$\alpha$ forest-Ly$\alpha$ emission
cross-correlation measurements. Instead we choose to compare the measurements
with the predictions from the model adopting the best-fit parameters,
$\beta_{\alpha}$ and $\langle\mu_{\alpha}\rangle$, from modeling the quasar-
Ly$\alpha$ emission correlation (Section 3.4.1). Such a consistency check is
shown in Figure 11.

The multipole measurements in Figure 11 indicate that there is no significant
detection of the Ly$\alpha$ forest-Ly$\alpha$ emission cross-correlation.
Quantitatively, a line of zero amplitude would lead to $\chi^{2}=19.8$ for a
total of 21 data points of the monopole, quadrupole, and hexadecapole in the
range of $4h^{-1}{\rm Mpc}<s<100h^{-1}{\rm Mpc}$. On the other hand, with the
large uncertainties in the data, our model predictions also appear to be
consistent with the measurements. The predictions from the bestfit model
(solid curves) give a value of $\chi^{2}=29.5$ for the above 21 data points,
within $\sim$1.3$\sigma$ of the expected mean $\chi^{2}$ value. We note that
the monopole is consistent with that in Croft et al. (2018), as long as the
uncertainties are taken into account (see their Fig.11). Our model has a much
lower amplitude than their galaxy-dominated model (model G), leading to a
closer match to the data. This is a manifestation of the lower
$\langle\mu_{\alpha}\rangle$ value inferred from our quasar-Ly$\alpha$
emission cross-correlation measurements.

##  5 Discussion: possible Ly$\alpha$ sources

Our quasar-Ly$\alpha$ emission cross-correlation measurements can be explained
by a model with Ly$\alpha$ emission associated with star-forming galaxies (§
3), and the Ly$\alpha$ forest-Ly$\alpha$ emission cross-correlation
measurements are also consistent with such an explanation (§ 4). The model,
however, does not provide details on the relation between Ly$\alpha$ emission
and galaxies, which we explore in this section.

As shown in Figure 8, the measured Ly$\alpha$ luminosity density, $\rho_{\rm
Ly\alpha}=6.6_{-3.1}^{+3.3}\times 10^{40}$ erg s-1 cMpc-3, computed from our
best-fit $\langle\mu_{\alpha}\rangle$ under the galaxy-dominated case. This
iceberg of Ly$\alpha$ emission can hardly be accounted for by Ly$\alpha$
emission from LAEs based on observed Ly$\alpha$ LFs, as shown in Figure 8 with
Ly$\alpha$ luminosity densities obtained from integrating the Ly$\alpha$ LF of
LAEs down to a low luminosity. For example, the value of $\rho_{\rm Ly\alpha}$
calculated by integrating the LAE LF at $z=2.5\pm 0.1$ in Sobral et al. (2018)
down to $1.75\times 10^{41}{\rm erg~{}s^{-1}}$ is $7.4^{+0.8}_{-0.7}\times
10^{39}{\rm erg~{}s^{-1}~{}cMpc^{-3}}$, only $\sim 12\%$ of our estimate. That
is, Ly$\alpha$ emission formally detected from LAEs is only the tip of the
iceberg.

Conversely, if we assume that all the Ly$\alpha$ photons detected in our work
are produced by star formation activities and neglect any dust effect on
Ly$\alpha$ emission, the implied SFRD $\rho_{\rm SFR}$ approximates the lower
bound of the dust-corrected cosmic $\rho_{\rm SFRD}$ determined by UV and IR
observations (see Figure 9).

There have to be some other sources responsible for the excessive Ly$\alpha$
emission. In this section, we explore two possible sources based on previous
observations and models: Ly$\alpha$ emission within an aperture centered on
star-forming galaxies, including LAEs and Lyman break galaxies (LBGs), with a
typical aperture of $2\arcsec$ in diameter in most NB surveys; Ly$\alpha$
emission outside the aperture usually missed for individual galaxies in NB
surveys, commonly called extended or diffuse Ly$\alpha$ halos. We name the two
components as inner and outer part of Ly$\alpha$ emission, respectively. For
the outer, diffuse Ly$\alpha$ halo component, we do not intend to discuss its
origin here (e.g., Zheng et al., 2011b; Lake et al., 2015) but adopt an
observation-motivated empirical model to estimate its contribution.

We argue that almost all star-forming galaxies produce Ly$\alpha$ emission,
and actually, significant emission may be originated from their halos. This
should contribute to the bulk of faint diffuse Ly$\alpha$ emission in the
Universe, as detected in this work.

###  5.1 Inner Part of Ly$\alpha$ Emission for UV-selected Star-forming
Galaxies

A large portion of LBGs exhibit Ly$\alpha$ emission, though their rest-frame
equivalent width (REW) might not satisfy the criteria for LAE selections
(Shapley et al., 2003; de La Vieuville et al., 2020) if measured with a
typical aperture of $2\arcsec$ in diameter in NB surveys. It is also detected
in deep stacks of luminous and massive LBGs (Steidel et al., 2011) and in
individual UV-selected galaxies in recent MUSE eXtremely Deep Field (MXDF)
observations (Kusakabe et al., 2022).

Dijkstra & Wyithe (2012) reported the Ly$\alpha$ REW distribution of $\sim
800$ $z\sim 3$ LBGs spectroscopically observed by Shapley et al. (2003) with
$1.4\arcsec$ slits, which can be described well by an exponential function.
This sample includes both Ly$\alpha$ emission (${\rm REW}>0$ Å) and Ly$\alpha$
absorption (${\rm REW}<0$ Å) within the central aperture. Combined with this
empirical model of Ly$\alpha$ ${\rm REW}$ distribution for star-forming
galaxies, we perform integration over the UV LF to obtain the corresponding
Ly$\alpha$ luminosity density,

$\displaystyle\rho^{\rm inner}_{\rm Ly\alpha}=$ $\displaystyle\int_{M_{\rm
UV,min}}^{M_{\rm UV,max}}[\langle L_{\alpha}(M_{\rm UV})\rangle\Phi^{e}_{\rm
UV}(M_{\rm UV})$ (25) $\displaystyle+\langle L^{\rm Abs}_{\alpha}(M_{\rm
UV})\rangle\Phi^{a}_{\rm UV}(M_{\rm UV})]dM_{\rm UV}$

where $\langle L_{\alpha}(M_{\rm UV})\rangle$ is the mean Ly$\alpha$
luminosity within the aperture of the ${\rm REW}>0$ Å population at a given UV
luminosity and $\langle L^{\rm Abs}_{\alpha}(M_{\rm UV})\rangle$ is the
absorption of the ${\rm REW}<0$ Å population making a negative contribution.
The function $\Phi^{e}_{\rm UV}$ is the UV LF for the ${\rm REW}>0$ Å
population, which is the overall UV LF $\Phi_{\rm UV}$ multiplied by the (UV
luminosity-dependent) fraction of such a population, and $\Phi^{a}_{\rm UV}$
for the ${\rm REW}<0$ Å population likewise. More details on the calculations
in our adopted model are presented in Appendix B.

We select five observed UV LFs around $z\approx 2.4$ from the literature
(Table 1), and calculate the corresponding Ly$\alpha$ luminosity densities,
which are shown in Table 2.

We note that the distribution of Ly$\alpha$ REW within the central aperture is
mainly determined by three factors: the intrinsic REW from photoionization and
recombination in the H II region of star-forming galaxies, the dust
extinction, and the scattering-induced escape fraction. The empirically
modelled Ly$\alpha$ REW distribution in Dijkstra & Wyithe (2012) we adopt
reflects the combination of the three factors.

Table 1: A compilation of the derived Schechter function parameters for the
galaxy UV LFs adopted in this work.

Source | $z$ |  $\lambda_{\rm UV}$1 (Å) | $M^{*}$ | $\Phi^{*}(10^{-3}{\rm cMpc}^{-3})$ | $\alpha$
---|---|---|---|---|---
Reddy & Steidel (2009) | 2.3 | 1700 | $-20.70\pm 0.11$ | $2.75\pm 0.54$ | $-1.73\pm 0.07$
Sawicki (2012) | 2.2 | 1700 | $-21.00\pm 0.50$ | $2.74\pm 0.24$ | $-1.47\pm 0.24$
Parsa et al. (2016) | 2.25 | 1700 | $-19.99\pm 0.08$ | $6.20\pm 0.77$ | $-1.31\pm 0.04$
Bouwens et al. (2015) | - | 1600 | $-20.89+0.12z$ | $0.48\times 10^{-0.19(z-6)}$ | $-1.85-0.09(z-6)$
extrapolation2 | 2.4 |  | -20.60 | 2.3 | -1.53
Parsa et al. (2016) | - | 1700 | $\frac{-35.4(1+z)^{0.524}}{1+(1+z)^{0.678}}$ | $-0.36z+2.8$ | $-0.106z-1.187$
extrapolation3 | 2.4 |  | -20.41 | 1.9 | -1.44

  * 1

Rest-frame UV wavelength where the UV LF is measured. Note that $\lambda_{\rm
UV}$ for Bouwens et al. (2015) is 1600 Å, while the empirical model in
Dijkstra & Wyithe (2012) as summarized in Appendix B.1 adopts 1700 Å. We just
assume that UV LFs are not sensitive to such a subtle difference in
$\lambda_{\rm UV}$.

  * 2

Extrapolation of the Schechter parameters of the UV LF to $z=2.4$ adopting the
best-fitting formula in Bouwens et al. (2015) for the redshift evolution.

  * 3

Extrapolation to $z=2.4$, based on the simple parametric fits to published
Schechter parameters in Parsa et al. (2016). Note that this fitting is meant
to illustrate the overall evolutionary trend, but not to indicate a best
estimate of true parameter evolution.

###  5.2 Outer Part of Ly$\alpha$ Emission from Galaxy Halos

As discussed before, many previous works have reported detections of extended
Ly$\alpha$ emission around high-redshift galaxies, either by discoveries of
Ly$\alpha$ halos/blobs around bright individual star-forming galaxies through
ultradeep exposures (Steidel et al., 2000; Matsuda et al., 2004, 2011;
Wisotzki et al., 2016; Leclercq et al., 2017; Kusakabe et al., 2022), or by
employing stacking analyses on large samples (Steidel et al., 2011; Matsuda et
al., 2012; Momose et al., 2014, 2016; Xue et al., 2017). Most extended
Ly$\alpha$-emitting halos are discovered around LAEs (Wisotzki et al., 2016;
Leclercq et al., 2017); they are also prevalent around non-LAEs, e.g., UV-
selected galaxies, due to a significant amount of cool/warm gas in their CGM
(Steidel et al., 2011; Kusakabe et al., 2022).

The cumulative fraction of the large-aperture Ly$\alpha$ flux, shown in Fig.10
of Steidel et al. (2011), indicates that a 2 arcsec aperture adopted by
typical deep narrow/medium-band LAE surveys could miss $\sim$50% Ly$\alpha$
emission for LBGs with net (positive) Ly$\alpha$ emission. Thus Equation (25)
could underestimate the total Ly$\alpha$ flux from ${{\rm REW}}>0$ Å galaxies
roughly by a factor of 2. For galaxies whose inner parts present net
Ly$\alpha$ absorption, the existence of extended Ly$\alpha$ halos has been
strongly confirmed by the sample with Ly$\alpha$ ${{\rm REW}}<0$ Å in Steidel
et al. (2011), whose radial SB profile outside 10 kpc is qualitatively similar
to that of the non-LAE sub-samples.

Given the above observational results, we adopt the reasonable model that all
star-forming galaxies, whether showing Ly$\alpha$ emission or absorption
within the central aperture, have Ly$\alpha$ emitting halos. Based on the
strong anti-correlation between Ly$\alpha$ luminosities of Ly$\alpha$ halos
and the corresponding UV magnitudes reported in Leclercq et al. (2017), we
assume that the Ly$\alpha$ luminosity from halos of galaxies with ${{\rm
REW}}<0$ Å depends on $M_{\rm UV}$ only. We further assume that it is equal to
the inner part originated from the ${{\rm REW}}>0$ Å galaxy population at a
given $M_{\rm UV}$ (Steidel et al., 2011). Therefore we express the total
contribution to the Ly$\alpha$ luminosity density from the outer part as

$\rho^{\rm outer}_{\rm Ly\alpha}=\int_{M_{\rm UV,min}}^{M_{\rm UV,max}}\langle
L_{\alpha}(M_{\rm UV})\rangle\Phi_{\rm UV}(M_{\rm UV})dM_{\rm UV},$ (26)

where $\Phi_{\rm UV}$ denotes the UV LF for the entire population (See
Appendix B).

Clearly, the total Ly$\alpha$ luminosity density should be $\rho_{\rm
Ly\alpha}^{\rm tot}=\rho_{\rm Ly\alpha}^{\rm inner}+\rho_{\rm Ly\alpha}^{\rm
outer}$. Note that the total Ly$\alpha$ luminosity density estimated from the
model is just a lower limit as discussed in Appendix B, since we (1) adopt a
constant scaling factor for the Dijkstra & Wyithe (2012) empirical model and
(2) use this empirical model that is designed for the ${\rm REW}>0$ Å
population to describe the ${\rm REW}<0$ Å one. A brief summary of the
estimated Ly$\alpha$ luminosity density is listed in Table 2 and Figure 12.
Revealed by Figure 12, the total Ly$\alpha$ luminosity density derived from
our model is consistent with our detection within 1$\sigma$ (or $\sim
1.3\sigma$ when using the $z=2.4$ UV LF in Parsa et al. 2016). We argue that,
star-forming galaxies, which contain the inner part of Ly$\alpha$ emission
that can be captured by the aperture photometry in deep NB surveys and the
outer part of Ly$\alpha$ emission from their halos, usually outside the
aperture, could produce sufficient Ly$\alpha$ emission to explain our
detection from the quasar-Ly$\alpha$ emission cross-correlation measurement.

Our derived $\rho_{\rm Ly\alpha}$ is higher than the result of Wisotzki et al.
(2018), who use MUSE observations of extended Ly$\alpha$ emission from LAEs to
infer a nearly 100% sky coverage of Ly$\alpha$ emision. The LAE sample they
use are selected from the Hubble Deep Field South (HDFS) and the Hubble Ultra
Deep Field (HUDF), a subset of LAEs whose Ly$\alpha$ LFs has been analyzed in
Drake et al. (2017a) and Drake et al. (2017b) (though the sample in Wisotzki
et al. (2018) contains a few additional LAEs). As shown in Figure 8,
$\rho_{\rm Ly\alpha}$ estimated in Drake et al. (2017a) is lower than ours,
too. Our result implies that Wisotzki et al. (2018) may underestimate the
Ly$\alpha$ sky coverage at a given SB level when simply focusing on LAEs and
ignoring the diffuse Ly$\alpha$ emission from faint UV-selected galaxies.

As shown in Figure 12, about half of the detected Ly$\alpha$ photons come from
the inner part of galaxies. By assuming that they all stem from star formation
activities, we estimate the escape fraction $f_{\rm esc}$ for these Ly$\alpha$
photons to be roughly $0.21_{-0.11}^{+0.21}$, where the cosmic intrinsic
Ly$\alpha$ luminosity density due to star formation is calculated based on the
cosmic SFRD shown in Figure 9, yielding $1.44^{+10.1}_{-6.1}\times 10^{41}{\rm
erg\,s^{-1}Mpc^{-3}}$. While the estimated $f_{\rm esc}$ appears consistent
with previous work within $1\sigma$ uncertainties (e.g., $\sim 10\%$ in Chiang
et al. 2019), we emphasize that the galaxy population involved in our
modelling is different from LAEs in typical NB surveys. We include galaxies
with low Ly$\alpha$ REW usually not identified as LAEs, which boost our
estimate for $f_{\rm esc}$ compared with LAE-derived ones.

Table 2: Model Ly$\alpha$ luminosity density $\rho_{\rm Ly\alpha}$ by
integrating UV LFs (from $M_{\rm UV,min}=-24$ to $M_{\rm UV,max}=-12$) based
on Schechter functions from various sources as in Table 1. See Section 5.1 and
5.2 for more details.

[t]  Source $z$ $\rho_{\rm Ly\alpha}$ ($10^{40}$ erg s-1 cMpc-3) inner1 outer2
total3 Reddy & Steidel (2009) 2.3 3.82 4.17 8.00 Sawicki (2012) 2.2 2.10 2.71
4.81 Parsa et al. (2016) 2.25 1.83 2.05 3.88 Bouwens et al. (2015) 2.4 1.61
1.87 3.49 extrapolation4 Parsa et al. (2016) 2.4 0.97 1.12 2.09 extrapolation

  * 1

Ly$\alpha$ luminosity density from emission that would be captured within an
aperture of $2\arcsec$ in diameter, computed from Equation (25). Galaxies with
Ly$\alpha$ ${{\rm REW}}>0$ Å contribute a positive part and the Ly$\alpha$
${{\rm REW}}<0$ Å population contribute a negative one.

  * 2

Ly$\alpha$ luminosity density from emission outside the 2$\arcsec$ aperture
for all galaxies, i.e., the diffuse Ly$\alpha$ halo component, computed from
Equation (26). At a given UV luminosity, we assume that the populations with
central ${\rm REW}>0$ Å and ${\rm REW}<0$ Å have the same diffuse halo
Ly$\alpha$ luminosity, which is set to be the same as that from the inner part
of the ${\rm REW}>0$ Å population in our model based on the results in Steidel
et al. (2011).

  * 3

Total Ly$\alpha$ luminosity density contributed by the three components
discussed above.

  * 4

Same as in Table 1.

Figure 12: Ly$\alpha$ luminosity density computed in our model by integrating
different observed UV LFs. Different colors denote the inner and outer
Ly$\alpha$ parts, as described in Table 2 in details. The model Ly$\alpha$
luminosity densities are compared with that inferred from our quasar-
Ly$\alpha$ emission cross-correlation measurements, $6.6\times 10^{40}$ erg
s-1 cMpc-3 (solid) with 1$\sigma$ errorbars (dashed) of $\pm 3.2\times
10^{40}$ erg s-1 cMpc-3.

##  6 Summary and conclusion

In this work, we have performed a cross-correlation analysis of the SDSS
BOSS/eBOSS LRG residual spectra at wavelengths $\lambda$ = 3647–5471Å and DR16
quasars at a redshift range of $2<z<3.5$. This enables a measurements of the
cross-correlation between quasar position and Ly$\alpha$ emission intensity
(embedded in the residual LRG spectra) at a median redshift $z\sim 2.4$. The
Ly$\alpha$ SB profile around quasars is obtained by projecting our cross-
correlation results into a pseudo-narrow band, and the truncated forms of the
monopole, quadrupole and hexadecapole of the quasar-Ly$\alpha$ emission cross-
correlation are computed by discarding small-scale signals within
$r_{\perp}<4h^{-1}{\rm cMpc}$.

Our work improves upon that in Croft et al. (2018) by making use of the final
SDSS-IV release of LRG spectra and quasar catalog. While our Ly$\alpha$ SB
profile measurements are consistent with that in Croft et al. (2018), our
inferred large-scale clustering amplitude is about 2.2 times lower. Although
the absolute uncertainty in our work is about 25% lower, the lower clustering
amplitude leads to a larger fractional uncertainty. This is a reflection of
our more rigorous treatment to possible contaminated fibers and our exclusion
of the small-scale signals in modelling the multipoles. With this lower
amplitude, our measured Ly$\alpha$ forest-Ly$\alpha$ emission cross-
correlation can also be consistently explained.

Like Croft et al. (2018), on sub-Mpc scales the obtained Ly$\alpha$ SB forms a
natural extrapolation of that observed from the luminous Ly$\alpha$ blobs on
smaller scales (Borisova et al., 2016; Cai et al., 2019). Unlike Croft et al.
(2018), we find that the amplitudes of the large-scale Ly$\alpha$ SB and
quasar-Ly$\alpha$ emission cross-correlation cannot result from the Ly$\alpha$
emission around quasars, as this would require the average Ly$\alpha$
luminosity of quasars to be about two orders of magnitude higher than observed
given their rather low number density.

To figure out the most possible sources that contribute to the detected
Ly$\alpha$ signals, we construct a simple analytical model, which combines the
SB profile and multipole measurements. The inferred Ly$\alpha$ luminosity
density, $6.6_{-3.1}^{+3.3}\times 10^{40}{\rm erg\,s^{-1}cMpc^{-3}}$, is much
higher than those from integrating the Ly$\alpha$ LFs of LAEs. We fix the
luminosity weighted bias of galaxies $b_{\alpha}$ to be 3 in our modelling,
which turns out to be a good estimate. But bear in mind that the luminosity
density scales with 3/$b_{\alpha}$ if $b_{\alpha}$ deviates from that value.
Our model rules out the possibility that the diffuse emission is due to
reprocessed energy from the quasars themselves, and support the hypothesis
that star-forming galaxies clustered around are responsible for the detected
signal. For the Ly$\alpha$ forest-Ly$\alpha$ emission cross-correlation, the
prediction from our model matches the measurement, although the current
measurement is consistent with a null detection given the low signal-to-noise
ratio. We argue that most star-forming galaxies exhibit Ly$\alpha$ emission.
These include galaxy populations with either Ly$\alpha$ emission or Ly$\alpha$
absorption at the center, while both populations have diffuse Ly$\alpha$
emitting halos, which are usually missed in individual LAEs from deep narrow-
band surveys. Our estimates based on the empirical model of Dijkstra & Wyithe
(2012) and the observed UV LFs of star-forming galaxies are able to match the
Ly$\alpha$ luminosity density inferred from our cross-correlation
measurements. The picture is supported by stacked analysis from NB surveys
(e.g., Steidel et al., 2011) and by the IFU observations of Ly$\alpha$
emission associated with UV-selected galaxies (e.g., Kusakabe et al., 2022).

Our work shows an enormous promise of Ly$\alpha$ intensity mapping as a probe
of large scale structure. One can also utilize this technique to explore the
intensity of other spectral lines, once a larger data set is provided. The
next-generation cosmological spectroscopic survey, the ongoing Dark Energy
Spectroscopic Instrument (DESI; DESI Collaboration et al. 2016), will enlarge
the galaxy/quasar survey volume at least by an order of magnitude compared to
SDSS BOSS/eBOSS. We expect the intensity mapping technique carried out in DESI
will bring us new insights into the Universe. Deep surveys of Ly$\alpha$
emission around star-forming galaxies, especially the UV-selected population
(e.g., Kusakabe et al., 2022) will shed light on the intensity mapping
measurements and provide inputs for building the corresponding model.
Moreover, more realistic modelling of physical processes such as radiative
transfer and quasar proximity effect should be considered to advance our
understanding of the Ly$\alpha$ emission iceberg in the Universe.

We thank Kyle Dawson, Rupert Croft, and Coast Zhang for useful discussions.
X.L. and Z.C. are supported by the National Key R&D Program of China (grant
No.2018YFA0404503) and the National Science Foundation of China (grant No.
12073014). Z.Z. is supported by NSF grant AST-2007499. Funding for the Sloan
Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the
U.S. Department of Energy Office of Science, and the Participating
Institutions. SDSS-IV acknowledges support and resources from the Center for
High Performance Computing at the University of Utah. The SDSS website is
www.sdss.org. SDSS-IV is managed by the Astrophysical Research Consortium for
the Participating Institutions of the SDSS Collaboration including the
Brazilian Participation Group, the Carnegie Institution for Science, Carnegie
Mellon University, Center for Astrophysics — Harvard & Smithsonian, the
Chilean Participation Group, the French Participation Group, Instituto de
Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the
Physics and Mathematics of the Universe (IPMU) / University of Tokyo, the
Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz
Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie
(MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-
Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical
Observatories of China, New Mexico State University, New York University,
University of Notre Dame, Observatário Nacional / MCTI, The Ohio State
University, Pennsylvania State University, Shanghai Astronomical Observatory,
United Kingdom Participation Group, Universidad Nacional Autónoma de México,
University of Arizona, University of Colorado Boulder, University of Oxford,
University of Portsmouth, University of Utah, University of Virginia,
University of Washington, University of Wisconsin, Vanderbilt University, and
Yale University.

##  Appendix A Correcting Measurement Systematics

Dealing with possible contamination is a difficult problem in all intensity
mapping (IM) experiments. Since the expected signals in our measurement have
gone beyond the detection capability of any current instruments, it is crucial
to remove possible systematics. In this section we discuss three main sources
of potential contamination: cross-talk effect among spectra in adjacent
fibers, correlation at $r_{\|}=0$ for pixel-pixel correlation, and spurious
signal on larger scales, and then demonstrate that we have removed them
carefully from our measurement.

###  A.1 Quasar Stray Light Contamination

The BOSS/eBOSS spectrograph has 1,000 fibers per plate, which disperse light
onto the same 4096-column CCD. Light from one fiber would possibly leak into
the extraction aperture for another fiber, but the level of this light
contamination is negligible in the SDSS data reduction pipeline. However, our
intensity-mapping technique reaches far beyond the instrument capability
($\sim 10^{-17}{\rm erg\,s^{-1}cm^{-2}\textup{\AA}^{-1}arcsec^{-2}}$), so this
light contamination should be treated cautiously. When cross-correlating
quasar-LRG spectrum pixels, the cross-correlation between quasars and its leak
into LRG spectra will lead to a contamination 3-4 orders of magnitude higher
than the targetted Ly$\alpha$ signals, due to the bright and broad Ly$\alpha$
features of quasars.

In Croft et al. (2016), quasar stray light contamination is removed through
discarding any quasar-LRG spectrum pixel pairs once on the CCD the quasar is
within five fibers apart from the LRG, i.e., $\Delta{\rm fiber}<5$. Moreover,
Croft et al. (2018) reported that the remaining quasar stray light would still
lead to contamination as a result of quasar clustering effect: an LRG fiber
with $\Delta_{\rm fiber}\geq 5$ away from a quasar fiber may be contaminated
by another quasar, and if not corrected for, the cross-correlation between
Ly$\alpha$ emission in the LRG fiber with the first quasar would have the
quasar clustering signal imprinted. They find that the quasar clustering
effect would reach 50% of the signal in Croft et al. (2016) on scales of
$|r_{\perp}|<10h^{-1}{\rm Mpc}$ and $|r_{\|}|<10h^{-1}{\rm Mpc}$. Croft et al.
(2018) correct such an effect by generating a set of mock spectra, which
contain quasar contaminating light only, and performing the same cross-
correlation procedure to measure the intensity of clustering. Then the
clustering signal from the mock is subtracted from their originally measured
signals.

The key of the algorithm in Croft et al. (2018) is to estimate the light
leakage fraction so that cross-correlation of mock spectra can precisely
reproduce the quasar clustering effect. The fraction measured in Croft et al.
(2016) is no longer applicable to our sample spectra, however, due to the
recent updates on DR16 optical spectra
pipeline333https://www.sdss.org/dr16/spectro/pipeline/. In our measurement, to
be conservative, instead, we exclude any LRG fiber once it is within 5 fiber
or less apart from a quasar fiber, and this fiducial sample selection will
remove both quasar stray light contamination and quasar clustering effect
simultaneously. We also repeat the algorithm introduced in Croft et al.
(2018), removing the quasasr clustering systematics by subtracting the cross-
correlation pattern produced by mock spectra, and then measure the
corresponding multipoles and SB profiles. To ensure the robustness of our
fiducial sample selection, i.e., excluding all LRG fibers of $\Delta{\rm
fiber}<5$, we perform the same fitting procedure mentioned in Section 3.4
under the galaxy-dominated assumption for test cases with various sample
selections. A comparison of results with differently selected samples is
demonstrated in Figure 13. The result of $\Delta{\rm fiber}<5$ LRG exclusion
are in fact consistent with that of $\Delta{\rm fiber}<8$ and $\Delta{\rm
fiber}<10$ within $1\sigma$, implying that Four fiducial selection can remove
the contamination well. It in general accords with the result of
$\Delta$fiber$<5$ pair exclusion, i.e., the method used in Croft et al.
(2018), though there is a tiny offset in best-fit $\beta_{\alpha}$ and the
uncertainties of the three parameters from the latter are smaller.

Figure 13: A test for the robustness of our fiducial sample selection, by
performing joint fits to the multipoles and SB profiles of quasar-Ly$\alpha$
emission cross-correlation after different sample selections. The model
parameters are the same as in Figure 7. Fiber Cut 5 refers to the sample that
any LRG fiber with $\Delta{\rm fiber}<5$ from a quasar fiber is excluded,
which is the fiducial sample for our measurements.; Fiber Cut 8 and Fiber Cut
10 refers to $\Delta{\rm fiber}<8$ and $\Delta{\rm fiber}<10$, respectively.
Pair Cut 5 refers to the sample in which a quasar-LRG spectrum pixel pair is
excluded if they satisfies $\Delta{\rm fiber}<5$, i.e. the method adopted in
Croft et al. (2018).

###  A.2 $r_{\|}=0$ Correlation for Pixel-pixel Pairs

In measuring the Ly$\alpha$ forest-Ly$\alpha$ emission cross-correlation, we
need to remove an artifact of correlation around $r_{\|}=0$, introduced by the
spectra pipeline.

In BOSS/eBOSS, each half-plate has $\sim$500 fibers (450 science fibers and
$\sim$ 40 sky fibers) with two spectrographs. The sky-subtraction for
individual spectra is done indepently for each spectrograph. Poisson
fluctuations in sky spectra will induce correlations in those spectra obtained
with the same spectrograph at the same observed wavelength (Bautista et al.,
2017; du Mas des Bourboux et al., 2019, 2020). Namely, a positive correlation
is expected for spectrum pixel pairs on the same half-plates at
$\Delta\lambda_{\rm obs}=0$, leading to an excess correlation in $r_{\|}=0$
bins. Furthermore, the continuum fitting procedure designed for Ly$\alpha$
forest transmission fields may smooth the excess correlation at $r_{\|}=0$,
extending it to larger $|r_{\|}|$. Therefore we reject Ly$\alpha$ forest-
Ly$\alpha$ emission pixel pairs once they are observed on the same-half plate.

To evaluate how this same-spectrograph induced systematics would contaminate
the signals and whether we have fully removed it, we perform a measurements of
cross-correlation between Ly$\alpha$ forest transmission pixels and Ly$\alpha$
emission pixels, as a function of their observed wavelength separations
$\Delta\lambda_{\rm obs}$ and transverse separations $\Delta\theta$,

$\xi_{f\alpha}(\Delta\lambda_{\rm
obs},\Delta\theta)=\frac{1}{\sum_{i=1}^{N}w_{ri,\alpha}w_{ri,f}}\sum_{i=1}^{N}w_{ri,\alpha}w_{ri,f}\Delta_{\mu,ri}(\lambda_{\rm
obs},\theta)\delta_{f,{ri}}(\lambda_{\rm obs}+\Delta\lambda_{\rm
obs},\theta+\Delta\theta),$ (A1)

where $\Delta\theta$ can be easily converted to the transverse comoving
separation at $z=2.41$ by $R_{\perp}=\Delta\theta\cdot D_{\rm C}(z=2.41)$. The
cross-correlation results of different-half-plate pixel pairs, same-half-plate
pixel pairs and all pixel pairs without selection preference are shown in
Figure 14. The contamination at $\lambda_{\rm obs}=0$ for same-half-plate
pairs reaches several times $10^{-21}{\rm
erg\,s^{-1}cm^{-2}\textup{\AA}^{-1}arcsec^{-2}}$ (middle panel), even stronger
than the targeted signals, stressing the necessity of rejecting same-half-
plate pixel pairs. While the contamination is largely removed when we only use
different-half-plate pixel pairs (left panel), there still appears to be a
residual weak correlation at $\Delta\lambda_{\rm obs}\sim 0$, not expected
from pure sky-subtraction effects. The exact source of such a weak correlation
at $\Delta\lambda_{\rm obs}\sim 0$ may be related to some details in the
processing procedure in the spectra pipeline. To proceed, we adopt a
conservative method to remove the effect of this weak correlation by
discarding any signal within $|r_{\|}|<4$cMpc, at the expense of slightly
reducing the signal-to-noise ratio of our measurement.

Figure 14: Ly$\alpha$ forest-Ly$\alpha$ emission pixel-pixel correlation as a
function of observed wavelength separations $\Delta\lambda_{\rm obs}$ and
transverse separations at a redshift of 2.41. Different half-plate pairs
refers to the sample in which the selected pixel pairs reside on different
half-plates; Same half-plate pairs refers to the sample in which all pairs are
on the same half-plates. All pairs refers to the sample without preference for
the plates/fibers. The bottom panels show the correlations averaged within
$2\textup{\AA}$ at $\Delta\lambda_{\rm obs}\sim 0$.

###  A.3 Large-scale Correction

As discussed in Croft et al. (2016), one may find non-zero cross-correlation
for large pair separation with no physical significance we concern about. We
correct this spurious signal by subtracting the average correlation over
80–400$h^{-1}{\rm Mpc}$ along both the line-of-sight and orthogonal
directions, following the method described in Croft et al. (2016).

##  Appendix B Model for Ly$\alpha$ Luminosity Density Contributed by Star-
forming Galaxies

In our model, star-forming galaxies dominate the Ly$\alpha$ luminosity
density. We first review the Dijkstra & Wyithe (2012) model for the REW
distribution of Ly$\alpha$ emission from an inner aperture around star-forming
galaxies. With such an REW distribution, we present our model of Ly$\alpha$
luminosity density from contributions of Ly$\alpha$ emission within the inner
aperture and from the outer halo.

###  B.1 Model for Ly$\alpha$ Rest-frame Equivalent Width Distribution of UV-
selected Galaxies

Dijkstra & Wyithe (2012) modelled the conditional probability density function
(PDF) for the REW of Ly$\alpha$ emission (from the central aperture around
LBGs) using an exponential function whose scaling factor ${{\rm REW}}_{c}$
depends on $M_{\rm UV}$ and $z$,

$P\left({{\rm REW}}\mid
M_{\mathrm{UV}}\right)=\begin{cases}\mathcal{N}\exp\left(-\frac{\mathrm{REW}}{\operatorname{REW}_{\mathrm{c}}}\right),&{{\rm
REW}}\in\left(x_{\min},x_{\max}\right)\\\ 0,&\text{ otherwise }\end{cases}$
(B1)

where $\mathcal{N}$ denotes a normalization constant. The choice of the
normalization factor $\mathcal{N}$ allows that all drop-out galaxies have
$x_{\min}\leq{{\rm REW}}\leq x_{\max}$,

$\mathcal{N}^{-1}=\operatorname{REW}_{c}\left[\exp\left(-\frac{x_{\min}}{\operatorname{REW}_{\mathrm{c}}\left(M_{\mathrm{UV}}\right)}\right)-\exp\left(-\frac{x_{\max}}{\operatorname{REW}_{\mathrm{c}}\left(M_{\mathrm{UV}}\right)}\right)\right].$
(B2)

To match the $M_{\rm UV}$-dependence of the observed fraction of LAEs (${{\rm
REW}}>50$Å) in drop-out galaxies, they fixed $x_{\max}=300$ and assumed
$x_{\min}\equiv-a_{1}$ (both in units of Å),

$a_{1}=\left\\{\begin{array}[]{lcl}20&&M_{\mathrm{UV}}<-21.5\\\
20-6\left(M_{\mathrm{UV}}+21.5\right)^{2}&&-21.5\leq M_{\mathrm{UV}}<-19\\\
-17.5&&\text{ other. }\end{array}\right.$ (B3)

In their fiducial model, ${{\rm REW}_{c}}$ evolves with $M_{\rm UV}$ and $z$,

${{\rm REW}}_{c}(M_{\rm UV},z)={{\rm REW}}_{c,0}+\mu_{1}(M_{\rm
UV}+21.9)+\mu_{2}(z-4),$ (B4)

where the best-fitting parameters are ${{\rm REW}}_{c,0}=23$Å, $\mu_{1}=7$Å,
$\mu_{2}=6$Å. Note that the fitting formula applies only in the observed range
of UV magnitudes and the evolution is frozen for $M_{\rm UV}>-19$. However, in
our analysis we adopt a constant ${{\rm REW}}_{c}=22$Å, which depicts the REW
distribution of the 400 brightest LBG sample of Shapley et al. (2003) well but
underpredicts the faint-end LAE fraction, as disccussed in Appendix A1 of
Dijkstra & Wyithe (2012). With this constant ${{\rm REW}}_{c}$, we would
underestimate the Ly$\alpha$ luminosity contributed by UV-faint galaxies, and
the total estimated Ly$\alpha$ emission would be a lower limit.

The Ly$\alpha$ luminosity at a given ${\rm REW}$ and ${{\rm UV}}$ luminosity
can be expressed as

$L_{\alpha}\left({{\rm REW}},M_{\rm UV}\right)=L_{\rm
UV,\nu}\left(\nu_{\alpha}/\lambda_{\alpha}\right)\left(\lambda_{\mathrm{UV}}/\lambda_{\alpha}\right)^{-\beta-2}\cdot{{\rm
REW}},$ (B5)

with the absolute AB magnitude $M_{\rm UV}=-2.5\log[{L_{\rm UV,\nu}}/({\rm
erg\,s^{-1}Hz^{-1}})]+51.6$. The parameter $\beta$ characterizes the slope of
the UV continuum, such that $L_{\rm UV,\lambda}=\nu L_{\rm
UV,\nu}/\lambda\propto\lambda^{\beta}$. We adopt $\lambda_{\rm UV}=1700$Å and
fix $\beta=-1.7$ as in Dijkstra & Wyithe (2012). The adopted wavelength is the
same as in the UV LF measurements (Table 1), except for the Bouwens et al.
(2015) UV LF (measured at 1600Å). In our calculation, we ignore the slight
wavelength shift in the Bouwens et al. (2015) UV LF, as the effect in the UV
luminosity computation is less than 2%.

###  B.2 Model for the Inner and Outer Ly$\alpha$ Emission Component

We separate star-forming galaxies into two populations based on the case of
Ly$\alpha$ radiation within the central 2$\arcsec$ aperture, one with
Ly$\alpha$ emission (${\rm REW}>0$) and one with Ly$\alpha$ absorption (${\rm
REW}<0$). We can express the corresponding UV LFs as

$\Phi^{e}_{\rm UV}(M_{\rm UV})=\frac{\int_{0}^{+\infty}P({\rm REW}\mid M_{\rm
UV})d{\rm REW}}{\int_{-\infty}^{+\infty}P({\rm REW}\mid M_{\rm UV})d{\rm
REW}}\,\Phi_{\rm UV}(M_{\rm UV})$ (B6)

for the ${\rm REW}>0$ population and

$\Phi^{a}_{\rm UV}(M_{\rm UV})=\frac{\int_{-\infty}^{0}P({\rm REW}\mid M_{\rm
UV})d{\rm REW}}{\int_{-\infty}^{+\infty}P({\rm REW}\mid M_{\rm UV})d{\rm
REW}}\,\Phi_{\rm UV}(M_{\rm UV})$ (B7)

for the ${\rm REW}<0$ population, where $P({\rm REW}\mid M_{\rm UV})$ is the
${\rm REW}$ distribution for galaxies with UV luminosity $M_{\rm UV}$.
Clearly, by construction, $\Phi^{e}_{\rm UV}+\Phi^{a}_{\rm UV}=\Phi_{\rm UV}$.
Note that we formally use $-\infty$ and $+\infty$ for clarity, while the true
cutoff thresholds are encoded in $P({\rm REW}\mid M_{\rm UV})$, which takes
the form of Equation (B1) if adopting the Dijkstra & Wyithe (2012) model.

The mean Ly$\alpha$ luminosity within the 2$\arcsec$ aperture of the ${\rm
REW}>0$ population at a given UV luminosity is

$\langle L_{\alpha}(M_{\rm UV})\rangle=\frac{\int_{0}^{+\infty}L_{\alpha}({\rm
REW},M_{\rm UV})P({\rm REW}\mid M_{\rm UV})d{\rm
REW}}{\int_{0}^{+\infty}P({\rm REW}\mid M_{\rm UV})d{\rm REW}}.$ (B8)

where $L_{\alpha}({\rm REW},M_{\rm UV})$ can be calculated through Equation
(B5). Figure 15 presents the evolution of $\langle L_{\alpha}(M_{\rm
UV})\rangle$, $\Phi_{\rm UV}^{a}$ and $\Phi_{\rm UV}^{e}$ with $M_{\rm UV}$ in
our model. We also show the expected Ly$\alpha$ luminosity for the SFR
associated with the UV luminosity, calculated through the relations that SFR
of 1$M_{\odot}{\rm yr^{-1}}$ corresponds to UV luminosity $L_{\nu}=1.4\times
10^{-28}{\rm erg~{}s^{-1}Hz^{-1}}$ and Ly$\alpha$ luminosity
$L_{\alpha}=1.1\times 10^{42}{\rm erg~{}s^{-1}}$. It is much higher than our
modelled Ly$\alpha$ luminosity, consistent with the measurements in Figure 9.

In addition, the net absorption from the ${\rm REW}<0$ population will also
make a negative contribution. The ‘absorbed’ luminosity could be descibed as

$\langle L_{\alpha}^{\rm Abs}(M_{\rm
UV})\rangle=\frac{\int_{-\infty}^{0}L_{\alpha}({\rm REW},M_{\rm UV})P({\rm
REW}\mid M_{\rm UV})d{\rm REW}}{\int_{-\infty}^{0}P({\rm REW}\mid M_{\rm
UV})d{\rm REW}},$ (B9)

which would yield a negative value.

The contribution to the Ly$\alpha$ luminosity density from the inner part
comes from the emission of the ${\rm REW}>0$ population and the absorption of
the ${\rm REW}<0$ population, which is

$\rho^{\rm inner}_{\rm Ly\alpha}=\int_{M_{\rm UV,min}}^{M_{\rm
UV,max}}\left[\langle L_{\alpha}(M_{\rm UV})\rangle\Phi^{e}_{\rm UV}(M_{\rm
UV})+\langle L^{\rm Abs}_{\alpha}(M_{\rm UV})\rangle\Phi^{a}_{\rm UV}(M_{\rm
UV})\right]dM_{\rm UV}.$ (B10)

In our model the negative absorption component is actually insignificant
compared to the emission one, with the former being about 1–4% of the latter
depending on the adopted UV LF.

Based on the finding in Steidel et al. (2011), we assume that the Ly$\alpha$
luminosity in the diffuse halo component is the same as that from the central
aperture in the ${\rm REW}>0$ population and that the diffuse component in the
${\rm REW}<0$ population takes the same value at a given UV luminosity. Then
the contribution from the outer part Ly$\alpha$ emission of the ${\rm REW}>0$
population has the same expression as in the above equation, while that from
the ${\rm REW}<0$ population is obtained by replacing $\Phi^{e}_{\rm UV}$ with
$\Phi^{a}_{\rm UV}$. The total outer part contribution from Ly$\alpha$ halos
is then

$\rho^{\rm outer}_{\rm Ly\alpha}=\int_{M_{\rm UV,min}}^{M_{\rm UV,max}}\langle
L_{\alpha}(M_{\rm UV})\rangle\Phi_{\rm UV}(M_{\rm UV})dM_{\rm UV}.$ (B11)

We adopt $M_{\rm UV,min}=-24$ and $M_{\rm UV,max}=-12$ in our calculation.

The outer part Ly$\alpha$ emission can have contributions from satellite
galaxies in high-mass halos (e.g., Momose et al., 2016; Lake et al., 2015;
Mitchell et al., 2021), while the UV LF used to compute the inner part
Ly$\alpha$ emission should already include the satellite population.
Therefore, in our model there is a possibility of double-counting the
contribution of Ly$\alpha$ emission from the satellites. From halo modelling
of LBG clustering, Cooray & Ouchi (2006) find that the contribution from
satellites to the UV LF is at a level of $\sim 10^{-3}$–$10^{-2}$ over a wide
luminosity range and that it becomes even lower at the faint end ($M_{\rm
UV}>-17$). A similar result is also obtained by Jose et al. (2013). These
empirical results suggest that the contribution of satellite galaxies to the
total cosmic Ly$\alpha$ luminosity density is negligible, and we simply ignore
the effect induced by possibly double-counting satellites here.

Note that our model is just a rough esitmate of the total Ly$\alpha$
luminosity, with systematics arising from both the Ly$\alpha$ REW PDF and UV
LFs. For example, the Dijkstra & Wyithe (2012) REW PDF may underpredict the
number of large REW systems, leading to an underestimate of the total
Ly$\alpha$ luminosity. On the other hand, the modelled REW PDF may not
describe the number of galaxies with net absorption very well. However, these
uncertainties would not change our main claim significantly. Future
observations of UV luminosity dependent Ly$\alpha$ REW distribution and
measurements of UV LFs are expected to improve the modelling.

Figure 15: Left: Mean Ly$\alpha$ luminosity $\langle L_{\alpha}(M_{\rm
UV})\rangle$ within the 2$\arcsec$ aperture of the ${\rm REW}>0$Å population
as a function of the UV magnitude $M_{\mathrm{UV}}$, as presented in Equation
B8. The gray dotted lines denote the turning points of $a_{1}$ as expressed in
Equation B3. The green dashed line denotes the expected Ly$\alpha$ luminosity
for the SFR associated with the UV luminosity. Right: UV LF of the ${\rm
REW}<0$Å population $\Phi_{\rm UV}^{a}$, the ${\rm REW}>0$Å population
$\Phi_{\rm UV}^{e}$, and the entire population $\Phi_{\rm UV}$ as a function
of $M_{\rm UV}$, as presented in Equation B6 and B7. We take the Reddy &
Steidel (2009) UV LF as an example. The gray dotted lines denotes one of the
turning points of $a_{1}$ (Equation B3), where $a_{1}=0$ and ${\rm REW}$ keeps
larger than 0 as $M_{\rm UV}$ increases. That is, we assume that there is no
${\rm REW}<0$Å population over this $M_{\rm UV}$ range, which will lead to an
underestimation of the total Ly$\alpha$ luminosity.

## References

  * Ahumada et al. (2020) Ahumada, R., Prieto, C. A., Almeida, A., et al. 2020, ApJS, 249, 3, doi: 10.3847/1538-4365/ab929e
  * Anderson et al. (2018) Anderson, C. J., Luciw, N. J., Li, Y. C., et al. 2018, MNRAS, 476, 3382, doi: 10.1093/mnras/sty346
  * Arrigoni Battaia et al. (2016) Arrigoni Battaia, F., Hennawi, J. F., Cantalupo, S., & Prochaska, J. X. 2016, ApJ, 829, 3, doi: 10.3847/0004-637X/829/1/3
  * Arrigoni Battaia et al. (2018) Arrigoni Battaia, F., Hennawi, J. F., Prochaska, J. X., et al. 2018, MNRAS, 482, 3162, doi: 10.1093/mnras/sty2827
  * Bacon et al. (2010) Bacon, R., Accardo, M., Adjali, L., et al. 2010, in Ground-based and Airborne Instrumentation for Astronomy III, ed. I. S. McLean, S. K. Ramsay, & H. Takami, Vol. 7735, International Society for Optics and Photonics (SPIE), 131 – 139, doi: 10.1117/12.856027
  * Bacon et al. (2021) Bacon, R., Mary, D., Garel, T., et al. 2021, A&A, 647, A107, doi: 10.1051/0004-6361/202039887
  * Bautista et al. (2017) Bautista, J. E., Busca, N. G., Guy, J., et al. 2017, A&A, 603, A12, doi: 10.1051/0004-6361/201730533
  * Behroozi et al. (2019) Behroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. 2019, MNRAS, 488, 3143, doi: 10.1093/mnras/stz1182
  * Blomqvist et al. (2019) Blomqvist, M., du Mas des Bourboux, H., Busca, N. G., et al. 2019, A&A, 629, A86, doi: 10.1051/0004-6361/201935641
  * Bolton et al. (2012) Bolton, A. S., Schlegel, D. J., Aubourg, É., et al. 2012, AJ, 144, 144, doi: 10.1088/0004-6256/144/5/144
  * Borisova et al. (2016) Borisova, E., Cantalupo, S., Lilly, S. J., et al. 2016, ApJ, 831, 39, doi: 10.3847/0004-637X/831/1/39
  * Bouwens et al. (2015) Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2015, ApJ, 803, 34, doi: 10.1088/0004-637X/803/1/34
  * Brodzeller & Dawson (2022) Brodzeller, A., & Dawson, K. 2022, AJ, 163, 110, doi: 10.3847/1538-3881/ac4600
  * Cai et al. (2017) Cai, Z., Fan, X., Bian, F., et al. 2017, ApJ, 839, 131, doi: 10.3847/1538-4357/aa6a1a
  * Cai et al. (2018) Cai, Z., Hamden, E., Matuszewski, M., et al. 2018, ApJ, 861, L3, doi: 10.3847/2041-8213/aacce6
  * Cai et al. (2019) Cai, Z., Cantalupo, S., Prochaska, J. X., et al. 2019, ApJS, 245, 23, doi: 10.3847/1538-4365/ab4796
  * Cannon et al. (2006) Cannon, R., Drinkwater, M., Edge, A., et al. 2006, MNRAS, 372, 425, doi: 10.1111/j.1365-2966.2006.10875.x
  * Cantalupo et al. (2014) Cantalupo, S., Arrigoni-Battaia, F., Prochaska, J. X., Hennawi, J. F., & Madau, P. 2014, Nature, 506, 63, doi: 10.1038/nature12898
  * Cantalupo et al. (2008) Cantalupo, S., Porciani, C., & Lilly, S. J. 2008, ApJ, 672, 48, doi: 10.1086/523298
  * Cantalupo et al. (2005) Cantalupo, S., Porciani, C., Lilly, S. J., & Miniati, F. 2005, ApJ, 628, 61, doi: 10.1086/430758
  * Chiang et al. (2019) Chiang, Y.-K., Ménard, B., & Schiminovich, D. 2019, ApJ, 877, 150, doi: 10.3847/1538-4357/ab1b35
  * Cooray & Ouchi (2006) Cooray, A., & Ouchi, M. 2006, MNRAS, 369, 1869, doi: 10.1111/j.1365-2966.2006.10437.x
  * 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
  * de La Vieuville et al. (2019) de La Vieuville, G., Bina, D., Pello, R., et al. 2019, A&A, 628, A3, doi: 10.1051/0004-6361/201834471
  * de La Vieuville et al. (2020) de La Vieuville, G., Pelló, R., Richard, J., et al. 2020, A&A, 644, A39, doi: 10.1051/0004-6361/202037651
  * 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
  * Dijkstra & Wyithe (2012) Dijkstra, M., & Wyithe, J. S. B. 2012, MNRAS, 419, 3181, doi: 10.1111/j.1365-2966.2011.19958.x
  * Drake et al. (2017a) Drake, A. B., Garel, T., Wisotzki, L., et al. 2017a, A&A, 608, A6, doi: 10.1051/0004-6361/201731431
  * Drake et al. (2017b) Drake, A. B., Guiderdoni, B., Blaizot, J., et al. 2017b, MNRAS, 471, 267, doi: 10.1093/mnras/stx1515
  * du Mas des Bourboux et al. (2019) du Mas des Bourboux, H., Dawson, K. S., Busca, N. G., et al. 2019, ApJ, 878, 47, doi: 10.3847/1538-4357/ab1d49
  * du Mas des Bourboux et al. (2020) du Mas des Bourboux, H., Rich, J., Font-Ribera, A., et al. 2020, ApJ, 901, 153, doi: 10.3847/1538-4357/abb085
  * Eisenstein et al. (2001) Eisenstein, D. J., Annis, J., Gunn, J. E., et al. 2001, AJ, 122, 2267, doi: 10.1086/323717
  * Elias et al. (2020) Elias, L. M., Genel, S., Sternberg, A., et al. 2020, MNRAS, 494, 5439, doi: 10.1093/mnras/staa1059
  * Font-Ribera et al. (2013) Font-Ribera, A., Arnau, E., Miralda-Escudé, J., et al. 2013, J. Cosmology Astropart. Phys, 2013, 018, doi: 10.1088/1475-7516/2013/05/018
  * Fumagalli et al. (2011) Fumagalli, M., O’Meara, J. M., & Prochaska, J. X. 2011, Science, 334, 1245, doi: 10.1126/science.1213581
  * Gallego et al. (2018) Gallego, S. G., Cantalupo, S., Lilly, S., et al. 2018, MNRAS, 475, 3854, doi: 10.1093/mnras/sty037
  * Gallego et al. (2021) Gallego, S. G., Cantalupo, S., Sarpas, S., et al. 2021, MNRAS, doi: 10.1093/mnras/stab796
  * Giavalisco et al. (2011) Giavalisco, M., Vanzella, E., Salimbeni, S., et al. 2011, ApJ, 743, 95, doi: 10.1088/0004-637X/743/1/95
  * Górski et al. (2005) Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759, doi: 10.1086/427976
  * Gould & Weinberg (1996) Gould, A., & Weinberg, D. H. 1996, ApJ, 468, 462, doi: 10.1086/177707
  * Hamilton (1992) Hamilton, A. J. S. 1992, ApJ, 385, L5, doi: 10.1086/186264
  * Hawkins et al. (2003) Hawkins, E., Maddox, S., Cole, S., et al. 2003, MNRAS, 346, 78, doi: 10.1046/j.1365-2966.2003.07063.x
  * Heneka et al. (2017) Heneka, C., Cooray, A., & Feng, C. 2017, ApJ, 848, 52, doi: 10.3847/1538-4357/aa8eed
  * Hennawi et al. (2015) Hennawi, J. F., Prochaska, J. X., Cantalupo, S., & Arrigoni-Battaia, F. 2015, Science, 348, 779, doi: 10.1126/science.aaa5397
  * Ho et al. (2021) Ho, M.-F., Bird, S., & Garnett, R. 2021, MNRAS, 507, 704, doi: 10.1093/mnras/stab2169
  * Hu et al. (2019) Hu, W., Wang, J., Zheng, Z.-Y., et al. 2019, ApJ, 886, 90, doi: 10.3847/1538-4357/ab4cf4
  * Jose et al. (2013) Jose, C., Subramanian, K., Srianand, R., & Samui, S. 2013, MNRAS, 429, 2333, doi: 10.1093/mnras/sts503
  * Kaiser (1987) Kaiser, N. 1987, MNRAS, 227, 1, doi: 10.1093/mnras/227.1.1
  * Kennicutt (1998) Kennicutt, Robert C., J. 1998, ARA&A, 36, 189, doi: 10.1146/annurev.astro.36.1.189
  * Kollmeier et al. (2010) Kollmeier, J. A., Zheng, Z., Davé, R., et al. 2010, ApJ, 708, 1048, doi: 10.1088/0004-637X/708/2/1048
  * Kovetz et al. (2017) Kovetz, E. D., Viero, M. P., Lidz, A., et al. 2017, arXiv e-prints, arXiv:1709.09066.  https://arxiv.org/abs/1709.09066
  * Kusakabe et al. (2022) Kusakabe, H., Verhamme, A., Blaizot, J., et al. 2022, arXiv e-prints, arXiv:2201.07257.  https://arxiv.org/abs/2201.07257
  * Lake et al. (2015) Lake, E., Zheng, Z., Cen, R., et al. 2015, ApJ, 806, 46, doi: 10.1088/0004-637X/806/1/46
  * Leclercq et al. (2017) Leclercq, F., Bacon, R., Wisotzki, L., et al. 2017, A&A, 608, A8, doi: 10.1051/0004-6361/201731480
  * Li et al. (2021) Li, Z., Steidel, C. C., Gronke, M., & Chen, Y. 2021, MNRAS, 502, 2389, doi: 10.1093/mnras/staa3951
  * Lujan Niemeyer et al. (2022) Lujan Niemeyer, M., Komatsu, E., Byrohl, C., et al. 2022, ApJ, 929, 90, doi: 10.3847/1538-4357/ac5cb8
  * Lyke et al. (2020a) Lyke, B. W., Higley, A. N., McLane, J. N., et al. 2020a, ApJS, 250, 8, doi: 10.3847/1538-4365/aba623
  * Lyke et al. (2020b) —. 2020b, ApJS, 250, 8, doi: 10.3847/1538-4365/aba623
  * 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
  * Matsuda et al. (2004) Matsuda, Y., Yamada, T., Hayashino, T., et al. 2004, AJ, 128, 569, doi: 10.1086/422020
  * Matsuda et al. (2011) —. 2011, MNRAS, 410, L13, doi: 10.1111/j.1745-3933.2010.00969.x
  * Matsuda et al. (2012) —. 2012, MNRAS, 425, 878, doi: 10.1111/j.1365-2966.2012.21143.x
  * McCarthy et al. (2019) McCarthy, K. S., Zheng, Z., & Guo, H. 2019, MNRAS, 487, 2424, doi: 10.1093/mnras/stz1461
  * Mitchell et al. (2021) Mitchell, P. D., Blaizot, J., Cadiou, C., et al. 2021, MNRAS, 501, 5757, doi: 10.1093/mnras/stab035
  * Momose et al. (2014) Momose, R., Ouchi, M., Nakajima, K., et al. 2014, MNRAS, 442, 110, doi: 10.1093/mnras/stu825
  * Momose et al. (2016) —. 2016, MNRAS, 457, 2318, doi: 10.1093/mnras/stw021
  * Morrissey et al. (2018) Morrissey, P., Matuszewski, M., Martin, D. C., et al. 2018, ApJ, 864, 93, doi: 10.3847/1538-4357/aad597
  * Moster et al. (2013) Moster, B. P., Naab, T., & White, S. D. M. 2013, MNRAS, 428, 3121, doi: 10.1093/mnras/sts261
  * 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
  * Noterdaeme et al. (2012) Noterdaeme, P., Petitjean, P., Carithers, W. C., et al. 2012, A&A, 547, L1, doi: 10.1051/0004-6361/201220259
  * Ouchi et al. (2008) Ouchi, M., Shimasaku, K., Akiyama, M., et al. 2008, ApJS, 176, 301, doi: 10.1086/527673
  * Ouchi et al. (2010) Ouchi, M., Shimasaku, K., Furusawa, H., et al. 2010, ApJ, 723, 869, doi: 10.1088/0004-637X/723/1/869
  * Pâris et al. (2018) Pâris, I., Petitjean, P., Aubourg, É., et al. 2018, A&A, 613, A51, doi: 10.1051/0004-6361/201732445
  * Parsa et al. (2016) Parsa, S., Dunlop, J. S., McLure, R. J., & Mortlock, A. 2016, MNRAS, 456, 3194, doi: 10.1093/mnras/stv2857
  * Percival & White (2009) Percival, W. J., & White, M. 2009, MNRAS, 393, 297, doi: 10.1111/j.1365-2966.2008.14211.x
  * Planck Collaboration et al. (2020) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6, doi: 10.1051/0004-6361/201833910
  * Reddy & Steidel (2009) Reddy, N. A., & Steidel, C. C. 2009, ApJ, 692, 778, doi: 10.1088/0004-637X/692/1/778
  * Reid et al. (2015) Reid, B., Ho, S., Padmanabhan, N., et al. 2015, MNRAS, 455, 1553, doi: 10.1093/mnras/stv2382
  * Robertson et al. (2015) Robertson, B. E., Ellis, R. S., Furlanetto, S. R., & Dunlop, J. S. 2015, ApJ, 802, L19, doi: 10.1088/2041-8205/802/2/L19
  * Ross et al. (2013) Ross, N. P., McGreer, I. D., White, M., et al. 2013, ApJ, 773, 14, doi: 10.1088/0004-637X/773/1/14
  * Sawicki (2012) Sawicki, M. 2012, MNRAS, 421, 2187, doi: 10.1111/j.1365-2966.2012.20452.x
  * Seljak (2012) Seljak, U. 2012, J. Cosmology Astropart. Phys, 2012, 004, doi: 10.1088/1475-7516/2012/03/004
  * Shapley et al. (2003) Shapley, A. E., Steidel, C. C., Pettini, M., & Adelberger, K. L. 2003, ApJ, 588, 65, doi: 10.1086/373922
  * Silva et al. (2016) Silva, M. B., Kooistra, R., & Zaroubi, S. 2016, MNRAS, 462, 1961, doi: 10.1093/mnras/stw1777
  * Silva et al. (2013) Silva, M. B., Santos, M. G., Gong, Y., Cooray, A., & Bock, J. 2013, ApJ, 763, 132, doi: 10.1088/0004-637X/763/2/132
  * Sobral et al. (2018) Sobral, D., Santos, S., Matthee, J., et al. 2018, MNRAS, 476, 4725, doi: 10.1093/mnras/sty378
  * Sobral et al. (2013) Sobral, D., Smail, I., Best, P. N., et al. 2013, MNRAS, 428, 1128, doi: 10.1093/mnras/sts096
  * Steidel et al. (2000) Steidel, C. C., Adelberger, K. L., Shapley, A. E., et al. 2000, ApJ, 532, 170, doi: 10.1086/308568
  * Steidel et al. (2011) Steidel, C. C., Bogosavljević, M., Shapley, A. E., et al. 2011, ApJ, 736, 160, doi: 10.1088/0004-637X/736/2/160
  * Tramonte & Ma (2020) Tramonte, D., & Ma, Y.-Z. 2020, MNRAS, 498, 5916, doi: 10.1093/mnras/staa2727
  * Tramonte et al. (2019) Tramonte, D., Ma, Y.-Z., Li, Y.-C., & Staveley-Smith, L. 2019, MNRAS, 489, 385, doi: 10.1093/mnras/stz2146
  * Véron-Cetty & Véron (2010) Véron-Cetty, M. P., & Véron, P. 2010, A&A, 518, A10, doi: 10.1051/0004-6361/201014188
  * Wisotzki et al. (2016) Wisotzki, L., Bacon, R., Blaizot, J., et al. 2016, A&A, 587, A98, doi: 10.1051/0004-6361/201527384
  * Wisotzki et al. (2018) Wisotzki, L., Bacon, R., Brinchmann, J., et al. 2018, Nature, 562, 229, doi: 10.1038/s41586-018-0564-6
  * Witstok et al. (2019) Witstok, J., Puchwein, E., Kulkarni, G., Smit, R., & Haehnelt, M. G. 2019, arXiv e-prints, arXiv:1905.06954.  https://arxiv.org/abs/1905.06954
  * Wold et al. (2017) Wold, I. G. B., Finkelstein, S. L., Barger, A. J., Cowie, L. L., & Rosenwasser, B. 2017, ApJ, 848, 108, doi: 10.3847/1538-4357/aa8d6b
  * Xue et al. (2017) Xue, R., Lee, K.-S., Dey, A., et al. 2017, ApJ, 837, 172, doi: 10.3847/1538-4357/837/2/172
  * Zheng et al. (2011a) Zheng, Z., Cen, R., Trac, H., & Miralda-Escudé, J. 2011a, ApJ, 726, 38, doi: 10.1088/0004-637X/726/1/38
  * Zheng et al. (2011b) Zheng, Z., Cen, R., Weinberg, D., Trac, H., & Miralda-Escudé, J. 2011b, ApJ, 739, 62, doi: 10.1088/0004-637X/739/2/62