© The Authors 2023

1. 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 intergalactic 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¹² − 10¹³ 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 h) 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⁻¹ 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¹² − 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 Wide-field 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⁻¹⁷ erg s⁻¹ cm⁻² arcsec⁻²) over > 100 physical kpc, with multiple embedded AGN and active galaxies (e.g., Hennawi et al. 2015; Arrigoni Battaia et al. 2018b, 2022; Chen et al. 2021), 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₀ = 70 km s⁻¹ Mpc⁻¹, Ω_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).

2. Observations

2.1. 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 vs. ∼7 h). 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 (FLASHLIGHT; Arrigoni Battaia et al. 2014), MApping the Most Massive Overdensity Through Hydrogen (MAMMOTH; e.g., Cai et al. 2017b), Quasar Snapshot Observations with MUse: Search for Extended Ultraviolet eMission (QSO MUSEUM; Arrigoni Battaia et al. 2019), and MUSE Analysis of Gas around Galaxies (MAGG; Fossati et al. 2021). The FLASHLIGHT and QSO MUSEUM surveys targeted so far 26 and 61 relatively bright (17.0 < i < 20.0 and 17.4 < i < 19.0) high-redshift (2 < z < 2.3 and 3.03 < z < 3.46) quasars with the aim to uncover surrounding extended Lyα emission and characterize its properties and those of the emitting gas (e.g., Arrigoni Battaia et al. 2016). 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₁₄₅₀ 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.

Fig. 1. Absolute magnitude at 1450 Å M1450 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.

2.2. SCUBA-2/JCMT observations and data reduction

The observations for the twenty fields in this study (Table 1) were conducted with the Submillimetre Common-User Bolometer Array 2 (SCUBA-2; Holland et al. 2013) on the James Clerk Maxwell Telescope (JCMT) during flexible observing under several programs from September 2, 2017, until July 17, 2021 (IDs: M17BP024, M18AP054, M18BP062, M19BP048, M21AP041, M21AP046), under good weather conditions (band 1 and 2, atmospheric opacity at 225 GHz τ_225 GHz ≤ 0.07). For these observations, a Daisy pattern was used to cover an area of ≃13.7′ in diameter, which was centered at the location of each known quasar or ELAN in those fields. To facilitate the scheduling, the observations were organized into scans/cycles of about 30 min. We obtained a total exposure time on source ranging from a minimum of 162 min up to 402 − 423 min for the three ELAN fields.

For the data reduction, we closely followed the method used in Chen et al. (2013a), Arrigoni Battaia et al. (2018a), and Nowotka 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 dimmconfig_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 years¹. 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⁻¹ for 450 μm and 537 Jy pW⁻¹ 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.07 mJy beam⁻¹ and 3.82–20.50 mJy 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.

3. 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 cross-checked 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. Figures 2 and C.1–C.6 show the final S/N maps at 850 and 450 μm for each targeted field with the discovered sources overplotted.

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.

3.1. 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 $\sqrt{t_1\times t_2}/(t_1+t_2)$, where t₁ and t₂ 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.

3.2. 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.

4. 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.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 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

$$\begin{array}{c}\hfill \frac{\mathrm{d}N}{\mathrm{d}S}=\frac{N_0}{S_0}{\left(\frac{S}{S_0}\right)}^{-\gamma }\exp (-\frac{S}{S_0}),\end{array}$$(1)

where N₀ (units of deg⁻²) is the normalization and S₀ (units of mJy) is the characteristic flux. To avoid degeneracy between N₀ and S₀, we keep the value of S₀ 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. Figures 3, 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 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).

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.

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₀ = 7180 ± 1220 deg⁻², S₀ = 2.5 mJy, and γ = 1.5 ± 0.4. Figures 3, G.1–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.

5. Results

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₀, 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 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 selection 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_{-1.0}^{+2.3}$ ($6_{-3.3}^{+5.8}$) sources for fluxes > 30.4 mJy (> 24.0 mJy), which are consistent with the average blank field values from the aforementioned works, $0.6_{-0.6}^{+1.3}$ ($1.5_{-1.4}^{+2.8}$). 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₈₅₀/f₄₅₀ 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.

5.2. 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 one-to-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.

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 S0, but with 2× and 4× higher normalization (N0), 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.

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 μm > 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¹² 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 10 000 (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 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 assuming 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.

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.

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 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.

5.3. 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 rest-frame 1450 Å, M₁₄₅₀ (middle panel), and the deboosted flux at 850 μm, $f_{850}^{\mathrm{Deboosted}}$ 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.

Fig. 6. Comparison of the obtained overdensity factors (Sect. 5.1) vs. some of the quasar’s properties: redshift (top), M1450 (middle), and $f_{850}^{\mathrm{Deboosted}}$ for the quasar or ELAN counterpart (bottom). Magenta dots indicate systems with radio-loud sources.

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.7 < M_i < −23.8 or about −27.4 < M₁₄₅₀ < −22.5 and in 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_{850}^{\mathrm{Deboosted}}\gtrsim 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}_{-17}^{+26}$%) 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 line-of-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), and Fossati et al. (2021). As done in Arrigoni Battaia et al. (2023), we chose the SB threshold 1.23 × 10⁻¹⁷ erg s⁻¹ cm⁻² arcsec⁻² at z = 2.0412 (Jackpot), and obtained the threshold at different redshifts by considering the cosmological SB dimming (1 + z)⁴. 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.

Fig. 7. Cumulative overdensity factor δcumul obtained in Sect. 5.1 vs. 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.

6. Discussion

6.1. Do ELANe sit in similarly overdense fields as other quasars?

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 follow-up ALMA survey targeting the brightest 850 μm detections in this work is already on-going (Wang et al., in prep.).

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.

6.2. 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° (angle α measured from east to north), and computed the rms distance of the sources to each direction². 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 10 000 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.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.

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.

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. We draw 10 000 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.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² 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.

Fig. 10. Stacked source distribution maps and their asymmetry profiles. Left: stacked maps of mean number of sources per cMpc2 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, $\sqrt{\mathrm{area}/\pi }$, 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.03 corresponding to stacks of all 20 fields where the positions of the sources have been randomized.

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³. From the Stokes parameters we can estimate the asymmetry as $(1-\sqrt{Q^2+U^2})/(1+\sqrt{Q^2+U^2})$. 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.03 predicted for randomly distributed sources. An additional 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 deg 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 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 $\lambda =3.{33}_{-0.14}^{+0.15}$ cMpc (Fig. 11). We assume a typical log-likelihood $\log (L)=-n\log (\lambda )-(\sum _{i=1}^n|\mathit{y}{|}_i)/(n\lambda )$ 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)⁴.

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 λ.

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 overdense areas (see, e.g., Fig. 40 in Cautun et al. 2014), which is reassuring given that on average the targeted fields host overdensities. 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.

6.3. 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.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 deg 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 deg, we find that 73% of them (or 11 out of 15) have their Lyα nebula major axis at less than 45 deg 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 ${\alpha }_{\mathrm{rms}}^{\mathrm{all}}$ and ${\alpha }_{\mathrm{inertia}}^{\mathrm{all}}$, as previously done for the individual fields. The left panel of Fig. 12 shows that ${\alpha }_{\mathrm{rms}}^{\mathrm{all}}={11}_{-15}^{+16}$ and ${\alpha }_{\mathrm{inertia}}^{\mathrm{all}}=-{18}_{-24}^{+23}$, ultimately proving that, on average, the Lyα nebulae are oriented along the same direction as the overdensities they inhabit.

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.

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 (${\alpha }_{\mathrm{rms}}^{\mathrm{ELAN}}={41}_{-17}^{+16}$ and ${\alpha }_{\mathrm{inertia}}^{\mathrm{ELAN}}={68}_{-34}^{+39}$). 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 (${\alpha }_{\mathrm{rms}}^{\mathrm{quasar}}=-3_{-16}^{+14}$ and ${\alpha }_{\mathrm{inertia}}^{\mathrm{quasar}}=-{20}_{-16}^{+15}$). 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 small-scale 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.

7. 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² (or ∼450 cMpc²) for each field, with an average central noise of 0.9 mJy beam⁻¹ at 850 μm and 11.7 mJy beam⁻¹ 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 expectations that place the effective radius of protoclusters associated with halo masses M_DM > 10¹² 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–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.

Acknowledgments

We thank the anonymous referee for their careful reading of the manuscript and useful suggestions. We thank Elisabeta Lusso and Ian Smail for providing comments to an earlier draft of this work, and Daniela Galárraga-Espinosa for fruitful discussions on filament sizes estimations from cosmological simulations. The James Clerk Maxwell Telescope is operated by the East Asian Observatory on behalf of The National Astronomical Observatory of Japan; Academia Sinica Institute of Astronomy and Astrophysics; the Korea Astronomy and Space Science Institute; the Operation, Maintenance, and Upgrading Fund for Astronomical Telescopes and Facility Instruments, budgeted from the Ministry of Finance (MOF) of China and administrated by the Chinese Academy of Sciences (CAS), as well as the National Key R&D Program of China (No. 2017YFA0402700). Additional funding support is provided by the Science and Technology Facilities Council of the United Kingdom and participating universities in the United Kingdom and Canada. A.O. is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – 443044596. C.-C.C. acknowledges support from the National Science and Technology Council of Taiwan (NSTC 109-2112-M-001-016-MY3 and 111-2112-M-001-045-MY3), as well as Academia Sinica through the Career Development Award (AS-CDA-112-M02). This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 757535) and by Fondazione Cariplo, grant No 2018-2329.

References

Appendix A: Flux conversion factors

In this work we adopt the recommended FCFs listed on the JCMT web pages⁵. For completeness, we list them here in Table A.1. Our conclusions are unchanged when using the old FCFs values from Dempsey et al. (2013).

Appendix B: S/N histograms

In this appendix we show for completeness the normalized histograms of the S/N values for the pixels within the effective area of the 850 and 450 μm maps for each field. For the 850 μm histograms we indicate the S/N threshold used for the counts analysis (dashed vertical line).

Fig. B.1. Normalized histograms of the S/N values for the pixels within the effective area of the 850 and 450 μm maps. The light and dark gray histograms indicate the distributions of the jackknife maps and of the data, respectively. The data, especially at 850 μm, show excesses at high S/N where sources contribute to the distribution. The dashed vertical line represents the cut used in the simulations of the number counts (see Table 2). The dotted line represents the threshold (4σ) above which we catalog sources.

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. 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).

Appendix C: S/N maps at 850 μm and 450 μm

In this appendix we give the obtained maps for the remaining 17 fields not shown in Fig. 2.

Fig. C.1. Same as for Fig. 2, but for two ELAN fields, Jackpot and MAMMOTH-1.

Fig. C.2. Same as for Fig. 2, but for three quasar fields from the QSO MUSEUM survey: J0525-2338, SDSSJ1209+1138, and SDSSJ1025+0452.

Fig. C.3. Same as for Fig. 2, but for three quasar fields from the QSO MUSEUM survey: SDSSJ1557+1540, SDSSJ1342+1702, and Q0115-30.

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. C.5. Same as for Fig. 2, but for three quasar fields from the MAGG survey: J033900-013318, J111113-080401, and J193957-100241.

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).

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.

Appendix E: Completeness tests

In this appendix we report the figure for the completeness tests of each field and band.

Fig. E.1. Completeness for sources above 4σ as a function of flux for each field and at 850 μm (top) and at 450 μm (bottom).

Appendix F: Source catalogs

In this appendix we report all the source catalogs for the twenty fields in this study, built following the procedure discussed in Sect. 3. Flux boosting and positional uncertainty are computed as done in Nowotka et al. (2022). All catalogs are ordered following the S/N of the sources. For the three ELAN fields, we kept the naming from Nowotka et al. (2022) for the sources in common. Detections at < 10″ from the central system (quasar or ELAN) are indicated by an asterisk next to their name. These sources are likely counterparts to the ELANe and quasars. We stress that caution has to be used when working with these catalogs for follow-up observations as some of these sources could be spurious, especially at 450 μm (Sect. 3). The tables of the catalogs F1 to F40 are available at the CDS.

Appendix G: Differential number counts and Monte Carlo simulations for all fields

In this appendix we show the differential number counts and the results of the Monte Carlo simulations for all the fields not shown in Fig. 3.

Fig. G.1. Same as Fig. 3, but for the two ELAN fields, MAMMOTH-1 and Jackpot.

Fig. G.2. Same as Fig. 3, but for the eight fields from the QSO MUSEUM (see Table 1).

Fig. G.3. Same as Fig. 3, but for the nine fields from the MAGG survey (the last field J233446-090812 is also from the QSO MUSEUM survey; see Table 1).

Appendix H: Different stacked maps for the ELAN fields

In Sect. 6.2 we show that the ELAN fields have a stacked map showing a secondary substructure and that their asymmetry profile decreases abruptly with respect to those of quasars (Fig. 10). This finding is evident when obtaining the stacked maps by rotating the fields using the direction of the overdensities estimated with both the rms distances (α_rms) or the 2D inertia tensor (α_inertia) because the two angles are consistent for these fields. In this appendix we show that the second substructure is due to chance superposition of structures in each field and that the decrease in the asymmetry profile does not change significantly depending on the alignment of the three available fields. In particular, we generated three additional stacked maps by rotating one of the ELAN fields by additional 180 degrees, after its first rotation with α_rms. Figure H.1 shows that the maps are all characterized by secondary peaks, while Fig. H.2 shows that the asymmetry profile for these different maps and the one in the main text are consistent. Therefore, ELAN fields are characterized by dense substructures in the surrounding of the central Lyα nebula, which remain evident in the stacked maps of a few fields.

Fig. H.1. Versions of the stacked ELAN fields maps in Fig. 10, where one of the three fields (indicated in each panel) is rotated by an extra 180 degrees after its first rotation by αrms. These maps demonstrate that the secondary contours in the central panels of Fig. 10 are due to the chance alignment of major substructures in the individual ELAN fields.

Fig. H.2. Asymmetry profiles (dashed, dotted, and dashed-dotted curves) extracted from the contours in Fig. H.1 compared with the corresponding profile for the stacked ELAN fields after rotating the individual fields by αrms. The extra rotation by 180 degrees for one of the three fields results in very similar asymmetry profiles. This means that the shape of these profiles is stable against the rotation of the three ELAN fields.

Appendix I: Preferred direction of the overdensities

For completeness, in this appendix we show the histograms for the calculation of the preferred direction of the overdensities α_rms and α_inertia, together with a visualization of the discovered directions overlaid on the source distribution and Lyα nebula for all the remaining fields not shown in the main text (see Sect. 6.2).

Fig. I.1. Same as Fig. 9, but for the MAMMOTH-1 and Jackpot fields.

Fig. I.2. Same as Figs. 9 and I.1, but for another four sources (see the labels).

Fig. I.3. Same as Figs. 9, I.1, and I.2, but for another four sources (see the labels).

Fig. I.4. Same as Figs. 9, I.1, I.2, and I.3, but for another four sources (see the labels).

Fig. I.5. Same as Figs. 9, I.1, I.2, I.3, and I.4, but for another three sources (see the labels).

All Tables

All Figures

	Fig. 1. Absolute magnitude at 1450 Å M1450 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.
In the text

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.

In the text

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.

In the text

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 S0, but with 2× and 4× higher normalization (N0), 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.

In the text

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.

In the text

	Fig. 6. Comparison of the obtained overdensity factors (Sect. 5.1) vs. some of the quasar’s properties: redshift (top), M1450 (middle), and $f_{850}^{\mathrm{Deboosted}}$ for the quasar or ELAN counterpart (bottom). Magenta dots indicate systems with radio-loud sources.
In the text

Fig. 7. Cumulative overdensity factor δcumul obtained in Sect. 5.1 vs. 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.

In the text

	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.
In the text

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.

In the text

Fig. 10. Stacked source distribution maps and their asymmetry profiles. Left: stacked maps of mean number of sources per cMpc2 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, $\sqrt{\mathrm{area}/\pi }$, 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.03 corresponding to stacks of all 20 fields where the positions of the sources have been randomized.

In the text

	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 λ.
In the text

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.

In the text

Fig. B.1. Normalized histograms of the S/N values for the pixels within the effective area of the 850 and 450 μm maps. The light and dark gray histograms indicate the distributions of the jackknife maps and of the data, respectively. The data, especially at 850 μm, show excesses at high S/N where sources contribute to the distribution. The dashed vertical line represents the cut used in the simulations of the number counts (see Table 2). The dotted line represents the threshold (4σ) above which we catalog sources.

In the text

	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.
In the text

	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).
In the text

	Fig. C.1. Same as for Fig. 2, but for two ELAN fields, Jackpot and MAMMOTH-1.
In the text

	Fig. C.2. Same as for Fig. 2, but for three quasar fields from the QSO MUSEUM survey: J0525-2338, SDSSJ1209+1138, and SDSSJ1025+0452.
In the text

	Fig. C.3. Same as for Fig. 2, but for three quasar fields from the QSO MUSEUM survey: SDSSJ1557+1540, SDSSJ1342+1702, and Q0115-30.
In the text

	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.
In the text

	Fig. C.5. Same as for Fig. 2, but for three quasar fields from the MAGG survey: J033900-013318, J111113-080401, and J193957-100241.
In the text

	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).
In the text

	Fig. E.1. Completeness for sources above 4σ as a function of flux for each field and at 850 μm (top) and at 450 μm (bottom).
In the text

	Fig. G.1. Same as Fig. 3, but for the two ELAN fields, MAMMOTH-1 and Jackpot.
In the text

	Fig. G.2. Same as Fig. 3, but for the eight fields from the QSO MUSEUM (see Table 1).
In the text

	Fig. G.3. Same as Fig. 3, but for the nine fields from the MAGG survey (the last field J233446-090812 is also from the QSO MUSEUM survey; see Table 1).
In the text

	Fig. H.1. Versions of the stacked ELAN fields maps in Fig. 10, where one of the three fields (indicated in each panel) is rotated by an extra 180 degrees after its first rotation by αrms. These maps demonstrate that the secondary contours in the central panels of Fig. 10 are due to the chance alignment of major substructures in the individual ELAN fields.
In the text

Fig. H.2. Asymmetry profiles (dashed, dotted, and dashed-dotted curves) extracted from the contours in Fig. H.1 compared with the corresponding profile for the stacked ELAN fields after rotating the individual fields by αrms. The extra rotation by 180 degrees for one of the three fields results in very similar asymmetry profiles. This means that the shape of these profiles is stable against the rotation of the three ELAN fields.

In the text

	Fig. I.1. Same as Fig. 9, but for the MAMMOTH-1 and Jackpot fields.
In the text

	Fig. I.2. Same as Figs. 9 and I.1, but for another four sources (see the labels).
In the text

	Fig. I.3. Same as Figs. 9, I.1, and I.2, but for another four sources (see the labels).
In the text

	Fig. I.4. Same as Figs. 9, I.1, I.2, and I.3, but for another four sources (see the labels).
In the text

	Fig. I.5. Same as Figs. 9, I.1, I.2, I.3, and I.4, but for another three sources (see the labels).
In the text