JCMT/SCUBA-2 uncovers an excess of 850 µ m counts on megaparsec scales around high-redshift quasars Characterization of the overdensities and their alignment with the quasars’ Ly α nebulae

We conducted a systematic survey of the environment of high-redshift quasars at submillimeter wavelengths to unveil and characterize the surrounding distribution of dusty submillimeter galaxies (SMGs). We took sensitive observations with the SCUBA-2 instrument on the James Clerk Maxwell Telescope for 3 enormous Lyman-alpha nebulae (ELANe) and 17 quasar fields in the redshift range 2 < z < 4 . 2 selected from recent Lyman alpha (Ly α ) surveys. These observations uncovered 523 and 101 sources at 850 µ m and 450 µ m, respectively, with signal-to-noise ratios (S / N) > 4 or detected in both bands at S / N > 3. We ran self-consistent Monte Carlo simulations to construct 850 µ m number counts and unveil an excess of sources in 75% of the targeted fields. Overall, regions around ELANe and quasars are overabundant with respect to blank fields by a factor of 3 . 4 ± 0 . 4 and 2 . 5 ± 0 . 2, respectively (weighted averages). Therefore, the excess of submillimeter sources is likely part of the megaparsec-scale environment around these systems. By combining all fields and repeating the count analysis in radial apertures, we find (at high significance, ≳ 5 σ ) a decrease in the overdensity factor from > 3 within ∼ 2 cMpc to ∼ 2 in the annulus at the edge of the surveyed field ( ∼ 10 cMpc), which suggests that the physical extent of the overdensities is larger than our maps. We computed preferred directions for the overdensities of SMGs from the positions of the sources and used them to orient and create stacked maps of source densities for the quasars’ environment. This stacking unveils an elongated structure reminiscent of a large-scale filament with a scale width of ≈ 3 cMpc. Finally, the directions of the overdensities are roughly aligned with the major axis of the Ly α nebulae, suggesting that the latter trace, on scales of hundreds of kiloparsecs, the central regions of the projected large-scale structure described by the SMGs on megaparsec scales. Confirming member associations of the SMGs is required to further characterize their spatial and kinematic distribution around ELANe and


Introduction
Quasars, accreting supermassive black holes at the center of galaxies, are one of the most luminous sources throughout cosmic time.Their extraordinary energy budget is expected to affect their host galaxies and the surrounding material out to intergalac-⋆ Tables F1 to F40 (the source catalogs) are only available in electronic form at the CDS via anonymous ftp to cdsarc.cds.unistra.fr(130.79.128.5) or via https://cdsarc.cds.unistra.fr/cgi-bin/qcat?J/A+A/ ⋆⋆ The reduced images are only available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr(130.79.128.5) or via http://cdsarc.ustrasbg.fr/viz-bin/qcat?J/A+A/ tic scales (e.g., Rees 1988;Silk & Rees 1998;Fabian 2012;King & Pounds 2015;Harrison 2017).Current galaxy formation models require this quasar feedback to explain several observables, especially at the massive end of the galaxy population and for intergalactic material in groups and clusters (e.g., Croton et al. 2006;Dubois et al. 2012;Crain et al. 2015;Weinberger et al. 2017;Costa et al. 2022).Bright quasars are relatively rare objects characterized by a strong evolution of their number density with redshift (e.g., Shen et al. 2020 and references therein).The substantial rise in the number of quasars from cosmic dawn to their peak at z ∼ 2 and the subsequent fall to z = 0 can be explained in the framework of the hierarchical growth of structures together with a decrease in fuel at late times (e.g., Efstathiou & Rees 1988;Kauffmann & Haehnelt 2000;Hopkins et al. 2007).Because of their paucity, bright quasars are usually assumed to pinpoint the location of the highest density peaks in the dark matter (DM) distribution.Indeed, clustering studies infer that high-redshift quasars should inhabit relatively massive DM halos in the currently probed redshift range, 0.5 < z ≲ 4 (M DM ∼ 10 12 − 10 13 M ⊙ ; e.g., White et al. 2012;Eftekharzadeh et al. 2015).
To confirm these expectations, several works have tried to directly constrain the density field around quasars by searching for associated unobscured companions (e.g., Lyman alpha emitters; e.g., Steidel et al. 2000) with different techniques and on different scales (from a few hundred kiloparsecs out to megaparsec scales), resulting in heterogeneous results (e.g., Willott et al. 2005;Kashikawa et al. 2007;Trainor & Steidel 2012;Kikuta et al. 2017;Mazzucchelli et al. 2017;García-Vergara et al. 2019).Possible causes of this lack of agreement could be related to (i) observational strategies or available instrumentation (i.e., small sample statistics for each study and/or shallow sensitivity) and/or (ii) a physical phenomenon: the galaxies in quasar environments are significantly dustier, resulting in a biased picture when looking only at unobscured tracers (e.g., García-Vergara et al. 2020).Recent efforts suggest that both effects could be at work and have started to confirm the expectations from clustering studies and, therefore, the association of quasars with rich environments.For example, Fossati et al. (2021) targeted 27 z ∼ 3 − 4.5 quasar fields with deep observations (∼ 4 hours) using the Multi Unit Spectroscopic Explorer (MUSE; Bacon et al. 2010) and found a significant overdensity of Lyman alpha emitters (LAEs) within 300 km s −1 of the quasars' systemic redshift and on scales of a few hundred kiloparsecs.From this overdensity and the kinematics of the LAEs, the authors inferred a DM halo mass in the range M DM ≈ 10 12 − 10 12.5 M ⊙ , which is in line with the aforementioned clustering measurements.Meanwhile, several studies have reported the presence of companions or overdense environments in submillimeter galaxies (SMGs; e.g., Smail et al. 1997) from tens of kiloparsecs up to a few hundred kiloparsecs around 2 < z ≲ 7 quasars, thanks especially to the Atacama Large Millimeter/submillimeter Array (ALMA; e.g., Decarli et al. 2017;Trakhtenbrot et al. 2017;Venemans et al. 2020;Bischetti et al. 2021;Chen et al. 2021;García-Vergara et al. 2022).Altogether, these works paint a coherent picture that stresses the importance of multiwavelength observations in determining the environment of quasars.
In this framework, submillimeter observations have targeted diverse active galactic nucleus (AGN) samples to constrain the richness of their environment on megaparsec scales.These works discovered an excess of submillimeter source counts around some high-redshift radio galaxies (e.g., Stevens et al. 2003;Dannerbauer et al. 2014), a few high-redshift quasars (Priddey et al. 2008), or IR-luminous AGN selected by the Widefield Infrared Survey Explorer (WISE) space telescope (Jones et al. 2017).For fields with large numbers of detections, it has also been assessed whether the SMGs trace large-scale structures with a preferential alignment with the radio galaxies' jets (Stevens et al. 2003).The latest work exploring this aspect targeted the fields of 17 high-redshift radio galaxies and found evidence for the alignment of SMGs with the perpendicular direction with respect to the radio jets (Zeballos et al. 2018).This misalignment has been interpreted as a signature that the SMGs trace a large-scale filament feeding the protocluster core.
Nonetheless, the overdensities are only found in a few fields out of larger parent samples, possibly suggesting that SMGs are not a universal tracer of overdensities around AGN or that they may preferentially coevolve around certain types of quasars (Rigby et al. 2014;Zeballos et al. 2018).To further investigate these hypotheses and extend the types of AGN environments studied, Arrigoni Battaia et al. (2018a) targeted one z ∼ 2 enormous Lyman alpha nebula (ELAN; e.g., Cai et al. 2017a).ELANe are rare systems characterized by an extreme Lyman alpha (Lyα) surface brightness (SB Lyα ∼ 10 −17 erg s −1 cm −2 arcsec −2 ) over > 100 physical kpc, with multiple embedded AGN and active galaxies (e.g., Hennawi et al. 2015;Arrigoni Battaia et al. 2018b;Chen et al. 2021;Arrigoni Battaia et al. 2022), making them good candidates for protoclusters (e.g., Hennawi et al. 2015).Arrigoni Battaia et al. (2018a) reported a factor of 4 overdensity for submillimeter sources on megaparsec scales, providing further evidence of the richness of ELAN environments.More recently, Nowotka et al. (2022) expanded that survey by targeting three additional ELAN fields and revealed a ubiquitous overabundance of SMGs.Assuming that all the submillimeter sources are associated with the targeted ELANe, the star formation rate density is found to be similar to that around other quasar samples or in protoclusters at similar redshifts (z ∼ 2−3; e.g., Clements et al. 2014;Dannerbauer et al. 2014;Kato et al. 2016).Therefore, ELANe seem to pinpoint regions where the coevolution of quasars and SMGs is favored.
Larger statistical surveys are needed to confirm these findings and assess whether ELANe inhabit very different environments than other bright quasars characterized by less extreme Lyα nebulae.In this paper, we report the results obtained by analyzing 20 fields centered on 17 quasars and 3 ELANe.The paper is structured as follows.Sections 2 and 3 describe the sample selection, the observations, their data reduction, and the source extraction, together with the criteria for making our catalogs.Section 4 reports on the calculation of the differential number counts.We highlight our results in Sect.5, including (i) an estimate of the overdensity factors for each field and on average; (ii) a study of the variation in the overdensity as a function of radial aperture; and (iii) a comparison of the overdensity factors with each system's properties (both quasar and Lyα nebula properties).Section 6 discusses (i) the implications of the similarity of overdensities found around all the studied systems; (ii) the presence of a characteristic scale for the overdensities determined by stacking maps of counts; and (iii) the alignment between overdensities and extended Lyα emission in these fields.Finally, we summarize our work in Sect.7. Throughout this paper, we adopt the cosmological parameters H 0 = 70 km s −1 Mpc −1 , Ω M = 0.3, and Ω Λ = 0.7.In this cosmology, 1 cMpc corresponds to about 0.5 arcmin at the mean redshift of our sample (z = 3.337).

Sample selection
We report on submillimeter data for 17 quasar fields in the redshift range z ∼ 3 − 4.2 and for 3 ELAN fields (Jackpot, MAMMOTH-I, Fabulous).The latter three objects have been already studied in Nowotka et al. (2022), but with shorter total exposure times than used here (∼3 versus ∼7 hours).All these systems have information on their Lyα emission on hundreds of kiloparsecs thanks to recent surveys with narrow-band filters or the MUSE instrument: Fluorescent Lyman-Alpha Survey of cosmic Hydrogen iLlumInated by hIGH-redshifT quasars (FLASH-LIGHT; Arrigoni Battaia et al. 2014), MApping the Most Mas- Notes. (a) ID from the QSO MUSEUM sample. (b) Systemic redshift of the quasar evaluated from the peak of the C IV line, correcting for the expected shift as estimated by Shen et al. (2016).The intrinsic uncertainty on this correction is ∼ 415 km s −1 and dominate the error badget (∆z ≈ 0.007). (c) Noise level at the center of each map. (d) Effective area assumed in this work and defined above the 3 times the central noise.The MAMMOTH survey identified overdensities at z = 2 − 3 by locating extremely rare, strong intergalactic H i Lyα absorption systems and absorption groups (Cai et al. 2017a).These fields were then observed with narrow-band filters targeting the Lyα line at those redshifts to characterize the population of LAEs.The results of this effort includes the discovery of the MAMMOTH-I ELAN (Cai et al. 2017b).The MAGG survey targeted 28 relatively bright (16.5 < r < 20.0) z = 3.2 − 4.5 quasars with the aim to study the halo gas of (i) z = 3 − 4.5 star-forming galaxies in the surroundings of optically thick H i absorbers (Lofthouse et al. 2020), (ii) z = 3 − 4.5 quasars (Fossati et al. 2021), and (iii) z ∼ 1 galaxies (Dutta et al. 2020).
The targets of our work were selected to cover the full range of Lyα nebulae (e.g., area) and quasars (e.g., luminosity) properties of the parent sample, after taking into account source visibility.
Five systems have a radio-loud central or companion quasar with unresolved morphology at the resolution of the Very Large Array (VLA) Faint Images of the Radio Sky at Twenty-Centimeters (FIRST) survey (∼ 5 arcsec; Becker et al. 1994) or of the NRAO VLA Sky Survey (NVSS; ∼ 45 arcsec; Condon et al. 1998).Figure 1 shows the absolute magnitude at 1450 Å M 1450 of the central AGN and the areas of the Lyα nebulae around the targeted systems in comparison to the aforementioned parent sample.In Table 1 we list the selected sources with their coordinates and redshifts.

SCUBA-2/JCMT observations and data reduction
The observations for the twenty fields in this study ( For the data reduction, we closely followed the method used in Chen et al. (2013a), Arrigoni Battaia et al. (2018a), andNowotka et al. (2022).Briefly, the data were reduced using the Dynamic Iterative Map Maker (DIMM) included in the SMURF package from the STARLINK software (Jenness et al. 2011;Chapin et al. 2013).The standard configuration file dimmcon-fig_blank_field.lis was adopted for our science purposes.Data were reduced for each scan and the MOSAIC_JCMT_IMAGES recipe in PICARD, the Pipeline for Combining and Analyzing Reduced Data (Jenness et al. 2008), was used to co-add the reduced scans into the final maps.
Point source detectability was increased by applying a standard matched filter to the final maps, using the PICARD recipe SCUBA2_MATCHED_FILTER.For flux calibration, we adopted the recommended flux conversion factors (FCFs) with 10% upward corrections obtained from the analysis of archival data spanning ten years1 .These FCFs are date dependent, and we list the relevant values for our work in Appendix A. We also generate calibrated maps using the standard FCFs (491 Jy pW −1 for 450 µm and 537 Jy pW −1 for 850 µm) with 10% upward corrections (Dempsey et al. 2013).These values have a difference of only 5 − 10% with respect to the new FCFs.We tested our results against such a change in FCFs and our findings are unchanged.
The final central noise levels for our data are in the range 0.53-1.07mJy/beam and 3.82-20.50mJy/beam at 850 µm and 450 µm, respectively.In the reminder of this work we focused on the regions of the data characterized by a noise level less than three times the central noise.We refer to this area as the effective area.Table 1 summarizes the observations log, including the map center, the effective area, the central noise, the total exposure time, the opacity, and the program ID for each field and band.

Source extraction and catalogs
We extracted the detections from all maps following Arrigoni Battaia et al. (2018a) and Nowotka et al. (2022).In brief, we focused on the effective area of each map (see Table 1) and found all sources with peak S/N > 3.This is done by looking recursively for the maximum pixel within the selected region, extracting the position and the information of the peak, and subtracting a scaled point spread function (PSF) centered at that position.The process was stopped when the peak S/N went below 3. We used these list of sources as preliminary catalogs, and crosschecked the catalogs at 850 and 450 µm for each field to select counterparts in the other band.A counterpart is defined as a detected source at 450 µm laying within the 850 µm beam.
In the final catalogs we list every > 4σ source, but also every > 3σ source characterized by a > 3σ counterpart in the other band.Considering all fields, we discovered 523 sources at 850 µm and 101 sources at 450 µm.In Tables F1 to F40 (available only in electronic form), we list the positions, positional error, S/N, fluxes, deboosted fluxes and counterparts for all these sources.In the same tables we also list the distances between detections at 850 µm and their counterparts at 450 µm.Figure 2 and  Figures C.1 to C.6 show the final S/N maps at 850 and 450 µm for each targeted field with the discovered sources overplotted.

Reliability of source extraction
The aforementioned catalogs could be affected by spurious sources.To quantify them, we use several methods as in Arrigoni Battaia et al. (2018a).First, we inverted the final maps and applied the source extraction algorithm.The median number of sources per field at > 4σ at 850 and 450 µm are found to be six and one, respectively.We stress that the matched filter technique introduce strong negative features in the vicinity of the brighter sources in each field.This results in a larger number of detections in the inverted maps for fields with very bright sources and therefore making this first method not reliable in estimating false detections.Second, to get a better estimate, we applied the source extraction algorithm to true noise maps, and checked the number of detections with > 4σ.The true noise maps were constructed by using the jackknife resampling technique.In brief, for each field and band, we obtained two maps by co-adding roughly half of the data.The true noise maps were then the result of the subtraction of these two maps.Indeed, any real source in the maps should be subtracted irrespective of its significance, resulting in residual maps that are source-free noise maps.Importantly, to account for the slight difference in exposure time, the true noise maps were scaled by a factor of √ t 1 × t 2 /(t 1 + t 2 ), where t 1 and t 2 indicate the exposure time of each pixel from the two maps.These jackknife maps have a central noise in agreement with the science data to < 4% and < 11% at 850 µm and 450 µm, respectively.In these maps, the median number of sources at > 4σ is one and two at 850 and 450 µm, respectively.We thus expect a similar number of spurious sources in our > 4σ source catalogs, which corresponds to a median false detection rate of ∼ 5% and ∼ 45% at 850 and 450 µm, respectively.Table D.1 lists the numbers for each field and band for both methods.
In addition, as done in Arrigoni Battaia et al. (2018a), we estimated the number of spurious detections for the sources with 3 <S/N< 4 identified as having a counterpart in the other bandpass (lower portion of catalog tables in Appendix F).We proceed as follows.First, sources between 3 and 4σ were extracted from the jackknife maps and cross-correlated with the detections in the real data at the other wavelength.An average of 0.2 and 1.0 spurious sources matched a detection in the real maps.This test was then repeated 1000 times using random positions for the spurious sources within the effective area of our observations.The average of spurious sources at 850 and 450 µm was found to be 0.03 and 0.25, stressing that > 3σ sources selected in both bandpasses are even more reliable than > 4σ sources selected in only one bandpass.Table D.1 lists the numbers for each field and method.
Altogether, the sources at 450 µm without a detection at 850 µm are most likely spurious for the shallower observations of the 17 quasar fields.For completeness, we report all the sources in our catalogs.Our analysis focuses on the 850 µm data and, therefore, our conclusions are not affected by the spurious 450 µm sources.1).The 850 and 450 µm maps for all other fields are shown in Appendix C.

Completeness tests
For each field and band, we tested at which flux the catalogs can be considered complete by populating the true noise maps discussed in Sect.3.1 with randomly placed mock sources of a given flux.We then use the extraction algorithm and considered as recovered sources those detected above 4σ and within the beam area.In more detail, sources were injected within the effective area with fluxes from 0.1 to 30.1 mJy (0.1 to 80.1 mJy) with a step of 0.5 mJy (1.0 mJy) for 850 (450) µm.For each step in flux, we iterated the extraction by introducing 1000 sources.The results of these tests are shown in Fig. E.1.As expected, the three ELAN fields (MAMMOTH-1, Jackpot, Fabulous) have deeper data with the 50% completeness being on average around 3.2 and 24.0 mJy at 850 and 450 µm, respectively, and the 80% being on average around 4.0 and 30.4 mJy, respectively.The other fields have the average 50% and 80% completeness at 850 µm around 5.4 and 6.9 mJy, respectively, while the completeness at 450 µm is much worse (≲ 50% at 55 mJy).SDSSJ1209+1138 is the deepest of the quasar fields with the 50% and 80% completeness at 850 µm around 3.9 and 4.9 mJy, respectively.

Differential number counts
In this section we outline the calculation of the differential number counts and the underlying counts model for each field at 850 µm.Our aim is to ascertain the presence of an excess of sources in these fields with respect to blank reference fields (Sect.5.1).We proceeded following the method of Chen et al. (2013a), already used in the case of ELANe by Arrigoni Battaia et al. (2018a) and Nowotka et al. (2022).First, we determine down to which S/N there is a significant excess signal with respect to a pure noise distribution for each field, by inspecting the S/N histograms of the signal maps and the jackknife maps (see Figs. B.1,B.2,and B.3).We consider as significant excess a factor of ≥ 3× more pixels at a given S/N.This is the case down to a median S/N = 2.5 (dashed vertical line in those figures; see Table 2 for the specific values).This positive excess signal is due to real astronomical sources, while the negative excesses are due to the negative troughs of the matched-filter PSF (Chen et al. 2013a).This effect can be properly accounted for during sources extraction with a correct PSF.
As a next step, we extracted sources lowering the detection threshold to the aforementioned smaller value for each field.This can be done as positional information is not relevant for number counts analyses (e.g., Chen et al. 2013b).In addition, we extracted sources with the same thresholds from the respective jackknife maps.We then computed the raw number counts in equally spaced logarithmic bins between the minimum and maximum flux densities of the detected sources.To do so, we summed the number of sources in each flux bin and divided by the width and detectable area of the bin.The detectable area defines where a source with a certain flux density can be detected above the S/N threshold.The raw number counts were then obtained by subtracting off the spurious counts obtained from the jackknife maps.The average fraction of spurious sources down to the thresholds considered here is 40%.
We have then corrected for the known observational biases that affect the raw counts (e.g., flux boosting, incompleteness, source blending) with the use of Monte Carlo simulations following the exact procedure by Nowotka et al. (2022).These simulations determine the underlying count model able to correctly match the raw number counts in each field, and have been shown to be consistent to previously published results for both a known Notes. (a) ID in QSO MUSEUM. (b) S/N threshold used for source extraction for computing the number counts.
blank field, the Chandra Deep Field-North (CDF-N; Chen et al. 2013a), and an overdense field, MAMMOTH-I (Arrigoni Battaia et al. 2018a).We briefly summarize how the simulations proceed.More details can be found in Nowotka et al. (2022).First, to ease comparison with previous works, and in particular with the chosen fiducial blank field count model (Geach et al. 2017), we parametrized the underlying count model as a Schechter function of the form where N 0 (units of deg −2 ) is the normalization and S 0 (units of mJy) is the characteristic flux.To avoid degeneracy between N 0 and S 0 , we keep the value of S 0 fixed to 2.5, as done in Geach et al. (2017) and Nowotka et al. (2022).We fit this function to the raw number counts to obtain the parameters of an initial model used to inject mock sources in the jackknife map at random positions.We injected sources down to 1 mJy, the same lowest flux adopted by the model in Geach et al. (2017), to which we compare our results to estimate overdensities.We constructed in this way 100 simulated maps for each set of model parameters, and obtain the mean recovered counts.A new set of counts was then created by multiplying the raw counts by the ratio between raw counts and mean recovered counts.At this point, the simulation starts a new iteration from the Schechter function fitting.The process was repeated until the recovered counts fall within the 1σ error of the raw counts.Table 2 lists the parameters of the underlying models found using these simulations.We confirmed that our Monte Carlo simulations properly account for observational biases by creating 500 simulated maps for each field with mock sources injected into the respective jackknife map following the estimated underlying models.the quasars or sources suspected to be associated with these systems given the low probability of a chance detection (Arrigoni Battaia et al. 2018a;Nowotka et al. 2022).Similarly, for the Fabulous ELAN and SDSSJ 1342+1702, we applied the convergence criteria but neglecting the last bin, which clearly deviates from a Schechter due to the presence of a few bright sources likely associated with the systems (see, e.g., Arrigoni Battaia et al. 2022 for the Fabulous ELAN).
Finally, to test whether this procedure is able to assess the presence of an overdensity, we create 500 additional simulated maps, but this time injecting sources according to the fiducial blank field model from Geach et al. (2017) with the parameters N 0 = 7180 ± 1220 deg −2 ,S 0 = 2.5 mJy, and γ = 1.5 ± 0.4.Figure 3,G.1,G.2,and G.3 show that the realizations of the blank field model (gray shaded region) produce recovered counts that are lower than the raw counts in most fields, highlighting the presence of overdensities.We compute overdensity factors in Sect.5.1.

Overdensity factors
Our analysis of the differential number counts and the identification of the underlying models point out that most of the targeted fields are overdense in comparison to the reference blank field model by Geach et al. (2017).In this section we quantify the overdensity factor for each field following the two approaches presented in Nowotka et al. (2022).
First, the raw counts in each flux bin are corrected for observational biases by dividing them by the ratios of the mean recovered counts and the obtained underlying model.The blank field model is then fit to this corrected counts with only the normalization, N 0 , as a free parameter.The overdensity factors are then estimated as the ratio between the obtained normalization and the blank field model normalization.In Table 3 we list the obtained overdensity factors together with their errors taking into account both Poisson noise and uncertainties from the simulations.We find a weighted average overdensity factor for ELAN fields of 1.9 ± 0.2, which agrees within 1σ of the value found in Nowotka et al. (2022) using the same method.For the overall quasar fields we instead obtain a weighted average overdensity Notes. (a) ID in QSO MUSEUM. (b) Overdensity factor estimated as ratio of model normalizations with respect to the reference blank fields. (c) Overdensity factor estimated as ratio of the cumulative number of sources with respect to the reference blank fields.
factor of 1.8 ± 0.1, with the QSO MUSEUM sample showing higher values with respect to the MAGG sample, 2.3 ± 0.2 and 1.6 ± 0.1, respectively.This first method is appropriate if the shape of the underlying model is consistent with that of the blank field model.However, for many of the fields there is evidence that the models are steeper, meaning that the overdensity could be flux dependent (Figs. 3,G.2,and G.3), and higher than estimated with the first approach (Nowotka et al. 2022).To check this, we computed the cumulative number of sources for each field by integrating the models in the flux range defined by the minimum and maximum source fluxes, and compare it with the blank sky value within the same range.As expected, we find on average larger overdensity factors of 3.4 ± 0.4, 2.5 ± 0.2, 3.1 ± 0.4, and 2.3 ± 0.3, for the ELAN, all 17 quasars, and the QSO MUSEUM and MAGG fields, respectively.The overdensity factors obtained with this second approach for all the fields are also listed in Table 3.Since this second approach does not suffer from any bias in the selec-tion of a fiducial model shape, we base our discussion on these values (see Sect. 6.1).
We stress that the data at 450 µm do not provide enough statistics for a robust analysis of the counts as done for the 850 µm counts.However, for the ELAN fields in which we have better sensitivity, we checked the number of detected sources for fluxes above the 80% and 50% completeness level and compare it to the available blank field values at 450 µm (e.g., Chen et al. 2013b;Geach et al. 2013;Zavala et al. 2017), taking into account the completeness and the deboosted flux for our sources.We found that the ELAN fields have on average 1 +2.3  −1.0 (6 +5.8 −3.3 ) sources for fluxes > 30.4 mJy (> 24.0 mJy), which are consistent with the average blank field values from the aforementioned works, 0.6 +1.3 −0.6 (1.5 +2.8 −1.4 ).The absence of a 450 µm excess provides further evidence that the 850 µm overdensities are likely at z > 2. Further, we note that most of the sources in our catalogs have a flux ratio f 850 / f 450 consistent with expectations from the SMG template of the ALESS survey (Swinbank et al. 2014;da Cunha et al. 2015) redshifted to our targeted fields.In particular, we found that only 1.5% of all the sources in our catalogs have the flux ratio lower than the ALESS expectations.It is therefore likely that most of the detected sources are at a redshift similar to the targeted fields.

Source counts in radial apertures
The number of fields and the large number of sources detected at high significance (S/N> 4) at 850 µm allow us to study the source counts in radial apertures around the targeted systems.A similar analysis has been performed by Zeballos et al. (2018), who targeted a sample of 17 high-redshift radio galaxies (HzRGs) at 0.5 < z < 6.3.They used radial apertures defined by radii every 1.5 arcmin, out to 4.5 arcmin, finding marginal evidence of an overdensity (∼ 3× the blank fields) for the central aperture (0 < r < 1.5 arcmin).Given the large spread in redshift of their sources, this result could be biased as overdensities span different physical scales at different epochs.
For this analysis we decided to follow two approaches: following the analysis in Zeballos et al. (2018) to allow a oneto-one comparison, but also repeating the analysis using apertures in comoving units.First, we compute the corrected source counts in apertures similar to Zeballos et al. (2018): a circle of radius r = 1.5 arcmin, and then three additional annuli with 1.5 ′ < r < 3 ′ , 3 ′ < r < 4.5 ′ , and 4.5 ′ < r < 6 ′ , basically covering most of the effective area of each map.We start from our catalogs and correct them for completeness and flux boosting using curves obtained as described in Sect.3.2.However, since the rms in our maps degrades toward their outskirts, we compute completeness and flux boosting for the different apertures in each map.Each source flux is then corrected by dividing by the obtained median boosting factor at its observed flux.The contribution of each source to the counts is computed by dividing its single occurrence by the completeness at its observed flux.We then computed the number counts in logarithmically spaced flux bins, covering the minimum and maximum flux in each aperture.The top panels in Fig. 4 show the differential and cumulative number counts obtained with this method.For the inner aperture, we compute the source counts with two approaches: including (red circles) or excluding (purple pentagons) the possible counterparts to the targeted systems (i.e., sources at r < 10 arcsec; see the tables in Appendix F).The second approach is binned differently as the brightest sources are neglected.For both approaches, it can be seen that the inner aperture is consistently more overdense than the others, having counts ∼ 3× higher than blank fields at S 850µm > 4 mJy.Interestingly, as can be better seen below the top-right panel of Fig. 4, there is a promising trend indicating a gradual decrease in the overdensity for S 850µ > 4 mJy from the inner to the outer aperture.The outer aperture seems still slightly overdense (∼ 2×) with respect to the field, probably meaning that the physical extent of the overdensities around these systems is larger than our maps.This finding would be in agreement with theoretical expectations by Chiang et al. (2013) that place the effective radius of protoclusters associated with halo masses M DM > 10 12 M ⊙ at physical scales similar to the size of our maps (∼ 10 cMpc).Similar predictions have been provided by additional simulations (e.g., Muldrew et al. 2015), which qualitatively agree with the large extents for protoclusters reported in observations (e.g., Casey 2016).Table 4 lists the differential and cumulative number counts for all apertures.For the inner aperture we only list the values including all sources, as the results of the two aforementioned approaches are consistent within the uncertainties.
We quantified the significance of the obtained overdensities in different apertures by using the 10000 (500 per system) synthetic maps for blank fields obtained in Sect. 4 using the jackknife maps.Indeed, the comparison between real and mock maps allows us to estimate the significance of the radial trend because in the mock maps, sources are injected at random positions.We focused on comparing the total number of sources detected in all 20 fields for S/N> 4, which is the threshold for our catalogs.Specifically, we first computed the total number of sources per aperture for our catalogs, and then repeated the same experiment 500 times for samples of 20 mock blank fields, one for each targeted system (i.e., with the noise properties of that field).Comparing the number of sources obtained in the observed maps and in the mock maps allows us to obtain a significance without the  2) around the targeted systems with respect to an estimate for blank fields (dashed gray; Geach et al. 2017).The dot-dashed and dotted lines show models with the same γ and S 0 , but with 2× and 4× higher normalization (N 0 ), respectively.The panel in the bottom-right corner shows the ratio between the number of sources for each aperture and the reference blank field model.For a better visualization of the data points at low fluxes, this plot is cut at a ratio of 50, and therefore it does not show the uncertain data points at large fluxes, which by far exceed the blank field model.A trend of decreasing overdensity from the inner to the outer aperture is evident.Bottom: Similar to the top panel, but for the differential number counts obtained in apertures defined in comoving megaparsecs (see the legend and Sect.5.2) around the targeted systems with respect to an estimate for blank fields (dashed gray; Geach et al. 2017), but translated to comoving megaparsecs for the median redshift of our sample.The inset plot compares the number of detected sources with S/N> 4 for the inner bin with the density distribution function obtained from 500 realizations of 20 mock blank fields, similar to what is done in Fig. 5 with apertures in arcminutes.need of deboosting or completeness corrections.Figure 5 shows the density distribution functions obtained in each aperture for the 500 realizations of blank fields (histograms), which can be well fit by a Gaussian, divided by their mean value.We then computed the overdensity in each aperture as the ratio between the total number of observed counts divided by the mean of the Gaussian distribution for the respective blank field realizations.The total number of observed counts (vertical line with 1σ Poisson uncertainty) in each aperture exceeds the blank field by a factor of at least ∼ 2 in the outer annuli and increases to ∼ 3 in the central aperture, in agreement with the previous analysis of number counts (top panels of Fig. 4).We estimated the significance of the observed overdensities in each aperture by as-suming that the uncertainty on our measurement is described by Poisson statistics in the presence of a Gaussian background (the distributions for the blank-field).The significance can be then obtained following Vianello (2018) (their Eq. 15).We find that the overdensities in each aperture are detected at ≥ 4.7σ.Table 5 lists the overdensities in each aperture, together with their significance and the parameters of the Gaussians.
Further, we tested whether the use of apertures defined using comoving megaparsecs rather than arcminutes would result in an even stronger detection of overdensities as a function of aperture.We redid the previous analysis using a circle of radius r = 2 cMpc and then three additional annuli with 2 < r < 5, 5 < r < 8, and 8 < r < 10 cMpc, covering the portion of effective areas in Fig. 5: Overdensities of detected sources with S/N> 4 in different apertures in arcminutes (colored vertical lines) compared to the density distribution functions of overdensities from 500 realizations of 20 mock blank fields with similar noise properties as the 20 science maps, each divided by their mean value (gray distributions).The overdensity in each aperture is detected at a high significance, with the strength of the overdensity decreasing from the central to the outer aperture (see Table 5).The errors on overdensities take into account the Poisson uncertainties on the counts.The significance of the overdensities also taking the Gaussian "background" into account is given in Table 5.
Table 5: Data for the significance of the overdensities in different apertures around the targeted systems.

Aperture
Overdensity Significance µ BF Gauss σ BF Gauss 0 < r < 1.5 arcmin 3.1 6.9 32.6 5.4 10 arcsec < r < 1.5 arcmin 2.7 5.7 31.2 5.4 1.5 < r < 3 arcmin 1.9 4.7 61.6 7.5 3 < r < 4.5 arcmin 1.9 5.2 72.2 7.6 4.5 < r < 6 arcmin 1.9 5.4 79.6 7.1 common for different redshifts.The differential and cumulative number counts within these apertures are shown in the bottom panels of Fig. 4. As expected, given the relatively small redshift range, they confirm the results of the previous analysis.In particular, we find that the use of comoving megaparsecs slightly increase the overdensity value.For example, for the inner region we find 3.5× higher counts at S 850µm > 4 mJy.Once again, we quantify the significance of the overdensity, and find that the overdensity in the inner region is detected with a significance of 6.1 (see the inset in the bottom-right panel of Fig. 4).We further characterize the detected overdensities in Sect.6.2.

Overdensity factors versus system properties
In this section we investigate the presence of any trend between the obtained overdensity factors and the properties of the quasars and their associated extended Lyα emission.First we compare with the quasar properties.Figure 6 shows the overdensity factors plotted against redshift (top panel), quasar magnitude at restframe 1450 Å, M 1450 (middle panel), and the deboosted flux at 850 µm, f Deboosted 850 for detections within 10 arcsec of the optical quasar position (bottom panel).Different colors and symbols denote different samples: ELANe targeted in this work (blue stars), Slug ELAN (green star; from Nowotka et al. 2022), QSO MUSEUM fields (orange circles), and MAGG fields (black diamonds).We note some interesting features in the data that need to be verified with larger samples and with spectroscopically confirmed overdensities.
In the top panel of Fig. 6, it is remarkable that the three fields at z ≥ 4 show the smallest source excess in the sample.This fact could be due to a weaker coevolution between SMGs and quasars at high redshifts, or to the known redshift distribution of SMGs in blank fields, which peaks at z ∼ 2 − 3 (e.g., Dudzevičiūtė et al. 2020).In other words, estimating overdensities with respect to blank fields at these high redshifts could result in a stronger dilution of the signal than at z ∼ 2 − 3. The middle panel is confirming that the luminosity of the central AGN is not linked to its large-scale environment.Quasars with similar absolute magnitudes can have a different excess of submillimeter sources.This is in agreement with clustering studies for quasars that show no dependence on the quasar luminosities in the range (−28.ELANe (Hennawi et al. 2015;Cai et al. 2017;Arrigoni Battaia et al. 2018) Slug ELAN (Cantalupo et al. 2014) MAGG (Fossati et al. 2021) QSO MUSEUM (FAB19) Fig. 7: Cumulative overdensity factor δ cumul obtained in Sect.5.1 versus distance from the luminosity-comoving area relation for quasars' Lyα nebulae, in units of its rms.For the luminosity-area relation we use the reference inliers fit in Arrigoni Battaia et al. (2023).Each data point has a circle that is logarithmically scaled following the comoving area of the Lyα nebulae.Magenta dots indicate systems with radio-loud sources.The dotted line with shaded area indicates the relation and one rms = 0.17 dex. the redshift range 2.2 ≤ z ≤ 3.4 ; Eftekharzadeh et al. 2015).The bottom panel shows that central radio loud quasars have counterparts at 850 µm with f Deboosted 850 ≳ 7 mJy, and all have 850 µm source excesses ≳ 2.8× the blank fields.Similarly high values are found also for the two ELANe with a radio-loud companion.Six out of 14 (or 43 +26 −17 %) radio-quiet fields show lower cumulative overdensities than the radio-loud objects.Therefore, being a radio-loud quasar seems to guarantee a relatively high 850 µm source excess.This ubiquitous excess around radio-loud quasars is in contrast with the more heterogeneous source counts found around HzRGs (e.g., Dannerbauer et al. 2014;Zeballos et al. 2018).This fact could be related to the different inclination of the AGN jets with respect to the line-of-sight, and how they are expected to align with the large-scale structure.Indeed, jets should be aligned perpendicular to the large-scale filaments traced by the SMGs' positions (e.g., Zeballos et al. 2018 and references therein).If this is indeed the case, for some HzRGs (as they have the jet on the plane of the sky), the filaments could be aligned along the line-of-sight, making the overdensity of galaxies more compact and difficult to detect, on average.Instead, radio-loud quasars with unresolved radio emission (as in our sample) are expected to have their jet (if any) almost aligned with our lineof-sight (e.g., Sbarrato et al. 2022), implying a filament direction roughly on the plane of the sky.We find a first piece of evidence of this effect when looking at the different results for the overdensities obtained in apertures between our study and Zeballos et al. (2018).On average, our fields are overdense out to larger distances than those around HzRGs, which only show a clear excess for < 1.5 arcmin.This hypothesis could be tested by spectroscopically confirming the association of galaxies and reconstructing in 3D the large-scale structure around these systems.
Second, we looked for relations between the overdensity factors and the properties of the Lyα nebulae.Specifically, as we are comparing extended structures at different redshifts, we computed luminosities and comoving areas for the Lyα nebulae at a fixed redshift-corrected common threshold in surface brightness (SB) using the maps from the works by Cantalupo et al. (2014), Hennawi et al. (2015), Cai et al. (2017b), Arrigoni Battaia et al. (2019), andFossati et al. (2021).As done in Arrigoni Battaia et al. (2023), we chose the SB threshold 1.23 × 10 −17 erg s −1 cm −2 arcsec −2 at z = 2.0412 (Jackpot), and obtained the threshold at different redshifts by considering the cosmological SB dimming (1 + z) 4 .At each redshift these values are well above the 2σ threshold usually used to select extended Lyα emission.
We then compare the obtained values against the reference inliers fit of the luminosity-area relation in Arrigoni Battaia et al. (2023), and look for any trend with respect to the cumulative overdensity factors.Figure 7 shows the results of this exercise.Specifically, the plot shows the overdensity as a function of the distance from the luminosity-area relation (without taking the absolute value).Also, the size of the data points indicates the comoving area of the nebulae.Importantly, the MAGG quasar's nebulae lie on the luminosity-area relation discovered in that work, which uses only 2 < z < 3.5 quasars (dotted line with rms region).Therefore, this fact confirms that the relation is in place for higher redshifts (up to z = 4.2).Regarding the overdensity factors, given their uncertainties, there is no clear trend with respect to the distance from the relation.
Nowotka et al. ( 2022) obtained 850 µm number counts for three ELAN fields, finding on average overdensities of 3.6 ± 0.6.As we show in Sect.5.1, the overdensities around our sample are similar to those around ELANe.We further emphasize this by obtaining average differential number counts for the ELAN and the other quasars, from QSO MUSEUM and the MAGG sample.
Figure 8 shows the comparison between these average curves and the blank field counts.All three curves show similar overdensities.Given the rarity of the ELANe, this fact could mean that they represent only a short transitory phenomenon in a massive system evolution and/or a peculiar illumination or activity.
To verify this picture and the ubiquity of overdensities around the current sample and in general around high-z quasars, we need to spectroscopically confirm the association of the detected sources.Indeed, literature studies start to provide evidence for CO-emitting companions and/or overdensities of CO emitters associated with high-z quasars, but on smaller scales than here probed (e.g., Decarli et al. 2017;Chen et al. 2021;Li et al. 2021a;García-Vergara et al. 2022)

A characteristic spatial scale and orientation of the overdensities
Section 5.2 shows that the overdensities around the targeted systems are more significant within a radius of 2 cMpc, and then smoothly drops for larger circular apertures, without reaching the blank field counts for the explored area (out to 10 cMpc).To go beyond the simple circular apertures, we used two methods to define preferred directions for each of the 20 fields, assuming that the detected overdensities are associated with the central sources.In the first method, we sampled the directions on the plane of the sky that pass through the main quasar position in steps of 0.5 o (angle α measured from east to north), and computed the rms distance of the sources to each direction 2 .Then, we took the direction α rms as the preferred orientation, which minimizes the rms.Straightforwardly, if there exists a well defined angle for which the rms distance is very small, it means that the sources are distributed along a projected large-scale structure or filament defined by that angle.To constrain well α rms , we draw 10000 random subsamples (90%) of sources for each field (bootstrap with replacement), and for each one of them find the angle minimizing the rms distance.The distributions of these angles for three fields (one ELAN, the most overdense quasar field, and the least overdense quasar field) are given by the blue histograms in the left panels of Fig. 9 (top, central and bottom).The corresponding distribution for the remaining 17 fields are shown in Figs.I.2, I.3, I.4, and I.5.As the histograms can be quite noisy, we smoothed them using the Gaussian kernel density estimate to find the value of α rms where the distribution peaks.We associated errors with this peak value by computing the 34 percentiles left and right of peak.
The second method involves computing the 2D inertia tensor from the position of the sources and taking as preferred direction α inertia the eigenvector corresponding to the smallest eigenvalue (Zeballos et al. 2018).As for the rms method, when computing the inertia tensor, we weight the sources by their completeness.
2 Each source in the field is weighted by its completeness.
We draw 10000 random subsamples (90%) of sources for each field to construct the distributions of α inertia , from which we find the values corresponding to the peak, and the 34 percentiles left and right of peak.The distributions of α inertia are given by the orange histograms in the left panels of Figs. 9, I.2, I.3, I.4, and I.5.The angles α inertia are also measured north of east.
At this point, we can rotate each field such that its preferred direction, given either by α rms or by α inertia , becomes the new x-axis, and stack them together.The maps of Fig. 10 show the results of this exercise, where all the 20 fields have been used to construct the maps on the left panels, while the central and the right panels show the stacked maps of mean number of sources (corrected for completeness) per cMpc 2 for the ELAN and quasar fields, respectively.For the maps on the top row the preferred direction is given by α rms , while for the ones on the bottom by α inertia .The black contours in each map enclose increasing fractions of sources from 0.1 for the innermost contour to 1.0 for the outermost one in steps of 0.1.From these stacks of surface number density maps, we can see that the distribution of sources is more asymmetric (elongated along the x-direction) for the ELAN fields than for the quasar ones when looking at the central part of the maps.Also, it is clear that for all three cases, contours become increasingly more circular when moving outward from the center.
We quantify the asymmetry in the distribution of sources within these contours in the right most panels of Fig. 10.In practice, we use all the sources within each main contour (we do not consider the small secondary ones) to compute the Stokes parameters Q and U, again weighting each source by its completeness w i 3 .From the Stokes parameters we can estimate the asymmetry as . If all the sources are collapsed along a single direction, the asymmetry reaches its maximum value of 1, while for randomly distributed sources the asymmetry would be 0. The colored points in the profiles of Fig. 10 correspond to the stacks obtained after rotating each field by its peak α angle (α rms and α inertia for the top and bottom panel, respectively).The errors on the points have been obtained as the standard deviation of the asymmetry at each level over 1000 stacks, where each stack has been constructed using rotation angles drawn (independent for each field) from the angle distributions around the peaks.The horizontal dashed black line gives the prediction for the average stack of all 20 fields when the positions of the sources are randomized (distributed uniformly in the field-of-view).
Irrespective of the method used to define a preferential angle (rms distance and inertia tensor for the top and bottom panels, respectively), there is a marked difference between the asymmetry profile of ELAN and quasar fields.The ELAN fields show a broken profile, with large asymmetries (from 0.5 to 0.8) for the region enclosing ∼50% of the sources (within ∼6 cMpc of the central object), and lower asymmetries (∼0.3) for regions enclosing more than half of the sources.In contrast, quasar fields have a slowly decreasing profile, from asymmetries of ∼0.5-0.6 in the inner regions to ∼0.3 for the full field-of-view.Given the low number statistics, it is important to judge these results in comparison with the asymmetry level corresponding to positions of the sources randomly distributed within the field-of-view (horizontal dashed black line in both right most panels of Fig. 10).For all three combinations of stacks (only ELAN, only quasars, or all fields), the profiles of asymmetry are clearly above the level of 0.06±0.03predicted for randomly distributed sources.An ad- In each plot, the direction of the overdensity is color-coded in blue or orange depending on the method used, rms or inertia tensor, respectively (see Sect. 6.2).We show these plots for one ELAN (Fabulous) and for the two quasar fields with the largest (SDSSJ0819+0823) and smallest (J020944+051713) overdensity factor in our sample.In the central panels, the source sizes and colors indicate their completeness, which is used as a weight in the calculation of the direction of the overdensity.The Lyα SB maps are shown with a logarithmic color scheme stretching from 10 −19 to 5 × 10 −17 erg s cm −2 arcsec −2 , and with contours at S/N= 2, 4, 10, 20, and 50.We stress that the Lyα maps are on much smaller scales (see their scale bars) than the 850 µm detections in the central panels.
ditional interesting result is that the stack of only ELAN fields reveals clearly the presence of companion overdensities to the overdensity centered on the main quasar position, at (0,0) in the maps.The presence of these likely merging large-scale structures explains the broken asymmetry profile in ELAN fields.Indeed, the companion structure circularizes the overdensity on smaller scales in comparison to quasar fields.In Appendix H we show that (i) the secondary peak in the ELAN stacked map is due to chance superposition of overdensities in the individual fields, and (ii) the asymmetry profile is robust against 180 degree individual rotations of each of the three ELAN fields.
Further, we quantify a scale width for the elongated structure evident in the stacked maps and reminiscent of a large-scale filament.We focus on the stacked map obtained by using all fields for all the fields (left), only the ELANe (center), and only the quasar fields (right).For the maps on the top row the preferred direction is given by α rms and the ones on the bottom by α inertia .The black contours in each map enclose increasing fractions of sources from 0.1 for the innermost contour to 1.0 for the outermost one in steps of 0.1.Right: Asymmetry profiles of the stacked source distribution maps.The x-axis gives the radius of the circularized area of the contours, √ area/π, in comoving megaparsecs.From left to right, the points correspond to the asymmetry within the contours enclosing increasing fractions of sources, from 0.1 to 1.0 in steps of 0.1.For the profiles in the top panel, the stacking has been done after rotating each field by its α rms , while for the ones in the bottom panel the fields have been rotated by α inertia .The horizontal dashed black line together with its gray band (in both panels) marks the level of asymmetry of 0.06±0.03corresponding to stacks of all 20 fields where the positions of the sources have been randomized.rotated according to their α rms , and compute the distribution of distances (|y|) from the x-axis, which we consider the major axis of the structure.From this analysis we exclude the sources within a radius of 0.5 cMpc from the center, which roughly corresponds to the virial radius of quasar host halos, resulting in a sample of n = 489 sources.The obtained empirical probability distribution function (ePDF) is then fitted with an exponential PDF(|y|) = e −|y|/λ /λ using EMCEE (Foreman-Mackey et al. 2013) to find the scale parameter λ = 3.33 +0.15  −0.14 cMpc (Fig. 11).We assume a typical log-likelihood log(L) = −n log(λ) − ( n i=1 |y| i )/(nλ) and a flat prior on the scale parameter λ (0.1 < λ < 10).Further, we use the Kolmogorov-Smirnov (KS) test (e.g., Kolmogorov 1933) to determine whether our ePDF is indeed consistent with an exponential.We find that the ePDF is drawn from an exponential function with 95% confidence (KS = 0.0598, p-value= 0.0605) 4 .
The obtained scale width is comparable with the diameter of thick filaments at high redshift measured in cosmological simulations from their DM and gas (e.g., Cautun et al. 2014;Zhu et al. 2021).Intriguingly, thick filaments are expected to live in 4 Initially we fitted our data with the more generic gamma function Γ(k, λ) = (|y| k−1 e −x/λ )/(λ k Γ(k)), which reduces to the exponential when the shape parameter k = 1.We found the shape and scale parameters to be k = 1.14 +0.07 −0.06 and λ = 2.92 +0.22 −0.20 , respectively.The KS test confirms that a gamma function is consistent with the data (KS= 0.0573 and p-value= 0.0806), but given the marginal difference with an exponential we prefer to rely on the fit with only one parameter.overdense areas (see, e.g., Fig. 40 in Cautun et al. 2014), which is reassuring given that on average the targeted fields host overden-sities.However, we stress that the filament thickness is computed in different ways in cosmological simulation (e.g., using the beta model; Zhu et al. 2021) and frequently only focus at the smooth DM and gas mass distribution along the direction perpendicular to the filament spine, while we use halo/galaxy counts.A more robust comparison with cosmological simulations will be the focus of future works.Notes. (a) ID in QSO MUSEUM. (b) Angle north of east for the preferred direction of the overdensity using the rms method. (c) Angle north of east for the preferred direction of the overdensity using the 2D inertia tensor. (d) Major axis of the Lyα nebulae.

The alignment of Lyα nebulae and overdensities
In this section we investigate whether the large-scale Lyα emission on hundreds of kiloparsecs around the targeted systems has a preferred direction with respect to the orientation of the 850 µm overdensities found on megaparsec scales.As the reference angle for the Lyα nebulae, we used their major axis direction, β, obtained from the Stokes parameters weighted by the Lyα SB within the 2σ isophote of each nebulae, as already done in Arrigoni Battaia et al. (2019) for the QSO MUSEUM sample.On the other hand, for the overdensities, we used the preferred directions α rms and α inertia obtained in Sect.6.2.All these angles are listed in Table 6 and the obtained directions are shown in the Figs. 9, I.2, I.3, I.4, and I.5, together with the SB map of each system.
Comparing the angles, we find that the majority of the Lyα nebulae have major axis at less than 45 degrees from the preferred directions of the overdensities.Specifically, 75% and 60% of the sample when comparing to α rms and α inertia , respectively.If we focus only on the fields with |α rms − α inertia | < 30 degrees, we find that 73% of them (or 11 out of 15) have their Lyα nebula major axis at less than 45 degrees from the direction of their overdensity determined by α rms or α inertia .The alignment between Lyα emission and overdensities suggests that the extended Lyα emission traces, for most of the systems, the large-scale structure in which the quasars are embedded.This finding is in agreement with the work by Costa et al. (2022), which studied the extended Lyα emission around high-redshift quasars using cosmological zoom-in simulation post processed with a Lyα radiative transfer code.Indeed, they showed that for quasars host galaxies close to face-on (as it should be more probable for quasars), the Lyα emission is more extended and traces the surrounding large-scale structure (see, e.g., their Fig. 1 for a visual example).On the other hand, more inclined galaxies would hinder the propagation of Lyα photons, resulting in more lopsided nebulae that do not trace the large-scale structure well.
To further test the discovered alignment and measure an average angle between Lyα major axis and the overdensity, we combined the information from all the fields.Specifically, we first consider the catalog detections for each field and rotate them so that the Lyα major axis align with the x-axis.Then, we use all the sources from all fields simultaneously to compute the preferred direction of the overdensity α all rms and α all inertia , as previously done for the individual fields.The left panel of Fig. 12 shows that α all rms =11 +16 −15 and α all inertia =-18 +23 −24 , ultimately proving that, on average, the Lyα nebulae are oriented along the same direction as the overdensities they inhabit.
We also did the same exercise for the separate stacks of rotated ELAN and quasar fields (central and right panels of Fig. 12).This revealed that the ELANe are actually fields that do not have the nebulae orientated along their host large-scale 850µm source distributions but show a preferred relative angle (α ELAN rms

=41 +16
−17 and α ELAN inertia =68 +39 −34 ).This is in contrast to the stack of rotated quasar fields, which clearly proves that these Lyα nebulae are orientated along the large-scale structure (α quasar rms

=-3 +14
−16 and α quasar inertia =-20 +15 −16 ).These results can be understood in terms of the particular evolutionary stage and environment that ELANe showcase.Indeed, several works have reported evidence that ELANe pinpoint the location of protocluster cores given the richness of their environment on small scales and on hundreds of kiloparsecs (e.g., Hennawi et al. 2015;Arrigoni Battaia et al. 2018a;Nowotka et al. 2022).The misalignment between the main large-scale structure and the extended Lyα could therefore be due to the latter tracing inspiraling halo accretion (a.k.a.mergers, interactions) of active galaxies within these protocluster cores (e.g., Arrigoni Battaia et al. 2018b;Chen et al. 2021;Li et al. 2021b).On the other hand, single quasars powering Lyα nebulae likely live in relatively less turbulent and active smallscale environment, and therefore the central quasar is probably better at illuminating their host large-scale structure without contamination of other substructures: a node at the intersection of various filaments.This hypothesis is further strengthen by the difference in the asymmetry profiles and the presence of the secondary overdensity peak for the stacked map for ELAN fields reported in Sect.6.2 (Fig. 10).This likely companion overdensity would be further evidence in favor of more disturbed environments around ELANe.
Further confirmation of our findings requires spectroscopic redshifts for the 850 µm detections.As previously said, this effort is ongoing using sensitive interferometers (Wang et al. in prep.).Also, it is necessary to enlarge the submillimeter surveyed area around these sources to pinpoint the 3D location of the large-scale filaments hosting them, and verify on which scales the overdensity reaches field values.
- Stack of QSO fields rotated by Lya Fig. 12: Distributions of preferred angles (measured with respect to the x-axis) for the stack of all fields (left), ELAN only (central), and quasars only (right), after each individual field has been first oriented with the x-axis along the direction of its Lyα nebula.Both the minimum rms distance and inertia methods result in angles compatible with 0, within their respective errors, for the case of all the field stacks and the case of only the quasar stack.The stack of rotated ELAN fields results in angles that are not compatible with 0.

Summary
We searched for submillimeter sources in the fields of 3 ELANe and 17 z ∼ 3 − 4 quasars (9 from QSO MUSEUM and 8 from MAGG) using the JCMT/SCUBA-2 instrument and surveying an area of ∼ 125 arcmin 2 (or ∼ 450 cMpc 2 ) for each field, with an average central noise of 0.9 mJy beam −1 at 850 µm and 11.7 mJy beam −1 at 450 µm.This work aimed to confirm and extend the effort by Arrigoni Battaia et al. (2018a) and Nowotka et al. (2022), who show that ELANe pinpoint the location of excess 850 µm counts, further indicating that they likely inhabit protoclusters.Our analysis of this large data set includes: -The discovery of 523 and 102 sources (S/N> 4 or detected in both bands at S/N> 3) at 850 µm and 450 µm, respectively.
For each field, we include catalogs in Appendix F. -The calculation of differential number counts at 850 µm and of their underlying count models with self-consistent simulations, which allowed us to compute overdensity factors with respect to blank fields.We find that submillimeter sources with S 850µm ≳ 3 mJy are overabundant in 75% of the targeted fields, by a factor of 2 − 7 within a radius of ∼ 10 cMpc (Sect.5.1).The significance of the overdensity in each field varies between 1σ and 6σ depending on the method used in computing the overdensity factors, as shown by Nowotka et al. (2022).When normalized to the fiducial blank-field counts, the weighted average overdensity factors are 1.9±0.2,1.8 ± 0.1, 2.3 ± 0.2, and 1.6 ± 0.1 for the ELANe, all quasars together, the QSO MUSEUM objects, and the MAGG sample, respectively.When using a cumulative method over the flux density range probed, we find larger values of 3.4 ± 0.4, 2.5 ± 0.2, 3.1 ± 0.4, and 2.3 ± 0.3 for the ELANe, all quasars together, the QSO MUSEUM objects, and the MAGG sample, respectively.This difference is due to underlying models that are steeper than those of blank fields.ELANe therefore show similar overdensities as other quasar fields (see also Sect.6.1).-The finding, at high significance (≳ 5σ), of a radial decrease in the overdensity through the calculation of differential and cumulative number counts in radial apertures using the whole sample (Sect.5.2).Specifically, the overdensity decreases from > 3 within 1.5 arcmin (or about 2 cMpc) to ∼ 2 at the edge of the field.This result suggests that the physical extent of the overdensities around these systems is larger than our maps, in agreement with theoretical expecta-tions that place the effective radius of protoclusters associated with halo masses M DM > 10 12 M ⊙ at scales similar to the size of our maps (∼ 10 cMpc; Chiang et al. 2013;Muldrew et al. 2015).-The comparison of the computed overdensity factors with quasars and Lyα nebula properties, from which we find interesting features in the data that should be explored further (Sect.5.3): (i) z ≥ 4 quasar fields show the smallest source excess, (ii) quasars with similar absolute magnitudes can have a different excess of submillimeter sources, and (iii) being a radio-loud quasar seems to guarantee a 850 µm source excess ≳ 2.8 times that of the blank field.-The discovery of preferred directions for the overdensities of 850 µm sources, which allowed us to stack the fields to unveil an elongated structure reminiscent of a large-scale filament with a scale width of ≈ 4 cMpc (Sect.6.2).We also determine an asymmetry profile for the overdensities.A significant difference is found between ELAN and quasar fields.Overdensities around ELANe show a broken asymmetry profile, with large asymmetries (0.5 to 0.8) for the region enclosing ∼ 50% of the sources (or a radius of ≲ 6 cMpc), and lower asymmetries (∼ 0.3) for outer regions.On the other hand, quasar fields are characterized by a slowly decreasing profile, from large asymmetries (∼ 0.6) in the inner regions to ∼ 0.3 in the outer scales probed.This finding is connected to the presence of a secondary overdensity peak or likely merging large-scale structures in ELAN fields.-The discovery of a trend for the major axis of the Lyα nebulae surrounding the central quasars to align along the direction of the 850 µm overdensities (Sect.6.3).This result suggests that the Lyα emission on scales of hundreds of kiloparsecs traces the central portion of the projected large-scale structure on megaparsec scales.While this is true for quasar fields, ELANe do not have the Lyα nebulae oriented along their host large-scale 850 µm source distributions, likely indicating that the small-scale (hundreds of kiloparsecs) turbulent and active environment of ELANe contaminates the signal from more diffuse gas with a signal associated with active merging substructures.
Overall, our work showcases the richness of ELAN and z ∼ 3 − 4 quasar megaparsec environments at submillimeter wavelengths and paves the way to study the connection between the extended Lyα emission and the surrounding large-scale structure for these massive nodes of the cosmic web.Confirming member associations of the current detections with follow-up observations is ongoing (ALMA and NOEMA).These datasets will allow us to further understand the spatial and kinematic distribution of dusty star-forming galaxies around ELANe and quasars and characterize their cosmic environment.

Appendix D: Spurious sources
For completeness, we report the detailed results for each field and band on the number of spurious sources that could affect our catalogs, as discussed in Sect.3.1.Notes. (a) Estimate of the number of spurious sources at S/N > 4 obtained from the inverted maps. (b) More reliable estimate of the number of spurious sources at S/N > 4 obtained from the jackknife maps. (c) Estimate of the number of spurious sources for the range 3 < S/N < 4 obtained from the jackknife maps.In brackets we list the estimate for the 1000 repetitions (see Sect. 3.1).

Fig. 1 :
Fig. 1: Absolute magnitude at 1450 Å M 1450 of the central AGN versus Lyα nebula areas of the targeted systems (red) in comparison to the parent sample (gray).Each histogram on the side is normalized separately.

Fig. 2 :
Fig. 2: SCUBA-2 S/N maps at 850 µm (left column) and 450 µm (right column) for three example fields: one ELAN (Fabulous), the most overdense quasar field (SDSSJ 0819+0823), and one of the least overdense quasar fields (J020944+051713).The size of each panel is about 15 ′ × 15 ′ .The red star in each panel indicates the location of the ELAN and quasars.The 850 and 450 µm detections with corresponding ID numbers are highlighted as black circles and magenta diamonds, respectively.For each map, we indicate the noise contours (black dashed) for levels of 1.5, 2, and 3 times the central noise, the last of which defines the outer boundary for source extraction (Table 1).The 850 and 450 µm maps for all other fields are shown in Appendix C.
G.2, and G.3  show the mean recovered counts together with their 90% confidence intervals.The simulated counts are in good agreement with the raw counts for all fields.The brightest bins for the MAMMOTH-I ELAN, J0525-2338, and J230301-093930 contain only one or two sources and their counts are not substantially different from the model.These bright sources are

Fig. 3 :
Fig. 3: Differential number counts recovered from 500 realizations of the underlying number count models determined through Monte Carlo simulations (yellow) and the fiducial model (gray) for three example fields: one ELAN (Fabulous) and two quasar fields, the ones with the largest (SDSSJ 0819+0823) and smallest (J020944+051713) overdensity factors in our sample.The shaded regions represent the 90% confidence intervals of the counts recovered from the simulated maps.The black points are raw differential number counts as defined in Sect. 4. The underlying count models and the fiducial model (Geach et al. 2017) are plotted in green and dashed gray, respectively.The effect of flux boosting can be observed as the shaded regions are lying above their input model.Note that the lower zero limits of the fiducial model realizations at faint fluxes reflect the nonuniform noise structure of our maps, where incompleteness increases with distance from the center of the map.Similar plots for all the other fields are shown in Appendix G.

Fig. 4 :
Fig.4: Differential and cumulative number counts obtained in apertures.Top: Comparison of the differential (left panel) and cumulative (right panel) number counts obtained in apertures (in arcminutes; see the legend and Sect.5.2) around the targeted systems with respect to an estimate for blank fields (dashed gray;Geach et al. 2017).The dot-dashed and dotted lines show models with the same γ and S 0 , but with 2× and 4× higher normalization (N 0 ), respectively.The panel in the bottom-right corner shows the ratio between the number of sources for each aperture and the reference blank field model.For a better visualization of the data points at low fluxes, this plot is cut at a ratio of 50, and therefore it does not show the uncertain data points at large fluxes, which by far exceed the blank field model.A trend of decreasing overdensity from the inner to the outer aperture is evident.Bottom: Similar to the top panel, but for the differential number counts obtained in apertures defined in comoving megaparsecs (see the legend and Sect.5.2) around the targeted systems with respect to an estimate for blank fields (dashed gray;Geach et al. 2017), but translated to comoving megaparsecs for the median redshift of our sample.The inset plot compares the number of detected sources with S/N> 4 for the inner bin with the density distribution function obtained from 500 realizations of 20 mock blank fields, similar to what is done in Fig.5with apertures in arcminutes.

Fig. 8 :
Fig. 8: Average differential number counts for ELANe (yellow), MAGG (brown), and QSO MUSEUM (red) quasars shown in comparison with those from blank fields (gray).We indicate the representative uncertainties on these average estimates for quasars and blank fields around the ELANe and Geach et al. counts, respectively.

3Fig. 9 :
Fig.9: Density distribution functions for the direction of the overdensity (left), source distributions with the direction of the overdensity overlaid (center), and Lyα SB maps (right) with the direction of the overdensity and the Lyα major axis overlaid (magenta).In each plot, the direction of the overdensity is color-coded in blue or orange depending on the method used, rms or inertia tensor, respectively (see Sect. 6.2).We show these plots for one ELAN (Fabulous) and for the two quasar fields with the largest (SDSSJ0819+0823) and smallest (J020944+051713) overdensity factor in our sample.In the central panels, the source sizes and colors indicate their completeness, which is used as a weight in the calculation of the direction of the overdensity.The Lyα SB maps are shown with a logarithmic color scheme stretching from 10 −19 to 5 × 10 −17 erg s cm −2 arcsec −2 , and with contours at S/N= 2, 4, 10, 20, and 50.We stress that the Lyα maps are on much smaller scales (see their scale bars) than the 850 µm detections in the central panels.

Fig. 10 :
Fig. 10: Stacked source distribution maps and their asymmetry profiles.Left: Stacked maps of mean number of sources per cMpc 2for all the fields (left), only the ELANe (center), and only the quasar fields (right).For the maps on the top row the preferred direction is given by α rms and the ones on the bottom by α inertia .The black contours in each map enclose increasing fractions of sources from 0.1 for the innermost contour to 1.0 for the outermost one in steps of 0.1.Right: Asymmetry profiles of the stacked source distribution maps.The x-axis gives the radius of the circularized area of the contours, √ area/π, in comoving megaparsecs.From left to right, the points correspond to the asymmetry within the contours enclosing increasing fractions of sources, from 0.1 to 1.0 in steps of 0.1.For the profiles in the top panel, the stacking has been done after rotating each field by its α rms , while for the ones in the bottom panel the fields have been rotated by α inertia .The horizontal dashed black line together with its gray band (in both panels) marks the level of asymmetry of 0.06±0.03corresponding to stacks of all 20 fields where the positions of the sources have been randomized.

Fig. 11 :
Fig. 11: Empirical distribution function for the distances to the x-axis of the sources in the stack of fields rotated according to α rms (gray histogram).The blue curve indicates the most likely exponential distribution that fits the data, while the shaded area represents the 1σ uncertainty on the scale parameter λ.

Fig
Fig. B.2: Same as Fig. B.1, but for eight quasar fields from the QSO MUSEUM survey.Note that for J0525-2338 the jackknife map is affected by the very bright detection associated with this quasar.This should capture the number of spurious sources at high S/N present in the real data and close to this detection.

Fig
Fig. B.3: Same as Fig. B.1, but for nine quasar fields from the MAGG survey (the last target J233446-090812 is also part of the QSO MUSEUM survey).

Fig
Fig. C.4: Same as for Fig. 2, but for one quasar field from the QSO MUSEUM survey, Q2355+0108, and two quasar fields from the MAGG survey, J015741-010629 and J024401-013403.

Fig
Fig. C.6: Same as for Fig. 2, but for three quasar fields from the MAGG survey: J221527-161133, J230301-093930, and J233446-090812 (this last field is also in the QSO MUSEUM survey).
).Each histogram on the side is normalized separately.

Table 2 :
Best-fit parameters for the underlying count models parametrized as Schechter functions.

Table 3 :
Overdensity factors estimated for each field.

Table 4 :
Differential and cumulative source counts computed in areas centered on the targeted systems.
. A follow-up ALMA survey targeting the brightest 850 µm detections in this work is already ongoing (Wang et al. in prep.).

Table 6 :
Preferred directions of the overdensities and major axis of the Lyα nebulae.

Table D .
1: Estimates of numbers of spurious sources in our catalogs.