An Extreme Protocluster of Luminous Dusty Starbursts in the Early Universe

DRC was observed with The Atacama Pathfinder Experiment (APEX) telescope's Large APEX BOlometer CAmera (LABOCA—Kreysa et al. 2003; Siringo et al. 2009) at $870\,\mu {\rm{m}}$ as part of a program aimed at looking for overdensities of dusty starbursts around ultrared DSFGs (Lewis et al. 2017). The APEX observations for DRC were carried out during 2013 October under projects M-092.F-0015-2013 (P.I. A. Weiss) and 191A-0748 (P.I. R. J. Ivison). A compact-raster scanning mode was used, whereby the telescope scans in an Archimedean spiral for ${t}_{\mathrm{int}}=35\,{\rm{s}}$ at four equally spaced raster positions in a $27^{\prime\prime} \times 27^{\prime\prime}$ grid. Each scan was approximately ${t}_{\mathrm{int}}\approx 7\,\mathrm{minutes}$ long such that each raster position was visited three times leading to a fully sampled map over the full $11^{\prime}$-diameter field of view of LABOCA. During the observations, we recorded typical precipitable water vapor (PWV) values between 0.4 and 1.3 mm, corresponding to a zenith atmospheric opacity of $\tau =0.2\mbox{--}0.4$. Finally, the flux density scale was determined to an rms accuracy of ${\sigma }_{\mathrm{calib}}\approx 7$% using observations of the primary calibrators, Uranus and Neptune, while pointing was checked every hour using nearby quasars and found to be stable to ${\sigma }_{\mathrm{point}}\approx 3^{\prime\prime}$ (rms). A total of $11.4\,\mathrm{hr}$ were spent integrating on our target, covering an area of $124\,{\mathrm{arcmin}}^{2}$. The final map was then beam-smoothed, to a resolution of $27^{\prime\prime}$. The average rms background noise is then $\sim 1.9\,\mathrm{mJy}\,{\mathrm{beam}}^{-1}$. The data were reduced using the Python-based BOlometer data Analysis Software package (BOA v4.1 —Schuller 2012), following the prescription outlined in Siringo et al. (2009) and Schuller et al. (2009). More details about the observations, data reduction and source extraction can be found in Lewis et al. (2017).

The ALMA data presented in this paper come from four different projects: 2013.1.00449.S (P.I. A. Conley), 2013.A.00014.S—a Director's Discretionary Time (DDT) proposal (P.I. R. J. Ivison), 2013.1.00001.S (P.I. R. J. Ivison), and 2016.1.01287.S (P.I. I. Oteo).

In project 2013.1.00449.S, we carried out spectral scans in the 3 mm band on a sample of eight ultrared starbursts at ${z}_{\mathrm{phot}}\gtrsim 4$ selected from HerMES (Oliver et al. 2012) and H-ATLAS surveys with the aim of measuring their redshift via multiple CO line detections (Asboth et al. 2016; Oteo et al. 2016b; Fudamoto et al. 2017; Riechers et al. 2017). The observations were taken between 2014 July 03 and August 28, in a total of five scheduling blocks (SBs) corresponding to the five tunings needed to cover most of ALMA band 3, between 84 and $114.88\,\mathrm{GHz}$. The data for each SB were reduced in CASA, following standard procedures. Imaging was carried out using natural weighting to improve sensitivity, resulting in an average rms sensitivity of $\sim 0.75\,\mathrm{mJy}\,{\mathrm{beam}}^{-1}$ in channels binned to $100\,\mathrm{km}\,{{\rm{s}}}^{-1}$. The fwhm synthesized beam ranged between $0\buildrel{\prime\prime}\over{.} 6$ and $1\buildrel{\prime\prime}\over{.} 2$ due to the different array configurations used in the different tunings.

In project 2013.A.00014.S, we observed DRC for about one hour, aiming to confirm what were thought to be two faint emission lines, ¹²CO(7–6) and [C i] (2–1), detected in DRC at around $98.4\,\mathrm{GHz}$ in data from the earlier project, 2013.1.00449.S.¹⁷ The observations were carried out on 2015 January 14, when the array was in a compact configuration. The data were reduced using the ALMA pipeline and imaged using natural weighting to improve sensitivity. The resulting beam size was $\sim 2\buildrel{\prime\prime}\over{.} 0$ fwhm and the rms sensitivity was $\sim 0.13\,\mathrm{mJy}\,{\mathrm{beam}}^{-1}$ in $100\,\mathrm{km}\,{{\rm{s}}}^{-1}$ channels.

Project 2013.1.00001.S consisted of high-spatial-resolution ($0\buildrel{\prime\prime}\over{.} 12$) continuum observations of the brightest DRC component, what we will call DRC-1, at $870\,\mu {\rm{m}}$. The aim of the observations was to study the morphology and extent of the dust emission in a subsample of the ultrared DSFGs presented in Ivison et al. (2016). Details on the observations, data calibration, and imaging can be found in Oteo et al. (2016b, 2017b). Briefly, the data were calibrated using the ALMA pipeline and imaging was done using Briggs weighting, which represents a good compromise between depth and spatial resolution. The resulting rms sensitivity is $\sim 0.1\,\mathrm{mJy}\,{\mathrm{beam}}^{-1}$ with an fwhm synthesized beam of $\sim 0\buildrel{\prime\prime}\over{.} 12$ or $\sim 830\,\mathrm{pc}$ at z = 4.002.

In project 2016.1.01287.S, we observed DRC at $2\,\mathrm{mm}$ in 14 SBs, each about $70\,\mathrm{minutes}$ long, with the aim of detecting one or more additional emission lines to determine the redshift of DRC unambiguously. Due to the smaller primary beam in band 4, a two-pointing mosaic was used, such that roughly the same area was covered as for the earlier 3 mm observations, without significant primary beam attenuation. The data were calibrated using the ALMA pipeline and imaging was done using natural weighting in the concatenated visibilities corresponding to the two pointings. The resulting rms sensitivity is $\sim 50\,\mu \mathrm{Jy}\,{\mathrm{beam}}^{-1}$ in $100\,\mathrm{km}\,{{\rm{s}}}^{-1}$ channels and $\sim 6\,\mu \mathrm{Jy}\,{\mathrm{beam}}^{-1}$ in the continuum map, with a synthesized beam size (fwhm) of $1\buildrel{\prime\prime}\over{.} 68\times 1\buildrel{\prime\prime}\over{.} 54$ in both pointings.

We used the Karl G. Jansky Very Large Array (VLA) to observe DRC, covering 27.68–28.58 GHz and 32.01–32.90 GHz using the 8-bit samples, in the CnB configuration, during 2015 January 24–25 (project VLA/14B-497; P.I. R. J. Ivison). The data were calibrated using the VLA pipeline, and the calibrated visibilities were imaged using natural weighting to improve the sensitivity. These observations were aimed at detecting low-J CO lines from DRC, but after the redshift was confirmed with the ALMA and ATCA observations, we knew that no CO lines were covered by the VLA spectral setup. Therefore, we utilize only the continuum map in this paper, which has an rms sensitivity of $\sim 6.3\,\mu \mathrm{Jy}\,{\mathrm{beam}}^{-1}$ and a synthesized beam (fwhm) of $1\buildrel{\prime\prime}\over{.} 07\times 0\buildrel{\prime\prime}\over{.} 66$.

DRC was observed on the nights of 2015 August 13, and September 5 and 7 with the Multi Unit Spectroscopic Explorer (MUSE) integral field spectrograph mounted on the European Southern Observatory's Very Large Telescope UT4 (Bacon et al. 2010; DDT program 295.A-5029, P.I. I. Oteo), with a seeing varying between 09 and 11. The $1^{\prime} \times 1^{\prime}$ MUSE field of view was centered on the protocluster core, covering the 11 components detected in our ALMA Cycle 4 observations. With the nominal wavelength range (475–930 nm), we performed a series of exposures, each $15\,\mathrm{minutes}$, for a total on-source time of $3\,\mathrm{hr}$. Between individual exposures, the spectrograph was rotated by 90° and a random dithering pattern was added. The data were reduced using the MUSE pipeline, which performed all the basic reduction steps (Bacon et al. 2015). The subtraction of the sky emission was improved in the reduced cubes using a set of custom scripts.

FORS2 broadband imaging in the I filter were carried out during $4\,{\rm{h}}$ on the night of 2014 December 17 (project 093.A-0705, P.I: R. J. Ivison) under good weather conditions with seeing $\sim 1^{\prime\prime}$. The data were reduced with the esorex pipeline following the standard procedures and the final images were astrometrically calibrated using the available VIKING z-band imaging in the field. The observations reached a $5\sigma$ limiting magnitude of $25.3\,\mathrm{mag}$.

Broadband near-IR observations in the K_s filter were carried out between 2014 July 18 and November 04 with FLAMINGOS-2 mounted in the Gemini-South telescope (project GS-2014A-Q-58, P.I. L. Dunne). DRC was observed for a total of $4.1\,{\rm{h}}$ with an average seeing of $0\buildrel{\prime\prime}\over{.} 72$. A classical dither pattern was used to remove the sky emission during the data reduction, performed with THELI (Schirmer 2013). First, the flat-field and dark correction was applied to each science frame, then the sky emission in each science frame was subtracted using a dynamical model of four images taken immediately before and after each frame. The flat-field correction was applied, cosmic rays were removed and the vignetted region of the image was masked. Finally, the images were combined using SWarp (Bertin 2010) and astrometrically corrected and flux calibrated using 2MASS (Skrutskie et al. 2006). The final image has a $3\sigma$ limiting magnitude of $25.03\,\mathrm{mag}$ in a $2^{\prime\prime}$ aperture.

DRC was observed with Spitzer/IRAC at 3.6 and 4.5 μm on 2015 September 14 (program ID: 11107; PI: Pérez-Fournon). A 36-position dither pattern with 30 s exposures per frame was used, totaling 1080 s integrations in each band. Data reduction was performed with the MOPEX package using standard procedures. Absolute astrometry was obtained relative to Gaia DR1, yielding rms accuracies of $0\buildrel{\prime\prime}\over{.} 04$ and $0\buildrel{\prime\prime}\over{.} 06$ in the 3.6 and 4.5 μm bands, respectively.

The ¹²CO(2–1) emission from DRC was observed with the Australia Telescope Compact Array (project C3185, P.I. I. Oteo) during 2017 July 18–19 in its most compact configuration, with the $64\,\mathrm{MHz}$ spectral mode. The setup of the observations was carried out using the calibrator 1921−293, which was also used for bandpass calibration. The absolute flux scale (which we estimate to have an uncertainty of $\sim 15 \%$) was determined by observing 1934−638, while 0104−408 was used as the phase calibrator. The observations were carried out in good weather conditions and covered almost full tracks, for a total on-source time of $\sim 14\,\mathrm{hr}$. The data were calibrated by using the standard techniques in MIRIAD, including manual flagging of bad data (note that antenna CA06 was used in the observations but these data were flagged after the calibration was complete due to the poor data quality with respect to the other antennas). The calibrated visibilities were then transformed into CASA format, where the cubes and continuum maps were created by using natural weighting to improve sensitivity. The resulting synthesized beam fwhm is $14\buildrel{\prime\prime}\over{.} 2\times 10\buildrel{\prime\prime}\over{.} 6$ and the rms sensitivity is $\sim 0.13\,\mathrm{mJy}\,{\mathrm{beam}}^{-1}$ in $\sim 800\,\mathrm{km}\,{{\rm{s}}}^{-1}$ wide channels.

Continuum observations with ATCA at 32.5, 9.0, and $5.5\,\mathrm{GHz}$ were carried out between 2013 October and 2015 July (project C2905, P.I. L. Dunne). The $7\,\mathrm{mm}$ observations were initially aimed at confirming the redshift of DRC via detection of ¹²CO, but no lines were detected and thus we only use the continuum maps here. The data were reduced in MIRIAD following standard procedures, then the continuum maps were created from the calibrated visibilities in CASA. At $32.5\,\mathrm{GHz}$ the continuum map reached an rms sensitivity of $\sim 7.6\,\mu \mathrm{Jy}\,{\mathrm{beam}}^{-1}$ with a synthesized beam (fwhm) of $1\buildrel{\prime\prime}\over{.} 3\times 0\buildrel{\prime\prime}\over{.} 7$. The observations at 9.0 and $5.5\,\mathrm{GHz}$ reached rms sensitivities of $\sim 6.6\,\mu \mathrm{Jy}\,{\mathrm{beam}}^{-1}$.

As an example of the difficulties sometimes encountered when trying to confirm redshifts via the detection of atomic and/or molecular lines in submillimeter and millimeter windows, it is worth mentioning that the redshift of DRC was only confirmed unambiguously after a large series of observations spanning several years. First, DRC was included in the sample of ultrared galaxies from the H-ATLAS survey (Ivison et al. 2016; Fudamoto et al. 2017) selected for spectral scans in the $3\,\mathrm{mm}$ band with ALMA. In those spectral scans of DRC, we detected only very faint emission in the center of the band, which was interpreted at ¹²CO(7–6) alongside C i(2–1). DRC was then observed through an ALMA DDT program to confirm these lines. However, the observations revealed a single, extremely broad emission line that we initially associated with CO, which then yielded several discrete possibilities for the redshift. We checked three different redshift options using the Jansky VLA, but no CO line was detected. Our ALMA Cycle 4 observations, aimed at confirming the redshift of DRC unambiguously, finally detected a second and a third emission line, neither were at the frequency we were expecting. These detections, of what could only be ¹²CO(6–5) and ${{\rm{H}}}_{2}{\rm{O}}({2}_{11}\mbox{--}{2}_{02})$, revealed that the initial line detected in the ALMA $3\,\mathrm{mm}$ and DDT observations was not CO, but instead [C i](1–0), and that the redshift of DRC is ${z}_{\mathrm{spec}}=4.002$. At this redshift, the ¹²CO(4–3) transition was covered by our $3\,\mathrm{mm}$ spectral scan, but in a noisy region that prevented the line profile from being detected. It should also be noted that because DRC lies close to the edge of the H-ATLAS image of the South Galactic Pole, it sometimes dropped out of the H-ATLAS catalogs, depending on how much padding was adopted, leading to some confusion.

In order to estimate the total SFR of the 10 spectroscopically confirmed DRC components (see Table 1), we use their observed flux densities at $2\,\mathrm{mm}$, corresponding to $\sim 400\,\mu {\rm{m}}$ in the rest frame, which is the only wavelength where the continuum emission from all the DRC components has been detected. The flux density at $2\,\mathrm{mm}$ of each DRC component is converted to an SFR (see values for individual components in Table 1) by rescaling the ALESS template and integrating it between rest-frames 8 and $1000\,\mu {\rm{m}}$. This method assumes that the SEDs (from 8 to $1000\,\mu {\rm{m}}$ in the rest frame) of all DRC components are well represented by the ALESS template, a fact that should be considered when comparing DRC to other protoclusters. Under the same assumption, the continuum sensitivity of our ALMA $2\,\mathrm{mm}$ observations indicates that we are sensitive to sources with $\mathrm{SFR}\gtrsim 40\,{M}_{\odot }\,{\mathrm{yr}}^{-1}\,(3\sigma )$. We include DRC-5 in Table 1, assuming that it belongs to the same protocluster structure, although no emission line has been detected from this component. This is because our observations only covered $v\gt -1000\,\mathrm{km}\,{{\rm{s}}}^{-1}$ and, consequently, we cannot be sure that DRC-5 does not belong to the protocluster.

Among the SED templates used in Ivison et al. (2016), which include a representative range, only the one reported by Pearson et al. (2013) yields lower SFRs than the ALESS SED template, by a factor of $0.66\times$. The Cosmic Eyelash SED (Ivison et al. 2010; Swinbank et al. 2010) and the SED template reported by Pope et al. (2008) yield similar SFRs to the ALESS template. The Arp 220 SED gives higher SFRs than the ALESS template by a factor of $1.36\times$. The same happens for the SED of the lensed source, G15.141 at z = 4.24 (Cox et al. 2011; Lapi et al. 2011), which yields SFRs higher than the ALESS template by a factor of $2.21\times$. Therefore, although the uncertainties can be significant, the SED template used in this work to measure SFRs does not give the highest SFR—we are being conservative (modulo the possibility of a profoundly different IMF in these objects—Romano et al. 2017). Therefore, the high SFRs measured for the DRC components are not artificially high because of the chosen template, but because they are truly luminous.

The total obscured SFR of DRC—which, recall, is only the core of our larger overdensity of DSFGs—is as extreme as $\mathrm{SFR}\sim 6500\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$, and about 75% of that is taking place in three DRC components: DRC-1, DRC-2, and DRC-3, those with ${L}_{\mathrm{IR}}\geqslant 5\times {10}^{12}\,{L}_{\odot }$. The obscured SFR of the protocluster core derived from the ALMA continuum emission at $2\,\mathrm{mm}$ is lower in a factor of $1.37\times$ than that derived from the flux density associated with its extended LABOCA 870 μm emission (components A and B in the left panel of Figure 1), ${S}_{870\mu {\rm{m}}}=64\pm 11\,\mathrm{mJy}$, which would imply $\mathrm{SFR}\,\sim 8900\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$. It is also lower in a $1.34\times$ factor than the obscured SFR associated to the total IR luminosity derived in Ivison et al. (2016) from the Herschel+LABOCA/SCUBA-2 photometry. One reason for this discrepancy could be that the ALESS SED template does not provide a good representation of the dust emission for all the DRC components, such that some of them would have higher SFRs than those reported in Table 1.

The total obscured SFR of our protocluster core is the highest reported so far for protoclusters whose members have been spectroscopically confirmed at $z\gt 4$. As a reference, the total obscured SFR of AzTEC-3, a massive protocluster at ${z}_{\mathrm{spec}}=5.3$ (Capak et al. 2011), is $\mathrm{SFR}\sim 1600\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ as estimated from the only DSFG in this system. For a comparison at slightly lower redshifts, the total obscured SFR in the core of the protocluster SSA22 at ${z}_{\mathrm{spec}}=3.09$ (Steidel et al. 1998; Yamada et al. 2012) is $\mathrm{SFR}\sim 3,820\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$, calculated from the observed flux density at $1.1\,\mathrm{mm}$ of the components and adopting the ALESS SED template, as with DRC). CL J001, a concentration of DSFGs at ${z}_{\mathrm{spec}}=2.5$ (Wang et al. 2016) has $\mathrm{SFR}\sim 3400\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ in the central $80\,\mathrm{kpc}$ region, much lower than DRC. One of the few comparable cases in terms of the high obscured SFR is the COSMOS protocluster at z = 2.10 reported by Yuan et al. (2014) and Hung et al. (2016), but its SFR was determined over much wider scales than the core of our protocluster and, additionally, the COSMOS structure is at much lower redshift than DRC.

As reported above, DRC is just the core of an overdensity of DSFGs. The total $870\,\mu {\rm{m}}$ of all these DSFGs (the eight LABOCA sources, including DRC—see left panel of Figure 1) is $\approx 106\,\mathrm{mJy}$. If all those DSFGs were at the same redshift as DRC, the total SFR of the system would be around ${\rm{14,400}}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$. The fact that the DSFGs found around DRC are at its same redshift is supported by their photometric redshifts measured from the Herschel plus LABOCA/SCUBA-2 photometry (Lewis et al. 2017), although further observations are required to confirm this.

The [C i](1–0) transition has been proposed to be a good tracer of the total molecular gas mass, even better than low-J CO lines in high-redshift galaxies (Papadopoulos et al. 2004) and it has been used to estimate the molecular gas of several populations of high-redshift DSFGs (e.g., Walter et al. 2011; Alaghband-Zadeh et al. 2013; Bothwell et al. 2017). In this work, we use [C i](1–0) to estimate the molecular gas mass for the DRC components detected in that transition (see Figure 3). The line fluxes have been derived from the moment-0 maps of [C i](1–0) for each component. For the DRC components whose [C i](1–0) line is not detected (see, for example, DRC-2, whose [C i](1–0) is not clearly detected in the spectrum—the line is expected to be extremely broad, judging by the ¹²CO(6–5) profile), we have assumed that the [C i](1–0) and ¹²CO(6–5) transitions have the same width. Then, we have calculated the ${{\rm{H}}}_{2}$ mass of each DRC component following Bothwell et al. (2017) and Alaghband-Zadeh et al. (2013), but see also Papadopoulos & Greve (2004), Papadopoulos et al. (2004), and Weiß et al. (2003, 2005):

$$\begin{eqnarray}{M}_{{{\rm{H}}}_{2}} & = & 1375.8\,{D}_{{\rm{L}}}^{2}{(1+z)}^{-1}{\left(\displaystyle \frac{{X}_{[{\rm{C}}{\rm{I}}]}}{{10}^{-5}}\right)}^{-1}\\ & & \times {\left(\displaystyle \frac{{A}_{10}}{{10}^{-7}{{\rm{s}}}^{-1}}\right)}^{-1}\times {Q}_{10}^{-1}\,{S}_{[{\rm{C}}{\rm{I}}]}\,{\rm{\Delta }}v,\end{eqnarray} \tag{ 1 }$$

where ${X}_{[{\rm{C}}{\rm{I}}]}$ is the [C i]/${{\rm{H}}}_{2}$ abundance ratio (assumed to be ${X}_{[{\rm{C}}{\rm{I}}]}=3\times {10}^{-5}$), A₁₀ is the Einstein coefficient (${A}_{10}=7.93\times {10}^{-8}\,{{\rm{s}}}^{-1}$), and Q₁₀ is the excitation factor, assumed to be ${Q}_{10}=0.6$ (Bothwell et al. 2017). Equation (1) does not take into account the effect of the cosmic microwave background (CMB) on the [C i](1–0) line strength (da Cunha et al. 2013; Zhang et al. 2016). The CMB at z = 4 has a temperature of ${T}_{\mathrm{CMB}}\sim 13.7\,{\rm{K}}$, which is not negligible compared to the upper level of the [C i](1–0) line, $\sim 23.6\,{\rm{K}}$. In order to estimate the magnitude of the effect of the CMB, we take into account that the upper level of the [C i] line lies between the upper level energies of ¹²CO(2–1) and ¹²CO(3–2). Following da Cunha et al. (2013), the effect of the CMB would then be a suppression of the velocity-integrated flux of the [C i](1–0) line by a factor of $\sim 2\times$, assuming local thermodynamic equilibrium (the precise factor depends on the kinetic temperature). This means that the molecular gas masses of the DRC components would be a factor of $2\times$ higher than those given by Equation (1). This a very crude estimation, since the effect of the CMB on the observed line fluxes depends on a number of factors. The molecular gas masses of the DRC components, derived with the above assumptions, are quoted in Table 1. It can be seen that our protocluster is formed by galaxies with massive molecular gas reservoirs, whose gas masses range from ${M}_{{{\rm{H}}}_{2}}\sim 1.1\times {10}^{11}\,{M}_{\odot }$ to ${M}_{{{\rm{H}}}_{2}}\sim 2.6\,\times {10}^{11}\,{M}_{\odot }$. The total molecular gas is at least ${M}_{{{\rm{H}}}_{2}}\,\sim 6.6\times {10}^{11}\,{M}_{\odot }$.

Our ALMA $3\,\mathrm{mm}$ scan also covered the ¹²CO(4–3) emission line. It was not detected because of its faintness and because the noise was higher near its frequency. Once the redshift of DRC was confirmed, we calculated the moment-0 map of the ¹²CO(4–3) emission, assuming that it is as wide as the ¹²CO(6–5) transition detected in our ALMA band 4 observations. ¹²CO(4–3) emission is detected in the moment-0 map only from DRC-1. Using this line flux to measure the molecular gas mass of this DRC component, assuming the average luminosity ratio of CO lines in DSFGs, ${L}_{\mathrm{CO}(4-3)}^{{\prime} }/{L}_{\mathrm{CO}(1-0)}^{{\prime} }=0.46$ (Ivison et al. 2011; Bothwell et al. 2013; Carilli & Walter 2013) and the ${\alpha }_{\mathrm{CO}}$ for local ULIRGs, ${\alpha }_{\mathrm{CO}}=0.8{M}_{\odot }{({\rm{K}}\mathrm{km}{{\rm{s}}}^{-1}{\mathrm{pc}}^{2})}^{-1}$ (Downes & Solomon 1998), the derived molecular gas mass for DRC-1 from ¹²CO(4–3) is ${M}_{{{\rm{H}}}_{2}}\sim 1.1\times {10}^{11}\,{M}_{\odot }$. This value is lower than the molecular gas mass of DRC-1 derived from [C i](1–0), but is consistent given the large uncertainties involved in such determinations.

The molecular gas mass for the other DRC components could be determined from the detected ¹²CO(6–5) line transitions. However, this would lead to a very uncertain estimates because of the conversion from ¹²CO(6–5) to ¹²CO(1–0) luminosity and the conversion from ¹²CO(1–0) luminosity to molecular gas mass. Alternatively, we could use the rest-frame $870\,\mu {\rm{m}}$ luminosity as a proxy for the molecular gas mass, following Scoville et al. (2016), Hughes et al. (2017), and Oteo et al. (2017c). For the DRC components with [C i](1–0) detections, the molecular gas masses derived from [C i](1–0) are $\sim 2\mbox{--}5\times$ lower than those derived from the dust continuum luminosity, where the rest-frame luminosities at $870\,\mu {\rm{m}}$ have been obtained by redshifting the ALESS template to z = 4.002 and rescaling to the observed flux density of each DRC component at $2\,\mathrm{mm}$. We report only the molecular gas mass of the sources detected in [C i](1–0), noting that the total molecular gas mass of our protocluster might be much higher (also considering that some of the LABOCA-detected DSFGs around DRC might be at the same redshift).

We can estimate the gas-depletion time of the DRC components detected in [C i](1–0) from the ratio between their molecular gas masses and their SFRs. These are shown in Table 1 and range between ∼90 and 230 Myr. These gas-depletion times are comparable to those obtained for other luminous DSFGs at $z\gtrsim 4$ (Riechers et al. 2013, 2017; Hodge et al. 2015; Oteo et al. 2016b; Fudamoto et al. 2017).

The probability of detecting 10 short-lived, physically associated DSFGs—we do not consider DRC-5 here due to the lack of spectroscopic confirmation—is extremely low (see, for example, Casey 2016), meaning that such structures are extraordinary systems. The fact that some DRC components are detectable as luminous DSFGs at the same redshift and with relatively low gas-depletion times might suggest the existence of a mechanism that is able to trigger star formation simultaneously in different sources distributed across hundreds of kiloparsecs. Alternatively, given that the baryon to dark-matter ratio in galaxies is much lower than the cosmic value, such that halos contain huge amounts of gas that can flow to the star-forming regions, one can argue that star formation can be sustained over longer times until it is swept out by feedback effects; though, the absence of "normal" SFGs is puzzling. Relatively long starburst lifetimes have been reported in the past (Granato et al. 2004; Lapi et al. 2011), while others have reported the presence of galaxies in high-redshift clusters whose scatter in age is only a few hundreds megayears (Andreon et al. 2014).

Using the ¹²CO(6–5) line detections, which provide the highest signal-to-noise among all the detected lines, we have measured the central velocity of all the DRC components (see Table 1), from which we can derive a velocity dispersion, ${\sigma }_{{\rm{v}}}=794\,\mathrm{km}\,{{\rm{s}}}^{-1}$. This relatively high velocity dispersion is mainly caused by the velocity of components DRC-7 and DRC-10 with respect to the others. The fact that these two components have velocity offsets as large as $\sim 2000\,\mathrm{km}\,{{\rm{s}}}^{-1}$ with respect to the average velocity of all the other components might mean that we are seeing two groups, as has been reported in protoclusters at lower redshifts such as MRC 0052-241 at z = 2.86 (Venemans et al. 2007). We note that the velocity dispersion has been obtained from only 10 spectroscopically confirmed sources in the core of the protocluster.

We can estimate the total mass of our protocluster using several different methods. First, we can use the velocity dispersion (${\sigma }_{{\rm{v}}}$) and the relation between the velocity dispersion and the total mass (${M}_{\mathrm{total}}$) derived by Evrard et al. (2008), which has also been used to estimate the total mass of lower-redshift protoclusters (Wang et al. 2016):

$$\begin{eqnarray}&&{\sigma }_{{\rm{v}}}(M,z)={\sigma }_{\mathrm{DM},15}{\left(\displaystyle \frac{h(z){M}_{\mathrm{total}}}{{10}^{15}{M}_{\odot }}\right)}^{\alpha },\end{eqnarray} \tag{ 2 }$$

where h(z) is the dimensionless Hubble parameter. Using ${\sigma }_{\mathrm{DM},15}\sim 1,083\,\mathrm{km}\,{{\rm{s}}}^{-1}$, and $\alpha \sim 0.336$ (see also Wang et al. 2016), we derive a total mass of ${M}_{\mathrm{total}}\sim 9.3\times {10}^{13}\,{M}_{\odot }$. We note that this estimate is highly uncertain for two main reasons. First, protoclusters at $z\sim 4$ are not virialized. Second, the velocity dispersion has been measured using only the spectroscopically confirmed sources, which are all located in the core (the most violent region, which is forming stars at a rate of thousands of ${M}_{\odot }\,{\mathrm{yr}}^{-1}$). We can also estimate the total mass of the cluster by assuming a uniform spherical distribution with line-of-sight velocity dispersion and radius corresponding to those for our protocluster core. In this case, the virial mass would be ${M}_{\mathrm{total}}\sim 3.2\times {10}^{13}\,{M}_{\odot }$, slightly lower than the previous estimate. Next, we estimate the total mass using the relation between the total IR luminosity and halo mass at z = 4 by Aversa et al. (2015), obtaining ${M}_{\mathrm{total}}\sim 4.4\times {10}^{13}\,{M}_{\odot }$.

The different methods used above to estimate the total mass of the protocluster core give a wide range of values, indicative of the uncertain nature of this task. Despite the uncertainties, we have determined that our protocluster core is very massive. The total mass of the full protocluster would be considerably higher still, considering that DRC is only the core. It could be even more massive than the most massive progenitor halos predicted by simulations. As a reference, Chiang et al. (2013) predicted $\lt 2\times {10}^{13}\,{M}_{\odot }$ at z = 4. Using the evolutionary tracks of Chiang et al. (2013), we conclude that DRC might become an ultra-massive cluster at z = 0, with a total mass in the region of $2\times {10}^{15}\,{M}_{\odot };$ thus, DRC could be the progenitor of a Coma-like (The & White 1986; Kubo et al. 2007) or even more massive cluster. However, we note that exploring the evolution of DRC is highly uncertain because we are only probing the core of a likely larger overdensity of DSFGs and halos with similar masses can evolve very differently. Further observations of the surroundings of DRC would be needed to have a better insight of the evolution (i.e., if DRC is truly the progenitor of a massive galaxy cluster or it is just a group of sources caught in a active stage during their evolution) of this extreme structure we have discovered.

Figure 4 shows a false-color image of the stellar emission in the components of our protocluster core, which has been created using R, K_s, and IRAC $4.5\,\mu {\rm{m}}$ imaging. A variety of colors can be seen, with the brightest components at $2\,\mathrm{mm}$ being associated to reddest sources, as expected. We note that our optical and near-IR imaging does not reveal any sign (like arcs), which might indicate that our protocluster is being gravitationally amplified by a foreground cluster.

Using our MUSE observations of the protocluster core, we have created a continuum-subtracted Lyα image, which is shown in the top panel of Figure 5. The Lyα image reveals the presence of a Lyα blob next to DRC-2 that extends over $60\,\mathrm{kpc}$; similar to what is found around other high-redshift DSFGs (Ivison et al. 1998; Chapman et al. 2001; Umehata et al. 2015; Geach et al. 2016; Cai et al. 2017) and radio galaxies (Swinbank et al. 2015). Actually, as discussed in Chiang et al. (2015), there is a tendency for there to be extended Lyα emission in overdense environments (Matsuda et al. 2009; Yang et al. 2009; Erb et al. 2011; Matsuda et al. 2012). We note that part of the emission in this Lyα blob might be contaminated by a spectroscopically confirmed [O ii] emitter located at the south the Lyα blob (there is an additional [O ii] emission northwest of the MUSE field. The lower panel of Figure 5 shows the spectrum of the Lyα blob, integrated over its full extent, in comparison to the ¹²CO(6–5) line profile of the DRC component right next to it. The Lyα emission is even wider than the ¹²CO(6–5), but is still comparable to other Lyα blobs at these redshifts.

Zoom In Zoom Out Reset image size

It can be seen in Figure 4 that most of the DRC components detected at CO and $2\,\mathrm{mm}$ continuum do not show significant Lyα emission, as expected for DSFGs (although there are exceptions to this—Chapman et al. 2003, 2005; Casey et al. 2012; Oteo et al. 2012a, 2012b; Sandberg et al. 2015). Furthermore, apart from the Lyα blob, no Lyα emitters are found in the $1^{\prime} \times 1^{\prime}$ central region of the protocluster (the emission-line galaxy detected in the northeast is a low-redshift [O ii] emitter, with strong continuum and Balmer absorption lines), indicating that the star formation in this most extreme region is dominated by dust-obscured star formation. This might be a consequence of an extra Lyα depletion in the core of high-redshift protoclusters, as discussed in Shimakawa et al. (2017b).