#  Cosmic Himalayas: The Highest Quasar Density Peak Identified in a 10,000
deg2 Sky
with Spatial Discrepancies between Galaxies, Quasars, and IGM Hi

Yongming Liang Institute for Cosmic Ray Research, The University of Tokyo,
5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8582, Japan  Masami Ouchi Institute for
Cosmic Ray Research, The University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa,
Chiba 277-8582, Japan  National Astronomical Observatory of Japan, 2-21-1
Osawa, Mitaka, Tokyo 181-8588, Japan  Dongsheng Sun Institute for Cosmic Ray
Research, The University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8582,
Japan  Nobunari Kashikawa Department of Astronomy, School of Science, The
University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo, 113-0033, Japan  Research
Center for the Early Universe, The University of Tokyo, 7-3-1 Hongo, Bunkyo-
ku, Tokyo 113-0033, Japan  Zheng Cai Department of Astronomy, Tsinghua
University, Beijing 100084, China  Sebastiano Cantalupo Department of Physics,
University of Milan Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy
Kentaro Nagamine Department of Earth and Space Science, Osaka University, 1-1
Machikaneyama, Toyonaka, Osaka 560-0043, Japan  Kavli IPMU (WPI), The
University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8583, Japan
Department of Physics and Astronomy, University of Nevada, Las Vegas, 4505 S.
Maryland Pkwy, Las Vegas, NV 89154-4002, USA  Hidenobu Yajima Center for
Computational Sciences, University of Tsukuba, 1-1-1 Tennodai, Tsukuba,
Ibaraki 305-8577, Japan  Takanobu Kirihara Kitami Institute of Technology, 165
Koen-cho, Kitami, Hokkaido 090-8507, Japan  Center for Computational Sciences,
University of Tsukuba, 1-1-1 Tennodai, Tsukuba, Ibaraki 305-8577, Japan
Haibin Zhang National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka,
Tokyo 181-8588, Japan  Mingyu Li Department of Astronomy, Tsinghua University,
Beijing 100084, China  Rhythm Shimakawa Waseda Institute for Advanced Study
(WIAS), Waseda University, 1-21-1 Nishi-Waseda, Shinjuku, Tokyo 169-0051,
Japan  Center for Data Science, Waseda University, 1-6-1 Nishi-Waseda,
Shinjuku, Tokyo 169-0051, Japan  Xiaohui Fan Steward Observatory, University
of Arizona, 933 N. Cherry Avenue, Tucson, AZ 85721, USA  Kei Ito Department of
Astronomy, School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku,
Tokyo, 113-0033, Japan  Masayuki Tanaka Department of Astronomical Science,
SOKENDAI (The Graduate University for Advanced Studies), Mitaka, Tokyo
181-8588, Japan  National Astronomical Observatory of Japan, 2-21-1 Osawa,
Mitaka, Tokyo 181-8588, Japan  Yuichi Harikane Institute for Cosmic Ray
Research, The University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8582,
Japan  J. Xavier Prochaska UCO/Lick Observatory, University of California,
1156 High Street, Santa Cruz, CA 95064, USA  Andrea Travascio Department of
Physics, University of Milan Bicocca, Piazza della Scienza 3, I-20126 Milano,
Italy  Weichen Wang Department of Physics, University of Milan Bicocca, Piazza
della Scienza 3, I-20126 Milano, Italy  Martin Elvis Harvard-Smithsonian
Center for Astrophysics, 60 Garden St., Cambridge, MA 02138, USA  Giuseppina
Fabbiano Harvard-Smithsonian Center for Astrophysics, 60 Garden St.,
Cambridge, MA 02138, USA  Junya Arita Department of Astronomy, School of
Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo, 113-0033,
Japan  Masafusa Onoue Kavli IPMU (WPI), The University of Tokyo, 5-1-5
Kashiwanoha, Kashiwa, Chiba 277-8583, Japan  John D. Silverman Kavli IPMU
(WPI), The University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8583,
Japan  Dong Dong Shi Purple Mountain Observatory, Chinese Academy of Sciences,
10 Yuan Hua Road, Nanjing, Jiangsu 210023, People’s Republic of China  Fangxia
An Purple Mountain Observatory, Chinese Academy of Sciences, 10 Yuan Hua Road,
Nanjing, Jiangsu 210023, People’s Republic of China  Takuma Izumi National
Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan
Kazuhiro Shimasaku Department of Astronomy, School of Science, The University
of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo, 113-0033, Japan  Research Center for
the Early Universe, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo
113-0033, Japan  Hisakazu Uchiyama National Astronomical Observatory of Japan,
2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan  Chenghao Zhu Institute for Cosmic
Ray Research, The University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba
277-8582, Japan

###### Abstract

We report the identification of a quasar overdensity in the BOSSJ0210 field,
dubbed Cosmic Himalayas, consisting of 11 quasars at $z=2.16-2.20$, the
densest overdensity of quasars ($17\sigma$) in the $\sim$10,000 deg2 of the
Sloan Digital Sky Survey. We present the spatial distributions of galaxies and
quasars and an Hi absorption map of the intergalactic medium (IGM). On the map
of 465 galaxies selected from the MAMMOTH-Subaru survey, we find two galaxy
density peaks that do not fall on the quasar overdensity but instead exist at
the northwest and southeast sides, approximately 25 $h^{-1}$ comoving-Mpc
apart from the quasar overdensity. With a spatial resolution of 15 $h^{-1}$
comoving Mpc in projection, we produce a three-dimensional Hi tomography map
by the IGM Ly$\alpha$ forest in the spectra of 23 SDSS/eBOSS quasars behind
the quasar overdensity. Surprisingly, the quasar overdensity coincides with
neither an absorption peak nor a transmission peak of IGM Hi but lies near the
border separating opaque and transparent volumes, with the more luminous
quasars located in an environment with lesser IGM Hi. Hence remarkably, the
overdensity region traced by the 11 quasars, albeit all in coherently active
states, has no clear coincidence with peaks of galaxies or Hi absorption
densities. Current physical scenarios with mixtures of Hi overdensities and
quasar photoionization cannot fully interpret the emergence of Cosmic
Himalayas, suggesting this peculiar structure is an excellent laboratory to
unveil the interplay between galaxies, quasars, and the IGM.

Quasars (1319) — Large-scale structure of the universe (902) — Intergalactic
medium (813) — High-redshift galaxies (734) — Protoclusters (1297)

##  1 Introduction

Quasars, or quasi-stellar objects (QSOs), represent some of the most luminous
and extreme entities in the cosmos, characterized by their stellar-like
appearance. This appearance is attributed to the intense luminosity emanating
from the active galactic nucleus (AGN), where material actively accretes onto
a supermassive black hole (SMBH) at the galactic center (e.g., Hopkins et al.,
2006). Such extreme systems are pivotal to the fields of astronomy and
cosmology, offering insights into the growth mechanisms of SMBHs in the early
universe and aiding in our understanding of AGN feedback to galaxy formation
and evolution (Fabian, 2012).

Quasars are recognized as remarkable tracers of the universe’s large-scale
structure (LSS). Typically consisting of at least five quasars and spanning
sizes greater than $100$ $h^{-1}$ comoving-Mpc (cMpc), Large quasar groups
(LQGs), such as the one first discovered by Webster (1982), have predominantly
been identified at relatively low redshifts ($z<2$) through both initial
surveys and later, more comprehensive redshift surveys like 2dF and SDSS
(Clowes, 1986; Crampton et al., 1989; Clowes & Campusano, 1991; Komberg et
al., 1996; Clowes et al., 2012, 2013; Park et al., 2015). The presence of
these structures poses questions regarding their alignment with predictions of
standard LSS formation theories (Clowes et al., 2013; Pilipenko & Malinovsky,
2013; Nadathur, 2013). Moreover, a consensus suggests the need for a smaller
linking length, under $50$ $h^{-1}$cMpc, to address the impact of random
statistics within the SDSS quasar survey (Park et al., 2015), with similar
structures at higher redshifts being seldom identified.

Beyond redshift $z=2$, in the period known as Cosmic Noon, environments of
high density are considered precursors to contemporary galaxy clusters, i.e.,
protoclusters (Cen & Ostriker, 2000; Overzier, 2016). Luminous quasars,
including those that are radio-loud, have been identified as ideal tracers for
protoclusters, facilitating their detection and examination (Wylezalek et al.,
2013; Onoue et al., 2018), while some opposite suggestions also persist
(Kashikawa et al., 2007; Uchiyama et al., 2018, 2020). However, the role of
fainter AGNs and their influence on the clustering properties of galaxies in
these early cosmic structures was initially overlooked. Recent studies now
indicate that AGNs, irrespective of their luminosity, significantly contribute
to the formation and governance of galaxy populations within these dense early
universe environments (Ito et al., 2022; Shimakawa et al., 2023).

A promising method to study the distribution of matter at these high redshifts
is three-dimensional (3D) intergalactic medium (IGM) tomography. By analyzing
the Ly$\alpha$ absorption in the spectra of background quasars, researchers
can map the three-dimensional distribution of hydrogen gas, tracing the
filaments and voids that constitute the LSSs (Lee et al., 2014, 2016, 2018;
Newman et al., 2020; Horowitz et al., 2022; Sun et al., 2023). This technique
offers an unparalleled view of the interplay between quasars, galaxy
formation, and the IGM, providing insights into the underlying processes that
have feedback to LSS (Newman et al., 2022; Dong et al., 2023).

In this study, we report the discovery of Cosmic Himalayas, a significant
quasar overdensity within the J0210 field at $z\approx 2.16-2.20$. Eleven
quasars are identified within a volume of (40 cMpc)3 from SDSS-IV/eBOSS.
Previous observations utilized narrowband and $g$-band imaging from the Subaru
Telescope’s Hyper Suprime-Cam (HSC; Miyazaki et al., 2018) to study the
galaxy–IGM correlation in this field (Liang et al., 2021, hereafter L21). They
target the strong Ly$\alpha$ absorption in IGM employing Sloan Digital Sky
Survey (SDSS)/eBOSS spectra for IGM absorbers to identify massive structures
within the Mapping the Most Massive Overdensity Through Hydrogen with Subaru
project (MAMMOTH-Subaru; Cai et al., 2016, 2017; Liang et al., 2021; Ma et
al., 2024; Zhang et al., 2024, M. Li et al. in prep., Y. Liang et al. in
prep.).

This paper is organized as follows: Section 2 introduces the data and
observations utilized in this study. In Section 3, we analyze the spatial
distribution of SDSS quasars and Ly$\alpha$ emitters (LAEs), highlighting the
identification of the Cosmic Himalayas and its significance within the LSS.
Section 4 presents our methodology and results from the 3D IGM Hi tomography.
Finally, Section 5 discusses the implications of our findings on the
association between the clustered quasars and IGM hydrogen, the possible
scenarios resolving a peculiar ionized structure, and the potential triggers
of Cosmic Himalayas. We conclude with a summary of our key results in the
final section. Throughout the paper, we employ cosmology with
$\Omega_{m}=0.3$, $\Omega_{\Lambda}=0.7$, and $h=0.7$. AB magnitudes are used
if there is no specific mention.

##  2 Data

###  2.1 SDSS Quasars

ID | RA (J2000) | DEC (J2000) | Redshift | Plate | MJD | Fiber |  $i$-Mag (AB) |  $L_{\rm bol}/10^{46}$ erg s-1 | $\delta_{\rm F}$
---|---|---|---|---|---|---|---|---|---
$[1]$ | $[2]$ | $[3]$ | $[4]$ | $[5]$ | $[6]$ | $[7]$ | $[8]$ | $[9]$ | $[10]$
CH-Q01 | 02:09:53.16 | +00:55:11.0 | 2.1891$\pm$0.0002 | 9383 | 58097 | 585 | 18.89$\pm$0.01 | 6.579$\pm$3.289 | 0.131
CH-Q02 | 02:09:53.48 | +00:51:13.5 | 2.1907$\pm$0.0007 | 702 | 52178 | 535 | 19.58$\pm$0.02 | 2.742$\pm$1.371 | 0.160
CH-Q03 | 02:09:21.92 | +00:51:48.9 | 2.1669$\pm$0.0005 | 7833 | 57286 | 122 | 21.01$\pm$0.05 | 2.022$\pm$1.011 | 0.186
CH-Q04 | 02:10:09.69 | +00:39:51.5 | 2.1691$\pm$0.0004 | 7834 | 56979 | 635 | 21.42$\pm$0.10 | 1.008$\pm$0.504 | 0.084
CH-Q05 | 02:10:00.50 | +00:47:27.0 | 2.1748$\pm$0.0003 | 7833 | 57286 | 126 | 20.47$\pm$0.04 | 0.950$\pm$0.475 | 0.123
CH-Q06 | 02:09:59.37 | +01:06:48.7 | 2.1962$\pm$0.0003 | 7833 | 57286 | 134 | 21.10$\pm$0.06 | 0.825$\pm$0.413 | 0.119
CH-Q07 | 02:09:52.55 | +00:43:08.0 | 2.1715$\pm$0.0004 | 7834 | 56979 | 624 | 21.33$\pm$0.08 | 0.791$\pm$0.396 | 0.184
CH-Q08 | 02:10:24.75 | +00:53:57.7 | 2.1693$\pm$0.0009 | 4235 | 55451 | 834 | 21.78$\pm$0.13 | 0.725$\pm$0.363 | 0.081
CH-Q09 | 02:10:22.06 | +00:57:45.0 | 2.2000$\pm$0.0006 | 4235 | 55451 | 832 | 21.69$\pm$0.10 | 0.501$\pm$0.251 | -0.096
CH-Q10 | 02:10:13.03 | +00:56:13.5 | 2.1997$\pm$0.0005 | 7834 | 56979 | 621 | 21.64$\pm$0.11 | 0.457$\pm$0.229 | 0.019
CH-Q11 | 02:10:14.29 | +00:53:13.5 | 2.1921$\pm$0.0008 | 7834 | 56979 | 623 | 21.32$\pm$0.08 | 0.291$\pm$0.146 | 0.017
Table 1:  SDSS-IV/eBOSS Quasar Information for the Cosmic Himalayas: $[1]$
Quasar ID; $[2]$ Right Ascension (RA), equinox J2000; $[3]$ Declination (DEC),
equinox J2000; $[4]$ Spectroscopic redshift from both PCA fitting in the SDSS-
IV/eBOSS pipeline and visual confirmation; $[5]$ SDSS plate ID; $[6]$ SDSS
Modified Julian Date (MJD); $[7]$ SDSS fiber ID; $[8]$ SDSS $i$-band
magnitude; $[9]$ Bolometric luminosity ($L_{\rm bol}$) derived from
$L_{1350}$; $[10]$ Ly$\alpha$ transmission fluctuation ($\delta_{F}$) at the
quasar’s location.

Figure 1:  Left (a): Sky distribution in part of SDSS quasars, with whiter
shades indicating higher quasar densities. Right (b): Number counts of SDSS
quasars within a $\left(40~{}{\rm cMpc}\right)^{3}$ volume across the entire
10,000 deg2 survey area with a step size of 20 cMpc along each dimension. The
average count is 0.36, with J0210 hosting 11 quasars, marking a density nearly
30 times the mean. For context, other fields with high SDSS quasar densities
are highlighted (Cai et al., 2017; Shi et al., 2021); even these are less than
half of the overdensity found for J0210.  Figure 2: Bolometric luminosity
($L_{\rm bol}$) of proximate quasars in the Cosmic Himalayas. All quasars
exhibit $L_{\rm bol}\gtrsim 10^{45.5}~{}{\rm erg}~{}{\rm s}^{-1}$. The gray
shaded area suggests the central 68% scatter of $L_{\rm bol}$ in SDSS DR16
quasars spectroscopically determined by Wu & Shen (2022).

The field BOSS J0210+0052, also called J0210, is selected from the $10,000$
deg2 SDSS/eBOSS database. It is simultaneously traced by a group of line-of-
sight hinting at the intergalactic Ly$\alpha$ absorbers and by a significant
quasar overdensity (L21). The quasars are identified as type-1 AGNs with
strong broad emission lines, based on spectra from the Baryon Oscillation
Spectroscopic Survey (BOSS) of SDSS-III (Dawson et al., 2013) and the
extended-BOSS (eBOSS) of SDSS-IV (Dawson et al., 2016). Moreover, spectra of
background quasars from SDSS/eBOSS were utilized in L21 for selecting
candidate fields for Subaru observations and analyzing Ly$\alpha$ absorption
in IGM Hi for galaxy–-IGM correlation studies. The SDSS/eBOSS survey
encompasses over 200,000 spectra of quasars down to $r=22$ across a survey
area exceeding 10,000 deg2, translating to a survey volume greater than 1 Gpc3
at $z\sim 2$ (Dawson et al., 2016). In focusing on the field J0210, this study
examines both the proximate quasars within the redshift range $2.15\leq z\leq
2.20$ and the background quasars at $z>2.20$ to construct 3D IGM Hi tomography
maps, as detailed in Section 4.

###  2.2 Subaru/Hyper Suprime-Cam Images

The J0210 field was observed in L21 as a part of the MAMMOTH-Subaru project.
Observations were conducted to identify LAEs using HSC, a high-performance
camera with a wide field of view (FoV; $\approx 1.5$ degrees in diameter),
mounted at the prime focus of the 8.2-meter Subaru telescope atop Mauna Kea,
Hawaii (Miyazaki et al., 2012, 2018). Broadband $g$-band and narrowband (NB)
NB387 images ($\lambda_{0}$ = $3,862$ Å, FWHM = 56 Å) were taken by Subaru/HSC
during queue observations between January 11 and January 20, 2018. The total
on-source exposure times were 165 minutes for NB387 and 40 minutes for the
$g$-band, divided into eleven 900-second and four 600-second exposures,
respectively. Frames observed under poor conditions, such as when the seeing
exceeded FWHM $>1.\arcsec 3$ or tracking failed due to low transparency, were
excluded.

Data reduction is detailed in L21. In summary, Subaru/HSC data were reduced
with hscPipe versions 5.4 and 6.6 (Bosch et al., 2018; Aihara et al., 2019).
Zero-point calibration for the $g$-band was performed using standard methods,
while for NB387, calibration involved realistic stars from SDSS to mitigate
potential biases from stellar metallicity in stellar templates (e.g., Pickles,
1998). As a result, the images achieved $5\sigma$ limiting magnitudes of NB387
= 24.36 and $g$ = 26.34. The FWHMs of the point spread function (PSF) are
1.$\arcsec 2$ and 0.$\arcsec 9$ for NB387 and $g$-band. During source
extraction with SExtractor (Bertin & Arnouts, 1996), regions saturated by the
brightest stars or exhibiting high root-mean-square (RMS) noise due to CCD
mosaicking failures were masked, shown as the gray shaded regions in the
following sky maps. The stacked NB387 image enables the detection of
Ly$\alpha$ emission at redshifts of $z=2.177\pm 0.023$, with the $g$-band
image assessing the continuum levels of objects. In total, 465 LAEs were
selected by color excess in NB387, indicative of a rest-frame equivalent width
for Ly$\alpha$ emission EW${}_{0}\geq$ 20 Å.

This field is covered in the HSC Subaru Strategic Program’s wide-field layers
(SSP; Aihara et al., 2018, 2019, 2022), offering access to multiple broadband
optical data sets with limits of $r<26.1$, $i<25.9$, $z<25.1$, and $Y<24.4$.
Alongside photometry from SDSS, we utilized the HSC SSP third data release
(DR3; Aihara et al., 2022) to analyze the UV slopes and bolometric
luminosities of quasars, as detailed in Section 3.1.

##  3 Identifying The Cosmic Himalayas

###  3.1 Quasar Overdensity

Figure 3: Overdensity map of LAEs, with black circles marking LAE positions,
sized according to UV luminosity based on $g$-band. The background contour
represents LAE overdensity, $\delta_{\rm LAE}$, with bluer areas indicating
higher densities. Yellow crosses show SDSS/eBOSS quasars at $z=2.15-2.20$.
Black squares and crosses highlight identified Ly$\alpha$ blobs (LABs) and the
LAEs with NUV counterparts from GALEX MIS survey (NUV-LAEs). The black pluses
mark the median centroids of member LAEs in Node-A (NW) and Node-B (SE),
respectively, which are the two overdensities traced by LAEs. The gray shaded
regions indicate the masks due to saturation by bright stars or high RMS noise
with CCD mosaicking failures.

We report on a unique cosmological structure in the J0210 field, identified
through SDSS quasars and dubbed the Cosmic Himalayas. Eleven proximate quasars
are situated in a volume of approximately $(40{\rm cMpc})^{3}$ within J0210.
We visually inspect all their spectra and confirm that they are indeed type-1
QSO with broad emission. These quasars represent a group of luminous AGNs,
with M${}_{\rm UV}\lesssim-23$ at $z=2.167-2.200$, derived from their SDSS
$u$-band magnitudes. Their redshifts are spectroscopically estimated by the
SDSS pipeline based on PCA fitting, and we have also visually confirmed their
reliability. Their spectroscopic redshifts and $i$-band magnitudes, a primary
reference for SDSS/eBOSS target selection (Dawson et al., 2016), are listed in
Table 1.

The quasar concentration in the Cosmic Himalayas represents the highest
density peak of SDSS quasars across the SDSS/eBOSS survey area of 10,000 deg2
sky in the redshift range $z=2.15-2.20$, as illustrated in Figure 1 (a). The
distribution of SDSS quasar number counts is depicted in Figure 1 (b), with an
average quasar number of 0.36 per $(40{\rm cMpc})^{3}$ volume across the
survey area and a standard deviation of 0.63 from a Gaussian fit to the
distribution. Thus, the quasar overdensity, which is defined as the number
density excess compared to the mean over the SDSS survey area $\delta_{\rm
QSO}=N_{\rm QSO}/\left<N_{\rm QSO}\right>-1$, is calculated to be $\delta_{\rm
QSO}\approx 30$. This $\delta_{\rm QSO}$ reaches a significance of
$16.9\sigma$, marking it as the most significant AGN overdensity discovered to
date at $z>2$ on a similar or larger scales. This concentration exceeds by a
factor of two to three the previously identified extreme cases, such as
BOSS1244 or BOSS1441, which contained 4 or 5 quasars at a similar rest-UV
detection limit of $M_{\rm UV}=-23$ within a comparable volume (Cai et al.,
2017; Shi et al., 2021).

For reference, we estimate the bolometric luminosities $L_{\rm bol}$ of the
eleven quasars. Utilizing the SDSS $ugriz$ and HSC SSP $grizY$ bands, their UV
continuum slopes $\beta_{\rm QSO}$ are estimated by fitting a power-law
continuum. The slope $\beta_{\rm QSO}$ aids in predicting the monochromatic
continuum luminosity at rest-frame 1350 Å, denoted as $L_{1350}$. With the
correction factor from Shen et al. (2011), we calculate the bolometric
luminosity as follows:

$L_{\rm bol}=3.81\times L_{1350}.$ (1)

These results, detailed in Table 1 and depicted in Figure 2, show that all
eleven quasars exhibit $L_{\rm bol}\gtrsim 10^{45.5}$ erg s-1, indicating
their status as part of a relatively luminous population among DR16 quasars
(Wu & Shen, 2022). We label the eleven SDSS quasars with IDs ranked by their
$L_{\rm bol}$, and the quasars from CH-Q01 to CH-Q11 are those from the most
to the least. We note that the measurements based on SDSS and HSC photometry
are consistent with previous literature using SDSS spectra within 1% for
CH-Q01 and CH-Q02; within 0.3 dex for all quasars, including the faintest ones
(Wu & Shen, 2022).

###  3.2 LAE Overdensities and Properties

In L21, Subaru/HSC NB387 observations were performed on J0210 to select the
coeval LAEs – a young galaxy population with strong Ly$\alpha$ emission, at
redshift $z=2.177\pm 0.023$. While covering the quasars in redshifts, these
LAEs also well overlap the SDSS quasars in the Cosmic Himalayas in the
projected sky, as shown in Figure 3. Notably, the LAE candidates construct a
cosmological filament with an end-to-end size of over $100$ cMpc. This
cosmological filament intersects the quasar overdensity and links two cosmic
nodes, i.e., the overdensity of LAEs in the projected sky. The northwestern
(NW) node is labeled as Node-A, and the southeastern (SE) node is labeled as
Node-B. The centroids of the two nodes, iteratively derived from the median
coordinates of member galaxies as defined below, are labeled as points A and
B, respectively.

The LAE overdensities at Node-B and Node-A exhibit significances of
approximately $4\sigma$ and $6\sigma$, respectively, with peak overdensities
of $\delta_{\rm LAE}\approx 1.5$ and $\approx 2.5$ on a 20 cMpc scale. If we
assume each node will eventually collapse to one single virialized structure
at $z=0$, their predicted halo masses at $z=0$ are M${}_{h}(z=0)\approx
10^{14.3}~{}{\rm M}_{\odot}$ for Node-A and M${}_{h}(z=0)\approx
10^{14.8}~{}{\rm M}_{\odot}$ for Node-B, respectively, according to
simulations based on galaxy overdensity on a comparable 25 cMpc scale (Chiang
et al., 2013). This suggests both nodes are likely progenitors of Virgo-type
galaxy clusters. Contrary to expectations, the quasar overdensity is
positioned midway along the filament, not near either LAE node but with an
offset of about 25 cMpc.

Within J0210, we further examine the properties of LAEs at Node-A and Node-B,
containing 78 and 76 member LAEs, respectively. These members are defined by
their separation from points A or B, being less than $15$ $h^{-1}$cMpc, or
approximately $22$ cMpc. We set the criterion as $15$ $h^{-1}$cMpc because
this is a typical size to enclose the entire outskirts of protoclusters for
meaningful comparison to simulations (e.g., Chiang et al., 2013), a potential
scale to have the most significant galaxy–IGM correlation signal (e.g., Cai et
al., 2016, and L21), and a comparable scale to the spatial resolution of IGM
tomography map constructed in Section 4.

Initially, we employ the observed $g$-band magnitudes as proxies for the UV
luminosity of LAEs to conduct comparisons between the two nodes. In Figure 3,
the sizes of circles are encoded by LAEs’ $g$-band magnitudes, with larger
circles indicating higher UV luminosities. Visually, LAEs in Node-A appear
systematically more luminous than those in Node-B. This observation is
statistically supported by a two-sample Kolmogorov–Smirnov (K-S) test,
yielding a $p$-value of approximately $10^{-4}$, which quantitatively confirms
the significant difference in UV luminosity distributions between the two
nodes.

Secondly, we cross-matched LAE candidates with near-UV (NUV; effective
wavelength $\lambda_{\rm eff}=2,314$ Å) detections down to $m_{\rm AB}\approx
22.7$ from the GALEX Release (GR) 6/7 of the Medium-depth Imaging Survey (MIS;
Bianchi et al., 2017). The matching aperture is set to 2″, a size smaller than
the PSF of GALEX NUV images. This selection enables the exclusion of false
matches with high-resolution HSC images upon visual inspection. Only sources
with NUV $<23$ and errors $<0.5$ mag are accepted. Given the Lyman break at
$z>2.15$ falls redward of the NUV wavelength range, no detection is typically
expected for the star-forming galaxies within the shallow GALEX NUV images for
our continuum-faint LAEs, as per the LBG selections (Ly et al., 2009;
Haberzettl et al., 2012). Out of 465 LAEs, 13 were identified with NUV
detection, referred to as NUV-LAEs, and are indicated by thin black crosses in
Figure 3. Of these, 6 belong to Node-A and 2 to Node-B, corresponding to
fractions of 8% and 3% member LAEs, respectively, with NUV detection in each
node. These NUV-LAEs could represent either low-$z$ interlopers or faint AGN
candidates, with a notable clustering in Node-A and a tendency to avoid
Node-B.

A systematic search for Ly$\alpha$ blobs (LABs), galaxies exhibiting luminous
and extended Ly$\alpha$ emission, was conducted using data from the MAMMOTH-
Subaru survey, including the J0210 field (Zhang et al., 2023, M. Li et al. in
prep.). The LABs are selected when their Ly$\alpha$ isophotal area $A_{\rm
iso}$ is $3\sigma$ above the sequence of point-like sources at the given
Ly$\alpha$ luminosity $L_{\rm Ly\alpha}$ and $A_{\rm iso}$ is at least over
about 15 arcsec2. From the LAB catalog, we find ten LABs in J0210, and they
are indicated as the black squares in Figure 3. Notably, six LABs distribute
eastward of the quasar overdensity towards Node-B, two of which (CH-Q05 and
CH-Q08) are also identified as LABs, while the west lobe of the LAE filament
towards Node-A lacks any proximate LABs. LABs are thought to inhabit gas-rich
environments, attributed to Ly$\alpha$ resonant scattering, cooling radiation,
or recombination radiation as potential mechanisms producing extended
Ly$\alpha$ emission (Tumlinson et al., 2017; Ouchi et al., 2020). This
suggests that Node-B may be at an earlier stage of stellar mass buildup
compared to Node-A.

These characteristics indicate that the two cosmic nodes of LAEs may represent
independent substructures with distinct evolutionary histories or at different
evolution stages. Specifically, Node-A appears to be a more mature large-scale
structure compared to Node-B. Intriguingly, the overdensity of proximate
quasars is situated at the transition point of LAE properties.

##  4 3D IGM Hi Tomography Map

###  4.1 Methodology

Figure 4: Examples of the MF-PCA fitting of quasar spectra to derive
Ly$\alpha$ transmission fluctuations for reconstructing 3D Hi tomography maps.
The lower left panel is a collapsed map of Hi absorption within redshift
$2.15<z<2.20$. Red stars mark the LoSs used for tomography, with labeled ones
corresponding to the examples in other panels. Sequentially from the upper
left to the lower right panels, MF-PCA fits of quasar spectra along the LAE
filament are illustrated. Each example highlights the NB387
redshift/wavelength coverage with blue dashed vertical lines. In each
spectrum’s lower panel, thin and thick black lines represent the original and
$15h^{-1}{\rm cMpc}$-smoothed spectra, respectively, and the magenta line
denotes the MF-PCA fit. The upper panel of each spectrum shows the
$15h^{-1}{\rm cMpc}$-smoothed Ly$\alpha$ transmission fluctuation $\delta_{F}$
within the focused wavelength range, with the orange line indicating the
cosmic mean Ly$\alpha$ transmission at the relevant redshift. Grey-shaded
regions are masked out for analysis, with criteria detailed in text.

The overdensities of neutral hydrogen in IGM can be revealed by 3D IGM
tomography, which is reconstructed by the Ly$\alpha$ forest imprinting in the
continuum of background sources. As an initial inspection, we utilize the SDSS
quasar spectra to reconstruct a coarse IGM tomography map with a resolution of
15 $h^{-1}$cMpc. The reconstruction procedures of IGM tomography follow Sun et
al. (2023), which is further based on the methodology proposed in the COSMOS
Lyman-Alpha Mapping And Tomography Observations (CLAMATO; Lee et al., 2016,
2018; Horowitz et al., 2022). We briefly review the essential workflow in the
following paragraphs.

The IGM map is made from Hi absorption left on the spectra of multiple
background quasars in the sky. First, damped Ly$\alpha$ absorption systems
(DLAs), characterized by high-column density hydrogen gas with $N_{\rm
HI}>2\times 10^{20}~{}{\rm cm}^{-2}$ likely associated with galaxies, are
masked in the spectra to avoid biasing the estimate of IGM Hi Ly$\alpha$
absorption. The masked wavelength range is defined by the equivalent width
(EW) of the DLAs and requires a continuum-to-noise ratio (C/N) greater than
1.0. Additionally, potential metal absorption lines, such as Siiv $\lambda
1062$, Nii $\lambda 1084$, Ni $\lambda 1134$, and Ciii $\lambda 1176$, are
masked in the intrinsic quasar spectrum using a $10$ Å mask in the observer’s
frame (Lee et al., 2012).

Figure 5: Collapsed Hi absorption maps from three perspectives over about $60$
cMpc: (a) Lower left: The R.A.–Decl. plane, collapsed across $2.15<z<2.20$;
(b) Upper left: The R.A.–$z$ plane, collapsed within $0.6<{\rm Decl.}<1.2$;
(c) Lower right: The $z$–Decl. plane, collapsed for $32.2<{\rm R.A.}<32.8$.
The colored contour in the background represents Ly$\alpha$ transmission
fluctuation $\delta_{F}$ after a median stacking over the specific dimensions,
with a bluer color indicating higher $\delta_{F}$, indicative of lower Hi
density. Yellow crosses are the proximate quasars, the same as Figure 3, with
numbers indicating their IDs. Blue bars represent the redshift range of LAEs
identified using NB387.

Second, mean-flux regulated principal component analysis (MF-PCA; Lee et al.,
2016; Suzuki et al., 2005) is utilized to model the intrinsic continuum in the
Ly$\alpha$ forest wavelength range (i.e., rest-frame 1,216 – 1,600 Å) from the
spectrum redward of the Ly$\alpha$ emission. Following a conventional least-
squares PCA fit to the spectra of background quasars, the method includes an
adjustment step, aligning the cosmic mean Ly$\alpha$ transmission to correct
the flux uncertainty of SDSS/eBOSS spectra blueward of the Ly$\alpha$ emission
(Lee et al., 2012; Becker et al., 2013). The cosmic mean Ly$\alpha$
transmission follows the empirical relation, $\langle F(z)\rangle_{\rm
cos}=\exp\left[-0.001845(1+z)^{3.924}\right]$, approximating $0.84$ at
$z\approx 2.2$ (Faucher-Giguère et al., 2008). Spectra (1) – (4) in Figure 4
exemplify the fitted spectra and continua for four lines-of-sight (LoSs) along
the LAE filament.

To reconstruct the Hi tomography maps using these LoSs, we define a grid-based
Ly$\alpha$ forest fluctuation as:

$\delta_{F}=F_{\rm Ly\alpha}/\langle F(z)\rangle_{\rm cos}-1.$ (2)

Here, $F_{\rm Ly\alpha}$ represents the Ly$\alpha$ forest transmission at each
pixel on the background quasar spectrum, calculated as $F_{\rm
Ly\alpha}=f_{\rm obs}/f_{\rm int}$, where $f_{\rm obs}$ and $f_{\rm int}$ are
the observed spectrum and the intrinsic continuum from MF-PCA, respectively.
Note that $\delta_{F,~{}{\rm lim}}\approx 0.2$ is the physical upper limit for
fully ionized IGM with $F_{\rm Ly\alpha}=1$, and the measured values higher
than this limit will be replaced by $\delta_{F,~{}{\rm lim}}$. A Wiener filter
is then applied to the LoSs to reconstruct the grid-based Hi tomography map
across the entire 3D comoving volume covering the J0210 field. This volume
spans $30.5<{\rm R.A.}<34.5$ and $-1.1<{\rm Decl.}<2.9$ in the redshift range
$2.10<z<2.25$, consisting of $52\times 52\times 27$ Cartesian grids with a
cell size of $5.114~{}h^{-1}{\rm cMpc}$ along each axis. The Hi tomography map
is created over a broader volume than the J0210 region of interest to mitigate
boundary effects.

Initially, 455 spectra of background quasars from SDSS/eBOSS DR18 were
available for reconstructing the Hi tomography map. Of these, 428 were chosen
as LoSs after visual inspection post MF-PCA fitting to remove those of poor
quality. Ultimately, 242 LoSs, with suitable C/N ratios and redshift $z>2.23$,
were utilized to reconstruct the entire volume. Among them, 23 LoSs fall
within the J0210 HSC FoV providing critical information for our study, and are
depicted as red stars in Figure 4. A summary of the reconstructed Hi
tomography, a collapsed map by stacking the $\delta_{\rm F}$ over redshifts
$z=2.15-2.20$, is presented as the background contour color-coded in Figures 4
and 5.

###  4.2 IGM Tomography Results

####  4.2.1 Average IGM Properties Over $2.15\leq z\leq 2.20$

Figure 6: Profiles of averaged Hi absorption and LAE overdensity along the
large-scale filament. Dotted vertical lines, labeled with numbers, mark the
positions of proximate quasars by their IDs. Red circles indicate averaged Hi
absorption, $\langle\delta_{F}\rangle$, and blue squares denote averaged LAE
overdensity, $\langle\delta_{\rm LAE}\rangle$. Error bars reflect the 32%–68%
scatter from 1,000 Monte Carlo simulations involving random placements. The
origin for the X-axis is designated as point B.

The IGM tomography map is constructed in 3D. For a quantitative analysis, we
begin by examining collapsed maps on the R.A.–Decl. plane, where $\delta_{\rm
F}$ is averaged over the redshift range $2.15\leq z\leq 2.20$ using median
stacking, as depicted in the main panel of Figure 5. The top and side views,
shown in the upper and right-hand panels, illustrate the quasar distribution
along the redshift space. It’s important to note that these collapsed maps
average over dimensions of approximately $60$ cMpc each, determined by the
redshift coverage of NB387–Ly$\alpha$. Therefore, the collapsed $\delta_{\rm
F}$ values do not imply a direct connection to individual quasars but rather
indicate a statistical correlation on a scale of about 60 cMpc.

The IGM tomography map indicates that the hydrogen gas in Node-B is more
neutral than the cosmic average, with $\delta_{\rm F}\approx-0.1$, whereas
Node-A is situated in an almost fully ionized region, with $\delta_{\rm
F}\approx 0.2$. Notably, a bimodal characteristic in the ionization state of
the IGM also manifests along the LAE filament, reinforcing the distinct galaxy
properties observed between the nodes.

The quasar overdensity exhibits substructures in redshift space, with the
proximate quasars divisible into two groups: one comprises (CH-)Q03, Q04, Q05,
Q07, and Q08 at $z\approx 2.17$, and the other includes (CH-)Q01, Q02, Q06,
Q09, Q10, and Q11 at $z\approx 2.19-2.20$. Each group occupies a volume of
$\left(20~{}{\rm cMpc}\right)^{3}$, with a separation of approximately $20$
cMpc between the two. A visual inspection reveals that most proximate quasars
align with the interface between ionized and neutral regions. The exception is
CH-Q03, which is situated in a fully ionized region with $\delta_{\rm
F}\approx 0.2$ across all dimensions, maintaining a distance of at least $15$
cMpc from the $\delta_{\rm F}=0$ boundary.

For a more quantitative analysis, a Ly$\alpha$ transmission fluctuation
profile is constructed along the LAE filament from Node-B to Node-A, using
point B as the origin. The filament A–B is inclined at 23.46 degrees. The
width of the calculating box, perpendicular to filament A–B, is set to 15
$h^{-1}$cMpc, aligning with the smoothing scale of the IGM tomography map.
This profile is depicted by red circles with error bars and a curve in Figure
6, while the locations of proximate quasars are indicated by vertical dotted
black lines with ID labels. Additionally, the LAE overdensity profile within
the same volume is represented by blue squares with error bars and curves.
Errors for both profiles are derived from the 32% and 68% ranks in Monte Carlo
simulations, which involve randomly positioning the calculating boxes 1,000
times with random orientations.

This profile quantitatively highlights the transition observed in our visual
inspections, depicting the ionization stage shift of the IGM hydrogen from
ionized gas in Node-A to neutral Hi in Node-B along the LAE filament. The
proximate quasars exhibit a $\sim 25$ cMpc offset from the peaks of either the
Ly$\alpha$ absorption or Ly$\alpha$ transmission in the IGM. These quasars are
situated near the ionization frontiers where $\langle\delta_{\rm
F}\rangle\approx 0$, predominantly within relatively ionized regions.

The signal indicating an ionization stage transition between Node-A and
Node-B, marked by a decrement in $\Delta\delta_{\rm F}\gtrsim 0.2\pm 0.1$, is
significant on the scale of $50~{}h^{-1}{\rm cMpc}$. This observation exceeds
the limitations imposed by the IGM tomography map’s spatial resolution, which
is reconstructed from background SDSS quasar spectra with a mean LoS sampling
separation of $D_{\perp}\approx 15~{}h^{-1}{\rm cMpc}$. The observed change
distinctly stands out from general fields, where $\Delta\delta_{\rm F}\leq
0.03$, utilizing similar SDSS LoS samples and methodology (Sun et al., 2023).

####  4.2.2 Inspection on Idividual Redshifts

Figure 7: Hi tomography maps across different redshift slices. The colored
contour retains the interpretation from Figure 4, applied to distinct redshift
bins. Sequentially from the upper left to the bottom right panel, redshifts
increase from $z=2.15$ to $z=2.20$, with intervals indicated in each panel.
Yellow crosses mark proximate quasars within the redshift bins, while red
dashed-dotted circles represent the photoionizing radius $r_{\rm phot}$.
White-shaded areas denote regions where ionizing photons from these luminous
SDSS quasars surpass those from the UVB. Due to the absence of spectroscopic
redshifts for LAEs, their overdensity is consistently depicted as dotted black
contours across slices.

As indicated in Section 4.2.1, collapsed maps do not suit the examination of
individual targets. However, they hint at a general trend of proximate quasars
congregating near ionizing frontiers. Therefore, we delve into individual
redshift slices derived directly from the LoSs, constrained by the cell size
along the redshift axis at $5.114~{}h^{-1}{\rm cMpc}$ (or $7.3$ cMpc). These
slicing maps are displayed in Figure 7, with the colored background contour
maintaining the same interpretation as in Figure 5, but applied to maps at
specific redshifts instead of a collapsed overview. From the upper left to the
lower right panel, the redshift increases from $z=2.150$ to $z=2.200$, with
intervals noted in each panel. Proximate quasars within each redshift slice
are marked as yellow crosses. Meanwhile, since LAE candidates are
photometrically identified within $z=2.15-2.20$ without precise redshifts,
their overdensity is depicted in every panel using black dashed contours.

The IGM Hi distribution exhibits variation across different redshifts,
revealing a 3D structure. In every redshift slice, IGM hydrogen in Node-A
appears nearly fully ionized, indicating a confidently Hi-poor environment
within this LAE overdensity. Conversely, the ionization state of IGM in Node-B
varies with redshift. Should LAEs within Node-B cluster around $2.17\lesssim
z\lesssim 2.19$, the IGM Hi overdensity might align with galaxy distributions
in a massive structure. This pattern also suggests a transition in the IGM
ionization stage along the filament. The spectroscopic follow-up to determine
precise LAE redshifts is crucial for accurately characterizing this vast
structure.

The narrative deepens for the proximate quasars. The slicing maps reiterate
that 10 out of 11 clustering quasars, with the exception of CH-Q03, are
situated near the boundary between neutral and ionized regions at every
redshift, reinforcing the result that quasars are at or just beyond the
ionizing frontiers. This result aligns with the $\delta_{\rm F}$ profile
described in the collapsed map discussed in Section 4.2.1. Thus, it is posited
that the clustering quasars may be intricately linked to large-scale IGM
ionization processes.

##  5 Discussion

###  5.1 The Association between Quasars and IGM Hi

Figure 8: Bolometric luminosity ($L_{\rm bol}$) of proximate quasars plotted
against the local Ly$\alpha$ transmission fluctuation, $\delta_{F}$. Error
bars reflect a 50% uncertainty in converting $L_{\rm 1350}$ to $L_{\rm bol}$.
The curve illustrates a power-law relation between $L_{\rm bol}$ and
$\delta_{\rm F}$. Confidence intervals are determined through Monte Carlo
simulations that vary each data point within a Gaussian distribution using the
stated uncertainty; the 32% and 68% ranks define the lower and upper bounds,
respectively.

Our earlier findings indicate that proximate quasars may be associated with
IGM Hi within the underlying structures in 3D space. Beyond the spatial
distribution, we examine the correlation between quasars’ bolometric
luminosity $L_{\rm bol}$ and $\delta_{\rm F}$ at their respective locations,
and the $\delta_{\rm F}$ is given by the nearest grid constructed on 15
h-1cMpc in Section 4. This analysis, depicted in Figure 8, reveals a
significant positive correlation between $L_{\rm bol}$ and $\delta_{\rm F}$,
as evidenced by a Spearman rank test yielding a correlation coefficient
$\rho_{\rm S}=0.75$ with a $P$-value of $0.007$. A power law fit between
$L_{\rm bol}$ and $\delta_{\rm F}$ establishes the relation:

$\text{log}\left(L_{\rm bol}\right)=2.42^{+0.39}_{-0.63}\times\delta_{\rm
F}+45.93^{+0.10}_{-0.12},$ (3)

with errors derived from the 32%–68% ranks in a Monte Carlo simulation by
fluctuating data points with 50% uncertainties and repeating 1,000 times. The
fitted relation and its uncertainty are illustrated in Figure 8 by the red
dashed-dotted curve and shaded area, respectively. This relation further
substantiates a strong correlation, indicating that more luminous quasars
reside in more transparent IGM environments, thereby suggesting a potential
physical linkage between quasars and the IGM on scales exceeding
$15~{}h^{-1}$cMpc.

Figure 9: Reconstructed maps of baryon matter fluctuation, $\Delta_{b}$,
across different redshift slices. Redder shades signify higher $\Delta_{b}$
values, representing denser baryonic matter, encompassing both ionized and
neutral hydrogen in the IGM. Panel arrangements and symbols used are
consistent with those in Figure 7. Figure 10: IGM tomography map for
$2.189<z<2.194$, akin to Figure 7. Red circles and white spheres denote
photoionization sources from both quasars and NUV-LAEs. Quasars are marked as
yellow crosses with IDs, while NUV-LAEs are represented by thin black crosses.
Figure 11: Star formation rate (SFR) distribution of LAEs, with blue hatched
and shaded histograms representing member LAEs in Node-A and Node-B,
respectively. The light and dark gray shaded band suggest the standard
deviation of SFR in the LAEs identified in general fields at $z\approx 2.1$
and 2.2 (Guaita et al., 2011; Kusakabe et al., 2018).

To discern the photoionizing influence of quasars, a comparison between their
contributions and that from UV background (UVB) can be insightful, e.g.,
Kashikawa et al. (2007). We estimate the ionizing photons from quasars in
their proximity Hi following Arrigoni Battaia et al. (2015); Pezzulli &
Cantalupo (2019). The photoionization rate for an individual quasar is given
by:

$\begin{split}\Gamma_{\rm HI,~{}i}^{\rm QSO}=8.8\times
10^{-7}\xi_{i}\left(\frac{r_{i}}{10~{}{\rm kpc}}\right)^{-2},\end{split}$ (4)

where the $\Gamma_{\rm HI,i}^{\rm QSO}$ is the photoionization rate of a
specific quasar $i$, $\xi_{i}$ is a factor scaled by this quasar’s luminosity
and continuum slope, and $r_{i}$ is the distance between the calculating
location and this quasar. We note that this $\Gamma_{\rm HI,~{}i}^{\rm QSO}$
estimate assumes that the quasar radiation is spherically symmetric and their
luminosity is time-invariant, which may not be a case in reality, but it still
provides instructive ideas.

Assuming a power-law continuum for wavelengths shorter than the Lyman limit,
$\xi_{i}$ is calculated by rescaling observed $i$-band magnitudes of reference
quasars from Borisova et al. (2016):

$\xi_{i}=\left(\frac{L_{\rm\nu,LL,i}}{5.2\times 10^{31}~{}{\rm
erg~{}s^{-1}~{}Hz^{-1}}}\right)\left(\frac{\beta_{\rm
QSO,i}+3}{4.7}\right)^{-1},$ (5)

where $L_{\rm\nu,LL,i}$ is the monochromatic continuum luminosity at the Lyman
limit (rest-frame 912 Å), and $\beta_{\rm QSO}$ is the UV continuum slope
determined in Section 3.1. By integrating Equations (4) and (5), we can
identify the radius at which a specified photoionizing rate occurs, defining a
photoionizing radius $r_{\rm phot}$ where for distances $r<r_{\rm phot}$:

$\Gamma_{\rm HI,~{}i}^{\rm QSO}>\Gamma_{\rm HI}^{\rm UVB},$ (6)

and $\Gamma_{\rm HI}^{\rm UVB}\sim 10^{-12}~{}{\rm s}^{-1}$ at $z\approx 2$
(Haardt & Madau, 2012). Thus, within regions closer than $r_{\rm phot}$ to
proximate quasars, the ionizing photons from these luminous SDSS quasars
significantly exceed those from the UVB. Figure 7 illustrates the results,
with the radius of red dotted circles centered on quasars representing $r_{\rm
phot}$. Consequently, photoionization within the white-shaded areas is
predominantly driven by quasar ionization. We note that the shape of these
spherical regions can be altered if the emission from quasars is anisotropic.

Clearly, the most luminous quasars like CH-Q01 and CH-Q02 are capable of
influencing the IGM ionization stage up to a distance of 20 cMpc. Although the
relatively faint quasars may have a more limited impact on their surroundings,
their affecting areas can be overlapped to make a coherent contribution to
shape the topology of nearby ionizing frontiers. However, even considering
their joint efforts, the $r_{\rm phot}$ of the known eleven quasars likely
have little effect on the ionization state of the distant Nodes A and B.

The quasar and IGM Hi spatial distribution, the correlation between quasar
$L_{\rm bol}$ and $\delta_{F}$, and the analysis of $r_{\rm phot}$
collectively indicate that the IGM Hi on scales of tens cMpc is possibly
modulated by the clustering quasars within the Cosmic Himalayas. Albeit the
most luminous quasars may have a more significant impact, as illustrated in
the $L_{\rm bol}$ – $\delta_{F}$ relation and $r_{\rm phot}$, a single quasar
cannot do the whole job to generate a large scale ionized IGM gas structure as
observed beyond 20 cMpc. Multiple quasars being switched on simultaneously are
essential to produce this peculiar ionizing structure.

###  5.2 The Potential Baryonic Matter Distribution

We have noted that the quasar overdensity does not coincide with the peaks of
LAE density or the peaks of IGM Hi absorption/transmission within the J0210
field. However, the distribution of total baryonic matter, including ionized
hydrogen, has yet to be clarified. Utilizing the IGM Hi tomography map, it is
possible to reconstruct the distribution of baryonic matter by accounting for
the photoionization effects attributed to quasars.

In the IGM, the Hi Ly$\alpha$ optical depth, $\tau_{\rm Ly\alpha}$, is related
to the Hi column density, NHI, through the equation:

$\tau_{\rm Ly\alpha}=\sigma_{\rm Ly\alpha}{\rm N}_{\rm HI},$ (7)

where $\sigma_{\rm Ly\alpha}$ is the Ly$\alpha$ cross section. The NHI can be
determined using:

$\begin{split}{\rm N}_{\rm HI}&=l_{\rm Ly\alpha}n_{\rm HI}\\\ &=l_{\rm
Ly\alpha}n_{\rm H}\Delta_{b}x_{\rm HI},\end{split}$ (8)

in which $l_{\rm Ly\alpha}$ represents the line-of-sight length of the medium
that Ly$\alpha$ photons pass through, $n_{\rm HI}$ and $n_{\rm H}$ are the
number density of neutral hydrogen and cosmic hydrogen, $\Delta_{b}$ signifies
the baryon fluctuation, and $x_{\rm HI}$ is the Hi fraction. As indicated in
Section 4.1, at $z\sim 2.2$, cosmic Hi is predominantly ionized with $\langle
F(z)\rangle_{\rm cos}\sim 0.84$. Consequently, under conditions where
photoionization reaches equilibrium, it is reasonable to consider case A
(optically thin) hydrogen recombination. Assuming the equilibrium can be
established on a very short time-scale, $x_{\rm HI}$ can be approximated
following Meiksin (2009):

$x_{\rm HI}=\frac{\Gamma_{\rm rec}}{\Gamma_{\rm HI}^{\rm tot}}=\frac{n_{\rm
H}\Delta_{b}\alpha(T)}{\Gamma_{\rm HI}^{\rm tot}},$ (9)

where $\Gamma_{\rm rec}$ and $\Gamma_{\rm HI}^{\rm tot}$ represent the
hydrogen recombination rate and the total Hi photoionization rate,
respectively. The $\alpha(T)$ is the Case A hydrogen recombination coefficient
at temperature $T$, and typically, $T\approx 10^{4}$ K in IGM. Merging
Equations (7), (8), and (9) leads to:

$\tau_{\rm Ly\alpha}=\frac{\sigma_{\rm Ly\alpha}l_{\rm Ly\alpha}n_{\rm
H}^{2}\Delta_{b}^{2}\alpha(T)}{\Gamma_{\rm HI}^{\rm tot}}.$ (10)

Hence, Equation (10) implies:

$\begin{split}\tau_{\rm Ly\alpha}&\propto\frac{\Delta_{b}^{2}}{\Gamma_{\rm
HI}^{\rm tot}}\\\ &\propto\frac{\Delta_{b}^{2}}{B\times\Gamma_{\rm HI}^{\rm
UVB}},\end{split}$ (11)

where the boost factor $B$ incorporates the quasar contribution:

$B=\left(\Gamma_{\rm HI}^{\rm UVB}+\Gamma_{\rm HI}^{\rm
QSO}\right)/\Gamma_{\rm HI}^{\rm UVB}.$ (12)

The $\Gamma_{\rm HI}^{\rm UVB}$ is the UV background photoionization rate and
$\Gamma_{\rm HI}^{\rm QSO}$ denotes the photoionization rate from proximate
quasars. Meanwhile, by definition, the Hi optical depth is

$\begin{split}\tau_{\rm Ly\alpha}&=-{\rm ln}\left(F\right)\\\ &=-{\rm
ln}\left[\left(1+\delta_{F}\right)\left<F\right>_{\rm cos}\right],\end{split}$
(13)

with $F$, $\left<F\right>_{\rm cos}$, and $\delta_{F}$ representing the local
Ly$\alpha$ transmission, cosmic mean Ly$\alpha$ transmission, and local
Ly$\alpha$ transmission fluctuation, respectively. By connecting Equations
(11) and (13), a proportional relation is established:

$\Delta_{b}\propto\left\\{-\Gamma_{\rm HI}^{\rm tot}\times{\rm
ln}\left[\left(1+\delta_{F}\right)\left<F\right>_{\rm
cos}\right]\right\\}^{1/2}.$ (14)

When $\delta_{F}=0$ and $\Gamma_{\rm HI}^{\rm tot}=\Gamma_{\rm HI}^{\rm UVB}$,
there should be $\Delta_{b}=\rho_{b}/\left<\rho_{b}\right>=1$, with $\rho_{b}$
and $\left<\rho_{b}\right>$ denoting the local density and cosmic mean density
of baryons. Hence, we get:

$\begin{split}\Delta_{b}&=\left\\{\frac{\Gamma_{\rm HI}^{\rm tot}\times{\rm
ln}\left[\left(1+\delta_{F}\right)\left<F\right>_{\rm cos}\right]}{\Gamma_{\rm
HI}^{\rm UVB}\times{\rm ln}\left<F\right>_{\rm cos}}\right\\}^{1/2}\\\
&=\left\\{B\times\left[1+\frac{{\rm ln}\left(1+\delta_{F}\right)}{{\rm
ln}\left<F\right>_{\rm cos}}\right]\right\\}^{1/2},\end{split}$ (15)

revealing that the baryonic fluctuation $\Delta_{b}$ can be deduced from
Ly$\alpha$ transmission fluctuation $\delta_{F}$ and the boost factor $B$ due
to quasar proximity. To compute $B$ across the IGM tomography map’s grids, the
collective photoionization rate from all eleven quasars is derived from
Equation (4):

$\begin{split}\Gamma_{\rm HI}^{\rm QSO}&=\sum_{i}^{11}\Gamma_{\rm HI,i}^{\rm
QSO}\\\ &=\sum_{i}^{11}8.8\times 10^{-7}\xi_{i}\left(\frac{r_{i}}{10~{}{\rm
kpc}}\right)^{-2}.\end{split}$ (16)

The 3D tomography map depicting the baryonic matter fluctuation, $\Delta_{b}$,
is presented in Figure 9, with a redder color indicating denser matter. This
aligns with the discussions in Section 5.1 about quasars of relatively low
luminosity having minimal impact on the IGM ionizing stage on large scales.
The reconstructed $\Delta_{b}$ distribution does not exhibit a significant
deviation from $\delta_{F}$. Nonetheless, due to the predominant
photoionization influence of CH-Q01 and CH-Q02 in their vicinity, $\Delta_{b}$
is substantially elevated around these luminous quasars compared to Hi,
hinting at that the most luminous quasars are possibly situated in dense
environments of matters, as indirectly traced by IGM Hi.

We note that although the reconstructed $\Delta_{b}$ is called baryon
fluctuation here, it is indeed reconstructed purely based on Hi distribution,
which may have a bias factor different from that of Hii and affect the result.
It’s also important to recognize that the $\Delta_{b}$ reconstruction
presented only considers the ionization effects of known quasars from
SDSS/eBOSS, and it’s possible that additional ionizing sources have not been
accounted for. This potential oversight and its implications will be further
explored in Section 5.3.

###  5.3 The Lack of IGM Hi in Node-A

The reconstructed $\Delta_{b}$ map has offered valuable insights into the
correlation between quasars and the IGM. Nonetheless, the structure of Hi
around LAE Node-A remains perplexing. As illustrated in Figures 5 and 7, the
$\delta_{F}$ within Node-A is consistently around 0.2, indicating a highly
transparent structure at $z\approx 2.2$ amidst dense LAEs. Conversely,
luminous quasars identified in SDSS/eBOSS are situated at least 25 cMpc away
from this region, rendering their influence on the IGM environment around
Node-A negligible, as quantified by $r_{\rm phot}$. Furthermore, the
$\Delta_{b}$ map in Figure 9 reveals a pronounced density decline between the
quasar overdensity and Node-A within $2.189<z<2.194$. This abrupt decrement
seems physically implausible, hinting at the presence of unaccounted ionizing
sources near Node-A.

One possibility is that the optically obscured AGNs responsible for the
ionizing photon budget are overlooked in the SDSS/eBOSS quasar selection
process. This brings attention to NUV-LAEs, considered potential AGN
candidates due to their NUV excess in the Lyman break regime. Assuming, in the
most extreme scenario, that all NUV-LAEs are indeed AGNs and cluster around
$2.189<z<2.194$, we can calculate the photoionization radius $r_{\rm
phot}^{\rm NUV-LAE}$ for these sources using a methodology akin to that
discussed in Section 5.1. The $\xi$ for NUV-LAEs were acquired through
extrapolating an equivalent $i$-band from the NUV band assuming a power-law
continuum with typical $\beta_{\rm QSO}$. The $r_{\rm phot}^{\rm NUV-LAE}$,
represented by the radius of red dashed-dotted circles centered on black
crosses, is illustrated in Figure 10. This analysis suggests that, even in the
most extreme case, the NUV-LAEs alone are insufficient to account for the
fully ionized IGM around Node-A.

Another potential contributor could be excess star-formation activity in star-
forming galaxies within these regions, if we assume the escape fraction of
ionizing photons can be comparable in the two structures. To explore this, we
assess the UV star formation rate (SFR) of member LAEs in Node-A and Node-B,
accounting for 71 and 74 member LAEs, respectively, after excluding the NUV-
LAEs. We derive the UV SFR based on the scaling relation given by Madau et al.
(1998):

${\rm SFR}/{\rm M_{\odot}~{}yr^{-1}}=\frac{L_{\rm 1500}/{\rm erg}~{}{\rm
s}^{-1}~{}{\rm Hz}^{-1}}{8.0\times 10^{27}},$ (17)

where $L_{\rm 1500}$ represents the UV flux at $\lambda_{0}=1,500$ Å. A
correction factor of $0.08$ is applied to the HSC $g$-band magnitude to
account for the wavelength offset, assuming a typical UV slope of LAEs as
$\beta=-1.6$ (Kusakabe et al., 2018). The SFR distributions for LAEs are
depicted in Figure 11, showing a noticeable difference, as confirmed by a two-
sample K-S test with a $p$-value $<10^{-4}$. Compared to the LAEs in general
fields, Node-B members share a SFR distribution almost similar to Guaita et
al. (2011), and slightly higher than Kusakabe et al. (2018), but Node-A
members have systemically higher star formation activity than both. As a
result, the ionizing photon production from star-forming activities in Node-A
is systematically higher than that in Node-B by approximately 0.3 dex. This
significant amount suggests that the elevated SFR around Node-A could partly
account for the ionized state of the IGM in its vicinity.

Furthermore, the combined effects of obscured AGNs and star-forming galaxies
might collaboratively enhance the photoionization within Node-A. The potential
for a more refined IGM tomography, aiming for a resolution closer to $\sim 3$
cMpc in future studies with the upcoming Subaru Prime Focus Spectrograph
(PFS), promises to delineate their contributions more precisely.

###  5.4 Possible Triggers of The Densest Quasar

Quasars are often considered the aftermath of gas-rich major mergers (Hopkins
et al., 2006). This hypothesis has garnered observational support at $z<1$
based on the statistical analysis based on galaxy morphology (e.g., Goulding
et al., 2018). Therefore, active quasar activities can be a proxy of rich
environments with galaxy overdensities. Nonetheless, the quasars within the
Cosmic Himalayas deviate from this pattern. They do not coincide with the
density peaks of LAEs. Instead, they appear prominently positioned at the
midpoint of a large-scale filament. This fact complicates the conventional
understanding of quasars’ origins in major mergers, suggesting a more nuanced
interplay between quasar activity and its surrounding cosmic structures.

To reconcile the observed quasars of the Cosmic Himalayas with the standard
model that links quasar activity to major mergers, two potential explanations
can help resolve the apparent discrepancy. First, major mergers, capable of
triggering quasar activity, are also expected to generate significant amounts
of dust, rendering the quasar’s environment highly dusty (Hopkins et al.,
2008). This could hinder the detection of resonant Ly$\alpha$ photons from
galaxies, leading to an observational bias where many LAEs proximate to
quasars might be overlooked in NB imaging (e.g., Shimakawa et al., 2017; Ito
et al., 2021). It is noteworthy that a clear discrepancy between dusty star-
forming galaxies and the more typical star-forming galaxies with H$alpha$
emission show out in two massive protoclusters also traced by multiple quasars
(Zhang et al., 2022). Alternatively, the influence of quasars on their
surroundings might manifest as strong feedback mechanisms that heat the cold
gas in their vicinity, thereby inhibiting the recent formation of galaxies
(Dong et al., 2023, 2024). This negative feedback could physically limit the
presence and detectability of young LAEs near quasars.

If these scenarios hold, the actual galaxy overdensity within the Cosmic
Himalayas could be significantly higher than what is currently observed.
Subsequently, the estimated present-day halo mass of $M_{h}(z=0)\sim
10^{14.8}~{}{\rm M}_{\odot}$ for Node-A, derived from the distribution of
detectable LAEs, would merely represent a lower limit of the true galaxy
density and mass within this region.

Another possibility of a collision between galaxy (proto-)clusters, extending
beyond the conventional galaxy merger origin, emerges from the structure
delineated by LAEs. The two nodes seem to be linked by a cosmic filament, with
quasars positioned perpendicularly to this filament both in the projected
plane and in redshift space, as illustrated in Figure 5. This configuration
bears resemblance to morphological analogs in the local universe, e.g., the
Bullet Cluster (1E 0657-558) at $z=0.296$ (Tucker et al., 1998; Clowe et al.,
2004) and the Sausage and Toothbrush (CIZA J2242.8+5301 and 1RXS J0603.3+4213;
Stroe et al., 2015). These structures represent the aftermath of a collision
between two galaxy clusters, leading to intensive interaction of their
intracluster gases. Consequently, the gas, having lost significant momentum,
was separated from the dark matter halos, manifesting as strong X-ray emission
in shock regions. Notably, an elevated proportion of X-ray luminous AGNs has
been observed within the Bullet Cluster (Puccetti et al., 2020), hinting at
mechanisms that might facilitate increased gas inflow to activate AGNs during
such galactic cluster collisions.

It is crucial to recognize the disparity in spatial scale between the Bullet
Cluster (approximately $2~{}{\rm cMpc}$) and the LAE filament in J0210 (around
$100~{}{\rm cMpc}$). Nonetheless, given the potential size evolution of
protoclusters, the LAE filament with its two nodes could represent a
progenitor of structures similar to the Bullet Cluster at $z=2.2$ (Chiang et
al., 2013). To validate this hypothesis, further investigations, such as deep
X-ray observations of the diffuse intergalactic gas and detailed kinematic
studies, are essential.

##  6 Summary

In this paper, we report the discovery of an extraordinary cosmological
structure in the J210 field, the Cosmic Himalayas, traced by both the presence
of a significant quasar concentration and the distribution of LAEs indicative
of large-scale cosmic filaments. Utilizing data from SDSS/eBOSS for quasar
identification and 3D IGM tomography and Subaru/HSC narrowband imaging for LAE
selection, we present a comprehensive analysis revealing the spatial
distribution of quasars, LAEs, and the IGM Hi. Important results include:

  * •

The discovery of an unprecedented concentration of eleven type-1 AGNs at
$z\approx 2.16-2.20$ within a compact $(40~{}\mathrm{cMpc})^{3}$ region,
marking the densest quasar overdensities $\delta_{\rm QSO}\approx 30$ with
$16.9\sigma$ significance at $z>2$, situated at the midpoint of a $100$ cMpc
LAE filament. This structure connects two protocluster nodes, but both show an
approximately 25 cMpc offset to the quasar density peak. The inspection of LAE
properties, including luminosity, potential AGN activity, and extended
Ly$\alpha$ emission, reveals significant differences between the two nodes.

  * •

A 3D IGM Hi tomography map is reconstructed from SDSS/eBOSS, which reveals a
dual ionization state along the LAE filament: nearly fully ionized IGM around
Node-A at NW, contrasting with a transition to more neutral states towards the
Node-B at SE. This variance corroborates a similarly distinctive
characteristic found by LAE properties between the two nodes.

  * •

The analysis with the 3D IGM Hi tomography map indicates that quasars are
preferentially located near the boundaries between ionized and neutral IGM
regions. This spatial distribution, alongside a significant positive
correlation between quasar bolometric luminosity ($L_{\rm bol}$) and local Hi
transmission fluctuation ($\delta_{F}$), underscores a potential physical
association between clustered quasar activity and large-scale IGM ionization.

  * •

The offset between the peaks of quasars and LAEs challenges conventional
theories of quasar activation solely through major mergers. It either suggests
that a significant amount of other galaxies are missing or killed in the
quasar proximity, or it hints at complex triggering mechanisms, potentially
including intense environmental interactions reminiscent of cosmic structures
like the Bullet Cluster in the local universe.

We have discussed some potential scenarios with mixtures of Hi overdensities
and quasar photoionization to explain a part of the observed results, but no
conclusive theory can fully resolve the emergence of Cosmic Himalayas. The
presented findings, supported by rigorous data analyses, not only shed light
on the relationships between quasars, galaxies, and the IGM in an extreme
cosmic structure, the Cosmic Himalayas, but also pave the way for future
multi-wavelength, high-resolution observational campaigns aimed at unraveling
the intricate evolution of such remarkable cosmic phenomena.

We thank Yechi Zhang, Yuki Isobe, Yi Xu, Hiroya Umeda, Khee-Gan Lee, Antonio
Pensabene, Titouan Lazeyras, Marta Galbiati, Matteo Viel, Xianzhong Zheng,
Veronica Strazzullo, Maurilio Pannella, Alex Saro, Stefano Borgani, Laura
Pentericci, Mark Dickinson, Lorenzo Napolitano, Paolo Cassata, Giulia
Rodighiero, Olga Cucciati, Roberto Decarli, Akio K. Inoue, Yuqing Lou, Guochao
Sun for helpful discussions. This work was supported by JSPS KAKENHI Grant
Number JP24K17084. Y.L. acknowledges the travel support from the program
International Leading Research “Comprehensive understanding of the formation
history of structures in the Universe”.

Funding for the Sloan Digital Sky Survey V has been provided by the Alfred P.
Sloan Foundation, the Heising-Simons Foundation, the National Science
Foundation, and the Participating Institutions. SDSS acknowledges support and
resources from the Center for High-Performance Computing at the University of
Utah. The SDSS web site is www.sdss.org.

SDSS is managed by the Astrophysical Research Consortium for the Participating
Institutions of the SDSS Collaboration, including the Carnegie Institution for
Science, Chilean National Time Allocation Committee (CNTAC) ratified
researchers, the Gotham Participation Group, Harvard University, Heidelberg
University, The Johns Hopkins University, L’Ecole polytechnique fédérale de
Lausanne (EPFL), Leibniz-Institut für Astrophysik Potsdam (AIP), Max-Planck-
Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für
Extraterrestrische Physik (MPE), Nanjing University, National Astronomical
Observatories of China (NAOC), New Mexico State University, The Ohio State
University, Pennsylvania State University, Smithsonian Astrophysical
Observatory, Space Telescope Science Institute (STScI), the Stellar
Astrophysics Participation Group, Universidad Nacional Autónoma de México,
University of Arizona, University of Colorado Boulder, University of Illinois
at Urbana-Champaign, University of Toronto, University of Utah, University of
Virginia, Yale University, and Yunnan University.

This research is based in part on data collected at the Subaru Telescope,
which is operated by the National Astronomical Observatory of Japan. We are
honored and grateful for the opportunity of observing the Universe from
Maunakea, which has the cultural, historical, and natural significance in
Hawaii. The Hyper Suprime-Cam (HSC) collaboration includes the astronomical
communities of Japan and Taiwan, and Princeton University. The HSC
instrumentation and software were developed by the National Astronomical
Observatory of Japan (NAOJ), the Kavli Institute for the Physics and
Mathematics of the Universe (Kavli IPMU), the University of Tokyo, the High
Energy Accelerator Research Organization (KEK), the Academia Sinica Institute
for Astronomy and Astrophysics in Taiwan (ASIAA), and Princeton University.
This paper makes use of software developed for the Large Synoptic Survey
Telescope. We thank the LSST Project for making their code available as free
software at http://dm.lsst.org.

Data analysis was in part carried out on the Multi-wavelength Data Analysis
System operated by the Astronomy Data Center (ADC), National Astronomical
Observatory of Japan.

Facility: Sloan, Subaru (HSC), GALEX

## References

  * Aihara et al. (2018) Aihara, H., Arimoto, N., Armstrong, R., et al. 2018, PASJ, 70, S4, doi: 10.1093/pasj/psx066
  * Aihara et al. (2019) Aihara, H., AlSayyad, Y., Ando, M., et al. 2019, PASJ, 71, 114, doi: 10.1093/pasj/psz103
  * Aihara et al. (2022) —. 2022, PASJ, 74, 247, doi: 10.1093/pasj/psab122
  * Arrigoni Battaia et al. (2015) Arrigoni Battaia, F., Hennawi, J. F., Prochaska, J. X., & Cantalupo, S. 2015, ApJ, 809, 163, doi: 10.1088/0004-637X/809/2/163
  * Becker et al. (2013) Becker, G. D., Hewett, P. C., Worseck, G., & Prochaska, J. X. 2013, MNRAS, 430, 2067, doi: 10.1093/mnras/stt031
  * Bertin & Arnouts (1996) Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393, doi: 10.1051/aas:1996164
  * Bianchi et al. (2017) Bianchi, L., Shiao, B., & Thilker, D. 2017, ApJS, 230, 24, doi: 10.3847/1538-4365/aa7053
  * Borisova et al. (2016) Borisova, E., Cantalupo, S., Lilly, S. J., et al. 2016, The Astrophysical Journal, 831, 39
  * Bosch et al. (2018) Bosch, J., Armstrong, R., Bickerton, S., et al. 2018, PASJ, 70, S5, doi: 10.1093/pasj/psx080
  * Cai et al. (2016) Cai, Z., Fan, X., Peirani, S., et al. 2016, ApJ, 833, 135, doi: 10.3847/1538-4357/833/2/135
  * Cai et al. (2017) Cai, Z., Fan, X., Bian, F., et al. 2017, ApJ, 839, 131, doi: 10.3847/1538-4357/aa6a1a
  * Cen & Ostriker (2000) Cen, R., & Ostriker, J. P. 2000, ApJ, 538, 83, doi: 10.1086/309090
  * Chiang et al. (2013) Chiang, Y.-K., Overzier, R., & Gebhardt, K. 2013, ApJ, 779, 127, doi: 10.1088/0004-637X/779/2/127
  * Clowe et al. (2004) Clowe, D., Gonzalez, A., & Markevitch, M. 2004, ApJ, 604, 596, doi: 10.1086/381970
  * Clowes (1986) Clowes, R. G. 1986, MNRAS, 218, 139, doi: 10.1093/mnras/218.2.139
  * Clowes & Campusano (1991) Clowes, R. G., & Campusano, L. E. 1991, MNRAS, 249, 218, doi: 10.1093/mnras/249.2.218
  * Clowes et al. (2012) Clowes, R. G., Campusano, L. E., Graham, M. J., & Söchting, I. K. 2012, MNRAS, 419, 556, doi: 10.1111/j.1365-2966.2011.19719.x
  * Clowes et al. (2013) Clowes, R. G., Harris, K. A., Raghunathan, S., et al. 2013, MNRAS, 429, 2910, doi: 10.1093/mnras/sts497
  * Crampton et al. (1989) Crampton, D., Cowley, A. P., & Hartwick, F. D. A. 1989, ApJ, 345, 59, doi: 10.1086/167881
  * Dawson et al. (2013) Dawson, K. S., Schlegel, D. J., Ahn, C. P., et al. 2013, AJ, 145, 10, doi: 10.1088/0004-6256/145/1/10
  * Dawson et al. (2016) Dawson, K. S., Kneib, J.-P., Percival, W. J., et al. 2016, AJ, 151, 44, doi: 10.3847/0004-6256/151/2/44
  * Dong et al. (2023) Dong, C., Lee, K.-G., Ata, M., Horowitz, B., & Momose, R. 2023, ApJ, 945, L28, doi: 10.3847/2041-8213/acba89
  * Dong et al. (2024) Dong, C., Lee, K.-G., Davé, R., Cui, W., & Sorini, D. 2024, arXiv e-prints, arXiv:2402.13568, doi: 10.48550/arXiv.2402.13568
  * Fabian (2012) Fabian, A. C. 2012, ARA&A, 50, 455, doi: 10.1146/annurev-astro-081811-125521
  * Faucher-Giguère et al. (2008) Faucher-Giguère, C.-A., Prochaska, J. X., Lidz, A., Hernquist, L., & Zaldarriaga, M. 2008, ApJ, 681, 831, doi: 10.1086/588648
  * Goulding et al. (2018) Goulding, A. D., Zakamska, N. L., Alexandroff, R. M., et al. 2018, ApJ, 856, 4, doi: 10.3847/1538-4357/aab040
  * Guaita et al. (2011) Guaita, L., Acquaviva, V., Padilla, N., et al. 2011, ApJ, 733, 114, doi: 10.1088/0004-637X/733/2/114
  * Haardt & Madau (2012) Haardt, F., & Madau, P. 2012, ApJ, 746, 125, doi: 10.1088/0004-637X/746/2/125
  * Haberzettl et al. (2012) Haberzettl, L., Williger, G., Lehnert, M. D., Nesvadba, N., & Davies, L. 2012, ApJ, 745, 96, doi: 10.1088/0004-637X/745/1/96
  * Hopkins et al. (2006) Hopkins, P. F., Hernquist, L., Cox, T. J., et al. 2006, ApJS, 163, 1, doi: 10.1086/499298
  * Hopkins et al. (2008) Hopkins, P. F., Hernquist, L., Cox, T. J., & Kereš, D. 2008, ApJS, 175, 356, doi: 10.1086/524362
  * Horowitz et al. (2022) Horowitz, B., Lee, K.-G., Ata, M., et al. 2022, ApJS, 263, 27, doi: 10.3847/1538-4365/ac982d
  * Ito et al. (2021) Ito, K., Kashikawa, N., Tanaka, M., et al. 2021, ApJ, 916, 35, doi: 10.3847/1538-4357/abfc50
  * Ito et al. (2022) Ito, K., Tanaka, M., Miyaji, T., et al. 2022, ApJ, 929, 53, doi: 10.3847/1538-4357/ac5aaf
  * Kashikawa et al. (2007) Kashikawa, N., Kitayama, T., Doi, M., et al. 2007, ApJ, 663, 765, doi: 10.1086/518410
  * Komberg et al. (1996) Komberg, B. V., Kravtsov, A. V., & Lukash, V. N. 1996, MNRAS, 282, 713, doi: 10.1093/mnras/282.3.713
  * Kusakabe et al. (2018) Kusakabe, H., Shimasaku, K., Ouchi, M., et al. 2018, PASJ, 70, 4, doi: 10.1093/pasj/psx148
  * Lee et al. (2014) Lee, K.-G., Hennawi, J. F., White, M., Croft, R. A. C., & Ozbek, M. 2014, ApJ, 788, 49, doi: 10.1088/0004-637X/788/1/49
  * Lee et al. (2012) Lee, K.-G., Suzuki, N., & Spergel, D. N. 2012, AJ, 143, 51, doi: 10.1088/0004-6256/143/2/51
  * Lee et al. (2016) Lee, K.-G., Hennawi, J. F., White, M., et al. 2016, ApJ, 817, 160, doi: 10.3847/0004-637X/817/2/160
  * Lee et al. (2018) Lee, K.-G., Krolewski, A., White, M., et al. 2018, The Astrophysical Journal Supplement Series, 237, 31
  * Liang et al. (2021) Liang, Y., Kashikawa, N., Cai, Z., et al. 2021, ApJ, 907, 3, doi: 10.3847/1538-4357/abcd93
  * Ly et al. (2009) Ly, C., Malkan, M. A., Treu, T., et al. 2009, ApJ, 697, 1410, doi: 10.1088/0004-637X/697/2/1410
  * Ma et al. (2024) Ma, K., Zhang, H., Cai, Z., et al. 2024, ApJ, 961, 102, doi: 10.3847/1538-4357/ad04da
  * Madau et al. (1998) Madau, P., Pozzetti, L., & Dickinson, M. 1998, ApJ, 498, 106, doi: 10.1086/305523
  * Meiksin (2009) Meiksin, A. A. 2009, Reviews of Modern Physics, 81, 1405, doi: 10.1103/RevModPhys.81.1405
  * Miyazaki et al. (2012) Miyazaki, S., Komiyama, Y., Nakaya, H., et al. 2012, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 8446, Ground-based and Airborne Instrumentation for Astronomy IV, ed. I. S. McLean, S. K. Ramsay, & H. Takami, 84460Z, doi: 10.1117/12.926844
  * Miyazaki et al. (2018) Miyazaki, S., Komiyama, Y., Kawanomoto, S., et al. 2018, PASJ, 70, S1, doi: 10.1093/pasj/psx063
  * Nadathur (2013) Nadathur, S. 2013, MNRAS, 434, 398, doi: 10.1093/mnras/stt1028
  * Newman et al. (2020) Newman, A. B., Rudie, G. C., Blanc, G. A., et al. 2020, ApJ, 891, 147, doi: 10.3847/1538-4357/ab75ee
  * Newman et al. (2022) —. 2022, Nature, 606, 475, doi: 10.1038/s41586-022-04681-6
  * Onoue et al. (2018) Onoue, M., Kashikawa, N., Uchiyama, H., et al. 2018, PASJ, 70, S31, doi: 10.1093/pasj/psx092
  * Ouchi et al. (2020) Ouchi, M., Ono, Y., & Shibuya, T. 2020, ARA&A, 58, 617, doi: 10.1146/annurev-astro-032620-021859
  * Overzier (2016) Overzier, R. A. 2016, A&A Rev., 24, 14, doi: 10.1007/s00159-016-0100-3
  * Park et al. (2015) Park, C., Song, H., Einasto, M., Lietzen, H., & Heinamaki, P. 2015, Journal of Korean Astronomical Society, 48, 75, doi: 10.5303/JKAS.2015.48.1.75
  * Pezzulli & Cantalupo (2019) Pezzulli, G., & Cantalupo, S. 2019, MNRAS, 486, 1489, doi: 10.1093/mnras/stz906
  * Pickles (1998) Pickles, A. J. 1998, PASP, 110, 863, doi: 10.1086/316197
  * Pilipenko & Malinovsky (2013) Pilipenko, S., & Malinovsky, A. 2013, arXiv e-prints, arXiv:1306.3970, doi: 10.48550/arXiv.1306.3970
  * Puccetti et al. (2020) Puccetti, S., Fiore, F., Bongiorno, A., et al. 2020, A&A, 634, A137, doi: 10.1051/0004-6361/201833601
  * Shen et al. (2011) Shen, Y., Richards, G. T., Strauss, M. A., et al. 2011, ApJS, 194, 45, doi: 10.1088/0067-0049/194/2/45
  * Shi et al. (2021) Shi, D. D., Cai, Z., Fan, X., et al. 2021, ApJ, 915, 32, doi: 10.3847/1538-4357/abfec0
  * Shimakawa et al. (2017) Shimakawa, R., Kodama, T., Hayashi, M., et al. 2017, MNRAS, 468, L21, doi: 10.1093/mnrasl/slx019
  * Shimakawa et al. (2023) Shimakawa, R., Pérez-Martínez, J. M., Koyama, Y., et al. 2023, arXiv e-prints, arXiv:2306.06392, doi: 10.48550/arXiv.2306.06392
  * Stroe et al. (2015) Stroe, A., Sobral, D., Dawson, W., et al. 2015, MNRAS, 450, 646, doi: 10.1093/mnras/stu2519
  * Sun et al. (2023) Sun, D., Mawatari, K., Ouchi, M., et al. 2023, ApJ, 951, 25, doi: 10.3847/1538-4357/accf88
  * Suzuki et al. (2005) Suzuki, N., Tytler, D., Kirkman, D., O’Meara, J. M., & Lubin, D. 2005, ApJ, 618, 592, doi: 10.1086/426062
  * Tucker et al. (1998) Tucker, W., Blanco, P., Rappoport, S., et al. 1998, ApJ, 496, L5, doi: 10.1086/311234
  * Tumlinson et al. (2017) Tumlinson, J., Peeples, M. S., & Werk, J. K. 2017, ARA&A, 55, 389, doi: 10.1146/annurev-astro-091916-055240
  * Uchiyama et al. (2018) Uchiyama, H., Toshikawa, J., Kashikawa, N., et al. 2018, PASJ, 70, S32, doi: 10.1093/pasj/psx112
  * Uchiyama et al. (2020) Uchiyama, H., Akiyama, M., Toshikawa, J., et al. 2020, ApJ, 905, 125, doi: 10.3847/1538-4357/abc47b
  * Webster (1982) Webster, A. 1982, MNRAS, 199, 683, doi: 10.1093/mnras/199.3.683
  * Wu & Shen (2022) Wu, Q., & Shen, Y. 2022, ApJS, 263, 42, doi: 10.3847/1538-4365/ac9ead
  * Wylezalek et al. (2013) Wylezalek, D., Galametz, A., Stern, D., et al. 2013, ApJ, 769, 79, doi: 10.1088/0004-637X/769/1/79
  * Zhang et al. (2023) Zhang, H., Cai, Z., Li, M., et al. 2023, arXiv e-prints, arXiv:2301.07359, doi: 10.48550/arXiv.2301.07359
  * Zhang et al. (2024) Zhang, H., Cai, Z., Liang, Y., et al. 2024, ApJ, 961, 63, doi: 10.3847/1538-4357/ad07d3
  * Zhang et al. (2022) Zhang, Y., Zheng, X. Z., Shi, D. D., et al. 2022, MNRAS, 512, 4893, doi: 10.1093/mnras/stac824