# Detection of companion galaxies around hot dust-obscured hyper-luminous galaxy W0410-0913

M. Ginolfi, E. Piconcelli, L. Zappacosta, G. C. Jones, L. Pentericci, R.

Maiolino, A. Travascio, N. Menci, S. Carniani, F. Rizzo, et al.

► To cite this version:

M. Ginolfi, E. Piconcelli, L. Zappacosta, G. C. Jones, L. Pentericci, et al. Detection of companion galaxies around hot dust-obscured hyper-luminous galaxy W0410-0913. Nature Communications, 2022, 13, 10.1038/s41467-022-32297-x. insu-03839651

HAL Id: insu-03839651

https://insu.hal.science/insu-03839651

Submitted on 4 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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 Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Detection of companion galaxies around hot dust-obscured hyper-luminous galaxy W0410-0913

Received: 12 January 2022

Accepted: 22 July 2022

Published online: 05 August 2022

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

Check for updates

M. Ginolfi<sup>1</sup>✉, E. Piconcelli<sup>2</sup>, L. Zappacosta<sup>2</sup>, G. C. Jones<sup>3,4,5</sup>, L. Pentericci<sup>2</sup>, R. Maiolino<sup>4,5</sup>, A. Travascio<sup>6</sup>, N. Menci<sup>2</sup>, S. Carniani<sup>7</sup>, F. Rizzo<sup>8,9</sup>, F. Arrigoni Battaia<sup>10</sup>, S. Cantalupo<sup>6</sup>, C. De Breuck<sup>1</sup>, L. Graziani<sup>11</sup>, K. Knudsen<sup>12</sup>, P. Laursen<sup>8,9</sup>, V. Mainieri<sup>1</sup>, R. Schneider<sup>11</sup>, F. Stanley<sup>13</sup>, R. Valiante<sup>2</sup> & A. Verhamme<sup>14</sup>

The phase transition between galaxies and quasars is often identified with the rare population of hyper-luminous, hot dust-obscured galaxies. Galaxy formation models predict these systems to grow via mergers, that can deliver large amounts of gas toward their centers, induce intense bursts of star formation and feed their supermassive black holes. Here we report the detection of 24 galaxies emitting Lyman-α emission on projected physical scales of about 400 kpc around the hyper-luminous hot dust-obscured galaxy W0410-0913, at redshift $z = 3.631$, using Very Large Telescope observations. While this indicates that W0410-0913 evolves in a very dense environment, we do not find clear signs of mergers that could sustain its growth. Data suggest that if mergers occurred, as models expect, these would involve less massive satellites, with only a moderate impact on the internal interstellar medium of W0410-0913, which is sustained by a rotationally-supported fast-rotating molecular disk, as Atacama Large Millimeter Array observations suggest.

Observations carried out with the Wide-field Infrared Survey Explorer (WISE)<fnref n="1" /> recently led to the discovery of a population of hyper-luminous dust-obscured quasars with bolometric luminosities of $L_{bol} \ge 10^{13}-10^{14}L_{\odot}^{2,3}$, characterized by hot dust temperatures ($T_{dust} > 60$ K) and very red mid-infrared (IR) colors. These quasars are thought to be powered by buried active galactic nuclei (AGNs), with

accretion rates close to the Eddington limit<fnref n="4" />. There is a growing consensus that this rare population of extreme objects, often referred to as Hot Dust-Obscured Galaxies (Hot DOGs)², represents a unique and short-lived stage of galaxy evolution in which both heating and ejection of dusty interstellar medium (ISM) through winds powered by the central AGN are caught in the act⁵ (the so-

<footnote n="1">¹European Southern Observatory, Karl-Schwarzschild-Str. 2, D-85748 Garching bei München, Germany. ²INAF—Osservatorio Astronomico di Roma, Via Frascati 33, I-00040 Monte Porzio Catone, Italy. ³Department of Physics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford OX1 4RH, UK.</footnote>

<footnote n="4">⁴Cavendish Laboratory, University of Cambridge, 19 J. J. Thomson Ave., Cambridge CB3 OHE, UK. ⁵Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge CB3 OHA, UK. ⁶Department of Physics, University of Milan Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy. ⁷Scuola Normale Superiore, Piazza dei Cavalieri 7, 56126 Pisa, Italy. ⁸Cosmic Dawn Center (DAWN), Copenhagen, Denmark. ⁹Niels Bohr Institute, University of Copenhagen, Jagtvej 128, DK-2200 Copenhagen N, Denmark. ¹⁰Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str 1, D-85748 Garching bei München, Germany.</footnote>

<footnote n="11">¹¹Dipartimento di Fisica, Sapienza, Universita di Roma, Piazzale Aldo Moro 5, I-00185 Roma, Italy. ¹²Department of Space, Earth, and Environment, Chalmers University of Technology, Onsala Space Observatory, SE-439 92 Onsala, Sweden. ¹³Sorbonne Université, UPMC Université Paris 6 & CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98b boulevard Arago, 75014 Paris, France. ¹⁴Observatoire de Genève, Université de Genève, 51 Ch. des Maillettes, 1290 Versoix, Switzerland. ✉e-mail: michele.ginolfi@eso.org</footnote>
called blowout phase⁶,⁷), leading to an evolutionary transition from dusty starburst galaxies to ultraviolet (UV)-bright unobscured quasars, and finally to giant elliptical galaxies⁶.

Models predict that the growth of these massive star-bursting objects is driven by galaxy mergers⁸,⁹, which can deliver a significant amount of gas toward their centers, boosting their star formation rate (SFR) up to thousands of $M_\odot$yr⁻¹¹, and feeding the growth of their supermassive black holes (SMBH), thereby powering the AGN activity⁸. However, while cumulative number counts in wide-field surveys suggest that Hot DOGs tend to lie in dense large-scale environments¹¹,¹², direct observations of the mechanisms responsible for their rapid growth are still sparse or only restricted to scales of a few tens of kpc¹³. On the other hand, recent works¹⁴⁻¹⁶ have demonstrated the potential of spectroscopic imaging studies across the environment of luminous quasars and protocluster galaxies, finding companion galaxies out to hundreds of kpc.

In this work, we report rest-frame UV and rest-frame far-infrared (FIR) spectroscopic imaging observations of the Hot DOG WISE J041010.60-0913056.2 (hereafter W0410-0913), obtained with the integral field spectrograph Multi Unit Spectroscopic Explorer (MUSE) at the Very Large Telescope (VLT), and with the Atacama Large Millimeter/Sub-millimeter Array (ALMA), respectively. W0410-0913 was originally selected from a WISE all-sky survey¹ and, with a total FIR luminosity of $L_{FIR} \sim 2 \times 10^{14}$ L$_\odot$$^{17}$, and with large stellar$^{18}$ ($M_{star}$) and molecular gas$^{19}$ ($M_{gas}$) masses of $\ge 10^{11}$ M$_\odot$, is one of the brightest, most massive and gas-rich galaxies in the distant Universe$^{19,20}$. Extensive rest-frame FIR observations, tracing several molecular lines, have shown W0410-0913 to be located at a systemic redshift of $z=3.631^{18,21}$ (see details on the redshift determination in the subsection ALMA data in Methods), about only 1.7 Gyr after the Big Bang. Through the combined analysis of the Lyman-$\alpha$ transition of atomic hydrogen (at a rest-frame wavelength of 1215.67 Å) and the rotational transitions of carbon monoxide, these observations allow us to study both the circumgalactic environment surrounding the hyper-luminous galaxy, on scales of $\ge 100$ kpc, as well as the morphological and kinematic structure of its ISM, on scales of a few kpc, providing information on the multi-scale processes that drive its evolution.

## Results

**MUSE observations of the circumgalactic environment of W0410-0913.** We have obtained deep MUSE observations of W0410-0913 (exploiting a total telescope time of 2 h) as part of an observational campaign aimed at obtaining spatially resolved rest-frame UV spectroscopy observations of a Hot DOG and its halo gas, also called circumgalactic medium (CGM). Due to its wide field of view (FOV) and high sensitivity, MUSE is a suitable instrument to explore the surrounding environment of luminous sources in the distant Universe, and it has carried out observations of both morphology and kinematics of the gas distribution around luminous distant quasars$^{22,23}$ through direct imaging of the Lyman-$\alpha$ emission, which is mainly produced through recombination radiation following photo-ionization powered by bright sources$^{24}$.

The spectrum of W0410-0913 observed with MUSE shows highly attenuated rest-frame UV emission. We detect very broad (full widths at half maxima of about FWHM-2800 km s⁻¹) blueshifted components of $N\gamma\lambda$ 1240 and $C\gamma\lambda$ 1548,1550, both offset by more than 2000 km s⁻¹ with respect to the systemic redshift, likely tracing fast AGN-driven outflows$^{5,13,25}$. Also, we find a narrow (FWHM ~ 400 km s⁻¹) Lyman-$\alpha$ line, which, similarly to the narrow component of the optical [OIII] $\lambda$ 5007 line⁵, is blueshifted by about 1400 km s⁻¹. We find evidence for a small Lyman-$\alpha$ nebula associated to this line, with a maximum projected linear size of about 30 kpc. This is at odds with recent observations of several Lyman-$\alpha$ nebulae showing sizes up to hundreds of kpc detected with MUSE around UV-bright quasars, with typical $L_{bol}$ even lower than W0410-0913 and shorter exposure times$^{22,23}$. Our finding is rather to be

ascribed to the large amount of dust in W0410-0913, which prevents Lyman-$\alpha$ and UV photons to escape the galaxy and illuminate the CGM, suggesting that extreme physical conditions of the ISM can have a profound impact on the properties of the CGM accessible through Lyman-$\alpha$.

The lack of extended diffuse Lyman-$\alpha$ emission dominating the field at wavelengths close to the systemic redshift of W0410-0913 makes it possible to efficiently search for Lyman-$\alpha$-emitting companion galaxies in the MUSE datacube, and thus provides well-suited conditions to explore the surrounding environment of a Hot DOG on circumgalactic scales, informing on the physics that drives the rapid evolution of this rare population of objects.

Our MUSE observations (reaching a $2\sigma$ surface brightness (SB) limit of $5 \times 10^{-19}$ erg s⁻¹ cm⁻² arcsec⁻² in a 1 arcsec² aperture and in a single spectral channel) reveal a high overdensity of 24 Lyman-$\alpha$-emitting satellite galaxies associated with W0410-0913 within the MUSE FOV, in a region of about 400 projected physical kpc (diameter-scale) centered around the central massive luminous object (see Fig. 1). A compilation of their Lyman-$\alpha$ flux maps is shown in Fig. 1a, providing information on their projected spatial distribution in the field. The adopted techniques and the analysis performed to search for these galaxies (as well as a statistical assessment of the reliability of their detections) are discussed in the subsections describing the MUSE data in Methods.

Archival images from the Hubble Space Telescope (HST) obtained with the WFC3 F160W filter (probing rest-frame emission at around 0.3 µm) show that, despite the modest exposure (only 2400 s), ten among these sources (i.e., approximately 40% of the total) have rest-frame optical counterparts detected with a signal-to-noise ratio (SNR) higher than $4\sigma$, and corresponding SFRs ranging between $\ge 12-100$ M$_\odot$ yr⁻¹ (computed using standard $L_{UV/optical}$-SFR calibrations; see subsection Measure of the SFR from HST in Methods). Most of the remaining objects show only tentative individual detections at SNR of $2-3\sigma$, but from a stacking analysis, we can infer an average SFR of about 3 M$_\odot$ yr⁻¹.

MUSE data show that the companion galaxies discovered around W0410-0913 have offset line-of-sight velocities, traced by the peak of Lyman-$\alpha$ in their spectra (see Methods, where resonant scattering effects are also discussed), in the broad range of about [-6.000; +10.000] km s⁻¹ with respect to the systemic redshift of the central galaxy (see Fig. 1b, c), traced by its molecular gas. 75% of them (18 out of 24) show offsets in a smaller velocity range of about [-2.500; +1.500] km s⁻¹ (see subsection The circumgalactic environment of W0410-0913 observed with MUSE in Methods; see Fig. 1c), suggesting that these satellites are experiencing a strong attraction by W0410-0913.

To quantify the physical significance of the revealed overdensity we compare the Lyman-$\alpha$ luminosity-weighted number density of the companion galaxies associated with W0410-0913, $n_{LYE}$, with the expected value for field galaxies according to the Lyman-$\alpha$ luminosity functions measured in wide and deep surveys with MUSE$^{26,27}$ (see subsection The physical significance of the overdensity in Methods for details on the computation and on the comparison sample). We find that the observed $n_{LYE}$ is a factor of $\rho_{CGM} \sim 14$ higher than the expected value ($\langle n_{LYE} \rangle$; see Fig. 2). This is higher than other reported observations of overdensities previously discovered through Lyman-$\alpha$ in the redshift range $z \sim 2-7$$^{28-31}$, including proto-clusters around radio-galaxies$^{32,33}$ and multiple-quasar objects$^{34}$, although a direct comparison with these systems is challenged by the differences in selection criteria, cosmological volumes, and sensitivity limits. Also, we note that the measured overdensity is consistent with cosmological predictions in the proximity of massive, luminous objects like W0410-0913, at the redshift of interest (see subsection The physical significance of the overdensity in Methods).
<figure><img src="image_4.png" /><figcaption>b</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>c</figcaption></figure>

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

**Fig. 1** | Overdensity of Lyman-α emitting galaxies in the circumgalactic environment around the hyper-luminous W0410-0913, at z = 3.631. **a** We show a collage of Lyman-α flux maps showing the projected spatial distribution of 24 Lyman-α-emitting companion galaxies revealed by MUSE in the field of W0410-0913. We applied a spatial Gaussian filtering of 1 pixel (0.2") to the cutouts and we set a lower bound of $SB = 5 \times 10^{-18} \text{ erg s}^{-1} \text{ cm}^{-2} \text{ arcsec}^{-2}$ in the color-scale to improve the visualization. The Lyman-α emitters are distributed on scales of about 400 kpc around the W0410-0913 (see the black line, which corresponds to 300 kpc). The size of each cutout varies between 2" and 4", and is chosen according to the extensions of the emitting sources. See Table 1 for the definition of IDs and a summary of the measured properties of the Lyman-α emitters. **b** We color the cutouts shown in panel **a** according to the offset velocities of the companion galaxies with respect to the systemic redshift of W0410-0913. The redshift of each emitter is traced by the center wavelength obtained through a Gaussian fit of the Lyman-α spectrum. For W0410-0913, at the center, we show the velocity as traced by its Lyman-α, which is blueshifted with respect to the systemic redshift. The offset velocities range between -6.000 km s⁻¹ and 10.000 km s⁻¹. A histogram of their distribution is shown in panel **c**. Most of the companion galaxies (18 out of 24) show offsets within a smaller velocity range of about [-2.500; +1.500] km s⁻¹. These are more likely subject to the gravitational effect of the central massive object.

Altogether, our findings suggest that the circumgalactic environment surrounding W0410-0913 (probed by MUSE on projected physical scales of approximately 400 kpc) is one of the densest regions of the observed Universe. The number of gravitational encounters is expected to be enhanced in such a dense environment<sup>6,17,35,36</sup>, suggesting that galaxy mergers could be an efficient channel for the mass build-up of the massive luminous galaxy.

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

**Fig. 2** | The Hot DOG W0410-0913 evolves in an very dense environment. We quantify the significance of the overdensity of companion galaxies discovered around W0410-0913, both in the whole volume defined by the offset velocity range [-6.000; +10.000] km s⁻¹ (green triangle) and in the denser volume defined by the offset velocity range [-2.500; +1.500] km s⁻¹ (blue triangle). We compare their Lyman-α luminosity-weighted number densities with the expected values for field galaxies at comparable redshift according to the Lyman-α luminosity functions measured with MUSE in wide<sup>27</sup> (gray line and yellow shaded region, representing the ±1σ uncertainty on the fitting parameters) and deep<sup>26</sup> surveys (red points and black lines, representing the ±1σ errors). The number density of Lyman-α-emitting galaxies, $n_{\text{LYE}}$, associated with W0410-0913 is a factor of about 14 higher than the value (<n<sub>LYE</sub>>) expected for field galaxies with comparable Lyman-α luminosity. The green and blue error bars represent the range of Lyman-α luminosities covered by companion galaxies in the whole and denser volume, respectively.

However, our MUSE observations cannot firmly confirm this scenario, since we do not detect direct unequivocal signs of mergers. A tentative indication of ongoing dynamical interactions is provided by the morphological and kinematic distribution of gas along the line-of-sight as traced by the Lyman-α emission in and around the companion galaxies, that we study by using commonly adopted flux-extraction techniques that enable one to recover flux of extended sources allowing for local kinematic variations<sup>37</sup>. Specifically, we find that the three-dimensional (3D) distribution of circumgalactic gas illuminated by the Lyman-α-emitters dynamically associated with the central luminous galaxy (i.e., with an offset velocity range of about [-2.500; +1.500] km s⁻¹) is preferentially distributed over bridge-like shapes on scales up to ≥100 kpc, some of which appear to be aligned along the common axis with the central massive object (see Fig. 3). While these filamentary structures might be induced by the strong gravitational influence of W0410-0913, projection effects combined with the resonant nature of Lyman-α challenge this interpretation.

Further studies of the morphology and kinematics of the distribution of gas in the halo, tracing e.g., bright rest-frame FIR lines (like [CII]) with ALMA and rest-frame optical/near-IR lines with JWST, are needed to discern whether W0410-0913 is in the process of accreting its neighbors, similarly to what is done in highly interacting dense groups of massive galaxies in the local Universe through maps of atomic hydrogen filaments<sup>38</sup>.

Furthermore, the merging scenario is at odds with high-resolution spectroscopic imaging data of the internal gas structure in the ISM of the hyper-luminous central object.
<figure><img src="image_8.png" /></figure>

**Fig. 3** | Morphologically disturbed and filamentary distribution of gas in the surroundings of W0410-0913. In panels a, b we show a composition of the optimally extracted Lyman-α flux maps and kinematic maps of galaxies in the overdensity, respectively, obtained by computing the moment-0 and moment-1 of the voxels detected in their 3D-segmentation-masks (see subsection A

study of the 3D distribution of gas around W0410-0913 in Methods). The three-dimensional distribution of gas around galaxies in the denser dynamic volume (with offset velocities in the range [-2.500; +1.500] km s⁻¹) is suggestive of connecting filamentary structures on scales up to ≥100 kpc (see white and black scales).

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

**Fig. 4** | ALMA CO(6-5) observations of W0410-0913 show that its molecular ISM is structured in a rotationally-supported fast-rotating disk. In panel a we show the moment-0 (velocity-integrated intensity), in panel b the moment-1 (velocity field), and in panel c the moment 2 (velocity dispersion field). The black solid line in the lower-left corner of the moment-0 panel shows a 5 kpc physical scale. The solid line in the moment-1 represents the kinematic major axis, while the dashed line represents the minor axis. The black contours in the moment-0 represent the ±[2, 4, 8, 16, 32]σ levels emission, where σ = 3.3 × 10⁻² Jy beam⁻¹ km s⁻¹ is the root mean

square (rms) of the velocity-integrated CO(6-5) flux map. The ALMA beam, colored in red, is shown in the lower-left corners. We also show, as a function of radius, the profiles of rotational velocity (blue squares and lines) and velocity dispersion (red squares and lines) in panel d and the profile of their ratio (green squares and lines) in panel e. The error bars represent the standard deviation derived by the model. The dotted (dashed) black line in the lower panel indicates the average (median) value of $v_{rot}/\sigma$.

**ALMA observations of the internal ISM structure.** W0410-0913 has been extensively observed by ALMA targeting several rest-frame FIR lines in different bands. While existing ALMA observations cannot provide information on the companion galaxies, they succeed in describing several properties on scales of the ISM, complementing the

information retrieved by MUSE on the CGM environment. Studying the emission from the rotational transition of carbon monoxide CO(6-5), observed in Band-4 with a high resolution of about 0.2" (corresponding to 1.4 kpc at this redshift), we find a smooth velocity gradient (see a map of the line-of-sight velocity in Fig. 4a) with velocities ranging
between approximately ±500 km s⁻¹. We perform a kinematic model of W0410-0913 using the 3D tilted ring model fitting code 3DBarolo (see details on the procedure in Methods) and obtain models of the rotational velocity ($v_{rot}$) and velocity dispersion ($\sigma$) radial profiles across the whole CO(6-5) emitting region (Fig. 4a). The $v_{rot}$ and $\sigma$ radial profiles are shown in Fig. 4b. We find that, after a declining trend in the inner 1–2 kpc, $v_{rot}(r)$ follows a flat curve (within the uncertainties) reaching a value of about 500 km s⁻¹ at a radius of ~6 kpc. The $\sigma(r)$ profile instead peaks at the center (about 200 km s⁻¹ in the inner -1–2 kpc) and slowly declines outward, exhibiting high values of $\sigma$ ≥ 100 km s⁻¹ up to about 5–6 kpc. In the lower panel of Fig. 4b we show, as a function of the radius, the ratio between rotational velocity and velocity dispersion, $v_{rot}/\sigma$, which is typically used as an indicator of whether a system can be gravitationally supported by rotation or by turbulence³⁹. We find that, at each radius, the rotational velocity is higher than the dispersion velocity, with an average (median) value of $v_{rot}/\sigma = 3.5 \pm 1.8$ ($v_{rot}/\sigma = 3.0 \pm 0.5$). While W0410-0913 appears to be more turbulent than a few recently reported extreme dynamically cold disks at $z > 3$³⁹⁻⁴¹ (that show a surprisingly large $v_{rot}/\sigma$ ~ 10), the $v_{rot}/\sigma$ values that we estimate are consistent with state-of-art galaxy formation models⁴²,⁴³, indicating that the molecular gas emission observed in our target arises from a massive rotationally supported fast rotating disk.

## Discussion

ALMA observations of W0410-0913, zooming into the kpc-scale properties of its ISM, interpreted in the context of the CGM-scales studied with MUSE, challenge the scenario, suggested by models, that repeated galaxy encounters in the dense environment might have been efficient in building up the central galaxy. Indeed, multiple mergers in luminous systems are often shown to have a disruptive impact on the kinematics of the ISM¹⁷,³⁵,⁴⁴, which, in the case of W0410-0913, appears to be settled into a smooth rotation-dominated structure.

We speculate that a possible mechanism capable of solving this tension is that the mass growth of W0410-0913 could be sustained by minor mergers with relatively less massive companions, that do not significantly perturb the disk-like structure⁴³,⁴⁵. In this case, the accretion of discrete small satellites might have played a role in originating the observed high-velocity dispersion in W0410-0913 through gravitational instabilities. Also, minor mergers are usually invoked by models as an efficient channel of bulge formation at the core of massive disk galaxies⁴⁶. In the case of W0410-0913, the inner peak (<1–2 kpc) of the CO(6-5) $v_{rot}(r)$ and the high central concentration of flux in the intensity map could be interpreted as a tracer of a bulge-like structure, that minor mergers might have contributed to form. However, we note the existence of an alternative (but not exclusive) scenario in which the energy injected by the prominent AGN activity into the surrounding gas might contribute to the observed turbulence throughout the disk, as suggested by similar observations of turbulent ISM revealed through spectroscopic imaging of molecular gas in another distant hyper-luminous obscured quasar¹³,⁴⁷.

Altogether, the multi-wavelength and multi-scale observations presented here, obtained with some of the most powerful telescopes available (including the first observations targeting a Hot DOG with MUSE) provide key information on the cosmic assembly of a distant hyper-luminous hot-dust obscured galaxy. Data support that the mass build-up of W0410-0913 takes place in a dense circumgalactic environment, but we cannot identify robust signs of merger that could sustain the growth of the massive galaxy. Instead, we find that gravitational encounters, if any, can have only a modest impact on shaping the morpho-kinematic properties of the ISM. Our results serve as a benchmark to calibrate models of the rare and short-lived populations of Hot DOGs and progress in the understanding of their evolutionary transition towards luminous quasars and giant galaxies at cosmic noon.

## Methods

### Cosmology

Throughout the paper, we adopt a standard $\Lambda$ cold dark matter ($\Lambda$CDM) cosmology, using the following parameters compatible with Planck results⁴⁸: $\Omega_{\Lambda 0} = 0.69$, $\Omega_{m 0} = 0.31$, $\Omega_{b 0} = 0.05$, $H_0 = 67.7$ km s⁻¹ Mpc⁻¹. Based on these parameters, at the redshift of W0410-0913 ($z = 3.631$), the angular scale is $1'' = 7.39$ kpc.

### MUSE observations

MUSE observations of W0410-0913 were carried out in February 2019 and January 2021, during the ESO Observing Periods P102 (Program 0102.A-0729) and P106 (Program 106.2184.001), respectively (see Data Availability Section). The collection of data in P102 was organized in one observing block (OB) composed of four exposures of 680 s each, for a total exposure time on target of 45 m (1 h of telescope time including overheads). Seeing conditions ranged between FWHM = 0.8"–1.1", with airmass values of 1.1–1.3. In P106, the data acquisition was split into three exposures of 900 s each, with time on target (45 m) and telescope time (1 h) similar to the P102 program. Observations were carried out under good seeing conditions, FWHM = 0.8", and airmass $\le$ 1.1. In both programs the exposures were slightly dithered by $\pm 0.3"$, and rotated with respect to each other by 90 degrees, to reduce systematic noise.

We have processed the individual exposures from both datasets with standard data reduction techniques using the ESO-MUSE Pipeline Version 2.8.4⁴⁹. In detail, for each exposure, we used the MUSE_BIAS, MUSE_FLAT, and MUSE WAVECAL pipeline recipes⁴⁹ to perform bias subtraction, flat fielding, twilight and illumination correction, and wavelength calibration. In order to ensure an accurate relative astronomy, we registered the individual exposures using the position of point sources in the field. This step is needed since the spatial shifts introduced by the derotator wobble⁵⁰ can cause offsets of a few arcsecs.

We performed the next stages of data reduction using the tools provided by the CubExtractor software package³⁷ (Version 1.8), which has been specifically developed to improve the data quality for the detection of diffuse, low-SB emission in MUSE datacubes²²,²³,⁵¹,⁵². In detail, we improved the flat fielding by using the tool CubeFix³⁷, which performs an additional correction to homogenize the illumination across the field and as a function of wavelength. We performed the sky subtraction using CubeSharp³⁷, which can reach a better accuracy than the MUSE pipeline²². We then combined the corrected and sky-subtracted individual exposures from both datasets using an average $3\sigma$-clipping. The final datacube samples the instrument FOV of $1' \times 1'$ in pixels of size 0.2". Each pixel contains a spectrum covering the wavelength domain 4750–9350 Å, and individual channels have a spectral resolution of 1.25 Å.

Our combined MUSE observations, comprising a total telescope time of 2h, reach a ($2\sigma$) SB limit of $5 \times 10^{-19}$ erg s⁻¹ cm⁻² arcsec⁻² in $1$ arcsec² aperture and in a single spectral channel (1.25 Å), at wavelengths adjacent to 5632 Å, corresponding to the expected Lyman-$\alpha$ in W0410-0913 (systemic redshift of $z = 3.631$; see next Sections). For reference, the SB limit on a pseudo-narrow band (NB) image of 30 Å (a common size for filters used in past NB surveys), obtained by collapsing 24 wavelength layers close to the location of Lyman-$\alpha$, is $2.8 \times 10^{-18}$ erg s⁻¹ cm⁻² arcsec⁻² in $1$ arcsec² ($2\sigma$).

### The rest-frame UV spectrum of W0410-0913

In Fig. 5a, we show the observed rest-frame UV spectrum of W0410-0913, covering the whole MUSE wavelength domain (-1000–2000 Å rest-frame), extracted from a circular aperture with a radius of $1''$ (which maximizes the SNR) centered at the position of the luminous galaxy, based on its MUSE continuum image (obtained collapsing the full datacube). The spectrum shows a very faint rest-frame UV emission, with an average flux of $4.0(\pm 3.2) \times 10^{-19}$ erg s⁻¹ cm⁻² Å⁻¹ computed
<figure><img src="image_10.png" /></figure>

**Fig. 5** | MUSE Spectrum of W0410-0913. **a** We show a spectrum of W0410-0913, covering the full MUSE wavelength domain, extracted from a circular aperture with a 1"-radius. Vertical red lines are located at the expected position of several rest-frame UV lines, according to its redshift (as traced by the rest-frame FIR lines). A sigma-clipping procedure was performed to remove the residuals of sky-lines and make the visualization easier. Nvλ 1240 and Cmλλ 1548,1550 show very broad line-

profiles (FWHM-2800 km s⁻¹) and large blueshifts (−2000 km s⁻¹) with respect to W0410-0913, indicative of fast nuclear outflows. A narrow Lyman-α line (FWHM-400 km s⁻¹) is detected with an offset of −1460 km s⁻¹ from the systemic redshift. **b** We show a zoom of the spectrum at wavelengths closer to the expected Lyman-α line.

at the rest-frame wavelength of 1300 Å. This corresponds to a mono-chromatic luminosity of log(λL<sub>λ</sub>/[erg s<sup>-1</sup>]) = 43.8<sup>+0.3</sup><sub>-0.7</sub>, which is orders of magnitudes lower than the expected value log(λL<sub>λ</sub>/[erg s<sup>-1</sup>]) = 47.3<sup>+11</sup><sub>-11</sub>, estimated through a standard bolometric correction relation for rest-frame UV luminosities<sup>53</sup>, using L<sub>bol</sub> = 1.68 × 10<sup>14</sup> L<sub>☉</sub> for W0410-0913, as computed from multiwavelength photometry<sup>20</sup>. This faint UV emission can be explained through dust attenuation caused by the large amount of obscuring ISM in the system<sup>17-19</sup>.

At the systemic redshift of W0410-0913, z = 3.631, traced by its molecular gas emission, we do not find any Lyman-α emission. The average flux, computed within a spectral window of Δv = 1500 km s⁻¹ centered around the expected line center at 5632 Å, is consistent with zero within the uncertainties, ~3(±30) × 10⁻²⁰ erg s⁻¹ cm⁻² Å⁻¹ (see Fig. 5b), and by producing pseudo-NB images at the expected position of Lyman-α (collapsing the continuum-subtracted datacube along the wavelength axis using different spectral widths, from 5 to 30 Å) we do not find any sign of emission. However, a narrow emission line (FWHM-400(±80) km s⁻¹), emitted at the same location of W0410-0913, is detected in the spectrum with a blueshift of about -1500(±150) km s⁻¹ with respect to z = 3.631 (see Fig. 5). We interpret this line as a Lyman-α emission from W0410-0913, emitted at different velocities than the systemic redshift traced by the molecular gas. This interpretation is supported by the blueshift of about -1400 km s⁻¹ (with respect to the systemic redshift adopted in this work) of the narrow component of the non-resonant line [Oiii] λ 5007 from a rest-frame optical spectrum of W0410-0913 observed with KECK/NIRES⁵. These findings suggest that the rest-frame UV and optical narrow emission lines in W0410-0913 are emitted at different velocities than the molecular gas. An alternative, but less likely

(because of the reasons discussed above) interpretation is that the narrow rest-frame UV emission line that we find in the MUSE spectrum traces Lyman-α emission produced by a nearby companion galaxy that is part of the overdense environment (see subsection The circumgalactic environment of W0410-0913 observed with MUSE), and is located along the line of sight of W0410-0913.

We compute the spatial extension of such emission using the tool CubEx (from the CubExtractor package³⁷). CubEx is extensively used to perform automatic detection and extraction of extended sources in MUSE datacubes²²,²³, and it is based on a 3D extension of the connected-labeling-component algorithm with union finding, usually used in classical binary image analysis²². We feed CubEx with a sub-cube that we extract to have a wavelength range of Δλ-100 Å centered around the peak of Lyman-α emission (4388 Å), and from which we remove continuum sources by adopting a median-filtering approach. We also apply a spatial Gaussian filtering of 2 pixels (0.4") and a spectral smoothing with a Gaussian filter size of 1 layer (1.25 Å) to bring out extended but narrow features. Running CubEx with a SNR threshold for detection of SNR > 3 in individual datacube elements (voxels) we detect a nebula consisting of approximately 1600 connected voxels and a maximum spectral width, as defined in the 3D-segmentation-mask, of Δλ ~ 15 Å. From the optimally extracted image of the nebula, derived by applying the 3D-segmentation-mask to the datacube (see subsection A study of the 3D distribution of gas around W0410-0913, for details on the advantage of this method), we measure a spatial extension, defined as the maximum projected linear size, of about 30 kpc (see Fig. 6). This indicates that, under the reasonable assumption that such emission could be considered a Lyman-α nebula associated with the Hot DOG, it would be smaller than several nebulae
<figure><img src="image_11.png" /></figure>

**Fig. 6** | A small Lyman-α nebula, possibly associated to W0410-0913. **a** We show an optimally extracted map of the spatial distribution of the narrow and spectrally-offset Lyman-α emission emitted at the same location of the Hot DOG (see Fig. 5). We measure a spatial extension, defined as the maximum projected linear size, of about 30 kpc (see white line). The white contours represent the [3, 5, 10]σ levels of significance, where the noise is estimated by propagating the variance and taking into account the number of layers contributing to each pixel in the 3D mask. The last contour at 3σ defines the spatial projection of the 3D mask. For display purposes, we have added to the map of the projected 3D-mask one wavelength layer of the cube corresponding to the central wavelength of the nebula. **b** We show the 3σ contour of the nebula overlaid onto the HST WFC3 F160W image. Black contours are located around galaxies in the cutout that are detected in MUSE continuum (obtained collapsing the full datacube) and are defined as the regions in which 80% of the flux is contained.

observed with MUSE (usually with lower exposure times) around UV-bright quasars with lower $L_{bol}^{22,23}$, pointing toward a significant effect of a dusty ISM on the detectability of Lyman-α nebulae.

In the MUSE spectrum of W0410-0913 we also clearly detect the $N\lambda\lambda 1240$ line emission and the $C\iota\lambda 1548,1550$ doublet, both showing very broad line-profiles with FWHM-2800(±200) km s⁻¹, computed through a non-linear least-squares fit of a Gaussian profile. Similarly to other luminous quasars and other Hot DOGs¹³,⁵⁴, by computing the center of their Gaussian profiles, we measure large blueshifts for both $N\lambda$ and $C\iota$ of about 2000(±150) km s⁻¹ (see Fig. 5), corresponding to a redshift of $z = 3.600 \pm 0.005$. This value is consistent with the first published estimate of the W0410-0913 optical redshift, based on shallow observations of $C\iota$ and $HeII\lambda 1640$ emission obtained with Keck/LRIS⁵⁷. We note that the MUSE spectrum does not show the low SNR $HeII$ emission present in the Keck/LRIS, as the line falls in a noisy region of the spectrum. Furthermore, similarly to the $N\lambda$ and $C\iota$ broad lines in the MUSE spectrum, the optical [OIII] observed with KECK/NIRES (see above) has a broad component that is also very blueshifted, showing an offset velocity of 3000(±1000) km s⁻¹ with respect to the narrow component⁵⁸. Large blueshifts in the rest-frame UV and optical broad emission lines, like the ones observed in W0410-0913, are often interpreted as an indication of fast nuclear outflows⁵⁵.

## The circumgalactic environment of W0410-0913 observed with MUSE

To study the environment of W0410-0913 over an area of about 450 kpc × 450 kpc (corresponding to the MUSE 1' × 1' FOV), we search for Lyman-α-emitting companion galaxies in a continuum-subtracted sub-cube (produced through median-filtering) of about 380 Å centered around the expected Lyman-α emission from W0410-0913, corresponding to offset velocities of [−10.000; +10.000] km s⁻¹ with respect to its systemic redshift.

Our method can be summarized as follows. We inspect the whole wavelength domain by producing mini-cubes with spectral widths ranging from a minimum allowed number of 2 channels (2.5 Å; to avoid single-channel spikes of noise) to a maximum of 20 channels (about 25 Å; to allow for ~3 × typical FWHMs of Lyα-emitters at this redshift). We collapse each mini-cube along the wavelength direction to obtain pseudo-NB images that we scan through their whole spatial domain by measuring the aperture photometry of 1"-sized regions centered

around each pixel. We consider as candidate detections the regions with a significance of their measured flux of SNR > 3 (where the noise is extracted from the local background), we extract a spectrum, and we save their position (RA and DEC of the central pixel), central wavelength and spectral width. For multiple candidate detections arising from the same physical source, namely two or more outputs found at about the same location (within 0.5", corresponding to about half of the size of the PSF) and wavelength (within 4 Å), we keep the candidate that has higher SNR in the spectrum. We then perform a visual inspection of the pseudo-NB images and spectra of the candidate detections to make sure that there is no contamination by sky residuals. Also, for each candidate detection, we visually inspect the original non-continuum subtracted full datacube (covering the whole MUSE spectral domain) to check whether the candidate detections show the clear presence of multiple (at least two) emission lines at wavelengths not included in the sub-cube used for studying the environment of W0410-0913. In these cases, the sources would be line emitters at lower redshifts (emitting e.g., [OII] $\lambda$ 3727,3729, $H\beta$, [OIII] $\lambda$ 4959,5007, and $H\alpha$) and would not be relevant for our analysis. Through this analysis, we exclude from the list two sources emitting [OII] at $z = 0.495$ ($H\beta$ and the [OIII] doublet are detected at redder wavelengths) and $H\beta$ at $z = 0.167$ (with [OII], [OIII] and $H\alpha$ detected in the full spectrum). We assume that the remaining sources, showing a single emission line in the sub-cube centered around the redshift of W0410-0913 are Lyman-α-emitting galaxies in the surroundings of W0410-0913.

We find 24 Lyman-α emitters within a box with a diameter of about 400 kpc centered around W0410-0913. This result is confirmed by an independent analysis that we perform using CubEx. Although CubEx is usually exploited to extract diffuse extended emission, it can be used to mimic a line finder by decreasing the threshold of the minimum allowed number of connected voxels and by increasing the threshold of SNR for the detection of individual voxels. In detail, by adopting a SNR > 4 threshold over a lightly smoothed cube (0.2" pixels in the spatial domain, and no filtering on the wavelength direction), we find the same number of sources of the method described above, with only tiny differences (≤0.3") in the spatial location of the peak of a few objects, likely due to the different pre-processing of the cube.

In Supplementary Figs. 1 and 2, we show the spectra of these sources (extracted from a region with a radius of 1") centered around the central wavelength of their Lyman-α, that we compute through a non-linear least-squares Gaussian fit. A selection of these spectra is shown in Fig. 7a. From the modeled Lyman-α peaks we measure the offset velocities with respect to the redshift of W0410-0913, and find values in the range [−6.000; +10.000] km s⁻¹, spatially distributed in the MUSE FOV as shown in Fig. 1b. We note that radiative transfer effects are known to shift the peak of Lyman-α emission in the spectra of star-forming galaxies even by a few hundreds of km s⁻¹ compared to the systemic velocity⁵⁶, though this effect is unlikely to have a dominant impact on the suggested interpretation of the observations, since we measure offset velocities up to thousands of km s⁻¹.

We also produce pseudo-NB images of each object by collapsing the datacube across the wavelength direction within spectral ranges defined by [$\lambda_{cen} - \text{FWHM}:\lambda_{cen} + \text{FWHM}$], where the central wavelength $\lambda_{cen}$ and FWHM are defined by the spectra shown in Supplementary Figs. 1 and 2. A composition of the derived pseudo-NB images, showing the projected spatial distribution of the Lyman-α-emitting galaxies around W0410-0913, is shown in Fig. 1a.

As shown in Fig. 1c, the distribution of offset velocities is centered around W0410-0913, confirming that the massive luminous galaxy is located at the dynamic center of its dense environment and dominates its gravitational potential. Most of the companions (18 out of 24) show offset velocities in a smaller range of about [−2.500; +1.500] km s⁻¹. With the caveat of partial information provided by the measurement of
<table><caption>Table 1 | Measured properties of the Lyman-α emitters in the overdensity</caption><thead><tr><th>ID</th><th>RA<br>[deg]</th><th>DEC<br>[deg]</th><th>z<sub>Lyα</sub></th><th>FWHM<sub>Lyα</sub><br>[km s<sup>-1</sup>]</th><th>L<sub>Lyα</sub><br>log(L<sub>Lyα</sub>/L<sub>☉</sub>)</th><th>SNR</th><th>SFR<br>[M<sub>☉</sub> yr<sup>-1</sup>]</th></tr></thead><tbody><tr><td>#A</td><td>62.5433</td><td>-9.21704</td><td>3.5425 ± 0.0003</td><td>250 ± 52</td><td>42.21 ± 0.28</td><td>7.1</td><td>21 ± 3</td></tr><tr><td>#B</td><td>62.5442</td><td>-9.21871</td><td>3.5947 ± 0.0001</td><td>215 ± 15</td><td>42.38 ± 0.41</td><td>17.7</td><td>57 ± 3</td></tr><tr><td>#C</td><td>62.543</td><td>-9.22171</td><td>3.6100 ± 0.0003</td><td>292 ± 50</td><td>42.02 ± 0.32</td><td>6.1</td><td>--</td></tr><tr><td>#D</td><td>62.5437</td><td>-9.22193</td><td>3.6103 ± 0.0003</td><td>382 ± 62</td><td>42.15 ± 0.33</td><td>6.9</td><td>--</td></tr><tr><td>#E</td><td>62.5418</td><td>-9.22037</td><td>3.6133 ± 0.0003</td><td>374 ± 66</td><td>42.11 ± 0.31</td><td>7.0</td><td>--</td></tr><tr><td>#F</td><td>62.5463</td><td>-9.22299</td><td>3.6141 ± 0.0003</td><td>389 ± 50</td><td>42.24 ± 0.39</td><td>8.9</td><td>12 ± 3</td></tr><tr><td>#G</td><td>62.544</td><td>-9.2156</td><td>3.6173 ± 0.0002</td><td>450 ± 32</td><td>42.55 ± 0.51</td><td>16.6</td><td>64 ± 3</td></tr><tr><td>#H</td><td>62.5433</td><td>-9.21782</td><td>3.6181 ± 0.0004</td><td>367 ± 75</td><td>41.98 ± 0.28</td><td>5.6</td><td>--</td></tr><tr><td>#I</td><td>62.5411</td><td>-9.22176</td><td>3.6236 ± 0.0002</td><td>280 ± 34</td><td>42.13 ± 0.39</td><td>9.7</td><td>--</td></tr><tr><td>#J</td><td>62.5432</td><td>-9.2196</td><td>3.6252 ± 0.0002</td><td>232 ± 41</td><td>41.95 ± 0.31</td><td>6.4</td><td>--</td></tr><tr><td>#K</td><td>62.5434</td><td>-9.21904</td><td>3.6255 ± 0.0002</td><td>195 ± 32</td><td>41.94 ± 0.29</td><td>6.4</td><td>--</td></tr><tr><td>#L</td><td>62.542</td><td>-9.21932</td><td>3.6256 ± 0.0002</td><td>127 ± 29</td><td>41.81 ± 0.25</td><td>4.9</td><td>--</td></tr><tr><td>#M</td><td>62.5513</td><td>-9.21365</td><td>3.6360 ± 0.0009</td><td>701 ± 170</td><td>42.15 ± 0.25</td><td>4.5</td><td>--</td></tr><tr><td>#N</td><td>62.5436</td><td>-9.21437</td><td>3.6426 ± 0.0004</td><td>498 ± 69</td><td>42.20 ± 0.35</td><td>7.7</td><td>--</td></tr><tr><td>#O</td><td>62.5456</td><td>-9.21504</td><td>3.6466 ± 0.0004</td><td>528 ± 69</td><td>42.16 ± 0.35</td><td>7.5</td><td>--</td></tr><tr><td>#P</td><td>62.544</td><td>-9.21937</td><td>3.6478 ± 0.0005</td><td>489 ± 99</td><td>42.08 ± 0.20</td><td>5.6</td><td>23 ± 3</td></tr><tr><td>#Q</td><td>62.5456</td><td>-9.21382</td><td>3.6495 ± 0.0002</td><td>555 ± 45</td><td>42.48 ± 0.47</td><td>14.0</td><td>--</td></tr><tr><td>#R</td><td>62.5432</td><td>-9.21843</td><td>3.6504 ± 0.0005</td><td>647 ± 90</td><td>42.32 ± 0.36</td><td>8.1</td><td>--</td></tr><tr><td>#S</td><td>62.5464</td><td>-9.21743</td><td>3.6514 ± 0.0003</td><td>242 ± 54</td><td>41.89 ± 0.27</td><td>5.5</td><td>--</td></tr><tr><td>#T</td><td>62.539</td><td>-9.21532</td><td>3.6760 ± 0.0005</td><td>399 ± 90</td><td>41.92 ± 0.24</td><td>4.8</td><td>40 ± 3</td></tr><tr><td>#U</td><td>62.5486</td><td>-9.21315</td><td>3.7305 ± 0.0002</td><td>451 ± 40</td><td>42.44 ± 0.46</td><td>13.0</td><td>16 ± 3</td></tr><tr><td>#V</td><td>62.5446</td><td>-9.22293</td><td>3.7326 ± 0.0001</td><td>420 ± 25</td><td>42.66 ± 0.55</td><td>18.9</td><td>34 ± 3</td></tr><tr><td>#W</td><td>62.5372</td><td>-9.22532</td><td>3.7328 ± 0.0003</td><td>454 ± 55</td><td>42.36 ± 0.4</td><td>9.8</td><td>33 ± 3</td></tr><tr><td>#X</td><td>62.5449</td><td>-9.21782</td><td>3.7880 ± 0.0002</td><td>327 ± 41</td><td>42.16 ± 0.39</td><td>8.5</td><td>99 ± 3</td></tr></tbody></table>

In the first, second, and third columns we report the IDs and the coordinates (Right Ascension and Declination) associated to the companion galaxies detected through Lyman-α in the surroundings of W0410-0913. The coordinates correspond to the centroid of Lyman-α emission in the pseudo-NB images. Redshifts, FWHMs, and L<sub>Lyα</sub> (in the fourth, fifth and sixth column, respectively) are derived from the MUSE spectra as described in the Methods, and their errors represent the uncertainties provided by the Gaussian models. The SNR values in the seventh column are computed as the ratio between the integrated flux in the MUSE spectra and the noise associated to the width of the line (for each source, we considered spectral regions as large as 2 × FWHM). In the last column, the SFRs and their uncertainties are based on the HST images and are computed as discussed in the subsection Measure of the SFR from HST in Methods.

Lyman-α luminosities of galaxies in both volumes, namely $\Delta \log(L_{Lyα}/L_{☉}) = 0.86$ for the whole volume and $\Delta \log(L_{Lyα}/L_{☉}) = 0.75$ for the dense volume. We obtain Lyman-α luminosity-weighted number densities of $\Phi_{\text{whole}} = 3.5 \times 10^{-2} \text{dlog}(L_{Lyα}/L_{☉})^{-1} \text{Mpc}^{-3}$ and $\Phi_{\text{dense}} = 1.1 \times 10^{-1} \text{dlog}(L_{Lyα}/L_{☉})^{-1} \text{Mpc}^{-3}$. These values are higher than the expected values for field galaxies (see Fig. 2).

In the following we focus on the dense volume, which is more likely to contain companion galaxies gravitationally associated with W0410-0913. In detail, at the median Lyman-α luminosity of galaxies in the dense volume, $\log(L_{Lyα}/L_{☉}) = 42.2$, the MUSE-wide survey$^{27}$, which observed Lyman-α-emitters over a large area of $> 20 \text{arcmin}^2$, yields an expected $\Phi_{\text{field}} = 10^{-(2.1 \times 10^{-3})} \text{dlog}(L_{Lyα}/L_{☉})^{-1} \text{Mpc}^{-3}$, where the uncertainty represents the errors on the best-fitting Schechter function parameters$^{27}$. Comparing this value with $\Phi_{\text{dense}}$, we find that the circumgalactic environment around W0410-0913 has a denser concentration of Lyman-α-emitting galaxies than the average volume, with an estimated overdensity factor of $\rho_{\text{CGM}} = \frac{\Phi_{\text{dense}}}{\Phi_{\text{field}}} = 14^{+16}_{-8}$. We note that this calculation is conservative since it is based on the assumption that offset velocities with respect to the central object are tracing relative distances, while it is more likely that the observed dispersion of offset velocities (at least for objects in the dense volume) is affected by gravitational attraction (see Article and previous subsections). Thus the effective values of $\Phi$ in our field are likely to be higher, e.g., up to a factor of about 20 times higher in the dense volume case under the extreme assumption that the 18 Lyman-α emitters are contained in a cube with a side of 60" centered around W0410-0913.

Despite its high value, we find that the measured overdensity of Lyman-α emitters in the dense volume, $\rho_{\text{CGM}} = 14^{+16}_{-8}$, is consistent with

cosmological predictions in the proximity of extreme, luminous objects like W0410-0913.

To check this, we first estimate how rare the fluctuation of matter density associated with $\rho_{\text{CGM}}$ would be. The first step is to consider that in general the spatial clustering of observable galaxies does not precisely mirror the clustering of the bulk of matter in the Universe. This caveat is often addressed through the definition of a galaxy clustering bias, which describes the relationship between the spatial distribution of galaxies and the underlying dark matter density field$^{57}$. Thus, by using a standard assumption of linear bias$^{58}$, we convert the observed overdensity of Lyman-α emitters $\rho_{\text{CGM}}$ into an estimate of the underlying matter overdensity, $\rho_M$, by dividing it by the clustering bias parameter of Lyman-α emitters, $b$. By using an updated value of $b = 2.8 \pm 0.4$, recently computed through a 3D clustering analysis of hundreds of Lyman-α emitters at $z \lesssim 4$ in the large MUSE-Wide survey$^{59}$, and consistent with previous similar studies$^{58,60}$, we obtain an estimate of the matter overdensity in the surroundings of W0410-0913 of $\rho_M = 5^{+7}_{-3}$. To estimate the rarity of such a fluctuation, we start from the computation of the linear rms density fluctuation$^{61}$, $\sigma_{\text{LN}}$, predicted by the standard $\Lambda$CDM cosmology. For the dense volume over which we compute $\rho_M$, i.e., about 210 Mpc$^3$, we find a value of $\sigma_{\text{LN}} \sim 1.6$. The next step is to compute the probability associated with observing an overdensity $\rho_M$ given a certain density fluctuation. On the physical scales of our system, the density field is non-linear, and a possible approximation for the true density distribution $P(\rho)$ is a log-normal model$^{62}$, for which the $\sigma_{\text{LN}} \sim 1.6$ can be interpreted as the rms in $(\rho)$, with a mean in this quantity of $-\frac{\sigma_{\text{LN}}}{2}$. Under this assumption, a matter overdensity of $\rho_M = 5^{+7}_{-3}$ corresponds to a probability of $P(\rho_M) = 10^{-1.4 \pm 0.6}$. Given the considered volume of 210 Mpc$^3$ and the estimated
$P(\rho_M)$, we find that the comoving density of objects like the central galaxy would need to be about $n_{halo} - \frac{P(\rho_M)}{210 \text{ Mpc}^3} \sim 10^{-3.7 \pm 0.6} \text{ Mpc}^{-3}$. This values is consistent, within the uncertainties, to the predicted comoving density of halos hosting central massive galaxies like W0410-0913, with maximum velocities of $v_{max} = 500 \pm 70 \text{ km s}^{-1}$ (see the radial profile of the rotational velocity in Fig. 4b) at $z = 3.631$, which we find to be $n_{halo} \sim 10^{-3.8 \pm 0.5} \text{ Mpc}^{-3}$, as we compute integrating the halo mass function provided by the large N-body simulation MDPL2$^{63}$, included in the open-access CosmoSim database (https://www.cosmosim.org/cms/simulations/mdl2/).

## Measure of the SFR from HST

We estimate the SFR of the Lyman-$\alpha$-emitting companion galaxies around W0410-0913 using archival HST data, that we retrieve from the public Barbara A. Mikulski Archive for Space Telescopes (MAST) portal (https://archive.stsci.edu/). A 2' x 2' field around our object was observed for a total exposure time of 2400 s in October 2012 during the HST Cycle 20 (Program ID: 12930) exploiting the filter F160W on the Wide Field Camera 3 (WFC3), centered at about 15436 Å, which corresponds to -0.3 µm at the redshift of interest.

We extract the fluxes in the HST WFC3/F160W image from circular apertures with 1"-sized diameters centered around the peak emission, that we search within a radius of 0.5" from the Lyman-$\alpha$ peak (obtained through a 2D Gaussian fit of the pseudo-NB images; see Fig. 1a), to allow for slight spatial offsets between Lyman-$\alpha$ and the rest-frame UV/optical continuum$^{64}$. We then derive the luminosity density $L_v$ and convert it in SFR by using a standard rest-frame UV/optical $L_v$-SFR calibration and assuming a Chabrier initial mass function$^{65}$. Although the available HST data are shallow (with a $3\sigma$ flux limit corresponding to a $SFR_{3\sigma} \sim 9\text{M}_\odot \text{yr}^{-1}$), we detect 10 galaxies at SNR > 4, with star formation rates ranging between SFR ~ 12–100$\text{M}_\odot \text{yr}^{-1}$.

In Supplementary Fig. 3, we show the HST cutouts for the companion galaxies detected with the WFC3/F160W filter, along with the Lyman-$\alpha$ contours from MUSE. A selection of HST maps is shown in Fig. 7b. The SFRs are reported in Table 1.

To obtain a better constraint of the average star formation rate in the 14 companion galaxies not detected by HST we perform a median stacking of their 5" x 5" cutouts centered at the Lyman-$\alpha$ peaks. Prior to stacking, we run a sigma-clipping algorithm (with a $3\sigma$ threshold on individual pixels) over each cutout to mask possible continuum sources outside the central circular area of 1" diameter. This is necessary to guarantee an uncontaminated estimate of the noise level in the final image, that would be otherwise affected by bright field objects. The noise level in the stacked image decreases down to a corresponding $3\sigma$ SFR limit of $SFR_{3\sigma} \sim 2.5\text{M}_\odot \text{yr}^{-1}$, enabling us to measure with a SNR ~ 3.6 an average SFR ~ 3$\text{M}_\odot \text{yr}^{-1}$, suggesting that a non-negligible star formation is taking place in the companion galaxies not detected by HST too.

## A study of the 3D distribution of gas around W0410-0913

Extended diffuse Lyman-$\alpha$ emission, mainly produced through recombination radiation following photoionization (often referred to as fluorescence) powered by UV-bright sources$^{24}$, is a commonly used tracer to observe the cool gas (with temperatures of approximately $T = 10^4$ K) in the CGM around galaxies. At $z \ge 2-3$, where the Lyman-$\alpha$ is redshifted into the optical domain, great advances have been recently made with wide-field integral field spectrographs at large telescopes (like MUSE, or the Keck Wide-Field Imager), which have identified ubiquitous Lyman-$\alpha$ nebulae around quasars, extending up to hundreds of kpc$^{22,23}$, and Lyman-$\alpha$ halos around star-forming galaxies, extending up to tens of kpc, between factors of 5–20 larger than the rest-frame UV sizes$^{66}$.

We search for extended Lyman-$\alpha$ emission around the companion galaxies of W0410-0913 to study the morphological and kinematic

distribution of gas in the circumgalactic environment of W0410-0913. We run CubEx over the sub-cube in which the 24 Lyman-$\alpha$-emitters are found, using the same pre-processing and setup of parameters adopted for analyzing the Lyman-$\alpha$ nebula (see previous subsections), namely a SNR > 3 threshold for voxel-detection in a smoothed cube with a Gaussian spatial filter of 0.4" and a Gaussian spectral filter of 1.25 Å. This setup is different from that exploited in the line-emitters searches, and it is tailored at extracting faint and narrow features. We produce optimally extracted Lyman-$\alpha$ flux maps of the sources by selecting the voxels detected in their 3D-segmentation-masks and integrating their fluxes along the wavelength axis. This technique is operationally equivalent to combining different pseudo-NB images in which the size of the spectral filter is optimized for each pixel. The advantage of this approach, enabled by integral field spectrographs, is that one can maximize the SNR of the whole emitting region by revealing both kinematically narrow and broad features of a source at the same time. In contrast, imaging through classical NB filters with a fixed width would produce either a noise-dilution or an underestimation of the signal emitted.

A collection of the optimally extracted images of the companion galaxies is shown in Fig. 3a. We find that some of the luminous companion galaxies in the denser dynamic volume around W0410-0913, with offset velocities in the range [-2.500; +1.500] km s$^{-1}$, show extended diffuse Lyman-$\alpha$ envelopes (not revealed, or barely detected, in the pseudo-NB images; see Fig. 1a) that are suggestive of connecting filamentary structures on scales up to ≥100 kpc (see Fig. 3). We use the 3D-segmentation-masks provided by CubEx to produce also the 2D maps of the first moment (moment-1) of the flux distribution in its spectral domain. These maps give an indication of the centroid offset velocities with respect to the systemic redshift of W0410-0913 at each spatial location. A collage of such kinematics maps is showed in Fig. 3b.

## ALMA data

In this work, we use archival ALMA observations of W0410-0913 targeting the emission-line spectra arising from rotational transitions of carbon monoxide CO(4-3), CO(6-5), and CO(7-6), and the atomic carbon transition [CI]$^3P_2-^3P_1$ (hereafter [CI]). CO(4-3) 461GHz, redshifted to 99.5GHz at $z=3.631$, was observed with ALMA in Band 3 on January 2018 (Program ID: 2017.1.00123.S), with a total integration time of 11.5 min. Deep observations of CO(6-5) 691.4GHz, redshifted to 149.3GHz, were taken in Band 4 in December 2017 (Program ID: 2017.1.00908.S) exploiting a total integration time of about 3.6h. CO(7-6) 806.6GHz and the adjacent [CI] 809.3GHz, redshifted to 174.2GHz and 174.75GHz respectively, were observed in Band 5 on September 2018 (Program ID: 2017.1.00123.S), with a total integration time of about 10 min. See the Data Availability Section for information on the data repositories. We performed a standard calibration of the data using the Common Astronomy Software Applications package (CASA$^{67}$), running different versions of the software per each dataset, according to the Quality Assurance report provided by the ALMA staff. The pipeline calibration outputs did not show any issue, and thus we did not perform any additional flagging. We obtained continuum-subtracted visibilities using the CASA task uvcontsub, and we generated the continuum-subtracted datacubes of the spectral windows containing the emission lines by running the CASA cleaning algorithm clean with 500 iterations and a SNR threshold of 3. In order to maximize the sensitivity, we adopted a natural weighting of the visibilities for all the three datasets. The resulting synthesized beam in the CO(4-3) datacube is 0.7" x 0.5" (position angle, PA = 83°), 0.22" x 0.18" (PA = -69°) for CO(6-5), and 0.6" x 0.4" (PA = 80°) in the datacube containing CO(7-6) and [CI]. We adopted pixel sizes of 0.1", 0.03", and 0.1", and spectral bins of 50 km s$^{-1}$, 60 km s$^{-1}$, and 50 km s$^{-1}$ (which enable us to have both a good sampling of the lines and enough SNR per individual channel) for CO(4-3), CO(6-5) and CO(7-6) + [CI]
<figure><img src="image_12.png" /><figcaption>b</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>c</figcaption></figure>

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

**Fig. 8 | ALMA FIR lines observations of the ISM in W0410-0913.** a CO(4-3), b CO(6-5), and c CO(7-6) (yellow bars) and [CI] (green bars) spectra of W0410-0913 observed with ALMA. Gaussian models of the line profiles are shown with blue lines.

respectively. The average sensitivities reached in the cubes (in spectral regions close to the line frequencies) are 0.7 mJy beam⁻¹ in a spectral channel of 50 km s⁻¹ for CO(4-3), 0.08 mJy beam⁻¹ in a spectral channel of 60 km s⁻¹ for CO(6-5) and 0.8 mJy beam⁻¹ in a spectral channel of 50 km s⁻¹ for CO(7-6) + [CI].

In Fig. 8 we show the CO(4-3), CO(6-5), CO(7-6) and [CI] spectra of W0410-0913. These spectra are extracted from a spatial region defined by the connected pixels with SNR > 2 in the moment-0 maps, that we first obtain by collapsing the cube spectral channels in the range [f_cen-FWHM: f_cen + FWHM], where f_cen (the central frequency) and FWHM are derived by a Gaussian fit of a spectrum extracted from a 1"-sized circular aperture centered at the expected position of the object. Such a 2-iterations procedure helps in case of resolved observations to ensure that the whole line-emitting region is accounted for in the spectrum while avoiding signal dilution. Through a non-linear least-squares Gaussian fit of the CO(6-5) line, which has the highest SNR among the available data, we derive a systemic redshift of z<sub>CO(6-5)</sub> = 3.631±0.001. This value is consistent with an independent

measurement of z<sub>CO(6-5)</sub> recently reported<sup>18</sup> (z = 3.6301), and the tiny difference is probably due to a slightly different definition of the extracting aperture. The velocity axes of the spectra in Fig. 8 are aligned to the common systemic redshift of z<sub>CO(6-5)</sub> = 3.631. From our Gaussian fit we derive large FWHMs of about 800 km s<sup>-1</sup> (with average associated uncertainties of 70 km s<sup>-1</sup>) in CO(4-3), CO(6-5) and CO(7-6). [CI] appears to be narrower (FWHM-600 ± 90km s<sup>-1</sup>), although this might be an effect of spectral blending with the more luminous CO(7-6). Broad molecular line profiles in W0410-0913 were also reported in previously published ALMA observations of CO(4-3)<sup>19</sup>, and [CI]<sup>18</sup>, although in both cases the instrument spectral setup was tuned on the optical redshift<sup>5</sup>, which is blueshifted with respect to the systemic redshift as traced by the molecular gas, and the lines, close to the edge of the spectral window, were only partially detected. Such broad emission lines in the integrated spectra extracted from the whole galaxy reflect the rapidly spinning nature of its disk (see subsection Kinematic modeling of the CO(6-5) emission in the ISM of W0410-0913).

## Study of the CGM environment with ALMA

Unfortunately, none of these datasets, designed to observe the cold gas emission from the hyper-luminous W0410-0913, provides information about the rest-frame FIR line properties of the companion Lyman-α emitters that we discover with MUSE in its CGM. However this is not surprising, and it is mainly due to a combination of insufficient sensitivity and too high resolution (that dilutes the flux over several beams), as we test by running simulations of the estimated telescope time on the ALMA Observing Tool (OT; version released for Cycle 8). In detail, we first estimate the total molecular gas mass (M<sub>mol</sub>) from the SFR measured with HST by using standard M<sub>mol</sub>SFR scaling relations calibrated at high-z<sup>68</sup>. We then convert M<sub>mol</sub> into total CO fluxes by adopting a conservative Milky-Way-like conversion factor between M<sub>mol</sub> and the CO(1-0) luminosity, α<sub>CO</sub><sup>MW</sup> = 4.3, and standard Milky-Way-like intensity ratios between CO(1-0) and the rotational transitions observed by ALMA<sup>69</sup>. We note that the adopted value of α<sub>CO</sub> has been calibrated at solar metallicities, and should be considered conservative in our case since α<sub>CO</sub> typically increases at lower metallicities<sup>70</sup>, and high-z galaxies tend to be less metal-rich at a fixed stellar mass<sup>71</sup>. Assuming projected sizes of about 1" and line-widths of 300 km s<sup>-1</sup>, we find that the estimated total ALMA time to observe the integrated CO spectrum in a galaxy with SFR = 30 M<sub>☉</sub> yr<sup>-1</sup> (the median SFR of the HST-detected companions), with at least 3σ per channel in at least 3 spectral bins, at the same angular resolution of the available observations, would be about 1.6h for CO(4-3) and CO(7-6), and about 24.5h for CO(6-5), between a factor of 6–8 longer than the observed datasets. Other additional limitations derive from the ALMA smaller fields of view (the primary beams range from about 30" in Band 5 to about 45" in Band 3), and the narrower spectral coverage (sidebands ranging from 5000 km s<sup>-1</sup> in Band 5 to 9000 km s<sup>-1</sup> in Band 3) with respect to MUSE data. The combination of these effects reduces the number of companion galaxies that could be possibly detected (up to 60%, in the best case) and, together with the sensitivity issue discussed above, limits the information that can be derived through stacking analyses. Indeed we find signals that are consistent with noise in the stacked images that we produce, for each dataset, by combining moment-0 maps centered at the expected positions of the galaxies (using the peak of Lyman-α emission as a prior), obtained by collapsing three different velocity ranges, 300, 400, 500 km s<sup>-1</sup>, centered around the expected systemic frequencies (using the redshifts derived from the Lyman-α as a prior).

## Kinematic modeling of the CO(6-5) emission in the ISM of W0410-0913

The deep and highly resolved CO(6-5) observations of W0410-0913 are best suited to study the molecular gas kinematics in W0410-0913, that
<figure><img src="image_15.png" /><figcaption>**Fig. 9** | Moment maps of ALMA CO(6-5) observations compared with the 3D**Barolo** models. In panel a, b, c, we show the moment-0, moment-1, and moment 2 of the CO(6-5) data. For each of these panels, the three columns denote (from left to right) the observed data, models (created with 3D**Barolo**), and the corresponding residuals. The black solid line in the lower-left corner of the moment-0 panel shows a 5 kpc physical scale. The solid lines in the second row represent the kinematic</figcaption></figure>

major axis, while the dashed lines represent the minor axis. The black contours in the moment-0 maps are defined as in Fig. 4. In panel d we show the major axis (left) and the minor axis (right) PVDs. For each PVD, the data is shown by the background color, while the overlaid contours represent the model at 20%, 50%, and 80% of its maximum value. The ALMA beam, colored in red, is shown in the lower-left corners.

we model using the 3D tilted ring model fitting code 3DBarolo<sup>72</sup> (3D-Based Analysis of Rotating Objects from Line Observations). This technique has the advantage that it directly models the 3D cube taking beam-smearing effects into account and it has been extensively applied to observations of high-redshift (z> 2) galaxies<sup>73</sup>. The fit routine is performed on the continuum-subtracted CO(6-5) datacube and can be summarized as follows.

First, we isolate the signal by producing a 3D signal mask of the emission line through the 3D**Barolo** SEARCH tool (based on DUCHAMP<sup>74</sup>). In detail, this algorithm identifies pixels in the datacube with intensities above a certain SNR (SNR<sub>upper</sub> = 3.5, in our case), and then it searches around these signal peaks for emission above a second SNR threshold (i.e., SNR<sub>lower</sub> = 2.0, in our case). The 3D**Barolo** best-fit

results are highly dependent on the overall geometry of the model (i.e., number and width of the rings) and on the initial estimates of the parameters. Thus, once the signal is isolated we create the moment maps needed for the initial guesses. Using the CASA toolkit task *image.moments* we collapse the cube over the spectral channels containing the line and produce a moment-0 map, which we then fit with a 2D Gaussian (using the CASA task *image.fitcomponents*) to obtain the FWHM of the major and minor axes, the position angle and the central position. We also produce a velocity dispersion map (i.e., moment-2) including only the pixels that belong to the 3D mask identified before. Then, we fit a tilted ring model to the line emission by using the 3D**Barolo** key function 3DFIT. We adopt a value of twice the FWHM of the major axis of the restoring beam as a maximum model
radius. As the value for the width of each ring, we adopt the FWHM of the minor axis, in order to not over-resolve the data<sup>73</sup>. This corresponds to a number of 12 model rings. The model central position is fixed to be coincident with the centroid of a 2D Gaussian fit of the moment-0 map, and we use the maximum value in the moment-2 map as the initial guess of the velocity dispersion. A thin disk, with a scale height of 0.01", is assumed. We adopt an initial guess for the inclination of 45°, and we allow it to vary between 10 and 80°. Starting from this set of guesses, <sup>3D</sup>Barolo produces a model of the first ring by randomly populating the volume with discrete emitting gas clouds in order to reproduce the position angle, the rotational velocity, etc of the initial estimates. The resulting model datacube is then convolved with the ALMA synthesized beam and rescaled so that the velocity-integrated intensity of each spaxel in the data and model cubes are equal. Then, the values of the pixels in the data and model datacubes are compared, and each parameter is varied until the absolute residual (i.e., [model-observation]) is minimized. When this process reaches a minimum, the next ring is analyzed. During the first run <sup>3DFIT</sup> provides best-fit results for the rotational velocity, velocity dispersion, systemic velocity, inclination, and position angle. The three latter variables are averaged across all the rings, and the whole fit routine is run again, with only the rotational velocity and velocity dispersion as variables. Among the final <sup>3D</sup>Barolo outputs, we obtain morphological parameters (e.g., central position), and kinematic parameters such as inclination, position angle, velocity dispersion, and rotation.

In Fig. 9 we show that our model is able to reproduce the kinematic properties of the target to an excellent level, with small residuals in the moment maps and a good agreement between the model and data position-velocity diagram (PVDs). The galaxy has a well-defined rotation, with rotational velocities that range between $v_{rot} \sim 500$ km s⁻¹, as can be seen from the moment-1 map, showing the spatial distribution of line-of-sight velocities, and from the PVDs extracted along the major and minor axes (see Fig. 9). Despite the fast coherent rotation, the velocity dispersion is still large, in analogy with other recent ALMA observations of Hot DOGs reporting turbulent ISM<sup>18</sup>. Nevertheless, as discussed in the Article, we find that $v_{rot}/\sigma > 1$ at each ring in our model (see Fig. 4). We find morphological and kinematic position angles of 133 ± 5° and 133.8 ± 4.9° respectively, and morphological and kinematic inclination on the plane of the sky of 49.5 ± 3.3° and 50.1 ± 13.1°.

## Data availability

The MUSE data used in this work (programs: 0102.A-0729, 106.2184-001) are available at http://archive.eso.org/scienceportal/home. The archival HST (program ID: 12930) and ALMA (programs: #2017.1.00123.S, #2017.1.00908.S) data are available at https://archive.stsci.edu/ and http://almascience.eso.org/aq/respectively. All data generated and/or analyzed during the current study are available from the corresponding author upon request.

## Code availability

The MUSE pipeline recipes, used to reduce the MUSE data, are publicly available and the software distribution can be obtained at https://www.eso.org/sci/software/pipelines/muse/. The software CASA, used to process the ALMA data, is publicly available and can be downloaded at https://almascience.eso.org/processing/science-pipeline. The code used to model the molecular gas kinematics, <sup>3D</sup>Barolo, is publicly accessible through github (https://editeodoro.github.io/Bbarolo/), with a full documentation hosted at https://bbarolo.readthedocs.io/en/latest/. The package used for the 3D extraction of faint diffuse emission in the MUSE datacube, CubExtractor, is not publicly available, though the reader interested in using this code can contact Sebastiano Cantalupo at the mail address sebastiano.cantalupo@unimib.it.

## References

1. Wright, E. L. et al. The Wide-field Infrared Survey Explorer (WISE): Mission description and initial on-orbit performance. Astron. J. 140, 1868–1881 (2010).

2. Wu, J. et al. Submillimeter follow-up of WISE-selected hyperluminous galaxies. Astrophys. J. 756, 96 (2012).

3. Eisenhardt, P. R. M. et al. The first hyper-luminous infrared galaxy discovered by WISE. Astrophys. J. 755, 173 (2012).

4. Assef, R. J. et al. Half of the most luminous quasars may be obscured: investigating the nature of WISE-selected hot dust-obscured galaxies. Astrophysical J. 804, 27 (2015).

5. Finnerty, L. et al. Fast outflows in hot dust-obscured galaxies detected with Keck/NIRES. Astrophys. J. 905, 16 (2020).

6. Hopkins, P. F., Hernquist, L., Cox, T. J. & Kereš, D. A Cosmological framework for the co-evolution of quasars, supermassive black holes, and elliptical galaxies. I. Galaxy mergers and quasar activity. Astrophys. J. Suppl. Ser. 175, 356–389 (2008).

7. Bridge, C. R. et al. A new population of High-z, dusty Lyα emitters and blobs discovered by WISE: feedback caught in the act? Astrophys. J. 769, 91 (2013).

8. Hopkins, P. F. et al. The relation between quasar and merging galaxy luminosity functions and the merger-driven star formation history of the universe. Astrophys. J. 652, 864–888 (2006).

9. Di Matteo, P. et al. On the frequency, intensity, and duration of starburst episodes triggered by galaxy interactions and mergers. Astron. Astrophys. 492, 31–49 (2008).

10. Frey, S., Paragi, Z., Gabányi, K. É. & An, T. Four hot DOGs in the microwave. Mon. Not. R. Astron. Soc. 455, 2058–2065 (2016).

11. Fan, L., Jones, S. F., Han, Y. & Knudsen, K. K. The SCUBA-2 850 µm follow-up of WISE-selected, luminous dust-obscured Quasars. Publ. Astron. Soc. Pac. 129, 124101 (2017).

12. Jones, S. F. et al. Submillimetre observations of WISE-selected high-redshift, luminous, dusty galaxies. Mon. Not. R. Astron. Soc. 443, 146–157 (2014).

13. Díaz-Santos, T. et al. The multiple merger assembly of a hyperluminous obscured quasar at redshift 4.6. Science 362, 1034–1036 (2018).

14. Bischetti, M. et al. The WISS quasars project. V. ALMA reveals the assembly of a giant galaxy around a z = 4.4 hyper-luminous QSO. Astron. Astrophys. 617, A82 (2018).

15. Hill, R. et al. Megaparsec-scale structure around the protocluster core SPT2349-56 at z = 4.3. Mon. Not. R. Astron. Soc. 495, 3124–3159 (2020).

16. García-Vergara, C. et al. ALMA reveals a large overdensity and strong clustering of galaxies in quasar environments at z 4. Astrophys. J. 927, 65 (2022).

17. Fan, L. et al. The most luminous heavily obscured quasars have a high merger fraction: morphological study of WISE-selected hot dust-obscured galaxies. Astrophys. J. 822, L32 (2016).

18. Díaz-Santos, T. et al. Kinematics and star formation of high-redshift hot dust-obscured quasars as seen by ALMA. Astron. Astrophys. 654, A37 (2021).

19. Fan, L., Knudsen, K. K., Fogasy, J. & Drouart, G. ALMA detections of CO emission in the most luminous, heavily dust-obscured quasars at z > 3. Astrophys. J. 856, L5 (2018).

20. Tsai, C.-W. et al. The most luminous galaxies discovered by WISE. Astrophys. J. 805, 90 (2015).

21. Stanley, F. et al. Detection of H₂O and OH⁺ in z > 3 hot dust-obscured galaxies. Astron. Astrophys. 646, A178 (2021).

22. Borisova, E. et al. Ubiquitous giant Lyα Nebulae around the brightest Quasars at z-3.5 R revealed with MUSE. Astrophys. J. 831, 39 (2016).

23. Arrigoni Battaia, F. et al. QSO MUSEUM I: a sample of 61 extended Ly-α-emission nebulae surrounding z-3 quasars. Mon. Not. R. Astron. Soc. 482, 3162–3205 (2019).
24. Cantalupo, S., Arrigoni-Battaia, F., Prochaska, J. X., Hennawi, J. F. & Madau, P. A cosmic web filament revealed in Lyman-α emission around a luminous high-redshift quasar. *Nature* **506**, 63–66 (2014).

25. Jun, H. D. et al. Spectral classification and ionized gas outflows in z ~ 2 WISE-selected hot dust-obscured galaxies. *Astrophys. J.* **888**, 110 (2020).

26. Drake, A. B. et al. MUSE deep-fields: the Ly α luminosity function in the Hubble Deep Field-South at 2.91 < z < 6.64. *Mon. Not. R. Astron. Soc.* **471**, 267–278 (2017).

27. Herenz, E. C. et al. The MUSE-Wide Survey: A determination of the Lyman α emitter luminosity function at 3 < z < 6. *Astron. Astrophys.* **621**, A107 (2019).

28. Hayashino, T. et al. Large-scale structure of emission-line galaxies at z = 3.1. *Astron. J.* **128**, 2073–2079 (2004).

29. Bădescu, T. et al. Discovery of a protocluster associated with a Lyα Blob Pair at z = 2.3. *Astrophys. J.* **845**, 172 (2017).

30. Bacon, R. et al. The MUSE Extremely Deep Field: The cosmic web in emission at high redshift. *Astron. Astrophys.* **647**, A107 (2021).

31. Hu, W. et al. A Lyman-α protocluster at redshift 6.9. *Nat. Astron.* **5**, 485–490 (2021).

32. Venemans, B. P. et al. Properties of Lyα emitters around the radio galaxy MRC 0316 257. *Astron. Astrophys.* **431**, 793–812 (2005).

33. Venemans, B. P. et al. Protoclusters associated with z>2 radio galaxies. I. Characteristics of high redshift protoclusters. *Astron. Astrophys.* **461**, 823–845 (2007).

34. Hennawi, J. F., Prochaska, J. X., Cantalupo, S. & Arrigoni-Battaia, F. Quasar quartet embedded in giant nebula reveals rare massive structure in distant universe. *Science* **348**, 779–783 (2015).

35. Farrah, D. et al. The role of the most luminous obscured AGNs in galaxy assembly at z ~ 2. *Astrophys. J.* **844**, 106 (2017).

36. Ginolfi, M. et al. The infrared-luminous progenitors of high-z quasars. *Mont. Not. R. Astron. Soc.* **483**, 1256–1264 (2019).

37. Cantalupo, S. et al. The large- and small-scale properties of the intergalactic gas in the Slug Ly α nebula revealed by MUSE He II emission observations. *Mon. Not. R. Astron. Soc.* **483**, 5188–5204 (2019).

38. Yun, M. S., Ho, P. T. P. & Lo, K. Y. A high-resolution image of atomic hydrogen in the M81 group of galaxies. *Nature* **372**, 530–532 (1994).

39. Rizzo, F. et al. A dynamically cold disk galaxy in the early Universe. *Nature* **584**, 201–204 (2020).

40. Neeleman, M., Prochaska, J. X., Kanekar, N. & Rafelski, M. A cold, massive, rotating disk galaxy 1.5 billion years after the Big Bang. *Nature* **581**, 269–272 (2020).

41. Lelli, F. et al. A massive stellar bulge in a regularly rotating galaxy 1.2 billion years after the Big Bang. *Science* **371**, 713–716 (2021).

42. Pillepich, A. et al. First results from the TNG50 simulation: the evolution of stellar and gaseous discs across cosmic time. *Mon. Not. R. Astron. Soc.* **490**, 3196–3233 (2019).

43. Kretschmer, M., Dekel, A. & Teyssier, R. On the origin of surprisingly cold gas discs in galaxies at high redshift. *arXiv e-prints* arXiv:2103.06882 (2021).

44. Engel, H. et al. Most submillimeter galaxies are major mergers. *Astrophys. J.* **724**, 233–243 (2010).

45. Kohandel, M. et al. Kinematics of z ≥ 6 galaxies from [C II] line emission. *Mon. Not. R. Astron. Soc.* **487**, 3007–3020 (2019).

46. Kannan, R. et al. From discs to bulges: effect of mergers on the morphology of galaxies. *Mon. Not. R. Astron. Soc.* **452**, 4347–4360 (2015).

47. Díaz-Santos, T. et al. The strikingly uniform, highly turbulent interstellar medium of the most luminous galaxy in the universe. *Astrophys. J.* **816**, L6 (2016).

48. Planck Collaboration et al. Planck 2018 results. VI. Cosmological parameters. *Astron. Astrophys.* **641**, A6 (2020).

49. Weilbacher, P. M. et al. The data processing pipeline for the MUSE instrument. *Astron. Astrophys.* **641**, A28 (2020).

50. Bacon, R. et al. The MUSE 3D view of the Hubble Deep Field South. *Astron. Astrophys.* **575**, A75 (2015).

51. Umehata, H. et al. Gas filaments of the cosmic web located around active galaxies in a protocluster. *Science* **366**, 97–100 (2019).

52. Fossati, M. et al. MUSE analysis of gas around galaxies (MAGG)-III. The gas and galaxy environment of z = 3-4.5 quasars. *Mon. Not. R. Astronomical Soc.* **503**, 3044–3064 (2021).

53. Runnoe, J. C., Brotherton, M. S. & Shang, Z. Updating quasar bolometric luminosity corrections. *Mon. Not. R. Astron. Soc.* **422**, 478–493 (2012).

54. Marziani, P. et al. Quasar massive ionized outflows traced by CIV λ1549 and [OIII]λ4959,5007. *Front. Astron. Space Sci.* **4**, 16 (2017).

55. Vietri, G. et al. The WISSH quasars project. IV. Broad line region versus kiloparsec-scale winds. *Astron. Astrophys.* **617**, A81 (2018).

56. Verhamme, A. et al. Recovering the systemic redshift of galaxies from their Lyman alpha line profile. *Mon. Not. R. Astron. Soc.* **478**, L60–L65 (2018).

57. Coil, A. L. The Large-Scale Structure of the Universe. In Oswalt, T. D. & Keel, W. C. (eds.) *Planets, Stars and Stellar Systems. Volume 6: Extragalactic Astronomy and Cosmology*, 6, 387, https://doi.org/10.1007/978-94-007-5609-0_8 (2013).

58. Kovač, K., Somerville, R. S., Rhoads, J. E., Malhotra, S. & Wang, J. Clustering of Lyα Emitters at z ~ 4.5. *Astrophys. J.* **668**, 15–22 (2007).

59. Herrero Alonso, Y. et al. The MUSE-Wide survey: Three-dimensional clustering analysis of Lyman-α emitters at 3.3 < z < 6. *Astron. Astrophys.* **653**, A136 (2021).

60. Ouchi, M. et al. Systematic identification of LAEs for visible exploration and reionization research using Subaru HSC (SILVER-RUSH). I. Program strategy and clustering properties of ~2000 Lyα emitters at z=6-7 over the 0.3–0.5 Gpc² survey area. *Publ. Astron. Soc. Jpn.* **70**, S13 (2018).

61. Kravtsov, A. V. & Borgani, S. Formation of galaxy clusters. *Annu. Rev. Astron. Astrophys.* **50**, 353–409 (2012).

62. Klypin, A., Prada, F., Betancort-Rijo, J. & Albareti, F. D. Density distribution of the cosmological matter field. *Mon. Not. R. Astron. Soc.* **481**, 4588–4601 (2018).

63. Klypin, A., Yepes, G., Gottlöber, S., Prada, F. & Heß, S. MultiDark simulations: the story of dark matter halo concentrations and density profiles. *Mon. Not. R. Astron. Soc.* **457**, 4340–4359 (2016).

64. Hoag, A. et al. Constraining Lyman-alpha spatial offsets at 3 < z < 5.5 from VANDELS slit spectroscopy. *Mon. Not. R. Astron. Soc.* **488**, 706–719 (2019).

65. Jr. Kennicutt, R. C. Star formation in galaxies along the hubble sequence. *Annu. Rev. Astron. Astrophys.* **36**, 189–232 (1998).

66. Leclercq, F. et al. The MUSE Hubble Ultra Deep Field Survey. XIII. Spatially resolved spectral properties of Lyman α haloes around star-forming galaxies at z > 3. *Astron. Astrophys.* **635**, A82 (2020).

67. McMullin, J. P., Waters, B., Schiebel, D., Young, W. & Golap, K. CASA Architecture and Applications. In Shaw, R. A., Hill, F. & Bell, D. J. (eds.) *Astronomical Data Analysis Software and Systems XVI*, 376 of Astronomical Society of the Pacific Conference Series, 127 (2007).

68. Tacconi, L. J. et al. PHIBSS: unified scaling relations of gas depletion time and molecular gas fractions. *Astrophys. J.* **853**, 179 (2018).

69. Carilli, C. L. & Walter, F. Cool gas in high-redshift galaxies. *Annu. Rev. Astron. Astrophys.* **51**, 105–161 (2013).

70. Hunt, L. K., Tortora, C., Ginolfi, M. & Schneider, R. Scaling relations and baryonic cycling in local star-forming galaxies. II. Gas content and star-formation efficiency. *Astron. Astrophys.* **643**, A180 (2020).

71. Maiolino, R. & Mannucci, F. De re metallica: the cosmic chemical evolution of galaxies. *Astron. Astrophys. Rev.* **27**, 3 (2019).
72. Di Teodoro, E. M. & Fraternali, F. 3D BAROLO: a new 3D algorithm to derive rotation curves of galaxies. *Mon. Not. R. Astron. Soc.* **451**, 3021–3033 (2015).

73. Jones, G. C. et al. The ALPINE-ALMA [C II] Survey: kinematic diversity and rotation in massive star-forming galaxies at z 4.4-5.9. *Mon. Not. R. Astron. Soc.* **507**, 3540–3563 (2021).

74. Whiting, M. T. DUCHAMP: a 3D source finder for spectral-line data. *Mon. Not. R. Astron. Soc.* **421**, 3242–3256 (2012).

## Acknowledgements

ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), NSC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The CosmoSim database used in this paper is a service by the Leibniz-Institute for Astrophysics Potsdam (AIP). The Cosmic Dawn Center (DAWN) is funded by the Danish National Research Foundation under grant No. 140. G.C.J. acknowledges funding from European Research Council (ERC) Advanced Grants 695671 "QUENCH" and 789056 "FirstGalaxies" under the European Union's Horizon 2020 research and innovation programme, as well as support by the Science and Technology Facilities Council (STFC). F.R. acknowledges support from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No. 847523 'INTERACTIONS'. S.CAN gratefully acknowledges support from the ERC under the European Union's Horizon 2020 research and innovation programme grant agreement No 864361. R.M. acknowledges support by the Science and Technology Facilities Council (STFC), by the ERC Advanced Grant 695671 "QUENCH" and also funding from a research professorship from the Royal Society. E.P. acknowledges support from PRIN MIUR project "Black Hole winds and the Baryon Life Cycle of Galaxies: the stone-guest at the galaxy evolution supper", contract #2017PH3WAT. L.G. and R.S. acknowledge support from the Amaldi Research Center funded by the MIUR program "Dipartimento di Eccellenza" (CUP:B81118001170001).

## Author contributions

M.G., E.P., L.Z. and R.M. designed the MUSE observing strategy and wrote the MUSE proposals, with the help of F.A.B., S. Cantalupo, S. Carniani, L.P., G.C.J., A.T., R.V., L.G. and R.S.; M.G. analyzed the MUSE data, with the help of A.T.; M.G., G.C.J., S.Carniani, and F.R. analyzed the ALMA data. M.G. and L.P. analyzed the HST data. C.D.B., K.K., P.L., V.M.,

F.S., and A.V. contributed to the interpretation of the MUSE data. N.M. led the calculation on the physical significance of the overdensity. All authors contributed to the interpretation and the discussion of the scientific results. M.G. wrote the manuscripts and all authors provided comments.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at

https://doi.org/10.1038/s41467-022-32297-x.

**Correspondence** and requests for materials should be addressed to M. Ginolfi.

**Peer review information** *Nature Communications* thanks John Peacock and the other, anonymous, reviewer for their contribution to the peer review of this work.

**Reprints and permission information** is available at

http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2022