# Tomography of the environment of the COSMOS/AzTEC-3 submillimeter galaxy at $z \sim 5.3$ revealed by $Ly\alpha$ and MUSE observations

Lucia Guaita, Manuel Aravena, Siddhartha Gurung-Lopez, Sebastiano

Cantalupo, Raffaele Marino, Dominik Riechers, Elisabete da Cunha, Jeff

Wagg, Hiddo S.B Algera, Helmut Dannerbauer, et al.

► To cite this version:

Lucia Guaita, Manuel Aravena, Siddhartha Gurung-Lopez, Sebastiano Cantalupo, Raffaele Marino, et al.. Tomography of the environment of the COSMOS/AzTEC-3 submillimeter galaxy at $z \sim 5.3$ revealed by $Ly\alpha$ and MUSE observations. Astronomy & Astrophysics - A&A, 2022, 660, pp.A137. 10.1051/0004-6361/202142384. [hal-04265439](https://hal.science/hal-04265439v1)

HAL Id: hal-04265439

https://hal.science/hal-04265439v1

Submitted on 30 Oct 2023

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Tomography of the environment of the COSMOS/AzTEC-3 submillimeter galaxy at z ~ 5.3 revealed by Ly$\alpha$ and MUSE observations*,**

L. Guaita<fnref n="1" />,<fnref n="2" />, M. Aravena<fnref n="2" />, S. Gurung-Lopez<fnref n="3" />,<fnref n="4" />, S. Cantalupo<fnref n="5" />,<fnref n="6" />, R. Marino<fnref n="6" />, D. Riechers<fnref n="7" />, E. da Cunha<fnref n="8" />, J. Wagg<fnref n="9" />, H. S. B. Algera<fnref n="10" />,<fnref n="11" />, H. Dannerbauer<fnref n="12" />,<fnref n="13" />, and P. Cox<fnref n="14" />

Received 7 October 2021 / Accepted 4 February 2022

## ABSTRACT

**Context.** Submillimeter galaxies (SMGs) have been proposed as the progenitors of massive ellipticals in the local Universe. Mapping the neutral gas distribution and investigating the gas accretion toward the SMGs at high redshift can provide information on the way SMG environments can evolve into clusters at $z = 0$.

**Aims.** In this work, we study the members of the protocluster around AzTEC-3, a submillimeter galaxy at $z = 5.3$. We use Ly$\alpha$ emission and its synergy with previous CO and [CII]158 µm observations.

**Methods.** We analyzed the data from the Multi Unit Spectroscopic Explorer (MUSE) instrument in an area of $1.4 \times 1.4 \text{ arcmin}^2$ around AzTEC-3 and derived information on the Ly$\alpha$ line in emission. We compared the Ly$\alpha$ profile of various regions of the environment with the zELDA radiative transfer model, revealing the neutral gas distribution and kinematics.

**Results.** We identified ten Ly$\alpha$ emitting sources, including two regions with extended emission: one embedding AzTEC-3 and LBG-3, which is a star-forming galaxy located 2'' (12 kpc) north of the SMG and another toward LBG-1, which is a star-forming galaxy located 15'' (90 kpc) to the southeast. The two regions extend for $\sim 27 \times 38 \text{ kpc}^2$ ($\sim 170 \times 240 \text{ ckpc}^2$) and $\sim 20 \times 20 \text{ kpc}^2$ ($\sim 125 \times 125 \text{ ckpc}^2$), respectively. The sources appear distributed in an elongated configuration of about 70'' (430 kpc) in extent. The number of sources confirms the overdensity around AzTEC-3. We study the MUSE spectra of the AzTEC-3+LBG-3 system and LBG-1 in detail. For the AzTEC-3+LBG-3 system, the Ly$\alpha$ emission appears redshifted and more spatially extended than the [CII] line emission. Similarly, the Ly$\alpha$ line spectrum is broader in velocity than [CII] for LBG-1. In the former spectrum, the Ly$\alpha$ emission is elongated to the north of LBG-3 and to the south of AzTEC-3, where a faint Ly$\alpha$ emitting galaxy is also located. The elongated structures could resemble tidal features due to the interaction of the two galaxies with AzTEC-3. Also, we find a bridge of gas, revealed by the Ly$\alpha$ emission between AzTEC-3 and LBG-3. The Ly$\alpha$ emission toward LBG-1 embeds its three components. The HI kinematics support the idea of a merger of the three components.

**Conclusions.** Given the availability of CO and [CII] observations from previous campaigns, and the Ly$\alpha$ information from our MUSE dataset, we find evidence of starburst-driven phenomena and interactions around AzTEC-3. The stellar mass of the galaxies of the overdensity and the Ly$\alpha$ luminosity of the HI nebula associated with AzTEC-3 imply a dark matter halo of $\sim 10^{12} \text{ M}_\odot$ at $z = 5.3$. By comparing this with semi-analytical models, the dark matter halo mass indicates that the region could evolve into a cluster of $2 \times 10^{13} \text{ M}_\odot$ by $z = 2$ and into a Fornax-type cluster at $z = 0$ with a typical mass of $2 \times 10^{14} \text{ M}_\odot$.

**Key words.** galaxies: high-redshift – galaxies: interactions – galaxies: evolution – galaxies: kinematics and dynamics – galaxies: starburst – galaxies: clusters: general

<footnote n="1">¹ Departamento de Ciencias Fisicas, Facultad de Ciencias Exactas, Universidad Andres Bello, Fernandez Concha 700, Las Condes, Santiago, Chile
e-mail: lucia.guaita@unab.cl</footnote>

<footnote n="2">² Núcleo de Astronomía, Facultad de Ingeniería y Ciencia, Universidad Diego Portales, Av. Ejército 441, Santiago, Chile</footnote>

<footnote n="3">³ Observatori Astronòmic, Universitat de València, C/ Catedrático José Beltran, 2, 46980 Paterna, València, Spain</footnote>

<footnote n="4">⁴ Departament d'Astronomia i Astrofísica, Universitat de València, 46100 Burjassot, València, Spain</footnote>

<footnote n="5">⁵ Dipartimento di Fisica, Università di Milano Bicocca, Piazza della Scienza 3, 20126 Milano, Italy</footnote>

<footnote n="6">⁶ Institute for Astronomy, ETH Zurich, 8093 Zurich, Switzerland</footnote>

<footnote n="7">⁷ I. Physikalisches Institut, Universität zu Köln, Zülpicher Strasse 77, 50937 Köln, Germany</footnote>

<footnote n="8">⁸ International Centre for Radio Astronomy Research, University of Western Australia, 35 Stirling Hwy, Crawley, WA 6009, Australia</footnote>

<footnote n="9">⁹ SKA Organization, Lower Withington Macclesfield, Cheshire SK11 9DL, UK</footnote>

<footnote n="10">¹⁰ Hiroshima Astrophysical Science Center, Hiroshima University, 1-3-1 Kagamiyama, Higashi-Hiroshima, Hiroshima 739-8526, Japan</footnote>

<footnote n="11">¹¹ National Astronomical Observatory of Japan, 2-21-1, Osawa, Mitaka, Tokyo, Japan</footnote>

<footnote n="12">¹² Instituto de Astrofísica de Canarias (IAC), 38205 La Laguna, Tenerife, Spain</footnote>

<footnote n="13">¹³ Universidad de La Laguna, Dpto. Astrofísica, 38206 La Laguna, Tenerife, Spain</footnote>

<footnote n="14">¹⁴ Sorbonne Université, UPMC Université Paris 6 and CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98b boulevard Arago, 75014 Paris, France</footnote>

<footnote>* The reduced mosaic of the MUSE observations is only available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsarc.u-strasbg.fr/viz-bin/cat/J/A+A/660/A137</footnote>

<footnote>** Based on observations made with ESO Telescopes at the Paranal Observatory, under programme ID 094.A-0487.</footnote>
# 1. Introduction

In the hierarchical theory of structure formation, initial small density fluctuations give rise to the formation of the first stars and galaxies. These structures subsequently grow larger and more massive via mergers and accretion (e.g., White & Rees 1978). This produces denser and denser regions with a strong gravitational potential that can influence the distribution and the evolution of galaxies.

The densest regions at high redshift can be pinpointed by searching for massive dusty galaxies, such as submillimeter galaxies (SMGs). Submillimeter galaxies (Blain et al. 2002; Casey et al. 2014) are characterized by strong star-formation events, typically associated with large amounts of gas and dust which often obscures rest-frame ultraviolet to optical wavelengths (Dudzevičiūtė et al. 2020). They typically present rapid gas consumption through high star-formation efficiencies that are associated with major mergers. Overdense regions around SMGs are interesting because SMGs are thought to be the progenitors of the most massive galaxies in the local Universe (Chapman et al. 2005; Swinbank et al. 2008; Stach et al. 2021).

By studying the distribution of the neutral gas around SMGs and within their overdense environment, we can investigate accretion events and obtain information on their evolution. An efficient method to reveal the distribution and study the kinematics of the neutral hydrogen (HI) gas at high redshift is the detection of Lyα emission (e.g., Cantalupo et al. 2014; Matthee et al. 2020a; Daddi et al. 2021).

It is important to note that Lyα is the strongest recombination line of neutral Hydrogen. It is mainly produced in star-forming regions where the ionizing radiation emitted by young, hot O and B stars ionize their surrounding gas which recombines in relatively short timescales, depending on its density (>1 atoms cm⁻³) and temperature (10⁴ K) conditions (Dijkstra 2016; Cantalupo 2017). Due to their short wavelength, Lyα photons are easily absorbed by dust and, due to the resonant nature of the Lyα transition, they are scattered by HI atoms. Also, the escape of Lyα photons out of the dense interstellar medium (ISM) gas can be favored by HI kinematics, which can produce an asymmetric profile where the red peak is the dominant one in case of an outflowing gas (Verhamme et al. 2006; Schaerer & Verhamme 2008; Laursen et al. 2009; Orsi et al. 2008; Gurung-López et al. 2019a). These phenomena produce an emission line with a red and a blue peak and some absorption at the resonant wavelength of Lyα. The relative intensity and the separation between the two peaks is related to the amount and distribution of dust and HI atoms, and the Lyα emission is typically spatially more extended than the UV continuum (e.g., Steidel et al. 2011; Wisotzki et al. 2016; Leclercq et al. 2017). Another process that can produce Lyα photons and that has been proposed to explain the very extended Lyα nebulae observed at high redshift is collisional excitation (Haiman et al. 2000; Fardal et al. 2001) with electrons. This process converts the thermal energy of the gas into radiation, and therefore cools the gas. It is typically less energetic than recombination, except in the case of a strong gravitational potential.

Together with Lyα, also [CII]158 µm has been detected in galaxies at $z > 4$ and used to investigate the distribution and the properties of their gas (e.g., Capak et al. 2015; Maiolino et al. 2015; Pentericci et al. 2016; Carniani et al. 2018). [CII]158 µm is the dominant cooling line of the ISM in star-forming galaxies, where it can carry up to 1% of the far-infrared luminosity (e.g., Israel et al. 1996). It can originate from different phases of the ISM (Hollenbach & Tielens 1999; Stacey et al. 1991; Goldsmith et al. 2012; Pineda et al. 2013); mainly in the photo-dissociation regions where the UV radiation from hot stars can dissociate molecules and ionize atoms (Stacey et al. 1991), also in the cold neutral atomic medium, but sometimes either in the warm or in the ionized medium (e.g., Pineda et al. 2013; Maiolino et al. 2015). The association with the molecular CO lines can help in understanding its origin from molecular clouds. Therefore, [CII]158 µm is an ideal tracer of star formation, but also of the distribution, dynamics, and enrichment of the ISM in star-forming galaxies and also of their circum-galactic medium (CGM).

We expect that the HI gas in dense environments affect the stage of evolution of the galaxies inside the dense regions and the shape of the emission of Lyα and possibly [CII]158 µm. For instance, the observations of protoclusters around radio galaxies had shown an inside-out picture in which dusty starburst galaxies are located in the cores and young Lyα emitters are distributed in the outskirts (Kuiper et al. 2011). However, this picture is not so straightforward because there are also dense regions at high redshift composed of more than one dense clump (e.g., Cortese et al. 2006; Zheng et al. 2016; Guaita et al. 2020) and also dense clumps could be traced by different galaxy populations. Shi et al. (2019), for example, showed an overdensity of Lyα emitters separated from a peak traced by more massive Lyman Break galaxies, maybe indicating the presence of a stream of gas falling into one side of the structure. The SSA22 protocluster (Steidel et al. 2000) at $z \sim 3$ is an example in which the overdensity is traced by both Lyα emitters and Lyman Break galaxies. The filamentary structure traced by the Lyα emitting sources shows the effect of the gravitational potential that seems to result in the formation of giant Lyα nebulae in the intersection of the filaments.

Furthermore, in dense regions, the major, gas-rich merger rate can be higher than in the field, given the possible galaxy encounters (e.g., Hine et al. 2016; Tacconi et al. 2013). In a merging process, the gas can be compressed either in the central and in the external regions of the progenitor galaxies, favoring strong episodes of star formation (e.g., Hernquist 1989). The properties of the merging gas and of the star-formation phenomenon could shape the Lyα emission of the merging system. Yajima et al. (2013) studied the formation of extended Lyα emission from interacting galaxies at high redshift using a combination of hydrodynamic simulations with three-dimensional radiative transfer calculations. They showed that at the first passage the triggered star-formation event produces two Lyα peaks coinciding with the two progenitor nuclei. The Lyα peaks decrease in intensity after the first passage and increase again in the second passage. The Lyα emission continues to decline and its distribution becomes more compact in the final stage of the merger. Recently, Romano et al. (2021) studied 75 main-sequence star-forming galaxies at $4.4 < z < 5.9$ through their [CII]158 µm emission and reported about the presence of significant merger activity at $z \sim 5$.

In this paper, we study the environment surrounding the submillimeter galaxy, AzTEC-3, at $z = 5.3$ (Capak et al. 2011). AzTEC-3 was discovered (Younger et al. 2007) in the AzTEC survey (Austermann et al. 2009) and it was the first SMG found at $z > 5$ with a large overdensity associated with it. Capak et al. (2011) discovered at least ten star-forming galaxies in a circular area of 2 comoving Mpc (cMpc) around AzTEC-3 and estimated that the SMG makes at least one-fourth of the mass of the entire system. Riechers et al. (2010) detected CO molecular gas emission from a compact ($<2.3$ kpc scale) region on top of the SMG, estimating a CO mass on the order of $5 \times 10^{10} M_{\odot}$ that
could be depleted in 50 Myr at a star-formation rate (SFR) of 1100 $M_{\odot}$ yr$^{-1}$. The detection of luminous CO emission implies enrichment with heavy elements in the material that fuels the star formation in AzTEC-3. The compactness of the CO-emitting region and the consequent high star-formation rate density make AzTEC-3 a special object to study, probably affected by its special environment (see also Riechers et al. 2020). In fact, a Lyman Break galaxy (LBG-3, nomenclature from Riechers et al. 2014) is observed at 2" (about 75 comoving kpc = ckpc) and another one (LBG-1, nomenclature from Riechers et al. 2014) is observed at 15" (about 580 comoving kpc) from AzTEC-3. LBG-1 is composed of three main clumps as seen in the rest-frame UV images from the Hubble Space Telescope (HST) and it makes a significant fraction of the mass of the protocluster (about half of the mass of AzTEC-3). From the spectral energy distribution point of view, LBG-1 can be considered as a typical star-forming galaxy at $z \sim 5$, with little dust obscuration and a young starburst in addition to an old underlying stellar population (Capak et al. 2011). No CO emission was detected in LBG-1 (Riechers et al. 2014; Pavesi et al. 2019) and this is consistent with its low metallicity. Riechers et al. (2014) used the Atacama Large Millimeter/submillimeter Array (ALMA) and detected [CII]158 µm emission from the same compact region of AzTEC-3 as CO. Also, they observed a hint of a tidal feature toward LBG-3. [CII]158 µm emission was also detected at the position of LBG-1 with a hint of velocity gradient, which could be related to a merging event within the three clumps.

Keeping the previous observations in mind, we study here deep VLT/MUSE (Multi Unit Spectroscopic Explorer) observations of the region around AzTEC-3 and LBG-1. The MUSE data cover the 800–1500 Å rest-frame wavelengths at $z \approx 5.3$, including $Ly\alpha$. We take advantage of the synergy between the previous [CII] information and the MUSE $Ly\alpha$ detections to study the properties of the galaxies in the protocluster and to investigate whether the environment could play a role in the stage of evolution of AzTEC-3. It is worth noting that many $z > 5$ SMGs are not detected at UV wavelengths in deep HST data (e.g., HDF850.1 and GN10 as explained in Riechers et al. 2020), so AzTEC-3 is a relatively rare case where MUSE UV studies are possible.

The paper is organized as follows. In Sects. 2 and 2.1, we explain the MUSE observation and the reduction of the data. In Sect. 2.2, we describe the method adopted to analyze the MUSE datacube and to detect $Ly\alpha$ emitting sources. In Sect. 3, we present the sample of $Ly\alpha$ emitting galaxies and star-forming galaxies in the field. In Sects. 4.1 and 4.2, we discuss in detail the spectroscopic properties of the $Ly\alpha$ emission around AzTEC-3 and LBG-1, the two galaxies with known systemic redshift from previous [CII] observations. In Sect. 5, we present the radiative transfer model used to fit the $Ly\alpha$ emission profiles. In Sect. 6, we provide a discussion about the AzTEC-3 environment coming from our MUSE $Ly\alpha$, the CO and [CII] observations. In Sect. 7, we summarize our work. Throughout the paper, we adopt a standard cosmology ($H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$, $\Omega_0 = 0.3$).

## 2. Observations

The field around the AzTEC-3 (RA = 150.0863, Dec = 2.589) submillimeter galaxy was selected to match the central region of the CO Luminosity Density at High-$z$ (COLDz) survey (Pavesi et al. 2018a; Riechers et al. 2019), which corresponds to the deepest field with CO(2–1) coverage at $z \approx 5.3$. It also matches the COSMOS-XS survey field (van der Vlugt et al.

2021; Algera et al. 2020) which has the deepest radio continuum data in the COSMOS field (see also Algera et al. 2021).

Data from the HST are available in the F606W, F814W, F105W, F125W, and F160W filters. The F814W image was obtained in Cycle 12 and 13 (Scoville et al. 2007), the F606W, F125W, and F160W images were obtained in 2014 (Riechers 2013), and the F105W image was obtained in Cycle 22 (PI: P.L. Capak) (see also Barisic et al. 2017). The environment around AzTEC-3 is quite rich. As we can see in Fig. 1, within a radius of 2" there are at least three sources at similar redshift. At 15" from AzTEC-3, there is also LBG-1 (Riechers et al. 2014).

During December 2014 and February 2015, we observed the AzTEC-3 field with the Multi Unit Spectroscopic Explorer (MUSE, Bacon et al. 2010, 2014, 2015). MUSE is a second generation instrument installed on the Nasmyth focus of UT4 at the Very Large Telescope (Chile). The spectral resolution of MUSE at 7000 Å is about 2700 and each resolution element is sampled by 2.5 pixels along the spectral direction. The spectral sampling is 1.25 Å per pixel.

The pixel scale of MUSE is 0.2" per pixel (about 8 ckpc per pixel) and the point spread function (PSF) of our MUSE data is 0.7", corresponding to about 4 kpc or 27 ckpc at $z = 5.3$. The field was observed in four MUSE pointings, two of which contain AzTEC-3 at the edges, and the total area is $1.4 \times 1.4 \text{ arcmin}^2$.

Each pointing was observed with 20 frames of 900 s, making a total exposure time of 20 h in the very center of the overlapping regions. For every observation night, we were able to use bias, dark, dome flat, sky flat frames from our run to perform the basic reduction steps, wavelength-calibration files, standard stars, and illumination-correction files close to each observing night were also used for proper calibration.

### 2.1. Reduction of the MUSE datacubes

We used the standard ESO MUSE pipeline (Weilbacher et al. 2014, v2.2) to perform the initial reduction of the MUSE datacubes. The pipeline carries out bias, dark and flat field corrections, calibrates the data in wavelength and astrometry, and applies a basic illumination correction. Afterwards, we improved the quality of the illumination correction and sky subtraction by applying the CubExtractor (CubEx) software following Cantalupo et al. (2019) (see also Borisova et al. 2016; Fumagalli et al. 2016, 2017; Mackenzie et al. 2019, for details).

We reduced the MUSE dataset in two ways. First, we reduced the datacubes corresponding to each pointing separately. In the final combinations, we generated a mean and median stack of the four pointings. Then, we produced a mosaic by combining the information contained in all the 80 MUSE OBJECT_PIXTABLEs and by generating individual cubes, all of the size of the final mosaic. Basically, NaN values are inserted in the areas of the individual cubes outside the observation field of view. We calculated the frame offset before combination, by cross-matching the coordinates of the sources identified with the SourceExtractor software (Bertin & Arnouts 1996). We made a mean combination of all the frames ("mosaic"), two mean combinations of the even ("even mosaic") and odd ("odd mosaic") numbered subexposures to have independent sets of data to be used as sanity checks of the detections, and a median combination of all the frames.

We performed the two reductions because the mosaic is more difficult to handle in the stage of the analysis due to its size, but it is fundamental to improve the signal-to-noise ration (S/N) of the overlapping areas among the pointings, in particular the region
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Left panel: mosaic of the entire area observed with MUSE, obtained collapsing the cube information in the 7600–7700 Å wavelength range (“mosaic”), which corresponds to the redshifted $Ly\alpha$ at $z \approx 5.3$. We identified six regions characterized by similar background level, the center of the four pointings (P1, P2, P3, P4), the central vertical and horizontal (black rectangles) stripes of overlapping regions. The white square corresponds to the area zoomed in the right panel. Right panel: zoom of the region containing AzTEC-3 (upper red circle) and LBG-1 (lower red circle) in the HST F160W image. We show the position of LBG-3 as black square, the position of LBG-2 (1447523 in Capak et al. 2011 and COSMOS2015_849887 in Laigle et al. 2016) as blue circle, and the position of COSMOS2015_848724 as black circle. These sources are all at $z \sim 5.3$. The yellow circle indicates the position of a knot of star formation which can also be part of AzTEC-3. It is worth noting that LBG-1 is located in the center of the mosaic, in the lower corner of the 20 h overlapping region.</figcaption></figure>

just to the south of AzTEC-3. The difficulty comes from the fact that the background level of the mosaic is different at the center of the pointings and in the overlapping areas and this could affect the measured S/N in the source detection. To overcome this difficulty we divided the mosaic into six regions, the center of each of the four pointings, the vertical overlapping area (vertical stripe), and the horizontal overlapping area (horizontal stripe), as shown in Fig. 1. However, these regions have sharp separations and the reduction of the four pointings separately are used to identify sources otherwise missed due to border effects. As expected, the detection limit increases in the overlapping regions (Table 1). For comparison, the $5\sigma$ detection limit of the white light image of the entire area is 26.2 when the central wavelength of the MUSE coverage is taken as the reference wavelength. We matched the astrometry of MUSE to that of the images of the Hubble Space Telescope (HST) and we reached an accuracy better than 0.1" on average. In Fig. 1, we show the narrow bands obtained collapsing the MUSE datacube in the 7600–7700 Å wavelength range, that contains the observed wavelength of the $Ly\alpha$ emission line at $z = 5.3$ ($\lambda_{\text{air}} \sim 7654$ Å).

## 2.2. Detection and identification of the line emitters

The inspection of the combined MUSE datacubes was performed with the CubEx software (Cantalupo et al. 2019). To maximize the S/N of the detections, we first ran a continuum-subtraction procedure, which is performed by the CubeBKG-Sub function of the CubEx package. The function subtracts the continuum sources in the field, likely to be low-$z$ objects. The continuum subtraction is performed through a median filtering along the spectral dimension, spaxel by spaxel. Following Marino et al. (2018), we chose the size of one cell of continuum in the wavelength direction equal to 50 and the continuum filter radius in the wavelength direction equal to 3.

<table><caption>**Table 1.** $5\sigma$ detection limits of the mean-combination narrow bands.</caption><thead><tr><td>NB</td><td>mag_7660</td></tr></thead><tbody><tr><td>P1</td><td>26.2</td></tr><tr><td>P2</td><td>26.0</td></tr><tr><td>P3</td><td>26.1</td></tr><tr><td>P4</td><td>26.1</td></tr><tr><td>Horizontal stripe</td><td>26.7</td></tr><tr><td>Vertical stripe</td><td>26.7</td></tr><tr><td>Entire mosaic</td><td>26.3</td></tr></tbody></table>

**Notes.** The narrow-band (NB) images correspond to the ones showed in Fig. 1, the center of the four pointings, the vertical and horizontal stripes of overlapping regions, and the entire mosaic. The $5\sigma$ detection limits in the second column are the magnitudes corresponding to a flux density at 7660 Å ($mag\_7660$) within an aperture of 1". The Point Spread Function of the entire mosaic is 0.7".

To avoid the subtraction of emission lines in the wavelength range of our interest, some wavelengths were masked in the CubeBKGSub procedure. For the specific case of the detections of the members of the AzTEC-3 protocluster, we masked the wavelength range corresponding to the $Ly\alpha$ emission line at the redshift of the AzTEC-3 submillimeter galaxy (2324–2336 pixels in the spectral direction, corresponding to $[-40, +550]$ km s$^{-1}$ from the AzTEC-3 systemic redshift at the $Ly\alpha$ wavelength).

Then, CubEx was used to perform extraction, detection, and photometry of sources with arbitrary spatial and spectral shapes directly in the datacubes. We rescaled the data variance (see Mackenzie et al. 2019) and we applied smoothing before extraction. The smoothing function is a Gaussian with a two-pixel radius in the spatial direction. The datacube elements, called
'voxels', were selected as detections if they matched the following criteria (see also Marino et al. 2018; Mackenzie et al. 2019; Cantalupo et al. 2019): signal-to-noise ratio of the individual voxels larger than 3; integrated S/N ratio in the 3D space equal to 7; integrated spectral S/N (using 1 ds optimally-extracted spectra) equal to 3; S/N of an individual voxel equal to 3 if the voxel is connected to a previously detected voxel; minimum number of connected voxels for detection equal to 30; minimum number of spectral pixels for detection equal to 3. The main sources of fake detections were sky line residuals.

We ran CubEx in the 1000–1100 pixel wavelength range (arbitrary much lower wavelength than 7660 Å) to immediately remove bright stars from the catalog. With the CubEx package we also generated optimally extracted images (OEimages) of the detections showing the extension of the detection emission line in the spatial direction (details in Borisova et al. 2016), weighted by signal to noise, and the signal-to-noise maps. We chose OEimages that are smoothed with a box car of a 2-pixel radius in the spatial direction for inspections.

To confirm a line detection at $z \sim 5.3$, we therefore, made sure that it did not show any counterpart in the HST F606W image; it had reasonable segmentation and S/N datacubes; it was not at the wavelength of a strong sky line residual; it had a redshift not more than $2000 \text{ km s}^{-1}$ from the AzTEC-3 systemic redshift; it was not associated to a detection at a different wavelength implying the emission line was not $Ly\alpha$; it had a reasonable S/N also in the even and odd mosaics (integrated S/N ratio larger than 5). This was also useful to exclude detections from cosmic rays (e.g., Lofthouse et al. 2020). Finally, we checked that its 1 d spectrum showed a reasonable $Ly\alpha$ shape at the detection wavelength, for example not too narrow like a spike of a few bad pixels. The criteria were applied in the form of a visual inspection of the OEimages and of the spectra.

To estimate the fidelity of our detections (i.e., down to which signal to noise the line detections start to be spurious), we applied the same set of criteria for the detection datacubes and for the datacubes multiplied by -1, that is the negative cubes. In Fig. 2, we show the ratio between the number of negative and positive detections as a function of signal to noise for the entire mosaic and for the four individual pointings. The signal to noise is calculated as the ratio between the isophotal flux integrated over the line detection by CubEx and its error. This signal-to-noise definition is different from the one based on pixel-to-pixel noise and that can produce lower values, but directly uses the output of CubEx. A S/N equal to 11 allowed that more than 60% of the detections were likely to be real detections in the mosaic. Therefore, we fixed to 11 the S/N of the line detections that made our “Main” sample. We propagated the 40% chance to have a fake detection in the estimation of the density of the protocluster (see the following sections), even if the visual inspection described above already removed 20% of the sources in the fidelity-cut catalog.

The choice of the detection parameters was guided by the extensive tests of the code performance made in the mentioned papers (Marino et al. 2018; Cantalupo et al. 2019; Mackenzie et al. 2019; Lofthouse et al. 2020), with the scope of exploiting CubEx for the detection of the faintest $Ly\alpha$ nebulae. However, we also tested the code with a variety of combinations of parameters before choosing the best set that allows a fidelity of 60% at S/N larger than 11.

We tested the detection of sources in the mosaic as well the detection of fake sources in the negative of the mosaic cube, by changing the spatial and spectral S/N of the individual voxels (values of 2, 3, 7, and 10), the integrated S/N (values of 5, 7, and

<figure><img src="image_4.png" /></figure>

**Fig. 2.** Detection fidelity as a function of signal-to-noise ratio. The fidelity is expressed as $1 - N_{\text{neg}}/N_{\text{pos}}$, where $N_{\text{neg}}$ is the number of the sources detected in the datacubes multiplied by -1 and $N_{\text{pos}}$ in the number of the detections in the original datacubes, when we apply the list of detection criteria described in the text. The red stars and line correspond to the detections in the “mosaic”, while the other lines correspond to the individual pointings, magenta for P1, blue for P2, cyan for P3, and green for P4. The signal-to-noise ratio in the x-axis correspond to the ratio between the isophotal flux and its error as provided by CubEx. The shapes of the lines reflect the uncertainty in the fidelity measurement.

10), and the minimum number of connected voxels (values of 5, 10, 30, and 50). We found that with a spatial S/N larger than 6, we only detect a few sources with very high isophotal S/N. For a S/N of connected voxels lower than 3, we detect more fake sources than real sources. In the case the integrated S/N is larger than 8, we miss 80% of the real sources for an isophotal S/N of 11. We found that the minimum number of voxels does not change the fidelity fraction more than 5%, however, for a value larger than 30 and for a spatial and spectral S/N of 3, we miss 15% of sources. Therefore, we fixed this parameter to 30.

As described in the following section, we kept a few lower S/N sources, with a counterpart in the HST F160W image and which $Ly\alpha$ emission line seemed convincing in the 1D spectrum. These sources make the “Supplemental” sample (Table 2). For one of the supplemental sources, the counterpart in the F160W image is also present in the COSMOS2015 catalog (Laigle et al. 2016) with a redshift consistent with that of the $Ly\alpha$ detection (see Table 2). However, given the high density of sources in the HST image, we calculated the probability of chance alignment for all the supplemental sources. We first ran Source Extractor (Bertin & Arnouts 1996) with the F160W image as detection and a detection threshold of $2\sigma$. In this run, we identified the counterparts of all the three sources in the Supplemental sample. Then, we calculated the surface density, $n(<\text{mag})$, of the sources detected in the F160W catalog, brighter than those counterparts, and obtained $n(m < m_{446}) = 0.09 \text{ arcsec}^{-2}$, $n(m < m_{414}) = 0.07 \text{ arcsec}^{-2}$, and $n(m < m_{199}) = 0.04 \text{ arcsec}^{-2}$. By following the discussion in Downes et al. (1986) (and references therein), we estimated the probability of chance alignment as $P = 1 - \exp(-n\pi R^2)$, where n is the surface density (sources per arcsec$^2$), $\pi R^2$ is the searching area of the counterparts, and R is the searching radius in arcsec given by the MUSE PSF. The probabilities resulted $P_{446} = 0.09$, $P_{414} = 0.07$, and $P_{199} = 0.04$. A chance probability of less than 10% supports the idea that the low signal-to-noise sources are unlikely to be chance
<table><caption>**Table 2.** Properties of the sources detected in the MUSE datacube and compatible with $z \approx 5.3$.</caption><thead><tr><th>ID<br>(1)</th><th>RA<br>(2)</th><th>Dec<br>(3)</th><th>λ<sub>Lya</sub><br>(Å)<br>(4)</th><th>F(Lya)<br>(10<sup>-18</sup> erg s<sup>-1</sup> cm<sup>-2</sup>)<br>(5)</th><th>S/N<br>(6)</th><th>zLya-z[CII]<sup>(c)</sup><br>(km s<sup>-1</sup>)<br>(7)</th><th>M<sub>*</sub><br>(10<sup>9</sup> M<sub>☉</sub>)<br>(8)</th><th>M<sub>CO</sub><br>(10<sup>9</sup> M<sub>☉</sub>)<br>(9)</th><th>M<sub>dys</sub>([CII])<br>(10<sup>9</sup> M<sub>☉</sub>)<br>(10)</th><th>M<sub>dust</sub><br>(10<sup>8</sup> M<sub>☉</sub>)<br>(11)</th><th>SFR<sup>(d)</sup><br>(M<sub>☉</sub> yr<sup>-1</sup>)<br>(12)</th><th>A<sub>V</sub><sup>(d)</sup><br>(13)</th><th>IAU name<br>(14)</th></tr></thead><tbody><tr><td colspan="15"><strong>Main</strong></td></tr><tr><td>mosaic_1513 <sup>(a)</sup></td><td>150.0864</td><td>2.5891</td><td>7663.2</td><td>16.92 ± 0.37</td><td>46</td><td>380 ± 5</td><td>20, 0.95</td><td>57 ± 5, –</td><td>97 ± 16, –</td><td>2.66 ± 0.76, –</td><td>1100 <sup>(e)</sup>, 18</td><td>0.8, 0.0</td><td>LAEJ100020.74+023520.8</td></tr><tr><td>mosaic_1496 <sup>(b)</sup></td><td>150.0897</td><td>2.5865</td><td>7662.7</td><td>5.63 ± 0.26</td><td>22</td><td>540 ± 10</td><td>14.8 ± 2.0 <sup>(f)</sup></td><td>–</td><td>19.0 ± 5.0</td><td>&lt;0.9</td><td>&lt;54 <sup>(e)</sup></td><td>0.2 <sup>(f)</sup></td><td>LAEJ100021.53+023511.4</td></tr><tr><td>mosaic_1548 <sup>(g)</sup></td><td>150.0909</td><td>2.5770</td><td>7666.5</td><td>5.42 ± 0.34</td><td>16</td><td></td><td>12.2 ± 2.2</td><td></td><td></td><td></td><td>15</td><td>0.0</td><td>LAEJ100021.82+023437.2</td></tr><tr><td>mosaic_1520 <sup>(h)</sup></td><td>150.0874</td><td>2.5938</td><td>7657.9</td><td>1.89 ± 0.17</td><td>11</td><td></td><td>2.8 ± 3.8</td><td></td><td></td><td></td><td>130</td><td>0.7</td><td>LAEJ100020.98+023537.7</td></tr><tr><td>mosaic_770</td><td>150.1012</td><td>2.5915</td><td>7634.5</td><td>2.60 ± 0.24</td><td>11</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>LAEJ100024.29+023529.4</td></tr><tr><td>mosaic_1035</td><td>150.0823</td><td>2.5959</td><td>7640.7</td><td>2.53 ± 0.23</td><td>11</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>LAEJ100019.75+023545.2</td></tr><tr><td>P3_547</td><td>150.0895</td><td>2.5842</td><td>7643.9</td><td>1.60 ± 0.14</td><td>11</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>LAEJ100021.48+023503.1</td></tr><tr><td colspan="15"><strong>Supplemental</strong></td></tr><tr><td>P2_446</td><td>150.0852</td><td>2.5877</td><td>7643.7</td><td>1.68 ± 0.19</td><td>9</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>LAEJ100020.45+023515.7</td></tr><tr><td>mosaic_414</td><td>150.0924</td><td>2.5752</td><td>7667.8</td><td>0.96 ± 0.10</td><td>10</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>LAEJ100022.18+023430.7</td></tr><tr><td>mosaic_199 <sup>(i)</sup></td><td>150.0860</td><td>2.5881</td><td>7662.7</td><td>5.14 ± 1.15</td><td>5</td><td></td><td>0.5 ± 0.4</td><td></td><td></td><td></td><td>23</td><td>0.2</td><td>LAEJ100020.64+023517.2</td></tr></tbody></table>

**Notes.** Lyα and physical properties of the galaxies detected in the MUSE datacube. In the first part of the table, we list the sources in the 'Main' sample, in the second part the "Supplemental" sample. The columns correspond to (1) ID of the sources, (2) their RA and (3) declination, (4) wavelength of the detected emission line which is interpreted as Lyα, (5) isophotal integrated flux of the detected Lyα emission line, (6) signal-to-noise ratio of the detected emission line, (7) velocity difference between the redshift implied by the detected Lyα central wavelength and the systemic redshift inferred by the [CII] emission line, (8) stellar, (9) molecular-gas, (10) dynamical, and (11) dust mass, (12) star-formation rate, (13) dust extinction in magnitude, (14) name of the sources according to the IAU standard. The mass estimates are reported only for the sources with CO and/or [CII] detections (Capak et al. 2011; Riechers et al. 2010, 2014), and for the sources with a counterpart in the COSMOS2015 catalog (Laigle et al. 2016). Star-formation rate and dust extinction are reported for the sources with a counterpart in the COSMOS2015 catalog. <sup>(a)</sup>System composed by AzTEC-3 and LBG-3. LBG-3 is named 1447526 in Capak et al. (2011), COLDz.COS.0 in Riechers et al. (2020), and COSMOS2015_849732 in the COSMOS2015 catalog. The physical parameters corresponding to AzTEC-3 (first number) and LBG-3 (second number) are taken from Riechers et al. (2014, 2020). <sup>(b)</sup>LBG-1 as defined in Riechers et al. (2014), 1447524 in Capak et al. (2011), HZ6 in Capak et al. (2015), and COSMOS2015_848185 in the COSMOS2015 catalog. The parameter values are taken from Riechers et al. (2014), Pavesi et al. (2019), Faisst et al. (2020) (see also Barovic et al. 2017, for HST and ALMA measurements). <sup>(c)</sup>Velocity reported only for the sources with [CII] detections. <sup>(d)</sup>SFR and A<sub>V</sub> from the best fit SED fitting that provides the photometric redshift in the COSMOS2015 catalog, except for AzTEC-3 and LBG-1. <sup>(e)</sup>SFR from FIR. <sup>(f)</sup>Capak et al. (2011) reported M<sub>*</sub> = 3.2×10<sup>9</sup> M<sub>☉</sub> and A<sub>V</sub> = 0.5. <sup>(g)</sup>It overlaps with COSMOS2015_841844 (z<sub>phot</sub> ≈ 5.4 ± 0.1). <sup>(h)</sup>It overlaps with COSMOS2015_852474 (z<sub>phot</sub> ≈ 5.4 ± 0.9). <sup>(i)</sup>It overlaps with COSMOS2015_848724 (z<sub>phot</sub> ≈ 5.1 ± 1.1).

detections. However, we treat them as candidates in the estimation of the density of the protocluster in Sect. 3.4.

## 3. Results

We present here the list of Lyα emitting sources found in our MUSE datacube with the method described in Sect. 2.2. In addition, we also discuss sources from the literature with photometric redshift consistent with $z \approx 5.3$ that could show a Lyα emission in the MUSE datacube, but that do not match our detection criteria.

### 3.1. Line detections compatible with Lyα at $z \approx 5.3$

By using the method described in Sect. 2.2, we identified ten sources with Lyα emission line compatible to be at $z \approx 5.3$. Eight of them are detected in the mosaic, one is detected in P2, and one is detected in P3. In Table 2, we report the properties of the ten sources, naming them based on the detection cube. The first seven sources make our Main list; mosaic_1513 corresponds to the system composed of AzTEC-3 and the Lyman break galaxy (LBG-3) located at 2" north of it; mosaic_1496 corresponds to LBG-1; mosaic_1513 and mosaic_1496 are the most extended we have in our list of Lyα emitting sources and are characterized by previous ALMA [CII] detections. They deserve a special discussion (see the following sections). In addition, mosaic_1548, mosaic_1520, and mosaic_199 overlap with sources in the COSMOS2015 catalog, where their photometric redshifts agree with the Lyα redshift.

In Fig. 3, we show the spatial distribution of the ten Lyα detections on the RGB image corresponding to the MUSE mosaic. The inserts show their Lyα emission line in velocity

space with respect to the systemic redshift of AzTEC-3. As we can see in the inserts, the ranges in spectral pixels of the CubEx detections are large for a few sources. In particular, for the mosaic_1548 source, CubEx identified 12 spectral pixels which correspond to a line with a tail.

In Fig. 4, we show the 1D MUSE spectra extracted with CubEx from an area corresponding to the Lyα detection toward AzTEC-3 and LBG-1. These are the only sources in the field with a known systemic redshift from [CII] (Riechers et al. 2014). In the next section and in the appendix, we show the optimally extracted narrow-band images of all the ten Lyα detections.

### 3.2. Sources at $z \approx 5.3$ with HST counterparts

To find additional Lyα emitting sources that might have escaped our detection threshold, we also selected sources with photometric redshifts consistent with that of AzTEC-3 ($4.9 < z_{phot} < 5.6$ which is $z = 5.3 \pm 3\sigma$ given the typical $z_{phot}$ uncertainty at $3 < z < 6$) in the COSMOS2015 catalog and extracted their MUSE spectra within a fixed aperture of 0.7" radius (about twice the mosaic PSF size). Five sources in addition to the ones mentioned in the previous section did not show any Lyα emission in the MUSE spectra.

We also investigated a method to identify a $z \sim 5$ star-forming galaxy from the HST F606W, F814W, F105W, F125W, and F160W images. As we can see in Fig. 5, such a source is characterized by a sharp decrement between the F606W and F814W filters. We considered representative spectra of star-forming galaxies with Lyα in absorption and in emission (Shapley et al. 2003) and we convolved them with the HST filter transmission curves. We found that a star-forming galaxy at $z \sim 5.3$ is characterized by a $V_{606} - I_{814}$ color of about 3, given the non detection in F606W, and a flat spectrum redder than the
<figure><img src="image_5.png" /></figure>

**Fig. 3.** RGB image composed of the HST images in F814W, F125W, and F160W filters covering the area of our MUSE observations. The overlaid contours represent the $3\sigma$ levels of the $Ly\alpha$ emission. Red contours correspond to the sources with counterparts in the COSMOS2015 catalog. Green contours correspond to the $Ly\alpha$ detections without counterparts. The inserts contain the normalized $Ly\alpha$ profiles (normalized flux density versus velocity) of the $Ly\alpha$ detections indicated in the title of the inserts. The wavelengths of the maximum of the line detected by CubEx is indicated with a vertical dashed blue line. The blue shaded areas correspond to the spectral pixel ranges of the lines detected by CubEx. The zero velocity corresponds to the systemic redshift of the [CII] detection of AzTEC-3 and it is indicated as a vertical red line in each insert.

<figure><img src="image_6.png" /></figure>

**Fig. 4.** Left panel: 1D spectrum of the AzTEC-3+LBG-3 system in the $Ly\alpha$ wavelength range. The original-sampling spectrum is shown in blue with error bars, the one smoothed by 3 × 3 wavelength channels in black. The error spectrum comes from the rescaled variance cube. The magenta spectrum is the theoretical sky spectrum from Rousselot et al. (2000) in arbitrary units. Right panels: same as the left panel for LBG-1.

F814W filter. We ran Source Extractor in dual mode, with the F105W image as a detection and the other HST-filter images as measurement images. We chose a detection minarea equal to the PSF, a detection threshold of 2, and an optimal extraction aperture given the image PSF. Also, we asked for a signal-to-noise ratio of 3 in the detection flux. We found one source with color consistent with a $z \sim 5.3$, that is also contained in the COSMOS2015 catalog with $z_{phot} = 5.2$ (COSMOS2015_842471). However, the MUSE spectrum does not show any significant emission line at the wavelength of $Ly\alpha$.

The best photometric filter to identify the members of the AzTEC-3 protocluster would be centered at ~7600 Å (the wavelength of $Ly\alpha$ at $z \sim 5.3$). Sobral et al. (2018) studied the galaxies detected in the I767 filter and compiled a sample of $Ly\alpha$ emitters at $z \approx 5.3$. However, all the sources in their sample are located outside our MUSE coverage.

Along the same line, we extracted the MUSE spectra of the galaxies discovered by Capak et al. (2011) in the AzTEC-3 overdensity. The ones showing $Ly\alpha$ emission at $z = 5.3$ in our MUSE datacube are LBG1447524 and LBG1447526 (LBG-1 and LBG-3).

### 3.3. Other catalogs from the literature

We also inspected other catalogs from the literature covering the AzTEC-3 area, the catalog of CO detections from
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** HST filter transmission curves (ACS F606W in yellow, ACS F814W in cyan, WFC3 F105W in magenta) and representative spectra of star-forming galaxies at $z \sim 5.3$ (Ly$\alpha$ absorber in red, Ly$\alpha$ emitter with a small, intermediate, and large EW(Ly$\alpha$) in green, blue, black, respectively from Shapley et al. 2003). We corrected all the representative spectra for the intergalactic-medium absorption at $z = 5.3$ by assuming the formalism in Madau (1995), before using them.</figcaption></figure>

Pavesi et al. (2018a) and the radio-continuum detection catalog from Smolčić et al. (2017) and van der Vlugt et al. (2021). We required non detection in the F606W image and we found only one source with Ly$\alpha$ at $z \sim 5.3$ in both catalogs, corresponding to the AzTEC-3 galaxy. In the Pavesi’s catalog there is a source at $z = 5.3$ (in addition to AzTEC-3) without Ly$\alpha$ emission from our MUSE observation that is shown in Fig. 6 (see next section).

### 3.4. Space distribution of all the candidates at $z \approx 5.3$

In Fig. 6, we show the location in the RA–Dec plane of the Ly$\alpha$ emitting sources detected in the MUSE datacube and of the sources from the literature with photometric redshift consistent with being at $z \approx 5.3$.

To verify if the group of sources around AzTEC-3 constitutes an overdensity of galaxies and confirms the presence of a protocluster, we compare them with the catalog of Ly$\alpha$ emitting sources detected with MUSE in the *Hubble* Ultra Deep Field (Inami et al. 2017, mosaic datacube). For a minimum Ly$\alpha$ luminosity of $4 \times 10^{41}$ erg s$^{-1}$, there are 77 galaxies at $4.9 < z < 5.6$ (redshift bin of 0.7) in the mosaic area of $3' \times 3'$ (30 800 sources per deg$^2$), which means 2200 sources per deg$^2$ in the redshift range corresponding to the peak of our Ly$\alpha$ line detections ($5.28 < z < 5.33$, redshift bin of 0.05). At the same limiting Ly$\alpha$ luminosity, we detected eight emitters in the Main sample, taking AzTEC-3 and LBG-3 as two separated galaxies, and 11 if we also consider the Supplemental sample. Therefore, we found $9 \pm 2$ emitters in our field ($=4 \times 10^{-4}$ deg$^2$), that can be translated into $9 \pm 3$ emitters if we take also into account the 40% chance of fake detections due to the fidelity-cut criterion in the uncertainty budget or $22\,500 \pm 7500$ sources per deg$^2$. This implies an overdensity of $10 \pm 3$. We, also, considered all the sources at $4.9 < z_{phot} < 5.6$ in the COSMOS2015 catalog located outside the area covered by our MUSE observations and found a density of $1862 \pm 300$ galaxies per deg$^2$. The overdensity of Ly$\alpha$ emitters in the MUSE area is $12 \pm 5$ times the density of all the sources in the COSMOS2015 catalog, with and without Ly$\alpha$ emission.

<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Location of the Ly$\alpha$ emitting sources detected in the AzTEC-3 environment. The ones with Ly$\alpha$ central wavelength larger (red) and smaller (blue) than that of the AzTEC-3 system are shown as small circles. The AzTEC-3 submillimeter and the LBG-3 galaxies are shown as yellow star and yellow diamond. Sources from the literature with photometric redshifts consistent with $z \sim 5.3$ are also shown, five COSMOS2015 sources with $4.9 < z_{phot} < 5.6$ (cyan big circles), Lyman break galaxies listed in Capak et al. (2011), but without Ly$\alpha$ emission in our MUSE datacube (green squares), one CO detection from Pavesi et al. (2018a) at $z = 5.3$ (gray circle). The black cross is the center of coordinates of all the sources shown in the figure and the gray cross the center obtained without considering AzTEC-3. The blue is the center of coordinates only of the sources with stellar mass larger than $10^9 M_{\odot}$ and the red cross is the weighted mass barycenter of the Ly$\alpha$ emitting sources detected in the field. The green cross is the weighted mass barycenter of the sources with a stellar mass estimation but without considering AzTEC-3.</figcaption></figure>

These estimations confirm the presence of an overdense region around AzTEC-3.

For comparison, the SSA22 protocluster (Steidel et al. 2000), which also contains a few SMGs, is an overdensity of about six at $z \sim 3.1$. The overdensity of Ly$\alpha$ emitters discovered by Shi et al. (2019) contains about four times more galaxies than the field. In the *Hubble* Deep Field, Walter et al. (2012), Calvi et al. (2021) found an overdensity around the submillimeter galaxy HDF850.1 at $z = 5.183$. They presented 23 spectroscopically confirmed star-forming galaxies, including Ly$\alpha$ emitters, at $z = 5.2$ with a surface density of at least a factor of two higher than the field. None of the sources, in addition to HDF850.1, are detected at submillimeter wavelengths. Pavesi et al. (2018b) discovered an overdensity at $z \sim 5.7$ similar to that around AzTEC-3. At the center of the overdensity there is a dusty starburst, probably in the phase of an on-going merger. A star-forming galaxy as massive as LBG-1 is located at $13''$ away from the starburst. An overdensity of Ly$\alpha$ emitters was discovered in the same area with an overdensity parameter of about ten. This and our area of observations are therefore ones of the regions at $z > 5$, characterized by the highest overdensities in the COSMOS field.

Other overdensities of Ly$\alpha$ emitters have been discovered in the literature at $4 < z < 6$ around radio galaxies or around quasars expected to trace dense regions of the Universe. For instance, Venemans et al. (2002) discovered that the region around the luminous radio galaxy TN J1338–1942 at $z = 4.1$ is 15 times more dense than the field; Venemans et al. (2004) found a density up to six times higher than in the field around
<figure><img src="image_9.png" /></figure>

**Fig. 7.** *Upper left:* signal-to-noise optimally extracted narrow-band image (Borisova et al. 2016) of the continuum-subtracted $Ly\alpha$ emission toward the AzTEC-3+LBG-3 system, rescaled to surface brightness units, and smoothed using a 2-pixel Gaussian kernel. The narrow band is built considering the 13 wavelength channels of the detection. The 5 white contours reveal the $1.3 \sim 3.8 \text{ erg s}^{-1} \text{ cm}^{-2} \text{ arcsec}^{-2}$ levels, which correspond to $3-10\sigma$ in this image. The black letters indicate the directions we refer to in the text (N = north, NW = northwest, S = south, SE = southeast, W = west of the SMG). The small black segment indicates the $1''$ scale. We also show the continuum-subtracted $Ly\alpha$ surface brightness contours overplotted to the HST ACS F606W (upper right), F814W (lower left), and F160W (bottom right) images. The $0.5''$-side square indicates the center of the F160W position of LBG-3. The PSF of our MUSE data is $0.7''$, corresponding to about 4 kpc or 27 ckpc at $z = 5.3$. The eight red contours correspond to the [CII] emission (from $3\sigma$ to $10\sigma$ where $1\sigma = 2.38 \times 10^{-4} \text{ Jy beam}^{-1}$) at 301.671–302.142 GHz, corresponding to $-230$ to $230 \text{ km s}^{-1}$ from the peak emission. The synthesized beam of the [CII] ALMA observations is $0.63'' \times 0.56''$ as reported in Riechers et al. (2014). The $3\sigma$ contour extends toward LBG-3.

the radio galaxy TN J0924–2201 at $z = 5.2$; Zheng et al. (2006) found an overdensity on the order of six around the radio-loud quasar SDSS J0836+0054 at $z = 5.8$.

Assuming that all the galaxies shown in Fig. 6 belong to the protocluster, we can see that they follow the location of our $Ly\alpha$ emitting sources. The center of coordinates of the most massive sources with stellar mass larger than $10^9 M_{\odot}$ and the weighted mass barycenter of the sources with a stellar mass estimation, even without considering AzTEC-3, are close to the submillimeter galaxy on the side of its $Ly\alpha$ emission peak. This implies that the region of the protocluster where the gravitational potential is the strongest is close to the SMG, which in turn is very close to the center of the protocluster.

# 4. Spectroscopic properties from the MUSE spectra

## 4.1. AzTEC-3

In Fig. 7, we show the continuum subtracted $Ly\alpha$ surface brightness toward AzTEC-3 corresponding to 7658–7670 Å. The original narrow band image and the continuum subtracted one are produced within CubEx optimizing the signal-to-noise ratio, as described in Borisova et al. (2016). The two main peaks of the emission roughly correspond to the positions of AzTEC-3 and the LBG-3. The emissions toward AzTEC-3 and LBG-3 are

blended at the $5\sigma$ level. The entire emission occupies an area of $4.4'' \times 6.3''$ (about $27 \times 38 \text{ kpc}^2$ and $170 \times 240 \text{ ckpc}^2$).

As expected for a source at $z \sim 5$, the sources are not visible in the F606W filter and the $Ly\alpha$ wavelength is contained in the F814W filter. The peak of the $Ly\alpha$ emission toward the SMG is offset toward the southeast with respect to the star-formation knots seen in the HST images and traced by the [CII] detection. This may be related to the presence of dust on the main knot of star formation of the SMG (see next section). The LBG peak is offset by $0.5''$ (20 ckpc) with respect to the source seen in the HST images and elongated to the north. Since the $Ly\alpha$ emission is very sensitive to the even small presence of dust, this is in agreement with the idea that on the main star-formation knot of LBG-3 some of the $Ly\alpha$ photons could be absorbed by dust (see next section).

Even if the $Ly\alpha$ emissions toward the SMG and the LBG seem connected, the spatial shape of the highest S/N contours indicates that the peak $Ly\alpha$ emission of the SMG is quite symmetrical. However, that of the LBG seems to be elongated and tilted toward the north (N position in the figure) and extends up to $2''$ (almost 80 ckpc). Also, the $Ly\alpha$ emission of the SMG shows extended components toward various directions. It is worth noting that a low S/N emission from the SMG toward the direction of LBG-3 was also observed in the [CII] observations, despite the main [CII] flux was found to be concentrated in a
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Lyα channel maps toward the AzTEC-3+LBG-3 system. Each panel corresponds to the indicated wavelength and velocities with respect to the systemic redshift inferred from the [CII] emission line. The white contours represent the isophotes at $S/N = 4, 5, \text{and } 6$. The red square corresponds to the F160W position of LBG-3 and the red circle to the location of the peak of the [CII] emission of AzTEC-3. The green contours in the fifth panel represent the contours of the continuum-subtracted Lyα surface brightness. The horizontal gray-color bar shows a signal-to-noise ratio from 3 to 5. We refer to Fig. 7 for the 1" scale.</figcaption></figure>

compact position on top of the main knot of star formation, as described in Riechers et al. (2014).

In Fig. 8, we show the Lyα channel maps of the AzTEC-3+LBG-3 system. Each channel is separated by approximately 1.5 Å. The figure shows that the Lyα emissions toward the two galaxies peak at similar wavelengths, 7662 Å for LBG-3 and 7664 Å for AzTEC-3, respectively. The Lyα emission associated to the SMG is visible in all the channels at $4\sigma$, while that probably associated to LBG-3 is significant at 7660–7665 Å. At 7665 Å, the emission is mainly concentrated on top of the SMG, but shows a kind of a bridge between SMG and LBG, which could represent a region of interaction between the two close-by galaxies (see Sect. 6). The Lyα emission of the bridge is observed from 7660 to 7666 Å, with a maximum at 7664 Å. It occupies a wide range in Lyα wavelengths, that could correspond to a region of a large range of gas velocities, together with large HI column densities. These gas characteristics could be the result of the interactions in the AzTEC-3+LBG-3 system.

At the position of the main star-forming knot of the SMG and of the peak of the Lyα emission associated to the SMG, the Lyα emission extends from 200 to 600 km s⁻¹. The velocity could be associated to a star-formation driven outflow of neutral gas, departing from the SMG and mixing with the interacting gas in the bridge between SMG and LBG. This interaction could produce some of the extended tail seen also in the Lyα 1D spectrum. The blue shift of low-ionization absorption lines, directly tracing the kinematics of the HI gas, together with the radiative-transfer modeling of the Lyα emission (Sect. 6), sensitive to kinematics and amount of HI gas, could provide support of this hypothesis.

At 7662 Å, the Lyα emission appears elongated to the north of LBG-3 (N position in Fig. 7) and to the southeast of AzTEC-3 (SE position in Fig. 7). Above LBG-3, the Lyα emission is com-

pact in wavelength. The Lyα emission is compact in wavelength (around 300 km s⁻¹) also in the region below the peak of the Lyα emission associated to the SMG. In this position, the gas kinematics could be associated with neutral gas flows between the SMG and the mosaic_199 source (Table 3).

Since Lyα photons can be scattered by HI gas and so reveal its presence, the elongated shape of the Lyα emission to the north of LBG-3 (N position in Fig. 7), in the bridge between SMG and LBG, and to the southeast of the SMG (SE position in Fig. 7) could be the result of the interaction between AzTEC-3, LBG-3, and possibly mosaic_199. As shown in the simulations described in Yajima et al. (2013), at the first passage of a merging event between star-forming galaxies, the Lyα emission coming from the progenitors can be intense due to the triggered star-formation event and its shape could also follow the gas distorted during the interaction in the outer regions.

In Fig. 4, we can see the integrated 1D spectrum of the AzTEC-3+LBG-3 system in the Lyα wavelength ranges, free of strong sky line residuals. Neither metal nor AGN-diagnostic lines are detected in the spectrum. As expected from Fig. 8, the main peak of the Lyα line is redshifted with respect to the systemic redshift by about 400 km s⁻¹. The entire line including its red tail occupies a velocity range up to 900 km s⁻¹. As shown in Verhamme et al. (2006), the shift of the Lyα red peak can be 2 or 3 times the velocity of the outflowing gas depending on the HI column density conditions. Riechers et al. (2014) measured that the central velocity of the OH163 µm doublet from the ALMA observations was blueshifted of about 100 km s⁻¹ with respect to the [CII] line and this could be associated to an outflow, maybe produced by the AzTEC-3 starburst.

The integrated flux of the Lyα emission of the entire system is (16.92 ± 0.37) × 10⁻¹⁸ s⁻¹ cm⁻². By estimating a continuum on the order of 3 × 10⁻¹⁹ erg s⁻¹ cm⁻² Å⁻¹ on the red side of
<table><caption>**Table 3.** Best fit parameters of the zELDA models.</caption><thead><tr><th>Spectrum</th><th>V<sub>exp</sub><br>(km s<sup>-1</sup>)</th><th>NHI<br>(10<sup>19</sup> atoms cm<sup>-2</sup>)</th><th>&chi;<sup>2</sup>/d.o.f.</th></tr></thead><tbody><tr><td>AzTEC-3+LBG-3<sup>1</sup></td><td>21[15–28]</td><td>19[18–22]</td><td>0.92</td></tr><tr><td>AzTEC-3+LBG-3<sup>2</sup></td><td>20[3–93]</td><td>1.3[0.1–7.8]</td><td>1.65</td></tr><tr><td>AzTEC-3+LBG-3<sup>3</sup></td><td>25[20–30]</td><td>19[18–20]</td><td>0.48</td></tr><tr><td>onlySMG Ly&alpha; <sup>(a)</sup></td><td>42[34–56]</td><td>16[14–18]</td><td>1.14</td></tr><tr><td>onlySMG Ly&alpha; <sup>(b)</sup></td><td>785[754–825]</td><td>32[22–46]</td><td>2.24</td></tr><tr><td>onlySMG UVknots <sup>(a)</sup></td><td>60[40–127]</td><td>0.3[0.1–0.8]</td><td>0.89</td></tr><tr><td>onlySMG UVknots <sup>(b)</sup></td><td>13[1–71]</td><td>0.1[0.02–2]</td><td>1.39</td></tr><tr><td>onlyLBG3 <sup>(a)</sup></td><td>15[6–23]</td><td>21[18–28]</td><td>0.93</td></tr><tr><td>bridgeSMGLBG <sup>(a)</sup></td><td>15[7–25]</td><td>19[17–28]</td><td>0.63</td></tr><tr><td>LBG-1 <sup>(a)</sup></td><td>12[4–30]</td><td>68[53–102]</td><td>1.09</td></tr><tr><td>LBG-1 <sup>(b)</sup></td><td>83[25–293]</td><td>24[2–92]</td><td>1.07</td></tr><tr><td>onLBG1horizontally <sup>(a)</sup></td><td>10[2–28]</td><td>70[48–107]</td><td>0.54</td></tr><tr><td>onLBG1tilted <sup>(a)</sup></td><td>14[4–30]</td><td>60[43–94]</td><td>0.47</td></tr><tr><td>belowLBG1 <sup>(a)</sup></td><td>9[1–58]</td><td>2.2[0.1–37]</td><td>0.38</td></tr><tr><td>aboveLBG1 <sup>(a)</sup></td><td>12[2–28]</td><td>69[45–103]</td><td>0.41</td></tr><tr><td>mosaic_1548 <sup>(a)</sup></td><td>30[20–46]</td><td>0.1[0.02–0.3](z = 5.304)</td><td>0.54</td></tr><tr><td>mosaic_1548 <sup>(inflow)</sup></td><td>10[4–30]</td><td>0.8[0.1–1.8](z = 5.308)</td><td>0.46</td></tr><tr><td>mosaic_1520 <sup>(a)</sup></td><td>26[6–61]</td><td>0.8[0.1–8](z = 5.298)</td><td>0.36</td></tr></tbody></table>

**Notes.** Expansion velocity, $V_{exp}$ and HI column density, NHI of the zELDA models that correspond to the best fits of the spectra in the first column. The last column provides an estimation of the goodness of the fit in terms of reduced $\chi^2$. The numbers in the second and third columns correspond to the 50th percentile, while the two numbers in parentheses are the variables corresponding to the 16th and 84th percentiles of the parameter space of the $V_{exp}$ and the NHI values of the models. AzTEC-3+LBG-3¹ corresponds to the fit performed at wavelength larger than the systemic redshift (green curve in Fig. 14a); AzTEC-3+LBG-3² at wavelength larger than the maximum of the main peak (blue curve in Fig. 14a); AzTEC-3+LBG-3³ corresponds to fit performed on the spectrum corrected for IGM absorption based on the prescription of Madau (1995) (red curve in Fig. 14a). <sup>(a)</sup>Corresponds to the fit of the main peak. <sup>(b)</sup>Corresponds to the fit of the extended tail. <sup>(inflow)</sup>Proposed best fit solution with inflow. The fit of the last three models was performed leaving the redshift as a free parameter. The best fit systemic redshift is indicated in parenthesis.

Ly$\alpha$ by linear fit, we calculate a rest-frame EW(Ly$\alpha$) = 13 ± 2 Å. The FWHM(Ly$\alpha$) of the main peak excluding the long tail is 220 km s$^{-1}$.

In Fig. 9, we show the spectroscopic comparison between [CII] and Ly$\alpha$ emission in velocity space. It is worth remembering that the Ly$\alpha$ emission contains contributions from AzTEC-3 and also from LBG-3, while the [CII] emission is only concentrated at the location of the SMG as seen in the HST images. The zero velocity is given by the [CII] systemic redshift. Both AzTEC-3 and LBG-3 contribute to the main peak and to the extended tail (at more than 500 km s$^{-1}$) of the integrated Ly$\alpha$ emission, even if the tail is mostly dominated by the SMG.

We separated the 1D Ly$\alpha$ spectrum of the SMG and of LBG-3 in two ways, with CubEx and extracting a fixed-aperture spectrum from the MUSE datacube. To make CubEx detect two separated (not blended) sources, we increased the S/N of the connected pixels from 3 to 5, the spatial and spectral S/N from 3 to 7, and decreased the minimum number of voxels from 30 to 5 in the detection parameters. This way CubEx is able to detect the brightest Ly$\alpha$ emission associated only to the SMG and that associated only to the LBG-3. Also, we extracted MUSE spectra within a 0.7"-radius aperture and obtained the spectra at the position of the brightest Ly$\alpha$ regions only coming from the SMG (onlySMG Ly$\alpha$ spectrum) and only from the LBG (onlyLBG spectrum). The profile of the Ly$\alpha$ emission just coming from the SMG presents a main peak which is fainter and broader than the main peak of the emission coming from the entire AzTEC-3+LBG-3 system (left panel of Fig. 9). The Ly$\alpha$ profile, obtained from an aperture including all the rest-frame UV knots of the SMG (onlySMG UV knots spectrum) shows a sharper slope on

the red side. However, the overall shapes of the Ly$\alpha$ emissions are comparable.

Even if we focus on the brightest part of the Ly$\alpha$ emission only associated to the SMG (as the [CII] emission is), we see that the Ly$\alpha$ main peak is redshifted with respect to the [CII] by the same amount as the entire AzTEC-3+LBG-3 system.

It is worth noting that Riechers et al. (2014) estimated a star-formation rate per unit of area of 530 $M_{\odot}$ yr$^{-1}$ kpc$^2$ (see also Riechers et al. 2020, for a discussion), which is larger than the limit of 0.1 $M_{\odot}$ yr$^{-1}$ kpc$^{-2}$, satisfied by local starbursts and high-z Lyman Break galaxies to sustain starburst-driven outflows (Heckman 2001). Starburst-driven outflows could favor the escape of Ly$\alpha$ photons even from dusty interstellar media (e.g., Kunth et al. 1998; Verhamme et al. 2006), as that of the SMG, and they could play an important role in producing the observed Ly$\alpha$ surface brightness from AzTEC-3. At 34 ckpc (0.9") from the star-forming region traced by the [CII] emission peak toward the southeast of the SMG (SE position in Fig. 7), we measure a Ly$\alpha$ surface brightness of ~4.4 × 10$^{-18}$ erg s$^{-1}$ cm$^{-2}$ arcsec$^{-2}$, indicating a $L(Ly\alpha) \sim SB(Ly\alpha) \times 4\pi 0.9^2 \times 4\pi D_L^2(z = 5.3) \sim 1.3 \times 10^{43}$ erg s$^{-1}$.

It is interesting to investigate what could be the main mechanism of the production of the Ly$\alpha$ photons that extend up to at least 30 ckpc. The molecular gas mass of 5.7 × 10$^{10}$ $M_{\odot}$ (Riechers et al. 2020) and the SFR$_{FIR}$ of 1100 $M_{\odot}$ yr$^{-1}$ imply that the starburst could be maintained for a time scale on the order of 50 Myr, as an upper limit assuming that no significant gas mass is lost due to the outflow itself. An outflow of a constant velocity of 800 km s$^{-1}$ would reach a maximum distance of 40 kpc in a time scale of 50 Myr and could have already reached
<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Left panel: Lyα (blue) and [CII] (red) profiles of the AzTEC-3+LBG-3 system in velocity space and normalized to the maximum. The [CII] spectrum is taken from Riechers et al. (2014) and it is associated to AzTEC-3 (see Fig. 7). Right panel: normalized [CII] profile (red), Lyα profile of the emission extracted at the position of the main Lyα peak associated to the SMG (blue and black curves), extracted from the position of the rest-frame UV knots of AzTEC-3 as seen in the F160W image (green), and extracted at the position of the main Lyα peak associated to LBG-3 (cyan). The extraction is done in an aperture of 0.7" radius for the blue, green, and cyan profiles and increasing the S/N of the connected pixels in CubEx for the black profile. For simplicity, we call the spectra as onlySMG Lyα, onlySMG UV knots, and onlyLBG. The normalization of these spectra is performed with the same factor used to normalize the spectrum of the AzTEC-3+LBG-3 system in the left panel to show the different intensities. The normalization of the [CII] profile is also performed accordingly.</figcaption></figure>

part of that distance in a shorter time scale. This outflow could channel the escape of ionizing radiation, ionized, and neutral gas.

Moreover, the Lyα luminosity of $1.3 \times 10^{43}$ erg s$^{-1}$ at a distance of 34 ckpc could be explained in terms of recombination of atoms ionized by a starburst of a thousand solar masses per year. In fact, following Eq. (7) in Cantalupo (2017), the luminosity could be translated into a ionization luminosity of $\sim 1 \times 10^{54}$ photons s$^{-1}$ consistent with that of a starburst of 1000 $M_{\odot}$ yr$^{-1}$ (e.g., Leitherer et al. 1999). Also, following Eq. (5) from the same paper, it can be shown that, for a volume of $4/3\pi(34 \text{ckpc})^3$, a $T = 10^4$ K, and a filling factor of the interstellar medium clouds of $10^{-4}$ (e.g., Cantalupo et al. 2014), a cloud density of $\sim 3$ atoms cm$^{-3}$ would imply the observed $L(Ly\alpha)$ at 34 ckpc. At this density and temperature, recombination is a plausible scenario to explain the observed $L(Ly\alpha)$ and surface brightness.

The luminosity at 90 ckpc would imply a lower cloud density of less than 1 atom cm$^{-3}$ that difficultly explains a recombination mechanism from the starburst. However, at this larger distance the HI scattering (free of dust absorption) could play a role in the case of $NHI > 10^{21}$ atoms cm$^{-2}$. The radiative-transfer-model fit of the 1D integrated Lyα profile of the SMG and of the AzTEC-3+LBG-3 system can give insight on the NHI quantity. The gas at the peak of the Lyα emission and at 90 ckpc from the star-forming region could be the result of the interaction between the SMG and the mosaic_199 source (see below).

## 4.2. LBG-1

Among the star-forming galaxies with Lyα in emission detected in the MUSE datacube, we focus here on LBG-1, because it is the only galaxy, apart from AzTEC-3, for which we have a systemic redshift, inferred from the [CII] emission line (Riechers et al. 2014). In the next section, we consider two additional Lyα detections with photometric redshifts from the COSMOS2015 catalog (Fig. 3).

In Fig. 10, we show the emission detected by CubEx at a wavelength around 7661 Å, corresponding to Lyα toward LBG-1.

This galaxy was studied in Capak et al. (2011), Riechers et al. (2014), Pavesi et al. (2019). Its known physical properties are reported in Table 2. Unlike the UV which clearly shows at least three components, the Lyα emission is smooth and encompasses the three UV components. The integrated $L(Ly\alpha)$ is equal to $(15.9 \pm 0.7) \times 10^{41}$ erg s$^{-1}$ ($(5.5 \pm 0.3) \times 10^{42}$ erg s$^{-1}$ if corrected by dust extinction, following Calzetti et al. 2000 and $A_V = 0.2$) and rest-frame EW(Lyα) = $2.0 \pm 0.2$ Å, by assuming the continuum magnitude ($z_p = 23.71$) from the COSMOS2015 database (the continuum is very noisy and affected by sky line residuals redwards of the Lyα wavelength in the MUSE spectrum).

The size of the entire Lyα emission is $3.2'' \times 3.2''$ (about $20 \times 20$ kpc$^2$ or $125 \times 125$ ckpc$^2$). As shown in Riechers et al. (2014), the [CII] emission at 301.974 GHz encompasses the three knots with a lower signal-to-noise tail toward the southwest (SW position in Fig. 10) as well. At MUSE resolution, we see that significant Lyα flux comes from the three components of LBG-1, maybe indicating that a low dust extinction is allowing Lyα photons to escape also close to the regions of star formation. In Fig. 11, we present the Lyα channel maps toward LBG-1. The figure shows that the Lyα emission is concentrated on the northern region of LBG-1 and that only at the longest wavelengths the emission comes only from the south (S and SW positions in Fig. 10). This could indicate that toward the south, LBG-1 is characterized by larger HI column densities and higher outflow velocities.

To increase the S/N and investigate the origin of the Lyα emission coming from the different directions of LBG-1, we slice the MUSE cube along the pseudo slits shown in the left side of Fig. 12, horizontal on the Lyα emission peak (first row), tilted by 110° on the Lyα emission peak (second row), tilted by ten degrees below the Lyα emission peak (third row), and tilted by ten degrees above the Lyα emission peak (fourth row). These slit orientations are chosen to isolate the Lyα emission coming from the star-forming knots of LBG-1 to that coming from above and below of them. We can see that at the position of the main Lyα peak, the emission extends from 7658 to 7672 Å. The highest velocity regions come from the center toward the southwest
<figure><img src="image_12.png" /><figcaption>**Fig. 10.** *Upper left:* signal-to-noise optimally extracted narrow-band image (Borisova et al. 2016) of the continuum-subtracted $Ly\alpha$ emission, rescaled to surface brightness units, and smoothed using a 2-pixel Gaussian kernel of LBG-1. The five white contours reveal the $0.8 \times 10^{-18} \text{ erg s}^{-1} \text{ cm}^{-2} \text{ arcsec}^{-2}$ levels, which correspond to 3 and $8\sigma$ in this image. The letters indicate reference positions described in the text. The small black segment indicates the 1" scale. We also show the $Ly\alpha$ surface brightness contours overplotted to the F606W (upper right), F814W (lower left) ACS, and F160W WFC3 (lower right) HST images. The PSF of our MUSE data is 0.7", corresponding to about 4 kpc or 27 ckpc at $z = 5.3$. The eight red contours correspond to the $3-10\sigma$ ($1\sigma = 2.38 \times 10^{-4} \text{ Jy beam}^{-1}$) of the [CII] emission at 301.847–302.093 GHz, corresponding to $-120$ to $120 \text{ km s}^{-1}$ from the peak emission. The synthesized beam size of the ALMA [CII] observations is $0.63'' \times 0.56''$ as reported in Riechers et al. (2014).</figcaption></figure>

<caption>**Fig. 11.** $Ly\alpha$ channel map of LBG-1. Each panel corresponds to the indicated wavelength and velocities with respect to the systemic redshift inferred from the [CII] emission line. The white contours represent the isophotes at $S/N = 4, 5, \text{ and } 6$. The red circle shows the F160W position of LBG-1 as in Fig. 1, drawn to drive the eye. The green contours in the *third panel* represent the contours of the continuum-subtracted $Ly\alpha$ surface brightness. The horizontal gray-color bar shows a signal-to-noise ratio from 3 to 5. We refer to Fig. 10 for the 1" scale.</caption>

(S and SW position in Fig. 10). Toward the extreme southwest (SW position in Fig. 10), the $Ly\alpha$ emission reaches $7672 \AA$, corresponding to $400 \text{ km s}^{-1}$ with respect to the main peak wavelength. This value can be related to the gas kinematics, in a wrap of the material of the three knots and/or to high HI col-

umn densities (see Sect. 6). The kinematics of the gas could be related to the interplay of the gas in the merger of the three knots, as suggested by Riechers et al. (2014), the high column density could suggest that the merger is “wet” and/or there is a large gas reservoir supplying the star formation in the entire region.
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** *Left panels:* pseudo slits oriented horizontally on the $Ly\alpha$ emission peak (first row), tilted by 110° on the $Ly\alpha$ emission peak (second row), tilted by ten degrees below the $Ly\alpha$ emission peak (third row), and tilted by ten degrees above the $Ly\alpha$ emission peak (fourth row). *Middle panels:* position-velocity diagrams obtaining slicing the continuum-subtracted MUSE datacube with a pseudo slit oriented as shown in the *left panels*. In the left y-axis, we indicate the wavelength and in the right y-axis the velocity with respect to the systemic redshift inferred from the [CII] emission line, the x-axis corresponds to the position along the pseudo slit in degrees (from the east to the west, except in the second row where it is from the south to the north). The white contours correspond to the 1, 2, 3σ level. *Right panels:* 1D spectra corresponding to the position-velocity diagrams of the *middle panels*. Wavelength is on the y-axis in the same range as in the *middle panels* and the flux normalized to the maximum is on the x-axis. The color coding of these spectra is like in Fig. 13.</figcaption></figure>

By slicing above the main peak, the $Ly\alpha$ emission is confined in wavelength.

In Fig. 13, we show the [CII] and $Ly\alpha$ profiles as a function of the velocity with respect to the [CII] systemic redshift. The integrated 1D spectrum of the $Ly\alpha$ emission line is composed by a narrower peak and a much broader component that reaches the continuum at 7680 Å. The FWHM of $Ly\alpha$ is larger than that of [CII]. This is not unusual in high-redshift $Ly\alpha$ emitting galaxies (e.g., Matthee et al. 2020a) and it could be related to the variety of HI kinematics properties that can condition the formation and escape of $Ly\alpha$ photons. Also, the tail of the $Ly\alpha$ emission could indicate that there is gas in a wide range of kinematic conditions and with different HI column densities around the three star-forming knots. It could be consistent with rotation among the three knots as well. In fact, the wavelength range of this tail corresponds to the feature visible in the second and third rows of Fig. 12 at large velocities and coming from the southwest. It is possible that this feature comes from the location of the [CII] emission and also from the region below the three knots. No metal lines are detected at enough signal to noise in the MUSE spectrum.

In Fig. 13b, we show the [CII] profile in velocity space together with the $Ly\alpha$ emission extracted slicing the $Ly\alpha$ image with pseudo slits oriented as in Fig. 12. These spectra are very noisy, but we can see that the highest peak of the integrated $Ly\alpha$ emission receives a contribution from the three knots, while the tail may not be originated from the northern region of LBG-1 (N position in Fig. 10). This could indicate the presence of a gas with higher velocity to the south and may suggest a merger

of the three components. Alternatively, the gas could be characterized by larger NHI. By studying the $Ly\alpha$ profile in terms of radiative transfer models could inform about the two possibilities (see next section).

## 5. $Ly\alpha$ radiative transfer model

We use the most updated version of the FLaREON radiative transfer model (Gurung-López et al. 2019b) to quantify the properties of the interstellar, circum-galactic, and intergalactic (IGM) medium that can condition the shape and intensity of the $Ly\alpha$ emission in the AzTEC-3 protocluster. FLaREON was based on LyART (Orsi et al. 2012), a radiative transfer Monte Carlo code of $Ly\alpha$ emission, and predicted the $Ly\alpha$ line profile escaping a galaxy through different outflow configurations. Several outflow geometries were implemented, such as thin shell of HI gas (see also Verhamme et al. 2006; Gronke & Dijkstra 2016), and a galactic wind (Fig. 1 in Gurung-López et al. 2019a). In the new version of FLaREON, called zELDA, Gurung-Lopez et al. (2022) developed a new thin shell model in which the intrinsic galaxy spectrum contains a continuum in addition to a (Gaussian) $Ly\alpha$ emission line. The variable set delivered by zELDA includes the systemic redshift, $z_{\text{in}}$, the outflow expansion velocity, $V_{\text{exp}}$, the HI column density, NHI, and the dust optical depth, $\tau$. Inflows are predicted such as the outflows, but with negative expansion velocities.

To be able to perform a quantitative comparison between our observed $Ly\alpha$ profiles and the ones predicted by zELDA, a Gaussian kernel representing the MUSE resolution is applied to the
<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Left panel: 1D profiles of the Lyα (blue) and [CII] (red) emissions toward LBG-1. The [CII] profile is taken from Riechers et al. (2014) and corresponds to the region outlined in Fig. 10 (red contours). Right panel: [CII] profile in red as in the left panel and Lyα profile obtained slicing the Lyα emission horizontally on top of LBG-1 (green), tilted by 110° on top of LBG-1 (cyan), tilted by ten degrees below LBG-1 (magenta), and tilted by ten degrees above LBG-1 (blue). The color coding of these spectra is like in Fig. 12. The two former spectra are the ones with higher signal to noise as can be seen by the smaller error bars. The thin black curve corresponds to the spectrum of the integrated Lyα emission as the blue curve of the left panel.</figcaption></figure>

predicted profiles and the systemic redshift is provided. Also, the observed-spectrum sampling and signal to noise are taken care of in the observed frame (see Sect. 6.1 of Gurung-López et al. 2021, for details). The quantitative comparison involves obtaining the model that produces the best fit of the observed spectrum. The best fit is obtained with a Markov chain Monte Carlo (MCMC) approach. In each step of the chain, a line profile corresponding to the variables mentioned above is computed and compared with the observed spectrum. The MCMC takes into account the uncertainties in the observed spectrum. After the MCMC run, we consider the model corresponding to the 50th percentile of the probability distribution function of the variable as the best-fit profile and the 16th (84th) percentile as the model corresponding to the lower (upper) limit of each variable.

As shown in Gurung-Lopez et al. (2022), there is a grid of parameters where the MCMC chains can look for the best combination of parameters to reproduce the observed spectrum. Expansion velocities from 0 to 1000 km s⁻¹ and HI column densities from 10¹⁷ to 10²¹.⁵ atoms cm⁻² are explored (their Sect. 2). However, the starting point of each chain is set by performing an initial optimization, which is an attempt to narrow the parameter space (their Sect. 4). The uncertainties of the best-fit zELDA parameters depend on the resolution, sampling, and signal to noise of the observed spectra (their Sect. 4). The zELDA performance is tested for a signal to noise larger than 5 at the wavelength of the main Lyα peak. The dust optical depth is poorly constrained even in the best signal-to-noise scenario, while the uncertainty on $V_{\text{exp}}$ can increase by a factor of 25% when the signal to noise changes from ten to 5 and the uncertainty on NHI is almost four times larger. It means that by fitting only the Lyα profile with zELDA, we are not able to provide a significant estimation of the dust absorption, but we can aim at a reliable estimation of $V_{\text{exp}}$ and NHI when the signal to noise of the observed Lyα spectrum is larger than ten.

According to the analysis in Gurung-Lopez et al. (2022), we can expect an anticorrelation between Lyα luminosity and NHI. This is explained considering that a gas with higher NHI would produce a larger number of scattering events and so a longer escaping path that can be associated to a larger dust attenuation and so a lower escape fraction of Lyα photons than in the case of lower NHI. However, the scattering events could also produce larger Lyα nebulae depending on the viewing angles.

We performed the zELDA fits of the spectra of the AzTEC-3+LBG-3 system (mosaic_1513), of LBG-1 (mosaic_1496), of the fixed-size aperture spectra extracted on top of the SMG Lyα emission peak (onlySMG Lyα), on top of the UV knots of the SMG (onlySMG UV knots), on top of LBG-3 (onlyLBG), on the bridge between SMG and LBG (bridgeSMGLBG), extracted from pseudo-slits oriented horizontally on LBG-1, tilted, tilted below, and tilted above LBG-1 (see Figs. 12 and 13). Also, we performed the fit for the sources mosaic_1548 and mosaic_1520 for which we have a signal-to-noise ratio larger than ten and an estimation of the physical parameters and photometric redshift from the COSMOS2015 catalog. We report the best fit values of the most reliable parameters, $V_{\text{exp}}$ and NHI, in Table 3, together with the $\chi^2$ of the best fit. We also report the ranges of the parameters contained within the 16th and 84th percentiles to show the extent of the parameter space.

## 5.1. Radiative transfer modeling of the AzTEC-3+LBG-3 system

In Fig. 14a, we show the models that best fit the observed spectrum of the AzTEC-3+LBG-3 system. The observed spectrum is shown in the vacuum framework as the zELDA models. By fixing the systemic redshift to that provided by the [CII] emission line, we estimated the MCMC model that best fits the data at wavelengths larger than the Lyα systemic wavelength (green curve in the figure). In fact, the zELDA models do not take into account the IGM absorption, but at $z \approx 5.3$ the IGM could condition the spectrum at wavelengths bluer that Lyα. Also, we tried to correct the bluer side of the spectrum for the IGM effect by using the average prescription by Madau (1995) and run an MCMC chain of the corrected spectrum of the AzTEC-3+LBG-3 system (red curve in the figure). The best fit zELDA models in these two cases provided similar combinations of parameters. However, since the IGM could affect the blue side of the Lyα emission line all the way through the maximum of its red peak, we tried also a zELDA fit only of the reddest part of the Lyα emission line, at $\lambda_{\text{vacuum}} > 7664$ Å (blue curve in the figure). The spectrum of the AzTEC-3+LBG-3 system is composed by one main peak and an extended tail at $\lambda_{\text{vacuum}} > 7668$ Å. By masking the wavelengths bluer than the maximum of the Lyα peak,
<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 14.** Left panel: observed-frame spectrum at the Lyα wavelength of the AzTEC-3+LBG-3 system in the vacuum. The observed spectrum is shown as data points with error bars and the smoothed spectrum as a black curve like in Fig. 9. The best-fit zELDA model that corresponds to the 50th percentile is shown in green and the green shaded area contains the models within the 16th and 84th percentiles obtained fitting the spectrum at wavelength larger than the systemic redshift. In blue, we show the models obtained fitting the spectrum only at the wavelength larger than the maximum. The red curve corresponds to the best fit in the wavelength range 7630–7700 Å for the AzTEC-3+LBG-3 spectrum corrected for the IGM absorption, by using the prescription of Madau (1995). Right panel: best-fit zELDA models for the spectra extracted in a 0.7" aperture located on the position of the main Lyα emission associated to AzTEC-3 (upper left), the position of the UV knots of the SMG (upper right), the position of the main Lyα emission associated to LBG-3 (lower left), and the bridge between SMG and LBG (lower right). Black dots with error bars are the observed-frame spectra as shown in Fig. 9, the black lines are the smoothed spectra, the green (blue) curves and shaded areas are the best-fit models and the models within the 16th and 84th percentile of the model parameter space obtained fitting the spectra at wavelength larger than the systemic redshift (larger than 7669 Å). Vertical blue dashed lines indicate the Lyα redshift given by the AzTEC-3 [CII] detection. The zELDA fits are performed fixing the redshift to the systemic inferred by the [CII] emission peak.</figcaption></figure>

zELDA finds the best compromise between the red peak and the extended tail of the line as the best fit, which corresponds to models with a wide range of HI column densities (see Table 3) and $V_{\text{exp}}$ up to 90 km s$^{-1}$.

In general, zELDA mainly accounts for the main peak of the Lyα emission line, which is dominant in intensity with respect to the extended tail. To identify the regions that most contribute to the extended tail at $\lambda_{\text{vacuum}} > 7668$ Å, we ran zELDA fits of the onlySMG Lyα, onlySMG UV knots, onlyLBG, and bridgeSMGLBG spectra (Fig. 14b). A main red peak and an extended tail are seen in the onlySMG spectra and we found a zELDA best fit for the main peak and one for the extended tail for them. For the onlySMG Lyα spectrum, the Lyα profile is consistent with zELDA models characterized by $V_{\text{exp}} \le 50$ km s$^{-1}$ and NHI on the order of $10^{20}$ atoms cm$^{-2}$ for the main red peak and by $V_{\text{exp}} \sim 800$ km s$^{-1}$ and NHI on the order of $3 \times 10^{20}$ atoms cm$^{-2}$ for the extended tail. The NHI value could be responsible for the scattering of Lyα photons up to 90 ckpc as mentioned in Sect. 4.1. On the contrary, in the region of the UV knots of star formation, the AzTEC-3 Lyα profile is consistent with zELDA models with $V_{\text{exp}} < 100$ km s$^{-1}$ and NHI $< 10^{19}$ atoms cm$^{-2}$.

The best fit zELDA models of the bridgeSMGLBG and onlyLBG spectra are both characterized by $V_{\text{exp}} \sim 20$ km s$^{-1}$ and NHI on the order of $10^{20}$ atoms cm$^{-2}$. Therefore, the HI column density could be responsible for the escape of Lyα photons in the distorted region above LBG-3 (N position in Fig. 7) and in the region of interaction between AzTEC-3 and LBG-3. In this region of interaction, the gas could be turned into stars, produce star formation under favorable conditions, and allow the formation of new Lyα photons as well. In Fig. 16a, we present a cartoon to better visualize the combination of parameters inferred by zELDA for the AzTEC-3+LBG-3 system.

## 5.2. Radiative transfer modeling of LBG-1

In Fig. 15a, we show the results of modeling the Lyα profile of LBG-1. We performed the fit at wavelengths larger than the systemic redshift and at $\lambda_{\text{vacuum}} > 7669$ Å to account only for the tail. In the LBG-1 spectrum, the main red peak has intensity more comparable to the extended tail than in the case of the AzTEC-3+LBG-3 spectrum, making the best fit zELDA models broader on average to account for the tail. The best fit of the main peak is consistent with a model with $V_{\text{exp}} < 30$ km s$^{-1}$ and NHI $= 5–10 \times 10^{20}$ atoms cm$^{-2}$. The extended tail is consistent with models with a large range of expansion velocities, up to 300 km s$^{-1}$, and HI column density up to $9 \times 10^{20}$ atoms cm$^{-2}$.

To investigate which regions of space and combination of parameters better reproduce the extended tail, we performed zELDA fits of the profiles obtained slicing the Lyα emission like in Fig. 12 and we show the results in Fig. 15b. Given the noisy spectra, we obtained a large uncertainty in the best fit models and their parameters. However, we can see that the best fit models of the slices containing Lyα emission coming from the N and E positions (included in the horizontal, tilted pseudo slit located on top of LBG-1, and the one above LBG-1, see Fig. 10) are consistent with HI column density up to $10^{21}$ atoms cm$^{-2}$ and $V_{\text{exp}} < 30$ km s$^{-1}$. The best fit is consistent with models of lower NHI and larger $V_{\text{exp}}$ in the other position. This could support the idea of random distortion of the gas due to the merger of the three components toward the south of LBG-1. The combination of parameters inferred by zELDA for the LBG-1 region are shown in Fig. 16b.

## 5.3. Radiative transfer modeling of the Lyα emitting sources without systemic redshift

We ran zELDA leaving the redshift as a free parameter for the Lyα emitting sources in the AzTEC-3 environment
<figure><img src="image_17.png" /></figure>

**Fig. 15.** Left panel: observed-frame spectrum at the Lyα wavelength of LBG-1. The observed spectrum in vacuum wavelengths is shown as black datapoints with error bar, the black curve is the smoothed spectrum. The green curve is the best-fit zELDA model that corresponds to percentile 50 and the green shaded area contains the models within percentiles 16 to 84. The blue curve corresponds to a zELDA fit at wavelengths larger than 7669 Å, obtained with the scope of finding a zELDA model for the extended tail at wavelengths larger than the main red peak. Right panels: zELDA fits of the spectra obtained collapsing the LBG-1 spectrum as in Fig. 13. In green, we show the best fit and the models within the 16th and 84th percentile obtained for wavelengths larger than the systemic redshift. The blue curves correspond to the best fit of the extended tail for the spectrum in the left panel. Vertical blue dashed lines indicate the Lyα redshift given by the LBG-1 [CII] detection, we fixed while performing the fits.

<figure><img src="image_18.png" /></figure>

**Fig. 16.** Qualitative representation of the combination of the zELDA parameters that best fit the Lyα emission of the AzTEC-3+LBG-3 system (left panel) and of the LBG-1 region (right panel). The size of the dots in the cartoon is proportional to $V_{exp}$ and the color scale corresponds to the NHI value in a way that fainter color indicates lower NHI. The smallest and darkest dots show $V_{exp} < 30$ km s$^{-1}$ and NHI ~ $10^{21}$ atoms cm$^{-2}$. The largest dot corresponds to $V_{exp}$ up to 800 km s$^{-1}$, the faintest color to NHI ~ $10^{19}$ atoms cm$^{-2}$.

without systemic redshift (Fig. 17). We focus on the sources with a counterpart in the COSMOS2015 catalog and with a photometric redshift consistent with $z \sim 5.3$. The best fit parameters are shown in Table 3, where we also report the best fit systemic redshift value. As expected, the best fit redshift is consistent with $z \sim 5.3$. The best fit zELDA models are consistent with galaxies with HI column densities up to $10^{20}$ atoms cm$^{-2}$ and expansion velocities on the order of a few tens of km s$^{-1}$. This combination of parameters is consistent with the fact that these sources are characterized by low masses and moderate star-formation rates (see Table 2). In particular, the mosaic_1520 source is characterized by a star-formation rate of more than $100 M_{\odot}$ yr$^{-1}$ and it is plausible that this rate is accompanied by an outflow with a velocity up to 60 km s$^{-1}$. Since the CubEx detection map shows a S/N larger than 5 from 7662 to 7675 Å, we propose an alternative zELDA solution for the spectrum of the mosaic_1548 source. In this solution, we found a combination of parameters that account for two peaks, a higher blue one at 7665 Å and a lower red one at 7672 Å. The best fit zELDA model suggests the presence of an inflow of a few tens of km s$^{-1}$ and a gas with a low column density up to $2 \times 10^{19}$ atoms cm$^{-2}$. If confirmed, this solution indicates the presence of an inflow at the outskirt of the protocluster. However, the presence of an inflow, visible as a Lyα blue

peak at $z \sim 5.3$, would imply the presence of an ionized bubble in the vicinity of the region of the production of Lyα photons and that those photons survive the IGM absorption from $z = 5$ to $z = 0$.

## 6. Discussion

In the $1.4 \times 1.4$ arcmin$^2$ region around the submillimeter galaxy AzTEC-3, we found ten Lyα emitting sources, including one extended emission toward AzTEC-3 itself and one toward the Lyman break galaxy, LBG-1. Significant [CII] detections were also observed at the positions of AzTEC-3 and LBG-1 (Riechers et al. 2014, 2020; Pavesi et al. 2019) and a significant CO molecular gas detection was detected at the positions of AzTEC-3 (Riechers et al. 2010). In the following section, we investigate the relation between Lyα, [CII], and CO emissions and we discuss the properties of the AzTEC-3 protocluster we can derive from them.

### 6.1. Dust and gas content revealed by Lyα, [CII], and CO detections in the AzTEC-3 environment

As described in Riechers et al. (2010, 2014, 2020), AzTEC-3 shows detections of the CO transitions from the Very Large
<figure><img src="image_19.png" /><figcaption>**Fig. 17.** Observed-frame spectrum at the Lyα wavelength of mosaic_1548 (**left**), mosaic_1520 (**right**). The observed spectrum in vacuum wavelengths is shown as black data points with error bars, the black curve is the smoothed spectrum. The green curve is the best-fit zELDA model of a redshifted Lyα main peak that corresponds to percentile 50 and the green shaded area contains the models within percentiles 16 to 84. The vertical blue dashed lines indicate the Lyα redshift given by the 50th percentile model. For the left hand side spectrum, we also show in magenta a zELDA fit of the two visible peaks at 7665 and 7672 Å. The systemic wavelength corresponding to this fit is at 7668.1 Å and the model implies a blue peak higher than the red peak as in the case of an inflow (see Table 2 for the best fit parameters).</figcaption></figure>

Array, the NOrthern Extended Millimeter Array, and the Plateau de Bure interferometer, and also rest-frame 1 mm continuum, [CII], and OH163 µm doublet emissions from an ALMA campaign and the COLDz survey.

The CO peak emission was observed to be offset from the HST ACS F814W image of the SMG by about 0.5" (3 pkpc). This offset could be partially due to dust obscuration in the rest-frame UV image (Riechers et al. 2010). The dust obscuration could also explain the difference between the UV- and FIR-derived star-formation rates of AzTEC-3 (Capak et al. 2011), suggesting that the regions of most intense star formation are highly dust obscured. In fact, a dust mass on the order of $3 \times 10^8 M_\odot$ was estimated from the SED fit of the optical-through-IR spectral energy distribution of the SMG (Riechers et al. 2014, 2020). Our MUSE observation shows a shift between the HST position of the SMG and the Lyα peak emission, which intensifies toward the southeast and the northeast. Lyα is not detected to the west. This could indicate that the dust on the main knots of star formation prevents the escape of Lyα photons on one side.

The molecular gas detections were in favor of a scenario of two gas components, a diffuse (1.3" = 50 ckpc), low-excitation component and a denser (1" = 40 ckpc), high-excitation component, where the diffuse gas has properties similar to those of normal high-z galaxies while the dense gas components has properties more similar to those of high-z FIR luminous quasars. However, as described in Riechers et al. (2010), the size of the CO emission was observed to be overall quite compact, 2.3 pkpc. Our MUSE observations show an even more diffuse HI gas, revealed by the Lyα emission, that extends up to 150 ckpc, embedding or blending with the HI gas of LBG-3. The zELDA modeling we performed on the Lyα profile informs about the presence of HI column densities up to $3 \times 10^{20}$ atoms cm$^{-2}$ above LBG-3 (Sect. 5).

The [CII] emission toward AzTEC-3 was also observed to be compact and its peak position to be consistent with that of CO (Riechers et al. 2014). The fact that the CO and [CII] peak emissions are spatially consistent implies that [CII] is tracing the position of the star-forming regions and the location of the molecular gas clouds rather than other gas phase regions of the SMG. The star-forming regions of AzTEC-3 would then contain some metals, unlike LBG-3 that lacks [CII] detection probably due to its low metallicity. The Lyα emission is much more extended than the star-forming regions traced by [CII], probably revealing the shape and distribution of the HI gas in the system,

elongated above LBG-3 and extending southern than the SMG. The elongated shape could be the result of the interaction between AzTEC-3 and LBG-3.

At the compact position of the [CII] emission, $F(Lyα) = (1.9 \pm 0.7) \times 10^{-18} \text{ erg s}^{-1} \text{ cm}^{-2}$, rest-frame EW(Lyα) = 6.7 ± 1.6 Å, and $L(Lyα) = (1.4 \pm 0.5) \times 10^8 L_\odot$, giving a ratio of $L(Lyα)/L([CII]) = (2.1 \pm 0.7) \times 10^{-2}$ (see also Riechers et al. 2020). Extracting the MUSE spectrum within a 0.7"-aperture on top of LBG-3, we estimate a $F(Lyα) = (3.8 \pm 0.6) \times 10^{-18} \text{ erg s}^{-1} \text{ cm}^{-2}$, a rest-frame EW(Lyα) = 9.5 ± 1.7 Å, and $L(Lyα) = (2.8 \pm 0.4) \times 10^8 L_\odot$. At this position no significant [CII] emission was detected and $L([CII]) < 0.17 \times 10^9 L_\odot$ (Riechers et al. 2014). So, on top of LBG-3, $L(Lyα)/L([CII]) > 0.6$.

We can now investigate the $L(Lyα)/L([CII])$ ratio in terms of its relation to star formation. A Lyα luminosity lower than the theoretical value and so a lower $L(Lyα)/L([CII])$ ratio can be usually explained by the presence of dust. We refer to the SFR versus $L([CII])$ empirical relation obtained by De Looze et al. (2014) to estimate the $L(Lyα)/L([CII])$ expected value for different kind of high-redshift star forming galaxies. We consider their Eq. (3) to relate the SFR and the luminosity of a SFR-sensitive emission line. Also, we assume that SFR is the total UV plus FIR star-formation rate and that it is equal to the SFR(Lyα) in the case all the emitted Lyα photons escape the galaxy (Fig. 18). We take the SFR(Lyα) expression from Kennicutt (1998), for which $SFR = 10 M_\odot \text{ yr}^{-1}$ implies an intrinsic $L(Lyα) = 10^{43} \text{ erg s}^{-1}$. By assuming the $L([CII])$ vs SFR relation of local starburst and high-z star-forming galaxies (Eqs. (7), (8), (11), (14), and (17) in De Looze et al. 2014), $L(Lyα)/L([CII]) \sim 0.3-0.7$ at a SFR(UV+FIR) of about $1200 M_\odot \text{ yr}^{-1}$ (yellow vertical line in Fig. 18). The $L(Lyα)/L([CII])$ theoretical ratio could decrease in the case of subsolar metallicity. A dust extinction corresponding to $A_V = 0.5$ could explain the $L(Lyα)/L([CII])$ ratio calculated for the spectrum extracted on top of the SMG, while the theoretical value is comparable with the ratio we calculate for the spectrum extracted on top of LBG-3.

Harikane et al. (2018) showed an anticorrelation between $L([CII])$ and rest-frame EW(Lyα) for galaxies at $5 < z < 7$ (see the right panel of Fig. 18). Based on the EW(Lyα) we calculate, either the $L([CII])$ is much stronger in AzTEC-3, implying that AzTEC-3 is a source more metal rich than the galaxies studied in that work, or Lyα is more suppressed due to the
<figure><img src="image_20.png" /><figcaption>**Fig. 18.** Left panel: total SFR vs [CII] luminosity from the empirical relations of De Looze et al. (2014). The red circle corresponds to the total SFR and $L([CII])$ of AzTEC-3. Middle panel: theoretical value of $L(Ly\alpha)/L([CII])$ versus total SFR obtained assuming that SFR(Ly$\alpha$) is equal to the total SFR and SFR(Ly$\alpha$) is calculated with the calibration from Kennicutt (1998). This calculation is intended to give an idea of the theoretical value of $L(Ly\alpha)/L([CII])$ at the location where the [CII] emission is concentrated. The vertical yellow line corresponds to the total SFR of AzTEC-3 for all the possible values of $L(Ly\alpha)/L([CII])$. Right panel: observed $L(Ly\alpha)/L([CII])$ versus rest-frame Ly$\alpha$ equivalent width as shown in Harikane et al. (2018). The average relation for all the galaxies in that work is shown as a red dashed line. The gray symbols are their data points for the sources with a Ly$\alpha$ detection, excluding the upper limits. The red and black symbols correspond to AzTEC-3 (red circle), LBG-3 (black facedown triangle), and LBG-1 (black square).</figcaption></figure>

higher dust content in AzTEC-3 (their Figs. 18–20) which is more likely the case given the nature of AzTEC-3. Also, a $\log(L([CII])/SFR_{FIR}) = 6.8$ is lower than the value they calculated for a galaxy with a rest-frame EW(Ly$\alpha$) lower than 30 Å. This is not surprising given the fact that the $L([CII])/L_{FIR}$ ratio suggested the presence of a stronger-than-average, not normal, far-UV radiation field (Riechers et al. 2014, 2020).

Riechers et al. (2014) showed that the [CII] emission is compact over the entire velocity range, but also showed a blue low-significance wing toward LBG-3, that could correspond to either a tidal feature or a [CII] outflow. A blueshift of the OH163 µm doublet also supported this hypothesis. The presence of an outflow is plausible due to the intense AzTEC-3 star-formation rate. In fact, the SFR value and the size of the star-formation region as given by the [CII] emission determine a star-formation rate per unit of area larger than the limit expected for supporting starburt-driven outflows (Heckman 2001; Riechers et al. 2020). Modeling the Ly$\alpha$ profile, we found best fit models that are consistent with the presence of outflows from the UV star-forming regions of the SMG and the presence of high NHI where the Ly$\alpha$ emission is brighter, where, therefore, HI scattering could be efficient. This suggests that the gas kinematics around AzTEC-3 are the result of its ongoing starburst, and may be further shaped by its possible interaction with LBG-3.

The only other galaxy in the AzTEC-3 neighborhood with a significant detection of [CII] is LBG-1. Riechers et al. (2014) found that the [CII] emission is spatially resolved and covers all the three optical knots. The infrared SED fitting gave a dust mass less than $9 \times 10^7 M_{\odot}$, which implies a dust extinction much lower (Capak et al. 2011) than in the case of AzTEC-3 and a $SFR_{FIR} < 54 M_{\odot}$ yr$^{-1}$, which makes LBG-1 consistent with being a main sequence galaxy at $z = 5$, composed of a young starburst and an underlying older population. At the compact position of the [CII] emission, $F(Ly\alpha) = (1.3 \pm 0.8) \times 10^{-18}$ erg s$^{-1}$ cm$^{-2}$, rest-frame EW(Ly$\alpha$) = $26 \pm 11$ Å (assuming a continuum flux density of $1 \times 10^{-20}$ erg s$^{-1}$ cm$^{-2}$ Å$^{-1}$), and $L(Ly\alpha) = (1.0 \pm 0.6) \times 10^8 L_{\odot}$, giving a ratio of $L(Ly\alpha)/L([CII]) = 0.06 \pm 0.04$ (Pavesi et al. 2019), that is consistent with the theoretical value by assuming

a dust reddening corresponding to $A_V = 0.2$ (see Table 2) and extinction law from Calzetti et al. (2000).

In addition to LBG-1, we detected eight other Ly$\alpha$ emitting galaxies in the AzTEC-3 field. They tend to be located in an elongated configuration (see Fig. 3). Some other objects with photometric redshift consistent with that of the SMG, but without a Ly$\alpha$ emission in our MUSE data, are also distributed around AzTEC-3 (Fig. 6). They could lack of Ly$\alpha$ emission due to an unfavorable amount and distribution of dust toward the line of sight, but they could also be lacking the HI gas reservoir that would shine in Ly$\alpha$. It could have been consumed by star formation or by tidal stripping within the AzTEC-3 environment. Instead, the elongated distribution of the Ly$\alpha$ emitting galaxies could follow the elongated distribution of HI gas falling toward the SMG. More sensitive data would be needed to confirm this hypothesis.

## 6.2. Evidence of merger phenomena in the AzTEC-3 environment

There is evidence in the literature that supports the idea that compact SMGs could be the result of major mergers (e.g., Tacconi et al. 2006; Swinbank et al. 2008; Engel et al. 2010; Simpson et al. 2014). In the protocluster of our study, the high molecular gas mass (roughly 5 times the stellar mass) of AzTEC-3 and its $SFR_{FIR}$ equal to $1100 M_{\odot}$ yr$^{-1}$ (Riechers et al. 2014) suggest the presence of a heavily obscured starburst, possibly triggered by a major merger. The richness of the gas reservoir in the AzTEC-3 system (both dense and diffuse molecular gas, and diffuse atomic gas revealed by Ly$\alpha$) supports the idea of a gas-rich merger.

The FIR luminosity surface density of AzTEC-3 was estimated to be $5 \times 10^{12} L_{\odot}$ kpc$^2$ (Riechers et al. 2014) and it is at the limit for a starburst to be supported by radiation pressure. The [CII] to FIR luminosity ratio ($L([CII])/L_{FIR}$) and the no-detection in X rays (Riechers et al. 2020) suggested the presence of a strong radiation field, but gives no direct indication for the presence of an obscured AGN (Riechers et al. 2014). The merger of the SMG progenitors could have triggered the starburst we
currently see. Also, the observed $Ly\alpha$ luminosity supports the idea that the $Ly\alpha$ photons are mainly produced by recombination of atoms excited by a starburst radiation field (see Sect. 6).

Based on the kinematics of our observations of $Ly\alpha$ emission, a fast outflow can depart from the starburst and extend to the south. We also detected a bridge of high column density gas between AzTEC-3 and LBG-3 that can be the result of the merger event between AzTEC-3 and LBG-3. The elongated $Ly\alpha$ emission north of LBG-3 suggests tidal interaction with the SMG gravitational potential. The gas could have been stripped and could have assumed a tidal structure as a consequence of a first passage between LBG-3 and AzTEC-3. The structure seen in the region to the south of the SMG suggests interaction between the SMG and the mosaic_199 source, a third phenomenon of interaction (Fig. 8) around AzTEC-3.

In addition to this, in the field there is LBG-1 that itself could be a merging system. The $Ly\alpha$ observations and modeling suggest the presence of a gas with expansion velocities of 50 up to $200 \text{ km s}^{-1}$ in the southern region of LBG-1 that could follow the merging rotation and an NHI overall on the order of $10^{21}$ atoms cm$^{-2}$ that could make a gas-rich merger. We do not have enough resolution to provide a more detailed description of the LBG-1 merger.

However, as proposed by Riechers et al. (2014), this system could be in an early stage of merger. The final stage could drive to a more intense star-formation rate than the one currently measured, but probably to a much fainter $Ly\alpha$ emission absorbed by the dust produced in the star formation merger-driven event (Yajima et al. 2013).

Other systems in the literature at $z \sim 6$, such as CR7 (Matthee et al. 2020b) and Himiko (Ouchi et al. 2013) share properties with LBG-1. For both sources from the literature, the $Ly\alpha$ emissions were found to be elongated along the direction connecting the multiple components and the $Ly\alpha$ peaks were seen offset with respect to the brightest UV emissions. The offset could be related to a merger event (see also Jiang et al. 2013) and the elongation of the $Ly\alpha$ emission along the axis of the multiple components could trace the underlying gas distribution (Matthee et al. 2020b), which seems distorted toward the southwest in the case of LBG-1. Riechers et al. (2014) measured a $L([CII])/L_{FIR}$ ratio of $3-9 \times 10^{-3}$, consistent with the value for normal star-forming galaxies (see also Pavesi et al. 2019), and a $L([CII])/L(\text{CO}) > 4600$ indicative of a moderate radiation field. Also, the FWHM and the EW($Ly\alpha$) of LBG-1 are consistent with the values of normal star-forming galaxies rather than AGN.

## 6.3. Study of AGN activity in AzTEC-3

Our MUSE observations do not provide strong constraints on the AGN activity in AzTEC-3. However, by comparing the $Ly\alpha$ emission we detected with observations in the literature that support AGN activity, we found differences that point more to a merger-driven starburst rather than to an influential AGN.

Borisova et al. (2016) and den Brok et al. (2020) showed a large sample of $Ly\alpha$ nebulae around typeI and typeII AGN at $z=3$, detected trough the CubEx software. Despite the surface brightness dimming, the radii of those nebulae were observed to be larger than 450 ckpc (more than 110 pkpc) and between 220 and 450 ckpc (50–100 pkpc and up to 300 pkpc with 480 pkpc for the Slug nebula), respectively, much larger than in our case. Also, the integrated $Ly\alpha$ luminosities are 1 order of magnitude brighter than the value of our entire system, even in the

cases of obscured tori. Also, their 1D spectra show $Ly\alpha$ peaks that are redshifted more than $1000 \text{ km s}^{-1}$. The FWHM and the EW of the $Ly\alpha$ emission line we measured, even for the spectrum extracted just on top of the SMG, are consistent with the values measured for star-forming galaxies rather than AGNs (Henry et al. 2015; Matthee et al. 2020b). However, the luminosity values we measured can still be dominated by the effect of the dust of the SMG. In fact, correcting the $Ly\alpha$ flux by the extinction of $A_V = 0.8$ estimated for the starburst, the $L(Ly\alpha)$ emitted at the compact position of the SMG becomes comparable to the faintest sources studied in Borisova et al. (2016).

Emissions of $Ly\alpha$, comparable in size to that of the AzTEC-3+LBG-3 system, were found around two star-forming galaxies at $z=6.5$ and $6.6$ (VR7 and CR7 Matthee et al. 2020a,b). However, their $Ly\alpha$ luminosities are more than 5 times brighter than in our system, implying that the SMG dust could allow the escape of fewer $Ly\alpha$ photons, even if the NHI conditions in the CGM could be similar and could produce an equally efficient scattering. For CR7, the AGN contribution was disfavored given the UV to IR luminosity ratio (Matthee et al. 2020b). In the case of AzTEC-3, Riechers et al. (2014) pointed out that a luminous AGN component was also disfavored by the non detection of the high excitation CO(16–15), even though not detecting this transition is not a hard constraint. Algera et al. (2021) found that, based on its radio emission, AzTEC-3 does not show a sign of AGN activity either and the MUSE spectrum does not present typical AGN features. Therefore, the data so far are not in favor of a scenario in which an AGN is powering the emission of AzTEC-3, even if a very obscured AGN, also a possible result of the major merger event, could still be present.

There are examples in the literature where other mechanisms of the production of the radiation field were proposed and that showed different properties in the $Ly\alpha$ emission in comparison to the AzTEC-3 system. Daddi et al. (2021) found evidence of gas accretion toward a massive system also containing an SMG at $z=2.91$. The phenomenon of gas accretion was supported by a blue shifted component of the $Ly\alpha$ emission line. However, the $Ly\alpha$ emission they observed peaks in an empty region located at the center of the halo potential well and it extends over 1000 ckpc (300 pkpc), much larger than what we detect.

The SSA22 protocluster is an overdensity with more than one extended $Ly\alpha$ emitting source. The overdensity was detected as an excess of LBGs (Steidel et al. 2000) and more than 200 $Ly\alpha$ emitters (Hayashino et al. 2004), more than 30 $Ly\alpha$ nebulae (Matsuda et al. 2004), and more than 50 SMGs (Umehata et al. 2014) were discovered in the entire area. Both the $Ly\alpha$ emitters and the SMGs trace the densest peak of the protocluster. The extended $Ly\alpha$ sources include two giant $Ly\alpha$ blobs with sizes larger than 100 pkpc (Steidel et al. 2000). They were discovered to be located in the intersection of the three filaments, traced by the $Ly\alpha$ emitters, that characterize the protocluster. In agreement with the theory that galaxy formation at high redshift occurs along large-scale filamentary overdense regions, the intersection of filaments could evolve in low-redshift clusters and the giant $Ly\alpha$ blobs could be revealing cluster progenitors. The $Ly\alpha$ luminosity of the two giant $Ly\alpha$ blobs on the order of $10^{44} \text{ erg s}^{-1}$ (about 20 times brighter than the $Ly\alpha$ emission associated to the AzTEC-3+LBG-3 system) and the association of one of the two to an SMG (LAB1, Chapman et al. 2001) suggested they may be progenitors of very massive galaxies near the center of a massive cluster (Matsuda et al. 2005) at $z=0$. The star-formation rate estimated for the SMG associated to LAB1 is ~1000 $M_{\odot} \text{yr}^{-1}$
and the CO emission line implies a large amount of molecular gas (Chapman et al. 2004), comparable to the values of AzTEC-3. However, cooling flow could be a possible scenario for the Lyα emission (Chapman et al. 2001) also in SSAA2, due to the lack of an optical counterpart in LAB1, lack of both an optical and a submillimeter source in the other blob, and for the large and bright Lyα emission, unlike what we observe in the AzTEC-3 system.

### 6.4. Estimation of the mass of the AzTEC-3 protocluster and its fate at lower redshift

We found that six of the Lyα emitting sources detected in the AzTEC-3 environment, including AzTEC-3, LBG-3, and LBG-1, have a counterpart in the COSMOS2015 catalog and so we have an estimation of their stellar masses (from $5 \times 10^8$ to $2 \times 10^{10} M_\odot$). The LBGs in the field studied in Capak et al. (2011) have stellar masses between $6.3 \times 10^7$ and $2.5 \times 10^9 M_\odot$. By summing all these mass values, we obtain a minimum value of the stellar mass in the protocluster of $(5.8 \pm 2.3) \times 10^{10} M_\odot$. If we assume that the COSMOS2015 sources shown in Fig. 6 also belong to the AzTEC-3 protocluster, the minimum value of the stellar mass in the protocluster is $(8.1 \pm 2.4) \times 10^{10} M_\odot$.

An estimation of the dark-matter halo mass of the protocluster can be obtained by assuming the relation between stellar mass and halo mass for the clusters studied in van der Burg et al. (2014). For a stellar mass of $5.8 \times 10^{10} - 8.1 \times 10^{10} M_\odot$, we can calculate a minimum dark-matter halo mass of $4 \times 10^{11} - 8 \times 10^{11} M_\odot$, which is consistent with the lower limit of dark-matter halo mass estimated in Capak et al. (2011). It is worth noting that the protocluster projected area can be as big as a 2 cMpc-radius circle (Capak et al. 2011) larger than the area probed by our MUSE observations. However, the stellar mass could be concentrated within 1 cMpc, since the center of mass of the protocluster is close to the location of AzTEC-3 and the two most massive sources of the field are about 600 ckpc apart. The relation from van der Burg et al. (2014) was obtained for clusters at $z=1$ and there could be evolution between the properties of $z \sim 5$ protoclusters and $z=1$ clusters. Cen & Zheng (2013) provided an approximate equation to relate the observed Lyα luminosity of a Lyα nebula and the dark matter halo mass of its central galaxies. By following their Eqs. (2) and (4), we can give an estimation of the dark matter halo mass in which the AzTEC-3+LBG-3 Lyα nebula is located. The total Lyα luminosity of the AzTEC-3+LBG-3 system is $(5 \pm 1) \times 10^{42} \text{ erg s}^{-1}$. This luminosity implies a dark matter halo mass of $(1.84 \pm 0.04) \times 10^{12} M_\odot$ and a size of the nebula on the order of $8 \text{ arcsec}^2$ rescaled to $z=5.3$, which is consistent with what we observe. Considering the simulation in Chiang et al. (2013), a protocluster with a halo mass of $10^{12} M_\odot$ could evolve into a cluster of $2 \times 10^{13} M_\odot$ by $z=2$ and into a Fornax-type cluster at $z=0$ with a typical mass of $2 \times 10^{14} M_\odot$. Chiang et al. (2013) also showed that $z > 5$ protoclusters with $10^{12} M_\odot$ halo masses span a region in the sky within an effective radius of about 5 cMpc. Therefore, we could be observing just the central, most massive part of the entire overdensity.

Submillimeter galaxies are thought to be high-$z$ progenitors of present-day ellipticals that could be located in the core of present-day galaxy clusters. Swinbank et al. (2008) showed that SMGs at $z=2$ could reside in a similar mass halo as SMGs as $z>4$, but could be characterized by twice their stellar mass. To achieve that stellar mass, major mergers are a plausible explanation (Swinbank et al. 2008; Engel et al. 2010; Simpson et al. 2014; Dudzevičiūtė et al. 2020; Stach et al. 2021) and we have evidence of them in the protocluster around AzTEC-3. Unlike

some other overdensites of galaxies at $z > 4$ from the literature, the AzTEC-3 environment contains only one submillimeter galaxy. As an example, the SPT 2349-56 protocluster contains 14 submillimeter galaxies at $z=4.3$ and it is expected to evolve into a very massive cluster of more than $10^{15} M_\odot$ at $z=0$ (Miller et al. 2018). Also, the overdensity discovered by Oteo et al. (2018), containing at least ten dusty star-forming galaxies at $z \approx 4$, is expected to evolve in a cluster as massive as the Coma cluster at $z=0$.

## 7. Conclusions

In this work, we have analyzed the MUSE data of the environment around the AzTEC-3 submillimeter galaxy. We have made use of the CubExtractor (Cantalupo et al. 2019) software to reduce and analyze the data, as shown in Sects. 2 and 3. We found ten Lyα emitting sources, including an extended emission toward AzTEC-3 and LBG-3 with a total Lyα luminosity of $(5 \pm 1) \times 10^{42} \text{ erg s}^{-1}$ and another extended emission embedding the three components of LBG-1. The sources appear distributed in an elongated configuration of about $70''$ in extent. Five of the sources have a counterpart in the COSMOS2015 catalog. The presence of the ten Lyα emitting sources around AzTEC-3 confirms that the region is an overdensity at $z \sim 5.3$ (Capak et al. 2011, and Sect. 5). We calculated the center of mass of the protocluster which is close to the AzTEC-3+LBG-3 system, the most massive of the protocluster.

We studied the spectroscopic properties of the Lyα emission in the AzTEC-3 environment in detail and we compared with the zELDA radiative transfer model (Gurung-Lopez et al. 2022). We note that Lyα photons are emitted toward the southeast of the SMG, probably due to dust obscuration to the west. The Lyα emission of the SMG is blended with that of LBG-3 and shows a bridge between the two galaxies, which could reveal the gas distribution due to their interaction. The Lyα emission is elongated to the north of LBG-3 resembling a tidal feature due to the interaction with AzTEC-3. An elongated feature is also seen toward the south of the SMG that could be produced by the interaction with the mosaic_199 source. The presence of an outflow departing from the SMG is supported by the size of the star-forming region, the star-formation rate estimated from the [CII], and FIR observations (Riechers et al. 2014, 2020). It is plausible that this outflow cleared channels for the escape of Lyα photons and that the Lyα photons are produced by recombination of HI atoms based on the energetics of the Lyα emission (Cantalupo 2017, and Sect. 4). The HI scattering can explain the Lyα luminosity we observe at more than 90 ckpc. Also, the modeling of the Lyα line did not require any contribution from an AGN.

Based on the resolution of our MUSE observation, we found that the Lyα emission comes from the three components of LBG-1. This could be related to the low dust content in this galaxy. However, the emission is not symmetrical and could indicate the interaction of the three components. The Lyα emission line profile shows a main peak and an extended tail which is consistent with models of galaxies characterized by gas in a wide range of velocities up to $200 \text{ km s}^{-1}$. Overall the Lyα emission is consistent with models of galaxies characterized by NHI on the order of $10^{21} \text{ atoms cm}^{-2}$ that could indicate that the interaction of the three components is a gas-rich merger. The emission is redshifted with respect to the [CII] systemic redshift and the Lyα emission line is broader than [CII] in the 1D spectrum. This is not unusual for high-$z$ Lyα emitters and could be related to the variety of gas kinematics and column densities that condition the escape of Lyα photons and instead do not affect the escape
of [CII] photons. The radiative transfer fit of the Lyα of two other star-forming galaxies indicate that they are galaxies with HI column densities up to $10^{20}$ atoms cm$^{-2}$ and HI expansion velocities on the order of a few tens of km s$^{-1}$.

Given the availability of CO and [CII] observations from previous campaigns (Capak et al. 2011; Riechers et al. 2010, 2014, 2020; Pavesi et al. 2018a, 2019) and the Lyα information from our MUSE dataset, we can discuss weather the environment can play a role in the stage of evolution of AzTEC-3. We have evidence of starburst-driven phenomena and interactions around AzTEC-3 that could support the idea that the gravitational potential of the SMG is accreting nearby galaxies (like LBG-3 and possibly mosaic_199) and accumulating gas. Due to the gravitational interactions, the submillimeter galaxy could increase its mass significantly by $z = 2$. With LBG-1, we could be assisting to another merging event in the AzTEC-3 protocluster, and the two galaxies could all merge by $z = 0$ (Riechers et al. 2020). We estimated that the dark matter halo mass of the protocluster is on the order of $10^{12} M_{\odot}$. This value is lower than that of other protoclusters studied in the literature that contain more than one submillimeter galaxy (e.g., Miller et al. 2018). However, a dark matter halo of $10^{12} M_{\odot}$ could evolve into a cluster of $2 \times 10^{13} M_{\odot}$ by $z = 2$ and into a Fornax-type cluster at $z = 0$ with a typical mass of $2 \times 10^{14} M_{\odot}$ (Chiang et al. 2013).

**Acknowledgements.** We thank Franz Bauer, Lodovico Coccato, Yara Jaffe, Evelyn Johnston, Sam Kim, Michael Maseda, Ian Smail, Frederic Vogt for useful discussions. L.G. is grateful to ETH for the use of computer facilities and to the CAS-CONICYT, CAS18016, fellowship. S.C. gratefully acknowledges support from Swiss National Science Foundation grant PP00P2_190092 and from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme grant agreement No 864361. H.D. acknowledges financial support from the Agencia Estatal de Investigación del Ministerio de Ciencia e Innovación (AEI-MCINN) under grant (La evolución de los cúmulos de galaxias desde el amanecer hasta el mediodía cósmico) with reference (PID2019-105776GB-I00/DOI:10.13039/501100011033) and from the ACIISI, Consejería de Economía, Conocimiento y Empleo del Gobierno de Canarias and the European Regional Development Fund (ERDF) under grant with reference PROID2020010107.

## References

* Algeria, H. S. B., van der Vlugt, D., Hodge, J. A., et al. 2020, ApJ, 903, 139

* Algeria, H. S. B., Hodge, J. A., Riechers, D., et al. 2021, ApJ, 912, 73

* Austermann, J. E., Aretxaga, I., Hughes, D. H., et al. 2009, MNRAS, 393, 1573

* Bacon, R., Accardo, M., Adjali, L., et al. 2010, in Ground-based and Airborne Instrumentation for Astronomy III, eds. I. S. McLean, S. K. Ramsay, & H. Takami, SPIE Conf. Ser., 7735, 773508

* Bacon, R., Vernet, J., Borisova, E., et al. 2014, The Messenger, 157, 13

* Bacon, R., Brinchmann, J., Richard, J., et al. 2015, A&A, 575, A75

* Barisic, I., Faisst, A. L., Capak, P. L., et al. 2017, ApJ, 845, 41

* Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393

* Blain, A. W., Smail, I., Ivison, R. J., Kneib, J. P., & Frayer, D. T. 2002, Phys. Rep., 369, 111

* Borisova, E., Cantalupo, S., Lilly, S. J., et al. 2016, ApJ, 831, 39

* Calvi, R., Dannerbauer, H., Arrabal Haro, P., et al. 2021, MNRAS, 502, 4558

* Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682

* Cantalupo, S. 2017, in Gas Accretion and Giant Lyα Nebulae, eds. A. Fox, & R. Davé, 430, 195

* Cantalupo, S., Arrigoni-Battaia, F., Prochaska, J. X., Hennawi, J. F., & Madau, P. 2014, Nature, 506, 63

* Cantalupo, S., Pezzulli, G., Lilly, S. J., et al. 2019, MNRAS, 483, 5188

* Capak, P. L., Riechers, D., Scoville, N. Z., et al. 2011, Nature, 470, 233

* Capak, P. L., Carilli, C., Jones, G., et al. 2015, Nature, 522, 455

* Carniani, S., Maiolino, R., Amorin, R., et al. 2018, MNRAS, 478, 1170

* Casey, C. M., Scoville, N. Z., Sanders, D. B., et al. 2014, ApJ, 796, 95

* Cen, R., & Zheng, Z. 2013, ApJ, 775, 112

* Chapman, S. C., Lewis, G. F., Scott, D., et al. 2001, ApJ, 548, L17

* Chapman, S. C., Smail, I., Windhorst, R., Muxlow, T., & Ivison, R. J. 2004, ApJ, 611, 732

* Chapman, S. C., Blain, A. W., Smail, I., & Ivison, R. J. 2005, ApJ, 622, 772

* Chiang, Y.-K., Overzier, R., & Gebhardt, K. 2013, ApJ, 779, 127

* Cortese, L., Gavazzi, G., Boselli, A., et al. 2006, A&A, 453, 847

* Daddi, E., Valentino, F., Rich, R. M., et al. 2021, A&A, 649, A78

* De Looze, I., Cormier, D., Lebouteiller, V., et al. 2014, A&A, 568, A62

* den Brok, J. S., Cantalupo, S., Mackenzie, R., et al. 2020, MNRAS, 495, 1874

* Dijkstra, M. 2016, in Understanding the Epoch of Cosmic Reionization: Challenges and Progress, ed. A. Mesinger, Astrophys. Space Sci. Lib., 423, 145

* Downes, A. J. B., Peacock, J. A., Savage, A., & Carrie, D. R. 1986, MNRAS, 218, 31

* Dudzevičiūtė, U., Smail, I., Swinbank, A. M., et al. 2020, MNRAS, 494, 3828

* Engel, H., Tacconi, L. J., Davies, R. I., et al. 2010, ApJ, 724, 233

* Faisst, A. L., Fudamoto, Y., Oesch, P. A., et al. 2020, MNRAS, 498, 4192

* Fardal, M. A., Katz, N., Gardner, J. P., et al. 2001, ApJ, 562, 605

* Fumagalli, M., Cantalupo, S., Dekel, A., et al. 2016, MNRAS, 462, 1978

* Fumagalli, M., Mackenzie, R., Trayford, J., et al. 2017, MNRAS, 471, 3686

* Goldsmith, P. F., Langer, W. D., Pineda, J. L., & Velusamy, T. 2012, ApJS, 203, 13

* Gronke, M., & Dijkstra, M. 2016, ApJ, 826, 14

* Guaita, L., Pompei, E., Castellano, M., et al. 2020, A&A, 640, A107

* Gurung-López, S., Orsi, Á. A., Bonoli, S., Baugh, C. M., & Lacey, C. G. 2019a, MNRAS, 486, 1882

* Gurung-López, S., Orsi, Á. A., & Bonoli, S. 2019b, MNRAS, 490, 733

* Gurung-López, S., Saito, S., Baugh, C. M., et al. 2021, MNRAS, 500, 603

* Gurung-Lopez, S., Gronke, M., Saito, S., Bonoli, S., & Orsi, A. A. 2022, MNRAS, 510, 4525

* Haiman, Z., Spaans, M., & Quataert, E. 2000, ApJ, 537, L5

* Harikane, Y., Ouchi, M., Shibuya, T., et al. 2018, ApJ, 859, 84

* Hayashino, T., Matsuda, Y., Tamura, H., et al. 2004, AJ, 128, 2073

* Heckman, T. M. 2001, in Gas and Galaxy Evolution, eds. J. E. Hibbard, M. Rupen, & J. H. van Gorkom, ASP Conf. Ser., 240, 345

* Henry, A., Scarlata, C., Martin, C. L., & Erb, D. 2015, ApJ, 809, 19

* Hernquist, L. 1989, Nature, 340, 687

* Hine, N. K., Geach, J. E., Alexander, D. M., et al. 2016, MNRAS, 455, 2363

* Hollenbach, D. J., & Tielens, A. G. G. M. 1999, Rev. Mod. Phys., 71, 173

* Inami, H., Bacon, R., Brinchmann, J., et al. 2017, A&A, 608, A2

* Israel, F. P., Bontekoe, T. R., & Kester, D. J. M. 1996, A&A, 308, 723

* Jiang, L., Egami, E., Mechtley, M., et al. 2013, ApJ, 772, 99

* Kennicutt, R. C., Jr. 1998, ARA&A, 36, 189

* Kuiper, E., Hatch, N. A., Venemans, B. P., et al. 2011, MNRAS, 417, 1088

* Kunth, D., Mas-Hesse, J. M., Terlevich, E., et al. 1998, A&A, 334, 11

* Laigle, C., McCracken, H. J., Ilbert, O., et al. 2016, ApJS, 224, 24

* Laursen, P., Sommer-Larsen, J., & Andersen, A. C. 2009, ApJ, 704, 1640

* Leclercq, F., Bacon, R., Wisotzki, L., et al. 2017, A&A, 608, A8

* Leitherer, C., Schaerer, D., Goldader, J. D., et al. 1999, ApJS, 123, 3

* Lofthouse, E. K., Fumagalli, M., Fossati, M., et al. 2020, MNRAS, 491, 2057

* Mackenzie, R., Fumagalli, M., Theuns, T., et al. 2019, MNRAS, 487, 5070

* Madau, P. 1995, ApJ, 441, 18

* Maiolino, R., Carniani, S., Fontana, A., et al. 2015, MNRAS, 452, 54

* Marino, R. A., Cantalupo, S., Lilly, S. J., et al. 2018, ApJ, 859, 53

* Matsuda, Y., Yamada, T., Hayashino, T., et al. 2004, AJ, 128, 569

* Matsuda, Y., Yamada, T., Hayashino, T., et al. 2005, ApJ, 634, L125

* Matthee, J., Sobral, D., Gronke, M., et al. 2020a, MNRAS, 492, 1778

* Matthee, J., Pezzulli, G., Mackenzie, R., et al. 2020b, MNRAS, 3043

* Miller, T. B., Chapman, S. C., Aravena, M., et al. 2018, Nature, 556, 469

* Orsi, A., Lacey, C. G., Baugh, C. M., & Infante, L. 2008, MNRAS, 391, 1589

* Orsi, A., Lacey, C. G., & Baugh, C. M. 2012, MNRAS, 425, 87

* Oteo, I., Ivison, R. J., Dunne, L., et al. 2018, ApJ, 856, 72

* Ouchi, M., Ellis, R., Ono, Y., et al. 2013, ApJ, 778, 102

* Pavesi, R., Sharon, C. E., Riechers, D. A., et al. 2018a, ApJ, 864, 49

* Pavesi, R., Riechers, D. A., Sharon, C. E., et al. 2018b, ApJ, 861, 43

* Pavesi, R., Riechers, D. A., Faisst, A. L., Stacey, G. J., & Capak, P. L. 2019, ApJ, 882, 168

* Pentericci, L., Carniani, S., Castellano, M., et al. 2016, ApJ, 829, L11

* Pineda, J. L., Langer, W. D., Velusamy, T., & Goldsmith, P. F. 2013, A&A, 554, A103

* Riechers, D. 2013, A Simultaneous Measurement of the Cold Gas, Star Formation Rate, and Stellar Mass Histories of the Universe, HST Proposal

* Riechers, D. A., Capak, P. L., Carilli, C. L., et al. 2010, ApJ, 720, L131

* Riechers, D. A., Carilli, C. L., Capak, P. L., et al. 2014, ApJ, 796, 84

* Riechers, D. A., Pavesi, R., Sharon, C. E., et al. 2019, ApJ, 872, 7

* Riechers, D. A., Hodge, J. A., Pavesi, R., et al. 2020, ApJ, 895, 81

* Romano, M., Cassata, P., Morselli, L., et al. 2021, A&A, 653, A111

* Rousselot, P., Lidman, C., Cuby, J.-G., Moreels, G., & Monnet, G. 2000, A&A, 354, 1134

* Schaerer, D., & Verhamme, A. 2008, A&A, 480, 369

* Scoville, N., Abraham, R. G., Aussel, H., et al. 2007, ApJS, 172, 38
Shapley, A. E., Steidel, C. C., Pettini, M., & Adelberger, K. L. 2003, ApJ, 588, 65

Shi, K., Huang, Y., Lee, K.-S., et al. 2019, ApJ, 879, 9

Simpson, J. M., Swinbank, A. M., Smail, I., et al. 2014, ApJ, 788, 125

Smolčić, V., Novak, M., Bondi, M., et al. 2017, A&A, 602, A1

Sobral, D., Santos, S., Matthee, J., et al. 2018, MNRAS, 476, 4725

Stacey, G. J., Geis, N., Genzel, R., et al. 1991, ApJ, 373, 423

Stach, S. M., Smail, I., Amvrosiadis, A., et al. 2021, MNRAS, 504, 172

Steidel, C. C., Adelberger, K. L., Shapley, A. E., et al. 2000, ApJ, 532, 170

Steidel, C. C., Bogosavljević, M., Shapley, A. E., et al. 2011, ApJ, 736, 160

Swinbank, A. M., Lacey, C. G., Smail, I., et al. 2008, MNRAS, 391, 420

Tacconi, L. J., Neri, R., Chapman, S. C., et al. 2006, ApJ, 640, 228

Tacconi, L. J., Neri, R., Genzel, R., et al. 2013, ApJ, 768, 74

Umehata, H., Tamura, Y., Kohno, K., et al. 2014, MNRAS, 440, 3462

van der Burg, R. F. J., Muzzin, A., Hoekstra, H., et al. 2014, A&A, 561, A79

van der Vlugt, D., Algera, H. S. B., Hodge, J. A., et al. 2021, ApJ, 907, 5

Venemans, B. P., Kurk, J. D., Miley, G. K., et al. 2002, ApJ, 569, L11

Venemans, B. P., Röttgering, H. J. A., Overzier, R. A., et al. 2004, A&A, 424, L17

Verhamme, A., Schaerer, D., & Maselli, A. 2006, A&A, 460, 397

Walter, F., Decarli, R., Carilli, C., et al. 2012, Nature, 486, 233

Weilbacher, P. M., Streicher, O., Urrutia, T., et al. 2014, in Astronomical Data Analysis Software and Systems XXIII, eds. N. Manset, & P. Forshay, ASP Conf. Ser., 485, 451

White, S. D. M., & Rees, M. J. 1978, MNRAS, 183, 341

Wisotzki, L., Bacon, R., Blaizot, J., et al. 2016, A&A, 587, A98

Yajima, H., Li, Y., & Zhu, Q. 2013, ApJ, 773, 151

Younger, J. D., Fazio, G. G., Huang, J.-S., et al. 2007, ApJ, 671, 1531

Zheng, W., Overzier, R. A., Bouwens, R. J., et al. 2006, ApJ, 640, 574

Zheng, Z.-Y., Malhotra, S., Rhoads, J. E., et al. 2016, ApJS, 226, 23
<figure><img src="image_21.png" /><figcaption>**Fig. A.4.** Same as Fig. A.1, but for the mosaic_1035 source. The two white contours correspond to $0.6 \times 10^{-18}$ and $0.8 \times 10^{-18}$ erg sec$^{-1}$ cm$^{-2}$ arcsec$^{-2}$.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Fig. A.5.** Same as Fig. A.1 for the P3_547 source. The two white contours correspond to $0.5 \times 10^{-18}$ and $1.0 \times 10^{-18}$ erg sec$^{-1}$ cm$^{-2}$ arcsec$^{-2}$.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Fig. A.6.** Same as Fig. A.1 for the P2_446 source. A source is visible in the HST F160W image in the center of the Ly$\alpha$ emission, but does not have a counterpart in the COSMOS2015 catalog. The two white contours correspond to $0.5 \times 10^{-18}$ and $0.8 \times 10^{-18}$ erg sec$^{-1}$ cm$^{-2}$ arcsec$^{-2}$.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Fig. A.7.** Same as Fig. A.1, but for the mosaic_414 source. A source is visible in the HST F160W image in the center of the Ly$\alpha$ emission, but does not have a counterpart in the COSMOS2015 catalog. The two white contours correspond to $0.5 \times 10^{-18}$ and $0.7 \times 10^{-18}$ erg sec$^{-1}$ cm$^{-2}$ arcsec$^{-2}$.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>**Fig. A.8.** Same as Fig. A.1, but for the mosaic_199 source. The two white contours correspond to $0.5 \times 10^{-18}$ and $0.7 \times 10^{-18}$ erg sec$^{-1}$ cm$^{-2}$ arcsec$^{-2}$. A source is visible in the F160W image which is COSMOS2015_848724. Some low SN emission could be blended with the AzTEC-3+LBG-3 system (yellow contours).</figcaption></figure>