A Multiwavelength Study of ELAN Environments (AMUSE²). Mass Budget, Satellites Spin Alignment, and Gas Infall in a Massive z ∼ 3 Quasar Host Halo

We used the Large APEX BOlometer CAmera (LABOCA; Siringo et al. 2009) on the APEX telescope to map a field of ∼68 arcmin² around the ELAN hosting QSO. The 295 bolometers of LABOCA operate at an effective frequency of 345 GHz (or a wavelength of 870 μm), and the instrument is characterized by a main beam of 19''. The observations were conducted in service mode in October 2016 (ID: 098.A-0828(B); PI: F. Arrigoni Battaia) with zenith opacities between 0.2 and 0.4 at 870 μm. The field has been covered with a raster of spiral scanning mode, which optimizes the sampling of the field of view with the LABOCA instrument. The total integration time on source resulted in 22 hours consisting of 176 scans of 7.5 minutes each. The observations have been acquired with regular standard calibrations for pointing, focus, and flux calibration (see, e.g., Siringo et al. 2009 for details).

The data reduction was performed with the Python-based BOlometer data Analysis Software package (BoA; Schuller 2012) following the steps indicated in Siringo et al. (2009) and Schuller et al. (2009). Specifically, BoA processes LABOCA data including flux calibration, opacity correction, noise removal, and despiking of the timestreams. We ran BoA using the default reduction script reduce-map-weaksource.boa, which also filters out the low-frequency noise below 0.3 Hz. The scans are then co-added after being variance-weighted. The final outputs are a beam-smoothed flux density and a noise map. The final map achieves an rms noise level of 2.6–3 mJy beam⁻¹ in its central part. We show the map for the full area covered in Appendix A.

The observations with the Submillimetre Common-User Bolometer Array 2 (SCUBA-2) for this ELAN field were conducted at JCMT during flexible observing on 2018 February 12 and March 29 (program ID: M18AP054; PI: M. Fumagalli) under good weather conditions (band 1 and 2; τ_225GHz ≤ 0.07). The SCUBA-2 instrument observes simultaneously the same field at 850 and 450 μm, with an effective beam FWHM of 146 and 98, respectively (Dempsey et al. 2013). The observations were performed with a Daisy pattern covering ≃137 in diameter, and were centered at the location of QSO (and thus the ELAN). To facilitate the scheduling we divided the observations in five scans/cycles of about 30 minutes, for a total of 2.5 hours.

For the data reduction we closely followed the procedures in Chen et al. (2013a) and Arrigoni Battaia et al. (2018b). In brief, we reduced the data using the Dynamic Iterative Map Maker (DIMM) included in the Sub-Millimetre User Reduction Facility (SMURF) package from the STARLINK software (Jenness et al. 2011; Chapin et al. 2013). We adopted the standard configuration file dimmconfig_blank_field.lis for our science purposes. We thus reduced each scan and the MOSAIC_JCMT_IMAGES recipe in the Pipeline for Combining and Analyzing Reduced Data (PICARD; Jenness et al. 2008) was used to coadd the reduced scans into the final maps.

We applied to these final maps a standard matched filter to increase the point-source detectability, using the PICARD recipe SCUBA2_MATCHED_FILTER. We adopted the standard flux conversion factors (FCFs; 491 Jy pW⁻¹ for 450 μm and 537 Jy pW⁻¹ for 850 μm) with 10% upward corrections for flux calibration. The relative calibration accuracy is shown to be stable and good to 10% at 450 μm and 5% at 850 μm (Dempsey et al. 2013).

The final noise level at the location of the ELAN is 1.01 mJy beam⁻¹ and 10.97 mJy beam⁻¹ at 850 μm and 450 μm, respectively. Appendix A presents the SCUBA-2 maps for the whole field of view. In the remainder of this work we focus only on the ELAN location.

We performed the SMA (Ho et al. 2004) observations of this ELAN (Project code: 2017A-S015; PI: F. Arrigoni Battaia) on 2017 June 21 (UTC 3:00-8:30), June 27 (UTC 3:00-7:30), and July 10 (UTC 3:30-6:30), in the compact array configuration. However, for the observations on June 27, we only used the data taken after UTC 6:30 because we noticed a large antenna pointing error before then. The atmospheric opacity at 225 GHz (τ_225GHz) was 0.1–0.15, ∼0.1, and ∼0.05 during these three tracks of observations.

The observations were carried out in the dual-receiver mode by tuning the 345 and 400 GHz receivers to the same observing frequencies. These two receivers took left and right circular polarization, respectively, and covered the observing frequency of 329–337 GHz in the lower sideband and 345–353 GHz in the upper sideband. Correlations were performed by the SMA Wideband Astronomical ROACH2 Machine (SWARM), which sampled individual sidebands with 16,384 × 4 spectral channels. The integration time was 30 s. Prior to data calibration, we binned every 16 spectral channels to reduce the file sizes. The observations on our target source cover the uv distance range of ∼8.5–88.5 k λ.

The target sources were observed in scans of 12 minutes in duration, which were bracketed by scans on the gain calibration quasar source 1058+015 with 3 minutes in duration. We observed Titan in the first two tracks, and observed Callisto in the last track for absolute flux calibrations. We followed the standard data calibration strategy of SMA. The application of Tsys information and the absolute flux, passband, and gain calibrations were carried out using the MIR IDL software package (Qi 2003). The absolute flux scalings were derived by comparing the visibility amplitudes of the gain calibrators with those of the absolute flux calibrators (i.e., Titan and Callisto). The derived and applied fluxes of 1058+015 were 2.5 Jy in the first two tracks, and 2.7 Jy in the last track. We nominally quote the ∼15% typical absolute flux calibration error of SMA.

After calibration, the zeroth-order fitting of continuum levels and the joint weighted imaging of all continuum data were performed using the Miriad software package (Sault et al. 1995). We performed zeroth-order multifrequency imaging combining the upper- and lower-sideband data, to produce a sensitive continuum image at the central observing frequency (i.e., the local oscillator frequency). Due to the different performance of the 345 and 400 GHz receivers at the same observing frequency, it would be incorrect to treat half of the difference of the parallel hand correlations (i.e., (LL − RR)/2) as the thermal noise map. Instead, we first smoothed the upper-sideband image to the angular resolution of the lower-sideband image and then took half of their difference as the approximated noise map. Using natural weighting, we obtained a θ_maj × ${\theta }_{{\rm{\min }}}$ = 24 × 20 (P.A. = 67°) synthesized beam, and an rms noise level of 1.4 mJy beam⁻¹.

We performed four epochs of ALMA observations toward this ELAN (Project code: 2017.1.00560.S; PI: F. Arrigoni Battaia), on 2018 March 23, 24, 26, and 27 (UTC) to constrain the CO(5–4) line emission (ν_rest = 576.267 GHz) and its underlying 2 mm continuum. The pointing and phase referencing center is R.A. (J2000) = 10^h20^m0942, and decl. (J2000) = 10°40'0871. Combining all existing data yields an overall uv distance range covered of 12–740 m.

The spectral setup of all our observations is identical. There were two 2 GHz wide spectral windows (channel spacing 15.625 MHz) centered at the sky frequencies 149.514 and 151.201, and two 1.875 GHz wide spectral windows (channel spacing 3.906 MHz) centered at the sky frequencies 137.784 and 139.472 GHz. The latter two spectral windows with channel width of about 8.5 km s⁻¹ are expected to encompass the CO(5–4) emission.

For all four epochs of observations, quasar J1058+0133 was chosen as the flux and passband calibrator. We assume that J1058+0133 has an absolute flux of 3.09 Jy and a spectral index of −0.46 at the reference frequency 144.493 GHz, based on interpolating the calibrator grid survey measurements taken in Band 3 (∼91 and 103 GHz) on 2018 March 25, and in Band 7 (∼343 GHz) on 2018 February 09. Based on the results of the calibrator grid survey, we expect a nominal absolute flux error of ∼10%, and an in-band spectral index error of ∼0.1, as the grid survey measurements are sparsely sampled in time. We observed quasar J1025+1253 approximately every 11 minutes for complex gain calibration.

We calibrated the data using the CASA software package (McMullin et al. 2007) version 5.1. The derived fluxes of the gain calibrator J1025+1253 were in the range 0.42–0.48 Jy. We fit the continuum baselines using the CASA task uvcontsub. We jointly imaged all continuum data using the CASA task clean, which produced the Stokes I image by averaging the parallel linear correlation data (i.e., I = (XX +YY)/2). Our target sources are presumably weakly or not polarized. Therefore, we regarded the (XX − YY)/2 image as an approximated thermal noise map. The Briggs Robust = 2 weighted image achieved a synthesized beam of ${\theta }_{\mathrm{maj}}\times {\theta }_{\min }$ = 095 × 094 (P.A. = −53), and an rms noise of 4.7 μJy beam⁻¹.

For the spectral windows including the CO(5–4) emission, we generated the continuum from the channels not affected by line emission, and subtracted it from the data. Continuum-subtracted data cubes were created with the CASA task tclean, using Briggs cleaning with a robustness parameter of 2 (corresponding to natural visibility weights). This approach maximizes the signal-to-noise ratio and it is frequently used in observations of high-z quasars (e.g., Decarli et al. 2018; Bischetti et al. 2021).

We observed this ELAN in the H band with the High Acuity Wide-field K-band Imager (HAWK-I; Casali et al. 2006) on the Unit Telescope 4 (UT4; Yepun) of the VLT in service mode under project 0102.C-0589(D) (PI: F. Vogt). In this work we only focus on the H-band observations acquired with clear weather, i.e., 2019 February 15, 23 and March 9, 22, 23. HAWK-I has a field of view of $7\buildrel{\,\prime}\over{.} 5\times 7\buildrel{\,\prime}\over{.} 5$ covered by an array of 2 × 2 Hawaii-2RG detectors separated by 15'' gaps. The observational strategy consisted of three fast 60 s H-band exposures per observing block (OB), applying a dithering within a jitter box of 15''. The ELAN system was always acquired in the fourth quadrant, Q4, of the detector array. The total on-source time for our clear weather observations consists of 12 OBs, i.e., 36 minutes on source for the H band.

We reduced the data with the standard ESO pipeline version 2.4.3 for HAWK-I. ²⁰ In brief, the data were corrected for dark current and were flat-fielded. The sky subtraction was performed using the algorithm pawsky _mask, which iteratively estimates the background by stacking with rejection the science frames and by constructing a mask for the objects in the data. The sky estimation ends once the number of masked pixels converges. The photometry of the images was calibrated with Two Micron All Sky Survey (2MASS) stars in the field of view of our observations, achieving a 1σ AB magnitude limit of 26.0 mag in a 1 arcsec² aperture. The intrinsic uncertainty on the photometric calibration is 0.1 mag. The astrometry was calibrated against the 2MASS catalog (about 20 stars), with an average error in the coordinates fit of ∼02. This astrometry calibration agrees well with the GAIA Data Release (DR) 2 catalog (Gaia Collaboration et al. 2018). The seeing in the final combined image is of 05.

The first data we acquired on this system, the 870 μm APEX/LABOCA data, revealed a 4.8σ detection of 12.5 ± 2.6 mJy at the position of the ELAN (see top-left panel in Figure 1). This surprisingly strong detection in an ELAN was then confirmed by the deeper 850 μm JCMT/SCUBA-2 observations, with a flux density of 12.7 ± 1.0 mJy (see bottom-left panel in Figure 1). We corrected these observed flux densities for flux boosting (see Appendix B), obtaining f_Deboosted = 10.5 ± 2.2, and 11.7 ± 0.9 mJy for the LABOCA and SCUBA-2 detections, respectively (Table 1).

Zoom In Zoom Out Reset image size

Given the radio-quietness of all the sources within the ELAN (Arrigoni Battaia et al. 2018a), the detection traces thermal dust emission from embedded starbusting galaxies. Following Cowie et al. (2017), the strong detected fluxes imply a star formation rate (SFR) of SFR = 1680 ± 338 M_⊙ yr⁻¹. These unexpectedly bright single-dish detections have been a fundamental stepping stone for the follow-up observations with interferometers.

The SCUBA-2 450 μm data are not deep enough to detect emission from the ELAN. We report in Table 1 the upper limit at 450 μm obtained at the location of the 850 μm detection.

We extracted continuum sources from the SMA continuum map (top-right panel in Figure 1) using the same algorithm described in, e.g., Arrigoni Battaia et al. (2018b), but working with the SMA beam of the current data set. Briefly, the algorithm iteratively searches for maxima in the S/N map (Figure 2) while subtracting (at their locations) mock sources normalized to those peaks. The iterations are stopped once S/N = 2 is reached. The S/N peaks found by the algorithm are included in a source candidate catalog. In the current case, the algorithm found seven sources. Subsequently, the same algorithm is applied to the negative data set down to the same S/N threshold to estimate the number of spurious sources and clean the aforementioned catalog. We find no spurious sources within a radius of R = 7'' from the center of the map, one such source in the annuli within 7'' < R < 105, and four for R > 105. Therefore, we consider reliable five of the seven source candidates: the three sources detected within a radius of R = 7'' (QSO2, AGN1, SMG1) and two sources for R > 105 with detections in the ALMA data (QSO and LAE1). These five sources are indicated as yellow diamonds in Figure 2. Instead, the two potentially spurious sources are indicated with cyan diamonds, and are located, respectively, in the 7'' < R < 105 and R > 105 regions. This analysis is confirmed by the absence of emission at these two locations in the ALMA continuum map (see Section 3.3), while all the other sources are very close to the positions of known sources associated with the ELAN or with ALMA detections (see Section 3.3). As said, the detected sources are QSO, QSO2, AGN1, LAE1, and a newly discovered source SMG1. The positions, S/N, and fluxes for the five detections are listed in Table 2, together with the deboosted fluxes estimated as explained in Appendix C. Summing up the deboosted fluxes of all detected sources, we find agreement within uncertainties with the detections in the single-dish data sets, 14.7 ± 2.8 mJy. Therefore, all the continuum emission detected by LABOCA and SCUBA-2 is ascribed to compact sources.

Zoom In Zoom Out Reset image size

Table 2. Continuum Measurements from ALMA, SMA, HAWK-I, and MUSE

ID	fALMA	S/N	${f}_{\mathrm{ALMA}}^{\mathrm{Deboosted}}$	R.A.	Decl.	fSMA	S/N	${f}_{\mathrm{SMA}}^{\mathrm{Deboosted}}$	H	i	r
	(mJy)		(mJy)	(J2000)	(J2000)	(mJy)		(mJy)	(AB)	(AB)	(AB)
QSO	0.24	34.6	0.23 ± 0.01	10:20:10.00	+10:40:02.7	6.8 b	4.9	5.7 ± 1.8	17.1 ± 0.1	17.98 ± 0.01	17.67 ± 0.01
QSO2	0.06	8.9	0.06 ± 0.01	10:20:09.56	+10:40:05.3	3.3 b	2.4	2.0 ± 1.0	23.5 ± 0.1	24.30 ± 0.02	24.10 ± 0.01
LAE1	0.17	23.9	0.17 ± 0.01	10:20:10.15	+10:40:10.6	3.6 b	2.6	2.4 ± 1.1	24.6 ± 0.6	25.43 ± 0.05	25.32 ± 0.05
LAE2	<0.01 a	⋯	⋯	⋯	⋯	<2.8 a	⋯	⋯	>26.0 d	>27.8 d	>27.8 d
AGN1	0.19	18.2	0.13 ± 0.01	10:20:09.83	+10:40:14.7	4.0 b	2.8	2.8 ± 1.2	24.7 ± 0.7	26.20 ± 0.20	26.30 ± 0.20
SMG1	0.03	4.4	0.02 ± 0.01	10:20:09.18	+10:40:13.4	3.1 b	2.2	1.8 ± 0.9	24.6 ± 0.3	26.14 ± 0.12	26.88 ± 0.23
S1	0.04	6.6	0.04 ± 0.01	10:20:09.41	+10:40:04.9	<2.8 c	⋯	⋯	24.8 ± 0.3	27.18 ± 0.31	>28.6 d
S2	0.03	4.1	0.03 ± 0.01	10:20:10.16	+10:40:08.6	<2.8 c	⋯	⋯	23.2 ± 0.1	24.77 ± 0.03	24.75 ± 0.03
S3	0.03	4.1	0.03 ± 0.01	10:20:08.78	+10:40:16.3	<2.8 c	⋯	⋯	>26.0 d	>27.8 d	>27.8d

Notes.

^a 2σ upper limit at the Lyα position from Arrigoni Battaia et al. (2018a). ^b The coordinates of the SMA detections differ from the ALMA detections due to their lower S/N. For completeness, we report them for each source: [10:20:09.9, +10:40:03], [10:20:09.5, +10:40:05], [10:20:10.1, +10:40:10], [10:20:09.7, +10:40:16], [10:20:09.2, +10:40:14] for QSO, QSO2, LAE1, AGN1, and SMG1, respectively. ^c 2σ upper limit at the ALMA position. ^d 2σ limit at the ALMA position within an aperture of 2'' diameter.

Download table as:  ASCIITypeset image

We finally stress that the alignment of each individual SMA S/N > 2 source with the location of an ALMA detection (see Section 3.3) strongly suggests that the reported sources are reliable. Within R < 12'' from the phase center, where all our detections are located, we expect to have 144 independent beam areas. Following Gaussian noise, the chance to have a >2σ positive noise peak is ∼2.25%, so there could be ∼3 noise peaks above >2σ. The probability of having a >2σ noise peak at close locations to an ALMA source is therefore ${P}_{\det }\sim 3/144=0.0225$. The probability of having five noise peaks (as the detected sources QSO, QSO2, LAE1, AGN1, SMG1) with S/N > 2 and at close locations to ALMA detections is therefore very low. We estimated it to be ${P}_{\det }^{5}\sim 5.8\times {10}^{-9}$ (confirmed using Monte Carlo simulations).

We extracted sources from the ALMA continuum map at 2 mm (bottom-right panel in Figure 1) following the same method as for the SMA data, but using the ALMA beam. We considered as reliable only sources with S/N > 3.7. Indeed, above this threshold we did not find any spurious source in the negative map. Using this threshold, we found eight detections, shown as black circles in Figure 3: (i) the known sources QSO, QSO2, LAE1, and AGN1, (ii) the source SMG1 discovered with SMA, and (iii) three additional sources, which we dubbed S1, S2, and S3. The coordinates, fluxes, and S/N of all the sources, as well as their deboosted fluxes estimated as explained in Appendix D, are listed in Table 2.

Zoom In Zoom Out Reset image size

As evident from Figure 3, some of the ALMA detections have shifts of ∼few arcseconds with respect to the Lyα locations and to the SMA locations. While the latter is likely due to the low S/N of the SMA detections, we discuss in detail the shifts with respect to the Lyα locations in Section 6.

We inspected the H-band HAWK-I data presented in Section 2.5 at the location of the sources so far discussed. As can be seen in Figure 4, we find clear detections for QSO, QSO2, and S2, and fainter emission at the locations of SMG1, AGN1, LAE1, and S1. LAE2 and S3 are undetected at the current depth. We extract magnitudes with apertures of 2'' diameter (4× the seeing) for all the sources except QSO. Indeed, its flux is better captured by a 3'' diameter aperture. In Table 2 we list the magnitudes obtained with apertures centered on the ALMA detections, if available.

Zoom In Zoom Out Reset image size

Further, we obtained the observed optical magnitude, i and r, of all sources within the aforementioned apertures from the MUSE data presented in Arrigoni Battaia et al. (2018a). These magnitudes are listed in Table 2. We note that the MUSE i-band magnitude for AGN1 is different from the magnitude listed in Arrigoni Battaia et al. (2018a) because in that study the authors rely on the compact Lyα emission for determining the source position. However, there is an important shift between Lyα and the near- and far-infrared continua detected in this work (Figure 4; see discussion in Section 6).

We rely on the publicly available code LineSeeker (González-López et al. 2017) to robustly identify sources with CO(5–4) line emission in the ALMA observations. Indeed LineSeeker has been developed to systematically search for line emissions in ALMA data and quantify their significance (e.g., González-López et al. 2019). The code looks for signal in the ALMA data cubes on a channel-by-channel basis after convolving the data along the spectral axis with Gaussian kernels of different spectral widths. A list of line candidates is obtained by joining detected signals on different channels using the DBSCAN algorithm (Ester et al. 1996). LineSeeker finally provides an S/N for each emission-line candidate. This S/N is defined as the maximum value obtained from the different convolutions. Further, it runs the search algorithm also on the negative datacube and on simulated cubes to estimate the significance of the line candidates' S/N, providing probabilities for the line to be false.

Here we focus on 100% fidelity sources, i.e., sources whose probability of being false is 0 and whose S/N is larger than any detection in the negative data. In this way, we obtained four line detections at >10σ from four sources detected also in the continuum, QSO, QSO2, AGN1, and LAE1. All the other known sources do not show evidence of CO(5–4) emission consistent with the system redshift. We stress that we looked for sources with LineSeeker down to S/N = 5, but all additional line-emission candidates are found at the edge of the primary beam with very low fidelity, and therefore they are not reliable.

We then extracted the spectrum for each detected source using the minimum aperture that maximizes the flux densities. An aperture 2× the synthetized beam fulfilled this criterion. Figure 5 shows the four spectra binned to channels of 23.4 MHz (or ∼51 km s⁻¹), with the emission above 2× the rms highlighted in each spectrum. The spectra are shown indicating the velocity shift with respect to the QSO systemic redshift estimated from C iv (z = 3.164; Arrigoni Battaia et al. 2018a) after correcting for the expected velocity shift between C iv and systemic (Shen et al. 2016). In the remainder of the paper we will consider the redshift of the CO emission as systemic redshift for each detected object.

Zoom In Zoom Out Reset image size

All lines detected show velocity widths >200 km s⁻¹, when estimated using their second moment, though their shapes are relatively boxy and the widths of the highlighted velocity ranges in Figure 5 are as high as ∼1000 km s⁻¹ (see Table 3). The CO(5–4) line profiles for QSO and AGN1 are complex: QSO has three tentative peaks, while AGN1 has a profile with higher flux densities at the edges of the line. The profile of AGN1 is suggestive of a molecular gas disk, similar to what is seen in other AGN host galaxies (e.g., in low-redshift radio galaxies; Ocaña Flaquer et al. 2010). We will further discuss the Lyα and CO(5–4) velocity shifts and line shapes in Section 6. Table 3 lists the rms for these spectra and all the lines properties extracted using the highlighted velocities, i.e., redshifts, line widths, fluxes, line luminosities.

Table 3. CO(5–4) Detections from ALMA and their Derived Galaxy Properties

	QSO	QSO2	AGN1	LAE1
rms noise per 23.4 MHz [μ Jy]	118	117	133	103
S/N a	22.2	12.4	11.3	18.0
zCO(5−4) b	3.1695 ± 0.0008	3.1582 ± 0.0006	3.1777 ± 0.0009	3.1727 ± 0.0007
CO(5–4) line width [km s−1] c	418 ± 101	211 ± 103	443 ± 202	246 ± 96
ICO(5−4) [Jy km s−1]	0.43 ± 0.03	0.22 ± 0.02	0.28 ± 0.03	0.34 ± 0.02
LCO(5−4) [107 L⊙]	4.5 ± 0.3	2.3 ± 0.2	2.9 ± 0.3	3.6 ± 0.2
${L}_{\mathrm{CO}(5--4)}^{{\prime} }$ [109 K km s−1 pc2]	7.5 ± 0.4	3.9 ± 0.3	4.8 ± 0.5	6.0 ± 0.3
${L}_{\mathrm{CO}(1--0)}^{{\prime} }$ [109 K km s−1 pc2] d	18.7 ± 1.1	9.6 ± 0.8	12.1 ± 1.3	14.9 ± 0.8
2 mm major axis [''] e	1.07 ± 0.04	1.68 ± 0.21	1.12 ± 0.07	1.06 ± 0.05
2 mm minor axis [''] e	0.98 ± 0.03	1.29 ± 0.14	1.04 ± 0.06	1.03 ± 0.05
2 mm dec. major axis [''] f	0.52 ± 0.09	1.39 ± 0.27	0.64 ± 0.16	0.52 ± 0.17
2 mm dec. minor axis [''] f	0.30 ± 0.15	0.89 ± 0.22	0.44 ± 0.19	0.43 ± 0.14
R2mm [kpc] g	2.0 ± 0.3	5.3 ± 1.0	2.4 ± 0.6	2.0 ± 0.6
CO(5–4) major axis [''] e	1.30 ± 0.11	1.67 ± 0.15	1.75 ± 0.18	1.39 ± 0.11
CO(5–4) minor axis [''] e	1.12 ± 0.09	1.36 ± 0.11	1.57 ± 0.15	1.29 ± 0.10
CO(5–4) dec. major axis [''] f	0.78 ± 0.20	1.33 ± 0.20	1.41 ± 0.23	0.95 ± 0.19
CO(5–4) dec. minor axis [''] f	0.49 ± 0.27	0.89 ± 0.18	1.21 ± 0.22	0.78 ± 0.17
RCO(5−4) [kpc] g	3.0 ± 0.8	5.1 ± 0.8	5.4 ± 0.9	3.6 ± 0.7
Lbol [erg s−1] h	(4.6 ± 0.2) × 1047	(8.8 ± 0.9) × 1045	(1.0 ± 0.5) × 1046	(1.3 ± 0.2) × 1045
LIR [1012 L⊙] i	4.7 ± 0.2	7.3 ± 1.9	3.0 ± 1.0	1.2 ± 0.2
${L}_{\mathrm{IR}}^{\mathrm{total}}$ [L⊙] l	(1.7 ± 0.1) × 1014	(7.8 ± 2.3) × 1012	(1.0 ± 0.3) × 1013	(1.2 ± 0.2) × 1012
SFR [M⊙ yr−1] m	521 ± 74	183 ± 49	104 ± 21	50 ± 11
SFRIR [M⊙ yr−1] n	700 ± 30	1087 ± 283	447 ± 149	179 ± 30
Mdust [M⊙] h	(6.7 ± 1.3) × 108	(1.5 ± 0.4) × 108	(5.7 ± 1.3) × 108	(1.6 ± 0.4) × 109
Mdust [M⊙] o	(9 ± 3) × 108	(3 ± 2) × 108	(5 ± 2) × 108	(4 ± 2) × 108
${M}_{{{\rm{H}}}_{2},\mathrm{CO}}$ [M⊙] p	(1.5 ± 0.1) × 1010	(7.7 ± 0.7) × 109	(1.0 ± 0.1) × 1010	(1.2 ± 0.1) × 1010
${r}_{{{\rm{H}}}_{2},\mathrm{dust}}$ q	22 ± 4	51 ± 14	18 ± 4	8 ± 2
Mstar [M⊙] r	(1.2 ± 0.8) × 1011	(1.4 ± 0.4) × 1010	(1.1 ± 0.3) × 1010	(4.0 ± 0.5) × 109
MDM [M⊙] s	$({6.9}_{-5.6}^{+19.4})\times {10}^{12}$	(6.8 ± 1.0) × 1011	$({6.0}_{-1.6}^{+0.7})\times {10}^{11}$	$({3.7}_{-0.2}^{+0.3})\times {10}^{11}$

Notes.

^a Integrated signal-to-noise ratio. ^b Redshift corresponding to the first moment of the detected line. ^c Line width computed from the second moment of each detected line following, e.g., Equation (3) of Birkin et al. (2021). Given the relatively boxy shape of the lines we report here for completeness the line width corresponding to the colored range in Figure 5: 924, 513, 976, and 667 km s⁻¹ for QSO, QSO2, AGN1, and LAE1, respectively. ^d CO(1–0) luminosities obtained assuming I_CO(5−4)/I_CO(1−0) = 10 (see Section 4.2 for details). ^e Observed (i.e., beam-convolved) sizes of the 2 mm continuum and CO(5–4) emitting region from 2D Gaussian fit of the ALMA maps (see Section 3.5 for details). ^f Beam-deconvolved sizes of the 2 mm continuum and CO(5–4) emitting region from 2D Gaussian fit of the ALMA maps (see Section 3.5 for details). ^g Effective radius of the 2 mm continuum and CO(5–4) emitting region, defined as their major semiaxis (see Section 3.5 for details). ^h Obtained with the CIGALE fit. ⁱ Luminosity obtained by integrating only the dust emission of the SED due to star formation, in the rest-frame wavelength range 8–1000 μm. ^l Luminosity obtained by integrating the total SED in the rest-frame wavelength range 8–1000 μm. ^m SFR obtained with the CIGALE fit. ⁿ SFR obtained from the IR, assuming the relation (SFR _IR/[M_⊙ yr⁻¹]) = 3.88 × 10⁻⁴⁴(L_IR/[erg s⁻¹]) of Murphy et al. (2011). ^o Dust mass obtained using a modified blackbody (see Section 4.1 for details). ^p Molecular gas mass derived from the CO(5–4) line, assuming (i) a ratio I_CO(5−4)/I_CO(1−0) = 10 (or r₅₁ = 0.4), and (ii) a luminosity to gas mass conversion factor of α_CO = 0.8 ${M}_{\odot }{\left({\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{pc}}^{2}\right)}^{-1}$. The reported error on these measurements only includes the error on ${L}_{\mathrm{CO}(5-4)}^{{\prime} }$. Large uncertainties are expected due to the unconstrained r₅₁ and the known uncertainties on α_CO (see Section 4.2 for details). ^q Molecular gas-to-dust mass ratio computed using the dust mass estimated by CIGALE. The ratios obtained with the other dust mass estimates are consistent within 2σ. The errors on ${r}_{{{\rm{H}}}_{2},\mathrm{dust}}$ do not include the large uncertainties expected for ${M}_{{{\rm{H}}}_{2},\mathrm{CO}}$ (see Section 4.2 for details). ^r Stellar mass estimated by CIGALE. ^s Dark matter halo mass estimated assuming the stellar mass—halo mass relation of Moster et al. (2018), and interpolating their models for the redshift of each source.

Download table as:  ASCIITypeset image

Figure 6 shows cutouts of 7'' × 7'' (or about 53 kpc × 53 kpc) of the moment zero map and first and second-moment maps, together with the continuum at each source position. The first-moment maps are computed with respect to each source CO(5–4) redshift, as listed in Table 3 and shown in Figure 5. As can be seen from Figure 6, the 2 mm continuum and the CO(5–4) emission are found at consistent sky locations within uncertainties. For this reason and to avoid confusion, we only report the coordinates for the continuum (Table 2). Following Gaussian deconvolution theory, the position accuracy that can be achieved is size_beam/(S/N), where size_beam is the beam size for our ALMA observations. Therefore, the faintest sources detected, those at S/N = 4.1, have a position accuracy of 024.

Zoom In Zoom Out Reset image size

The sizes of the 2 mm continuum and CO(5–4) emitting regions are estimated from a 2D fit of the continuum and CO(5–4) moment zero maps. The fit is obtained using the task imfit within CASA, selecting a rectangular region of 4'' × 4'' around each source. We fit a 2D Gaussian profile with the centroid, major and minor axes, position angle, and integrated flux as free parameters.

All the observed sizes from the fits are in the range 1.1×–1.8× of the synthesized beams, with all the observed CO(5–4) emitting regions on scales ≳1.3 of the beam size. Given the high S/N of the detections, all the CO(5–4) emissions are therefore resolved (e.g., Decarli et al. 2018). The effective radius of each CO(5–4) emitting region, defined as the major semiaxis, is found to be R_CO(5−4) = 3.0 ± 0.8, 5.1 ± 0.8, 5.4 ± 0.9, 3.6 ± 0.7 kpc for QSO, QSO2, AGN1, and LAE1, respectively. ²¹ Hence QSO has likely the most compact molecular reservoir down to the current depth of the observations. QSO2 has also the continuum resolved on comparable sizes to its CO(5–4) emission, R_2mm = 5.3 ± 1.0 kpc. All these size measurements are reported in Table 3, and in Section 8 we discuss them in comparison with values from the literature. We stress that we also tested several Sérsic profile fits, finding that Gaussian profiles, i.e., Sérsic profiles with n ≈ 0.5, are preferred at the current spatial resolution.

In addition, there are hints for resolved kinematics within each source. Indeed, there are symmetric blue and redshifts within the first-moment maps of QSO2, AGN1, and LAE1 at the location of the highest S/N in the zero-moment maps. Similar kinematic features have been reported in other high-z quasars and have been interpreted as rotation (e.g., Bischetti et al. 2021). The major axis (same as the line of nodes) of these tentative rotation-like features was constrained by fitting a simple rotational curve to each object. Specifically, we perform chi-squared minimization to estimate the position angle of the major axis, defined as the angle taken in the counterclockwise direction between the north direction in the sky and the major axis of the galaxy. The rotational curve was assumed to follow the simplest function, the arctan (Courteau 1997), which is flexible enough to reasonably describe z ≳ 1 rotating galaxies (e.g., Miller et al. 2011; Swinbank et al. 2012). We follow the procedures described in Chen et al. (2017) to project the 1D arctan function to 2D, and run a Markov Chain Monte Carlo with the EMCEE Python package (Foreman-Mackey et al. 2013) to fit the velocity maps and obtain posterior probability distributions. Because the asymptotic velocity and inclination angle are essentially degenerate for our data quality, we treat these two as a single parameter in this S/N weighted fit. As a result, the position angles of the major axis are ${14}_{-5}^{+4}$, ${50}_{-22}^{+22}$, and ${120}_{-40}^{+9}$ degrees for QSO2, AGN1, and LAE1, respectively. We highlight the obtained line of nodes (magenta) in Figure 6 and list in Table 4 the angles ϕ defining these directions in the reference frame east of north, together with their uncertainties. Table 4 also lists the angles θ between the spin directions of QSO2, AGN1, and LAE1 with respect to the direction to QSO on the projected plane. The spin directions are assumed to be perpendicular to the major axis. For each companion source we found that its line of nodes (spin) is almost perpendicular (parallel) to its direction to QSO, though with large uncertainties.

Table 4. Angle Measurements (Reference Frame East of North)

ID	ϕa	θb	ϕoff c
	(deg)	(deg)	(deg)
QSO2	${14}_{-5}^{+4}$	${8}_{-5}^{+4}$	14 ± 39
AGN1	${50}_{-22}^{+22}$	${28}_{-22}^{+22}$	50 ± 76
LAE1	${120}_{-40}^{+9}$	${15}_{-40}^{+9}$	96 ± 28

Notes.

^a Major axis angle obtained from the first-moment map as described in Section 3.5. ^b Angle between the spin vector of each object and its direction to QSO. ^c Major axis angle obtained from the offset position in the zero-moment maps at negative and positive velocities, as described in Section 3.5.

Download table as:  ASCIITypeset image

To further inspect these velocity shifts, we produced zero-moment maps within several velocity ranges. Figure 7 shows an example of these maps, highlighting the location of the CO(5–4) emission at negative (blue contours), positive (red), and around zero (yellow) velocities with respect to the 2 mm continuum emission. This test further confirms our previous analysis. In particular, QSO does not show a spatial offset between the emission at positive and negative velocities. Therefore, the three peaks in its integrated CO(5–4) spectrum do not seem to be associated with three distinct components at this spatial resolution. On the other hand, we find significant offsets between emission at negative and positive velocities of 04 ± 02, 04 ± 02, and 05 ± 01 for QSO2, AGN1, and LAE1, respectively. As a sanity check, the directions of these shifts are consistent with the major axis computed by fitting the moment maps. However, they have larger uncertainties as they are not based on a fit to the full data information. Specifically, the angles between the negative and positive peaks are found to be 14° ± 39°, 50° ± 76°, and 96° ± 28° for QSO2, AGN1, and LAE1, respectively. We also list these angles in Table 4.

Zoom In Zoom Out Reset image size

Finally, Figure 6 compares the location of the centroid of the Lyα emission of each source with the millimeter observations. While the QSO and QSO2 locations are consistent at different wavelengths, AGN1 and LAE1 show measurable shifts between the 2 mm continuum (or CO(5–4) emission) and the Lyα emission. We estimate these to be 15 (or ∼11.4 kpc) and 09 (or ∼6.8 kpc), respectively, for AGN1 and LAE1. These shifts are significant for both the ALMA and MUSE (066 seeing) observations. We stress that we have verified the astrometric calibration of the MUSE observations presented in Arrigoni Battaia et al. (2018a) against the two available sources (one is QSO) in the GAIA DR2 catalog (Gaia Collaboration et al. 2018) within the observations' field of view. We found agreement between the astrometric calibration done using the Sloan Digital Sky Survey DR12 catalog (Alam et al. 2015) and the few GAIA sources, confirming the precision of astrometry of about one pixel of the integral field unit data (∼02). We further note that the offsets of AGN1 and LAE1 are in opposite directions with respect to each other, which would require a rather weird distortion pattern throughout the data to cancel them out. Therefore, we are confident that the aforementioned shifts between millimeter observations and Lyα are real. We discuss these shifts in Section 6.

The dust and stellar masses, as well as the SFRs, are estimated by fitting the spectral energy distribution (SED) of each source, as is usually done in the literature (e.g., Calistro Rivera et al. 2016; Circosta et al. 2018). The SEDs are built using the data described in the previous sections, together with the information at 3.4, 4.6, 12, 22 μm from the All Wide-field Infrared Survey Explorer (AllWISE) source Catalog, ²² and at 1.4 GHz from Very Large Array (VLA) Far-Infrared and Submillimetre Telescope (FIRST; Becker et al. 1994). These additional data points are listed in Appendix E (Table 5). The data set therefore includes at least 11 data points for each source. However, given their faintness, QSO2, AGN1, and LAE1 are detected in only 5 bands out of the 11 currently available.

We rely on the SED fitting code CIGALE (v2018.0; Boquien et al. 2019), which covers the full range of the current data sets, from rest-frame UV to radio emission. CIGALE fits simultaneously all this wavelength range imposing energy balance between the UV and the infrared (IR) emission (reprocessed dust emission), while decomposing the SED into different physically motivated components. This energy balance is critical for getting meaningful stellar masses with few data points. For our specific case, we select (i) an AGN component (accretion disk plus hot dust emission; Fritz et al. 2006), (ii) dust emission from star-forming regions (Draine & Li 2007; Draine et al. 2014), (iii) radio synchrotron emission, and (iv) stellar emission from the host galaxy, which is modeled by an exponentially declining star formation history (SFH), the simple stellar population models of Bruzual & Charlot (2003), a Chabrier initial mass function (Chabrier 2003), and a modified starburst attenuation law (based on Calzetti et al. 2000 and Leitherer et al. 2002). Details on these specific models and a comparison with other models implemented in CIGALE are discussed in Boquien et al. (2019). The parameters available for each model using the code's notation and the ranges explored by our fit are:

• 1.  
  AGN emission: this model has seven parameters, five of which are left free to explore all the values allowed by CIGALE. The remaining parameters are the AGN fraction (fracAGN; defined as the ratio of the AGN luminosity to the sum of the AGN and dust luminosities) and the angle between the equatorial axis and the line of sight (psy). We let fracAGN vary between 0 and 1 in steps of 0.05. psy is allowed to vary between 0.001 and 40.100 for type-2 AGN, and between 50.100 and 89.990 for type-1 AGN. In the case of LAE1, for which no AGN signature is present in the MUSE data (Arrigoni Battaia et al. 2018a), we neglect the AGN component during the fit. However, as its data points are very similar to those of AGN1, for completeness, we report in Appendix F the fit including a type-2 AGN component.
• 2.  
  dust emission: this model has four parameters (mass fraction of polycyclic aromatic hydrocarbon, minimum radiation field, power-law slope, fraction of dust illuminated), which are left free to explore all the values allowed by CIGALE.
• 3.  
  synchrotron emission: this model has two parameters, the value of the far-IR (FIR)-to-radio coefficient (Helou et al. 1985) and the slope of the synchrotron power law. Given the absence of tight constraints in the radio for any of the sources we fixed the slope to −1.0 (as observed in sources within other ELANe; e.g., Decarli et al. 2021) and the ratio to an arbitrary value satisfying the VLA FIRST upper limits. This portion of the SED has to be considered simply as illustrative.
• 4.  
  stellar emission: the SFH is modeled with two parameters, the age and e-folding time τ. The age is allowed to vary between 0.1 Gyr and 2 Gyr (about the age of the universe at z = 3.164) in steps of 0.1 Gyr, while τ can vary between 0.1 and 10 Gyr in steps of 0.1 Gyr. The attenuation model is set up so that the final E(B − V) is between 0 and 3. All the other eight parameters are kept to the default values.

In addition, for high-redshift sources, CIGALE applies a correction to rest-frame UV data for the attenuation from the foreground intergalactic medium following Meiksin (2006).

In the continuum fit we do not include the nebular emission component, for which CIGALE has built-in templates. Indeed, we find that the nebular Ly α line emission from some of these sources is displaced with respect to the continuum (e.g., Section 3.5). The best-fit models obtained following this procedure using the pdf_analysis module in CIGALE are shown in Figure 8, together with their χ² values and the observed data points. The likelihood-weighted output dust and stellar masses, and the SFRs together with their likelihood-weighted uncertainties are listed in Table 3. We stress that these uncertainties do not include systematic errors due to the models used, a priori assumptions on the nature of the sources (which are key in determining the use or not of an AGN component given the small number of data points for the faint sources), and the discrete coverage of the parameter space.

Zoom In Zoom Out Reset image size

Specifically, we find dust masses in the range M_dust = 1.5–16 × 10⁸ M_⊙, with LAE1 being the most dust-rich object in the system. As an additional check, we computed the dust masses using a modified blackbody model, assuming (i) a dust temperature T_dust = 40 K in the range for high-redshift quasars (e.g., Carilli & Walter 2013), (ii) a dust opacity at 850 μm of κ_d = 0.43 cm² g⁻¹ (Li & Draine 2001), and (iii) a fixed dust emissivity power-law spectral index β derived from the SMA and ALMA continuum. We find β = 2.4 ± 0.1, 2.9 ± 0.2, 2.3 ± 0.2, 1.8 ± 0.2, for QSO, QSO2, AGN1, and LAE1, respectively. The values for QSO, QSO2, and AGN1 are on the high side of the values usually found for high-redshift quasars (e.g., β = 1.95 ± 0.3, Priddey & McMahon 2001; β = 1.6 ± 0.1, Beelen et al. 2006) or dusty star-forming galaxies (e.g., β = 2.0 ± 0.2 , Magnelli et al. 2012), and are consistent with the value of β = 2.5 used to fit SEDs of HzRGs (Falkendal et al. 2019). The dust masses derived with this method are M_dust = (9 ± 3) × 10⁸ M_⊙, (3 ± 2) × 10⁸ M_⊙, (5 ± 2) × 10⁸ M_⊙, and (4 ± 2) × 10⁸ M_⊙, respectively, for QSO, QSO2, AGN1, and LAE1. Hence, they agree with the CIGALE output (LAE1 within 2σ). We stress that in this calculation we assumed the dust to be optically thin in all four sources. Given the current source sizes estimated at 2 mm, β, and the assumed dust opacity, this assumption is confirmed except for QSO, for which τ_dust ≥ 1 at λ < 145 μm. Nevertheless, we do not have any information on the source sizes at λ < 2 mm, and we decided to quote for QSO the M_dust value for the optically thin case for ease of comparison with the other sources. An optically thick calculation for QSO would give lower dust masses (e.g., Spilker et al. 2016; Cortzen et al. 2020) in better agreement with the CIGALE fit. The obtained dust masses are (i) in the range usually derived for high-redshift quasars hosts from z ∼ 2 up to z ∼ 7 (e.g., Weiß et al. 2003; Schumacher et al. 2012; Venemans et al. 2016), (ii) within the typical range for high-redshift dusty, star-forming galaxies (e.g., Casey et al. 2014; Dudzevičiūtė et al. 2020), and (iii) similar to what is reported for HzRGs (few times 10⁸ M_⊙ assuming T ∼ 50 K; e.g., Archibald et al. 2001).

The stellar masses are found to be in the range M_star = 4.0 × 10⁹ M_⊙ to 1.2 × 10¹¹ M_⊙, with the host of QSO being the most massive galaxy in this system, as expected (Arrigoni Battaia et al. 2018a). The other associated galaxies are about 10× less massive (Table 3). Given that QSO greatly outshine its host galaxy, its stellar mass derived with CIGALE has to be taken with caution. However, we show in Section 4.3 that the dynamical mass derived from CO(5–4) is consistent with the presence of such a massive galaxy. The obtained stellar masses agree well with literature values for quasar hosts (e.g., ∼10¹¹ M_⊙, Decarli et al. 2010) and the more massive and IR-detected LAEs (e.g., 10⁹–10^10.5 M_⊙; Ono et al. 2010) found at similar redshifts. The QSO host has a stellar mass consistent with the median stellar mass for Spitzer 1 < z < 5.2 HzRGs (∼10¹¹ M_⊙; De Breuck et al. 2010).

With these estimates, the dust-to-stellar mass ratios are M_dust/M_star = 0.006 ± 0.004, 0.010 ± 0.004, 0.05 ± 0.02, 0.4 ± 0.1 for QSO, QSO2, AGN1, and LAE1, respectively. These values are in agreement within uncertainties with observations of main-sequence and starburst galaxies at these redshifts (in the range 0.001 – 0.1; e.g., da Cunha et al. 2015; Donevski et al. 2020), except LAE1, which has a larger value. Given the similarity with the SED of AGN1, this tension could be solved by including an AGN component in the SED fit for LAE1 (M_dust/M_star = 0.2 ± 0.1; see Appendix F). The current data set does not show a clear evidence for AGN activity in LAE1, but it could be obscured. Deep infrared and X-ray observations are therefore needed to verify the nature of this source.

Further, CIGALE determined the instantaneous SFR for each source, indicating SFR = 521 ± 74, 183 ± 49, 104 ± 21, 50 ± 11 M_⊙ yr⁻¹ for QSO, QSO2, AGN1, LAE1, respectively. These SFRs are in line with values from the literature. In particular, the SFR in the QSO host agrees well with values usually found in z ∼ 2 − 3 quasars (e.g., Harris et al. 2016). We also computed the SFR from the total infrared (IR) emission using the rest-frame wavelength range 8 − 1000 μm for the obtained SED. For this purpose, we used the relation (SFR _IR/[M_⊙ yr⁻¹]) = 3.88 × 10⁻⁴⁴(L_IR/[erg s⁻¹]) of Murphy et al. (2011), usually employed for high-redshift quasars (e.g., Venemans et al. 2017). This relation has been computed using Starburst99 (Leitherer et al. 1999) with a fixed SFR and a Kroupa initial mass function (Kroupa 2001), and assumes that the entire Balmer continuum is absorbed and reradiated by dust in the optically thin regime. Applying this relation only to the IR luminosity L_IR computed excluding the AGN contribution (see Table 3), results in SFR_IR = 700 ± 30, 1087 ± 283, 447 ± 149, 179 ± 30 M_⊙ yr⁻¹ for QSO, QSO2, AGN1, LAE1, respectively. These values are higher than the instantaneous SFRs, reflecting the longer timescales probed by the IR tracers.

Further, in Figure 9 we show the locations of QSO, QSO2, AGN1, and LAE1 in the L_IR versus ${L}_{\mathrm{CO}(5-4)}^{{\prime} }$ plot in comparison to submillimeter galaxies (SMGs) and high-redshift quasars. The CO(5–4) line is known to be linked to L_IR through a relation of the form log ${L}_{\mathrm{IR}}=\alpha \mathrm{log}{L}_{\mathrm{CO}(5-4)}^{{\prime} }+\beta$ from low to high redshift (e.g., Greve et al. 2014; Daddi et al. 2015; dotted and dashed lines). Our sources are consistent with the scatter of known high-redshift sources (e.g., Valentino et al. 2020; blue line with intrinsic scatter), ensuring that the AGN-corrected L_IR obtained from the SED fit are reasonable.

Zoom In Zoom Out Reset image size

Summarizing, we find values for dust and stellar masses, and SFRs within the scatter of observations reported in the literature. Follow-up deep observations in the near-IR (NIR) and millimeter regimes are needed to better constrain the SED of the targeted sources (especially for QSO2, AGN1, and LAE1), and therefore verify their nature and the current estimates.

It is common practice to obtain the molecular mass through the equation ${M}_{\mathrm{gas}}=\alpha {L}_{\mathrm{CO}(1-0)}^{{\prime} }$, where α is the CO luminosity to gas mass conversion factor, and ${L}_{\mathrm{CO}(1-0)}^{{\prime} }$ is the CO(1–0) luminosity in units of K km s⁻¹ pc² (e.g., Carilli & Walter 2013; Aravena et al. 2019). To use this equation, we assume α_CO = 0.8 M_⊙ ${\left({\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{pc}}^{2}\right)}^{-1}$. This value has been derived for local ultraluminous infrared galaxies (ULIRGs; e.g., Downes & Solomon 1998), and it is commonly used to calculate molecular gas masses in high-redshift quasars (z ∼ 2–7; e.g., Riechers et al. 2006; Coppin et al. 2008; Carilli & Walter 2013; Venemans et al. 2017). As only the CO(5–4) line flux is available, we have to further assume a CO spectral line energy distribution (CO SLED) to derive how strong the CO(1–0) transition is. Current statistics show that the CO SLED of high-redshift quasars peaks at high-J transition, CO(6-5) and CO(7-6), with a minimum flux ratio CO(5–4)/CO(1–0) ∼ 10 (Weiss et al. 2007; Carilli & Walter 2013). Therefore, we derive the corresponding ${L}_{\mathrm{CO}(1-0)}^{{\prime} }$ assuming this ratio. The values obtained for the three AGN (QSO, QSO2, AGN1) are considered as possible upper limits (i.e., their CO SLED could be more excited; e.g., Weiß et al. 2007), while for LAE1 it is possibly a lower limit. For completeness we list the inferred CO(1–0) luminosities in Table 3. We note that the adopted CO(5–4)/CO(1–0) corresponds to ${r}_{51}\,={L}_{\mathrm{CO}(5-4)}^{{\prime} }/{L}_{\mathrm{CO}(1-0)}^{{\prime} }=0.4$ (the range reported by Carilli & Walter (2013) is r₅₁ = 0.69 ± 0.3). This value is also consistent within 2σ with: the values reported for the spectrum obtained for z > 2 gravitationally lensed dusty star-forming galaxies by Spilker et al. (2014; r₅₁ = 0.67 ± 0.20), the median value reported for 32 z ∼ 1.2–4.1 luminous submillimeter galaxies (r₅₁ = 0.32 ± 0.05; Bothwell et al. 2013), the average value for 8 z > 2 star-forming galaxies (r₅₁ = 0.44 ± 0.11; Boogaard et al. 2020), and for the compilation of SMGs in Birkin et al. (2021 ; r₅₁ = 0.35 ± 0.08).

Using the derived CO(1–0) luminosities and the assumed α_CO, the resulting gas masses M_gas are (1.5 ± 0.1) × 10¹⁰, (7.7 ± 0.7) × 10⁹, (1.0 ± 0.1) × 10¹⁰, and (1.2 ± 0.1) × 10¹⁰ M_⊙ for QSO, QSO2, AGN1, and LAE1, respectively (see also Table 3). We stress that the errors on the molecular gas masses do not include systematics due to the uncertain α_CO, the CO excitation, and the possibility of having some CO-dark gas in these systems (e.g., Balashev et al. 2017). For example, if the conditions in the molecular gas are more similar to the Milky Way, the molecular gas masses could be up to a factor of 5 larger, i.e., α_CO ∼ 4 M_⊙ ${\left({\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{pc}}^{2}\right)}^{-1}$ (e.g., Bolatto et al. 2013). Nevertheless, the obtained molecular gas masses are within the ranges reported in the literature for high-redshift quasars and dusty star-forming galaxies (M_gas = few × 10¹⁰(α_CO/0.8) M_⊙; e.g., Bothwell et al. 2013; Carilli & Walter 2013), and are also at the low-end of masses reported for HzRGs (10¹⁰ − 10¹¹ M_⊙; e.g., Miley & De Breuck 2008).

When comparing the obtained molecular gas masses to the dust masses derived in Section 4.1, we find molecular-to-dust mass ratios ${r}_{{{\rm{H}}}_{2},\ \mathrm{dust}}$ in the range 8–51 (see Table 3), which are very low in comparison to the usually assumed gas-to-dust mass ratios for local galaxies (∼100; e.g., Draine et al. 2007; Galametz et al. 2011) and high-redshift massive star-forming galaxies (∼100; e.g., Riechers et al. 2013), even when correcting them for the fraction of gas in molecular form (∼80%; e.g., Riechers et al. 2013). Interestingly, these values are more similar to what is seen in SMGs (e.g., ${r}_{{{\rm{H}}}_{2},\ \mathrm{dust}}={28}_{-6}^{+7};$ Kovács et al. 2006). Therefore, our measurement could be due to a real molecular gas deficiency or efficient dust absorption (i.e., larger dust opacities) in these sources, or could be related to the assumptions made to determine the molecular gas masses. Assuming a Milky Way value of α_CO ∼ 4 M_⊙ ${\left({\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{pc}}^{2}\right)}^{-1}$, the molecular-to-dust mass ratios increase, but still two sources (AGN1 and LAE1) have values reflecting a depletion of molecular gas. Less excited CO ladders would then be needed to increase our molecular mass estimates and thus alleviate the tension for these sources. It is clear that our measurements need to be refined with follow-up observations of additional CO transitions (especially at lower J; ALMA, Northern Extended Millimeter Array, J-VLA) or other molecular gas tracers (e.g., [C i]), to at least remove the uncertainties on the CO excitation ladder. For completeness, the ${r}_{{{\rm{H}}}_{2},\ \mathrm{dust}}$ values are also listed in Table 3.

In this section we outline rough estimates of the dynamical masses for QSO, QSO2, AGN1, and LAE1, using the kinematics and sizes of the CO(5–4) emitting region. In turn, these dynamical masses can be compared to the galaxies' mass budget derived in the previous sections. As rotation-like signatures are present in most of the CO(5–4) maps (Section 3.5), the bulk of the molecular gas is assumed to be in a disk with an inclination i. This approach is common practice in the study of unresolved line tracers of molecular gas associated with high-redshift quasars (e.g., Decarli et al. 2018). In this framework, the dynamical mass can be obtained as ${M}_{\mathrm{dyn}}={G}^{-1}{R}_{\mathrm{CO}(5-4)}{\left(\mathrm{FWHM}/\sin (i)\right)}^{2}$ (Willott et al. 2015), where G is the gravitational constant, R_CO(5−4) is the size of the CO(5–4) emitting region, and FWHM is its line width. We omitted the frequently used 0.75 factor to scale the line FWHM to the width of the line at 20% because the integrated CO(5–4) line shape is not a simple Gaussian.

As the quality of our data does not allow an estimate of the inclination angle (Section 3.5), we assume the mean inclination angle for randomly oriented disks $\langle \sin (i)\rangle =\pi /4$ (e.g., Law et al. 2009), obtaining M_dyn = (1.9 ± 0.8) × 10¹¹, (8.5 ± 6.7) ×10¹⁰, (3.9 ± 2.9) × 10¹¹, (8.2 ± 5.2) × 10¹⁰ M_⊙ for QSO, QSO2, AGN1, and LAE1, respectively. These dynamical masses are consistent within their large uncertainties with the sum of the different galaxy components, i.e., molecular, stellar, dust, and DM. We computed the DM mass within R_CO(5−4) assuming a Navarro–Frenk–White (NFW; Navarro et al. 1997) profile and the concentration–halo mass relation of Dutton & Macciò (2014). We found DM masses of $({2.4}_{-1.1}^{+1.3})\times {10}^{10}$, (2.6 ± 0.2) × 10¹⁰, $({2.8}_{-0.4}^{+0.1})\times {10}^{10}$, (1.2 ± 0.1) × 10¹⁰ M_⊙ for QSO, QSO2, AGN1, and LAE1, respectively.

For QSO2, AGN1, and LAE1, the dynamical masses including some pressure support can be further assessed, using the observed asymptotic rotational velocities ${v}_{\mathrm{rot}}^{\mathrm{obs}}$ obtained from the 2D fit of their first-moment maps described in Section 3.5, and the velocity dispersion σ within their effective radii in their second-moment maps. In this framework, ${M}_{\mathrm{dyn}}=2{R}_{\mathrm{CO}(5-4)}({v}_{\mathrm{rot}}^{2}+{\sigma }^{2})/G$ (e.g., Smit et al. 2018), where ${v}_{\mathrm{rot}}={v}_{\mathrm{rot}}^{\mathrm{obs}}/\sin (i)$, assuming again $\langle \sin (i)\rangle =\pi /4$. Using the computed values of ${v}_{\mathrm{rot}}^{\mathrm{obs}}={170}_{-60}^{+30},{190}_{-60}^{+60},{190}_{-40}^{+50}$ km s⁻¹ and σ = 119 ± 89, 259 ± 153, 156 ± 138 km s⁻¹, we obtain ${M}_{\mathrm{dyn}}={1.4}_{-0.7}^{+0.6}\times {10}^{11},{3.1}_{-2.1}^{+2.1}\times {10}^{11},{1.4}_{-0.8}^{+0.8}\times {10}^{11}$ M_⊙ for QSO2, AGN1, and LAE1, respectively. These masses are consistent with the previously determined values.

Notwithstanding the aforementioned fair agreement between the estimated dynamical masses and the mass budget for each galaxy, we note that the obtained dynamical masses are usually providing larger mass estimates, especially for AGN1. This tentative mismatch can however be evidence of turbulence injection in the molecular reservoir due to different physical processes expected in galaxy evolution: infall of gas at velocities of hundreds of km s⁻¹, stellar and AGN feedback. In other words, turbulence due to these processes could explain the large velocity dispersions seen in the four CO(5–4) detected objects (Figure 6). Higher-resolution observations, exploiting the ALMA longest baselines, are required to firmly assess the dynamical masses of these sources.

In this section we present our estimate for the mass budget of the whole system and compare that to the cosmic baryon fraction. To compute the DM and baryonic (stars, molecular gas, dust, atomic gas) components, we rely on the previously obtained values and on several additional assumptions which are needed to overcome the limitations of the current observations.

We assume that the sources detected by ALMA, i.e., QSO, QSO2, AGN1, LAE1, are the most massive objects in this system, and neglect the contributions from additional sources. It will be clear that additional satellite masses are well accommodated within the final error estimates.

4.4.1. The Dark Matter Component

4.4.2. The Baryonic Components

For the baryon budget, we proceed by simply adding up the masses of each component for QSO, QSO2, AGN1, and LAE1 and propagating their errors, finding ${M}_{\mathrm{star}}^{\mathrm{total}}=(1.5\pm 0.8)\times {10}^{11}$ M_⊙, ${M}_{\mathrm{dust}}^{\mathrm{total}}=(3.0\pm 0.4)\times {10}^{9}$ M_⊙, ${M}_{{{\rm{H}}}_{2}}^{\mathrm{total}}=({4.5}_{-0.2}^{+17.9})\,\times \,{10}^{10}$ M_⊙, for the total stellar, dust, and molecular masses. In the error budget for the molecular mass we include the large uncertainty (a factor of 5) on α_CO. This large uncertainty should also include the possibility of molecular gas extending on scales larger than the body of galaxies as seen, for example, in HzRGs environments (e.g., Emonts et al. 2016; see Section 8 for discussion). From the mass budget we then miss the atomic gas components at different temperatures, i.e., cold (∼100 K), cool (∼10⁴ K), and warm/hot (>10⁵ K).

We can predict the amount of cold atomic gas by assuming that the interstellar-medium molecular gas fraction is ∼80% at high redshift (e.g., Riechers et al. 2013), and in turn that the cold atomic gas fraction is therefore f_cold ∼ 20 %. This is also consistent with current estimates for such massive halos at z ∼ 3 from semiempirical models of galaxy evolution (e.g., Popping et al. 2015). We therefore include in the budget a total cold atomic mass of ${M}_{\mathrm{HI}}^{\mathrm{cold}}=({0.9}_{-0.1}^{+3.6})\times {10}^{10}$ M_⊙. This prediction could be tested by targeting the [C ii] emission at 158 μm with, e.g., ALMA (e.g., Fujimoto et al. 2020).

To derive the total cool gas mass, we can instead rely on the large-scale Lyα emission detected with VLT/MUSE in Arrigoni Battaia et al. (2018a). Given that the projected maximum distance of the Lyα emission is comparable with the obtained virial radii, we assume that all the Lyα nebula sits within the halo. This is also in agreement with the discussion in Arrigoni Battaia et al. (2018a), who argued that the Lyα emission traces the motions of substructures accreting within the bright quasar massive halo. We then assume that the visible Lyα emitting gas is the densest cool gas in the halo, and thus the one contributing to most of its mass. The cool gas mass can then be estimated as ${M}_{\mathrm{cool}}^{\mathrm{total}}=A\,({m}_{{\rm{p}}}/X)\,{f}_{C}\,{N}_{{\rm{H}}}$ (Hennawi & Prochaska 2013), where A is the projected area on the sky covered by the ELAN in cm² (609.36 arcsec² within the 2σ isophote in Arrigoni Battaia et al. 2018a), m_p is the proton mass, X = 0.76 is the hydrogen mass fraction (e.g., Pagel 1997), f_C is the cool gas covering factor within the ELAN, and N_H is the total hydrogen column density of the emitting gas. We assume (i) f_C = 1 as it has been shown that the observed morphology of extended Lyα nebulae can be reproduced if f_C ≳ 0.5 (Arrigoni Battaia et al. 2015b), and (ii) a constant $\mathrm{log}({N}_{H}/{\mathrm{cm}}^{-2})=20.5\pm 1.0$, which is the median value found by Lau et al. (2016) for optically thick absorbers in z ∼ 2 quasar halos (see their Figure 15). For this latter value we assume an error, which encompasses the large uncertainties in some of those authors data points. Inserting these values in the aforementioned relation gives ${M}_{\mathrm{cool}}^{\mathrm{total}}={1.2}_{-1.0}^{+10.5}\times {10}^{11}$ M_⊙. We stress that this calculation neglects additional cool gas within the halo not emitting Lyα above the sensitivity of the observations in Arrigoni Battaia et al. (2018a). However, the current area of the nebula covers 42% or 25% of the projected halo for the Moster et al. (2018) or the Tempel et al. (2014) calculation, respectively. Even if we assume the full halo to be covered by Lyα emission the total mass would increase by a factor of 2.4 or 4, respectively, thus falling within the errors of the previous measurement. Interestingly, the obtained ${M}_{\mathrm{cool}}^{\mathrm{total}}$ agrees well with cool gas masses reported for z ∼ 2–3 quasars halos (${M}_{\mathrm{cool}}^{\mathrm{total}}\gt {10}^{10}$ M_⊙; e.g., Prochaska et al. 2013; Lau et al. 2016), and reported for other ELANe (1.0 × 10¹¹ ${M}_{\odot }\,\lt {M}_{\mathrm{cool}}^{\mathrm{total}}\lt 6.5\times {10}^{11}$ M_⊙; Hennawi et al. 2015).

Figure 10 summarizes the discussed baryonic components as fractions of the total mass of the system, which has been derived by assuming a halo baryon fraction equal to the cosmic value (15.6%). It is clear that the stars, dust, molecular, cold and cool atomic components add up to a small fraction of the cosmic value, with 21% or 10% of the baryons in these constituents depending on the DM mass considered, Moster et al. (2018) or Tempel et al. (2014), respectively. These values, though uncertain, are lower than the estimated value reported for z ∼ 2 quasars (56%; Lau et al. 2016). We can easily explain these differences with the larger halo masses derived in this work with respect to the assumed halo mass in Lau et al. (2016 ; M_DM = 10^12.5 M_⊙). In other words, we find similar baryonic masses but in a halo which is 2.7× or 5.9× more massive. If all the quasars in Lau et al. (2016) inhabit DM halos as massive as the one of this ELAN, they would have a similar baryon budget.

Zoom In Zoom Out Reset image size

As expected from galaxy formation theories (e.g., Dekel & Birnboim 2006), our analysis suggests that the rest of the baryonic mass is in a warm/hot phase which permeates the halo of this massive system. Assuming a halo baryon fraction equal to the cosmic value, this warm/hot phase would represent a reservoir as massive as ${M}_{\mathrm{warm}-\mathrm{hot}}^{\mathrm{total}}=({1.3}_{-0.2}^{+1.2})\times {10}^{12}$ M_⊙ or ${M}_{\mathrm{warm}-\mathrm{hot}}^{\mathrm{total}}=({3.1}_{-0.1}^{+1.1})\times {10}^{12}$ M_⊙, for Moster et al. (2018) and Tempel et al. (2014) DM calculations, respectively (Figure 10). The warm/hot phase together with the cool phase would then represent 87% or 94% of the baryon fraction.

We can further gain some intuition on the halo gas physical properties by assuming the cool and warm/hot phases coexist in pressure equilibrium. This assumption is likely not valid in turbulent massive halos (e.g., Nelson et al. 2020), but it is useful as a first-order approximation. We can therefore derive the physical properties of the two phases, namely, the temperature (T_cool, T_warm/hot), volume density (${n}_{{\rm{H}}}^{\mathrm{cool}}$, ${n}_{{\rm{H}}}^{\mathrm{warm}-\mathrm{hot}}$), and volume filling factors (${f}_{{\rm{V}}}^{\mathrm{cool}}$, ${f}_{{\rm{V}}}^{\mathrm{warm}-\mathrm{hot}}$). To do so, the following three relations have to be considered simultaneously: (i) the Lyα surface brightness in an optically thin scenario ${\mathrm{SB}}_{\mathrm{Ly}\alpha }\propto {\alpha }_{{\rm{A}}}({T}_{\mathrm{cool}}){N}_{{\rm{H}}}^{\mathrm{cool}}{n}_{{\rm{H}}}^{\mathrm{cool}}{f}_{{\rm{C}}}$ (Hennawi & Prochaska 2013), where α_A(T_cool) is the temperature-dependent coefficient for case A recombinations (e.g., Hui & Gnedin 1997); (ii) the pressure balance ${n}_{{\rm{H}}}^{\mathrm{cool}}{T}_{\mathrm{cool}}={n}_{{\rm{H}}}^{\mathrm{warm}-\mathrm{hot}}{T}_{\mathrm{warm}-\mathrm{hot}};$ (iii) the mass ratio of the two phases ${M}_{\mathrm{warm}-\mathrm{hot}}/{M}_{\mathrm{cool}}=({n}_{{\rm{H}}}^{\mathrm{cool}}/{n}_{{\rm{H}}}^{\mathrm{warm}-\mathrm{hot}})({f}_{{\rm{V}}}^{\mathrm{warm}-\mathrm{hot}}/{f}_{{\rm{V}}}^{\mathrm{cool}})$. Using the observed average SB_Lyα = 6.08 × 10⁻¹⁸ erg s⁻¹ cm⁻² arcsec⁻², T_warm/hot = T_virial = GMm_p/(3k_B R_vir) (e.g., White & Rees 1978), and the mass ratio obtained by assuming a halo baryon fraction equal to the cosmic value, we find T_cool = 10^4.4 (10^4.7) K, ${n}_{{\rm{H}}}^{\mathrm{cool}}=3.1\ (3.1)$ cm⁻³, ${n}_{{\rm{H}}}^{\mathrm{warm}-\mathrm{hot}}\,={10}^{-2.3}\ ({10}^{-2})$ cm⁻³, ${f}_{{\rm{V}}}^{\mathrm{cool}}=2.6\times {10}^{-4}\ (1.2\times {10}^{-4})$, where we quote in brackets the value for the DM calculation following the formalism in Tempel et al. (2014). A scale length for the structures in the cool gas responsible for the Lyα emission can then be computed as ${l}_{\mathrm{cool}}={N}_{{\rm{H}}}^{\mathrm{cool}}/{n}_{{\rm{H}}}^{\mathrm{cool}}=56\ (33)$ pc. This simple calculation agrees with previously reported properties for cool dense gas in ELANe (Cantalupo et al. 2014; Arrigoni Battaia et al. 2015a; Hennawi et al. 2015; see also discussion in Pezzulli & Cantalupo 2019).

Several recent works have studied the survival of cool clouds against hydrodynamic instabilities while moving throughout the hot halo with velocities of the order of few hundreds kilometers per second (e.g., Gronke & Oh 2018; Kanjilal et al. 2021). In particular, it has been shown that if the cool dense gas falls out of pressure balance (e.g., due to radiative processes), it could shatter in smaller structures (McCourt et al. 2018) and entrain in the warm/hot medium without being destroyed by Kelvin–Helmholtz instabilites if the original cloud sizes are larger than a critical scale (Equation (2) in Gronke & Oh 2018). The gas properties we found (scale length, temperature, density), if translated to cloud properties, fulfill this survival criterion and, given the density contrast, might lead to those clouds whose fragments are able to coagulate into larger cloudlets and therefore survive within the harsh environment of a quasar hot halo (Gronke & Oh 2020). Therefore, these small-scale processes could be the reason why there is enough dense gas resulting in the bright ELAN emission (see also Mandelker et al. 2020).

The aforementioned pressure balance calculation assumes that all Lyα emission is due to recombination, but as we will show in Section 6 there is evidence for resonant scattering at least on scales of ∼15 kpc (or 2'') around compact sources. For this reason, we recompute the aforementioned quantities by assuming that only a fraction of the observed SB_Lyα is due to recombination. As expected, lowering the SB_Lyα signal due to recombination decreases ${n}_{{\rm{H}}}^{\mathrm{cool}}$ (and therefore increases ${f}_{{\rm{V}}}^{\mathrm{cool}}$), while increasing T_cool. For example, assuming a recombination fraction of 50%, we find T_cool = 10^4.6 (10^4.9) K, ${n}_{{\rm{H}}}^{\mathrm{cool}}\,=1.9\ (1.9)$ cm⁻³, ${f}_{{\rm{V}}}^{\mathrm{cool}}\,=\,4.1\,\times \,{10}^{-4}\ (1.9\,\times \,{10}^{-4})$, l_cool = 89 (53) pc, where the values in brackets correspond to the DM calculation following the formalism in Tempel et al. (2014). The warm/hot densities are not affected because they are linked to the mass ratio of the two phases. Also in this case, the aforementioned cloud survival scenario holds.

We further conduct this calculation assuming the possibility that the halo baryon content is only a small fraction of the cosmic value. Indeed, current cosmological hydrodynamical simulations of structure formation implement strong feedback (supernova and AGN) recipes to match the observed properties of galaxies (e.g., Schaye et al. 2015; Springel et al. 2018). These feedbacks are able to eject large fractions of baryons from the halos of interest here, making the halo baryon fraction only ∼1/3 of the cosmic value (e.g., Davé 2009; Wright et al. 2020). In this framework, the hot reservoir within the virial radius will decrease, resulting in M_hot/M_cool = 2.9 or 7.6 for the Moster et al. (2018) and Tempel et al. (2014) calculations, respectively. In other words, the hot phase is 74% or 88% of the halo baryons, which would be in agreement with recent results from cosmological simulations (∼80%; e.g., Gabor & Davé 2015; Correa et al. 2018). Assuming 50% Lyα emission from recombination, we then find T_cool = 10^4.2 (10^4.5) K, ${n}_{{\rm{H}}}^{\mathrm{cool}}=0.7\ (1.2)$ cm⁻³, ${n}_{{\rm{H}}}^{\mathrm{warm}-\mathrm{hot}}={10}^{-2.9}\ ({10}^{-2.6})$ cm⁻³, ${f}_{{\rm{V}}}^{\mathrm{cool}}=6.5\times {10}^{-4}\ (2.9\times {10}^{-4})$, l_cool = 141 (83) pc, where again the values in brackets correspond to the DM calculation following the formalism in Tempel et al. (2014). Also in this case, the aforementioned cloud survival scenario holds. Direct observations of the warm/hot phase are therefore crucial for constraining the warm/hot fraction, and ultimately galaxy formation models.

Here we briefly discuss the fate of this ELAN system in light of the ensemble of our findings presented in previous sections. We first focus on the halo mass. Following the expected evolution of DM halos in cosmological simulations (e.g., van den Bosch et al. 2014), the obtained M_DM values at z ∼ 3 would then result in halo masses M_DM > 10¹⁴ M_⊙ at z = 0. This result is a first evidence that this ELAN could be considered as the nursery of a local elliptical galaxy. Using then the molecular depletion timescale t_depl = M_H2/SFR (e.g., Tacconi et al. 2020 and references therein), we can assess how long the current star formation can be sustained in the ELAN system without any additional fuel from gas recycling or infall. Assuming the SFRs obtained with CIGALE for the sources within the ELAN, we find ${t}_{\mathrm{depl}}^{\mathrm{QSO}}=29\mbox{--}144$ Myr, ${t}_{\mathrm{depl}}^{\mathrm{QSO}2}=42\mbox{--}210$ Myr, ${t}_{\mathrm{depl}}^{\mathrm{AGN}1}\,=96\mbox{--}480$ Myr, ${t}_{\mathrm{depl}}^{\mathrm{LAE}1}=240\mbox{--}1200$ Myr, where these conservative ranges take into account the uncertainty on α_CO (i.e., a factor of 5). Without the help of recycling and infall, these objects are thus not able to sustain their current SFR for long periods, with the longest t_depl allowing to possibly reach z ∼ 2 if the merger between QSO and LAE1 does not happen by then.

To fuel the system for longer periods, a net mass infall is therefore required. To compute a rough estimate for the mass inflow rate, we assume gas infall velocities constant with the radius and given by the first-order approximation v_in = 0.9v_vir (Goerdt & Ceverino 2015), and use the Lyα emission as a mass tracer. Assuming the total cool gas mass ${M}_{\mathrm{cool}}^{\mathrm{total}}=1.2\times {10}^{11}$ M_⊙ (Section 4.4.2) and that all the mass will end up onto the central object, we then find a mass inflow rate of ${\dot{M}}_{\mathrm{in}}^{\mathrm{cool}}(z=3.164)=320$ M_⊙ yr⁻¹ for both DM halos obtained in Section 4.4.1. This fresh fuel for star formation is able to delay the depletion of the central object by a factor of 2.6 at a fixed SFR, but is not able to keep up with the star formation rate of the central object. As the gas accretion rate is expected to decrease at lower redshifts (e.g., McBride et al. 2009), a corresponding decrease in the SFR is expected with a certain delay, with the system activity almost shut down by z ∼ 1.

We can indeed compute a z = 0 stellar mass by assuming (i) the gas accretion rate as a function of z in McBride et al. (2009; i.e., ${\dot{M}}_{\mathrm{in}}\propto {\left(1+z\right)}^{2.5}$ for z ≥ 1 and ${\dot{M}}_{\mathrm{in}}\propto {\left(1+z\right)}^{1.5}$ for z < 1) normalized to ${\dot{M}}_{\mathrm{cool}}(z=3.164)$, and (ii) that all cold and cool material and satellites now in the QSO halo will end up in the central object by then. We find that the halo accretion down to z = 0 contributes 7.2 × 10¹¹ M_⊙, while the latter 2 × 10¹¹ M_⊙, if outflows are not effective in removing mass from such a massive galaxy/halo (i.e., the wind material rains back onto the galaxy; e.g., Oppenheimer & Davé 2008). This is certainly true at low redshifts, while at high-redshift winds could push material out of the halo virial radius. Nonetheless, this material may have fall back into the central or accreted satellites by z = 0. We further assume that the accreted mass is translated to stellar mass with a star formation efficiency per freefall time that scales as 1+z (Scoville et al. 2017). For this rough calculation we use an average freefall time of 10 Myr for molecular clouds (Chevance et al. 2020). By z = 0 the stellar mass due to large-scale accretion is 1.5 × 10¹¹ M_⊙. Taking into account all these assumptions and adding together (i) the stellar mass due to large-scale accretion, (ii) the mass from the satellites, and (iii) the z = 3.164 stellar mass of the QSO host, we obtain a z = 0 stellar mass M_star(z = 0) = 4.7 × 10¹¹ M_⊙, 94% of which has been built before z ∼ 1. The obtained M_star(z = 0) is similar to the masses of local giant elliptical galaxies, e.g., NGC 4365 and NGC 5044 (${M}_{\mathrm{star}}^{\mathrm{NGC}\,4365}=(4.3\,\pm 0.7)\times {10}^{11}$ M_⊙ and ${M}_{\mathrm{star}}^{\mathrm{NGC}\,5044}=(5.7\pm 0.7)\times {10}^{11}$ M_⊙; e.g., Spavone et al. 2017). Computing the expected halo mass using the relations of Moster et al. (2018), we find M_DM = 10^14.7 M_⊙, which could represent a local galaxy cluster. This calculation is in agreement with the evolutionary tracks in Chiang et al. (2013; see, e.g., their Figure 2).

This rough calculation once again points to the fact that this ELAN should be part of a large-scale protocluster, as found for other ELANe (e.g., Hennawi et al. 2015; Cai et al. 2017; see also discussion in Arrigoni Battaia et al. 2018a). Wide-field coverage is therefore needed to confirm this hypothesis and further pin down the evolution of this system.