ALMA survey of a massive node of the Cosmic Web at z~3. I. Discovery of a large overdensity of CO emitters

Sub-mm surveys toward overdense regions in the early Universe are essential to uncover the obscured star-formation and the cold gas content of assembling galaxies within massive dark matter halos. In this work, we present deep ALMA mosaic observations covering an area of $\sim 2'\times2'$ around MQN01 (MUSE Quasar Nebula 01), one of the largest and brightest Ly-$\alpha$ emitting nebulae discovered thus far surrounding a radio-quiet quasar at $z\simeq3.25$. Our observations target the 1.2- and the 3-mm dust continuum, as well as the carbon monoxide CO(4-3) transition in galaxies in the vicinity of the quasar. We identify a robust sample of eleven CO line-emitting galaxies (including a closely-separated quasar companion) which lie within $\pm 4000\,{\rm km\,s^{-1}}$ relatively to the quasar systemic redshift. A fraction of these objects are missed in previous deep rest-frame optical/UV surveys thus highlighting the critical role of (sub-)mm imaging. We also detect a total of eleven sources revealed in their 1.2-mm dust continuum with six of them having either high-fidelity spectroscopic redshift information from rest-frame UV metal absorptions, or CO line which place them in the same narrow redshift range. A comparison of the CO luminosity function (LF) and 1.2-mm number count density with that of the general fields points to a galaxy overdensity of $\delta>10$. We find evidence of a systematic flattening at the bright-end of the CO LF with respect to the trend measured in blank fields. Our findings reveal that galaxies in dense regions at $z\sim3$ are more massive and significantly richer in molecular gas than galaxies in fields, hence enabling a faster and accelerated assembly. This is the first of a series of studies to characterize one of the densest regions of the Universe found so far at $z>3$.


Introduction
The formation and evolution of galaxies and active galactic nuclei (AGN) is believed to occur within a network of diffuse intergalactic gas (IGM) distributed along filaments and sheets on Mpc-scale structures (the so-called "Cosmic Web"; Bond et al. 1996).Cosmological simulations predict that galaxy formation takes place in the densest regions of this structure where the assembly of galaxies is regulated by the complex interplay between accretion of gas from the Cosmic Web and gas ejection from galaxies into the IGM triggered by feedback mechanisms acting on galactic scales (such as, star-formation and AGNdriven outflows, see, e.g., McNamara & Nulsen 2012;Pike et al. 2014;Wilkinson et al. 2018).However, the details of these processes are still poorly understood.The collection of many lines of evidence over the years led to a growing consensus concerning the key role of the large-scale environment in shaping galaxies during their evolution, at least at z ≲ 1. Elliptical galaxies are preferably found in clusters (Oemler 1974;Davis & Geller 1976;Dressler 1980;Postman & Geller 1984;Dressler et al. 1997;Boselli & Gavazzi 2006) and tend to be passive (with red colors) while galaxies in fields are predominantly blue spirals exhibiting substantial star-formation activity (see, e.g., Lewis et al. 2002;Gómez et al. 2003;Balogh et al. 2004;Hogg et al. 2004;Kauffmann et al. 2004;Tanaka et al. 2004;Park et al. 2007;Peng et al. 2010).The situation at higher redshifts is however unclear.While on the one hand, some results suggest that the aforementioned trend holds up to z ∼ 2 (see, e.g.Postman et al. 2005;Muzzin et al. 2012;Quadri et al. 2012;Darvish et al. 2016;Fossati et al. 2017;Kawinwanichakij et al. 2017;Pérez-Martínez et al. 2023;Taylor et al. 2023), others support the evidence of enhanced star formation activity in protoclusters at z > 2 compared to the field galaxies (e.g., Elbaz et al. 2007;Cooper et al. 2008;Ideue et al. 2009;Tran et al. 2010;Shimakawa et al. 2018;Ito et al. 2020, but see, also, Domínguez-Gómez et al. 2023), which is consistent with the hypothesis that galaxies in dense environments assemble their mass more rapidly and at earlier times (see, Alberts & Noble 2022, for a recent review and further dis-cussion).In order to get further insights on how galaxies form and evolve in connection with their large-scale environment, it is therefore crucial to obtain a comprehensive view of overdense regions of galaxies at early epochs, especially during the peak of galaxy assembly and AGN activity at z ∼ 2 − 3 (Madau & Dickinson 2014).
At high redshifts, there is an increasing contribution from dusty star-forming galaxies to the cosmic star-formation rate density (see, e.g., Casey et al. 2014;Dunlop et al. 2017;Hodge & da Cunha 2020) meaning that a larger fraction of galaxies might remain undetected even in deep rest-frame optical/UV surveys, and possibly even in near-infrared (NIR) observations (e.g.Williams et al. 2019;Yamaguchi et al. 2019;Smail et al. 2021Smail et al. , 2023;;Manning et al. 2022).Such galaxies can however be uncovered in the rest-frame far-infrared (FIR) band where the thermal emission by dust grains dominates the galaxy spectral energy distribution (SED).At such wavelengths, atomic finestructure and molecular emission lines are the major coolants of the gas-phase galaxy interstellar medium (ISM), and they can therefore be targeted to trace the cold gas in galaxies (see, Carilli & Walter 2013, for a review).
With the advent of sensitive facilities in the (sub-)mm such as the Atacama Large (sub-)mm Array (ALMA), astronomers have started to map dense galaxy environments such as (proto-)clusters of galaxies at increasingly high redshift to understand how cold gas -the ultimate fuel of star-formation -and dustwhich is a proxy of the galaxy metal enrichment -are affected in galaxy-rich environments with respect to galaxies living in isolation.Such studies have been mainly conducted by targeting carbon monoxide (CO) rotational lines, the singly-ionized atomic carbon transition [CII] 158 µm , as well as the FIR dust continuum at z ∼ 1 − 2 (e.g.Hayashi et al. 2017;Noble et al. 2017Noble et al. , 2019;;Rudnick et al. 2017;Stach et al. 2017;Coogan et al. 2018;Williams et al. 2022), z ∼ 2 − 3 (e.g., Wang et al. 2016Wang et al. , 2018;;Lee et al. 2017;Castignani et al. 2019;Gómez-Guijarro et al. 2019;Tadaki et al. 2019;Champagne et al. 2021;Jin et al. 2021;Aoyama et al. 2022), up to z ∼ 3 − 4 (see, e.g., Hodge et al. 2013;Miller et al. 2018;Oteo et al. 2018;Umehata et al. 2019;Hill et al. 2020;Polletta et al. 2022).These investigations led to the discovery of numerous gas-rich galaxies in (the core of) galaxy (proto-)clusters.Nevertheless, despite all these efforts, the emerging picture is still unclear and contradictory with tentative evidence of an enhanced star-formation rate, gas and dust fraction, and molecular gas excitation in clustered galaxies, at least at 1 < z < 2.
Crucially, the aforementioned works highlight the importance of sub-mm (pseudo-)blind surveys toward galaxy dense regions at high-z to probe galaxy CO luminosity functions (LFs) and the spectral-line energy distribution, as well as the mm number counts which can provide us with key clues on the physical processes that are acting in the nodes of the Cosmic Web.In this work, we present ALMA observations toward the MUSE Quasar Nebula 01 (MQN01) field.This field hosts a giant Ly-α emitting nebula initially discovered via a blind survey of bright radioquiet quasars (or QSOs) at 3 < z < 4 (Borisova et al. 2016) using the Multi Unit Spectroscopic Explorer (MUSE) mounted on the Very Large Telescope (VLT).The MQN01 nebula surrounding the QSO CTS G18.013 at z = 3.25 is one of the largest nebulae (∼ 30 ′′ corresponding to ∼ 230 physical kpc) discovered in this survey that also exhibits a filamentary morphology.The largest Ly-α nebulae discovered so far are often associated with an overdensity of AGN and massive (dusty) star-forming galaxies (see, e.g., Hennawi et al. 2015;Cai et al. 2017;Cantalupo 2017;Arrigoni Battaia et al. 2018a,b;Umehata et al. 2019;Nowotka et al. 2022).MUSE follow-up observations extending both the covered area and sensitivity limit mapped a large ∼ 4 arcmin 2 area around MQN01, revealing a high concentration of Lyman Break Galaxies (LBGs) embedded in an extended Ly-α emitting structure (Cantalupo et al., in prep., Galbiati et al., in prep.).To obtain a full picture of this exceptional field we used ALMA to perform mosaic observations targeting the mm dust continuum and the CO(4-3) rotational transition in galaxies embedded in MQN01.This work is part of an extensive multiwavelength survey of the MQN01 field which has been conducted in the FIR to the X-rays regime using multiple facilities.Here, we report galaxy detections and field statistics obtained via our mm observations using ALMA.The full characterization of individual sources and their correlation with the Ly-α emitting gas will be presented in future works.
This paper is structured as follow: in Sect. 2 we present our survey design, the acquired observations, the reduction of the data, and the ancillary datasets.In Sect.3, we discuss the source extraction and the measurements of continuum fluxes and line luminosities of the selected candidates.In Sect.4, we present our results which include the analysis of the CO LF, and the mmcontinuum source number count density.We dedicate Sect. 5 to the interpretation and discussion of our results and we compare them with previous works putting our findings into a more general context.Finally, in Sect.6 we draw our conclusions.

Survey design
We observed the MQN01 field with ALMA 12-m array using band 3 and 6 in Cycle 8 (Program ID. 2021.1.00793.S, PI: S. Cantalupo).The observations were designed to cover the entire field of view (FoV) of the MUSE mosaic (∼ 4 arcmin 2 , corresponding to ≃ 16 cMpc 2 at z = 3.25) by performing a Nyquistsampled mosaics following the standard hexagonal pattern to achieve a uniform sensitivity across the entire field.
The band 3 mosaic consists of 27 pointings each with a Half Power Beam Width (HPBW) of ≃ 53 ′′ resulting in a covered rectangular sky area of ≃ 3. ′ 0×3.′ 2. Observations were carried out in the Frequency Division Mode (FDM).The frequency setup consists of four 1.875 GHz-wide spectral windows (SPWs).We tuned two adjacent SPWs in the Upper Side Band (USB) centered respectively at 107.20 GHz and 109.00GHz such that they encompass the CO(4-3) transition (rest-frame frequency ν rest = 461.041GHz), as well as the underlying 3-mm dust continuum in a contiguous redshift bin of ∆z ≃ 0.15 corresponding to ∆ = (−4000, +6100) km s −1 around z = 3.25.We tuned the other two SPWs in the Lower Side Band (LSB) covering a continuous bandwidth in the frequency range 94.06 − 97.74 GHz.The total effective bandwidth of ALMA band 3 observations is 7.354 GHz.The native spectral resolution of the acquired data is 1.95 MHz (∼ 5.4 km s −1 ).Observations were carried out in eighteen Execution Blocks (EBs) during the period 2022 January 16 -May 19 employing a total on source exposure time of 16.4 hours and maximum antenna baseline of 976.6 m.During the executions, the quasars J0025-4803 and J2357-5311 were observed as phase and flux calibrator, respectively.
(1) Central wavelength of the entire frequency setting. (2)Half Power Beam Width (i.e., ∼ 1.13 × λ/D, where λ is the observed wavelength, and D is the ALMA antenna diameter, see Remijan et al. 2019) at the representative frequency of 109.00GHz (band 3) and 245.00 GHz (band 6). ( 3 Observed sky angular area with primary beam sensitivity ≥ 50%. ( 4 Synthesized beam size and position angle (PA) at the representative frequency in case of a natural weighting scheme of the visibilities. (5)Mean Precipitable Water Vapor during observations. (6)The representative bandwidth is 100 km s −1 (corresponding to 36.358MHz) at the representative frequency of band 3, and 6.89 GHz (aggregate continuum) for band 6 observations.jacent 1.875 GHz-wide SPWs that we disposed to cover the 1.2-mm dust continuum together with the CO(9-8) (ν rest = 1036.912GHz), as well as the adjacent transitions of the Hydroxyl Ion OH + (1 1 -0 1 ) (ν rest = 1033.119GHz) that are redshifted in the ALMA band 6 at z ≃ 3.25.The two SPWs in the USB are centered respectively at 243.20 GHz, and 245.00 GHz sampling a contiguous CO(9-8) line redshift bin of ∆z ≃ 0.06, corresponding to ∆ = (−2400, +2100) km s −1 around z = 3.25.The SPWs in the LSB cover the frequency range 228.26 − 231.94 GHz.The effective total bandwidth is 7.35 GHz with a native frequency sampling of 3.9 MHz (∼ 4.8 km s −1 ).The observations were carried out in four EBs during the period 2022 April 3 -13 employing a total on source observation time of 3.3 hours and a maximum baseline of 500.2 m.For such observations, the quasars J0025-4803 and J2258-2758 were used as phase and flux calibrator, respectively.
In Table 1 we report the details of the observations presented in this work.In Fig. 1 we show the combined primary beam (PB) response of the mosaics in the two different bands together with the disposition of the pointings.The last EB of band 6 observations has been affected by increased noise during the last two scans, impacting the sensitivity of the mosaic pointings no.77-114.As a result, the sensitivity in the Northern part of the FoV is about 30% lower relatively to the mosaic center.

Data reduction
We performed data reduction using the Common Astronomy Sofware Application (CASA;McMullin et al. 2007;Hunter et al. 2023).We calibrated the data of both band 3 and 6 observations by running the pipeline scripts (scriptForPI) delivered alongside with the raw Measurement Sets (MSs) by using the CASA pipeline version 6.2.1.In the case of band 6 data, the calibration includes the renormalization correction related to the ALMA amplitude normalization strategy.We imaged the ALMA band 3 and 6 visibilities by using the CASA task tclean by adopting a natural weighting scheme, and by using "mosaic" as a gridding convolution function.We set the phase center at the coordinates ICRS 00:41:32.27 -49:36:19.60.In order to Nyquist sample the longest baselines we set pixel sizes of 0. ′′ 2 and 0. ′′ 15 for band 3 and 6 data, respectively.For each observation in ALMA band 3 and 6, we obtained two sets of "dirty" cubes (i.e., without performing the cleaning process) with a channel width of 25 km s −1 , and 40 km s −1 , respectively in the LSB and USB.During the "cleaning" procedure, we set "cube" as spectral definition mode (specmode) and niter = 0. We also obtained "dirty" continuum images by aggregating all the four SPWs in each frequency setting and we performed the Fourier transform with the tclean by setting specmode="mfs".The resulting beam sizes of the ALMA band 3 data are 1.′′ 4 × 1. ′′ 3, and 1. ′′ 5 × 1. ′′ 4 for data cube at the reference frequency of 109.00GHz, and for continuum image, respectively.The data processing of ALMA band 6 observations yields to beam sizes of 1. ′′ 3 × 1. ′′ 1, and 1. ′′ 2 × 1. ′′ 0 for the continuum image, and the data cube at the reference frequency of 245.00 GHz.After measuring the RMS (root-mean-square) of the "dirty" data and continuum images, we obtained "cleaned" data cubes and continuum images by setting tight circular aper-   tures on ≥ 5σ sources in their continuum and we performed the cleaning using the tclean task down to 2σ (nsigma=2).

Ancillary data
As part of our ALMA program, the presented mosaics have been complemented with an additional single high-resolution (∼ 0.25 ′′ ) pointing encompassing the central region of MQN01 field using ALMA band 3. We designed these observations in order to spatially resolve the CO(4-3) line emission of the quasar host galaxy CTS G18.01, as well as that of possible closelyseparated galaxies.Full details and analysis of this data will be presented in a future paper.In the current work we benefit from this data to disentangle the QSO line and continuum emission which appear blended with a nearby companion located at ∼ 1 ′′ sky-projected distance in the lower-resolution mosaics (see, Sect.3).
Our ALMA observations are part of an extensive multiwavelength observational campaign ranging from the (rest-frame) FIR to the X-ray regime using both ground-based and space telescopes.These data will be presented in full details in future works.In this work we benefit from some of the acquired data.
Here we summarize the main characteristics of such observations.
In this work, we make use of deep VLT/MUSE spectroscopic observations toward the MQN01 field (Cantalupo et al., in prep., and Galbiati et al., in prep.).Such observations consist in four 10 hours-pointing (40 hours of total exposure time) using the MUSE-WFM (Wide Field Mode) integral field spectrograph with adaptive optics covering a FoV of ≃ 2 ′ × 2 ′ which is fully sampled by our ALMA data.The MUSE observations provide us with integral field spectroscopy between 4650 − 9300 Å with a spectral resolution of R ∼ 2000 − 3500.
In this field, galaxies are identified via their rest-frame far-UV (FUV) continuum emission in the white-light image.The final sample comprises only those sources with a high-confidence measurement of the spectroscopic redshift from interstellar absorption lines (e.g., SiII 1260, 1526Å, OI 1303Å, CII 1334Å, and CIV 1548, 1550 Å).The catalog includes 38 galaxies with secure redshift between 2.7 < z < 4.3, out of which 22 lie within ±1000 km s −1 with respect to the QSO CTS G18.01 systemic redshift (z = 3.25; Galbiati et al., in prep., for full details).
We imaged the MQN01 field in the optical by using the VLT/FORS2 (Focal Reducer and low-dispersion Spectrograph 2; Appenzeller & Rupprecht 1992) instrument by acquiring broadband mosaics (≃ 13×13 arcmin 2 ) in U (central wavelength λ 0 = 361 nm; ∼ 7 hours of exposure time), B (λ 0 = 440 nm; ∼ 5 hours of exposure time), and R (λ 0 = 655 nm; ∼ 45 minutes of exposure time) filters.Such observations cover a 36× larger sky area with respect to the MUSE mosaic.This allows us to extend the census of the z = 3.0 − 3.5 population of LBGs well beyond the MUSE FoV via color-color selection tested against the spectroscopic information available in the MUSE FoV (Galbiati et al. in prep.).
In this work, we additionally benefit from NIR photometric images taken with NIRCam (NIR Camera) instrument on board the James Webb Space Telescope (JWST; Rigby et al. 2023).Observations were acquired by using the extra-wide filters F150W2, and F322W2 in the short-and long-wavelength channel, respectively, with an on-source exposure time of 1632 sec per filter.These images cover a total FoV of 2 × 5 arcmin 2 across the two detectors of the camera, with one detector encompassing the sky area observed with MUSE.

2.8740
Notes.Upper part: high-fidelity (F ≥ 0.9) candidate list extracted via blind search using LineSeeker.Lower part: source candidates extracted by cross-matching ALMA 1.2-mm continuum-selected sources with F > 0.2 with the MUSE catalog of z ≥ 2.5 sources.
(1) Identifier of ALMA 1.2-mm continuum-selected candidates.In case a source is also detected in its CO(4-3) line (see, Sect.3.2), the corresponding identifier from Table 3 is reported.In case a source has been identified as low-z interloper, its identifier is also reported (see, Sect 3.5).The QSO and the nearby companion are labeled with the Calán-Tololo Survey (CTS, Maza et al. 1995) identifier, and with "Object B", respectively.Object B has been selected by LineSeeker together with the QSO as a single source.Here we report, separately the coordinate of the continuum peak of the QSO CTS G18.01 and Object B obtained via a visual inspection.Their S/N and the fidelity are set to the values provided by the code that are referred to the continuum peak of the Object B which appears brighter than the QSO at 1.2 mm. (2)Right Acension (ICRS). ( 3 Declination (ICRS). ( 4 Integrated flux density over the ≥ 2σ isophote.This quantity is not reported for compact sources at low S/N and for the QSO and Object B which are partially blended. (5)Continuum peak flux density at 1.2 mm. (6)Continuum peak flux density at 3 mm. (7)Signal-to-noise ratio of the 1.2mm-selected sources. (8)Fidelity and its uncertainties estimated by LineSeeker using negative detections. (9)High-confidence redshift estimate from MUSE spectroscopic data.For those sources for which the redshift is not provided they either have an uncertain redshift measurement or they do not have any FUVcontinuum counterpart from MUSE. ( †) Source C06 has been unambiguously identified as a low-z interloper ("i1") showing a bright emission line in the LSB of 3-mm ALMA band.For this galaxy we therefore report the low-confidence redshift provided by the analysis of the MUSE spectrum.

Source search and characterization
We performed a blind search of continuum-and line-emitting sources in both our ALMA band 3 and 6 images and cubes by using the Python-based code LineSeeker 1 (see, González-López et al. 2017b, 2019, for full details).This code was originally developed to search for line and continuum emission of galaxies in the ALMA Frontier Fields survey (see González-López et al. 2017a,b), and was subsequently employed in various other surveys such as the ALMA Spectroscopic Survey in the Hubble Ultra Deep Field (ASPECS; see, e.g., Decarli et al. 2019;González-López et al. 2019, 2020), the Multiwavelength Study of ELAN Environments (AMUSE2 ; see, e.g.Chen et al. 2021;Arrigoni Battaia et al. 2022), and the Northern Extended Millimeter Array (NOEMA) Molecular Line Scan of the Hubble Deep Field North (Boogaard et al. 2023).Here we summarize the operation and the main features of the code.LineSeeker adopts a matched-filter approach.The code combines adjacent spectral channels by convolving the data cube along the spectral axis using Gaussian kernels with a range of widths.The RMS of the resulting images is then estimated via a sigma clipping at 5σ to remove the spurious increases of the noise due to possibly bright lines or continuum emission within the convolved channels.Then all the voxels above a given signalto-noise ratio 2 (S/N) threshold are stored for each convolution kernels.Finally, a list of line (or continuum3 ) emitter candidates is generated by grouping voxels from the different channels cor-responding to unique sources by using the Density-Based Spatial Clustering of Application with Noise algorithm (DBSCAN; Ester et al. 1996) included in the Python package Scikit-learn (Pedregosa et al. 2011).The final S/N of the sources is then selected as the maximum value obtained through the different convolutions.DBSCAN is also able to recover extended sources traced by S/N ≥ 2 pixels, however, a visual inspection is needed in order to verify if the single extended source is actually composed by multiple blended sources along the line of sight.
For each source candidate selected by LineSeeker, the code automatically estimates the probability of false-positive detection based on the source S/N.To this purpose, the code is run on the negative (i.e., multiplied by −1) cube or image.In a deep extragalactic field in the (sub-)mm, the majority of the surveyed area is expected to be empty sky; hence, any "negative" peak is a realization of noise.The statistics of negative detections is then compared to that of the positive ones.The fidelity (or reliability) of a positive detection as a function of its S/N is therefore computed as where N neg and N pos are the number of negative and positive detections at a given S/N, respectively.To compute the fidelity at any S/N following Eq.(1), LineSeeker assumes that the noise in the data is Gaussian distributed and compute the best fit model of the cumulative distribution of negative detections using a function of the form N[1 − erf(SN/ √ 2 σ)], with erf that is the error function and N, σ are free parameters.
Other similar source-finding algorithms are available in the literature such as FindClump (Walter et al. 2016), and MF3D (Pavesi et al. 2018), that mostly differ on details (such as, e.g., the adopted spectral filter function).Comparisons between the codes yield consistent results to within ∼ 10% (see, González-López et al. 2019).

1.2-mm continuum-selected candidates
We performed a source search of 1.2-mm continuum-emitting galaxy candidates in MQN01 field by running LineSeeker on the "dirty" continuum band 6 image, excluding the region with PB response < 50% in which the low telescope sensitivity enhances the fraction of spurious candidates.The "dirty" data are preferred over the "cleaned" ones since in the former the intrinsic properties of the noise are preserved.Also, we do not correct our "dirty" image for the PB response to preserve the spatial homogeneity of the noise across the FoV.We therefore extracted all S/N ≥ 3 detections, and we selected the source candidates with estimated fidelity of F ≥ 90% corresponding to S/N ≥ 4.7.With this method, we retrieve a total of nine sources including the QSO CTS G18.01, and a closely (on-sky) separated source (hereafter, Object B) partially blended with the QSO.
We complemented our blind search of 1.2-mm continuum candidates in the field by cross-matching our MUSE catalog of high-z sources (z > 2.5, see Sect.2.3) with the low-fidelity (F > 20%) sources selected by LineSeeker.In this process, we cross-matched the on-sky spatial position of the MUSE and ALMA continuum sources within a separation limit of 0. ′′ 6.We chose this separation since it is about one half of the angular resolution of the ALMA image thus accounting for possible spatial offset between the ALMA low-S/N FIR-and MUSE FUVcontinuum peak4 .This separation also corresponds to the maximum observed angular distance between our ALMA candidates selected in the blind search and their MUSE counterparts.By doing so, we recovered two additional sources in the field.In Table 2 we report the final catalog of the eleven ALMA 1.2mm continuum-selected sources.In Fig. 2, we show the location of the sources detected in MUSE within ±1000 km s −1 with respect to QSO CTS G18.01 and the ALMA 1.2-mm continuumselected sources in the MQN01 field.We labeled the latter as C01 -C09.

CO(4-3) line-emitting candidates
We used LineSeeker to blindly search for CO(4-3) emission lines that are expected to be redshifted in the USB of the ALMA band 3 datacube.We opted to extract sources in the "dirty" cube not corrected for the PB response over the area where the combined mosaic sensitivity is ≥ 50% (see also, Sect.3.1).We run the line-search algorithm on the datacube that we spectrally binned at 25 km s −1 using Gaussian kernels with widths ranging from 0 to 18 channels.This range enables the code to match the typical linewidths of CO lines observed in high-z galaxies (FWHM ∼ 50 − 1000 km s −1 ; see, e.g., Carilli & Walter 2013).We extracted all line-emitting candidates with estimated S/N ≥ 3. We then selected those sources with estimated fidelity of F ≥ 90% (see, Sect.3).Similarly to the search of 1.2-mm continuum emitters, we cross-matched the location of the MUSE LBGs with the low-fidelity (F > 20%) line-emitting candidates from LineSeeker within a separation limit of 0. ′′ 6.As a result of this procedure, we extracted one additional source.The CObased redshift of this candidate is within ±200 km s −1 with that based on the Ly-α emission line derived from MUSE spectro-scopic data.The difference is consistent with the typical shift observed between the Ly-α line and the systemic redshift of high-z Ly-α emitters (LAEs; see, e.g., Guaita et al. 2013;Muzahid et al. 2020aMuzahid et al. ,b, 2021)).Within the selected sample of 13 galaxies we identified two interlopers at lower redshift (dubbed as "i3" and "i4") which we therefore excluded from the final sample (see, Sect.3.5, for a detailed discussion).In Table 3 we report the final catalog of the eleven ALMA CO(4-3) line-emitting candidates.In Fig. 2, we draw the location of these galaxies (L01 -L09) as well as that of the identified low-z interlopers within the field.

3-mm continuum-selected candidates
Similarly to what described in Sect.3.1, we complemented our source extraction by searching candidates detected in continuum at 3 mm.For this purpose, we run LineSeeker on the "dirty" aggregate (LSB+USB) continuum ALMA band 3 image.This blind search results in six continuum-detected sources (including the QSO and Object B) with fidelity F ≥ 90% .We also carefully inspected the ALMA 3mm-continuum image in the position of secure sources (F = 100%) revealed in the ALMA band 6 at 1.2 mm.We therefore included two additional sources in the final sample corresponding to C04 and C05 (see, Table 2) which show convincing continuum emission at 3 mm.We verified our conclusion by cross-matching the catalog of secure positive continuum detections at 1.2 mm with that at 3 mm provided by Line-Seeker.Finally, the cross-match between the catalog of high-z MUSE LBGs and 3 mm-continuum candidates did not provide us with any additional source.The locations of these detections in the MQN01 field are indicated in Fig. 2 by green squares.

Completeness and Flux Boosting
In order to compute the source number counts and the CO(4-3) LF it is crucial to determine the probability of detecting sources in our blind search.Such information is enclosed in the continuum and CO line completeness functions.In addition, we need to estimate how the measured fluxes are affected by the noise in the real data.Indeed, sources with low S/N have higher probability to be recovered with higher flux than the intrinsic value because of the boosting effect produced by the noise fluctuations (the so-called "flux-boosting effect"; see, e.g., Hogg & Turner 1998;Scott et al. 2002;Coppin et al. 2006).We computed the completeness of our survey by following a common approach widely adopted in the literature (see, e.g., Hatsukade et al. 2016Hatsukade et al. , 2018;;Umehata et al. 2017;González-López et al. 2019, 2020;Béthermin et al. 2020;Boogaard et al. 2023).We injected artificial sources in the real data and we then performed the source search by using LineSeeker.More details about this exercise and how it is used to perform the corrections to the LFs are described in the following.
To estimate the completeness function of 1.2 mm-continuum detections in the ALMA band 6 we created artificial point-like sources by rescaling the synthesized beam model and we injected them into the dirty ALMA Band 6 continuum map at random positions within the source-search area of our survey (i.e., where the combined mosaic PB response is ≥ 50%).To take into account the effect of the variation of sensitivity across the mosaic field, we rescaled the source fluxes by the PB response at the input location of each source.We then run the source-search algorithm and we computed the source detection rate.We considered a source recovered if it is extracted within 1 ′′ from its input location with a fidelity ≥ 90%.We repeated this procedure Article number, page 7 of 25 A&A proofs: manuscript no.mqn01-I_arXiv Table 3. ALMA CO(4-3) line-selected sources.Notes.Upper part: high-fidelity (F ≥ 0.9) candidate list extracted via blind search using LineSeeker.Lower part: source candidates extracted by cross-matching ALMA CO(4-3) line-emitting candidates with F > 0.2 with the MUSE catalog of z ≥ 2.5 sources.
(1) Identifier of ALMA CO(4-3) line-emitting (1.2-mm continuum-selected) candidates.The QSO and the nearby companion (Object B) are labeled as in Table 2. Object B has been selected by LineSeeker together with the QSO as a single source.Here we report, separately the coordinate of the continuum peak of the QSO CTS G18.01 and Object B obtained via a visual inspection.Their S/N and the fidelity are set to the values provided by the code that are referred to the CO(4-3) peak of the QSO.
(2) Right Acension (ICRS). ( 3 Declination (ICRS). ( 4 Our best estimate of the CO(4-3) flux derived from the Gaussian fit of the line emission (see Sect 3.6). ( 5) CO(4-3) luminosities estimated via Eq.( 2) and (3), respectively. (8)Maximum value of the Signal-to-noise ratio obtained through the different convolution (see Sect. 3). ( 9 Fidelity and its uncertainties estimated by LineSeeker using negative detections. (10)Redshift estimate from the CO line centroid. (11)Velocity shift with respect to the QSO CTS G18.01 redshift. (12)High-confidence redshift estimate from MUSE spectroscopic data.For those sources for which the redshift is not provided they either have an uncertain redshift measurement or they do not have any FUV-continuum counterpart from MUSE.
by injecting 20 sources simultaneously in the image and iterating for 100 times for each 1.2-mm flux density value within the range 0.02 − 0.46 mJy in steps of 0.02 mJy.The total number of injected sources in the simulation is 46000.During the simulation we also evaluated the effect of flux boosting by computing the ratio of measured-to-injected flux of the artificial sources as (F meas 1.2 mm − F inj 1.2 mm )/F inj 1.2 mm .In Fig. 3 we report the output of our simulation.As a result, the completeness in the flux range of the detected sources is in the range 50 − 100%, while the flux boosting effect does not affect significantly the measured 1.2-mm flux density of our selected sample of galaxies.
We evaluated the completeness of our blind survey of lineemitting sources and the effect of the flux boosting on the observed line emission.On the basis of the line measurements of the selected CO(4-3) emitters, we injected point-like artificial emission lines using Gaussian spectral profile with FWHM and velocity-integrated flux ranging within 100−1000 km s −1 in steps of 100 km s −1 , and 0.02−0.5 Jy km s −1 in steps of 0.02 Jy km s −1 , respectively.We injected the artificial sources at random positions (within the volume where the mosaic PB response is ≥ 50%) in the USB of the dirty ALMA band 3 cube binned at 25 km s −1 (i.e., where we search for the CO(4-3)-line emitting candidates; see Sect.3.2).We then rescaled the signal from the artificial line emitters in each channel of the cube by the PB response at the input position of the sources.We then ran Line-Seeker and we computed the source detection rate following the same criteria adopted for the artificial sources in the continuum.To estimate the line flux boosting effect for the recovered sources, we obtained the line-velocity integrated map using channels within ±2σ with respect to the line centroid provided by LineSeeker5 .We then measured the total source fluxes by performing fit of the moment-0 map by using a 2D Gaussian model6 .For each values of the line FWHM and line-velocity integrated flux, we injected simultaneously 50 sources in the cube measurements of our selected sample of CO(4-3) line-emitting galaxies.

Identification of interlopers across the redshift range
For the purpose of our analysis it is crucial to identify possible interlopers at different redshifts among the ALMA-selected galaxies.In the ALMA 3-mm band, sources at any redshift are expected to be revealed in their line emission primarily via CO rotational lines, or the fine-structure transitions of the atomic neutral carbon [CI] 609, 370 µm (albeit the latter are expected to be much fainter than low-J CO lines; see, e.g., González-López et al. 2019;Decarli et al. 2020, but see, also, Gullberg et al. 2016).In Fig. 5, we draw the redshifted frequency of such transitions entering in our ALMA band 3 SPWs up to z = 6.The various transitions probe galaxies within different redshift intervals and cosmological volumes depending on the surveyed sky area and the encompassed frequency range (see, Table 4).By using the available best-fit model to the CO LFs from Boogaard et al. (2023), we computed the number of expected galaxies at various redshifts within the cosmological volume probed by our observations above the limiting luminosity (see, Table 4).As a result, ∼ 2 sources at z ∼ 1.1 are expected to be detected via their CO(2-1) line in the USB of our ALMA band 3 survey, ∼ 1 at z ∼ 2.2 through CO(3-2) line, and ∼ 0.5 CO(5-4) line emitters at z ∼ 4.3.For the targeted CO(4-3) line the number of expected sources is ∼ 2.5, dropping to ∼ 0.5 when considering the volume within ±1000 km s −1 around the QSO CTS G18.01.However, these estimates suffer from large uncertainties given the poor sampling of the CO LFs (see, e.g., Decarli et al. 2019Decarli et al. , 2020;;Boogaard et al. 2023), and therefore must be taken with caution.In order to unambiguously identify the emission lines observed in our survey, we made use of our multiwavelength datasets of the MQN01 field.As described in Sect.3.2, our blind search for CO(4-3) line-emitting candidates yielded an initial sample of 13 sources.Within the candidates belonging to this sample, 4 of out 13 (≃ 30%; L03, L07, L08, i4, excluding Object B) do not correspond to any high-z MUSE-selected LBGs/AGN.In addition, three candidates (L05, L06, i3) are located outside the MUSE mosaic footprint.For the first group of galaxies, the analysis of their MUSE spectra enabled us to identify i4 as a low-z interloper.The spectrum of this galaxy shows MgII] 2796, 2803 Å absorptions, and [OII] 3726, 3729 Å line emission, the wavelengths of which unambiguously place this object at z ≃ 1.170, hence the observed line emission in our ALMA band 3 observation is the CO(2-1) at the redshifted frequency of 106.2 GHz.In the case of sources outside the area surveyed by MUSE, the analysis of UBR colors obtained from VLT/FORS2 observations (see, Sect.2.3; Galbiati et al., in prep., for full details) allowed us to identify, at high confidence, i3 as another one low-z galaxy interloper.The line emission of i3 in our ALMA band 3 is observed at 109.2 GHz, the frequency of which possibly corresponds to either the CO(2-1) at z ≃ 1.11 or CO(1-0) z ≃ 0.06.Future spectroscopic follow-up observations will possibly pin-down the precise redshift of i3.As a result of these analysis, we excluded the aforementioned two sources from the final sample of CO(4-3) line-emitting candidates.
(1) CO and [CI] emission lines entering in our ALMA band 3 observations. (2)Redshift range in which the line can be detected.
(3) Cosmological volume in the redshift range within the effective area of our survey (PB ≥ 0.5). ( 4 Limiting 5σ line luminosity assuming FWHM = 300 km s −1 . which belongs to the sample of the 1.2-mm continuum-selected candidates (C06), while the other one ("i2") which is detected in its 3-mm continuum emission, both showing a very bright emission line in the LSB of the ALMA band 3 datacube.Via the analysis of the MUSE spectra, we accurately determined the redshift of i2 to be z ≃ 1.42, therefore identifying the detected line as the CO(2-1).Regarding i1, its redshift determination is highly uncertain given the lack of clear absorption or emission features in the MUSE datacube.However, the carefully inspection of the rest-frame UV spectrum of i1 points to a tentative redshift measurement of z ≃ 2.54, thus suggesting that the emission line we detected in our ALMA observations corresponds to the CO(3-2) transition.
We then assessed the nature of our ALMA continuum selected sources at 1.2 mm by inspecting their MUSE spectra.Among these, 6 out of 11 (QSO, C01, C02, C06, C08, and C09) have high-confident spectroscopic redshift derived from their rest-frame UV spectra, while 5 out of 11 (C03, C04, C05, and C07) are either not detected in MUSE or the available spectrum does not allow to pin-down a precise spectroscopic redshift.For the first group, the spectroscopic information available places four of them in the proximity of the QSO CTS G18.01 redshift, with three of them (QSO, C01, C02) having CO(4-3) detection, and one (C08) which is only detected in its 1.2-mm dust continuum.The other two sources for which redshift measurements from MUSE are available are C06, and C09.C06 is the 1.2mm continuum counterpart of the interloper i1, possibly detected via the CO(3-2) in the LSB of our ALMA band 3 observations.The MUSE spectra of C09, instead, exhibits SiIV 1394, 1403 Å, and CIV 1548, 1550 Å absorption lines thus determining that this source is an interloper located at z ≃ 2.874.At this redshift, we do not expect to detect any bright emission line in our ALMA datacubes.We dubbed this source as "i5".
In summary, among all the sources extracted with high fidelity from our ALMA data, we unambiguously identified five interlopers located to different redshifts with respect to the QSO, namely i1 (z ≃ 2.54), i2 (z ≃ 1.42), i3 (either z ≃ 1.11 or z ≃ 0.06), i4 (z ≃ 1.170), and i5 (z ≃ 2.874).Interestingly, these findings are consistent to our predictions based on the CO LFs of blank fields.
However, the subsample of secure sources with two independent redshift measurements from both MUSE and the CO(4-3) line is composed by QSO, L01 (C01), L02 (C02), L04, and L09, all lying within ±1000 km s −1 with respect to the QSO systemic redshift.Regarding the remaining sources, both Object B and L03 (C04) are either revealed in the mm dust continuum with ALMA or have a counterpart in the optical/NIR (see, Appendix A).The analysis of VLT/FORS2 UBR colors of L03 (C04) suggests that this source belong to z ∼ 3 − 3.5 thus supporting the conclusion that L03 is actually detected via CO(4-3) at z ≃ 3.25.Finally, the recently acquired spectrum of the Object B with the NIR spectrograph on board of JWST definitely confirms that this source is located in the proximity of the QSO (Pensabene et al., in prep.).On the other hand, L05, L06, L07, and L08 are only detected in line with ALMA.Therefore, we cannot rule out that the latters sources actually represent either false-positive detections or interlopers located at different redshift.Future deep NIR observations are needed to assess the nature of such objects.Interestingly, these sources exhibits a large velocity shift and significant spatial separation from the QSO.In what follows, we took the aforementioned consideration into account in the computation of the CO LF.

Source fluxes and luminosities
The majority of CO(4-3) and continuum emitters detected in this work appear as compact spatially-unresolved sources.For such objects, the total continuum or line flux can be measured through a standard single-pixel analysis of the data.However, in the case of partially-resolved objects or extended sources with complex morphology, this simple method yields to significant flux underestimation.In this work, for such sources we therefore performed source flux measurements by applying the 2σ-clipped photometry7 (see, e.g., Béthermin et al. 2020) as described below.
We measured the 1.2-and 3-mm flux of the sources from the cleaned ALMA band 6 and 3 continuum images, respectively, not corrected for the combined mosaic PB response.For each source, we sum the flux density in mJy beam −1 enclosed in the contiguous area around the source including all pixels with S/N ≥ 28 .We then divided the total flux density per beam by the synthesized beam area (in pixel units) and we rescaled the flux for the PB response at the location of the source.We computed the flux uncertainty by rescaling the noise by the square root of the number of independent beams within the integration area.In Table 2 we report our source flux measurements as well as their peak flux in the continuum.To understand which sources can be considered point-like, we computed the uncertainty-normalized difference between the two flux estimates as ∆ F = (F 2σ − F peak )/ σ 2 2σ + σ 2 peak , where F 2σ , F peak , σ 2σ , and σ peak are the 1.2-mm continuum flux obtained via the 2σ-clipped photometry, the peak flux, and their uncertainties, respectively.Within this formalism, we expect sources which are significantly resolved to have ∆ F > 1.Sources C01 and C02 are such cases; for them we therefore adopted S 2σ 1.2 mm as our fiducial 1.2-mm continuum flux density measurement 9 .We however note that the proper flux measurement of QSO CTS G18.01, and its closely-separated companion (Object B) is challenging due to the partial blending of the sources at the current resolution of the continuum data (see Sect. 2.2).In order to minimize the flux contamination, for such sources we adopted their 1.2-mm continuum flux peak F peak 1.2 mm .We measured the CO line fluxes by following the iterative process presented in Béthermin et al. (2020).As for the continuum, the CO line emission of the QSO CTS G18.01 and Object B are partially blended.To measure their fluxes in what follows, we employed the ALMA band 3 high-resolution (∼ 0.25 ′′ ) data (see, Sect.2.3) where the line emission of the sources are spatially-resolved.For each CO(4-3)-line emitter we extracted the spectrum at the peak position of the source from the "cleaned" ALMA band 3 datacube binned at 40 km s −1 .We then scaled the source spectrum by the PB response, and we performed a fit using a Gaussian profile for the line and a constant for the underlying 3-mm dust continuum using the curve_fit task included in the SciPy package (Virtanen et al. 2020).We then produced the line-velocity integrated map using all the channels within ±2σ from the line centroid.Subsequently, we re-extracted the source spectrum by summing all the spectra in pixels showing S/N ≥ 2 in the moment-0 map, and we rescaled the channel fluxes and their uncertainties respectively by the synthesized beam area, and the square root of the number of independent beams within the integration area.In this process, we masked all the pixels below the chosen threshold which we confidently believe to be not related to any real emission from the source.This new spectrum is more informative in the case of (partially-)resolved sources since it includes the signal from the entire source line-emitting region.We therefore performed a new spectral fit using the best-fit parameters of the previous iteration as starting point for the fitting code.We repeated the aforementioned steps a few times until convergence.All the extracted source spectra remain stable after a few (< 10) iterations.Similarly to continuum flux estimates discussed above, we computed the quantity ∆ F for each source.As a result, L01, L02, the QSO, and the Object B10 all exhibit ∆ F > 1. Accordingly, this analysis determined our final source spectra.We finally performed a finer fit of the final source spectrum by sampling the parameter space through the Python Monte Carlo Markow Chain (MCMC) ensemble sampler emcee (Foreman-Mackey et al. 2013).We employed flat priors on the basis of the best-fit parameters derived from the last iteration and we assume Gaussian uncertainties in the definition of the likelihood.We finally derived the line luminosities as (see, e.g., Solomon et al. 1997): where S ∆ is the velocity-integrated line flux in Jy km s −1 , ν obs is the observed central frequency of the line in GHz, z is the source redshift measured from the CO line centroid, and D L is the luminosity distance in Mpc.The relation between Eq. ( 2) and ( 3) is L CO = 3 × 10 −11 ν 3 rest L ′ CO , where ν rest is the line rest frequency in GHz.
As our observations probe the (rest-frame) FIR wavelengths of high-z galaxies, the Cosmic Microwave Background (CMB) might have an impact on our continuum and line measurements by increasing both the dust and the line excitation temperature.In addition, the CMB provides a strong background signal (see, da Cunha et al. 2013;Zhang et al. 2016).
By assuming a modified black body with typical values for SMGs at z ∼ 3 (i.e., spectral index 1.6, and dust temperature of 35 K; see, e.g., Kovács et al. 2006), the CMB affects our 1.2 mm dust continuum measurements by ≲ 5%, thus compatible with the typical uncertainties on our flux estimates.The CMB effect further decreases assuming higher dust temperature.
Regarding the CO(4-3) line measurements, under the assumption of local thermodynamic equilibrium (LTE; i.e., the line excitation temperature equaling the gas kinetic temperature), the effect of CMB in reducing the recovered line fluxes is always < 20% for any excitation temperature > 40 K at z ≲ 4 (see, Decarli et al. 2020).However, due to the unknown excitation temperature of the gas and the unverified assumption of LTE for our sources, in this work we opted to not apply the CMB-related corrections (see, e.g., Decarli et al. 2020;Boogaard et al. 2023).
We list the measured fluxes and derived quantities of the lines in Table 3, and we report the final source spectra in Fig. 6.In Fig 7 we report the source maps of the combined sample sources within ∆ QSO < 4000 km s −1 either detected in their CO(4-3) line or continuum at 1.2 mm.

CO Luminosity Function Analysis
We computed the CO(4-3) LF by following the approach described in Decarli et al. (2016Decarli et al. ( , 2019Decarli et al. ( , 2020)); Riechers et al. (2019); Boogaard et al. (2023).We define the CO LF as where Φ is the number of sources per comoving Mpc 3 in the luminosity interval log L ′ CO ± 0.5 ∆(log L ′ CO ), ∆V is the comoving volume of our survey, and F i and C i are the fidelity and completeness associated to each source, respectively.We computed two different CO LFs, one in ∆V corresponding to the redshift range within ±4000 km s −1 with respect to the systemic redshift of QSO CTS G18.01 (∆ QSO ), and another one within ∆ QSO < 1000 km s −1 .The corresponding cosmological volumes are respectively ∆V 4000 = 2395 cMpc 3 , and ∆V 1000 = 599 cMpc 3 .While ∆V 4000 contains all our eleven CO(4-3) lineemitting candidates, ∆V 1000 encompasses only sources having a optical/NIR counterparts (see also, Appendix A) with a spectroscopic-confirmed redshift, thus including secure sources (F = 1).In this regard, we use sources in ∆V 1000 to obtain a "raw" CO LF without applying the completeness correction.
To obtain the CO LF in the MQN01 field, we performed a Monte Carlo simulation of 1000 independent realizations of the LF.In each iteration, we varied the CO line luminosity, line FHWM, and the source fidelity of our CO line emitters within their uncertainties.For secure sources, we fixed the fidelity value to F = 1, for all the other sources we treated the fidelity as upper limit.This approach provides a conservative treatment of our fidelity estimates attempting to include the systematic uncertainties associated with sources without any clear multiwavelength counterparts (see, e.g., Pavesi et al. 2018;Decarli et al. 2019;Riechers et al. 2019).In each iteration, we extracted a number (P) from a random uniform distribution in the interval [0; 1].We then selected those sources entering in the realization having P ≤ F. In this way, sources with higher fidelity have a larger chance to be selected.We computed the number of sources and associated 1σ Poisson confidence intervals (Gehrels 1986) in Flux density (mJy) Fig. 6.Spectra of the CO(4-3) line-emitting sources in our ALMA survey of the MQN01 field.The red lines are the best fit models to the spectra (yellow bins).The green bins indicate the channels we used to compute the line-velocity integrated maps (see Fig. 7) defined within ±2σ of the best-fitting Gaussian line.The horizontal gray bands are the rms noise in each channel.The blue vertical line indicate the QSO CTS G18.01 systemic velocity.Sources are labeled according to the ID CO reported in Table 3. 0.5 dex bins.We then rescaled the resulting counts and uncertainties by the completeness corrections11 , and then divided them by the effective volume of the survey and by the luminosity bin width.Finally, we averaged over the realizations.We repeated the entire simulation five times by shifting the luminosity bins of 0.1 dex in order to mitigate the dependence of the result on the bin definitions and to expose the intra-bin variations.For bins with a less than one count on average, we report a 1σ upper limit.The resulting CO(4-3) LF is shown in Fig. 8 (panel a; red and gold symbols, see caption for a detailed explanation).In the 2.0" 0.0" -2.0" 2.0" 0.0" -2.0" 2.0"0.0"-2.0" 2.0"0.0"-2.0"2.0" 0.0" -2.0"  3. Sources with spectroscopic-redshift confirmation are marked with a violet square in the upper-left corner.
same figure, we also show the LF derived from blank fields at the same redshift reported from the literature.
Our aim is to compare the CO LF in MQN01 with that of blank fields.To do so, we combined data from ASPECS (The ALMA Spectroscopic Survey in the Hubble Ultra Deep Field Walter et al. 2016;González-López et al. 2019;Decarli et al. 2019Decarli et al. , 2020)), and Boogaard et al. (2023) who performed a blind search of molecular lines in the Hubble Deep Field North (HDF-N) using the Northern Extended Millimeter Array (NOEMA).These works provide us with the most up-to-date CO(4-3) LF at z ≃ 3.5.To enable a direct comparison with our result in the MQN01 field, we recomputed LFs consistently as described above adopting common luminosity bins.Additionally, we obtained the CO(4-3) LF from the SIDES simulation (Simulated Infrared Dusty Extragalactic Sky; Béthermin et al. 2022;Gkogkou et al. 2023) at z = 3 − 3.5.As a comparison, in Fig. 8, we also rescaled by a factor of 8× the best-fit Schechter function to the blank fields (see, Boogaard et al. 2023) as well as the CO LF predicted from SIDES.Interestingly, the CO LF in MQN01 differs significantly not only in normalization but also in shape with respect to that of blank fields showing a flattening at its bright end.We further discuss this point in Sect. 5.
We then computed the cumulative source number counts per comoving volume n CO (> L ′ CO ) by integrating the five different CO LFs separately, each of which has uncorrelated luminosity bins.We show our results in Fig. 8 (panel c).We also obtained cumulative source number counts from both the best-fit of blank fields and SIDES by integrating the corresponding CO LFs.In Fig. 8, we also rescaled the latter functions by a factor of 25× as a comparison to our results.
Finally, we evaluated the galaxy overdensity δ(> L ′ CO ) as traced by CO(4-3) by computing the bin-to-bin ratios between the cumulative number counts in MQNQ01 and those in blank fields.The measured cumulative overdensity of CO(4-3) line emitters is a strong function of luminosity, increasing from ∼ 8 (at the lowest luminosities) to ∼ 200 (at the highest luminosities).Above the limiting CO luminosity of our survey we estimated δ 4000 CO (> L ′ lim ) = 14 +9 −6 , and δ 1000 CO (> L ′ lim ) = 18 +14 −8 , in ∆V 4000 , and ∆V 1000 , respectively.The overdensity plot is also shown in Fig. 8 (panel d).4-3) LF in MQN01 and in blank fields.The (dark/light) red dots and boxes represent the CO LF in MQN01 field (fidelity (F) + completeness (C)-corrected/uncorrected, respectively) in ∆V 4000 , while gold hatched boxes are the uncorrected (raw) CO LF within ∆V 1000 .The blue and green boxes with black dots are data from ASPECS and HDF-N, respectively, which are representative of blank fields at these redshifts (e.g., Decarli et al. 2019Decarli et al. , 2020;;Boogaard et al. 2023).The downward arrows indicate 1-σ upper limits.The dashed black line is the best-fit Schechter function to the ASPECS+HDF-N data reported in Boogaard et al. (2023).The orange line is the CO(4-3) LF prediction from SIDES simulation (Béthermin et al. 2022;Gkogkou et al. 2023)

1.2-mm continuum source number counts
We investigate the galaxy overdensity in the MQN01 field as traced by dust continuum emission at 1.2 mm.Thanks to the combination of our deep spectroscopic VLT/MUSE and ALMA surveys we can pin-down the redshift of a sufficiently large sample of continuum-selected galaxies in the MQN01 field.We computed the source count density of our continuum-selected galaxies in a given cosmological volume around the CTS QSO G18.01 quasar by adapting the recipe from (e.g.,) Hatsukade et al. (2013Hatsukade et al. ( , 2016Hatsukade et al. ( , 2018)); Carniani et al. (2015); Aravena et al. (2016);Fujimoto et al. (2016Fujimoto et al. ( , 2023)); Umehata et al. (2017Umehata et al. ( , 2018)); González-López et al. (2019, 2020).We defined the differential source number count density per flux bin S 1.2 mm as or, In the computation of the source number density, we included all our 1.2 mm-continuum selected sources having spectroscopic redshift information either from the analysis of VLT/MUSE spectrum, or from their CO(4-3) line which lie within ∆ QSO < 1000 km s −1 (∆V = 430 cMpc 3 )12 .This selection yielded to a sample of six sources (including the QSO and the Object B; see Table 2, and 3).
To compare our results to blank fields, we employed the AS-PECS publicly-available catalogs13 (Aravena et al. 2019(Aravena et al. , 2020;;Boogaard et al. 2019Boogaard et al. , 2020;;González-López et al. 2019, 2020) of the ALMA band 6 survey from which we select all 1.2 mmcontinuum detected galaxies with spectroscopic redshift information.To increase the statistics of the sample we considered all sources with 2.5 ≤ z ≤ 3.5.This selection effectively yielded to a sample of ten galaxies in the redshift range z = 2.543 − 2.981 (median redshift z = 2.685).Conservatively, we employed such range to compute the corresponding cosmological volume.
Additionally, we compared our source number count density with that observed toward the SSA22 field which is found to show a large overdensity of Ly-α emitters and sub-millimeter Fig. 9.The overdensity of dusty star forming galaxies in the MQN01.Panel a): Differential 1.2 mm continuum source number count density in MQN01 (red points and boxes; including galaxies within ∆ QSO < 1000 km s −1 ); SSA22 (black dots and green boxes; Umehata et al. 2017Umehata et al. , 2019, within ±1000 km s −1 with respect to the median source redshift), and in blank fields as derived from ASPECS large program (black points and blue boxes; Aravena et al. 2019Aravena et al. , 2020;;Boogaard et al. 2019Boogaard et al. , 2020;;González-López et al. 2019, 2020)  galaxies (SMGs) at z ∼ 3 (see, Steidel et al. 1998Steidel et al. , 2000;;Yamada et al. 2012;Umehata et al. 2015Umehata et al. , 2017Umehata et al. , 2018Umehata et al. , 2019)).To this purpose, we cross-matched the catalog of the CO(3-2) emitters with that of continuum-selected galaxies at 1.1 mm as revealed by the ALMA deep field survey in the SSA22 field A (ADF22A; see, Umehata et al. 2017Umehata et al. , 2019)).The final sample comprises ten sources with redshift measurement within ±1000 km s −1 around the median redshift of the sources z = 3.0951.
The sources we selected from ASPECS surveys are observed at 1.2 mm but lie in a different redshift range, while that from Umehata et al. (2017Umehata et al. ( , 2019) are observed at 1.1 mm.In order to mitigate possible biases introduced in the selection of sources at different redshifts and which are observed at different wavelengths, we employed a "k-correction" to translate the observed monochromatic flux S λobs (z) of a source at redshift z to S 1.2 mm (z 0 ) defined as where ν 0 = c(1 + z 0 )/λ 1.2 mm and ν = c(1 + z)/λ obs .Here we estimated k λobs by assuming a modified black body emis-sion with typical dust temperature ranging in T dust = 20 − 45 K, and spectral index β = 1.5 − 2.0.We additionally take into account the contrast effect produced by the CMB at high-z (see, da Cunha et al. 2013) the temperature of which ranges between T CMB = 9.5 − 12.3 K at z = 2.5 − 3.5.Under such assumptions and by using a reference redshift of z 0 = 3.25, we found k 1.2 mm (z = 2.5) = 1.11 +0.06 −0.10 , k 1.2 mm (z = 3.5) = 0.97 +0.03 −0.02 , and k 1.1 mm (z = 3.1) = 0.796 ± 0.009, the uncertainties of which we took into account in converting the observed fluxes and computing the source number counts.
We therefore computed the differential and cumulative source number count densities, as well as source overdensity in MQN01, by adopting the same approach described in Sect.4.1.We used common bins for all the different fields corresponding to ∆ log S = 0.4, and we repeated the simulation two times by shifting the flux bins of 0.1 dex.Since all the sources involved in this computation are detected via multiple tracers, we set the fidelity of all sources to 1.In addition, we opted to ignore the completeness corrections which we expect to be not relevant due to the large uncertainties involved in such analysis.We report our results in Fig. 9. Overall, the 1.2-mm source number count density in MQN01 shows higher normalization and similar shape with respect that of blank fields at similar redshift without any evidence of a flattening at the bright end.Interestingly, our result resembles the source number count density measured in SSA22.Overall, the overdensity of sources detected in continuum at 1.2 mm in the MQN01 field within ∆ QSO < 1000 km s −1 is consistent with that estimated in the same redshift range by using CO emitters.We further discuss this point in Sect. 5.

The overdensity in MQN01
In Sect.4.1 we obtained the CO LF in the MQN01 field within ∆ QSO < 4000 km s −1 , and ∆ QSO < 1000 km s −1 .The analysis we performed on the data allowed us to estimate the presence of a galaxy overdensity traced by the molecular gas content.We found clear evidence of an overdensity in MQN01 which is a factor of ∼10 -100 higher than the field, depending on the luminosity cut.The luminosity dependence of the overdensity indicates that the most luminous objects are also the most overabundant.Taking into account all the CO line-detected sources (above the limiting CO(4-3) luminosity L ′ lim ), we estimated δ 4000 CO (> L ′ lim ) = 14 +9 −6 in ∆V 4000 , and counting only secure sources within ∆V 1000 without applying the completeness correction.Our results suggest that galaxies in MQN01 field are possibly part of a structure extending within a cosmological volume of at least ∼ 600 cMpc 3 as probed by our ALMA survey.
In addition, we can investigate the normalization and shape of the CO LF in MQN01 which encloses key information on the assembly of galaxies in such an overdense region.However, an accurate analysis of the LF trend and its interpretation is challenging given the low statistics of the sample and large uncertainties.To compare the CO LF in MQN01 with that measured in blank fields, we considered the one computed in ∆V 4000 corrected for both the source fidelity and the survey completeness.
The observed CO LF in the MQN01 field can be produced by a combination of an excess in galaxy number counts (thus increasing the overall normalization), and an enhanced CO luminosity of galaxies in the field (which translates to a shift of the source counts toward higher luminosities).However, a pure luminosity shift seems at odds with the data.Indeed, assuming that the CO LF in MQN01 follows the trend expected for blank fields (typically fit by a Schechter function; Schechter 1976), there is evidence for a further excess at the high-luminosity end producing a flattening at log L ′ CO /(K km s −1 pc 2 ) 10, which is also reflected in an increase of the cumulative overdensity values δ(> L ′ CO ) reported in Fig. 8. Interestingly, this excess is even more evident when considering the "raw" CO LF within ∆V 1000 .To demonstrate this fact, in Fig. 8 we report the rescaled version of the best-fit model to blank field data and the result obtained from the SIDES simulation.If the observed flattening at the high-luminosity tail of the LF is ascribed to a systematical increase of galaxy CO(4-3) luminosities, this can be either due to a higher galaxy molecular gas budget, or enhanced molecular gas excitation due to a significant star-formation activity and/or higher AGN fraction in the field which can increase the fraction of molecules populating higher-J CO states without necessary increasing the bulk of the CO-traced molecular gas mass.Interestingly, four out of seven (≃ 57%) of the CO-detected sources in ∆ QSO < 1000 km s −1 are also detected in the Xrays thus suggesting an intense AGN activity in the field (Travascio et al., in prep.).In the first scenario, if such galaxies follow the Kennicutt-Schmidt (KS) relation (Schmidt 1959;Kennicutt 1998;Kennicutt & Evans 2012, for a review), this would imply that galaxy assembly is accelerated in this dense field with more massive galaxies experiencing considerable episodes of star-formation.Regarding the second possible scenario, the sole detection of CO(4-3) is not sufficient to determine the dominant physical mechanism responsible for the molecular gas excitation.A proper sampling of multiple CO rotational ladders is therefore needed to understand if there is an important contribution to the enhancement of the CO(4-3) luminosity from high star-formation activity or X-ray radiation from AGN (see, e.g., Wolfire et al. 2022, for a review).Indeed, highly excited molecular gas and high star-formation efficiency have been found in the core of protoclusters at high-z (see, e.g., Lee et al. 2017;Coogan et al. 2018).However, since the aforementioned effects are interlaced, we are hampered in drawing strong conclusions here.Further information can be obtained by evaluating the SFR of the selected CO line-emitting galaxies, which is possible in the mm regime by measuring the galaxy FIR luminosities (see, e.g., De Looze et al. 2014).However, our observations provide us with at best two different photometric measurements in the Rayleigh-Jeans (RJ) tail of the dust SED, therefore preventing a proper estimation of the L FIR without an a-priori assumption on the galaxy dust temperature (T dust ).
The comparison of the CO LFs is based on the evaluation of differential number counts which are likely to be affected by Poissonian statistical oscillations.Cumulative luminosity function n CO (> L ′ ) can mitigate the uncertainties providing a more robust result.Here, the flattening that occurs at log L ′ CO /(K km s −1 pc 2 ) 10 is less sharp.However, there is a hint of a shallower decline of n CO with respect to blank fields.
We can then ask if the overdensity revealed in MQN01 via CO line emission is different from that revealed by other tracers.In Sect.4.2, we analyzed the 1.2-mm source number count density in MQN01 within ∆ QSO < 1000 km s −1 , and we compared with that of blank fields by selecting a sample of galaxies between 2.5 ≤ z ≤ 3.5 from ASPECS large program (see, Aravena et al. 2019Aravena et al. , 2020;;Boogaard et al. 2019Boogaard et al. , 2020;;González-López et al. 2019, 2020).This comparison yields to an overdensity signal which is also of the order of 10 -100 depending on the flux cut (see, Fig. 9).We stress that such analysis is affected by large uncertainties due to low statistics of the selected samples.In addition, the selection of galaxies which are detected in their rest-frame sub-mm dust continuum emission and for which spectroscopic redshift information is available from different lines can introduce biases which are difficult to assess (see, Sect.4.2 for further details).Any interpretation here should therefore be taken with caution.Given the above limitations, an analysis of the source number count trend as function of continuum flux density cannot be performed.Notwithstanding, we argue that overdensity traced by dust continuum is consistent with what we found from the analysis of CO LFs.The source number counts shown of MQN01 in Fig. 9 do not exhibit any evidence of a flattening as observed in the CO LF, and can therefore be interpreted as a pure normalization effect with respect to that of blank fields.However, the 1.2-mm flux density at z ∼ 3 (corresponding to ∼ 300 µm rest-frame wavelength) probes the RJ tail of the FIR dust emission which scales as ∼ M dust T dust (see, e.g., Hodge & da Cunha 2020).Therefore, any possible contribution to increasing the observed 1.2 mm flux density in galaxies located in the MQN01 volume cannot be directly attributed to either warmer dust temperature or a larger amount (M dust ) of cold dust in galaxies.We also compare the observed 1.2-mm number counts with that toward the SSA22 protocluster at z ∼ 3 (Umehata et al. 2017(Umehata et al. , 2019) ) obtained via a survey with similar design as we present here.The spectroscopic redshift information of galaxies in SSA22 are obtained through the measurement of CO(3-2) line thus providing a fair comparable sample.At zeroth order, the differential and cumulative number counts in MQN01 field are consistent with what has been found in SSA22 field.This evidence possibly suggests that similar physical mechanisms are caught in actions in these environments (we further discuss this point in Sect.5.3).

Molecular gas density in MQN01 field
We can estimate the total molecular gas content in the MQN01 structure by summing up the molecular gas masses of galaxies within the structure.To do so, we followed the recipe from (e.g.,) Decarli et al. (2016); Jin et al. ( 2021) where M H 2 is the galaxy molecular gas mass.To estimate M H 2 for each galaxy, we have to assume a CO(4-3)-to-CO(1-0) conversion factor (L ′ CO(1−0) = L ′ CO(4−3) /r 41 ), and an Bolatto et al. 2013, for a review).Such parameters are unknown for our CO(4-3) emitters in the MQN01 field.We adopted r 41 = 0.61 ± 0.13 obtained via modeling of the CO ladders of sources at z = 2.0 − 2.7 in the AS-PECS field (see, Boogaard et al. 2020, and references therein).The typical values of α CO adopted in the literature vary from 0.8 M ⊙ (K km s −1 pc 2 ) −1 for local ULIRGs (ultra luminous infrared galaxies), SMGs, and quasar hosts (Downes & Solomon 1998, but see, also, Montoya Arroyave et al. 2023 who measured a higher average value of 1.7 ± 0.5 M ⊙ (K km s −1 pc 2 ) −1 for local (U)LIRGs), to ∼ 4 M ⊙ (K km s −1 pc 2 ) −1 for the giant molecular clouds (GMCs) in the Milky Way (see,e.g., Bolatto et al. 2013;Carilli & Walter 2013, for further discussion).To take into account all the uncertainties, we computed ρ(H 2 ) by performing a Monte Carlo Simulation by varying all the measurements within their uncertainties.For the source fidelity and completeness, we adopted the same approach as in the evaluation of the CO LFs, and continuum number counts (see, Sects.4.1 and 4.2).During the simulation, we also accounted for the uncertainty on r 41 , and we varied α CO /(M ⊙ (K km s −1 pc 2 ) −1 ) uniformly in the range 0.8 − 4.3.This procedure yields ρ(H 2 ) = 3.7 +1.0 −0.9 × 10 8 M ⊙ cMpc −3 , where the nominal values and uncertainties are computed by taking the 50th, 5th, and 95th percentile, respectively.Similarly, we computed the molecular gas density for galaxies within ∆ QSO < 1000 km s −1 ignoring completeness corrections and we obtained ρ(H 2 ) = (1.0 ± 0.3) × 10 9 M ⊙ cMpc −3 .
In Fig. 10, we compare our results with the molecular gas density across cosmic time in blank fields as reported by several works in the literature.We report results from the CO Luminosity Density at High Redshift (COLDz; Pavesi et al. 2018;Riechers et al. 2019), ASPECS Large Program (Decarli et al. 2019(Decarli et al. , 2020)), the Very Large Array (VLA) -ASPECS (VLASPECS; Riechers et al. 2020), the IRAM (Institute for Radio Astronomy in the Millimeter Range) Plateau de Bure High-z Blue Sequence Survey 2 (PHIBBS2; Freundlich et al. 2019;Lenkić et al. 2020), the NOEMA survey in the HDF-N (Boogaard et al. 2023), and the xCOLD GASS survey (Saintonge et al. 2011(Saintonge et al. , 2017;;Fletcher et al. 2021).We also report measurements of the cosmic molecular gas density obtained via the dust continuum survey (see, Berta et al. 2013;Scoville et al. 2017;Magnelli et al. 2020), as well as the best-fit to data obtained from Walter et al. (2020).We found that the molecular gas budget residing in MQN01 is ∼ 10× the expected molecular gas density at z ∼ 3 consistent with what .Molecular gas density in the MQN01 field for galaxies within ∆ QSO < 1000 km s −1 (red star) compared to results in blank fields from HDF-N+ASPECS-LP (Boogaard et al. 2023), ASPECS-LP (Decarli et al. 2019(Decarli et al. , 2020)), PHIBBS2 (Lenkić et al. 2020), COLDz (Riechers et al. 2019), VLASPECS (Riechers et al. 2020), and the xCOLD GASS survey (Fletcher et al. 2021).We also report measurements of the cosmic molecular gas density obtained via dust continuum (gray bars; Berta et al. 2013;Scoville et al. 2017;Magnelli et al. 2020), and the best-fitting function from Walter et al. (2020).The green symbol is the molecular gas density in the Spiderweb protocluster (Jin et al. 2021).
we found from the analysis of the CO LF.This corroborates our findings indicating that galaxies in MQN01 are evolving in a gasrich overdense environment.
We also compare our result with that obtained in the Spiderweb protocluster at z = 2.16 reported by Jin et al. (2021, ρ(H 2 ) = 9.0 +3.6  −3.4 × 10 8 M ⊙ cMpc −3 ) as traced by CO(1-0).The emission from this line is a more direct proxy (because it does not require assumptions on the CO line excitation) of the molecular gas mass in galaxies (albeit it still requires an assumption on the α CO luminosity-to-H 2 mass conversion factor).Jin et al. (2021) adopted two typical values of α CO /(M ⊙ (K km s −1 pc 2 ) −1 ) for their sources equal to 0.8 for starburst-like objects (Emonts et al. 2018) and 3.6 for disk-like galaxies (Daddi et al. 2010).We note here that mm observations toward the Spiderweb protocluster probed a much larger volume (∼ 6600 cMpc 3 ) with respect to our ALMA survey in the MQN01 field .However, our measurements point to a molecular mass density which is comparable to that found in one of the most massive protoclusters at later times.

Comparison with (sub-)mm surveys of high-z protoclusters
Comparing our findings with results obtained in other fields hosting high-z galaxy-rich environments would allow us to put our results into a wider context.However, this is somehow difficult since a fair comparison requires at least a similar survey design at comparable redshift.Indeed, the observed galaxy population strongly depends on the observed tracers, selection method, sensitivity, and area of the survey.In our work we found a peculiar feature in the CO LF suggesting a flattening at the bright end relatively to the trend expected in blank fields.Despite numerous works in the literature reporting CO surveys in overdense environments (see, Sect.1), a systematic analysis of the CO LF in such environments across the redshift range is still missing.Jin et al. (2021) reported a large area (∼ 21 arcmin 2 ) survey of CO(1-0) with the Australia Telescope Compact Array (ATCA) toward the Spiderweb protocluster at z = 2.16.By using 46 identified CO emitters, they put a first constraint of the CO(1-0) LF in a protocluster environment.By fitting the data using a Schechter function with a fixed slope as determined by studies in blank fields (such as COLDz; see Riechers et al. 2019), they found a CO(1-0) LF normalization which is ∼ 1.5 dex higher than that in fields.Interestingly, the CO(1-0) LF reported by Jin et al. (2021), does not seems to show any signature of flattening at its bright end, even when considering a small volume (∼ 1650 cMpc 3 ) centered around the starbursting radio galaxy MRC1138-262 (see, e.g.Miley et al. 2006;Emonts et al. 2016) where the highest concentration of galaxies has been found.
Another attempt to study the LF in high-z protocluster field is reported by Hill et al. (2020).The latter present ALMA followup observations of SPT2349-56, a star-forming protocluster core at z = 4.3, targeting the CO(4-3) as well as the singly-ionized carbon fine-structure line [CII] 158 µm .They revealed 24 lineemitting sources in ∼ 7.2 arcmin 2 through which they obtained the FIR and [CII] 158 µm luminosity differential number counts.Intriguingly, they found that the FIR LF in the core of SPT2349-56 is biased toward bright galaxies compared to what is predicted for field galaxies, thus suggesting an enhanced star-formation activity possibly triggered by ongoing mergers (see, e.g., Tacconi et al. 2008;Engel et al. 2010;Luo et al. 2015) which are expected to be more common in overdense regions (see, e.g., Lotz et al. 2013;Hine et al. 2016).On the basis of the almost-linear CO-FIR luminosity relation (see, e.g.Liu et al. 2015;Kamenetzky et al. 2016) this would directly translate in an increased L ′ CO(4−3) luminosity at the bright end of the CO LF as we observe in the MQN01 field.
A comparable example to our observations is the SSA22 field, a protocluster at z ≃ 3.1 surveyed by Umehata et al. (2015Umehata et al. ( , 2017) ) with ALMA achieving a sensitivity of ∼ 60 µJy beam −1 over an area of ∼ 2 ′ × 3 ′ probing the 1.1 mm continuum.Therefore this survey is fairly comparable to our work both in terms of survey design and redshift.Umehata et al. (2018Umehata et al. ( , 2019) ) followed-up this field by extending both the probed area and including CO(3-2) line observations using ALMA band 3.In the combined survey, Umehata et al. (2018) report the detection of 35 SMGs and the 1.1-mm number counts which reveal an excess by a factor of several with respect to blank field.In this work we combined data from Umehata et al. (2018Umehata et al. ( , 2019)), and we rescaled the observed fluxes to 1.2 mm at z = 3.25 (see Sect. 4.2).The source number counts in MQN01 and SSA22 appear to be overlapped following a similar trend (albeit with large uncertainties).This might indicate that dust-obscured star-formation and metal production in galaxies in MQN01 field are enhanced in the protocluster.In addition, we conclude that the galaxy density as traced by the 1.2 mm dust continuum emission in MQN01 is comparable to that revealed in SSA22 field.

Galaxy velocity distribution
In Fig. 11 we report the distribution of ∆ QSO of the combined ALMA+MUSE sample, including sources with high-fidelity spectroscopic redshift within ∆ QSO < 4000 km s −1 .Interestingly, we can identify a "core" population of galaxies having ∆ QSO < 1000 km s −1 which appear to be Gaussian distributed Fig. 11.Right axis: velocity distribution of galaxies in MQN01 field either detected in their CO(4-3) lines in our ALMA band 6 survey (orange bins), or having high-fidelity spectroscopic redshift within ∆ QSO < 4000 km s −1 measured from MUSE data (violet bins).The blue hatched bins are the distribution of galaxy velocities in the "core".Left axis: the black and blue lines are the cumulative distribution of the whole sample, and that of the "core", respectively.The light blue curve is the theoretical prediction for a Gaussian distribution of velocities with central value and dispersion as reported in the legend.The gray band represents the 1-σ confidence interval.as expected for cluster (i.e., virialized) galaxies (see, e.g., Yahil & Vidal 1977), with the distribution median value appearing slightly shifted with respect to the QSO systemic redshift of ∼ −200 km s −1 .To test this hypothesis, in Fig. 11 (blue bins) we additionally report the (cumulative) velocity distribution of galaxies belonging to ∆ QSO < 1000 km s −1 , within the estimated core virial radius of R vir ∼ 100 kpc (see, Sect.5.4.2).We therefore tested the Gaussianity of the sample by performing a Lilliefors test (Lilliefors 1967) 14 and we obtained a p-value of 0.94.Therefore, we do not have significant evidence for rejecting the null hypothesis, but the small number statistics prevents here any strong conclusion.The mean velocity and the standard deviation of the "core" galaxies is 0 = −137 +84 −86 km s −1 , and σ = 208 +45 −59 km s −1 , where uncertainties are computed using bootstrap resampling.In Fig. 11 we report the theoretical cumulative Gaussian.The latter fairly reproduces the observed velocity distribution of the "core" galaxies.

Halo mass estimate
Under the assumption that the observed galaxies in the "core" of MQN01 field are virialized within a large dark matter (DM) halo (see, e.g., Miller et al. 2018;Hill et al. 2020), we can estimate the total halo mass on the basis of the observed velocity dispersion of the galaxies.To do so, we employed the scaling relation found by Evrard et al. (2008, see, also Rines et al. 2013) which is based on a suite of N-body simulations with various cosmologies: where h(z) = H(z)/100 km s −1 Mpc −1 , and σ r is the line-of-sight velocity dispersion of the halo members.
To define the halo members we proceeded iteratively.We took all galaxies within ∆ QSO < 1000 km s −1 , and we computed the dispersion of the line-of-sight velocity of galaxies (assuming they are due to purely peculiar motion), and the average position of the galaxies (i.e., the geometrical barycenter projected on the sky plane)15 .We therefore computed M 200 using Eq. 9, and the virial radius (R 200 ≡ R vir ) as R vir = [GM 200 /(100 H(z) 2 )] 1/3 .We then iterated the procedure by selecting all galaxies within R vir until convergence.By doing so, we determined a total dynamical halo mass of M 200 = 2.2 +4.0 −1.8 × 10 12 M ⊙ , and a virial radius of R vir = 94 +37 −41 kpc, where the nominal values and uncertainties are computed with a Monte Carlo sampling and by taking the 5th, 50th, and 95th percentile of the distribution.However, we note that such estimates involve a set of assumptions and therefore has to be taken with caution.In particular, we have assumed that the halo is spherical and that thus the velocity dispersion can be simply derived from the line-ofsight velocity.However, numerical simulations show that halos are likely to be triaxial spheroids whose axis ratios depends on halo mass and that the line-of-sight velocity, on average, underestimate the dispersion but has significant scatter (see, e.g., Elahi et al. 2018a,b).Also, the kinematical analysis of cosmological simulations presented by de Beer et al. ( 2023) suggest significant variations of the line-of-sight velocity of substructures within halos, due to asymmetries in matter accretion.Taken at face value, the mass estimate above would suggest that the host halo of the MQN01 quasar should be similar to the average quasar population at similar redshifts, (i.e., M 200 ≃ (1.5 ± 0.5) × 10 12 M ⊙ , see de Beer et al. 2023, and references therein).This is somewhat in contrast with the large overdensity of CO emitters found around MQN01 on scales of a few cMpc, implying that similar overdensities should have been commonly found around the typical z ∼ 3.5 quasars, unless the MQN01 is a structure on larger scales still in the process of merging.In order to test this hypothesis, we estimated the clustering signal in our field as discussed in detail in Appendix B. By fixing the slope of the cross-correlation function to γ = 1.8 as commonly done in the literature at these redshifts (see, e.g., Ouchi et al. 2004;Diener et al. 2017;Fossati et al. 2021;García-Vergara et al. 2022), we find a cross-correlation length for our CO emitters around the QSO CTS G18.01 of r 0,CO = 16 +13 −7 cMpc.However, we stress that the low number counts limits the statistical significance of our results in absolute terms.There are only very few other studies in the literature regarding the correlation length of CO emitters around typical quasars at z > 3.By looking at the fields of 17 quasars at z ∼ 3.8 over an area of ≃ 66 ′′ radius per field, García-Vergara et al. ( 2022) built a sample of five CO emitters in a velocity window of ±1000 km s −1 above the limiting luminosities (at 5.6σ) that are dependent on distance from the quasars and that typically vary from L ′ CO ≃ 4.5 × 10 9 K km s −1 pc 2 at the center of the field to ≃ 3.1 × 10 10 K km s −1 pc 2 at a distance of ∼ 1.9 cMpc (rescaled following the primary-beam response of ALMA observations).We notice that our ALMA band 3 mosaic CO(4-3) emitters MUSE-selected galaxies Fig. 12.Estimated three-dimensional galaxy velocities (including a statistical correction factor √ 3 assuming dynamical symmetry) as function of the projected separation with respect to the estimated halo center (also shown in Fig. 2).The CO(4-3)-line emitting galaxies, and those with high-fidelity spectroscopic redshift from MUSE within ∆ QSO < 4000 km s −1 are indicated by orange circles, and violet crosses, respectively.The arrows indicate sources that lie outside the reported velocity range.The vertical dashed black and green lines mark the estimated virial (R vir ) and "turnaround" radius (R ta ).The blue and red curves are the predicted escape velocity corresponding to a central point-like object and a truncated NFW mass profile, respectively with a total mass equal to the estimated M 200 .The shaded areas report the 1σ uncertainties.Galaxies which fall within the virial radius and the envelope are expected to be gravitationally bound.observation does not suffer from the same radial sensitivity dependence.In particular, García-Vergara et al. (2022) find a crosscorrelation length of r 0,CO = 12 +4 −3 cMpc (using the same fixed slope of γ = 1.8) and a total overdensity of δ = 17 +12 −8 .However, these estimates do not take into account a possible effect on the overdensity on the assumed luminosity cut.Our result shows that the overdensity increases with the CO emitter line luminosity (see, Fig. 8).Thus, we expect that a radial dependence of the sensitivity limits would artificially increase the overdensity as a function of radius and thus increasing the correlation length.We can directly test these expectations using our dataset.In particular, we have simulated in our data a similar sensitivity radial dependence as described above, obtaining a correlation length of r 0,CO = 23 +22 −10 cMpc which is significantly larger with respect to the previous estimate with uniform sensitivity and with respect to the value found around typical quasars at z ∼ 3.8 by García-Vergara et al. ( 2022) (albeit consistent within large uncertainties due to the low number statistics).This result suggests that either the MQN01 quasar host halo is much more massive than suggested by the kinematical analysis discussed above or that the quasar is surrounded by a rich structure consisting in multiple massive halos, or both.Ongoing surveys tracing the galaxy population on scales of tens of cMpc (Galbiati et al., in prep.) and aimed at characterizing the halos of other galaxies in the field will enable us to solve this dichotomy.
In the assumptions that the halo mass is consistent with the kinematical analysis above, we can also estimate which galaxies are bound to the core by comparing their line-of-sight velocity as a function of the projected distance from the estimated center of the halo (i.e., the geometrical barycenter of galaxies within the "core"; see, Fig. 12).In particular, we shifted the galaxy velocities by their mean value within the virial radius and rescaled them by √ 3. We also report the escape velocity assuming a central point-like mass esc = √ 2GM 200 /R 200 as well as a more realistic model which is a truncated NFW mass profile (Navarro et al. 1996)  16 .Galaxies that fall within the envelope are expected to be gravitationally bound.However, sources which are located well beyond the virial radius are expected to be dragged by the Hubble flow, unless they have significant tangential velocity component toward the halo center.More precisely, galaxies that are expected to escape from the DM halo potential are those located beyond the so-called "turnaround radius" which is estimated from numerical simulations to be R ta ≃ 3.3R vir (see, e.g., Gunn & Gott 1972;Busha et al. 2003Busha et al. , 2005)).This is the situation for a large fraction of galaxies, including 7 out of 11 COline emitters in our selected sample.Interestingly, we can identify a galaxy group (which includes sources L01, L02, L03, and L09) which lie within ∼ 200 − 500 kpc and |∆ | < 1000 km s −1 from the inferred halo center.Such symmetrical configuration cannot be due to selection effect and it could suggest the existence of a causal connection between the "core", and a possible larger structure extending a few hundreds of kpc.In particular, L01 -L03, and L09 show small line-of-sight velocity and are located close to the "turnaround" radius thus possibly suggesting that we are witnessing the moment in which these galaxies are decoupling from the Hubble flow and about to start their infall toward the core of the protocluster.However, we note that skyprojected distances of galaxies from the halo center has to be considered as lower limit to the true distances.Indeed, we expect to see in projection some galaxies having R sky ∼ R ta with threedimensional distance R > R ta .These galaxies are physically in the Hubble flow and dynamically disconnected from the protocluster.This possibly explains why some galaxies (in particular several MUSE-selected sources), have a projected distance within (or comparable to) the "turnaround" radius while exhibiting high line-of-sight velocities consistent with the Hubble flow.Finally, we note that galaxies having large distance from the halo center also tend to show large velocity shifts.This is the case of L05 -L08, which are all among the most distance (in the phase space) sources from the protocluster core and for which we do not have an independent redshift confirmation.Hence, in case future studies will determine that such sources are actually interlopers at different redshift or false-positive detections, this will not affect the dynamical analysis presented in this work.In conclusion, our analysis suggests that while the "core" appears to be bound, the entire system is not virialized.This conclusion is further supported by the measured galaxy overdensity considering that a virialized system should have a DM overdensity of around ∼ 200 which is always smaller than that traced by galaxies (see, e.g., Lacey & Cole 1994;Desjacques et al. 2018).
Finally, using the sources associated with the core, we can set a tentative limit of the baryon fraction in molecular form within the halo.We summed up the molecular mass of the individual CO-detected sources within the virial radius (i.e., QSO, Object B, and L04) which we estimated following the procedure described in Sect.5.2, assuming α CO = 0.8 M ⊙ (K km s −1 pc 2 ) −1 typical of quasar hosts and SMGs (e.g., Downes & Solomon 1998).By doing so, we found M core H2 = (8.8+0.9 −0.8 ) × 10 10 M ⊙ which 16 We truncated the NFW gravitational potential at r = R is about ∼ 25% of the total baryon mass of the halo assuming the maximum cosmic baryon fraction of Ω b /Ω m = 0.156 (Planck Collaboration et al. 2020).The molecular baryon fraction is therefore f mol ≡ M core H2 /M 200 < 0.04 +0.07 −0.02 (assuming the estimated DM halo mass above as a lower limit).

Summary and Conclusions
We presented an ALMA survey of the MQN01 field hosting a giant Ly-α nebula around a radio-quiet quasar at z = 3.25 (Borisova et al. 2016).Our survey is designed to primarily target the CO(4-3) transition as well as the 1.2-mm (rest-frame ∼ 300 µm) dust continuum of galaxies in the entire MUSE FoV (∼ 4 arcmin 2 ).The combination of our FIR observations with the multiwavelength information collected using multiple facilities, allowed us to unveil the molecular gas content and the cold dust emission of galaxies in this field.Below, we summarize our main findings: -We identified a robust sample of eleven CO(4-3) emitters within ∆ QSO < 4000 km s −1 , including a closely-separated companion (Object B) ∼ 1 ′′ south to the QSO CTS G18.01.Object B was not detected previously in optical data due to its proximity to the QSO.Nine sources (including the QSO) are in the area covered by the MUSE observations, five of which (≃ 56%, including Object B) do not have any MUSE counterparts thus highlighting the crucial role of (sub-)mm surveys to obtain a complete census of galaxy population at high-z.-Our observations revealed eleven sources through their 1.2mm dust continuum emission.Five of them (including the QSO, and Object B), are also detected in the CO(4-3) line, and another one having MUSE counterpart in the rest-frame UV.This implies the presence of six 1.2-mm-continuum selected galaxies in a narrow redshift range of ∆ QSO < 1000 km s −1 .-We analyzed the CO(4-3) LFs in MQN01 field within 2395 cMpc 3 ( ∆ QSO < 4000 km s −1 ), and 599 cMpc 3 ( ∆ QSO < 1000 km s −1 ) and we compare them with that of blank fields.We found a systematic excess of galaxy number counts which points to a galaxy overdensity of δ 4000 CO (> L ′ lim ) = 14 +9 −6 and δ 1000 CO (> L ′ lim ) = 18 +14 −8 above the limiting CO(4-3) luminosity of our ALMA survey.Notably, we found evidence of a flattening at the bright-end of the LF with respect to the expected trend in blank fields at similar redshift.Despite the limited statistics, this results suggest that massive galaxies in dense environments at z ∼ 3 are richer in molecular gas with respect to the field allowing them to grow faster than their counterparts in average environments.
-We obtained the number counts density in MQN01 for sources detected in their 1.2 mm dust continuum which have high-fidelity spectroscopic information from rest-frame UV metal absorptions from MUSE, or CO line.By building a comparison sample for blank fields exploiting ASPECS catalogs, and applying proper "k-corrections", we found evidence for a systematic excess of source number counts.The overdensity revealed by our ALMA 1.2-mm survey is consistent (within the large uncertainties) with that traced by the molecular gas.Our source number count density is also consistent with that in the SSA22 protocluster at z = 3.1 suggesting that obscured star-formation and consequently the metal production are similarly enhanced in such overdense environments.
Pensabene et al.: The ALMA view of MQN01 field -We studied the velocity distribution of galaxies in the field by using sources with high-fidelity spectroscopic redshift information from the MUSE and ALMA survey and we found evidence for a "virialized" structure which might represent the "core" of a protocluster in formation.Our dynamical analysis led to a total halo mass of 2.2 +4.0 −1.8 × 10 12 M ⊙ which, however, should be considered as a lower limit on the true halo mass due to geometrical reason.
Further efforts are needed to deeply understand how galaxies are shaped during their evolution in dense environments and to shed light on the galaxy-large-scale structure connection.In particular, in future works we will fully characterize the galaxy population in MQN01 by exploiting our rich multiwavelength dataset, and we will compare the physical properties of galaxies (such as, e.g., their stellar and dust mass, SFRs, gas fraction, morphology, AGN contribution) with that of field galaxies at similar redshift.This will enable us to better assess the environmental effect on the galaxy assembly in one of the densest region of the Universe discovered so far at z ∼ 3. Through our high-resolution ALMA observations we will also study the "core" of the MQN01 structure by characterizing in details the galaxy morphology and kinematics, in particular those of the QSO host and the closely-separated Object B. In future studies, we will also explore the correlation between galaxy properties and the diffuse ionized gas on large-scale as probed by the Ly-α line emission.All these works will help us to dissect the galaxy assembly processes to address the question about how galaxies get their gas and how this affects the galaxy gas content.

Fig. 1 .
Fig. 1.Combined primary beam response for our ALMA mosaics at the representative frequency of 109.0GHz (band 3, left panel) and for the aggregate continuum mosaic image of ALMA band 6 (right panel).The white contours correspond to a primary beam response of 30% (solid line) and 50% (dashed line).The circles show the disposition of the pointings with diameter equal to the Half Power Beam Width (HPBW) of the ALMA 12m antennas at the reference frequency of the setup.

Fig. 2 .
Fig.2.Footprints of our ALMA and MUSE observations toward the MQN01 field.The orange and blue contours encircle the area where the combined mosaic PB response is ≥ 0.5 for ALMA band 3 and 6, respectively.Within these areas we performed the source search.Violet contour draws the MUSE footprint.The background is a composition of the JWST NIRCam F322W2 image (center and south-est corner) complemented with the VLT/FORS2 R-band observation (north-est, north-and south-west corners, and the south-est gap between the JWST NIRCam detectors).The PSF-like emission of QSO CTS G18.01 has been subtracted in the JWST image revealing a nearby south-est quasar companion (Object B).Orange circles and blue squares indicate the locations of the ALMA CO(4-3) line-emitting and 1.2-mm continuum-selected sources in our ALMA band 3 and 6 observations, respectively.Green squares pin-point galaxy candidates detected in continuum at 3 mm in ALMA band 3. Grey hexagons denote sources corresponding to low-z counterparts, all but i5 show bright emission line in ALMA band 3. Violet crosses are the MUSE-selected LBGs belonging within ±1000 km s −1 with respect to the QSO CTS G18.01 systemic redshift.The dashed red circle shows the inner 100 pkpc around the estimated center of the protocluster core.

Fig. 3 .Fig. 4 .
Fig. 3. Completeness (red squares, left axis) and flux boosting (blue circles and density plot, right axis) corrected for the PB response as a function of 1.2-mm flux density (bottom axis) and S/N (top axis) of injected sources in the ALMA band 6 continuum image.The errorbars are derived by computing the 16th and 84th percentile of the distributions of completeness and flux boosting measurements in each flux bin.The red solid line is the best-fitting function to the completeness values modeled as C(F 1.2 mm ) = {1 + erf[(F 1.2 mm − A)/B]}/2 with A = 0.197 ± 0.004 and B = 0.083 ± 0.007.The density plot shows the values of the flux boosting as a function of the measured flux of the injected sources.

MQN01Fig. 5 .
Fig. 5.The observed frequency of CO and [CI] 609 µm transitions as a function of redshift.The boxes encircle both the frequency and redshift ranges covered by our ALMA band 3 observations.Boxes are colorcoded by the cosmological volume probed by each transitions.Grey bands show the SPW coverage in the lower (LSB) and the upper sideband (USB).The mean redshift of transitions entering in the USB are also reported.The blue arrow points to the box representing the volume within ±1000 km s −1 around the QSO CTS G18.01 systemic redshift.

Fig. 7 .
Fig. 7. Maps of the ALMA-selected galaxies in MQN01 field either detected in their CO(4-3) line or in continuum at 1.2-mm, or both.The color scale at the top refers to line-velocity integrated maps for all sources except C08, which is only detected at 1.2-mm in the dust continuum (bottom right color scale).Grey and black solid contours correspond to the velocity-integrated and 1.2-mm dust continuum maps, respectively and scales as [2, 3, 2 N ] × σ with N > 1 integer number.Dotted contours are the −2σ level.The yellow ellipse drawn in the bottom right corner represents the FWHM of the synthesized beam of the ALMA band 3 (yellow fill and gray line), and 6 (black line) observations.Sources are labeled according to the ID CO and/or ID 1.2 mm reported in Table3.Sources with spectroscopic-redshift confirmation are marked with a violet square in the upper-left corner.

Fig. 8 .
Fig.8.The overdensity of CO emitters in the MQN01 field.Panel a): CO(4-3) LF in MQN01 and in blank fields.The (dark/light) red dots and boxes represent the CO LF in MQN01 field (fidelity (F) + completeness (C)-corrected/uncorrected, respectively) in ∆V 4000 , while gold hatched boxes are the uncorrected (raw) CO LF within ∆V 1000 .The blue and green boxes with black dots are data from ASPECS and HDF-N, respectively, which are representative of blank fields at these redshifts (e.g.,Decarli et al. 2019Decarli et al. , 2020;;Boogaard et al. 2023).The downward arrows indicate 1-σ upper limits.The dashed black line is the best-fit Schechter function to the ASPECS+HDF-N data reported inBoogaard et al. (2023).The orange line is the CO(4-3) LF prediction from SIDES simulation(Béthermin et al. 2022;Gkogkou et al. 2023) at z = 3 − 3.25.Grey lines represent the rescaled versions of such models.Panel b): Number of sources in each luminosity bin averaged over the iterations.The horizontal dashed line reports the total average number of sources entering in each bin of the LFs.Panel c): Cumulative source number density as a function of the CO(4-3) line luminosity obtained by integrating the CO LFs.Panel d): Source overdensity in MQN01 field as a function of the CO(4-3) line luminosity obtained by a bin-to-bin ratio of the cumulative source number densities.
Fig.8.The overdensity of CO emitters in the MQN01 field.Panel a): CO(4-3) LF in MQN01 and in blank fields.The (dark/light) red dots and boxes represent the CO LF in MQN01 field (fidelity (F) + completeness (C)-corrected/uncorrected, respectively) in ∆V 4000 , while gold hatched boxes are the uncorrected (raw) CO LF within ∆V 1000 .The blue and green boxes with black dots are data from ASPECS and HDF-N, respectively, which are representative of blank fields at these redshifts (e.g.,Decarli et al. 2019Decarli et al. , 2020;;Boogaard et al. 2023).The downward arrows indicate 1-σ upper limits.The dashed black line is the best-fit Schechter function to the ASPECS+HDF-N data reported inBoogaard et al. (2023).The orange line is the CO(4-3) LF prediction from SIDES simulation(Béthermin et al. 2022;Gkogkou et al. 2023) at z = 3 − 3.25.Grey lines represent the rescaled versions of such models.Panel b): Number of sources in each luminosity bin averaged over the iterations.The horizontal dashed line reports the total average number of sources entering in each bin of the LFs.Panel c): Cumulative source number density as a function of the CO(4-3) line luminosity obtained by integrating the CO LFs.Panel d): Source overdensity in MQN01 field as a function of the CO(4-3) line luminosity obtained by a bin-to-bin ratio of the cumulative source number densities.
Fig.9.The overdensity of dusty star forming galaxies in the MQN01.Panel a): Differential 1.2 mm continuum source number count density in MQN01 (red points and boxes; including galaxies within ∆ QSO < 1000 km s −1 ); SSA22 (black dots and green boxes;Umehata et al. 2017Umehata et al. , 2019, within ±1000 km s −1 with respect to the median source redshift), and in blank fields as derived from ASPECS large program (black points and blue boxes;Aravena et al. 2019Aravena et al. , 2020;;Boogaard et al. 2019Boogaard et al. , 2020;;González-López et al. 2019, 2020)  including sources with spectroscopic redshift between 2.5 ≤ z ≤ 3.5.The downward arrows indicate 1-σ upper limits.Panel b): Number of sources in each luminosity bin averaged over the iterations.The horizontal gray line reports the total average number of sources entering in each bin for the computation of the differential number counts.Panel c): Cumulative source number density as function of 1.2 mm continuum flux obtained by integrating the differential number count density.Panel d): Cumulative source overdensity in MQN01 as a function of S 1.2 mm .The horizontal lines represent the overdensity range of CO(4-3) emitters in ∆V 1000 .

Table 1 .
Characteristics of ALMA observations toward the MQN01 field.

Table 4 .
Redshift bin, cosmological volume, and limiting luminosity of the main emission lines entering in the upper sideband of our ALMA band 3 survey.