#  Photometric IGM tomography with Subaru/HSC: the large-scale structure of
Ly$\bm{\alpha}$ emitters and IGM transmission in the COSMOS field at
$\bm{z\sim 5}$

Koki Kakiichi,1 Joseph F. Hennawi,1,2 Yoshiaki Ono,3 Akio K. Inoue,4,5 Masami
Ouchi,6,3,7 Richard S. Ellis,8 Romain A. Meyer,9 and Sarah I. Bosman9
1Department of Physics, Broida Hall, University of California, Santa Barbara,
CA 93106-9530, USA
2Leiden Observatory, Leiden University, Niels Bohrweg 2, 2333 CA Leiden,
Netherlands
3Institute for Cosmic Ray Research, The University of Tokyo, 5-1-5
Kashiwanoha, Kashiwa, Chiba 277-8582, Japan
4Waseda Research Institute for Science and Engineering, Faculty of Science and
Engineering, Waseda University, 3-4-1, Okubo, Shinjuku, Tokyo 169-8555, Japan
5Department of Physics, School of Advanced Science and Engineering, Faculty of
Science and Engineering, Waseda University, 3-4-1, Okubo, Shinjuku, Tokyo
169-8555, Japan
6National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo
181-8588, Japan
7Kavli Institute for the Physics and Mathematics of the Universe (WPI),
University of Tokyo, Kashiwa, Chiba 277-8583, Japan
8Department of Physics and Astronomy, University College London, Gower Street,
London WC1E 6BT, UK
9Max Planck Institut für Astronomie, Königstuhl 17, D-69117, Heidelberg,
Germany E-mail<EMAIL_ADDRESS>(KK)

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We present a novel technique called “photometric IGM tomography” to map the
intergalactic medium (IGM) at $z\simeq 4.9$ in the COSMOS field. It utilizes
deep narrow-band (NB) imaging to photometrically detect faint Ly$\alpha$
forest transmission in background galaxies across the Subaru/Hyper-Suprime Cam
(HSC)’s $1.8\rm\,sq.\,deg$ field of view and locate Ly$\alpha$ emitters (LAEs)
in the same cosmic volume. Using ultra-deep HSC images and Bayesian spectral
energy distribution fitting, we measure the Ly$\alpha$ forest transmission at
$z\simeq 4.9$ along a large number ($140$) of background galaxies selected
from the DEIMOS10k spectroscopic catalogue at $4.98<z<5.89$ and the SILVERRUSH
LAEs at $z\simeq 5.7$. We photometrically measure the mean Ly$\alpha$ forest
transmission and achieve a result consistent with previous measurements based
on quasar spectra. We also measure the angular LAE-Ly$\alpha$ forest cross-
correlation and Ly$\alpha$ forest auto-correlation functions and place an
observational constraint on the large-scale fluctuations of the IGM around
LAEs at $z\simeq 4.9$. Finally, we present the reconstructed 2D tomographic
map of the IGM, co-spatial with the large-scale structure of LAEs, at a
transverse resolution of $11\,h^{-1}\rm cMpc$ across $140\,h^{-1}\rm cMpc$ in
the COSMOS field at $z\simeq 4.9$. We discuss the observational requirements
and the potential applications of this new technique for understanding the
sources of reionization, quasar radiative history, and galaxy-IGM correlations
across $z\sim 3-6$. Our results represent the first proof-of-concept of
photometric IGM tomography, offering a new route to examining early galaxy
evolution in the context of the large-scale cosmic web from the epoch of
reionization to cosmic noon.

###### keywords:

methods: observational – intergalactic medium – dark ages, reionization, first
stars – large-scale structure of Universe

††pubyear: 2015††pagerange:  Photometric IGM tomography with Subaru/HSC: the
large-scale structure of Ly$\bm{\alpha}$ emitters and IGM transmission in the
COSMOS field at $\bm{z\sim 5}$–Online supplementary material: Table of the
estimated Ly$\alpha$ forest transmission

##  1 Introduction

Cosmography, i.e. the science of mapping the Universe, is a fundamental pillar
of astronomy. Mapping the distribution of objects on the sky has been pivotal
for the discovery of the large-scale structure of the Universe. Our modern
cosmological model is largely based on the maps of cosmic microwave background
fluctuations (e.g Planck Collaboration et al., 2020), the large-scale
distribution of galaxies (e.g. eBOSS Collaboration: Alam et al., 2021), and
gravitational lensing (e.g DES Collaboration: Abbott et al., 2022). Mapping
the structure of the intergalactic medium (IGM) with 21-cm tomography has
great promise in advancing our understanding of the epoch of reionization and
cosmic dawn (e.g. Pritchard & Loeb, 2012). However, while significant progress
has been made in searching for the 21-cm power spectrum (Mertens et al., 2020;
Trott et al., 2020; Abdurashidova et al., 2022) and the global signal (Bowman
et al., 2018; Singh et al., 2022), there are still many challenges to overcome
before IGM tomography can be achieved with the 21-cm line.

Meanwhile, the Ly$\alpha$ forest remains the best probe of the IGM available
to date (e.g Becker et al. 2015a; McQuinn 2016 for reviews). Recent
measurements of effective optical depth indicate an ending of reionization as
late as $z\sim 5.3-5.7$ (Becker et al., 2015b; Eilers et al., 2018; Bosman et
al., 2022). The presence of long Gunn-Peterson troughs extending $\sim
10-100\,h^{-1}$ comoving Mpc (cMpc) (Becker et al., 2015b; Zhu et al., 2022)
and transmission spikes (Barnett et al., 2017; Yang et al., 2020) in the same
redshift range indicates large spatial variation in the IGM opacity at the
tail end of reionization. Simulations suggest the large-scale fluctuations of
Ly$\alpha$ forest transmission could be caused by the UV background
fluctuations from galaxies (Becker et al., 2015b; D’Aloisio et al., 2018;
Davies et al., 2018) or luminous rare sources such as AGN (Chardin et al.,
2015, 2017; Meiksin, 2020), thermal fluctuations in the IGM (D’Aloisio et al.,
2015; Keating et al., 2018), islands of neutral hydrogen due to the late end
of reionization (Kulkarni et al., 2019; Keating et al., 2020; Nasir &
D’Aloisio, 2020), and/or the spatially-varying distribution of self-shielded
absorbers which modulate the mean free path of the ionizing radiation (Davies
& Furlanetto, 2016; D’Aloisio et al., 2018). However, without directly
observing the sources of ionizing radiation, it is difficult to understand how
the interplay of these various physical processes and how they shape the
physical state of the IGM during the reionization process.

Establishing the direct spatial correlation between galaxy populations and
Ly$\alpha$ forest transmission of the IGM is key to understanding how
reionization proceeded. Previous ground-based surveys have mapped the
distribution of galaxies both photometrically (Becker et al., 2018; Kashino et
al., 2020; Christenson et al., 2021; Ishimoto et al., 2022) and
spectroscopically (Kakiichi et al., 2018; Meyer et al., 2019; Meyer et al.,
2020; Bosman et al., 2020) along sightlines to luminous $z\gtrsim 6$ quasars
where high quality Ly$\alpha$ forest spectra are available. By using
Subaru/Hyper-Suprime Cam (HSC) narrow-band (NB) imaging in a total of six
quasar fields, Becker et al. (2018); Christenson et al. (2021); Ishimoto et
al. (2022) find that $\sim 20\,h^{-1}\rm cMpc$ scale galaxy underdensities
(overdensities) at $z\simeq 5.7$ correlate with opaque (transmissive) regions
of the IGM on $\sim 50-100\,h^{-1}\rm cMpc$ scale. These observations favour
the scenario where ionizing radiation from galaxies drives the large-scale UV
background fluctuations at the tail end of reionization and/or completely
neutral islands still exist in the IGM at $z<5.7$. The spectroscopic survey
using Keck/DEIMOS and VLT/MUSE (Kakiichi et al., 2018; Meyer et al., 2020)
supports a similar picture. Surveying a total of eight quasar fields, Meyer et
al. (2020) find a large scale excess transmission of Ly$\alpha$ forest around
$z\simeq 5.8$ galaxies on scales of $\sim 10-40\,h^{-1}\rm cMpc$ at $\sim
2-3\sigma$. By modelling the observed galaxy-Ly$\alpha$ forest cross-
correlation signal, they interpreted that this excess transmission is caused
by the UV background fluctuations driven by the faint unseen population of
galaxies clustered around luminous galaxies, with an average Lyman continuum
(LyC) escape fraction of $\langle f_{\rm esc}\rangle\simeq 14\,\%$ at $z\simeq
5.8$.

Recent fully-coupled cosmological radiation hydrodynamic simulations also show
the excess transmission in the large-scale galaxy-Ly$\alpha$ forest cross-
correlation and suggest that the cross-correlation signal contains important
information about the timing of reionization (Garaldi et al., 2022). The
large-scale UV background fluctuations in the Ly$\alpha$ forest have long been
recognised to contain valuable information about the nature of LyC sources
(e.g. host halo mass) through their clustering properties (Pontzen, 2014;
Gontcho A Gontcho et al., 2014; Meiksin & McQuinn, 2019; Wolfson et al.,
2022). In summary, establishing the spatial correlation between galaxies and
Ly$\alpha$ forest at $5<z<7$ offers a smoking gun test of the reionization
process and is one of the key science goals of ongoing JWST quasar field
surveys (ID 2078, PI: Wang et al. 2021; ID 1243, PI: Lilly et al. 2017, see
also Kashino et al. 2022). However, due to the rarity of high-redshift
quasars, the connection between galaxies and the IGM in these surveys will
remain limited to one-dimensional skewers.

Ultimately we seek to map both galaxies and the IGM in three dimensions. At
intermediate redshifts $z\sim 2-3$, deep spectroscopic samples of background
galaxies have enabled the construction of 3D Ly$\alpha$ forest tomographic
maps of the IGM (Lee et al., 2014a, b, 2018; Newman et al., 2020; Horowitz et
al., 2021) and the measurements of the galaxy-Ly$\alpha$ forest cross-
correlation (Steidel et al., 2010; Chen et al., 2020). However, extending this
approach to higher redshifts is extremely challenging because the diminishing
Ly$\alpha$ forest transmission demands a much larger investment of telescope
time. Spectroscopically detecting the UV continua and Ly$\alpha$ forest
transmission of $z\simeq 5-6$ galaxies would require 30-m class telescopes
such as the Thirty-Meter Telescope (TMT), Giant Magellan Telescope (GMT), and
Extremely Large Telescope (ELT) (Japelj et al., 2019).

An alternative approach for IGM tomography is to utilize ultra-deep NB imaging
to detect the Ly$\alpha$ forest transmission at a fixed redshift against the
backdrop of more distant galaxies. As the throughput of an imager is much
higher than a typical spectrograph, this photometric measurement of the
Ly$\alpha$ forest transmission can be more sensitive than one based on
spectroscopy. For the $\sim 60\,\%$ throughput of the HSC imager (cf. $\sim
10-20\,\%$ throughput of a typical spectrograph), the NB measurement of the
Ly$\alpha$ forest with the Subaru 8.2-m can be comparable to undertaking a
spectroscopic IGM survey with a 14-20 m telescope. The typical NB filter width
is $\simeq 100\,\mbox{\AA}$ which corresponds to a line-of-sight distance of
$\sim 30\,h^{-1}\rm cMpc$ at $z\sim 5-6$. This matches the scale of
fluctuations in the Ly$\alpha$ forest transmission seen by previous quasar
field surveys (Becker et al., 2018; Kakiichi et al., 2018; Kashino et al.,
2019; Meyer et al., 2019; Meyer et al., 2020; Christenson et al., 2021;
Ishimoto et al., 2022). Since up to a few hundred background galaxies can be
identified using extant spectroscopic catalogues and NB-selected Ly$\alpha$
emitters (LAEs) in well-studied extragalactic fields, photometric IGM
tomography can provide a $\times$ 100 increase in the number of galaxy-
Ly$\alpha$ sightline pairs compared to quasar surveys - a huge boost in
statistical power. In an earlier article, Kakiichi et al. (2022) outlined the
strategy for photometric IGM tomography which provides a path forward to map
to map the Ly$\alpha$ forest transmission of the IGM and Ly$\alpha$ emitting
galaxies in the same cosmic volume using only imaging data.

In this paper, we apply this photometric IGM tomography technique to the well-
studied extragalactic COSMOS field and present measurements of the LAE-
Ly$\alpha$ forest cross-correlation and the auto-correlation of the Ly$\alpha$
forest at $z\simeq 4.9$. The wealth of deep multi-wavelength imaging and
spectroscopic data makes COSMOS an ideal field to demonstrate the method. We
provide the first large-scale 2D tomographic map of the IGM at $z\simeq 4.9$
across the $1.8\rm\,deg^{2}$ ($\sim 140\,h^{-1}\rm cMpc$ in diameter) field of
view of Subaru/HSC, enabling us to directly visualise the spatial connection
between galaxies and the IGM. After the reionization process is complete, we
expect that the galaxy-Ly$\alpha$ forest cross-correlation will evolve from
positive (Kakiichi et al., 2018; Meyer et al., 2020; Garaldi et al., 2019)
owing to the large-scale UV background fluctuations and/or ionized bubbles to
a negative (i.e. anti-correlation) signal (e.g. Turner et al., 2017; Nagamine
et al., 2021) due to the increasing impact of gas overdensities around
galaxies at lower redshifts. Our study at $z\simeq 4.9$ will provide a clue
for how the cross-correlation evolves across cosmic time.

In Section 2 we describe the imaging data used in this analysis. Section 3
describes the galaxy catalogues and the selection for the photometric IGM
tomography. Section 4 presents the method to estimate the Ly$\alpha$ forest
transmission along background galaxies using a Bayesian spectral energy
distribution (SED) fitting framework. In Sections 5-8, we present our main
results including the measurements of mean Ly$\alpha$ forest transmission
(Section 5), LAE-Ly$\alpha$ forest cross-correlation (Section 6), auto-
correlation of the Ly$\alpha$ forest (Section 7), and the reconstruction of
the 2D IGM tomographic map (Section 8). In Section 9 we discuss the
requirement to improve the photometric IGM tomography and various science
applications. We summarise our results and conclusions in Section 10.
Throughput this paper we assume cosmological parameters
$(\Omega_{m},\Omega_{\Lambda},\Omega_{b},h,\sigma_{8},n_{s})=(0.3089,0.6911,0.0486,0.6774,0.8159,0.9667)$
(Planck Collaboration et al., 2016). We use cMpc (pMpc) to indicate distances
in comoving (proper) units. All magnitudes in this paper are quoted in the AB
system (Oke & Gunn, 1983).

##  2 Data

Figure 1: $5\sigma$ limiting magnitudes of the NB718, $z$, and $y$-band images with fixed $1.5^{\prime\prime}$ aperture in the COSMOS field (tract 9813). Table 1: Summary of the foreground NB images Filter | Ly$\alpha$ redshift | bkg. source redshift |  $5\sigma$ depth (1.5″) |  $3\sigma$ depth (1.5″) |  $1\sigma$ depth (1.5″) | Ref
---|---|---|---|---|---|---
|  |  | [AB mag] | [AB mag] | [AB mag] |
NB718 | 4.90 $(4.85\,<\,z_{\rm Ly\alpha}\,<\,4.94)$ | $4.98<z<5.89$ | 26.30 | 26.86 | 28.05 | Inoue et al (2020)
b computed from the median of random sky objects in masked region (i.e.
excluding near bright stars and artefacts) in each patch of tract 9813.
Table 2: Summary of the BB images from HSC-SSP DR3 in the UD-COSMOS field (tract 9813) Median $5\sigma$ depth (1.5″) [mag] | Reference
---|---
$g$ | $r$ | $i$ | $z$ | $y$ |
27.85 | 27.39 | 27.22 | 26.86 | 26.23 | Aihara et al (2022)

We use the public release of a Subaru HSC NB718 image from CHORUS DR1 (Inoue
et al., 2020) and broad-band (BB) $grizy$ and NB816 images from HSC-SSP DR3
(Aihara et al., 2022) in the ultra-deep layer of COSMOS field (tract 9813).
These images cover an area of approximately $1.67\times 1.67\rm\,deg^{2}$
centred at $\rm(RA,DEC)=(10h01m00s,+2d14m00s)$. Co-added images are retrieved
from the public data release website.111HSC-SSP DR3: https://hsc-
release.mtk.nao.ac.jp/doc/index.php/data-access__pdr3/,222CHORUS DR1:
https://hsc-release.mtk.nao.ac.jp/doc/index.php/chorus/ The ancillary data
including PSF FWHM and limiting magnitudes for each patch is also downloaded
from here. We use NB718 as a foreground NB filter to measure the Ly$\alpha$
forest transmission at $z\simeq 4.9$ (the central wavelength of
$7170.5\,\mbox{\AA}$ corresponds to Ly$\alpha$ redshift of $z_{\rm
Ly\alpha}=4.898$) covering the range of $4.85<z<4.94$ and $z$\- and $y$-bands
to measure the UV continua of the background galaxies.

We mask regions contaminated by artefacts (bright star haloes, ghosts,
blooming, channel-stop, dip). Since we measure the residual transmitted fluxes
in the NB718 image, particular care is needed for this filter as artefacts
could potentially cause false positive detections. To address this we first
flag all the masked pixels reported by CHORUS PDR1 (Inoue et al., 2020). This
includes the pixels with the hscPipe flags: pixelflags_bright_object=True
(pixels affected by bright objects) or pixelflags_saturatedcenter=True (pixels
affected by count saturation). Pixels affected by the haloes of bright stars
are also masked. After visual inspection of the NB718 image with the CHORUS
PDR1 mask overlaid, we find that there are still some regions affected by the
outer ghosts of bright stars, which extend to approximately $\sim
270\rm\,arcsec$ in radius. To mask these, we conservatively follow the
procedure adopted by HSC-SSP DR3 (Aihara et al., 2022). We select bright stars
with $G<10\rm\,mag$ from GAIA DR3
catalogue333https://gea.esac.esa.int/archive/ in the footprint of tract 9813.
We then mask regions inside 320 arcsec radius around these Gaia stars. 320
arcsec corresponds to the size of the outermost ghost identified by Aihara et
al. (2022). They find that the outer ghost is significant for a star brighter
than $\sim 7\rm\,mag$ and the inner ghost is dominant at $\sim 7-9\rm\,mag$.
We apply the same mask for the broad-band images.

The limiting magnitudes of the NB and BB images are estimated using synthetic
apertures randomly distributed in the blank sky regions of the image. For
NB718 we use the published limiting magnitude map from the CHORUS PDR1. For
the BB images, we retrieve the synthetic apertures located in the empty
regions of the sky (hereafter sky objects) from HSC database and perform
photometry with a fixed $1.5\arcsec$ aperture. As the sensitivity varies
across the field of view, we estimate the limiting magnitudes for each patch.
In each patch, there are typically $\sim 40$ sky objects and we compute the
limiting magnitude from the standard deviation of the photometric measurements
of the fluxes for the sky objects. Figure 1 shows the limiting magnitudes of
NB718, $z$, and $y$-bands for each patch in the field. The NB limiting
magnitudes vary by about $\sim 0.18$ mag. The median depths of the NB and BB
images are summarised in Tables 1 and 2.

A rule-of-thumb for the required depth for IGM tomography is given by Kakiichi
et al. (2022). Assuming the flat UV continuum slope of a background galaxy,
the NB magnitude needs to reach

$m_{\rm NB}=m_{\rm BB}-2.5\log_{10}e^{-\tau_{\rm eff}(z)}\approx m_{\rm
BB}+\tau_{\rm eff}(z),$ (1)

to detect the Ly$\alpha$ forest transmission with an effective optical depth
$\tau_{\rm eff}(z)$. Here, the NB magnitude corresponds to NB718 and the BB
magnitude corresponds to $z$-band which covers the UV continuum of a
background galaxy. Assuming the effective optical depth of $\tau_{\rm
eff}(z)=1.5$ at $z=4.9$ (Becker et al., 2013; Eilers et al., 2018; Bosman et
al., 2022), the required NB718 depth for a $z=25.0\rm\,mag$ background source
is $26.5\,\rm mag$. The existing NB718 depth meets this requirement at
$>3\sigma$ (Table 1). For a fainter background source with $z=26.5\rm\,mag$,
the existing NB718 depth still has a $\sim 1\sigma$ sensitivity to the mean
Ly$\alpha$ forest transmission. The HSC imaging of the COSMOS field thus has
sufficient sensitivity for photometric IGM tomography.

###  2.1 Photometry

We measure the $grizy$ and NB718 photometry from HSC-SSP DR3 and CHORUS PDR1
data using a fixed $1.5\arcsec$ aperture for the background sources. The zero
points of all photometric bands is $27\,\rm mag/DN$ according to the HSC-SSP
data release. We ignore a few percent level correction arising from aperture
corrections during the photometric calibration stage. This is negligible
compared with the other photometric errors described below. We assign the
photometric error of an object based on the limiting magnitude of the patch
where the object is located.

##  3 Catalogues

The redshift range for the background sources is chosen such that the
Ly$\alpha$ forest range between Ly$\beta$ ($\lambda_{\beta}=1026$ Å) and
Ly$\alpha$ ($\lambda_{\alpha}=1216$ Å) lines is covered by the NB718 filter.
The lower and upper redshifts are set by $z_{\rm min}=\lambda_{\rm
NB,max}/\lambda_{\alpha}-1$ and $z_{\rm max}=\lambda_{\rm
NB,min}/\lambda_{\beta}-1$ where $\lambda_{\rm NB,max}$ $\lambda_{\rm NB,min}$
are the maximum and minimum wavelengths of the filter. We define the minimum
and maximum wavelengths as a range where the NB718 filter transmission is
$>50\,\%$. The appropriate background source redshift for $z=4.9$ IGM
tomography is thus $4.98<z<5.89$.

To locate foreground galaxies at $z=4.9$ in the same redshift slice
corresponding to that for which our IGM transmission is being measured, we
employ the SILVERRUSH catalogue of $z\simeq 4.9$ LAEs (ver20210224, Ono et al.
2021).

To identify background galaxies at $4.98<z<5.89$, we use both the SILVERRUSH
catalogue of $z\simeq 5.7$ LAEs (Ono et al., 2021) and the spectroscopic
redshift (spec-z) catalogue compiled with HSC-SSP DR3 (Aihara et al., 2022).
For the latter, we find that DEIMOS10k (Hasinger et al., 2018) was the primary
source of the background galaxies as we will describe below. Thus in the
remainder of the paper, we refer the background galaxies derived from the
catalogues to as the LAE and DEIMOS10k samples, respectively.

###  3.1 Spectroscopic redshift catalogue

The spec-z catalogue associated with HSC-SSP DR3 is a compilation of public
spectroscopic redshifts from numerous previous redshift surveys including
2dFGRS (Colless et al., 2003), 3D-HST (Skelton et al., 2014; Momcheva et al.,
2016), 6dFGRS (Jones et al., 2004, 2009), C3R2 DR2 (Masters et al., 2017;
Masters et al., 2019), DEEP2 DR4 (Davis et al., 2003; Newman et al., 2013),
DEEP3 (Cooper et al., 2011, 2012), DEIMOS10k (Hasinger et al., 2018), FMOS-
COSMOS (Silverman et al., 2015; Kashino et al., 2019), GAMA DR2 (Liske et al.,
2015), LEGA-C DR2 (Straatman et al., 2018), PRIMUS DR1 (Coil et al., 2011;
Cool et al., 2013), SDSS DR16 (Ahumada et al., 2020), SDSS IV QSO catalog
(Pâris et al., 2018), UDSz (Bradshaw et al., 2013; McLure et al., 2013),
VANDELS DR1 (Pentericci et al., 2018), VIPERS PDR1 (Garilli et al., 2014),
VVDS (Le Fèvre et al., 2013), WiggleZ DR1 (Drinkwater et al., 2010), and
zCOSMOS DR3 (Lilly et al., 2009).

Spectroscopic sources are matched to the HSC photometry by position, thus the
catalogue only includes the objects detected by HSC-SSP DR3. Using the CAS
search, we retrieve the HSC-SSP DR3 spec-z catalogue after it was cross-
matched using the object_id. We re-measure the fixed $1.5^{\prime\prime}$
aperature photometry in $grizy$ and NB718 for all the objects to ensure
consistent measurements of the fluxes across all bands.

Figure 2: An example background galaxy from the DEIMOS10k catalogue. (Left):
Postage stamps show $5\times 5$ arcsec cutouts of $grizY$ images around the
object marked with a white crosshairs. (Right): DEIMOS 1D spectrum (flux:
black, red: noise). The vertical dotted line indicates the Ly$\alpha$ line.
Figure 3: An example background LAE from the SILVERRUSH catalogue. (Left):
Postage stamp showing $5\times 5$ arcsec cutouts of $grizY$ images with the
object marked with a white cross. (Right): NB816 image. The NB colour excess
is clearly detected.

Using the spec-z catalogue, we search for background source candidates at
$4.98<z<5.89$ in the UD-COSMOS field (tract 9813). The catalogue contains a
spec-z flag (specz_flag_homogeneous=True for secure and False for insecure)
after homogenising the quality flags444https://hsc-
release.mtk.nao.ac.jp/doc/index.php/catalog-of-spectroscopic-redshifts__pdr3/
of spectroscopic redshifts from the above surveys. Selecting only objects with
specz_flag_homogeneous=True, we find 236 candidates in the required redshift
range, of which 60 belong to DEIMOS10k (Hasinger et al., 2018) and 176 belong
to 3D-HST (v4.1.5, Momcheva et al. 2016). The majority of bright candidates
with $z\lesssim 25.5$ comes from DEIMOS10k whereas fainter candidates are
mostly from 3D-HST.

In order to visually confirm the spectroscopic redshifts, we downloaded the
original DEIMOS10k spectra from NASA/IPAC Infrared Science Archive
(IRSA)555https://irsa.ipac.caltech.edu/data/COSMOS/overview.html, and the
3D-HST spectra from MAST archive666https://archive.stsci.edu/prepds/3d-hst/.
For 54 of the 60 DEIMOS10k sources, we confirmed an emission line feature
(mostly single Ly$\alpha$ line, one Lyman break only, one quasar). We rejected
6 objects because either (1) no published spectrum is available (DEIMOS10k ID:
L234173) or (2) we could not visually confirm the reported redshift (L420065,
L430951, L378903, C563716, L442206).

Out of the remaining 54 DEIMOS10k sources, we removed 11 residing in the
masked regions. Furthermore, in order to secure a reliable UV continuum
detection in each background sources, we applied a $5\sigma$ detection cut in
the $z$-band using the limiting magnitude appropriate for the relevant patch.
3 candidates fail to meet this criterion in the $z$-band of HSC-SSP DR3. We
also require a $3\sigma$ non-detection in $g$-band in order to reject low-
redshift interlopers. This leads to the removal of a further 4 sources. As a
result we finally have 36 DEIMOS10k background sources for our IGM tomography.
We present the postage stamp image and DEIMOS spectrum of a representative
background source in Figure 2.

For the 3D-HST sources, the catalogued redshifts are determined from either
photometric and/or grism spectroscopic data. The ACS/G800L grism spectra cover
Ly$\alpha$ in our desired redshift range. We downloaded the 3D-HST catalogue
and find that all the relevant sources have only photometric redshifts; we
could not find any sources with a convincing Ly$\alpha$ line or Lyman break in
the grism spectra. In principle, we can use background objects with photo-z’s
whose 95% confidence interval (i.e. z_best_l95, z_best_u95 from
cosmos_3dhst.v4.1.5.zbest.dat) lie within our required range. This would
ensure that their Ly$\alpha$ forest region is appropriately covered by the
NB718 filter. However, since it is unclear how catastrophic photo-z errors
might affect the quality of our IGM tomography, we decided to remove all the
3D-HST objects from our final background sources in this paper. Nonetheless,
in future work, it will be interesting to examine the utility of the photo-z
background sources for IGM tomography.

To summarise, our final catalogue of background spec-z sources contains 36
objects from DEIMOS10k. Their $z$-band magnitudes and distribution in the UD-
COSMOS field are shown in Figures 4 and 5, respectively. Because the previous
DEIMOS spectroscopic campaigns are focused near the central region of the
COSMOS field, their distribution reflects their survey footprints.

Figure 4: Magnitude distribution of background sources for NB718 IGM
tomography at $z\simeq 4.9$. Our final background source catalogue contains
151 sources in total (red: 36 spec-z objects from DEIMOS10k, blue: 115 LAE
objects from SILVERRUSH). Figure 5: Sky distribution of background sources for
NB718 IGM tomography (blue circles: $z=5.7$ LAEs from SILVERRUSH, red squares:
spec-z from DEIMOS 10k). Masked regions and those outside the field-of-view
are indicated by the gray shaded regions.

###  3.2 LAE catalogue

We use the LAE catalogue from Ono et al. (2021) constructed as part of the
SILVERRUSH programme. The catalogue is based on the data from CHORUS survey
(Inoue et al., 2020) and from the HSC-SSP internal data release of S18A which
is basically identical to the Public Data Release 2 (Aihara et al., 2019). To
ensure homogeneous photometric measurements for the IGM tomography, we re-
measure the fixed $1.5^{\prime\prime}$ aperture photometry at the coordinates
of the SILVERRUSH LAEs using HSC-SSP DR3 $grizy$ and CHORUS NB718 images. We
use the SILVERRUSH catalogue to select both background and foreground LAEs at
$z=5.7$ and $z=4.9$ located by NB816 and NB718 colour excess, respectively.

####  3.2.1 Background LAE selection: $z\simeq 5.7$

To select $z\simeq 5.7$ background LAEs, we draw a sample from the SILVERRUSH
catalogue (Ono et al., 2021) which applies a NB colour excess $i-{\rm
NB816}\geq 1.2$ and $5\sigma$ detection in NB816. Such a NB816 colour excess
can locate LAEs at a redshift $z=5.726$ with an $\Delta z\simeq 0.1$ accuracy
(Ono et al., 2021) sufficient to ensure that the Ly$\alpha$ forest
transmission is covered by the NB718 filter. Details of the LAE catalogue
construction are described in Ono et al. (2021). The SILVERRUSH catalogue
contains 378 $z\simeq 5.7$ LAEs in the UD-COSMOS field. This double NB
technique (Kakiichi et al., 2022) allows us to efficiently assemble a large
number of background sources for IGM tomography.

In addition to the standard NB selection, we require a $5\sigma$ detection in
$z$-band. This criterion is met by 176 objects out of 378 LAEs. We also
removed 35 further objects which reside in the masked regions. While Ono et
al. (2021) already applied masks in constructing the original catalogue, our
revised mask in NB718 is more conservative. As before, we also require a
$3\sigma$ non-detection in the $g$-band to avoid possible low-redshift
interlopers; this removes a further 26 objects.

Thus, our final background LAE catalogue contains 115 objects. We visually
inspected all of these objects in HSC DR3 $grizy$ and NB816 images. The
$z$-band magnitudes and the spatial distribution of the background LAEs are
shown in Figures 4 and 5. The background LAEs are typically fainter than the
DEIMOS10k sample, but distributed more evenly across the entire UD-COSMOS
field as they are selected homogeneously via NB816 colour excess. We present
an example postage stamp of a background LAE in Figure 3.

####  3.2.2 Foreground LAE selection: $z\simeq 4.9$

In order to cross-correlate foreground LAEs with the Ly$\alpha$ forest
transmission, we also use the HSC NB718 data to select LAEs at $z\simeq 4.9$.
The SILVERRUSH catalogue applies the selection criteria: $ri-{\rm NB718}>0.7$
and $r-i>0.8$ and $ri-{\rm NB718}>(ri-{\rm NB718})_{3\sigma}$ and
$g>g_{2\sigma}$ where $ri$ is calculated by the linear combination of the
fluxes in $r$\- and $i$-bands, $f_{r}$ and $f_{i}$, following
$f_{ri}=0.3f_{r}+0.7f_{i}$ and the $2\sigma$ and $3\sigma$ subscripts denote
$2\sigma$ and $3\sigma$ limiting magnitudes (Ono et al., 2021). Further detail
is described in Ono et al. (2021). This gives 280 $z\simeq 4.9$ LAEs in the
UD-COSMOS field. We remove 17 objects lying in our updated masked regions. For
these foreground LAEs, unlike the background galaxies, we do not apply any
$z$-band detection cut and use all 263 NB718-selected LAEs for our subsequent
analysis. The average luminosities of the foreground LAEs are summarised in
Table 3.

Table 3: Average physical properties of foreground LAEs. Redshift |  $\log_{10}\langle L_{\alpha}\rangle$ (2.0″) |  $\langle M_{\rm UV}\rangle$ (2.0″)
---|---|---
| [$\rm erg\,s^{-1}$] | [AB mag]
$z=4.89$ | 42.62 | $-20.09$

##  4 IGM Ly$\alpha$ forest transmission

A measurement of the IGM Ly$\alpha$ forest transmission $T_{\rm IGM}$ using
the foreground NB718 filter requires us to infer the intrinsic spectral energy
distribution (SED) of each background galaxy in the absence of any IGM
absorption. The NB-integrated Ly$\alpha$ forest transmission is then
determined by the ratio between the observed and intrinsic NB fluxes,

$T_{\rm IGM}=\frac{f^{\rm obs}_{\rm NB}}{f^{\rm intr}_{\rm NB}}.$ (2)

There are several ways to perform this measurement. One obvious way to
estimate $f_{\rm NB}^{\rm intr}$, similar to the approach employed by Mawatari
et al. (2017), is to first fit a SED to each background galaxy using broad-
band photometry redward of the Ly$\alpha$ emission line $>1216$ Å , and then
extrapolate the continuum to the relevant rest-frame range of the Ly$\alpha$
forest between $1026$ Å and $1216$ Å covered by the foreground NB718 filter.
While intuitive, it is difficult to rigorously propagate the photometric
errors and systematic uncertainties of the SED modelling into the final
measurement of $T_{\rm IGM}$. Also, it is hard to quantify likely degeneracies
between the SED parameters (e.g. UV continuum slope, or age and dust
attenuation law) and the Ly$\alpha$ forest transmission.

A better way is to simultaneously fit both the galaxy SED and the Ly$\alpha$
forest transmission of the IGM in a fully Bayesian framework. This allows a
rigorous propagation of photometric and systematic errors in the $T_{\rm IGM}$
estimate for each background galaxy and characterises the full posteriors
including the degeneracy with the SED model parameters.

Figure 6: Representative examples illustrating results from the Bayesian SED
fitting framework for the DEIMOS10k (top panels) and LAE (bottom panels)
samples. The examples show a case for the detection of the transmitted
Ly$\alpha$ forest flux in NB718 (e.g. DEIMOS_2018_519281) and one for a non-
detection (e.g. SILVERRUSH_2021_7140). (Left): The best-fit power-law SEDs
(black solid) with the $14-86\%$ and $5-95\%$ confidence intervals (dark and
light grey shaded regions) are overlaid on the measured NB718 (blue), $z$
(yellow), and $y$ (red) band fluxes of the background source. The wavelength
coverage of each NB718, $z$ and $y$-band filter is indicated by the
transparent filled curves. $5^{\prime\prime}\times 5^{\prime\prime}$ postage
stamps show the images of $g$, NB718, $z$, and $y$ fluxes smoothed by a
Gaussian kernel with a standard deviation of one pixel. The colourbar scales
are the same for all images. (Right): MCMC corner plots of the key parameters
($T_{\rm IGM},M_{\mbox{\tiny UV}},\beta$).

###  4.1 Bayesian SED fitting framework

We apply a Bayesian SED fitting framework to measure $T_{\rm IGM}$. We forward
model the observed photometric fluxes in narrow- and broad-band filters using
realistic HSC filter transmission curves $t_{\rm NB}(\nu)$ and $t_{\rm
BB}(\nu)$ including the CCD quantum efficiency, the transmittance of the dewar
window and the Primary Focus Unit of the HSC. For our IGM tomography, $\rm
NB=NB718$ and ${\rm BB}=z,y$ since we use NB718 to measure the Ly$\alpha$
forest transmission and $z$\- and $y$-bands to constrain the intrinsic galaxy
SED.

We denote a model SED of a background galaxy by $L_{\nu}(\nu_{\rm
e}|\bm{\Theta})$ (in unit of $\rm erg\,s^{-1}\,Hz^{-1}$) where $\nu_{e}$ is
the emitted frequency at the rest-frame of the galaxy and $\bm{\Theta}$ is a
set of the SED parameters. We assume a power-law SED with
$L_{\nu}=L_{\nu}(1500\mbox{\AA})(\nu_{e}/\nu_{1500})^{-(2+\beta)}$ where
$L_{\nu}(1500\mbox{\AA})$ and $\nu_{1500}$ represent the luminosity and
frequency at $1500\,\mbox{\AA}$ respectively, and $\beta$ is the UV continuum
slope. Thus our SED parameters are $\bm{\Theta}=\\{M_{\rm UV},\beta\\}$. At
the known redshift of the object (either from spectroscopy or NB detection of
Ly$\alpha$), the observed flux is
$f_{\nu}(\nu|\bm{\Theta})=(1+z)L_{\nu}[\nu_{e}=\nu(1+z)|\bm{\Theta}]/(4\pi
D_{L}(z)^{2})$ where $\nu$ is the observed frequency and $D_{L}(z)$ is the
luminosity distance.

As the foreground NB718 filter covers a portion of Ly$\alpha$ forest of a
background galaxy, the observed NB718 flux is attenuated by
$e^{-\tau_{\alpha}}$ where $\tau_{\alpha}$ is the the Ly$\alpha$ optical depth
of the IGM, therefore,

$f_{\rm NB}(T_{\rm IGM},\bm{\Theta})=\frac{\int
e^{-\tau_{\alpha}}f_{\nu}(\nu|\bm{\Theta})\,t_{\rm NB}(\nu)d\nu}{\int t_{\rm
NB}(\nu)d\nu}\approx T_{\rm IGM}f_{\rm NB}^{\rm intr}(\bm{\Theta}).$ (3)

We define the NB-integrated Ly$\alpha$ forest transmission of the IGM as
$T_{\rm IGM}=\left.\int e^{-\tau_{\alpha}}\,t_{\rm NB}(\nu)d\nu\right/\int
t_{\rm NB}(\nu)d\nu$. In the absence of the IGM, the NB718 flux is $f_{\rm
NB}^{\rm intr}(\bm{\Theta})=\left.\int f_{\nu}(\nu|\bm{\Theta})\,t_{\rm
NB}(\nu)d\nu\right/\int t_{\rm NB}(\nu)d\nu$. The BB fluxes redward of
Ly$\alpha$ are not affected by the IGM. They are thus modelled as

$f_{\rm BB}(\bm{\Theta})=\frac{\int f_{\nu}(\nu|\bm{\Theta})\,t_{\rm
BB}(\nu)d\nu}{\int t_{\rm BB}(\nu)d\nu}.$ (4)

We assume the observed photometric noise follows a Gaussian distribution and
that noise levels in the various filters do not correlate with one another.
Therefore, the likelihood can be written as the sum of Gaussian likelihoods,

$\ln\mathcal{L}=-\frac{1}{2}\left[\frac{f_{\rm NB}^{\rm obs}-f_{\rm NB}(T_{\rm
IGM},\bm{\Theta})}{\sigma_{\rm NB}}\right]^{2}-\frac{1}{2}\sum_{{\rm
BB}=z,y}\left[\frac{f_{\rm BB}^{\rm obs}-f_{\rm BB}(\bm{\Theta})}{\sigma_{\rm
BB}}\right]^{2}.$ (5)

Using Bayes theorem, we can express the posterior as a product of prior
$P(T_{\rm IGM},\bm{\Theta})$ and likelihood $\mathcal{L}$,

$P(T_{\rm IGM},\bm{\Theta}|f_{\rm NB}^{\rm obs},\bm{f}_{\rm BB}^{\rm
obs})\propto P(T_{\rm IGM},\bm{\Theta})\mathcal{L}(f_{\rm NB}^{\rm
obs},\bm{f}_{\rm BB}^{\rm obs}|T_{\rm IGM},\bm{\Theta}).$ (6)

Thus the measurement of the Ly$\alpha$ forest transmission $T_{\rm IGM}$ along
each background galaxy is given by the marginalized posterior over the SED
parameters $\bm{\Theta}$,

$P(T_{\rm IGM}|f_{\rm NB}^{\rm obs},\bm{f}_{\rm BB}^{\rm obs})=\int P(T_{\rm
IGM},\bm{\Theta}|f_{\rm NB}^{\rm obs},\bm{f}_{\rm BB}^{\rm obs})d\bm{\Theta}.$
(7)

We implement this Bayesian SED fitting framework using a Markov Chain Monte
Carlo method emcee (Foreman-Mackey et al., 2013). We use flat priors in the
range of $-100<T_{\rm IGM}<100$, $-30<M_{\mbox{\tiny UV}}<-15$, and
$-5<\beta<2$ as our default. We justify the very wide range of flat priors
(instead of imposing a physical range between 0 and 1 for individual
measurements for $T_{\rm IGM}$) in Section 5.

When fitting the observed SEDs, we find occasional cases where the
marginalized posterior peaks at an unphysically large value $T_{\rm IGM}>1$.
This indicates there may be an unaccounted systematic error in our data. In
order to assess and remove such objects, we compute the probability that the
estimated $T_{\rm IGM}$ is greater than 1 using the posterior, i.e. $P(T_{\rm
IGM}>1|f_{\rm NB}^{\rm obs},\bm{f}_{\rm BB}^{\rm
obs})=\int_{1}^{\infty}P(T_{\rm IGM}|f_{\rm NB}^{\rm obs},\bm{f}_{\rm BB}^{\rm
obs})dT_{\rm IGM}$. We then flag objects with $P(T_{\rm IGM}>1|f_{\rm NB}^{\rm
obs},\bm{f}_{\rm BB}^{\rm obs})>50\,\%$. This choice is motivated by the fact
that a Gaussian posterior centred at $T_{\rm IGM}>1$ gives $>50\,\%$
probability that $T_{\rm IGM}$ is greater than 1, suggesting an unmodelled
systematic error while SED fitting. We then visually check the original images
and confirm such cases originate from systematic errors such as the under-
subtraction of the sky background in NB718 and/or contamination from nearby
objects in the photometric aperture. We apply this validation procedure in
both the DEIMOS10k and LAE samples. We find no such cases (out of 36) in the
DEIMOS10k catalogue. However, 11 cases (out of 115) are flagged in the LAE
sample. Here, our visual inspection confirms likely contamination from nearby
objects or a diffuse NB718 image compared to that in the BB $z$-band,
suggesting that the region may be affected by the faint halo or ghost of a
nearby bright star or by incorrect sky subtraction. Thus these 11 objects are
removed and we use the remaining 104 background LAEs for our subsequent
analysis.

###  4.2 Individual IGM Ly$\alpha$ forest transmission measurements

Figure 7: (Top): Distribution of estimated $T_{\rm IGM}$ values from the mean
of each posterior of individual background galaxies (blue: LAEs, red:
DEIMOS10k) using a flat prior (filled histogram) and Gaussian prior (step
histogram) on the $\beta$ UV slope. The median standard deviation of the
posteriors is indicated as the typical error of the measurement. (Bottom): As
above but for the UV slope $\beta$. The Gaussian fit to the distribution of
$\beta$ slopes from Bouwens et al. (2014) is indicated by the dashed line.

Figure 8: Comparison of the intrinsic galaxy SED without the IGM absorption
from BPASS+CLOUDY and power-law fits to the LAE (left) and DEIMOS10k (right)
samples. The dark (light) gray regions indicate the $16-84$ ($5-95$)
percentiles of the range of the best-fit power-law models with the white line
indicating the population-averaged $\beta$ slope. The default parameters of
the BPASS+CLOUDY model are stellar age $t_{\rm age}=10\rm\,Myr$, metallicity
$Z=0.20Z_{\odot}$, ionization parameter $\log_{10}U=-2.5$ and Calzetti (2001)
dust attenuation $E(B-V)=0.1$. The filter transmission curves for NB718
(blue), $z$ (yellow), and $y$ (red) bands are indicated by the shaded regions.

By applying the Bayesian SED fitting framework, we measure the Ly$\alpha$
forest transmission of the IGM along sightlines to individual background
galaxies. The method simultaneously returns the constraints on the Ly$\alpha$
forest transmission $T_{\rm IGM}$ and the SED parameters ($M_{\mbox{\tiny
UV}}$ and $\beta$) for each source. In Figure 6 we show two representative
examples of the derived constraints on these physical parameters for the
DEIMOS10k and LAE samples. When the transmitted Ly$\alpha$ forest flux is
detected in NB718 as demonstrated by the DEIMOS10k sample (e.g.
DEIMOS_2018_L519281), we can clearly measure the Ly$\alpha$ forest
transmission. In the case of a non-detection in NB718 (e.g.
SILVERRUSH_2021_3832) the derived constraint on $T_{\rm IGM}$ is consistent
with zero within the photometric uncertainty. The table of the MCMC results
for the full DEIMOS10k and LAE samples is available as the supplementary
online material.

Figure 7 shows the distribution of the estimated Ly$\alpha$ forest
transmission. The SED fitting provides $|\delta T_{\rm IGM}/T_{\rm IGM}|\sim
49\%$ and $76\%$ determinations of the Ly$\alpha$ forest transmission for the
DEIMOS10k and LAE samples with median $\delta T_{\rm IGM}\sim 0.10$ and
$0.42$, where $\delta T_{\rm IGM}$ and $T_{\rm IGM}$ are given by the standard
deviation and mean of the posterior. This includes errors from both
photometric noise and the uncertainty in the intrinsic UV continuum. Defining
the signal-to-noise ratio (SNR) to be the inverse of the relative error ${\rm
SNR}=|\delta T_{\rm IGM}/T_{\rm IGM}|^{-1}$, the median SNR of the individual
$T_{\rm IGM}$ measurements is thus ${\rm SNR}\simeq 2.0$ and $1.3$ for
DEIMOS10k and LAE samples, respectively.

The uncertainty in $T_{\rm IGM}$ for the LAEs is larger than the DEIMOS10k
sample. Since the LAE sample is typically fainter than the DEIMOS10k sample,
it is more severely affected by photometric noise because of (i) the reduced
contrast between the UV continuum level ($z$-band) and the Ly$\alpha$ forest
flux (NB718) and (ii) a less precise determination of the UV continuum slope
($z-y$ colour).

The uncertainty in the UV continuum slope introduces a degeneracy in the
estimate of the Ly$\alpha$ forest transmission. As discussed in Kakiichi et
al. (2022), for power-law spectra the uncertainty in the continuum slope
enters as $T_{\rm IGM}^{\rm estimated}=(\lambda_{\rm NB}/\lambda_{\rm
BB})^{\beta_{\rm true}-\beta_{\rm temp}}T_{\rm IGM}^{\rm true}$ where
$\beta_{\rm true}$ and $\beta_{\rm temp}$ are the continuum slopes of true and
template spectra and the central wavelengths of the filters are $\lambda_{\rm
NB}=7170.5\,\mbox{\AA}$ for NB718 and $\lambda_{\rm BB}=8912.6\,\mbox{\AA}$
for the $z$-band. The typical relative error $\delta\epsilon_{\rm
cont}\equiv|T_{\rm IGM}^{\rm true}-T_{\rm IGM}^{\rm estimated}|/T_{\rm
IGM}^{\rm true}$ is then estimated as:

$\langle\delta\epsilon_{\rm cont}\rangle=\int|1-(\lambda_{\rm NB}/\lambda_{\rm
BB})^{\beta_{\rm true}-\beta_{\rm temp}}|P(\beta_{\rm
true}|\bar{\beta},\sigma_{\beta})d\beta_{\rm true},$ (8)

resulting in $\langle\delta\epsilon_{\rm cont}\rangle\approx 11\,\%$
($21\,\%$) error for $\beta_{\rm temp}=-1.8$ ($-1.8\pm 1$) and assuming that
true $\beta$ slopes follows a Gaussian PDF with mean $\bar{\beta}=-1.8$ and
$\sigma_{\beta}=0.7$ consistent with Bouwens et al. (2014). In comparison, the
relative error due to the photometric noise in the Ly$\alpha$ forest
transmission can be estimated by approximating $T_{\rm IGM}\approx f_{\rm
NB}/f_{\rm z}$ as

$\delta\epsilon_{\rm phot}=\sqrt{(\delta f_{\rm NB718}/f_{\rm
NB718})^{2}+(\delta f_{\rm z}/f_{\rm z})^{2}},$ (9)

where $\delta f_{\rm NB718}$ and $\delta f_{\rm z}$ are the limiting fluxes in
NB718 and $z$. The median relative error of the DEIMOS10k and LAE samples are
$\sim 46\,\%$ and $\sim 89\,\%$. These estimates indicate that for our current
HSC depth, the photometric noise dominates the continuum slope uncertainty.

This point is further reinforced by the fact that the choice of flat or
Gaussian priors on the continuum slope has little impact on the resulting
distribution of individual $T_{\rm IGM}$ (Figure 7). A two-sample Kolmogorov-
Smirnov test indicates the difference is not statistically significant with
p-values of $p=0.392$ and $0.999$ for the LAE and DEIMOS10k samples
respectively.

####  4.2.1 Population synthesis vs power-law SEDs

We now test whether a power-law spectrum is sufficient to accurately represent
the intrinsic galaxy spectrum for the purposes of IGM tomography. We compare
the power-law spectrum with the intrinsic galaxy spectrum model without the
Ly$\alpha$ forest absorption from the stellar population synthesis code BPASS
(v2.2.1, Eldridge et al. 2017; Stanway & Eldridge 2018) processed with the
photionization code CLOUDY (c17.01 Ferland et al., 2017). We generate a grid
of model spectra with metallicities, $Z=0.1,0.2,0.3,0.5,1.0\rm\,Z_{\odot}$ and
stellar ages, $\log_{10}t_{\rm age}/{\rm yr}=6.0,6.5,7.0,7.5,8.0,8.5$, with a
Salpeter initial mass function with the upper mass limit set to $300\,\rm
M_{\odot}$ including binary stars. We assume an instantanous starburst to
model the LAEs and a continuous star formation history to model the Lyman-
break galaxies (LBGs) in the DEIMOS10k sample. We assume the stellar
population is spherically surrounded by gas with an electron density
$n_{e}=200\rm\,cm^{-3}$ and ionization parameter $\log_{10}U=-2.5$. We then
apply the Calzetti (2001) dust attenuation curve with $E(B-V)=0.0,0.1,0.2,0.4$
to the BPASS+CLOUDY outputs to model the intrinsic galaxy spectra.

In Figure 8 (left) we compare the LAE power-law spectra with a range of the
best-fit UV continuum slopes with intrinsic spectra calculate from the
BPASS+CLOUDY model with varying metallicities, ages, and dust extinctions. All
the spectra are normalized at $\sim 1330\,\mbox{\AA}$ corresponding to the
rest-frame wavelength coverage of the $z$-band. For a typical range of
metallicities $Z\sim 0.1-0.3\,Z_{\odot}$, stellar ages $t_{\rm age}\sim
1-30\rm\,Myr$, and dust extinction $E(B-V)\sim 0.0-0.3$ for LAEs (e.g. Ono et
al. 2010; Guaita et al. 2011; Nakajima et al. 2012; Hagen et al. 2014; Trainor
et al. 2016, see also reviews by Hayes 2019; Ouchi et al. 2020), the power-law
template approximates the continuum shape of the BPASS+CLOUDY spectra at $\sim
1000-1600\,\mbox{\AA}$ very well. The relative error in the estimated $T_{\rm
IGM}$ due to adopting different SEDs, $\delta\epsilon_{\rm SED}\equiv|(T_{\rm
IGM}^{\rm power-law}-T_{\rm IGM}^{\rm BPASS+CLOUDY})/T_{\rm IGM}^{\rm
BPASS+CLOUDY}|$, is given by

$\delta\epsilon_{\rm SED}=|1-f_{\rm NB}^{\rm intr,BPASS+CLOUDY}/f_{\rm
NB}^{\rm intr,power-law}|.$ (10)

Comparing the power-law template with median $\beta$ slope and the
BPASS+CLOUDY spectrum with typical LAE parameters of $Z=0.20\rm\,Z_{\odot}$,
$\log_{10}t_{\rm age}/{\rm yr}=7.0$, and $E(B-V)=0.10$, the relative error is
$\delta\epsilon_{\rm SED}\sim 6\,\%$ ($\sim 15$ and $17\,\%$ for $E(B-V)=0.2$
and $\log_{10}t_{\rm age}/{\rm yr}=7.5$ respectively). This is much smaller
than the error from photometric noise and be captured by the continuum slope
error budget of the power-law template. We conclude the power-law SED fit is
sufficient to predict the intrinsic flux at the Ly$\alpha$ forest region for
the current HSC depth.

We expect a similar uncertainty for the DEIMOS10k sample, which should
primarily consist of LBGs. Figure 8 (right) shows the same comparison
normalized at $\sim 1400\,\mbox{\AA}$ corresponding to the $z$-band coverage
at the mean redshift of the DEIMOS10k sample. LBGs typically span a range of
dust extinction $E(B-V)\sim 0.0-0.4$ (de Barros et al., 2014; Reddy et al.,
2016), metallicities $Z\sim 0.2-1.0\rm\,Z_{\odot}$ (Steidel et al., 2014), and
stellar ages $t_{\rm age}\sim 50-500\rm\,Myr$ (Stark et al., 2009; Curtis-Lake
et al., 2013; de Barros et al., 2014). Being more mature star-forming galaxies
than LAEs, their older stellar ages or higher dust extinctions introduce a
downward trend towards shorter wavelengths. For typical LBG parameters of
$Z=0.5\rm\,Z_{\odot}$, $\log_{10}\,t_{\rm age}/{\rm yr}=8.0$, and
$E(B-V)=0.2$, the impact on the estimated $T_{\rm IGM}$ between the
BPASS+CLOUDY and power-law SEDs is $\delta\epsilon_{\rm SED}\sim 27\rm\,\%$
($\sim 39$ and $31\rm\,\%$ for $E(B-V)=0.4$ and $\log_{10}\,t_{\rm age}/{\rm
yr}=8.5$ respectively). While larger than for the LAEs, this is still within
the photometric error for the individual measurements of $T_{\rm IGM}$.
However, if this range of stellar ages and dust extinctions is representative
of the true intrinsic spectra of background spec-z LBGs, it could introduce a
systematic bias in stacked measurements. The use of power-law SEDs for
background spec-z sample could systematically underestimate the measured
Ly$\alpha$ forest transmission because it predicts an intrinsic UV continuum
level larger than the true value. To quantify and eliminate the possible
impact of this limitation, we would need to extend our analysis (which
currently only uses $z$\- and $y$-bands) including near-infrared data to
better constrain the ages and dust extinctions of background galaxies. We
discuss this strategy further in Section 9, but in the following analysis we
only discuss the effect of this potential bias.

Another possible systematic arising from the assumption of a power-law
spectrum is the presence of stellar photospheric and interstellar absorption
lines in the Ly$\alpha$ forest region of a background galaxy. Prominant
absorption lines between Ly$\beta$ and Ly$\alpha$ lines
($1026-1216\,\mbox{\AA}$) include ${\rm C\,{\scriptstyle II}}\,\lambda 1036$,
${\rm S\,{\scriptstyle IV}}/{\rm Fe\,{\scriptstyle II}}\,\lambda 1063$, ${\rm
N\,{\scriptstyle II}}\,\lambda 1084$, ${\rm N\,{\scriptstyle I}}\,\lambda
1134$, ${\rm C\,{\scriptstyle III}}\,\lambda 1176$, ${\rm Si\,{\scriptstyle
II}}\,\lambda\lambda 1190,\\!1193$, ${\rm Si\,{\scriptstyle III}}\,\lambda
1207$ (Reddy et al., 2016). Previous spectroscopic IGM tomographic surveys
mitigated this issue by masking the $\pm 2-5\,\mbox{\AA}$ regions around each
absorption line (Lee et al., 2014b, 2018; Newman et al., 2020). For $z\sim
5.7$ background LAEs, the NB718 filter covers the rest-frame wavelength
between 1050 and 1082 Å, which coincides with the ${\rm S\,{\scriptstyle
IV}}/{\rm Fe\,{\scriptstyle II}}\,\lambda 1063$ line. According to a stacked
galaxy spectrum (Newman et al., 2020), the typical rest-frame equivalent width
of the absorption line is ${\rm EW}_{\rm abs}\sim 0.8\,\mbox{\AA}$. The effect
of an absorption line on the NB718 flux is thus $(1+z){\rm EW}_{\rm
abs}/\Delta\lambda_{\rm NB718}\approx 4.8\,\%$ reduction in flux integrated
over the NB718 filter. Accordingly, the resulting bias in the estimated
$T_{\rm IGM}$ is only

$\delta\epsilon_{\rm abs}=(1+z){\rm EW}_{\rm abs}/\Delta\lambda_{\rm NB718},$
(11)

i.e. $4.8\,\%$, which is negligible compared to the other sources of error
discussed above. For the DEIMOS10k sample, the spectroscopic redshifts span
the range of $4.98<z<5.89$. As the contamination in NB718 flux by the
absorption lines is expected to be randomized, the effect of absorption lines
on the overall DEIMOS10k sample should be negligible.

##  5 Mean Ly$\alpha$ forest transmission

###  5.1 Estimating mean Ly$\alpha$ forest transmission

The mean Ly$\alpha$ forest transmission is simply the mean in a representative
volume $V$, $\overline{T}_{\rm IGM}=V^{-1}\int T_{\rm IGM}(\bm{x})dV$ where
$\bm{x}$ is a 3D spatial position in the Universe. For IGM tomography, we can
only sample $T_{\rm IGM}(\bm{x})$ along sightlines to background galaxies. For
each background galaxy at a location $\bm{\theta}_{i}$, we sample $T_{\rm
IGM}(\bm{\theta}_{i})=T_{{\rm IGM},i}$ where $T_{{\rm IGM},i}$ is the measured
Ly$\alpha$ forest transmission integrated along the line of sight to an $i$th
galaxy over the width of the NB filter. Since the background galaxies are
uncorrelated with the foreground IGM structure, providing random sampling of
$T_{\rm IGM}(\bm{x})$, this is equivalent to performing the Monte Carlo
integration of the mean Ly$\alpha$ forest transmission, i.e.

$\overline{T}_{\rm IGM}=\frac{1}{N_{\rm bg}}\sum_{i=1}^{N_{\rm bg}}T_{{\rm
IGM},i}.$ (12)

Since we have noisy measurements of $T_{{\rm IGM},i}$ characterised by the
posterior $P(T_{{\rm IGM},i}|f_{{\rm NB},i}^{\rm obs},\bm{f}_{{\rm BB},i}^{\rm
obs})$ (equation 7) for a set of background galaxies $i=1,\dots,N_{\rm bg}$,
our estimate of the mean Ly$\alpha$ forest transmission is also noisy. Thus,
to characterise the uncertainty in the mean Ly$\alpha$ forest transmission, we
need to know the posterior probability of $\overline{T}_{\rm IGM}$, that is,
$P(\overline{T}_{\rm IGM}|\\{f_{{\rm NB},i}^{\rm obs},\bm{f}_{{\rm BB},i}^{\rm
obs}\\}_{i=1,\dots,N_{\rm bg}})$. Assuming individual measurements of $T_{{\rm
IGM},i}$ are statistically independent, the joint posterior probability is
simply the product of all the individual posteriors, i.e. $P(\\{T_{{\rm
IGM},i}\\}_{i=1,\dots,N_{\rm bg}}|\\{f_{{\rm NB},i}^{\rm obs},\bm{f}_{{\rm
BB},i}^{\rm obs}\\}_{i=1,\dots,N_{\rm bg}})=\prod_{i=1}^{N_{\rm bg}}P(T_{{\rm
IGM},i}|f_{{\rm NB},i}^{\rm obs},\bm{f}_{{\rm BB},i}^{\rm obs})$. Formally,
the posterior probability of the mean Ly$\alpha$ forest transmission can then
be written in terms of the convolution of the individual posteriors of
$T_{{\rm IGM},i}$ of all background galaxies (e.g. Sivia & Skilling, 2006,
Section 3.6),

$\displaystyle P(\overline{T}_{\rm IGM}|\\{f_{{\rm NB},i}^{\rm
obs},\bm{f}_{{\rm BB},i}^{\rm obs}\\}_{i=1,\dots,N_{\rm bg}})=$
$\displaystyle\int\prod_{i=1}^{N_{\rm bg}}dT_{{\rm IGM},i}P(T_{{\rm
IGM},i}|f_{{\rm NB},i}^{\rm obs},\bm{f}_{{\rm BB},i}^{\rm obs})\delta_{\rm
D}\left(\overline{T}_{\rm IGM}-\frac{1}{N_{\rm bg}}\sum_{i=1}^{N_{\rm
bg}}T_{{\rm IGM},i}\right),$ (13)

where $\delta_{\rm D}(x)$ is the Dirac delta function. Numerically, it is easy
to generate the posterior of $\overline{T}_{\rm IGM}$, i.e.
$P(\overline{T}_{\rm IGM}|\\{f_{{\rm NB},i}^{\rm obs},\bm{f}_{{\rm BB},i}^{\rm
obs}\\}_{i=1,\dots,N_{\rm bg}})$, by computing the histogram of many
$\overline{T}_{\rm IGM}$’s using random draws of $T_{{\rm IGM},i}$ from
$P(T_{{\rm IGM},i}|f_{{\rm NB},i}^{\rm obs},\bm{f}_{{\rm BB},i}^{\rm obs})$.

Equation (13) illustrates an important, but subtle point on the choice of
prior on $T_{{\rm IGM},i}$. When the data has no constraining power, i.e.
$P(T_{{\rm IGM},i}|f_{{\rm NB},i}^{\rm obs},\bm{f}_{{\rm BB},i}^{\rm
obs})\rightarrow P(T_{{\rm IGM},i})$, the posterior of $\overline{T}_{\rm
IGM}$ becomes a convolution of multiple priors of $T_{{\rm IGM},i}$. Assuming
any prior with mean $\overline{T}_{\rm prior}$ and standard deviation
$\sigma_{\rm prior}$ for all the individual measurements, by the virtue of
central limit theorem, the posterior of $\overline{T}_{\rm IGM}$ approaches a
Gaussian distribution with mean $\overline{T}_{\rm prior}$ and standard
deviation $\sigma_{\overline{T}_{\rm IGM}}=\sigma_{\rm prior}/\sqrt{N_{\rm
bg}}$. In our case, if we were to impose a flat prior between 0 and 1 (i.e.
$\overline{T}_{\rm prior}=0.5$ and $\sigma_{\rm prior}=0.29$) for the 104
individual measurements, the end result would be a Gaussian posterior on
$\overline{T}_{\rm IGM}$ with mean $\overline{T}_{\rm prior}=0.5$ and
$\sigma_{\overline{T}_{\rm IGM}}=0.29/\sqrt{104}=0.028$ even for a completely
uninformative dataset. This demonstrates that a reasonable prior on individual
measurements propagates into an unreasonably tight constraint on the mean
Ly$\alpha$ forest transmission. While imposing a flat prior with $0<T_{\rm
IGM}<1$ seems an innocent assumption, when we are interested in the mean we
should use a maximally non-informative prior on the individual measurement to
avoid an artificial constraint on the final estimate of $\overline{T}_{\rm
IGM}$. We thus use a flat prior on individual measurement allowing a very wide
range between $-100<T_{\rm IGM}<100$.

Given the full posterior probability of $\overline{T}_{\rm IGM}$, it is
natural to take our best estimate of the mean Ly$\alpha$ forest transmission
as the expectation value $\langle\overline{T}_{\rm
IGM}\rangle=\int\overline{T}_{\rm IGM}P(\overline{T}_{\rm IGM}|\\{f_{{\rm
NB},i}^{\rm obs},\bm{f}_{{\rm BB},i}^{\rm obs}\\}_{i=1,\dots,N_{\rm
bg}})d\overline{T}_{\rm IGM}$, which simplifies as the average of the
expectation values of individual measurements of $T_{{\rm IGM},i}$ (see
Appendix A for derivation),

$\langle\overline{T}_{\rm IGM}\rangle=\frac{1}{N_{\rm bg}}\sum_{i=1}^{N_{\rm
bg}}\langle T_{{\rm IGM},i}\rangle,$ (14)

where $\langle T_{{\rm IGM},i}\rangle=\int T_{{\rm IGM},i}P(T_{{\rm
IGM},i}|f_{{\rm NB},i}^{\rm obs},\bm{f}_{{\rm BB},i}^{\rm obs})dT_{{\rm
IGM},i}$. Similarly, the variance of the mean Ly$\alpha$ forest transmission,
$\sigma^{2}_{\overline{T}_{\rm IGM}}=\int(\overline{T}_{\rm
IGM}-\langle\overline{T}_{\rm IGM}\rangle)^{2}P(\overline{T}_{\rm
IGM}|\\{f_{{\rm NB},i}^{\rm obs},\bm{f}_{{\rm BB},i}^{\rm
obs}\\}_{i=1,\dots,N_{\rm bg}})d\overline{T}_{\rm IGM}$, is given by

$\displaystyle\sigma^{2}_{\overline{T}_{\rm IGM}}=\frac{1}{N_{\rm
bg}^{2}}\sum_{i=1}^{N_{\rm bg}}{\rm Var}[T_{{\rm IGM},i}].$ (15)

This again simply follows from the sum of the variances of individual
measurements, ${\rm Var}[T_{{\rm IGM},i}]=\int(T_{{\rm IGM},i}-\langle T_{{\rm
IGM},i}\rangle)^{2}P(T_{{\rm IGM},i}|f_{{\rm NB},i}^{\rm obs},\bm{f}_{{\rm
BB},i}^{\rm obs})dT_{{\rm IGM},i}$. Note that because the final variance
$\sigma^{2}_{\overline{T}_{\rm IGM}}$ scales like $1/N_{\rm bg}$ times the
average variance of the individual measurement, the final error scales as
$\sigma_{\overline{T}_{\rm IGM}}\propto 1/\sqrt{N_{\rm bg}}$. All these
quantities are easy to compute using the MCMC sample of $T_{{\rm IGM},i}$ of
the individual posteriors from the Bayesian SED fitting framework.

The central limit theorem guarantees that for a large number of background
galaxies $P(\overline{T}_{\rm IGM}|\\{f_{{\rm NB},i}^{\rm obs},\bm{f}_{{\rm
BB},i}^{\rm obs}\\}_{i=1,\dots,N_{\rm bg}})$ approaches a Gaussian
distribution. Thus the expectation value $\langle\overline{T}_{\rm
IGM}\rangle$ is equivalent to the the maximum a posteriori estimate of the
mean Ly$\alpha$ forest transmission, $\overline{T}_{\rm IGM}^{\rm MP}$, which
is given by

$\overline{T}_{\rm IGM}^{\rm
MP}=\underset{\scriptscriptstyle{\overline{T}}_{\rm IGM}}{\rm
arg\,max}~{}P(\overline{T}_{\rm IGM}|\\{f_{{\rm NB},i}^{\rm obs},\bm{f}_{{\rm
BB},i}^{\rm obs}\\}_{i=1,\dots,N_{\rm bg}}).$ (16)

Figure 9 shows the full posterior of $\overline{T}_{\rm IGM}$ directly
computed from random realizations of individual measurements, which explicitly
confirms the equivalence between the expectation value and the maximum a
posteriori estimate of $\overline{T}_{\rm IGM}$, i.e. $\overline{T}_{\rm
IGM}^{\rm MP}=\langle\overline{T}_{\rm IGM}\rangle$. As the full posterior is
Gaussian, the variance (equation 15) completely characterises the total error
in the estimated mean Ly$\alpha$ froest transmission as $\langle\bar{T}_{\rm
IGM}\rangle\pm\sigma_{\overline{T}_{\rm IGM}}$ at one sigma level. This error
is fully propagated, including both photometric noise and UV continuum
uncertainties, from the Bayesian SED fitting procedure for the individual
measurements of $T_{{\rm IGM},i}$.

Figure 9: Full posterior of the mean Ly$\alpha$ forest transmission
$\overline{T}_{\rm IGM}$, $P(\overline{T}_{\rm IGM}|\\{f_{{\rm NB},i}^{\rm
obs},\bm{f}_{{\rm BB},i}^{\rm obs}\\}_{i=1,\dots,N_{\rm bg}})$, for the LAE
(blue histogram) and DEIMOS10k (red histogram) samples. The histogram is
computed using the mean Ly$\alpha$ forest transmission from 1,000,000 random
realizations of a set of $T_{{\rm IGM},i}$ from the individual posteriors. The
blue and red solid curves are the Gaussian distributions with the expectation
value and variance computed via (14) and (15). Note that both the expectation
value and the maximum a posteriori estimate of $\overline{T}_{\rm IGM}$ are
equivalent.

###  5.2 Results

Using the DEIMOS10k and LAE samples, the mean Ly$\alpha$ forest transmission
at $z\simeq 4.9$ are estimated to be

$\langle\overline{T}_{\rm IGM}\rangle=0.179\pm 0.028~{}~{}~{}~{}\mbox{(for
DEIMOS10k)}$ (17)

and

$\langle\overline{T}_{\rm IGM}\rangle=0.293\pm 0.045~{}~{}~{}~{}\mbox{(for
LAE)}$ (18)

This represents the first $\sim 15\,\%$ photometric measurement of the mean
Ly$\alpha$ forest transmission of the IGM using background galaxies. The
statistical precision is comparable to the $\sim 8\,\%$ determinations based
using quasar spectra (Becker et al., 2013, 2015b; Eilers et al., 2018; Yang et
al., 2020; Bosman et al., 2018; Bosman et al., 2022). Our precision arises
from a large number of background galaxies roughly consistent with the
$\propto 1/\sqrt{N_{\rm bg}}$ scaling law, from which we expect $\delta T_{\rm
IGM}/T_{\rm IGM}\sim 0.49/\sqrt{36}$ $(0.76/\sqrt{104})$ $=8\,\%$ $(7\,\%)$
error on the mean transmission for the DEIMOS10k (LAE) sample. Note that the
error budget quoted above using equation (15) takes into account only the
photometric noise and UV continuum uncertainties, but not the error from
cosmic or patch-to-patch variance. To quantify its impact, we empirically
estimate the total error using Jackknife resampling (described in Section 6).
We find that the Jackknife errors are $\sigma_{\rm JK}=0.024$ and $0.047$ for
the DEIMOS10k and LAE samples respectively, comparable to our analytic
estimate of the (photometric noise + continuum) error. Thus, the additional
error from cosmic variance is not significant.

To visually confirm that the measured Ly$\alpha$ forest transmission is truly
representative of the physical value at $z\simeq 4.9$, we applied a sigma-
clipped mean stacking of the NB718 and BB images of the DEIMOS10k and LAE
sample in Figure 10. The signal is clearly detected in NB718 as well as via a
clear UV continuum detection in $z$. The signal is also detected in the median
stack. The non-detection in mean $g$ stack reinforces that the detected signal
originates from the Ly$\alpha$ forest transmission towards the background
galaxies and not due to low-redshift interlopers.

Figure 10: $5\arcsec\times 5\arcsec$ cutout images of the sigma-clipped mean
stack of $g$ (left), NB718 (middle), $z$ (right) images for the DEIMOS10k
(top) and LAE (bottom) samples.

####  5.2.1 Comparison with the literature

In Figure 11 we compare our mean Ly$\alpha$ forest transmission with the
measurements using quasars (Becker et al., 2013; Eilers et al., 2018; Bosman
et al., 2022) and galaxy spectra (Thomas et al., 2017, 2020, 2021) in the
literature. Using high signal-to-noise quasar spectra, the former measures the
mean Ly$\alpha$ forest transmission in bins of $50\,h^{-1}\rm cMpc$ length.
This is comparable to the line-of-sight comoving length $34\,h^{-1}\rm cMpc$
of the NB718 filter. Thomas et al. (2017, 2020, 2021) used a large
spectroscopic sample of galaxies from the VANDELS and VUDS surveys and
measured the Ly$\alpha$ forest transmission from the rest-frame
$1070-1170\,\mbox{\AA}$ region of background galaxies by using various IGM
templates in their spectral fitting method (see also Monzon et al. (2020) who
used a stacking method).

Our measurement using the DEIMOS10k sample is in excellent agreement with the
quasar studies, in particular with the latest high signal-to-noise measurement
of $\langle\overline{T}_{\rm IGM}\rangle=0.171\pm 0.014$ based on the XQR-30
quasar sample (Bosman et al., 2022). Although our measurement using background
LAEs is slightly higher, it is still broadly in agreement with the previous
values from the quasar- and galaxy-based measurements.

Our conclusion disagrees with the claim by Thomas et al. (2020) who argued
that photometric data is insufficient to constrain the Ly$\alpha$ forest
transmission. This is because their SED fitting used broad-band photometry to
determine the Ly$\alpha$ forest transmission which is contaminated by the
Ly$\alpha$ emission line and UV continuum and covers a region below the
Ly$\beta$ line depending on the redshift of a background galaxy. This renders
the resulting measurement of the Ly$\alpha$ forest transmission uncertain and
is a likely source of their $\sim 20\,\%$ discrepancy between their
photometric and spectroscopic measurements. In contrast, our method uses a NB
filter precisely covering the appropriate Ly$\alpha$ forest region enabling a
clean photometric measurement of the transmitted Ly$\alpha$ forest flux. We
argue that there is no fundamental limitation to the photometric approach when
a carefully chosen combination of a NB filter and background galaxy redshifts
is used.

Figure 11: Comparison of the mean Ly$\alpha$ forest transmission measured in
this work using the DEIMOS10k (red) and LAE (blue) samples with the previous
measurements based on quasar (Becker et al., 2013; Eilers et al., 2018; Bosman
et al., 2022) and galaxy spectra (Thomas et al., 2020, 2021).

###  5.3 Systematics: low-redshift interlopers and SEDs

There is a $\sim 2\sigma$ tension in our measurements of the mean Ly$\alpha$
forest transmission between the DEIMOS10k and LAE samples. The statistical
significance is calculated from the difference between the posterior means
divided by the quadrature sum of the errors (Lemos et al., 2021),
$(0.293-0.179)/\sqrt{0.045^{2}+0.028^{2}}=2.15\sigma$. While the discrepancy
is statistically insignificant, it could indicate systematic effects outside
our statistical error budget.

####  5.3.1 Low-redshift interlopers

The first possibility is contamination by the low-redshift interlopers in the
background $z\simeq 5.7$ LAE sample. A low-redshift interloper could bias the
result by introducing a fictitious transmissive sightline. Possible
interlopers in the LAE sample selected by a NB816 excess include low-redshift
galaxies with strong emission lines from $z\simeq 0.25$ H$\alpha$, $z\simeq
0.63$ $[{\rm O\,{\scriptstyle III}}]\,\lambda 5008$, and $z\simeq 1.19$ $[{\rm
O\,{\scriptstyle II}}]\,\lambda\lambda 3727,3729$, as well as slightly lower
redshift AGN at $z\simeq 4.28$ with ${\rm C\,{\scriptstyle IV}}\,\lambda 1549$
emission (e.g. Ouchi et al., 2008; Shibuya et al., 2018; Sobral et al., 2018).
Their rest-frame optical or UV continua could be mistaken in the NB718 filter
as the transmitted Ly$\alpha$ forest flux at $z\simeq 4.9$, artificially
increasing the estimated mean Ly$\alpha$ forest transmission. The effect of a
low-redshift interloper can be written as

$\langle\overline{T}_{\rm IGM}\rangle=(1-f_{\rm
bg.int})\langle\overline{T}_{\rm IGM}\rangle^{\rm true}+f_{\rm
bg.int}\langle\overline{T}_{\rm IGM}\rangle^{\rm bg.int},$ (19)

where $f_{\rm bg.int}$ is the contamination rate by low-redshift interlopers
in the background LAE sample and $\langle\overline{T}_{\rm IGM}\rangle^{\rm
bg.int}$ is the fictitious mean Ly$\alpha$ forest transmission along the
sightlines of the interlopers. The interloper fraction in the LAE selection is
typically $\sim 20\%$ (Shibuya et al., 2018). Assuming that
$\langle\overline{T}_{\rm IGM}\rangle^{\rm bg.int}=0.7$ and the measured mean
Ly$\alpha$ forest transmission from DEIMOS10k sample is the true value
$\langle\overline{T}_{\rm IGM}\rangle^{\rm true}=0.179$, the observed value
using NB-selected background LAEs will become $\langle\overline{T}_{\rm
IGM}\rangle=(1-0.2)\times 0.179+0.2\times 0.7=0.283$, being consistent with
the measurement from our LAE sample. This would resolve the tension between
DEIMOS10k and LAE samples.

In order to explore this further, we cross-matched our background LAE
catalogue with the spectroscopic catalogue compiled with HSC-SSP DR3 (Aihara
et al., 2022) containing 70,358 objects in the UD-COSMOS field (tract 9813).
We find no interloper in our background LAE catalogue while 10 objects are
spectroscopically-confirmed to be at $z\simeq 5.7$. We also applied a stricter
$g$-band non-detection cut compared to our default $3\sigma$ threshold and
find that for a $1\sigma$ ($2\sigma$) threshold the estimated mean Ly$\alpha$
forest transmission becomes $\langle\overline{T}_{\rm IGM}\rangle=0.277\pm
0.055$ $(0.298\pm 0.048)$ for LAE sample; the values are consistent with our
result from the $3\sigma$ threshold within $1\sigma$ error. Thus we find no
obvious evidence for low-redshift interlopers in our LAE sample.

Note that as we require $>5\sigma$ detection in the $z$-band
($<26.86\rm\,mag$), potential interlopers, if any, need to have red
$g-z\gtrsim 2.0$ colours and be fainter than $29.60-28.84$ mag $(1-2\sigma)$
in $g$-band. Such interlopers could be low-redshift dusty red galaxies or
Balmer break galaxies with H$\alpha$, $[{\rm O\,{\scriptstyle III}}]$, or
$[{\rm O\,{\scriptstyle II}}]$ doublet emission lines at $z\simeq
0.25,0.63,1.19$. Ultimately, spectroscopic follow-up of the background LAEs
would be necessary to fully reject systematic bias from low-redshift
interlopers. Following up a random subset would determine the interloper
fraction $f_{\rm bg.int}$ and by measuring the mean Ly$\alpha$ forest
transmission along the interlopers, one can also determine
$\langle\overline{T}_{\rm IGM}\rangle^{\rm bg.int}$. Then the observed
Ly$\alpha$ forest transmission using the parent photometric background LAE
sample, $\langle\overline{T}_{\rm IGM}\rangle^{\rm obs}$, could be
statistically corrected via

$\langle\overline{T}_{\rm IGM}\rangle^{\rm
corrected}=\frac{\langle\overline{T}_{\rm IGM}\rangle^{\rm obs}-f_{\rm
bg.int}\langle\overline{T}_{\rm IGM}\rangle^{\rm bg.int}}{1-f_{\rm bg.int}}.$
(20)

####  5.3.2 SED templates

As discussed in Section 4.2.1, model galaxy SEDs may introduce a systematic
error in $T_{\rm IGM}$. The effect is expected to be larger in the DEIMOS10k
sample which contains more mature galaxies with dust or more complex stellar
populations for which the UV continuum cannot precisely be determined by only
$z$ and $y$-band photometry. If the intrinsic continua were systematically
overestimated by $\sim 30\,\%$ (for example, because of an underestimated
$E(B-V)$), then a more flexible galaxy SED template would lead to
$\langle\overline{T}_{\rm IGM}\rangle=0.233\pm 0.036$ relaxing the tension
between the DEIMOS10k and LAE samples to $\sim 1\sigma$. However, this
explanation would weaken the agreement between both our measures and those
determined using quasars.

Our choice of the Ly$\alpha$ forest wavelength range $(1026-1216\,\mbox{\AA})$
is more generous compared to the $1040-1190\,\mbox{\AA}$ range commonly used
for spectroscopic IGM tomography (Lee et al., 2014b, 2018; Newman et al.,
2020). Although the effect of absorption lines is small (Section 4.2.1), the
wider range means that a power-law continuum might neglect the effect of
Ly$\beta$+${\rm C\,{\scriptstyle II}}\,\lambda 1036$ and Ly$\alpha$+${\rm
Si\,{\scriptstyle III}}\,\lambda 1207$ absorption lines. Also, neutral gas in
the circumgalactic medium of a background galaxy could contribute to the
Ly$\alpha$ absorption blueward of the line centre (Rudie et al., 2013;
Kakiichi et al., 2018; Bassett et al., 2021). We can test this effect by
adopting a more restrictive range of $1040-1190\,\mbox{\AA}$ and limiting the
redshift range of background galaxies to $5.072<z<5.841$ for the DEIMOS10k
sample. We find $\langle\overline{T}_{\rm IGM}\rangle=0.174\pm 0.035$,
consistent with our main result. Thus contamination by absorption lines is
negligible and cannot explain the tension.

####  5.3.3 Cosmic variance

Finally, cosmic variance may cause the tension between our two subsamples.
Although the Jackknife resampling error should include the effects of cosmic
variance, the small sample size may underestimate the effect. However, since
both the DEIMOS10k and LAE samples probe similar regions of the sky in the
COSMOS field, we consider cosmic variance an unlikely source of the tension.

In conclusion, we believe that the combination of a modest contribution from
low-redshift interlopers and some variation in the SED templates is the mostly
likely source of the $2\sigma$ tension between the DEIMOS10k and LAE samples.
This can only be tested and corrected by spectroscopic follow-up of the
background LAEs and inclusion of near-infrared photometry to better constrain
the background galaxy SEDs.

##  6 LAE-Ly$\alpha$ forest cross-correlation

Figure 12: The full posterior of the mean Ly$\alpha$ forest transmission
around $z\simeq 4.9$ LAEs, $P(\overline{T}_{\rm IGM}(\theta)|\\{f_{{\rm
NB},i}^{\rm obs},\bm{f}_{{\rm BB},i}^{\rm obs}\\}_{i=1,\dots,N_{\rm bg}})$,
using background $z\simeq 5.7$ LAE (blue histogram) and DEIMOS10k (red
histogram) samples. The histogram is computed via $\overline{T}_{\rm
IGM}(\theta)$ from 10,000 random realizations of a set of $T_{{\rm IGM},i}$
from the individual posteriors. The innermost (left), middle (middle),
outermost (rightl) angular bins are shown. The blue and red solid curves are
the Gaussian distributions with the expectation value and variance computed
via equations (22) and (23). The equivalence between the expectation value and
the maximum a posteriori estimate of the mean Ly$\alpha$ forest transmission
around LAEs is verified.

###  6.1 Estimating the mean Ly$\alpha$ forest transmission around LAEs

The NB718 dataset can be used to reveal $z\simeq 4.9$ LAEs in the same
redshift slice as the IGM Ly$\alpha$ forest transmission and thus how the
transmission varies as a function of angular separation $\theta$ from the
foreground LAEs. As before, we can use a Monte Carlo sampling to evaluate the
angular mean Ly$\alpha$ forest transmission around LAEs, $\overline{T}_{\rm
IGM}(\theta)=V(\theta)^{-1}\int T_{\rm IGM}(\bm{x})dV(\theta)$. By angular
averaging the transmission for all pairs of foreground LAEs and background
galaxy sightlines, we obtain

$\overline{T}_{\rm IGM}(\theta)=\frac{1}{N_{\rm
pair}(\theta)}\sum_{j=1}^{N_{\rm fg}}\sum_{i=1}^{N_{\rm bg}}T_{{\rm
IGM},i}\mathcal{I}(|\theta-\theta_{ij}|),$ (21)

where $N_{\rm pair}(\theta)=\sum_{j=1}^{N_{\rm fg}}\sum_{i=1}^{N_{\rm
bg}}\mathcal{I}(|\theta-\theta_{ij}|)$ is the number of pairs in each angular
bin and $\mathcal{I}(|\theta-\theta_{ij}|)$ is an indicator function equal to
unity when the angular separation
$\theta_{ij}=|\bm{\theta}_{i}-\bm{\theta}_{j}|$ between $i$ and $j$ is in the
angular bin specified by $\theta$ with width $\Delta\theta$, i.e.
$\mathcal{I}(|\theta-\theta_{ij}|)=1$ if $|\theta-\theta_{ij}|<\Delta\theta$,
and zero otherwise. The sums run over all foreground LAEs $j=1,\dots,N_{\rm
fg}$ and all background galaxy sightlines $i=1,\dots,N_{\rm bg}$.

Again as each background galaxy sightline provides a noisy measurement of
$T_{{\rm IGM},i}$, analogous to the argument made for the mean Ly$\alpha$
forest transmission, the full posterior probablity of the mean Ly$\alpha$
forest transmission around LAEs, $P(\overline{T}_{\rm IGM}(\theta)|\\{f_{{\rm
NB},i}^{\rm obs},\bm{f}_{{\rm BB},i}^{\rm obs}\\}_{i=1,\dots,N_{\rm bg}})$,
can be expressed in terms of the posteriors of individual $T_{{\rm IGM},i}$
measurements, which can be numerically computed by randomly sampling the
individual posteriors. The expectation value of the angular mean Ly$\alpha$
forest transmission around LAEs is given by (see Appendix A)

$\displaystyle\langle\overline{T}_{\rm IGM}(\theta)\rangle=\frac{1}{N_{\rm
pair}(\theta)}\sum_{j=1}^{N_{\rm fg}}\sum_{i=1}^{N_{\rm bg}}\langle T_{{\rm
IGM},i}\rangle\mathcal{I}(|\theta-\theta_{ij}|).$ (22)

The variance of the estimated mean Ly$\alpha$ forest transmission around LAEs
at each angular bin is given by

$\displaystyle{\rm Var}[\overline{T}_{\rm IGM}(\theta)]=\frac{1}{N_{\rm
pair}(\theta)^{2}}\sum_{j=1}^{N_{\rm fg}}\sum_{i=1}^{N_{\rm bg}}{\rm
Var}[T_{{\rm IGM},i}]\mathcal{I}(|\theta-\theta_{ij}|).$ (23)

Figure 12 verifies that the full posterior $P(\overline{T}_{\rm
IGM}(\theta)|\\{f_{{\rm NB},i}^{\rm obs},\bm{f}_{{\rm BB},i}^{\rm
obs}\\}_{i=1,\dots,N_{\rm bg}})$ follows a Gaussian distribution, which can be
fully characterised by an expectation value (22) and variance (23). The
maximum a posteriori estimation,

$\overline{T}^{\rm MP}_{\rm
IGM}(\theta)=\underset{\overline{T}_{\rm\scriptscriptstyle IGM}(\theta)}{\rm
arg\,max}~{}P(\overline{T}_{\rm IGM}(\theta)|\\{f_{{\rm NB},i}^{\rm
obs},\bm{f}_{{\rm BB},i}^{\rm obs}\\}_{i=1,\dots,N_{\rm bg}}),$ (24)

is therefore equivalent to the expectation value $\langle\overline{T}_{\rm
IGM}(\theta)\rangle$. The error estimated by the square root of the variance
(23) in the mean Ly$\alpha$ forest transmission around LAEs scales as $\propto
1/\sqrt{N_{\rm pair}(\theta)}$ and includes the full uncertainties from
photometric noise and continuum error in our individual $T_{\rm IGM}$
measurements from the Bayesian SED fitting framework.

###  6.2 The LAE-Ly$\alpha$ forest cross-correlation function

While the “mean Ly$\alpha$ forest transmission around LAEs” is well defined
given the observed distribution of LAEs, in order to examine the angular
“cross-correlation”, which is given by,

$\omega_{\rm g\alpha}(\theta)=\langle\overline{T}_{\rm
IGM}(\theta)\rangle/\langle\overline{T}_{\rm IGM}\rangle-1,$ (25)

we must quantify whether the excess probability of finding galaxies in the
environment of high or low Ly$\alpha$ forest transmission is statistically
significant. An additional uncertainty arises from the Poisson sampling of
foreground galaxies in the survey footprint, including the effect of mask
regions and the edge of the field-of-view.

To understand the scatter, we generate a random galaxy catalogue. We populate
$N_{\rm rand}$ objects at random locations $\\{\bm{\theta}_{i}^{\rm
rand}\\}_{i=1,\dots,N_{\rm rand}}$ within the survey footprint excluding the
actual masked regions. We then repeat the measurement of the mean Ly$\alpha$
forest transmission around the random objects using the actual background
sightlines,

$\langle\overline{T}_{\rm IGM}^{\rm random}(\theta)\rangle=\frac{1}{N_{\rm
pair}(\theta)}\sum_{j=1}^{N_{\rm rand}}\sum_{i=1}^{N_{\rm bg}}\langle T_{{\rm
IGM},i}\rangle\mathcal{I}(|\theta-\theta_{ij}^{\rm rand}|),$ (26)

where $\theta_{ij}^{\rm rand}$ is the angular separation between $i$-th random
galaxy position and the observed location of $j$-th background galaxy
sightline, $\theta_{ij}^{\rm rand}=|\bm{\theta}_{i}^{\rm
rand}-\bm{\theta}_{j}|$. This de-correlates the real cross-correlation signal
between LAEs and Ly$\alpha$ forest transmission and should converge to the
mean Ly$\alpha$ forest transmission, $\langle\overline{T}^{\rm random}_{\rm
IGM}(\theta)\rangle\approx\langle\overline{T}_{\rm IGM}\rangle$.

One can also de-correlate by randomly shuffling the observed values of
$T_{{\rm IGM},i}$ among the background galaxy sightlines but keeping their
angular locations fixed,

$\langle\overline{T}_{\rm IGM}^{\rm shuffle}(\theta)\rangle=\frac{1}{N_{\rm
pair}(\theta)}\sum_{j=1}^{N_{\rm fg}}\sum_{i=1}^{N_{\rm bg}}\langle T_{{\rm
IGM},i}^{\rm shuffle}\rangle\mathcal{I}(|\theta-\theta_{ij}|),$ (27)

where $\langle T_{{\rm IGM},i}^{\rm shuffle}\rangle$ is a random draw from a
set of real measurements $\\{\langle T_{\rm IGM,1}\rangle,\langle T_{\rm
IGM,2}\rangle,\dots,\langle T_{{\rm IGM},N_{\rm bg}}\rangle\\}$ without
replacement. The shuffled approach is convenient as it does not require us to
model the galaxy selection function. This serves to verify that the observed
LAE-Ly$\alpha$ forest cross-correlation is uncontaminated by the particular
distribution of background galaxies on the sky. Both the shuffled and random
measurements should be statistically identical, $\langle\overline{T}_{\rm
IGM}^{\rm shuffle}(\theta)\rangle\approx\langle\overline{T}_{\rm IGM}^{\rm
random}(\theta)\rangle$, and should converge to the mean Ly$\alpha$ forest
transmission $\langle\overline{T}_{\rm IGM}\rangle$ within the statistical
error.

Figure 13: Angular regions (coloured patches) used to compute the Jackknife
covariance matrix. Our survey footprint of the UD-COSMOS field are subdivided
into $N_{\rm JK}=20$ regions. Figure 14: Correlation coefficients of the
Jackknife covariance matrix for the LAE (left) and DEIMOS10k (right) samples.
The covariance matrix at the innermost bin for the DEIMOS10k sample is not
determined due to the small sample size in that bin.

Figure 15: Mean Ly$\alpha$ forest transmission profiles around $z=4.9$ LAEs
(black circles) using the $z=5.7$ LAE (left) and DEIMOS10k (right) samples as
background sources. The same profiles around random foreground objects (red
squares) and around foreground LAEs but using shuffled $T_{\rm IGM}$ along the
background sources (blue triangles) are offset by $\pm 0.02$ dex offset along
the x-axis for clarity. The gray shaded region indicates the mean Ly$\alpha$
forest transmission and its $1\sigma$ error. Top panels indicate the number of
foreground LAE - background sightline pairs for the LAE and DEIMOS10k samples
respectively.

We estimate the error on the cross-correlation using the Jackknife estimator
(e.g. Norberg et al., 2009). The Jackknife covariance matrix is given by

$\displaystyle{\rm Cov}_{\rm JK}\left[\omega_{\rm g\alpha}(\theta),\omega_{\rm
g\alpha}(\theta^{\prime})\right]=$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}\frac{N_{\rm JK}-1}{N_{\rm
JK}}\sum_{k=1}^{N_{\rm JK}}\left[\omega^{k}_{\rm
g\alpha}(\theta)-\overline{\omega}_{\rm g\alpha}^{\rm
JK}(\theta)\right]\left[\omega_{\rm
g\alpha}^{k}(\theta^{\prime})-\overline{\omega}_{\rm g\alpha}^{\rm
JK}(\theta^{\prime})\right],$ (28)

where

$\overline{\omega}_{\rm g\alpha}^{\rm JK}(\theta)=\frac{1}{N_{\rm
JK}}\sum_{k=1}^{N_{\rm JK}}\omega^{k}_{\rm g\alpha}(\theta),$ (29)

is the average of the cross-correlation functions from Jackknife resampling.
Jackknife regions are obtained using the k-means clustering
algorithm777https://github.com/esheldon/kmeans_radec (Kwan et al., 2017) on a
random galaxy catalogue with $N_{\rm rand}=50,000$. This algorithm subdivides
the observed survey area into $N_{\rm JK}$ regions of a roughly equal area as
shown in Figure 13. To compute the Jackknife covariance, we omit foreground
LAEs and $T_{\rm IGM}$ along background galaxy sightlines located in each
Jackknife region at a time and compute $N_{\rm JK}$ Jackknife re-sampled
cross-correlation functions $\omega^{k}_{\rm g\alpha}(\theta)$,
$k=1,\dots,N_{\rm JK}$, using the remaining objects. We use the same procedure
to compute the Jackknife error for the other summary statistics (the mean
Ly$\alpha$ forest transmission and Ly$\alpha$ forest auto-correlation
function) in this paper.

The correlation coefficients of the Jackknife covariance matrix are shown in
Figure 14. The covariance matrix of the innermost bins of the DEIMOS10k sample
could not be determined due to the small sample size. For both LAE and
DEIMOS10k samples, there are significant off-diagonal correlations between
angular bins as the same sightlines contribute multiple foreground LAE-
background sightline pairs.

###  6.3 Result

Figure 15 shows the angular mean Ly$\alpha$ forest transmission around LAEs at
$z\simeq 4.9$. We find no excess transmission or absorption in the NB-
integrated Ly$\alpha$ forest around the LAEs. The result is consistent with
the global mean $\langle\overline{T}_{\rm IGM}\rangle$ within the $2\sigma$
error both for the LAE and DEIMOS10k samples. We compare our result with the
random and shuffled measurements using the same number of foreground LAEs and
background galaxies in the real data. Both measurements show similar
fluctuations with the observed values, confirming that our result is
consistent with no spatial correlation.

####  6.3.1 Correcting for contamination by low-redshift interlopers

Figure 15 shows a mean offset between $\langle\overline{T}_{\rm
IGM}(\theta)\rangle$ measured using the LAE and DEIMOS10k samples. As
discussed in Section 5.3, this is likely caused by the low-redshift
interlopers in both the foreground and background LAE samples (Grasshorn
Gebhardt et al., 2019; Farrow et al., 2021). Since the distribution of any
low-redshift interlopers would be random relative to structures in the
tomographic slice of interest, the interlopers will dilute the observed cross-
correlation. Assuming foreground and background contamination fractions
$f_{\rm fg.int}$ and $f_{\rm bg.int}$, the observed angular mean Ly$\alpha$
forest transmission around the foreground LAEs can be expressed as (see
Appendix B)

$\displaystyle\langle\overline{T}_{\rm IGM}(\theta)\rangle=$
$\displaystyle\,(1-f_{\rm fg.int})(1-f_{\rm bg.int})\langle\overline{T}_{\rm
IGM}(\theta)\rangle^{\rm true}$ $\displaystyle+f_{\rm fg.int}(1-f_{\rm
bg.int})\langle\overline{T}_{\rm IGM}\rangle^{\rm true}+f_{\rm
bg.int}\langle\overline{T}_{\rm IGM}\rangle^{\rm bg.int},$ (30)

where $\langle\overline{T}_{\rm IGM}\rangle^{\rm true}$ is the true mean
Ly$\alpha$ forest transmission and $\langle\overline{T}_{\rm IGM}\rangle^{\rm
bg.int}$ is the fictitious mean Ly$\alpha$ forest transmission measured along
the low-redshift interlopers in the background LAE sample. The second and
third terms indicate contaminations from the low-redshifts interlopers in the
foreground and background LAE samples.

Following our definition of the observed LAE-Ly$\alpha$ forest cross-
correlation $\omega_{\rm g\alpha}(\theta)=\langle\overline{T}_{\rm
IGM}(\theta)\rangle/\langle\overline{T}_{\rm IGM}\rangle-1$, we can similarly
find that low-redshift interlopers dilute the cross-correlation amplitude by

$\displaystyle\omega_{\rm g\alpha}(\theta)=\frac{(1-f_{\rm fg.int})(1-f_{\rm
bg.int})\langle\overline{T}_{\rm IGM}\rangle^{\rm true}}{(1-f_{\rm
bg.int})\langle\overline{T}_{\rm IGM}\rangle^{\rm true}+f_{\rm
bg.int}\langle\overline{T}_{\rm IGM}\rangle^{\rm bg.int}}\omega^{\rm
true}_{\rm g\alpha}(\theta).$ (31)

The offset in the mean Ly$\alpha$ forest transmission around LAEs between the
background LAE and DEIMOS10k samples can be explained by this effect. As in
Section 5.3, we set $f_{\rm fg.int}=f_{\rm bg.int}=0.2$ both for foreground
and background LAE samples. The contamination fraction for the DEIMOS10k
sample is $f_{\rm bg.int}=0.0$ as all are confirmed spectroscopically. We
assume that the fictitious mean Ly$\alpha$ forest transmission along the
interlopers in the background LAE sample is $\langle\overline{T}_{\rm
IGM}\rangle^{\rm bg.int}=0.7$. As before, LAE interlopers can explain the
offset between the background LAE and DEIMOS10k samples.

These interlopers depress the observed LAE-Ly$\alpha$ forest cross-correlation
by $\omega_{\rm g\alpha}(\theta)\approx 0.40\,\omega^{\rm true}_{\rm
g\alpha}(\theta)$ and $\omega_{\rm g\alpha}(\theta)\approx 0.80\,\omega^{\rm
true}_{\rm g\alpha}(\theta)$ for background LAE and DEIMOS10k samples,
respectively. The cross-correlation measurement using the DEIMOS10k sample is
also affected because of the interloper contamination in the foreground LAEs.

All the terms in the damping pre-factor in Equation 31 can be determined and
statistically corrected a posteriori by spectroscopic follow up of a random
subset of the foreground and background LAE samples as described in Section
5.3.

####  6.3.2 Limit on the IGM fluctuations around LAEs

Figure 16 shows the observed LAE-Ly$\alpha$ forest cross-correlation at
$z\simeq 4.9$ after correcting for possible low-redshift interloper
contamination. To place an empirical constraint, we assume a simple power-law
form,

$\omega_{\rm g\alpha}^{\rm model}(\theta)=A_{0}(\theta/\theta_{0})^{-\gamma},$
(32)

where the fluctuations are characterised by the amplitude $A_{0}$ at angular
distance $\theta_{0}$ and the power-law slope $\gamma$. We assume a Gaussian
likelihood with the measured Jackknife covariance matrix and a fixed slope of
$\gamma=0.5$. The resulting $3\sigma$ bounds are shown in Figure 16. The
$3\sigma$ lower and upper limits are

$-0.29\left(\frac{r_{\perp}}{10\,h^{-1}\rm cMpc}\right)^{-0.5}<\omega_{\rm
g\alpha}^{\rm model}<0.07\left(\frac{r_{\perp}}{10\,h^{-1}\rm
cMpc}\right)^{-0.5}$ (33)

for the DEIMOS10k sample, and

$-0.58\left(\frac{r_{\perp}}{10\,h^{-1}\rm cMpc}\right)^{-0.5}<\omega_{\rm
g\alpha}^{\rm model}<0.40\left(\frac{r_{\perp}}{10\,h^{-1}\rm
cMpc}\right)^{-0.5}$ (34)

for the LAE sample. The derived lower and upper limits are consistent for both
the LAE and DEIMOS10k samples. While the bound from the DEIMOS10k sample is
slightly shifted to the negative cross-correlation, this is likely due to an
underestimated Jackknife error at small angular bins. The directly propagated
error (equation 23) from the Bayesian SED fitting framework indicate the error
from photometric noise and UV continuum uncertainty in the DEIMOS10k sample at
the inner bins should be larger than the empirical estimate from the Jackknife
method.

Our result indicates that the angular fluctuations of the Ly$\alpha$ forest
transmission around LAEs should be $\lesssim 58\,\%$ at $10\,h^{-1}\rm Mpc$
relative to the global mean at $z\simeq 4.9$. The physical interpretation of
the result will be discussed in a companion paper (Kakiichi et al in prep).

Figure 16: Observed LAE-Ly$\alpha$ forest cross-correlation function for the
DEIMOS10k (red squares) and LAE (blue circles) samples after correcting for
likely interloper contamination. The derived $3\sigma$ lower and upper limits
of the cross-correlation assuming the power-law with slope $\gamma=0.5$ are
shown with the red and blue shaded regions for the DEIMOS10k and LAE sample,
respectively. The Jackknife covariance matrices scaled by the interloper
correction factors are used to estimate the error.

##  7 Ly$\alpha$ forest auto-correlation

###  7.1 Estimating the Ly$\alpha$ forest auto-correlation

Figure 17: Ly$\alpha$ forest transmission auto-correlation function at $z=4.9$
using the $z=5.7$ LAE (left) and DEIMOS10k (right) samples. Values computed
using only shuffled background sources (red squares) and the cross-correlation
between data and shuffled samples (blue triangles) are shown offset $\pm 0.02$
dex along the x-axis for clarity. The gray shaded region indicates the
estimate in the case of no correlation based on the mean Ly$\alpha$ forest
transmission and its $1\sigma$ error. Top panels indicate the number of
sightline pairs for the LAE and DEIMOS10k samples.

We now turn our attention to examine the spatial fluctuations of Ly$\alpha$
forest transmission. These are expected to spatially correlate across
different sightlines due to large-scale fluctuations of the IGM. Unlike the
measurement of the 3D Ly$\alpha$ forest auto-correlation from spectra (e.g.
Slosar et al., 2011), photometric IGM tomography measures the angular auto-
correlation of Ly$\alpha$ forest transmission integrated over the line-of-
sight width of the NB filter ($\simeq 34\,h^{-1}\rm cMpc$).

In order to estimate the Ly$\alpha$ forest angular auto-correlation function,
using the pairs of NB-integrated Ly$\alpha$ forest transmission measurements,
we first compute, at each angular bin,

$T_{\rm IGM}T_{\rm IGM}(\theta)=\frac{1}{N_{\rm
pair}(\theta)}\sum_{i=1}^{N_{\rm bg}}\sum_{j>i}^{N_{\rm bg}}T_{{\rm
IGM},i}T_{{\rm IGM},j}\mathcal{I}(|\theta-\theta_{ij}|),$ (35)

where $N_{\rm pair}(\theta)=\sum_{i=1}^{N_{\rm bg}}\sum_{j>i}^{N_{\rm
bg}}\mathcal{I}(|\theta-\theta_{ij}|)$ is the number of pairs in each angular
bin. For independent measurements of $T_{{\rm IGM},i}$, the expectation value
of the Ly$\alpha$ forest angular auto-correlation is given by

$\langle T_{\rm IGM}T_{\rm IGM}(\theta)\rangle=\frac{1}{N_{\rm
pair}(\theta)}\sum_{i=1}^{N_{\rm bg}}\sum_{j>i}^{N_{\rm bg}}\langle T_{{\rm
IGM},i}\rangle\langle T_{{\rm
IGM},j}\rangle\mathcal{I}(|\theta-\theta_{ij}|).$ (36)

The error is computed from the Jackknife covariance matrix. The Ly$\alpha$
forest auto-correlation function is estimated by

$\omega_{\rm\alpha\alpha}(\theta)=\langle T_{\rm IGM}T_{\rm
IGM}(\theta)\rangle/\langle\overline{T}_{\rm IGM}\rangle^{2}-1.$ (37)

We can understand the scatter of $\langle T_{\rm IGM}T_{\rm
IGM}(\theta)\rangle$ in the absence of any spatial correlation. As the
Ly$\alpha$ forest auto-correlation has no complication from the window
function or the survey geometry, $\langle T_{\rm IGM}T_{\rm
IGM}(\theta)\rangle$ should be equal to $\langle\overline{T}_{\rm
IGM}\rangle^{2}$ if there is no spatial correlation. To test this, we de-
correlate the observed correlation by shuffling either one or both of the
measured values of $T_{\rm IGM}$ between the observed locations, i.e.

$\displaystyle\langle T_{\rm IGM}T_{\rm IGM}^{\rm shuffle}(\theta)\rangle=$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\frac{1}{N_{\rm
pair}(\theta)}\sum_{i=1}^{N_{\rm bg}}\sum_{j>i}^{N_{\rm bg}}\langle T_{{\rm
IGM},i}\rangle\langle T_{{\rm IGM},j}^{\rm
shuffle}(\theta)\rangle\mathcal{I}(|\theta-\theta_{ij}|),$ (38)

or

$\displaystyle\langle T_{\rm IGM}^{\rm shuffle}(\theta)T_{\rm IGM}^{\rm
shuffle}(\theta)\rangle=$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\frac{1}{N_{\rm
pair}(\theta)}\sum_{i=1}^{N_{\rm bg}}\sum_{j>i}^{N_{\rm bg}}\langle T_{{\rm
IGM},i}^{\rm shuffle}(\theta)\rangle\langle T_{{\rm IGM},j}^{\rm
shuffle}(\theta)\rangle\mathcal{I}(|\theta-\theta_{ij}|).$ (39)

Using the real set of $\\{T_{{\rm IGM},i}\\}_{i=1,\dots,N_{\rm bg}}$, we
generated a randomized set of $T_{\rm IGM}$ values keeping the angular
positions of the sightlines the same. This artificially de-correlates the
possible correlation. If there is no systematic, this should approach
$\simeq\langle\overline{T}_{\rm IGM}\rangle^{2}$.

###  7.2 Result

Figure 17 shows the observed auto-correlation function of the Ly$\alpha$
forest transmission at $z\simeq 4.9$. The observed auto-correlation is
consistent with the square of the mean and the shuffled results within
$2\sigma$ error, indicating the observed signal is consistent with no auto-
correlation. This null detection can be interpreted as the observed limit on
the Ly$\alpha$ forest transmission fluctuations at $z\simeq 4.9$.

####  7.2.1 Correcting for the contamination by low-redshift interlopers

Similar to the angular mean Ly$\alpha$ forest transmission around LAEs, Figure
17 shows an offset between $\langle T_{\rm IGM}T_{\rm IGM}(\theta)\rangle$
measured using the background LAE and DEIMOS10k samples, which is likely
caused by the low-redshift interlopers. The effect of the lower-redshift
interlopers in $\langle T_{\rm IGM}T_{\rm IGM}(\theta)\rangle$ can be
expressed as (see Appendix B)

$\displaystyle\langle T_{\rm IGM}T_{\rm IGM}(\theta)\rangle=(1-f_{\rm
bg.int})^{2}\langle T_{\rm IGM}T_{\rm IGM}(\theta)\rangle^{\rm true}+$
$\displaystyle 2(1-f_{\rm bg.int})f_{\rm bg.int}\langle T_{\rm
IGM}\rangle^{\rm true}\langle T_{\rm IGM}\rangle^{\rm bg.int}+(f_{\rm
bg.int}\langle T_{\rm IGM}\rangle^{\rm bg.int})^{2}.$ (40)

The second and third terms indicate contaminations from cross-correlation
between low-redshift interlopers and true background galaxies and the auto-
correlation of low-redshift interlopers, assuming there is no spatial
correlation. In terms of the Ly$\alpha$ forest angular auto-correlation
function, the true auto-correlation function $\omega^{\rm
true}_{\rm\alpha\alpha}(\theta)=\langle T_{\rm IGM}T_{\rm
IGM}(\theta)\rangle^{\rm true}/(\langle\overline{T}_{\rm IGM}\rangle^{\rm
true})^{2}-1$ is diluted by the interlopers such that

$\displaystyle\omega_{\rm\alpha\alpha}(\theta)=\left[\frac{(1-f_{\rm
bg.int})\langle\overline{T}_{\rm IGM}\rangle^{\rm true}}{(1-f_{\rm
bg.int})\langle\overline{T}_{\rm IGM}\rangle^{\rm true}+f_{\rm
bg.int}\langle\overline{T}_{\rm IGM}\rangle^{\rm
bg.int}}\right]^{2}\omega^{\rm true}_{\rm\alpha\alpha}(\theta).$ (41)

Again, all factors can be determined a posteriori using spectroscopic follow-
up of the background galaxy sample. Assuming $f_{\rm bg.int}=0.20$ and
$\langle\overline{T}_{\rm IGM}\rangle^{\rm bg.int}=0.7$ for the background LAE
sample and taking $\langle\overline{T}_{\rm IGM}\rangle^{\rm true}$ to be the
value from the DEIMOS10k sample can explain the observed offset in $\langle
T_{\rm IGM}T_{\rm IGM}(\theta)\rangle$. This corresponds to the damping of
$\omega_{\alpha\alpha}(\theta)=0.39\,\omega_{\alpha\alpha}^{\rm true}(\theta)$
for the observed Ly$\alpha$ forest angular auto-correlation function from the
LAE sample. There is no damping factor for the DEIMOS10k sample as the
interloper contamination for the spectroscopically confirmed sample is zero.

##  8 IGM tomographic map

Figure 18: (Left): Reconstructed 2D tomographic map of the Ly$\alpha$ forest
transmission of the IGM at $z\simeq 4.9$ (colour map). The locations of the
background ($z\simeq 5.7$) LAEs are shown with coloured circles with the
colours indicating the measured difference between the Ly$\alpha$ forest
transmission along the sightline and the global mean. (Middle): The standard
deviation of the reconstructed map including the uncertainties from the
photometric and continuum errors in the Bayesian SED fitting. The locations of
the background LAE sightlines are shown with open circles. (Right): The
signal-to-noise ratio of the reconstructed map. The masked regions are left
blank.

###  8.1 Reconstruction method

Finally, we present a reconstructed tomographic map of the IGM. This is
arguably the most unique aspect of photometric IGM tomography, since it
enables us to directly visualise the large-scale structures of the IGM and
galaxies in the same cosmic volume. To accomplish this we use the Nadaraya-
Watson estimator for the 2D tomographic map of the IGM Ly$\alpha$ forest
transmission fluctuations (Kakiichi et al., 2022),

$\Delta T_{\rm IGM}^{\rm 2D}(\bm{\theta})=\frac{\sum^{N_{\rm bg}}_{\rm
i=1}K_{R}(\bm{\theta}-\bm{\theta}_{i})(T_{{\rm
IGM},i}-\langle\overline{T}_{\rm IGM}\rangle)}{\sum^{N_{\rm bg}}_{\rm
i=1}K_{R}(\bm{\theta}-\bm{\theta}_{i})},$ (42)

where $K_{R}(\bm{\theta})=(2\pi R^{2})^{-1/2}\exp[-\bm{\theta}^{2}/(2R^{2})]$
is a Gaussian kernel with a smoothing length $R$. In practice, we create a 2D
map on pixelised map of $4096\times 4096$ pixels. Following a similar argument
as in previous sections, given the posterior of Ly$\alpha$ forest
transmission, the expectation value of the 2D tomographic map is given by

$\langle\Delta T_{\rm IGM}^{\rm
2D}(\bm{\theta})\rangle=\frac{\sum_{i=1}^{N_{\rm bg}}\displaystyle
K_{R}(\bm{\theta}-\bm{\theta}_{i})(\langle T_{{\rm
IGM},i}\rangle-\langle\overline{T}_{\rm IGM}\rangle)}{\sum^{N_{\rm bg}}_{\rm
i=1}K_{R}(\bm{\theta}-\bm{\theta}_{i})}.$ (43)

At each point $\bm{\theta}$, the estimator is simply the weighted sum of
(independent) individual $T_{{\rm IGM},i}$ measurements. Thus the variance can
be computed as

$\displaystyle{\rm Var}[\Delta T_{\rm IGM}^{\rm
2D}(\bm{\theta})]=\frac{\sum_{i=1}^{N_{\rm bg}}\displaystyle
K_{R}(\bm{\theta}-\bm{\theta}_{i})^{2}{\rm Var}\left[T_{{\rm
IGM},i}\right]}{\left[\sum^{N_{\rm bg}}_{\rm
i=1}K_{R}(\bm{\theta}-\bm{\theta}_{i})\right]^{2}}.$ (44)

This includes errors from photometric noise and the UV continuum uncertainty.
We define the SNR map as the ratio between the observed Ly$\alpha$ forest
transmission fluctuations and the standard deviation,

${\rm SNR}(\bm{\theta})=\frac{\left|\langle\Delta T^{\rm 2D}_{\rm
IGM}(\bm{\theta})\rangle\right|}{\sqrt{{\rm Var}[\Delta T^{\rm 2D}_{\rm
IGM}(\bm{\theta})]}}.$ (45)

Figure 19: Reconstructed 2D tomographic map of the Ly$\alpha$ forest
transmission of the IGM (coloured map) overlaid with the distribution of the
$z=4.9$ LAEs (coloured circles). The angular resolution of the reconstructed
map is $0.12\rm\,deg$ corresponding to $11\,h^{-1}\rm cMpc$. The colour of
each circle indicates the Ly$\alpha$ luminosity of the LAE. Masked regions are
left blank. Although the map is dominated by photometric noise, it represents
the first IGM tomographic map co-spatial to the large-scale structure of LAEs
in the same cosmic volume close to the end of cosmic reionization.  Figure 20:
Fluctuations of the reconstructed 2D Ly$\alpha$ forest tomographic map (left),
LAE number density field centre), and Ly$\alpha$ luminosity-weighted LAE
number density field (right) with the angular resolution of $0.12\rm\,deg$
($11\,h^{-1}\rm cMpc$) at $z\simeq 4.9$.

###  8.2 Result

In Figure 18 we show the reconstructed 2D tomographic map of the Ly$\alpha$
forest transmission at $z\simeq 4.9$. We only apply the map reconstruction to
the background LAE sample because their distribution spans the entire field of
view. The smoothing length is chosen as the mean inter-sightline separation,
$R=1/\sqrt{\Sigma_{\rm LAE}}\simeq 0.118\rm\,deg\,(7.1\arcmin)$ where
$\Sigma_{\rm LAE}$ is the surface number density of the background LAEs. We
can visually see large-scale fluctuations of the Ly$\alpha$ forest
transmission with a median contrast of $\left|\langle\Delta T^{\rm 2D}_{\rm
IGM}(\bm{x})\rangle\right|\sim 0.1$. The typical standard deviation in the
reconstructed map is $\sqrt{{\rm Var}[\Delta T^{\rm 2D}_{\rm
IGM}(\bm{x})]}\sim 0.14$. We find the mean SNR of the reconstructed map as
$\langle{\rm SNR}(\bm{x})\rangle=0.71$, indicating that our tomographic map is
still noisy with contributions from photometric errors and the UV continuum
uncertainty.

There are several tentative regions of transmissive and opaque transmission in
the map located at ${\rm(RA,DEC)}\simeq(150.2^{\circ},2.45^{\circ})$ and
${\rm(RA,DEC)}\simeq(150.0^{\circ},1.95^{\circ})$ respectively, with the peak
${\rm SNR}\simeq 1.5-2.0$. The Ly$\alpha$ forest transmission in these regions
is reasonably coherent. We require deeper NB imaging data to secure the
statistical significance. If confirmed, these opaque and transmissive regions
of the IGM may represent a protocluster and a highly ionized region of the IGM
by the enhanced UV background.

Although there is a higher transmissive region towards the edge at
${\rm(RA,DEC)}\simeq(150.8^{\circ},1.80^{\circ})$, as there is no background
galaxy here this is likely an artefact from the map reconstruction method. At
the edge of the field of view, the reconstruction method is more affected by
boundary effects. Although our estimator corrects for boundary effects by
incorporating the sightline density, the estimator is more sensitive to the
$T_{\rm IGM}$ values of individual background galaxies, whereas at the centre
of the field smoothing corrects outlier values of $T_{\rm IGM}$.

Figure 19 overlays the distribution of the $z=4.9$ LAEs on the reconstructed
tomographic map of the Ly$\alpha$ forest transmission. This represents the
highest redshift 2D tomographic map of the IGM with the galaxy distribution at
the present time and demonstrates the potential of the NB tomographic
technique to examine the galaxy-IGM connection closer to the reionization
epoch.

We briefly examine the spatial correlation between the $z=4.9$ LAE
distribution and the reconstructed Ly$\alpha$ forest transmission map. In
order to apply the spatial correlation analysis at the map level, we first
reconstruct LAE density map at the same smoothing scale using the Gaussian
kernel density estimator with the boundary correction,

$n_{\rm LAE}(\bm{\theta})=\frac{\sum^{N_{\rm bg}}_{\rm
i=1}w_{i}K_{R}(\bm{\theta}-\bm{\theta}_{i})}{\int
m(\bm{\theta}^{\prime})K_{R}(\bm{\theta}-\bm{\theta}^{\prime})d\bm{\theta}^{\prime}},$
(46)

where $m(\bm{\theta}^{\prime})$ represents the mask and the denominator is the
correction factor $C^{-1}=\int
m(\bm{\theta}^{\prime})K_{R}(\bm{\theta}-\bm{\theta}^{\prime})d\bm{\theta}^{\prime}$
for boundary effects including masked regions around bright stars. We
introduce weights $w_{i}$, where $w_{i}=1$ for the ordinary LAE density field
and $w_{i}=L_{\alpha,i}$ for the Ly$\alpha$ luminosity-weighted LAE density
field. The galaxy density fluctuation map is then $\delta_{\rm
LAE}(\bm{x})=n_{\rm LAE}/\bar{n}_{\rm LAE}-1$ where the mean density
$\bar{n}_{\rm LAE}$ is computed by excluding the masked regions.

In Figure 20 we show the comparison of the 2D tomographic map with the LAE
density and Ly$\alpha$ luminosity-weighted LAE density fields at the same
smoothing length. Using all unmasked regions, the Pearson correlation
coefficient $\rho$ indicates that there is a negligible correlation between
the reconstructed 2D tomographic map and the (luminosity-weighted) LAE density
map with $\rho=0.09$ ($0.06$). Within the precision of existing photometric
data, it appears that $z\simeq 4.9$ LAEs do not occupy extreme Ly$\alpha$
transmissive or opaque regions of the IGM. This is consistent with the two-
point cross-correlation analysis between LAE and Ly$\alpha$ forest
transmission.

To avoid confusion from boundary effects, we have also limited the map-level
analysis within the region with a low correction factor $C<2.0$. We then find
that the Pearson correlation coefficient becomes $\rho=0.22$ for the LAE
density field ($0.20$ for luminosity-weighted field), indicating a possible
weak positive correlation between the LAE number density and the Ly$\alpha$
forest transmission map of the IGM. Although this weak correlation is
intriguing, as noted above, our average SNR of the map is still low to
conclude.

###  8.3 Systematics

Unlike the cross- and auto-correlation functions between LAEs and Ly$\alpha$
forest transmission, the reconstruction of the 2D tomographic map demands a
higher purity of the background galaxy sample. A fictitious Ly$\alpha$ forest
transmission due to a low-redshift interloper would produce a fake
transmissive IGM region. While having many background galaxies within a
smoothing length of the reconstruction reduces the effect of interlopers, we
have not found a way to statistically correct for this effect as is possible
statistically using a spectroscopic subset in the case of the correlation
functions. It is difficult to quantify the level of contamination in each
transmissive or opaque region of the IGM identified with photometric IGM
tomography without directly confirming the redshifts of all background
galaxies spectroscopically.

##  9 Discussion

###  9.1 Improving photometric IGM tomography

####  9.1.1 Extremely-deep NB imaging and spectroscopic campaign

As discussed in Section 4.2, our estimate of $T_{\rm IGM}$ from individual
background galaxies is dominated by photometric noise. This error propagates
into our measurements of LAE-Ly$\alpha$ forest cross-correlation (Section 6)
and Ly$\alpha$ forest auto-correlation functions (Section 7) and the
reconstruction of the 2D tomographic map of the IGM (Section 8). As the
measured Ly$\alpha$ forest transmission depends on the contrast between the
foreground NB flux and the BB flux of a background galaxy, the noise scales
approximately as $\delta T_{\rm IGM}\approx\delta f_{\rm NB}/f_{\rm BB}$. To
improve the signal-to-noise ratio of photometric IGM tomography, we therefore
require (i) deeper imaging in the foreground NB filter and/or (ii) enlarge the
sample of bright (spectroscopically-confirmed) background galaxies for which
we can more accurately measure $T_{\rm IGM}$ given a great contrast between
the NB and BB filters.

While the current NB718 depth ($26.86\rm\,mag,\,3\sigma$) provides $3\sigma$
sensitivity to a typical mean value of the Ly$\alpha$ forest transmission when
using $25.2\,\rm mag$ ($z$-band) background sources, the majority of our
background galaxies are fainter. Secure ($>3\sigma$) detection of the
Ly$\alpha$ forest transmitted flux along individual sightlines is currently
only possible for rare bright background galaxies ($M_{\rm
UV}\lesssim-21.4\rm\,mag$). Extremely-deep NB imaging reaching $27.6\rm\,mag$
($28.2\rm\,mag$) at $3\sigma$ depth would allow us to detect typical
Ly$\alpha$ forest transmitted fluxes to fainter $26.0-26.4\rm\,mag$
($26.6-27.0\rm\,mag$) background sources at $2-3\sigma$ significance level
which comprises $\approx 50\,\%$ ($90\,\%$) of our background LAE+DEIMOS10k
sample. Based on the existing depth of NB718 from CHORUS PDR1 after $t_{\rm
exp}=7.7$ hour exposures (Inoue et al., 2020), and assuming a factor of
$\propto t_{\rm exp}^{-1/2}$ reduction of photometric noise, such an extremely
deep observation would require a total exposure of $\simeq 28$ (100) hours in
NB718. As the reconstructed tomographic map of the IGM is dominated by the
photometric noise, the improvement in the NB image quality by longer
integration will directly increase the signal-to-noise ratio of the IGM
tomographic map. While this may seem a significant investment of the telescope
time, given a large number of potential science applications of photometric
IGM tomography as we will discuss later, such an investment would be of great
interest.

Concerning an increase in the number of bright spectroscopically-confirmed
background sources, while we used a large compilation of spectroscopic
catalogues, previous surveys have focused primarily on the central part of the
COSMOS field, providing only 36 background objects with suitable spectroscopic
redshifts for our IGM tomography. According to the Bouwens et al. (2021) UV
luminosity function, there should be numerous star-forming galaxies brighter
than $m_{\rm UV}<25.5\rm\,mag$ ($M_{\rm UV}\lesssim-21$) in the appropriate
redshift range ($4.98<z<5.89$) with the surface density of $\Sigma_{\rm
g}\approx 726\rm\,deg^{-2}$. This corresponds to a total of $\approx 1280$
sources that can be in principle accessed across the HSC’s $1.76\,\rm deg^{2}$
field. Uncovering this population would provide a large boost in the number of
background galaxies (cf. the surface density of our background LAEs of
$\Sigma_{\rm LAE}\simeq 71.8\rm\,deg^{2}$). The bright UV continua will ensure
$\sim 2-3\sigma$ detection of the Ly$\alpha$ forest transmission with the
current NB718 depth. There are a number of photometric catalogues with dropout
selection (e.g. GOLDRUSH: Harikane et al., 2022) and photometric redshifts
(COSMOS2020: Weaver et al., 2022) in the COSMOS field. We expect at least
$10-20\,\%$ of such UV continuum selected objects will show observable
Ly$\alpha$ emission (Stark et al., 2010; Stark et al., 2011; Mallery et al.,
2012; Cassata et al., 2015; Arrabal Haro et al., 2018; Kusakabe et al., 2020).
A wide-field multi-object spectroscopic (MOS) follow-up campaign in the ultra-
deep HSC footprint of the COSMOS field can locate $\approx 128-256$ background
sources ($\Sigma_{\rm specz}\simeq 72.6-145.2\rm\,deg^{2}$), providing
$2-3\times$ increase in the total background sample (which is currently
dominated by $z\simeq 5.7$ LAEs with typical $\sim 26.4\rm\,mag$ continua). As
discussed in Section 5.3, including a subset of the background LAEs in the
spectroscopic follow-up campaign is also important as it allows us to
statistically correct for the systematic bias in the correlation functions by
low-redshift interlopers. A follow-up spectroscopic survey using wide-field
MOS instruments such as Keck/DEIMOS and upcoming Subaru/Prime Focus
Spectrograph (PFS) and VLT/MOONS is required to improve the significance and
angular resolution of the photometric IGM tomography.

Alternatively, it should also be possible to uncover large numbers of bright
star-forming galaxies with secure redshifts using rest-frame optical emission
line such as H$\beta+[{\rm O\,{\scriptstyle III}}]$ and H$\alpha$ using a
wide-field redshift survey with the NIRCam wide-field slitless spectrograph
(WFSS) on board JWST. Recently Sun et al. (2022b, a) (see also Matthee et al.,
2022a) suggest that $\sim 88\,\%$ of bright $z\sim 6$ star-forming galaxies
show strong $\rm H\beta+[{\rm O\,{\scriptstyle III}}]$ and H$\alpha$ emission
lines detectable with a shallow ($\sim 20\rm\,min$) integration. If this holds
true at $4.98<z<5.89$, the shallow wide-field NIRCam/WFSS survey tiling the
HSC COSMOS field can uncover a factor of $\sim 4-9$ larger number of
background galaxies ($\Sigma_{\rm JWST}\simeq 638\rm\,deg^{-2}$) than a
ground-based wide-field MOS survey targeting Ly$\alpha$ lines.

Figure 21: Comparison of the error budget in the cross-correlation estimated
from the Jackknife covariance matrix (black) and analytic variance (red) for
the LAE (filled circles) and DEIMOS10k (open squares) samples. The analytic
variance includes the error from photometric noise and continuum slope
uncertainty and follows the expected scaling $\sim(\mbox{(typical $T_{\rm
IGM}$ error)}/\langle\overline{T}_{\rm IGM}\rangle)/\sqrt{N_{\rm
pairs}(\theta)}$ as the number of pairs per bin increases.

####  9.1.2 Reducing the errors & systematics in the statistical analysis

For the statistical measurement of the correlation functions, there is patch-
to-patch variance in addition to the photometric noise and UV continuum
uncertainty. In Figure 21 we compare the Jackknife variance with our
propagated errors from photometric noise and UV continuum uncertainty
(Equation 23) from the Bayesian SED fitting framework in the angular mean
Ly$\alpha$ forest transmission around LAEs $\langle T_{\rm
IGM}(\theta)\rangle$.

We find that the photometric noise is the dominant source of uncertainties at
$\theta<3\rm\,arcmin$ ($<5\,h^{-1}\rm cMpc$). The Jackknife variance sometimes
underestimates the error due to the small number of pairs in the inner angular
bins. The photometric noise is comparable for both the LAE and DEIMOS10k
samples because, while the individual error in $T_{\rm IGM}$ is larger in the
LAE sample, the larger sample size reduces the error in the cross-correlation
function.

At larger radii $\theta>3\rm\,arcmin$ ($>5\,h^{-1}\rm cMpc$), the Jackknife
error becomes larger than the propagated error from photometric+continuum
errors, indicating that the patch-to-patch variance in the field becomes the
dominant uncertainty. While the photometric error is the dominant uncertainty
in individual $T_{\rm IGM}$, the error scale as $\propto\delta T_{\rm
IGM}/\sqrt{N_{\rm pair}(\theta)}$ where $\delta T_{\rm IGM}$ is the typical
photometric uncertainty in individual $T_{\rm IGM}$’s. The patch-to-patch
variance may arise from cosmic variance or systematics such as imperfect sky
subtraction or a coherent error in the photometric colours across the field.
We believe that the current dominant source of the patch-to-patch variance is
systematics. We checked this by comparing the Jackknife variances between the
LAE and DEIMOS10k samples. As the former is sampled from the larger area, if
the large-scale patch-to-patch variance is due to cosmic variance, we expect
the Jackknife variance to decrease. However, this is not the case. We have
also compared the Jackknife variance using a smaller number of angular bins to
reduce the photometric error. The photometric error decreases as expected. In
the both cases, the large-scale Jackknife variance remains roughly constant,
suggesting that the observed excess of the Jackknife variance compared to the
photometric+continuum error is likely caused by systematics such as sky
background subtraction or possibly reflected lights within the optical units
during the NB imaging. An unaccounted coherent change in the photometric
colours (e.g. ${\rm NB718}-z$ colour) across the field, e.g. due to imperfect
PSF matching or Galactic dust extinction corrections, could also be a cause of
the systematics. More careful data reduction, background subtraction, and
photometric calibration will be required to reduced the error in the large
angular bins. We emphasise that these issues can be resolved with additional
procedures during data reduction and calibration steps.

If the large-scale Jackknife variance is caused by the cosmic variance, we
would benefit by enlarging the survey area to increase the sampling of the
large-scale modes. Including other pointings such as the SXDS field would
increase the constraining power on the large-scale correlation functions.
Whether such an investment is worthwhile depends on the theoretically expected
scale of the LAE-Ly$\alpha$ forest cross-correlation and Ly$\alpha$ forest
auto-correlation functions (see companion paper, Kakiichi et al in prep).
Further theoretical studies on photometric IGM tomography are necessary in
order to understand the physical information contained on the different scales
of the correlation functions.

####  9.1.3 Combining multi-wavelength data

Our present analysis employed only two broad-band filters (HSC $z$ and $y$) to
constrain the SEDs of background galaxies. The SED uncertainties (i.e
continuum slope and SED template) can be mitigated by including multi-
wavelength datasets available in the COSMOS field. Inclusion of near-infrared
data such as UltraVISTA $JHK$ (McCracken et al., 2012) will provide a longer
baseline to better characterise the intrinsic SEDs. Constraining the dust and
age of the background galaxies would eliminate the $\sim 6-27\,\%$ error in
the measured $T_{\rm IGM}$ from SED fitting (Section 4.2). The soon-available
NIRCam imaging with F115W, F150W, F277W, F444W filters from COSMOS-Web
(Kartaltepe et al. 2021, ID 1727, PI:, see also Casey et al. 2022) can
precisely determine the rest UV-to-optical SEDs of background galaxies.

###  9.2 Science applications

####  9.2.1 LyC escape fraction and the nature of ionizing sources

The statistical analysis of the LAE-Ly$\alpha$ forest angular cross-
correlation from photometric IGM tomography provides a measure of the gas
overdensity, temperature, and UV background fluctuations of the IGM around
foreground LAEs. The angular cross-correlation is the NB-averaged version of
the underlying galaxy-Ly$\alpha$ forest 3D cross-correlation that can be
measured from the spectroscopic galaxy survey in quasar fields (Kakiichi et
al., 2018; Meyer et al., 2020). We can thus adopt a similar approach to
constrain the population-averaged LyC escape fraction from the LAE-Ly$\alpha$
forest cross-correlation. We will present the analysis in a companion paper
(Kakiichi et al in prep).

The observed limit on the LAE-Ly$\alpha$ forest cross-correlation could also
limit the contribution of bright LAEs to the UV background. Matthee et al.
(2022b); Naidu et al. (2022) claimed that bright LAEs with $L_{\rm
Ly\alpha}\gtrsim 10^{42}\rm\,erg\,s^{-1}$ could contribute significantly to
the total ionizing budget at $z\gtrsim 4$. Such bright LAEs may produce a
large proximity zone which could be observed by photometric IGM tomography.
The measurement of LAE-Ly$\alpha$ forest cross-correlation can be used to test
this scenario. We will discuss the implication of our result on on the
relative contribution of bright and faint galaxies to the total ionizing
budget in a following paper (Kakiichi et al in prep).

####  9.2.2 Search for ionized bubbles and protoclusters

Besides the statistical analysis, the IGM tomography can also be used to
search for ionized bubbles and/or protoclusters at high redshifts. While our
present analysis focused on NB718 filter ($z\simeq 4.9$), one can perform
tomographic analyses at any redshift where NB filters are available, including
NB527 ($z\simeq 3.31$) and NB816 ($z\simeq 5.72$) (see Kakiichi et al., 2022,
for the full list), assuming adequate deep NB imaging and background galaxies
are available.

Ionized bubbles or the high UV background regions are expected to show
transmissive Ly$\alpha$ forest with effective optical depth of $\tau_{\rm
eff}\approx 2$ at $z\simeq 5.7$ (Davies et al., 2018; Keating et al., 2020).
This corresponds to a NB-BB magnitude decrement of $m_{\rm NB816}-m_{\rm
UV}=-2.5\log_{10}e^{-\tau_{\rm eff}}\approx 2.2\rm\,mag$. With a NB816 depth
of 27.7 mag, it is possible to search for such regions using 25.5 mag
background sources at $5.81<z<6.86$.

The present $2-3\sigma$ depth of $26.89-27.33\rm\,mag$ of the NB816 image
already allows us to search for extreme transmissive IGM regions with
$\tau_{\rm eff}\sim 1-2$ ($m_{\rm NB816}-m_{\rm UV}\approx 1.1-2.2\rm\,mag$).
One such region ($\tau_{\rm eff}\simeq 1.2$) has been already discovered
serendipitously by Bosman et al. (2020) around a quasar using NB816 imaging.
While current cosmological simulations (Davies et al., 2018; Keating et al.,
2020) do not predict such extreme values, they are performed under the
assumption that reionization is driven by a large population of relatively
faint galaxies. Highly transmissive regions may exist if the contribution from
quasars (Chardin et al., 2017) or luminous galaxies with high LyC escape
fraction and ionizing photon production efficiency (e.g. Endsley et al., 2021;
Topping et al., 2022; Marques-Chaves et al., 2022) are important.

Similarly, one can conduct a search for protoclusters as coherently opaque
regions of the IGM in the tomographic map (Lee et al., 2016; Newman et al.,
2020). Indeed, Mawatari et al. (2017) used the photometric IGM tomographic
technique to examine the protocluster region with NB497 ($z\simeq 3.1$). The
$10-40\,h^{-1}\rm cMpc$ scales of coherently strong Ly$\alpha$ absorption (Cai
et al., 2016) is shown to be associated with an overdensity at $z\sim 2-3$
(Cai et al., 2017; Shi et al., 2021) whose scale matches closely with the
line-of-sight width of a NB filter. LAEs in the tomographic slice permit an
immediately confirmation of whether the coherently opaque IGM region is
associated with a galaxy overdensity without a separate spectroscopic or
dedicated imaging campaign. This provides an opportunity to test the relation
between galaxy overdensities and the IGM environments (e.g. Momose et al.,
2021; Newman et al., 2022) at higher redshifts without a need of extremely-
deep spectroscopy of background galaxies.

####  9.2.3 Quasar light-echoes and past AGN activity

The relatively short variability timescale of a quasar $t_{\rm Q}\sim
10^{6-8}\rm\,yr$ compared to the light crossing time of the IGM tomographic
map means that the radiation from the quasar could leave an imprint on the
ionization state of the IGM well after the non-thermal activity ceases
(Adelberger, 2004). Schmidt et al. (2019); Kakiichi et al. (2022) examined the
prospect of examining the lifetime/past AGN activity of an active quasar using
this light echo signal. Our photometric IGM tomography demonstrates that
searching for quasar light echoes is possible if an appropriate quasar field
is targeted.

In principle, the search for quasar light-echoes is not limited to the region
around a luminous quasar, but can be applied to any region around a massive
galaxy where quasar activity might have occurred during its previous $\sim
10^{8}\rm\,yr$. Searches for fossil light-echoes around LAEs, or LBGs at
redshift within the region of influence of the NB tomographic slice, could be
used to constrain the past luminous ionizing activity of a source as a
function of the travel time between the source and a position of the IGM.
Bosman et al. (2020) used the photometric detection of the Ly$\alpha$ forest
transmission in the NB816 filter located slightly foreground of $z\simeq 5.8$
quasar to show that the quasar was active for at least $\sim 2\times
10^{4}\rm\,yr$ in the past. Bowler et al. (2015, 2020); Ono et al. (2018);
Harikane et al. (2022) suggest that the double-power law luminosity function
at $z\gtrsim 4$ may be a sign of inefficient quenching by AGN feedback acting
on luminous galaxies. Finding a lack of highly transmissive regions around
luminous galaxies would imply that these systems could not release a
significant ionizing radiation by quasar activity during their history.

####  9.2.4 Correlation with direct LyC, ${\rm He\,{\scriptstyle II}}$, and
${\rm C\,{\scriptstyle IV}}$ emitters

Multiple HSC NB imaging available in the COSMOS field provides a NB
photometric measure of LyC leakage and the identification of strong ${\rm
He\,{\scriptstyle II}}$ and ${\rm C\,{\scriptstyle IV}}$ emission lines for
$z\simeq 4.9$ LAEs via the CHORUS survey (Inoue et al., 2020). The correlation
of these populations with photometric IGM tomographic map will have a number
of applications.

For example, direct search for LyC emission along the line-of-sight of
galaxies is increasingly difficult at $z\gtrsim 3.5$ as the IGM becomes
increasingly opaque on average. Bassett et al. (2021, 2022) showed that any
LyC detection could be biased towards the rare transmissive regions of the IGM
and the bias introduced by assuming an average IGM transmission is severe at
$3<z<5$. While LyC leaking candidates have reported at this redshift range
(e.g. Shapley et al., 2016; Vanzella et al., 2018; Ji et al., 2020; Meštrić et
al., 2020; Prichard et al., 2022; Rivera-Thorsen et al., 2022; Marques-Chaves
et al., 2022), the uncertain IGM transmission value against LyC photons make
it difficult to convert the observed LyC flux to an absolute LyC escape
fraction and sometimes results in unphysical values $f_{\rm esc}>1$. Fletcher
et al. (2019) also noted the possibility of spatial variations in the inferred
LyC escape fraction due to the uncertain IGM transmission. Photometric IGM
tomography provides an useful independent measure of the IGM transmission and
could be used to alleviate the issue of uncertain IGM LyC opacities along the
lines-of-sight to high-redshift galaxies.

Combining the direct detection of LyC leakage and the indirect measurement
from galaxy-Ly$\alpha$ forest cross-correlation, one can test the relative
contributions of luminous and faint galaxies to the ionising budget. The
former measures the ionising contribution of galaxies above the detection
limit, while the latter provides a population-averaged LyC escape fraction
including those below the detection limit. NB measurements of ${\rm
C\,{\scriptstyle IV}}$ and ${\rm He\,{\scriptstyle II}}$ from LAEs are also
valuable for examining hard ionising sources like AGN and X-ray binaries. The
IGM tomographic map enables us to connect the properties of these populations
with their large-scale IGM environments.

##  10 Conclusions

We present a novel technique called photometric IGM tomography to map the
large-scale structure of the IGM at $z\sim 5$ in the COSMOS field. The
technique utilizes ultra-deep NB718 imaging to detect the Ly$\alpha$ forest
transmission along sightlines to various background galaxies including
spectroscopically-confirmed DEIMOS10k sources and NB816-selected LAE catalogue
from SILVERRUSH. In this paper, we describe a science verification of this new
technique using public HSC data including HSC-SSP DR3 and CHORUS PDR1.

We have developed a Bayesian SED fitting framework to measure the Ly$\alpha$
forest transmission along background galaxies. This allows us to accurately
propagate the error from photometric noise and uncertainty from the assumed
SED template into the final measurement of the Ly$\alpha$ forest transmission.
At the current imaging depths, photometric noise from NB718 imaging dominates
the total error in the estimated Ly$\alpha$ forest transmission. Uncertainties
from the UV continuum slopes and SED templates are subdominant in the present
analysis.

Using a total sample of 140 background sources, we have photometrically
measured the NB718-integrated mean Ly$\alpha$ forest transmission at $z\simeq
4.9$. We find that our result is consistent with the previous measurement