11institutetext: Kapteyn Astronomical Institute, University of Groningen,
Postbus 800, 9700 AV Groningen, The Netherlands 22institutetext: SRON
Netherlands Institute for Space Research, Landleven 12, 9747 AD, Groningen,
The Netherlands 33institutetext: Imperial College London, Blackett Laboratory,
Prince Consort Road, London, SW7 2AZ, UK 44institutetext: Department of
Physics and Astronomy, University of Hawai’i, 2505 Correa Road, Honolulu, HI
96822, USA 55institutetext: Institute for Astronomy, 2680 Woodlawn Drive,
University of Hawai’i, Honolulu, HI 96822, USA 66institutetext: Department of
Space and Astronautical Science, Graduate University for Advanced Studies,
SOKENDAI, Shonankokusaimura, Hayama, Miura District, Kanagawa 240-0193, Japan
77institutetext: Institute of Space and Astronautical Science, Japan Aerospace
Exploration Agency, 3-1-1 Yoshinodai, Chuo-ku, Sagamihara, Kanagawa 252-5210,
Japan

# Probing the Mpc-scale environment of hyperluminous infrared galaxies at
$2<z<4$

F. Gao Probing the Mpc-scale environment of hyperluminous infrared galaxies at
$2<z<4$Probing the Mpc-scale environment of hyperluminous infrared galaxies at
$2<z<4$Probing the Mpc-scale environment of hyperluminous infrared galaxies at
$2<z<4$Probing the Mpc-scale environment of hyperluminous infrared galaxies at
$2<z<4$ L. Wang Probing the Mpc-scale environment of hyperluminous infrared
galaxies at $2<z<4$Probing the Mpc-scale environment of hyperluminous infrared
galaxies at $2<z<4$Probing the Mpc-scale environment of hyperluminous infrared
galaxies at $2<z<4$Probing the Mpc-scale environment of hyperluminous infrared
galaxies at $2<z<4$ A. F. Ramos Padilla Probing the Mpc-scale environment of
hyperluminous infrared galaxies at $2<z<4$Probing the Mpc-scale environment of
hyperluminous infrared galaxies at $2<z<4$Probing the Mpc-scale environment of
hyperluminous infrared galaxies at $2<z<4$Probing the Mpc-scale environment of
hyperluminous infrared galaxies at $2<z<4$ D. Clements Probing the Mpc-scale
environment of hyperluminous infrared galaxies at $2<z<4$Probing the Mpc-scale
environment of hyperluminous infrared galaxies at $2<z<4$ D. Farrah Probing
the Mpc-scale environment of hyperluminous infrared galaxies at $2<z<4$Probing
the Mpc-scale environment of hyperluminous infrared galaxies at $2<z<4$Probing
the Mpc-scale environment of hyperluminous infrared galaxies at $2<z<4$Probing
the Mpc-scale environment of hyperluminous infrared galaxies at $2<z<4$ T.
Huang Probing the Mpc-scale environment of hyperluminous infrared galaxies at
$2<z<4$Probing the Mpc-scale environment of hyperluminous infrared galaxies at
$2<z<4$Probing the Mpc-scale environment of hyperluminous infrared galaxies at
$2<z<4$Probing the Mpc-scale environment of hyperluminous infrared galaxies at
$2<z<4$

###### Abstract

Context. Protoclusters are progenitors of galaxy clusters and are important
for studying how halo mass and stellar mass assemble in the early universe.
Finding signposts of such over-dense regions, for example bright dusty star-
forming galaxies (DSFG), is a popular method to identify protocluster
candidates.

Aims. Hyperluminous infrared galaxies (HLIRGs), which are ultra-massive and
show extreme levels of dusty star formation and/or black hole accretion, are
expected to reside in overdense regions with massive halos. We study the Mpc-
scale environment of the largest HLIRG sample to date (526 HLIRGs over 26
deg2) and investigate whether they predominantly live in overdense regions.

Methods. We first explore the surface density of Herschel 250 $\mu$m sources
around HLIRGs and compare with that around random positions. Then, we compare
the spatial distribution of neighbours around HLIRGs with that around randomly
selected galaxies using a deep IRAC-selected catalogue with good-quality
photometric redshifts. We also use a redshift-matched quasar sample and
submillimeter galaxy (SMG) sample to validate our method, as previous
clustering studies have measured the host halo masses of these populations.
Finally, we adopt a Friends of Friends (FOF) algorithm to seek (proto)clusters
that host HLIRGs.

Results. We find that HLIRGs tend to have more bright star-forming neighbours
(with 250 $\mu$m flux density $>$10 mJy) within 100$\arcsec$ projected radius
($\sim 0.8$ Mpc at $2<z<4$) than a random galaxy at a $3.7\sigma$
significance. In our 3D analysis, we find relatively weak excess of IRAC-
selected sources within 3 Mpc around HLIRGs compared with random galaxy
neighbours, mainly influenced by photometric redshift uncertainty and survey
depth. We find a more significant difference (at a 4.7$\sigma$ significance)
in the number of Low Frequency Array (LOFAR)-detected neighbours in the
deepest ELAIS-N1 (EN1) field. HLIRGs at $3<z<4$ show stronger excess compared
to HLIRGs at $2<z<3$ ($0.13\pm{0.04}$ and $0.14\pm{0.01}$ neighbours around
HLIRGs and random positions at $2<z<3$ respectively, and $0.08\pm{0.04}$ and
$0.05\pm{0.01}$ neighbours around HLIRGs and random positions at $3<z<4$
respectively), consistent with cosmic downsizing. Finally, we select and
present a list of 30 most promising protocluster candidates for future follow-
up observations.

###### Key Words.:

Galaxies: clusters: general – Galaxies: groups: general –Galaxies: high-
redshift

##  1 Introduction

In the local universe, the morphology-density relation is well established:
galaxies that reside in dense regions such as clusters are more likely to be
massive ellipticals and have lower star-formation rates (SFR), while galaxies
in less crowded regions are more likely to be less massive spirals with
ongoing star formation (Dressler 1980; Lewis et al. 2002; Gómez et al. 2003;
Kauffmann et al. 2004; Thomas et al. 2005; Peng et al. 2010; Cappellari et al.
2011; Scoville et al. 2013). Numerous studies have also revealed that massive
early-type galaxies formed most of their stellar mass in the early universe
(galaxy down-sizing, e.g., Cowie et al. 1996; Thomas et al. 2005; De Lucia et
al. 2006; Cimatti et al. 2008; Pérez-González et al. 2008; Ilbert et al.
2010). Therefore, overdense regions at the present day such as galaxy clusters
that host massive early-type galaxies must have had rapid star-formation
activity at higher redshifts, when clusters are less virialized (often
referred as protoclusters; Overzier 2016). For example, Elbaz et al. (2007)
used data from the Great Observatories Origins Deep Survey(GOODS) and
demonstrated that SFR increases as galaxy density increases at z$\sim 1$,
opposite to what have been found for local galaxies. By studying galaxy
members in a Spitzer-selected cluster at z = 1.62, Tran et al. (2010) found
that the fraction of star-forming galaxies (SFGs) triples from the lowest to
the highest density regions. A recent work by Lemaux et al. (2020) which
studied 6730 spectroscopically confirmed SFGs at $2<z<5$ found a positive
trend between average SFR and local environment.

Protoclusters are progenitors of today’s most massive structures. They play an
important role in answering many key questions, such as how halo and stellar
mass assemble across cosmic time, how environment affects the evolution of
massive galaxies and so on. To study the nature of protoclusters, we first
need an efficient method to find them. However, identifying protoclusters is
more challenging than identifying their descendant clusters, which can be
found by searching for overdensity of red massive galaxies or hot intracluster
medium (ICM) through the Sunyaev–Zel’Dovich effect. In contrast, the lack of
red massive galaxies and mature ICM make it difficult to trace protoclusters
(see review Overzier 2016, and references therein).

By definition protoclusters can be selected by tracing overdensity of galaxies
in a small region (on the scale corresponding to the expected size of massive
halos scale of several Mpc), ideally spectroscopically confirmed. However,
good-quality photometric redshifts have also been used (e.g., Capak et al.
2011; Chiang et al. 2014; Franck & McGaugh 2016; Lemaux et al. 2018). An
alternative method is to trace surface overdensity of galaxies that occupy a
narrow redshift slice. For example, Lee et al. (2014) reported a large-scale
structure containing three protoclusters at $z=3.78$ in the Bo$\rm\ddot{o}$tes
field by searching for overdensity of Lyman $\alpha$ emitter (LAE) candidates.
In their work, deep imaging survey revealed 65 LAEs within a comoving volume
of $72\times 72\times 25$ Mpc3. These protoclusters may evolve into clusters
with a total stellar mass of a few times of 10${}^{14}\,\mathrm{M_{\sun}}$ to
10${}^{15}\,\mathrm{M_{\sun}}$ at the present day. Toshikawa et al. (2018)
selected 179 protocluster candidates at $z\sim 4$ by tracing surface
overdensity of Lyman break galaxies (LBGs), based on data from the Hyper
SuprimeCam Subaru strategic program (HSC-SSP) covering a wide area of 121
deg2. These methods provide a simple and direct way to identify potential
protoclusters, and thus contribute to a population census of protoclusters.
However, they usually require a large-area sky survey to conduct a systematic
search for overdense regions due to their rarity. Moreover, as redshift
increases, protoclusters occupying a given volume will cover a larger sky
area, making it more challenging to find them in a survey with limited sky
coverage. In addition, selections based on overdensities of LAEs or LBGs may
not be able to uncover overdensities of dusty star-forming galaxies which
dominate the star-forming population at high redshifts (e.g., see Casey et al.
2014, for a review) due to severe obscuration in the optical bands. Owing to
these difficulties, identification of protoclusters at high redshifts is so
far largely dependent on serendipitous findings, by targeting signposts which
are expected to reside in dense regions. Popular signposts include radio-loud
galaxies (e.g., Pentericci et al. 2000; Kurk et al. 2001; Hayashi et al. 2012;
Shimakawa et al. 2014), enormous Lyman $\alpha$ nebulae (e.g., Arrigoni
Battaia et al. 2018; Cai et al. 2018; Li et al. 2021; Nowotka et al. 2022) and
quasars (e.g., Steidel et al. 2005; Balmaverde et al. 2017; García-Vergara et
al. 2021).

Given the fact that larger SFR is observed in denser region at high redshifts
and the predominant role of dusty star-formation activity at high redshifts,
it is widely adopted to search for protoclusters either through overdensities
of dusty star-forming galaxies (DSFGs) or using bright DSFGs as a signpost.
DSFG can be selected as bright far-infrared (FIR) or submillimeter detections.
Their star-formation activity is obscured by dust which absorbs UV/optical
photons and re-emits them in the longer wavelength range. One efficient method
to find overdensity of DSFGs is using large FIR/submillimeter surveys covering
a wide area such as surveys with Planck and Herschel. The total flux density
within the beam of a low-resolution instrument is usually a combination of
multiple sources located within the same beam. Follow-up high-resolution
imaging can help confirm the overdense nature for regions with exceedingly
large flux densities (e.g., Negrello et al. 2005). For example, by observing
Planck pre-selected protocluster candidates through Herschel, a number of
protoclusters have been confirmed (e.g., Clements et al. 2014; Flores-Cacho et
al. 2016; Greenslade et al. 2018; Martinache et al. 2018). A famous
protocluster selected in this way is SPT2349-56 at $z=4.3$ which was first
found in the 2,500 deg2 South Pole Telescope (SPT) submillimeter survey.
Follow-up Atacama Large Millimeter/submillimeter Array (ALMA) observations
reveal that this protocluster consists of at least 14 gaseous galaxies in a
region of 130 kpc diameter, and is very likely to eventually grow into one of
the most massive structures in the local universe (Miller et al. 2018; Hill et
al. 2020). This method of targeting bright sources detected in low-resolution
single dish surveys can be very effective at pinpointing the most intense
star-forming cores of protoclusters, which may be evolutionarily connected to
the bright cluster galaxies (BCGs) in the local Universe (e.g., Rotermund et
al. 2021).

The success of tracing protoclusters through DSFGs naturally leads us to ask
whether the most extreme DSFG live in the densest regions and trace the most
massive protoclusters. Hyper luminous infrared galaxies (HLIRGs) are extremely
luminous in the infrared (IR) band (rest-frame 8-1000 $\mu$m), reaching total
IR luminosity $L_{IR}>10^{13}L_{\odot}$ (Rowan-Robinson 2000). This huge
amount of IR luminosity implies a SFR of a few thousand
$\mathrm{M_{\sun}}\,\mathrm{yr^{-1}}$. In addition, active galactic nucleus
(AGN) activity can also contribute significantly to their extreme IR
luminosity. HLIRGs are expected to live in overdense environment as they are
extreme DSFGs and/or show vigorous black hole accretion. Therefore, they are
good targets to search for protoclusters at higher redshifts. To date, there
are only a few studies focused on the environment of HLIRGs due to their
rarity. Farrah et al. (2004) studied six HLIRGs selected using data from the
Infrared Astronomical Satellite (IRAS) at $0.44<z<1.55$ and evaluated their
environment through the amplitude of the spatial cross-correlation function,
based on a quantitative measurement presented in Longair & Seldner (1979).
They found a larger average clustering level of their HLIRG sample compared
with local IR luminous galaxies. Jones et al. (2014) selected 10 hot dust-
obscured galaxies (hot DOGs) that are detected in WISE 12 or 22 $\mu m$ but
undetected in 3.4 and 4.6 $\mu m$, six of which are HLIRGs. They found a
factor of $\sim$3 overdensity of SMGs around these hot DOGs compared with
blank-field submillimeter surveys. Assef et al. (2015) counted galaxies having
red IRAC colors within 1′of 90 HLIRGs selected from WISE. They found that the
number of red IRAC sources is significantly higher than those around random
pointings.

Following the largest HLIRG sample constructed in Wang et al. (2021b) and
further studied in Gao et al. (2021), in this work we investigate whether they
live in dense regions as expected. The structure of this paper is as follows.
We describe our data in Section 2. We first introduce our HLIRG sample, then
describe the deblended Herschel 250 $\mu$m catalogs and the deep IRAC-selected
photometric redshift catalogs used to search for HLIRG neighbours. We also
include a quasar sample and a SMG sample for cross-checks. We study the number
of HLIRG neighbours within different projected separations as well as spatial
volumes and compare with that around randomly chosen sources in Section 3. In
addition, we discuss the influence due to photometric redshift uncertainty and
survey depth. We also adopt a Friends of Friends (FOF) algorithm to search for
overdense regions that are associated with HLIRGs. We summarize our findings
and provide a list of 30 most promising protocluster candidates for future
observations in Section 5. Throughout this paper, we assume a flat
$\Lambda$CDM universe with $\Omega_{\text{M}}=0.286$ and $H_{0}=69.3\,\rm
kms^{-1}Mpc^{-1}$ (Nine-year Wilkinson Microwave Anisotropy Probe (WMAP)
results; Hinshaw et al. 2013). Unless otherwise stated, we adopt a Salpeter
(1955) initial mass function (IMF). For magnitudes, we use standard AB
magnitude system.

##  2 Data and Method

In this section, we introduce the various samples used in this work. We first
describe our HLIRG sample in Section 2.1. They are selected in Wang et al.
(2021b) and further studied in Gao et al. (2021). In order to search for
neighbours of HLIRGs and study their distributions, we need deep source
catalogs which trace the galaxy density field. Source catalogs without
redshift information can help find neighbours within certain projected
separations, and source catalogs containing redshift information can assist in
seeking neighbours within certain spatial volumes. We first introduce the
deblended Herschel source catalogs in Section 2.2 used to search for SFGs
around HLIRGs within different projected separations. We then describe the
deep IRAC-selected source catalogs and their good-quality photometric
redshifts in Section 2.3 (referred as deep photo-z catalogs hereafter) to find
HLIRG neighbours that reside within different spatial volumes. We characterize
the overdense nature of HLIRG environment by comparing the distributions of
HLIRG neighbours with that of random galaxy neighbours. In order to validate
our method, we make use of a quasar sample and an SMG sample which have
estimates of host dark matter halo mass through clustering studies. We briefly
summarize our samples in Section 2.6.

Figure 1: Distributions of stellar mass (upper), SFR (middle) and AGN
luminosity (bottom; only for HLIRGs with AGN fraction¿0.1) as a function of
redshift for all HLIRGs. The normalized histograms are also inserted. The
median uncertainty is indicated in the right bottom corner of each panel.

###  2.1 The HLIRG sample

Our HLIRG sample is first identified in Wang et al. (2021b), based on a parent
sample of IR luminous galaxies selected from the Herschel blind source
catalogs. These blind Herschel source catalogs are then cross-matched with the
Low Frequency Array (LOFAR) 150 MHz source catalog, specifically the LOFAR
Two-metre Sky Survey (LoTSS) Deep Fields First data release (Duncan et al.
2021; Kondapally et al. 2021; Sabater et al. 2021; Tasse et al. 2021) in three
deep fields, Bo$\rm\ddot{o}$tes, ELAIS-N1 (EN1) and Lockman-Hole (LH). The
cross-matching between Herschel and LOFAR catalogs is secure thanks to the
well-known tight correlation between the far-IR and radio (FIRC; e.g.,
Calistro Rivera et al. 2017; Wang et al. 2019). The angular resolution ($\sim
6\arcsec$), positional accuracy (0.2″) and the superb sensitivity (average 71
$\upmu$Jy beam-1) of LOFAR imaging help to associate the multi-wavelength
counterparts of these radio sources. Only 2-3% of LOFAR sources do not have
multi-wavelength counterparts in these three fields (see Kondapally et al.
2021). The extraction of multi-wavelength photometry in the LOFAR deep fields
and the process of cross-matching between the LOFAR radio sources and their
multi-wavelength counterparts are fully described in Kondapally et al. (2021).
As our blind Herschel sources have already been matched to the LOFAR radio
sources, it is then straightforward to find their multi-wavelength
counterparts.

We selected Herschel blind sources with 250 $\mu$m fluxes above 45, 35 and 40
mJy in Bo$\rm\ddot{o}$tes, EN1 and LH respectively. At these flux density cuts
$>90\%$ Herschel sources are matched with LOFAR sources. The majority of
Herschel sources are matched with only one LOFAR source (the unique sample)
and around 16% are matched with at least two LOFAR sources (the multiple
sample). We deblended the Herschel 250, 350 and 500 $\mu$m flux densities for
the multiple sample using XID+, a probabilistic de-blender tool (Hurley et al.
2017) and taking advantage of the aforementioned FIRC to calculate flux
density priors. For more details we refer the reader to Wang et al. (2021b).

In Wang et al. (2021b) we used the spectral energy distribution (SED) fitting
code Code Investigating GALaxy Emission (CIGALE Burgarella et al. 2005; Noll
et al. 2009) to obtain physical properties of these IR luminous galaxies such
as the total IR luminosity. We then further analyzed the HLIRGs (69, 198, and
259 HLIRGs in Bo$\rm\ddot{o}$tes, EN1, and LH respectively) in Gao et al.
(2021), adopting two different SED fitting codes CIGALE and CYprus Models for
Galaxies and their NUclear Spectral (CYGNUS; Efstathiou & Rowan-Robinson 1995;
Efstathiou et al. 2000; Efstathiou & Rowan-Robinson 2003; Efstathiou et al.
2013) and various AGN models (Efstathiou & Rowan-Robinson 1995; Fritz et al.
2006; Stalevski et al. 2012), to explore in detail how their physical
properties depend on model assumptions. The distributions of stellar mass
estimates, SFR estimates and AGN luminosity estimates as a function of
redshift is shown in Figure 1. We find that HLIRGs are extremely massive (with
a median stellar mass of $10^{12}\,\mathrm{M_{\sun}}$) and have a co-moving
volume density higher than what is expected from previous studies of the
global stellar mass functions. They also undergo active star-forming
activities with a median SFR of
$10^{3.3}\,\mathrm{M_{\sun}}\,\mathrm{yr^{-1}}$. Moreover, many HLIRGs are
associated with vigorous AGN activity. There are $30-50\%$ HLIRGs that have
AGN fraction $>0.3$ depending on which SED fitting code and AGN model being
used.

###  2.2 The deblended Herschel 250 $\mu$m source catalogs

To explore whether HLIRGs live in dense environments such as protoclusters
with intense on-going star-formation activity, we need deep catalogs of SFGs.
The Herschel Extragalactic Legacy Project (HELP; Shirley et al. 2019) is a
large program which combines and homogenises a wide range of multi-wavelength
surveys in fields observed by Herschel. HELP presents a catalog of around 170
million sources covering 1270 deg2 selected in optical-near IR (NIR)
wavelength range.

We use the deblended Herschel 250 $\mu$m source catalogs (DMU26) from HELP to
search for star-forming neighbours around HLIRGs. The blind Herschel source
catalogs we used to build the parent sample of IR luminous galaxies are
selected by finding peaks in the matched filtered (MF) maps (Chapin et al.
2011) and thus limited to bright sources. Therefore, they are not suitable as
tracers of the general SFG population. In contrast, the deblended catalogs are
generated by deblending the Herschel maps based on positions of known sources
in the 24 $\mu$m maps with higher resolution. They are therefore deeper than
the blind Herschel source catalogs and thus better for studying the
environment of HLIRGs. We use the deblended 250 $\mu$m source catalogs to
search for DSFGs around HLIRGs within different projected separations and
compare with those around random positions. In order to remove less reliable
detections of low significance, we apply a flux density cut of 10 mJy,
resulting in 42 406, 42 694 and 41 651 sources in Bo$\rm\ddot{o}$tes, EN1 and
LH respectively.

Table 1: Statistics of the HLIRGs, quasars, SMGs, the deblended 250 $\mu$m sources and deep IRAC-seletced photo-z catalogs in the three fields. Since all HLIRGs are all above $z=1$ (only one HLIRG is at $z<1$), we select sources at $z>1$ in all catalogs. No redshift cut to the HELP 250 $\mu$m deblended catalogs has been made as we only use them in the surface density analysis. field | limit area (deg2) | HLIRG | quasar | SMG | deep photo-z catalog | deblended 250 $\mu$m catalog
---|---|---|---|---|---|---
Bo$\rm\ddot{o}$tes | 5.06 | 31 | 235 | – | 150 261 | 19 025
EN1 | 5.75 | 85 | 272 | – | 161 845 | 16 255
LH | 6.5 | 86 | 346 | 93 | 130 907 | 14 593
Note: SMGs only locate in the LH field.

Figure 2: Sky distributions of the various samples used in this study. Black
dots are 1000 randomly selected sources from the deep IRAC-selected photo-z
catalogs at $2<z<3$ for clarity. We require all samples to be located inside
the central black boxes in order to avoid searching for neighbours outside of
the sky area coverage. HLIRGs at $2<z<3$, redshift-matched quasars and SMGs
(only in the LH field) are plotted as red circles, blue squares and cyan
squares respectively.

###  2.3 Deep photo-z source catalog

In order to find neighbours of HLIRGs in a certain spatial volume and
characterize their environment in 3D, we need deep source catalogs that
contain good-quality redshift information. SED fitting to the multi-wavelength
photometry can help determine galaxy properties such as stellar mass for these
sources.

We use deep source catalogs in the three deep fields compiled in Kondapally et
al. (2021) to search for neighbours around our HLIRG sample. The authors
generated deep source catalogs for the LH and EN1 fields and combined existing
catalogs for the Bo$\rm\ddot{o}$tes field. Source detection is carried out
using SExtractor (Bertin & Arnouts 1996) and multi-wavelength photometry is
measured first through forced, matched aperture then corrected to total fluxes
using aperture corrections across all 20 and 16 bands for EN1 and LH fields
respectively.

Photometric redshifts of the deep source catalogs are determined in Duncan et
al. (2021), adopting a hybrid method that combines template fitting and
machine learning. They achieve a good agreement when comparing with
spectroscopic redshifts, with a 1.6-2% scatter from spectroscopic redshifts
for galaxies and a 1.5-1.8% of outlier fractions (defined as
—$\Delta$z—/(1+spec-z)¿0.15).

We first apply quality flags in the three source catalogs to eliminate
unreliable and duplicate sources. As elaborated in Kondapally et al. (2021),
the FLAG_CLEAN flag indicates regions masked by bright stars and FLAG_OVERLAP
flag indicates multiwavelength coverage in multiple surveys. For
Bo$\rm\ddot{o}$tes, an additional flag FLAG_DEEP represents duplicate sources
in the $I$-band catalog of Brown et al. (2008). We then require the signal-to-
noise (S/N) ratio in the Spitzer IRAC 3.6 $\mu$m band to be above 3$\sigma$ to
further reduce spurious sources, reaching a limiting magnitude of 23.63. Since
our HLIRGs are all at high redshifts beyond $z=1$, 3.6 $\mu$m is close to
rest-frame NIR band whose emission comes mostly from evolved stellar
populations. Hence the 3.6 $\mu$m emission provides a good indicator of
stellar mass (e.g., Kauffmann & Charlot 1998; Cole et al. 2001) and it suffers
less dust obscuration. After applying the quality flags and S/N cuts, the
total number of galaxies above $z=1$ reduces to 331 032, 275 672, and 369 729
sources in the Bo$\rm\ddot{o}$tes, EN1 and LH fields respectively.

Table 2: The number of HLIRGs, redshift-matched quasars and SMGs summed over the three fields. We select the same number of quasars in each redshift bin as the number of HLIRGs. Numbers in parenthesis are the number of HLIRGs in the LH field. We select SMGs with half the number of HLIRGs due to small number statistics of the SMG sample. redshift bin | HLIRG | quasar | SMG (LH field)
---|---|---|---
2.0-2.2 | 13(7) | 13 | 3
2.2-2.4 | 18(8) | 18 | 4
2.4-2.6 | 19(8) | 19 | 4
2.6-2.8 | 22(6) | 22 | 3
2.8-3.0 | 15(7) | 15 | 3
3.0-3.2 | 8(6) | – | 3
3.2-3.4 | 11(5) | – | 2
3.4-3.6 | 11(6) | – | 3
3.6-3.8 | 13(4) | – | 2
3.8-4.0 | 8(2) | – | 1

###  2.4 Quasar sample

Quasars are among the brightest galaxies that are powered by accretion onto
the central super massive black holes (Salpeter 1964; Lynden-Bell 1969). Since
the advent of large sky surveys such as the Sloan Digital Sky Survey (SDSS;
York et al. 2000) and the 2dF Quasi-Stellar Object (QSO) redshift survey (2QZ;
Croom et al. 2004), it has long been well established that quasar clustering
becomes stronger as redshift increases (Porciani et al. 2004; Croom et al.
2005; Myers et al. 2006; Shen et al. 2007). They are expected to be hosted by
dark matter halos with characteristic halo masses around a few times of
$10^{12}\,\mathrm{M_{\sun}}$ with little to no dependence on quasar luminosity
(Croom et al. 2005; Myers et al. 2007; Shen et al. 2007; Ross et al. 2009;
White et al. 2012; Eftekharzadeh et al. 2015; Wang et al. 2015). According to
the well studied relationship between stellar mass and halo mass (SMHM
relationship; Guo et al. 2010; Moster et al. 2010; Behroozi et al. 2013; Wang
et al. 2013; Rodríguez-Puebla et al. 2015; Behroozi et al. 2019), halos around
a few times $10^{12}\,\mathrm{M_{\sun}}$ host galaxies with stellar masses
around a few times of $10^{10}\,\mathrm{M_{\sun}}$ over a wide range of
redshifts.

We use the latest SDSS sixteenth quasar catalog (Lyke et al. 2020) which
consists of the largest number of spectroscopically confirmed quasars to date,
with 750 414 quasars over 7 500 deg2, reaching a 99.8% completeness with a
contamination fraction of 0.3-1.3%. We make use of this quasar catalog due to
several reasons. First, estimates of the host dark matter halo mass for quasar
already exists thanks to extensive clustering analyses. Comparing the
neighbours around HLIRGs with those around quasars can give us an indication
of the environment of HLIRGs relative to the quasars. Moreover, these quasars
partly overlap with our HLIRG sample in the redshift distributions as later
described in Section 2.6. In summary, the brightness of quasars, good
understanding of their host halo environment, and overlapping redshift range
make them a good sample to serve as a cross-check.

Figure 3: Redshift distribution of our HLIRG (red), quasar (blue) and SMG
(cyan) samples in the three deep fields after location restriction. Most
HLIRGs are locacted at $z>2$. SMGs have a broadly similar redshift
distribution as HLIRGs, while quasars only partly overlap with the HLIRGs at
$2<z<3$. Due to small number statistics, we only study HLIRGs at $2<z<4$ and
use quasars at $2<z<3$ for a cross-check.

###  2.5 Sub-millimeter Galaxy sample

Sub-millimeter Galaxies (SMG) are rapidly star-forming galaxies at high
redshifts. Their rest-frame FIR emission from obscured star-forming activity
shifts towards longer wavelengths and make them detectable in the
(sub-)millimeter band. HLIRGs and SMGs are broadly similar in terms of their
power source (i.e., dusty star formation, and AGN activity in some cases)
although SMGs typically have IR luminosity more comparable to ultra luminous
IR galaxies (ULIRGs), reaching a few times of 10${}^{12}L_{\odot}$ (e.g.,
Chapman et al. 2005; Magnelli et al. 2012; Swinbank et al. 2014). In addition,
SMGs in general are biased towards colder dust temperatures compared to
Herschel-selected galaxies with similar IR luminosities (e.g., Swinbank et al.
2014; da Cunha et al. 2015).

Clustering measurements of SMGs report a strong clustering and typical host
halos with $M_{halo}\sim 10^{12}-10^{13}\,\mathrm{M_{\sun}}$ (Blain et al.
2004; Hickox et al. 2012; Chen et al. 2016; Wilkinson et al. 2017; Stach et
al. 2021). Such halos will host galaxies with stellar masses a few times of
$10^{10}\,\mathrm{M_{\sun}}$ based on the aforementioned SMHM relationship at
$2<z<4$ where the distribution of SMGs peaks. Bright SMGs can be used as a
signpost to trace high-redshift protocluster candidates and several studies
have successfully found SMGs that live in protoclusters (e.g., Smail et al.
2003; Daddi et al. 2009; Riechers et al. 2014; Cheng et al. 2019; Wang et al.
2021a). In our work, we also select a SMG sample to compare their neighbours
with random galaxy neighbours. Similar to the quasar sample described in
Section 2.4, by comparing HLIRG neighbours with SMG neighbours, we can find
out whether HLIRGs reside in similar or more massive halos. Moreover, SMGs
typically have a broadly similar redshift range (as later described in Section
2.6) and number distribution as HLIRGs as they both trace dusty star-forming
activity in the early universe.

The SMG sample we use is taken from Geach et al. (2017), which contains $\sim$
3000 SMGs above 3.5 $\sigma$ at 850 $\mu$m over $\sim 5$ deg2, as part of the
James Clerk Maxwell Telescope (JCMT) SCUBA-2 Cosmology Legacy Survey (S2CLS).
This survey reached an average 1 $\sigma$ depth of 1.2 mJy/beam, close to the
SCUBA-2 confusion limit of 0.8 mJy/beam. The catalog covers six fields in
total, but only the Lockman-Hole North (LH-N) field falls into the regions we
study in this paper. We cross-match the 850 $\mu$m source catalog with LOFAR
radio-optical cross-identified catalog in the LH field, in order to associate
SMG sources with their multi-wavelength counterparts. In contrast to the
quasar sample, the SMG sample suffers from small sample size and only has
photometric redshifts. Therefore, when comparing both SMG neighbours and
quasar neighbours with random galaxy neighbours, we can also have an idea of
the influence due to small number statistics and photometric redshift
uncertainty, given the fact that both samples have been well studied through
clustering analyses, and the quasar sample is larger and has spectroscopic
redshifts. The effect from photometric redshift uncertainty will be further
discussed in Section 4.2.

Figure 4: Upper: Number distributions of deblended 250 $\mu$m galaxies around
HLIRGs within 100, 200 and 400$\arcsec$ respectively. The black symbols and
errorbars are the median value and 16th/84th percentiles calculated from 100
realizations of random galaxies. Bottom: the cumulative distribution function
(CDF) of the number of HLIRG neighbours and random neighbours. HLIRGs on
average have more 250 $\mu$m neighbours than random positions. This difference
becomes weaker as radius increases, reducing from 3.7$\sigma$ within
100$\arcsec$ to 2.7 $\sigma$ within 200$\arcsec$ and 1.6$\sigma$ within
400$\arcsec$ respectively.

###  2.6 Summary

Our HLIRG sample is selected from the blind Herschel catalogs. To find star-
forming galaxies around HLIRGs we take advantage of the HELP 250 $\mu$m
deblended source catalogs. We also use a deep IRAC-selected source catalog
with good-quality photometric redshifts in order to study the 3D environment
of the HLIRGs within different spatial volumes. We also include a quasar
sample (at $2<z<3$) and a SMG sample (only in the LH field) as cross-checks.

We first restrict all samples to a limited sky area within the coverage of the
deblended Herschel 250 $\mu$m source catalogs and the deep photo-z catalogs,
in order to avoid searching for neighbours outside of the sky coverage. In
Figure 2,the sky areas covered by the deep photo-z catalogues are illustrated
using 1000 random sources over the redshift range 2¡z¡3 for clarity. We
confine all samples to be located in the central part of the regions shown by
the black boxes. The sky coverage of the deblended Herschel 250 $\mu$m source
catalogs is almost the same as that covered by the deep photo-z catalogs,
therefore we adopt the same limited zones. The sky area and the total number
of each sample in the three fields are listed in Table 1.

When exploring the star-forming neighbours of HLIRGs, we study all HLIRGs in
the limited zones and compare the number of their neighbours with the number
of random galaxy neighbours. To characterize the 3D environment of HLIRG
neighbours, we compare the distribution of HLIRG neighbours from the deep
photo-z catalogs with the distribution of redshift-matched random galaxy
neighbours. We also study the distribution of redshift-matched quasar
neighbours and SMG neighbours as cross-checks. After restriction in position
according to the sky coverage of the deep photo-z catalogs, the redshift
distributions of HLIRGs, quasars and SMGs are shown in Figure 3. Our HLIRG
sample lies mostly in $2<z<5$ and quasar samples are abundant below $z<3$. In
spite of a small number, SMGs occupy roughly the same redshift range as
HLIRGs. Based on the size of each sample, we restrict our analysis to HLIRGs
at $2<z<4$, quasars at $2<z<3$, and SMGs at $2<z<4$. To obtain a redshift-
matched comparison, we split our HLIRG sample into 10 redshift bins with 0.2
interval over $2<z<4$. For each redshift bin, we randomly select the same
number of quasars as the number of HLIRGs. The number of the selected SMGs is
half of the number of HLIRGs, due to their small sample size. The numbers of
each sample in each redshift bin are listed in Table 2.

Figure 5: The distribution of HLIRG neighbours (red empty circles), quasar
neighbours (blue squares) and random galaxies neighbours (black stars) as a
function of stellar mass. The red solid circles represent the subset of HLIRGs
which are the most promising protocluster candidates (see Section 5). Three
columns represent neighbours found within 3 Mpc, 6 Mpc and 10 Mpc
respectively. Two rows display results in two redshift bins. The dashed lines
are stellar mass completeness limits in each redshift bins. We combine Poisson
error, photometric redshift uncertainty and standard deviation in the three
fields for the HLIRG neighbours. We only consider sampling uncertainty for
quasar neighbours and random neighbours. We find no excess in quasar
neighbours at $2<z<3$ and a weak excess in HLIRG neighbours at $3<z<4$
compared with random neighbours. This excess disappear as searching radius
increases. For the most promising protocluster candidates, we observe an
enhanced excess signal.

Figure 6: Similar to Figure 5, but only in the LH filed. Neighbours of SMGs
are plotted as cyan squares. We only include Poisson error for the SMG
neighbours. We also find weak excess in the SMG neighbours and this excess
disappears as searching radius increases.

##  3 Results

###  3.1 The 2D surface distribution of 250 $\mu$m deblended star-formation
galaxies around HLIRGs

High-redshift protoclusters are expected to be overdense regions with vigorous
on-going star-formation activity. Therefore, it is reasonable to target
overdensities of highly star-forming galaxies in order to identify
protocluster candidates. This can be achieved by using high-resolution
millimeter/submillimeter surveys to search for overdensity of DSFGs around
protocluster candidates, which can be pre-selected as sources with exceedingly
large flux density in low-resolution surveys. For example, Planck
Collaboration et al. (2016) presented a total of 2151 high-redshift sources
whose flux density at 545 GHz are above 500 mJy and colors indicate $z>2$.
Optical to submillimeter observations reveal that only a small fraction are
strong gravitational lensed sources and the vast majority of them ($\sim
97\%$) are overdense regions of DSFGs. In order to investigate further the
nature of these potential protoclusters, follow-up observations have been
made. For example, SCUBA-2 observations at 850 $\mu$m (MacKenzie et al. 2017;
Cheng et al. 2019, 2020) have found a significant overdensity of 850 $\mu$m
sources compared with blank-field distributions and a rapid on-going star-
formation activity. Another example is observing unresolved bright sources
detected by SPT111Due to the smaller beam size of SPT (1′compared to 3′of
Planck), SPT-detected bright sources are mostly lensed galaxies., a single-
dish submillimeter telescope with a beam size of one arcminute (Carlstrom et
al. 2011), through interferometric facilities such as ALMA (e.g, Chapman et
al. 2020; Wang et al. 2021a).

We search for deblended 250 $\mu$m sources surrounding our HLIRGs within
different projected separations, and compare with that surrounding random
galaxies. Those which have more star-forming neighbours than random positions
are likely to be potential protoclusters and can be used as targets for future
observations to distinguish real protocluster members from chance projections.
We first study the total number of HLIRGs in the central confined region and
count the number of their neighbours within 100$\arcsec$. We visually inspect
them and remove a few HLIRGs due to their closeness to the boundary
respectively, resulting in 27, 77 and 74 HLIRGs in Bo$\rm\ddot{o}$tes, EN1 and
LH respectively. We randomly select the same number of random galaxies and
count their neighbours within 100$\arcsec$. These randoms galaxies are
required to be located at least 200$\arcsec$ away from any HLIRG, to make sure
that there is no overlap between HLIRG environments and random environments.
We run random selections 100 times and calculate the average distribution of
the number of their neighbours. We find that HLIRGs have a median value of
$13.0\pm{5.0}$ neighbours while random galaxies have a median value of
$11.4\pm{0.5}$ neighbours (calculated as the mean value and the standard
deviation of median number of neighbours from 100 realizations). HLIRGs live
in an overdense environment at a 3.2$\sigma$ significance level. As radius
increases, the difference in environment between HLIRGs and random positions
becomes weaker (with 2.9$\sigma$ significance level within 200$\arcsec$) or
even disappear (within 400$\arcsec$).

To reduce contaminants from low-redshift sources, we adopt the Herschel color
cut which requires $S_{350\mu m}/S_{250\mu m}>0.7$ and $S_{500\mu m}/S_{350\mu
m}>0.6$ to effectively select sources at $1.5<z<3$ (see Planck Collaboration
et al. 2016). The numbers of HLIRGs occupying this redshift range are 16, 39,
and 45 in Bo$\rm\ddot{o}$tes, EN1 and LH respectively, and reduce to 16, 34,
and 40 after visual inspection to remove objects near the boundary. In the
upper panel of Figure 4, we present the number distributions of star-forming
neighbours around each HLIRGs within 100$\arcsec(\sim$ 1 Mpc at $1.5<z<3$),
200$\arcsec$ and 400$\arcsec$ respectively. The black symbols and errorbars
are the median values and 16th/84th percentiles calculated from 100
realizations which select the same number of random galaxies and count their
neighbours. We also require no overlap between HLIRG neighbours and random
neighbours by requiring random galaxies locate at least 2 times of searching
radius away from every HLIRG. The cumulative distribution function (CDF) of
the number of HLIRG neighbours and random neighbours are shown in the lower
panel of Figure 4. We observe a clear difference in the number of star-forming
neighbours within 100$\arcsec$ as HLIRGs at $1.5<z<3$ tend to have more
neighbours, with a median value of $4.0\pm{2.3}$, compared with a median value
of $2.9\pm{0.3}$ for random positions (calculated as the mean value and the
standard deviation of median number of neighbours from 100 realizations). This
difference implies that HLIRGs live in overdense regions at a 3.7$\sigma$
significance level. This difference becomes weaker as the search radius
increases, reducing to 2.7 and 1.6 $\sigma$ for 200” and 400”, respectively.
We also adopt IRAC color cut which requires $m_{3.6\mu m,AB}-m_{4.5\mu
m,AB}>-0.1$ to select sources at $z>1.3$ (Papovich 2008). The similar
discrepancy in the number of star-forming neighbours to what has been found
using the Herschel color cut also supports our expectation in the overdense
nature of HLIRG environments.

###  3.2 Spatial distribution of IRAC-selected sources around HLIRGs

To probe the small-scale environment of HLIRGs in 3D, we search for their
neighbours using the deep IRAC-selected photo-z catalogs. In this section, we
search for neighbours around HLIRGs, and calculate the co-moving volume
density of these neighbours as a function of stellar masses. We also study the
neighbours around quasars and SMGs to validate our method and to explore how
they compare with HLIRG neighbours. We adopt three fixed radius, 3 Mpc, 6 Mpc
and 10 Mpc (typical scale of protoclusters, see e.g., Cai et al. 2016) to find
nearby galaxies.

Following Wang et al. (2021b) and Gao et al. (2021), we use CIGALE to estimate
stellar mass for the deep photo-z source catalogs. We use a delayed-$\tau$
plus a starburst star formation history (SFH) and Bruzual & Charlot (2003)
single stellar populations (SSPs), with a Salpeter (1955) IMF and solar
metallicity. We adopt a double power-law dust attenuation law based on Charlot
& Fall (2000) and dust emission models from Draine et al. (2014). We use the
Fritz et al. (2006) AGN models but include fewer parameters for simplicity. We
refer to Gao et al. (2021) for a detailed description of the parameter space.

Figure 5 shows the distribution of stellar mass estimates of HLIRG neighbours
as well as redshift-matched quasar neighbours and random galaxy neighbours.
Figure 6 shows the same distribution but only for HLIRGs in LH and redshift-
matched SMGs. Dashed lines represent the stellar mass completeness limit
calculated using the method described in Pozzetti et al. (2010). Briefly, we
calculate the limiting stellar mass $M_{limit}$ for each galaxy which is the
value a galaxy will have if its apparent 3.6 $\mu$m magnitude equals to the
limiting magnitude of the entire catalog. Then, we select the 20% of the
faintest galaxies in each redshift bin and assign the completeness limit as
the $M_{limit}$ value below which 90% of these faint galaxies lie. The stellar
masses completeness limit derived in this way is
$10^{10.2}\,\mathrm{M_{\sun}}$ and $10^{10.5}\,\mathrm{M_{\sun}}$ in $2<z<3$
and $3<z<4$ respectively.

In Figures 5 and 6, the three columns represent neighbours found within
different radius as indicated at the top. The two rows display two redshift
bins as indicated in the upper-left corner. For the uncertainty of the HLIRGs
neighbours, we combine the Poisson error, standard deviation in three fields
(for Figure 5 only) and uncertainty from photometric redshifts in quadrature.
The latter comes from the standard deviation from 100 realizations. In every
realization the photometric redshift of each HLIRG is randomly determined from
a Gaussian distribution with a mean value of its photometric redshift
$z$-photo and a standard deviation of (1+$z_{\rm photo})\times$ 0.15 (typical
boundaries used for counting outliers in photometric redshift fitting, see
Section 2.3). Data points and lower/upper boundaries for the quasar neighbours
and random neighbours are median value and 16th/84th percentiles respectively
of 100 realizations. In each realization we randomly selected redshift-matched
quasars and random positions to take uncertainty due to sampling into account.
We do not consider Poisson error, standard deviation in three fields or
uncertainty from photometric redshifts, because uncertainty due to random
sampling dominates. For the uncertainty of SMG neighbours, we only consider
Poisson errors due to their small sample size.

As Figures 5 and 6 show, there is a weak excess of HLIRG neighbours at $3<z<4$
compared with random positions. This excess of HLIRG neighbours become weaker
as searching radius increases, implying that large-scale environment of HLIRGs
is no different from random positions. The environment of quasars and SMGs
within different scales is similar or marginally denser than random positions.
We attribute the lack of strong excess when compared with random positions to
the selection criterion of the deep photo-z catalogs we use. As explained in
Section 2.4 and 2.5, quasars are typically hosted in dark matter halos around
a few times $10^{12}\,\mathrm{M_{\sun}}$, while SMGs are typically hosted in
slightly more massive halos with masses around $10^{13}\,\mathrm{M_{\sun}}$.
Such halos will host galaxies with stellar masses a few time of
$10^{10}\,\mathrm{M_{\sun}}$. The completeness limits of our deep photo-z
catalogs reach $10^{10.2}\,\mathrm{M_{\sun}}$ and
$10^{10.5}\,\mathrm{M_{\sun}}$ in $2<z<3$ and $3<z<4$ respectively after
applying a 3$\sigma$ cut at 3.6 $\mu$m. It is expected that the environments
of SMGs and quasars are similar to the environment of these massive random
positions. We will discuss further about the influence of survey depth in
Section 4.1. In addition, there are some studies pointing out that quasars are
not good beacons of overdense structures (e.g., Trainor & Steidel 2012;
Fanidakis et al. 2013; Bañados et al. 2013; Uchiyama et al. 2018).

In order to test whether the stellar mass estimates are reliable, we carry out
several SED fitting runs with different models and parameters. For example, we
tried different SFHs, including double power-laws and delayed-$\tau$ models
without a later burst, and dust attenuation models with different slopes. We
find these changes do not introduce significant systematic difference to the
stellar mass estimates. We also adopt a machine learning random forest (RF)
method to estimate stellar mass in Appendix A. When comparing HLIRG neighbours
with random galaxy neighbours as a distribution of RF derived stellar mass, we
find a similar picture to the one presented in Figures 5 and 6.

###  3.3 Search for clusters using Friend of Friend algorithm

In this section, we take advantage of the method used in Huang et al. (2021),
which discovered 88 cluster candidates that consist of 4390 member galaxies at
$z<1.1$ in the 5.4 deg2 AKARI North Ecliptic Pole (NEP) field, based on
photometric redshifts. The authors calculated the local density for each
galaxy through the angular separation to the 10th nearest neighbour (Santos et
al. 2021), then adopted a FOF algorithm to overdense galaxies in order to
determine cluster members. The FOF algorithm is improved in that it considers
neighbour number as a criterion for grouping friends so that the chain-shaped
linking is avoided. Also, using redshift-dependent linking length makes the
FOF algorithm valid at higher redshift.

We make use of this method to find (proto)cluster candidates in deep IRAC-
selected photo-z catalogs after both quality flag cuts and 3.6 $\mu$ m
detection significance cuts are applied (see Section 2.3) to reduce spurious
sources which could influence the search for overdense regions. The Normalized
Median Absolute Deviation (NMAD) of photometric redshift uncertainty are
0.022, 0.029, and 0.046 for Bo$\rm\ddot{o}$tes, EN1, and LH fields,
respectively. There are 73 974, 61 281, and 72 639 overdensities selected as
having local density greater than 2 in three fields respectively. Then we
apply the FOF algorithm on the overdensities. The projected linking length is
determined by fitting the distances to the tenth neighbour for all galaxies.
Having at least 10 friends within the linking length is required to continue
making new friends. Clusters are selected as having more than 30 members. We
find 52, 49 and 57 clusters above z ¿ 1 in Bo$\rm\ddot{o}$tes, EN1 and LH
fields respectively with 1, 2, 1 of them hosting HLIRGs respectively.

The majority of our HLIRG sample are not associated with (proto)cluster
candidates identified by FOF algorithm. One explanation could be due to the
missing overdense regions traced by the IRAC-selected photo-z catalogs. Dust-
obscured overdense regions (i.e, 250 $\mu$m bright overdense regions in
Section 3.1) may have low S/N ratios or be totally missed in the photo-z
catalogs. For example, Kubo et al. (2019) found that at $z\sim 3.8$,
protoclusters may host obscured AGNs missed by optical selection. Therefore,
the potential overdense environment which host HLIRGs cannot be identified by
the FOF algorithm. Another explanation could be that HLIRGs reside in various
environments and only a small fraction of them are hosted by protoclusters. We
will discuss further in Section 4.3.

Figure 7: The normalized 3.6 $\mu$m magnitudes distribution of the deep IRAC-
selected photo-z catalogs (black), HLIRG neighbours (red hatched), quasar
neighbours (blue) and SMG neighbours (cyan; only in the LH field) within 3
Mpc. We find that HLIRG neighbours are predominantly faint in the 3.6 $\mu$m
band.

Figure 8: Number distribution of LOFAR-detected neighbours around HLIRGs
within 100$\arcsec$ compared with that around random positions in
Bo$\rm\ddot{o}$tes and EN1 fields respectively. The minimum 150 MHz flux
densities are 0.15 and 0.07 mJy respectively. We find a more significant
difference (4.7$\sigma$ in comparison with 1.9$\sigma$) between HLIRG
neighbours and random neighbours in the EN1 field which is the deepest LOFAR
field.

##  4 Discussion

###  4.1 The effect of survey depth

We do not find a strong excess in the quasar neighbours and SMG neighbour
compared to random neighbours. The main explanation is the similar halo mass
range that host quasars, SMGs and random galaxies within our deep photo-z
catalogs. In contrast, our HLIRGs are ultra-massive and thus are expect to
reside in dark matter halos with masses in a range of a few times of
$10^{13}\,\mathrm{M_{\sun}}$ to $10^{14}\,\mathrm{M_{\sun}}$ at $2<z<3$ where
most our HLIRG lie. Such extremely massive halos should host a large number of
member galaxies. However, we only observe a relatively weak excess in HLIRG
spatial neighbours compared to random neighbours.

We attribute this lack of excess to the fact that some HLIRG neighbours are
not detected due to the limited survey depth. We select the deep photo-z
catalogs requiring 3-$\sigma$ detection in the 3.6$\mu$m band which is
sensitive to emission from evolved stars. Some star-forming neighbours of
HLIRGs may be too dusty to be selected. Therefore, they are faint and of low
significance in the 3.6$\mu$m band, resulting in non-detection in the deep
photo-z catalogs. In Figure 7 we display the normalized probabilistic
distribution of 3.6 $\mu$m magnitudes for all sources in the deep photo-z
catalogs, as well as neighbours of HLIRGs, quasars and SMGs (only in the LH
field). We find that HLIRG neighbours and SMG neighbours typically have
fainter 3.6 $\mu$m magnitudes. It is highly possible that some neighbours are
too faint in the 3.6 $\mu$m bands and are not included in our deep photo-z
catalogs after applying a 3$\sigma$ significance cut. We may need deeper
catalogs in order to find these fainter neighbours and study the environment
of HLIRGs.We anticipate that with deeper catalogs, more HLIRG neighbours will
be revealed relative to the random neighbours, leading to a stronger excess
signal in the HLIRG environments when comparing with random positions.

We also explore the effect due to survey depth by seeking LOFAR-detected
neighbours around HLIRGs. The sensitivity of LOFAR radio survey varies in
different fields, reaching 32, 20 and 22 $\upmu$Jy beam-1 in
Bo$\rm\ddot{o}$tes, EN1 and LH fields respectively. We search for LOFAR-
detected neighbours within 100$\arcsec$ around 64 and 198 HLIRGs in
Bo$\rm\ddot{o}$tes and EN1 respectively as these two fields are the most
different: After requiring a 3$\sigma$ cut in the signal-to-noise ratio of the
150 MHz detections and confining all sources to be located in the central 14
and 24 degree2 respectively, the minimum 150 MHz flux densities are 0.16 and
0.07 mJy respectively. We keep most of HLIRGs in the Bo$\rm\ddot{o}$tes field
and all HLIRGs in the EN1 field as LOFAR observations cover a large sky area.
In each field we follow the same procedure as described in Section 3.1,
searching for neighbours around HLIRGs and counting neighbours around random
positions. The distribution of the number of LOFAR-detected neighbours around
HLIRGs in two fields are shown in Figure 8.

As Figure 8 shows, the difference in the abundance between HLIRG neighbours
and random neighbours is larger in the EN1 field, at a 4.7$\sigma$
significance level compared to 1.9$\sigma$ significance level in the
Bo$\rm\ddot{o}$tes field. When requiring the same depth of LOFAR detections
with 150 MHz flux density $>0.5$ mJy, this discrepancy in different fields
disappears, reaching 1.0$\sigma$ and 0.9$\sigma$ significance levels
respectively. We hypothesize that deeper surveys can play an important role in
studying small-scale environments as they can include fainter and less massive
galaxies which trace more typical environments with lower host halo mass
values. The excess signal may become more significant when comparing HLIRG
environments with these less massive environments.

###  4.2 The need for spectroscopic follow-up

Another potential reason for the relatively weak excess of spatial HLIRG
neighbours compared with random neighbours is due to photometric redshift
uncertainty. At a given sky position, a HLIRG at $2<z<4$ will move into a new
location that is $70-150$ Mpc away which is beyond the scale of its host halo
when its redshift is shifted with a small value of 0.1. Also true neighbours
of HLIRGs which live in the same overdense structure could be moved outside of
the searching radius due to photometric redshift uncertainty.

It is challenging to obtain spectroscopic redshifts for HLIRGs due to severe
dust obscuration. So far, we have managed to obtain spectroscopic redshifts
for a few HLIRGs. First, there are 28 HLIRGs (6, 8, and 14 in
Bo$\rm\ddot{o}$tes, EN1, and LH respectively) that have spectroscopic
redshifts collected and compiled in the deep source catalogs (Duncan et al.
2021). Second, two HLIRGs, LH_16049 and LH_23905 are observed by Northern
Extended Millimeter Array (NOEMA) and their spectroscopic redshifts are
consistent with photometric redshifts. LH_16049 has a $z_{\rm phot}$ of 2.47
and a $z_{\rm spec}$ of 2.464. LH_23905 has a $z_{\rm phot}$ of 2.66 and a
$z_{\rm spec}$ of 2.558 (private communications with Sharon Chelsea and Axel
Weiss). These two HLIRGs also demonstrate the good quality of photometric
redshifts. We also obtained SUBARU observations for an extreme starburst
candidate EN1_19770. Details of the observation campaign and data processing
are described in the Appendix B. We retrieve a spectroscopic redshift of 3.74
which is in good agreement with the photometric redshift estimate of 3.91 in
Duncan et al. (2021).

###  4.3 Are HLIRGs good tracers of protoclusters?

We find that HLIRGs in general have more bright Herschel sources (with
$250\mu$m flux density above 10 mJy) within 100$\arcsec$ than random positions
at a $3.7\sigma$ significance level. However, when seeking 3D neighbours with
IRAC-selected photo-z catalog, we only find relatively weak excess compared to
random positions. Apart from the influence due to survey depth and photometric
redshift uncertainty, another potential explanation could be attributed to
HLIRG themselves. Are HLIRGs good tracers of overdense regions such as
protoclusters?

From Figure 4 we find that around 68-87% of our HLIRGs (11 of 16, 23 of 34, 35
of 40 in Bo$\rm\ddot{o}$tes, EN1 and LH fields respectively) have more
Herschel-detected neighbours within 100$\arcsec$ than the median number of
neighbours of random positions while the remaining 13-32% have fewer Herschel-
detected neighbours. We select 30 HLIRGs as the most promising protocluster
candidates based on the number of Herschel neighbours and FOF algorithm (see
Section 5) and plot the stellar mass distribution of their neighbours as red
solid circles in Figure 5. We observe a slightly enhanced excess compared to
the stellar mass distribution of neighbours of all HLIRGs. We conclude that
only some HLIRGs are more likely than others to reside in overdense regions,
both in 2D and 3D environments. We note here that this conclusion is reached
under the above-mentioned limitations of survey depth and photometric redshift
uncertainty. HLIRGs that show overdensity features under such limitations can
be used as the most promising candidates to identify protoclusters in follow-
up observations given that spectroscopic campaigns are time-consuming (see
next section). However, the possibility that HLIRGs showing no signs of
overdensity in the current study may actually reside in overdense regions
cannot be ruled out completely.

Miller et al. (2015) investigated a sample of mock SMGs from the Bolshoi
cosmological simulation (Riebe et al. 2013; Klypin et al. 2011; Rodríguez-
Puebla et al. 2016) and found that SMGs are incomplete tracers of the most-
massive structures as the most majority of the most-massive structures do not
host SMGs. Since the enhanced starbursting mode is short-lived, SMGs can only
be observed during a relatively short period. Their rarity results in some of
the most-massive structures hosting no SMGs. Moreover, due to the cosmic
downsizing, the most-massive structures are expected to cease star-formation
activity at an earlier epoch. Miller et al. (2015) also found that dark matter
halos at $z<2.5$ are less likely than their high-redshift counterparts to host
SMGs. This is broadly consistent with our findings that HLIRGs at $3<z<4$ show
a relatively stronger excess than HLIRGs at $2<z<3$ when comparing with random
positions.

Observationally,Chapman et al. (2015) carried out observational campaign for
one of the most luminous SMGs HS1700.850.1 at $z=2.82$ and concluded that it
resides in relatively voids. Laporte et al. (2015) studied the environment
around an extreme luminous IR galaxy HFLS3 at $z=6.34$ and found no evidence
for an overdensity of bright sources.222HFLS3 is at much higher redshift than
our HLIRG sample. Both Chapman et al. (2015) and Laporte et al. (2015)
searched for overdensity of LBGs which may not be a good tracer of dusty
overdensities at high redshifts. Similarly, Trainor & Steidel (2012) studied a
sample of 15 most luminous QSOs at $z\sim 2.7$. They found that these
extremely luminous QSOs reside in host halos with masses similar to their less
luminous counterparts. The possibility that extreme luminous SMGs, HLIRGs or
QSOs are merely rare events and do not necessarily trace the most massive
structures can not be ruled out.

HLIRGs maybe a more incomplete tracers as being rarer and more extreme than
typical SMGs. This is supported by the (proto)clusters found using FOF
algorithm. Of all (proto)clusters candidates identified, $>96\%$ of them do
not host one of our HLIRGs. Nonetheless, the reason why some HLIRGs do not
show an overdense nature of MPC-scale environment is beyond the scope of this
paper. Our work aims to study the environment of the largest HLIRG sample and
select the most promising protocluster candidates that benefit further
research.

##  5 Conclusions

Table 3: A list of 30 HLIRGs which have been selected as the most promising signposts of potential protoclusters. These candidates are ranked by possibility of residing in overdense regions, accounting for factors such as the number of Herschel-detected neighbours, spatial neighbours found in the deep photo-z catalogs, the quality of their redshifts (specscopically confirmed or the reduced $\chi^{2}$ in the fitting for photometric redshifts), whether found in FOF algorithm and the q value in FIRC (defined as log$(\frac{L_{IR}/(3.75\times 10^{12}Hz)}{L_{150MHz}})$; smaller values indicating more radio-loud). LOFAR id | field | RA | DEC | redshift | reduced $\chi^{2}$ for $z_{phot}$ | FOF | q value
---|---|---|---|---|---|---|---
spectroscopic redshift
21243 | LH | 162.33080 | 57.39226 | 1.68 | – | F | 1.63
19726 | LH | 162.53003 | 58.08506 | 2.84 | – | T | 1.68
EN1
40215 | EN1 | 242.58979 | 54.03271 | 2.29 | 0.27 | F | 1.28
37547 | EN1 | 242.83600 | 53.76636 | 2.90 | 0.26 | F | 1.34
38464 | EN1 | 242.75260 | 56.52186 | 2.50 | 0.37 | F | 1.75
34758 | EN1 | 243.09614 | 54.26354 | 2.45 | 0.48 | F | 2.16
54459 | EN1 | 241.18199 | 54.34068 | 2.69 | 0.32 | F | 1.49
44576 | EN1 | 242.15551 | 56.47998 | 2.45 | 0.43 | F | 1.63
60669 | EN1 | 240.40852 | 54.98878 | 2.40 | 0.48 | F | 1.72
40404 | EN1 | 242.57690 | 54.00540 | 2.59 | 0.13 | F | 1.75
50542 | EN1 | 241.61670 | 54.22248 | 4.92 | 0.32 | T | 1.65
46515 | EN1 | 242.00465 | 54.43105 | 3.69 | 0.23 | F | 1.57
23891 | EN1 | 244.16897 | 55.90912 | 4.11 | 1.54 | F | 0.94
39610 | EN1 | 242.64656 | 54.60312 | 4.72 | 0.55 | F | 1.98
35782 | EN1 | 243.00559 | 54.72401 | 4.02 | 0.17 | F | 1.71
LH
20490 | LH | 162.44004 | 58.13847 | 2.49 | 0.48 | F | 1.28
27801 | LH | 161.50650 | 59.87516 | 2.02 | 0.34 | F | 1.72
23594 | LH | 162.03957 | 56.36271 | 2.05 | 0.31 | F | 1.54
21892 | LH | 162.26817 | 58.46467 | 2.37 | 0.29 | F | 1.58
41471 | LH | 159.46355 | 59.17206 | 2.10 | 0.26 | F | 1.85
20655 | LH | 162.41359 | 57.91304 | 3.08 | 0.21 | F | 1.25
13561 | LH | 163.39872 | 58.22669 | 4.06 | 0.12 | F | 1.21
39722 | LH | 159.81428 | 58.54783 | 3.59 | 0.37 | F | 1.67
23157 | LH | 162.11348 | 58.56543 | 3.42 | 1.10 | F | 1.52
23101 | LH | 162.11968 | 58.42583 | 3.65 | 0.63 | F | 1.74
16853 | LH | 162.86970 | 56.62488 | 3.93 | 2.19 | T | 1.54
Bo$\rm\ddot{o}$tes
13499 | Bo$\rm\ddot{o}$tes | 218.59504 | 35.48697 | 2.18 | 1.22 | T | 1.31
26351 | Bo$\rm\ddot{o}$tes | 217.12640 | 35.52119 | 4.25 | 0.06 | F | 1.11
30484 | Bo$\rm\ddot{o}$tes | 216.57340 | 33.02774 | 3.83 | 0.06 | F | 1.52
13358 | Bo$\rm\ddot{o}$tes | 218.59538 | 33.50608 | 5.54 | 0.01 | F | 1.52
Note: EN1_23891 is 18.4$\arcmin$ away from the protocluster candidate PCCS1
857 G085.48+43.36 in Greenslade et al. (2018), in spatial volume of 7.78 Mpc
at $z=4.11$. Bo$\rm\ddot{o}$tes_13499 is 4.39$\arcmin$ away from the
protocluster candidate PLCKERC545 G060.36+66.56 in Greenslade et al. (2018),
in spatial volume of 2.23 Mpc at $z=2.18$.

In this work, we probe the Mpc-scale environment of the largest sample of
HLIRGs in three fields, in order to investigate whether they generally trace
overdense regions as expected. We first search for Herschel-detected star-
forming neighbours around HLIRGs and compare with that around random
positions. We apply a Herschel color cut to reduce chance projections. We find
that on average HLIRGs do indeed have more bright (with 250 $\mu$m flux
density $>10$ mJy) star-forming neighbours at a 3.7$\sigma$ significance level
with a median value of $4.0\pm{2.3}$ in their close vicinity (100$\arcsec$)
while the median number of random neighbours is $2.9\pm{0.3}$.

We then search for neighbours in 3D around HLIRGs using deep IRAC-selected
photo-z catalogs. We include a quasar sample and a SMG sample as cross-checks
because clustering analyses of these two samples have derived host dark matter
halo masses to be around around a few times of 10${}^{12}\,\mathrm{M_{\sun}}$
and 10${}^{13}\,\mathrm{M_{\sun}}$ respectively. They also (partly) overlap in
redshifts with the HLIRGs. We study the comoving volume density of HLIRGs
neighbours as a function of stellar mass and compare with that of random
neighbours. We find the environment of quasars and SMGs are similar or
marginally denser than the environment of random positions, mainly due to the
fact that random galaxies in the IRAC-selected photo-z catalogs live in
similar massive halos based on their stellar mass estimates.

We only find relatively weak excess of HLIRG 3D neighbours compared with
random neighbours. This is due to a number of factors such as the influence of
survey depth, photometric redshift uncertainty which dilutes the excess
signal, small number statistics, and so on. The IRAC-selected photo-z catalogs
may not be able to uncover all HLIRG neighbours as they are predominantly
faint at 3.6 $\mu$m. We investigate the influence due to survey depth by
seeking LOFAR-detected neighbours. In the deepest EN1 field, HLIRGs have more
LOFAR-detected neighbours at a 4.7 $\sigma$ significance level while in the
shallowest Bo$\rm\ddot{o}$tes field the excess signal is reduced to
1.9$\sigma$ significance level. After applying a higher flux density cut in
both fields, the significance levels drop to 1.0$\sigma$ and 0.9$\sigma$
respectively, which is in agreement with our expectation. In terms of the
influence of photometric redshift uncertainty, a small displacement (0.1 at
$2<z<4$) in redshift will result in a change of tens to hundreds of Mpc in
location where is far away from host halos and thus no neighbours may be
found. We expect that with spectroscopic redshifts, we would be able to
observe a much stronger excess of massive neighbours around our HLIRGs
compared with random positions. So far we have managed to assemble 31
spectroscopic redshifts for our HLIRG sample from literature and new
optical/submillimeter observations. We also find that HLIRGs at $3<z<4$ show
stronger excess when comparing with random positions than HLIRGs at $2<z<3$,
suggesting cosmic downsizing, as the most massive halos may stop star-forming
activity at $2<z<3$ and hence cannot be traced by extreme SFR.

Finally, we present a list of 30 HLIRGs in Table 3 that are ranked by the
degree of overdensity in its Mpc-scale environment. We also display the
cumulative distribution of the 250 $\mu$m flux densities of their neighbours
within 100$\arcsec$ as well as the overdensity parameters measured at
different radius in Figure 10. We take into account factors such as the
quality of their spectroscopic redshifts, the number of Herschel detected
neighbours, whether or not being associated with protocluster candidates found
in the FOF algorithm and so on. Two candidates are close to the protocluster
candidates identified in Greenslade et al. (2018). The authors cross-matched
Herschel catalogs to Planck compact sources and selected candidates that have
at least 3$\sigma$ overdensities in either 250, 350, or 500 $\mu$m sources.

Unlike previous studies based on bright unresolved sources detected in Planck
or SPT survey that have been pre-selected as protocluster candidates, our work
first select the most extreme dusty star-forming galaxies based on Herschel
blind survey and then investigate their environment. We find on average HLIRG
have more star-forming neighbours in their close vicinity. Those having more
bright neighbours are very likely to be hosted in extreme massive halos and
thus potentially live in overdense environments such as protoclusters. By
targeting these most promising protocluster candidates with follow-up
observations, we can spectroscopically confirm them and study their physical
properties such as the host halo mass, total SFR, gas mass, gas depletion time
and star-formation efficiency to probe the evolution path of massive halos in
an early universe.

## Acknowledgement

We thank Scott Trager for his help with the Subaru FOCAS spectrum analysis for
the HLIRG EN1_19770. We thank Matthieu Bethermin and Giulia Rodighiero for
their enlightening comments and suggestions.

## References

  * Arrigoni Battaia et al. (2018) Arrigoni Battaia, F., Chen, C.-C., Fumagalli, M., et al. 2018, A&A, 620, A202
  * Assef et al. (2015) Assef, R. J., Eisenhardt, P. R. M., Stern, D., et al. 2015, ApJ, 804, 27
  * Bañados et al. (2013) Bañados, E., Venemans, B., Walter, F., et al. 2013, ApJ, 773, 178
  * Balmaverde et al. (2017) Balmaverde, B., Gilli, R., Mignoli, M., et al. 2017, A&A, 606, A23
  * Behroozi et al. (2019) Behroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. 2019, MNRAS, 488, 3143
  * Behroozi et al. (2013) Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57
  * Bertin & Arnouts (1996) Bertin, E. & Arnouts, S. 1996, A&AS, 117, 393
  * Blain et al. (2004) Blain, A. W., Chapman, S. C., Smail, I., & Ivison, R. 2004, ApJ, 611, 725
  * Brown et al. (2008) Brown, M. J. I., Zheng, Z., White, M., et al. 2008, ApJ, 682, 937
  * Bruzual & Charlot (2003) Bruzual, G. & Charlot, S. 2003, MNRAS, 344, 1000
  * Burgarella et al. (2005) Burgarella, D., Buat, V., & Iglesias-Páramo, J. 2005, MNRAS, 360, 1413
  * Cai et al. (2016) Cai, Z., Fan, X., Peirani, S., et al. 2016, ApJ, 833, 135
  * Cai et al. (2018) Cai, Z., Hamden, E., Matuszewski, M., et al. 2018, ApJ, 861, L3
  * Calistro Rivera et al. (2017) Calistro Rivera, G., Williams, W. L., Hardcastle, M. J., et al. 2017, MNRAS, 469, 3468
  * Capak et al. (2011) Capak, P. L., Riechers, D., Scoville, N. Z., et al. 2011, Nature, 470, 233
  * Cappellari et al. (2011) Cappellari, M., Emsellem, E., Krajnović, D., et al. 2011, MNRAS, 416, 1680
  * Carlstrom et al. (2011) Carlstrom, J. E., Ade, P. A. R., Aird, K. A., et al. 2011, PASP, 123, 568
  * Casey et al. (2014) Casey, C. M., Narayanan, D., & Cooray, A. 2014, Phys. Rep, 541, 45
  * Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763
  * Chapin et al. (2011) Chapin, E. L., Chapman, S. C., Coppin, K. E., et al. 2011, MNRAS, 411, 505
  * Chapman et al. (2020) Chapman, S. C., Aravena, M. A., Ashby, M. L. N., et al. 2020, A massive protocluster at z=7 selected by the South Pole Telescope, HST Proposal
  * Chapman et al. (2015) Chapman, S. C., Bertoldi, F., Smail, I., et al. 2015, MNRAS, 453, 951
  * Chapman et al. (2005) Chapman, S. C., Blain, A. W., Smail, I., & Ivison, R. J. 2005, ApJ, 622, 772
  * Charlot & Fall (2000) Charlot, S. & Fall, S. M. 2000, ApJ, 539, 718
  * Chen et al. (2016) Chen, C.-C., Smail, I., Ivison, R. J., et al. 2016, ApJ, 820, 82
  * Cheng et al. (2019) Cheng, T., Clements, D. L., Greenslade, J., et al. 2019, MNRAS, 490, 3840
  * Cheng et al. (2020) Cheng, T., Clements, D. L., Greenslade, J., et al. 2020, MNRAS, 494, 5985
  * Chiang et al. (2014) Chiang, Y.-K., Overzier, R., & Gebhardt, K. 2014, ApJ, 782, L3
  * Cimatti et al. (2008) Cimatti, A., Cassata, P., Pozzetti, L., et al. 2008, A&A, 482, 21
  * Clements et al. (2014) Clements, D. L., Braglia, F. G., Hyde, A. K., et al. 2014, MNRAS, 439, 1193
  * Cole et al. (2001) Cole, S., Norberg, P., Baugh, C. M., et al. 2001, MNRAS, 326, 255
  * Cowie et al. (1996) Cowie, L. L., Songaila, A., Hu, E. M., & Cohen, J. G. 1996, AJ, 112, 839
  * Croom et al. (2005) Croom, S. M., Boyle, B. J., Shanks, T., et al. 2005, MNRAS, 356, 415
  * Croom et al. (2004) Croom, S. M., Smith, R. J., Boyle, B. J., et al. 2004, MNRAS, 349, 1397
  * da Cunha et al. (2015) da Cunha, E., Walter, F., Smail, I. R., et al. 2015, ApJ, 806, 110
  * Daddi et al. (2009) Daddi, E., Dannerbauer, H., Stern, D., et al. 2009, ApJ, 694, 1517
  * De Lucia et al. (2006) De Lucia, G., Springel, V., White, S. D. M., Croton, D., & Kauffmann, G. 2006, MNRAS, 366, 499
  * Draine et al. (2014) Draine, B. T., Aniano, G., Krause, O., et al. 2014, ApJ, 780, 172
  * Dressler (1980) Dressler, A. 1980, ApJ, 236, 351
  * Duncan et al. (2021) Duncan, K. J., Kondapally, R., Brown, M. J. I., et al. 2021, A&A, 648, A4
  * Efstathiou et al. (2013) Efstathiou, A., Christopher, N., Verma, A., & Siebenmorgen, R. 2013, MNRAS, 436, 1873
  * Efstathiou & Rowan-Robinson (1995) Efstathiou, A. & Rowan-Robinson, M. 1995, MNRAS, 273, 649
  * Efstathiou & Rowan-Robinson (2003) Efstathiou, A. & Rowan-Robinson, M. 2003, MNRAS, 343, 322
  * Efstathiou et al. (2000) Efstathiou, A., Rowan-Robinson, M., & Siebenmorgen, R. 2000, MNRAS, 313, 734
  * Eftekharzadeh et al. (2015) Eftekharzadeh, S., Myers, A. D., White, M., et al. 2015, MNRAS, 453, 2779
  * Elbaz et al. (2007) Elbaz, D., Daddi, E., Le Borgne, D., et al. 2007, A&A, 468, 33
  * Fanidakis et al. (2013) Fanidakis, N., Macciò, A. V., Baugh, C. M., Lacey, C. G., & Frenk, C. S. 2013, MNRAS, 436, 315
  * Farrah et al. (2004) Farrah, D., Geach, J., Fox, M., et al. 2004, MNRAS, 349, 518
  * Flores-Cacho et al. (2016) Flores-Cacho, I., Pierini, D., Soucail, G., et al. 2016, A&A, 585, A54
  * Franck & McGaugh (2016) Franck, J. R. & McGaugh, S. S. 2016, ApJ, 833, 15
  * Fritz et al. (2006) Fritz, J., Franceschini, A., & Hatziminaoglou, E. 2006, MNRAS, 366, 767
  * Gao et al. (2021) Gao, F., Wang, L., Efstathiou, A., et al. 2021, A&A, 654, A117
  * García-Vergara et al. (2021) García-Vergara, C., Rybak, M., Hodge, J., et al. 2021, arXiv e-prints, arXiv:2109.09754
  * Geach et al. (2017) Geach, J. E., Dunlop, J. S., Halpern, M., et al. 2017, MNRAS, 465, 1789
  * Gómez et al. (2003) Gómez, P. L., Nichol, R. C., Miller, C. J., et al. 2003, ApJ, 584, 210
  * Greenslade et al. (2018) Greenslade, J., Clements, D. L., Cheng, T., et al. 2018, MNRAS, 476, 3336
  * Guo et al. (2010) Guo, Q., White, S., Li, C., & Boylan-Kolchin, M. 2010, MNRAS, 404, 1111
  * Hayashi et al. (2012) Hayashi, M., Kodama, T., Tadaki, K.-i., Koyama, Y., & Tanaka, I. 2012, ApJ, 757, 15
  * Hickox et al. (2012) Hickox, R. C., Wardlow, J. L., Smail, I., et al. 2012, MNRAS, 421, 284
  * Hill et al. (2020) Hill, R., Chapman, S., Scott, D., et al. 2020, MNRAS, 495, 3124
  * Hinshaw et al. (2013) Hinshaw, G., Larson, D., Komatsu, E., et al. 2013, ApJS, 208, 19
  * Huang et al. (2021) Huang, T.-C., Matsuhara, H., Goto, T., et al. 2021, MNRAS, 506, 6063
  * Hurley et al. (2017) Hurley, P. D., Oliver, S., Betancourt, M., et al. 2017, MNRAS, 464, 885
  * Ilbert et al. (2010) Ilbert, O., Salvato, M., Le Floc’h, E., et al. 2010, ApJ, 709, 644
  * Jones et al. (2014) Jones, S. F., Blain, A. W., Stern, D., et al. 2014, MNRAS, 443, 146
  * Kauffmann & Charlot (1998) Kauffmann, G. & Charlot, S. 1998, MNRAS, 297, L23
  * Kauffmann et al. (2004) Kauffmann, G., White, S. D. M., Heckman, T. M., et al. 2004, MNRAS, 353, 713
  * Klypin et al. (2011) Klypin, A. A., Trujillo-Gomez, S., & Primack, J. 2011, ApJ, 740, 102
  * Kondapally et al. (2021) Kondapally, R., Best, P. N., Hardcastle, M. J., et al. 2021, A&A, 648, A3
  * Kubo et al. (2019) Kubo, M., Toshikawa, J., Kashikawa, N., et al. 2019, ApJ, 887, 214
  * Kurk et al. (2001) Kurk, J. D., Pentericci, L., Röttgering, H. J. A., & Miley, G. K. 2001, Astrophysics and Space Science Supplement, 277, 543
  * Laigle et al. (2016) Laigle, C., McCracken, H. J., Ilbert, O., et al. 2016, ApJS, 224, 24
  * Laporte et al. (2015) Laporte, N., Pérez-Fournon, I., Calanog, J. A., et al. 2015, ApJ, 810, 130
  * Lee et al. (2014) Lee, K.-S., Dey, A., Hong, S., et al. 2014, ApJ, 796, 126
  * Lemaux et al. (2020) Lemaux, B. C., Cucciati, O., Le Fèvre, O., et al. 2020, arXiv e-prints, arXiv:2009.03324
  * Lemaux et al. (2018) Lemaux, B. C., Le Fèvre, O., Cucciati, O., et al. 2018, A&A, 615, A77
  * Lewis et al. (2002) Lewis, I., Balogh, M., De Propris, R., et al. 2002, MNRAS, 334, 673
  * Li et al. (2021) Li, Q., Wang, R., Dannerbauer, H., et al. 2021, ApJ, 922, 236
  * Longair & Seldner (1979) Longair, M. S. & Seldner, M. 1979, MNRAS, 189, 433
  * Lyke et al. (2020) Lyke, B. W., Higley, A. N., McLane, J. N., et al. 2020, ApJS, 250, 8
  * Lynden-Bell (1969) Lynden-Bell, D. 1969, Nature, 223, 690
  * MacKenzie et al. (2017) MacKenzie, T. P., Scott, D., Bianconi, M., et al. 2017, MNRAS, 468, 4006
  * Magnelli et al. (2012) Magnelli, B., Lutz, D., Santini, P., et al. 2012, A&A, 539, A155
  * Martinache et al. (2018) Martinache, C., Rettura, A., Dole, H., et al. 2018, A&A, 620, A198
  * Miller et al. (2018) Miller, T. B., Chapman, S. C., Aravena, M., et al. 2018, Nature, 556, 469
  * Miller et al. (2015) Miller, T. B., Hayward, C. C., Chapman, S. C., & Behroozi, P. S. 2015, MNRAS, 452, 878
  * Moster et al. (2010) Moster, B. P., Somerville, R. S., Maulbetsch, C., et al. 2010, ApJ, 710, 903
  * Myers et al. (2007) Myers, A. D., Brunner, R. J., Nichol, R. C., et al. 2007, ApJ, 658, 85
  * Myers et al. (2006) Myers, A. D., Brunner, R. J., Richards, G. T., et al. 2006, ApJ, 638, 622
  * Negrello et al. (2005) Negrello, M., González-Nuevo, J., Magliocchetti, M., et al. 2005, MNRAS, 358, 869
  * Noll et al. (2009) Noll, S., Burgarella, D., Giovannoli, E., et al. 2009, A&A, 507, 1793
  * Nowotka et al. (2022) Nowotka, M., Chen, C.-C., Battaia, F. A., et al. 2022, A&A, 658, A77
  * Overzier (2016) Overzier, R. A. 2016, A&A Rev., 24, 14
  * Papovich (2008) Papovich, C. 2008, ApJ, 676, 206
  * Peng et al. (2010) Peng, Y.-j., Lilly, S. J., Kovač, K., et al. 2010, ApJ, 721, 193
  * Pentericci et al. (2000) Pentericci, L., Kurk, J. D., Röttgering, H. J. A., et al. 2000, A&A, 361, L25
  * Pérez-González et al. (2008) Pérez-González, P. G., Rieke, G. H., Villar, V., et al. 2008, ApJ, 675, 234
  * Planck Collaboration et al. (2016) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 596, A100
  * Porciani et al. (2004) Porciani, C., Magliocchetti, M., & Norberg, P. 2004, MNRAS, 355, 1010
  * Pozzetti et al. (2010) Pozzetti, L., Bolzonella, M., Zucca, E., et al. 2010, A&A, 523, A13
  * Riebe et al. (2013) Riebe, K., Partl, A. M., Enke, H., et al. 2013, Astronomische Nachrichten, 334, 691
  * Riechers et al. (2014) Riechers, D. A., Carilli, C. L., Capak, P. L., et al. 2014, ApJ, 796, 84
  * Rodríguez-Puebla et al. (2015) Rodríguez-Puebla, A., Avila-Reese, V., Yang, X., et al. 2015, ApJ, 799, 130
  * Rodríguez-Puebla et al. (2016) Rodríguez-Puebla, A., Behroozi, P., Primack, J., et al. 2016, MNRAS, 462, 893
  * Ross et al. (2009) Ross, N. P., Shen, Y., Strauss, M. A., et al. 2009, ApJ, 697, 1634
  * Rotermund et al. (2021) Rotermund, K. M., Chapman, S. C., Phadke, K. A., et al. 2021, MNRAS, 502, 1797
  * Rowan-Robinson (2000) Rowan-Robinson, M. 2000, MNRAS, 316, 885
  * Sabater et al. (2021) Sabater, J., Best, P. N., Tasse, C., et al. 2021, A&A, 648, A2
  * Salpeter (1955) Salpeter, E. E. 1955, ApJ, 121, 161
  * Salpeter (1964) Salpeter, E. E. 1964, ApJ, 140, 796
  * Santos et al. (2021) Santos, D. J. D., Goto, T., Kim, S. J., et al. 2021, MNRAS, 507, 3070
  * Scoville et al. (2013) Scoville, N., Arnouts, S., Aussel, H., et al. 2013, ApJS, 206, 3
  * Shen et al. (2007) Shen, Y., Strauss, M. A., Oguri, M., et al. 2007, AJ, 133, 2222
  * Shimakawa et al. (2014) Shimakawa, R., Kodama, T., Tadaki, K. I., et al. 2014, MNRAS, 441, L1
  * Shirley et al. (2019) Shirley, R., Roehlly, Y., Hurley, P. D., et al. 2019, MNRAS, 490, 634
  * Smail et al. (2003) Smail, I., Ivison, R. J., Gilbank, D. G., et al. 2003, ApJ, 583, 551
  * Stach et al. (2021) Stach, S. M., Smail, I., Amvrosiadis, A., et al. 2021, MNRAS, 504, 172
  * Stalevski et al. (2012) Stalevski, M., Fritz, J., Baes, M., Nakos, T., & Popović, L. Č. 2012, MNRAS, 420, 2756
  * Steidel et al. (2005) Steidel, C. C., Adelberger, K. L., Shapley, A. E., et al. 2005, ApJ, 626, 44
  * Swinbank et al. (2014) Swinbank, A. M., Simpson, J. M., Smail, I., et al. 2014, MNRAS, 438, 1267
  * Tasse et al. (2021) Tasse, C., Shimwell, T., Hardcastle, M. J., et al. 2021, A&A, 648, A1
  * Thomas et al. (2005) Thomas, D., Maraston, C., Bender, R., & Mendes de Oliveira, C. 2005, ApJ, 621, 673
  * Toshikawa et al. (2018) Toshikawa, J., Uchiyama, H., Kashikawa, N., et al. 2018, PASJ, 70, S12
  * Trainor & Steidel (2012) Trainor, R. F. & Steidel, C. C. 2012, ApJ, 752, 39
  * Tran et al. (2010) Tran, K.-V. H., Papovich, C., Saintonge, A., et al. 2010, ApJ, 719, L126
  * Uchiyama et al. (2018) Uchiyama, H., Toshikawa, J., Kashikawa, N., et al. 2018, PASJ, 70, S32
  * Wang et al. (2021a) Wang, G. C. P., Hill, R., Chapman, S. C., et al. 2021a, MNRAS, 508, 3754
  * Wang et al. (2013) Wang, L., Farrah, D., Oliver, S. J., et al. 2013, MNRAS, 431, 648
  * Wang et al. (2021b) Wang, L., Gao, F., Best, P. N., et al. 2021b, A&A, 648, A8
  * Wang et al. (2019) Wang, L., Gao, F., Duncan, K. J., et al. 2019, A&A, 631, A109
  * Wang et al. (2015) Wang, L., Viero, M., Ross, N. P., et al. 2015, MNRAS, 449, 4476
  * White et al. (2012) White, M., Myers, A. D., Ross, N. P., et al. 2012, MNRAS, 424, 933
  * Wilkinson et al. (2017) Wilkinson, A., Almaini, O., Chen, C.-C., et al. 2017, MNRAS, 464, 1380
  * York et al. (2000) York, D. G., Adelman, J., Anderson, John E., J., et al. 2000, AJ, 120, 1579

##  Appendix A Estimate stellar mass using Random Forest

To explore if there is any dependence on the method by which stellar masses
are estimated, we repeat the analysis on the distribution of stellar mass of
neighbours with an alternative method to derive stellar masses for the IRAC-
selected sources. We adopt the Random Forest (RF) method to derive stellar
mass estimates in EN1 and trained a RF estimator with COSMOS2015 data (Laigle
et al. 2016). We use the RandomForestRegressor included in Python package
scikit-learn to fit the decision trees in the Random Forest, which creates a
regression model that can be applied to the EN1 field data. We first selected
the COSMOS2015 galaxies to have $>3-\sigma$ detection in Spitzer IRAC1 band at
3.6 $\mu$m, which is the same criterion we use to reduce spurious sources in
Section 2.2. Galaxies with stellar masses below $10^{7.5}\,\mathrm{M_{\sun}}$
are discarded as they may provide us with incorrect low values. We divided the
COSMOS2015 data into three redshift-matched sub-samples: train, test and
validation, containing 180 931, 54 240 and 488 281 galaxies respectively. We
use seven overlapping bands between COSMOS2015 catalog and EN1 source catalog
($u$ band data from Canada-France Hawaii Telescope; $i,r,z$ and $y$ from
SUBARU; 3.6 $\mu$m and 4.5 $\mu$m from Spitzer) together with the photometric
redshift information, as the features for the regression model.

We use different combinations of values for the parameters max_depth (maximum
depth of the three), max_features (number of features to be consider for the
best split) and n_extimators (number of trees in the forest), inside
RandomForestRegressor to select the best parameters in order to train the
data. We start with a random search to narrow down the range and then we use
the grid search to select the best parameters. We found that the best values
for the RandomForestRegressor were max_depth = 31, max_features =
log2(n_features) and n_extimators = 300. This led to an accuracy of $\sim$
99.08% and average error of $\sim$ 0.084 dex for the stellar mass in both the
test and validation samples. With the trained regression model, we estimate
the stellar masses for 275,672 galaxies in the EN1 field (see Section 2.2). We
compare them with our stellar mass estimates using CIGALE SED fitting code,
finding a mean difference and standard deviation of $0.41\pm{0.27}$ dex.
However, COSMOS2015 adopted the Chabrier (2003) IMF while we use the Salpeter
(1955) IMF in this paper. Considering this difference, the mean difference
reduces to $\sim$0.18 dex. We also compare the RF stellar masses of HLIRG
neighbours with that of random galaxy neighbours and only find weak excess as
discussed in Section 3.2. Thus, we conclude that the lack of significant
excess of HLIRG neighbours compared to random galaxy neighbours is not due to
the method used to derive stellar masses.

##  Appendix B Spectroscopic follow-up observation of the HLIRG EN1_19770

Follow-up observation of EN1_19770 was carried out on June 20th of 2021 on the
Subaru 8m Telescope using the Faint Object Camera and Spectrograph (FOCAS)
instrument in clear sky conditions under long slit mode, the B300 grism and a
0.5” slit. Two 900 second integrations were taken along with associated
acquisition and calibration observations. We reduce FOCAS observations for
EN1_19770, a galaxy near the field of view of galaxy LEDA 2514599. We follow
the cookbook v1.0.3 of the instrument, where we perform the different steps
from the raw images using IRAF version 2.16 and the FOCASRED package. First,
we subtract the offset of the images (i.e. bias) from the read-out process of
the charge-coupled device (CCD). Second, we correct the intrinsic performance
of the CCDs by diving the images with the flat-field taken before the
observations. For this step, we use the flatnorm task to normalise the
spectrum using an specified region of the average detector image. Then, we
correct the distortion pattern on the FOCAS data using the distcalib task. We
calibrate the wavelength using the night sky lines for the object and the
emission light from a ThAr lamp for the standard star (HZ44). We perform the
sky-subtraction using the background task in IRAF by fitting the sky region
close to the two galaxies located in the same slit (EN1_19770 and LEDA
2514599). We calibrate the flux with HZ44 to them combine the different
exposures to remove the cosmic rays. Finally, we correct by interstellar
extinction and heliocentric velocity and extract the final 1D spectrum, as
shown in Figure 9. We derive a spectroscopic redshift of z = 3.74 based on the
detection of Si IV and Si II absorption lines as shown in the figure. This
source has a photometric redshift of 3.91 with a minimum value 3.82 in Duncan
et al. (2021). Thus, the photometric redshift estimate agrees well with the
derived spectroscopic redshift.

Figure 9: The optical spectrum for EN1_19770. We retrieve a spectroscopic
redshift of 3.74 based on the detection of Si IV and Si II absorption lines.

##  Appendix C Overdensity of the most promising protocluster candidates

In the left panel of each subfigure in Figure 10, we present the 250 $\mu$m
flux density distribution of neighbours within 100$\arcsec$ of the most
promising protocluster candidates. The median value and uncertainty are drawn
from 100 realizations. In the right panel of each subfigure, we plot the
overdensity parameter calculates as $\frac{N_{HLIRG}-N_{random}}{N_{random}}$
where $N_{HLIRG}$ is the number of neighbours aroud HLIRGs within difference
radius and $N_{random}$ is the mean number of neighbours around random
positions drawn from 100 random selections. Values above zero indicate an
overdensity.

Figure 10: Left: Cumulative 250 $\mu$m flux density distribution of HLIRG
neighbours within 100$\arcsec$ compared with random neighbours. Right:
overdensity parameters measured within different radius.

Figure: Continued

Figure: Continued