© ESO 2020

1. Introduction

The discovery of bright submm-selected galaxies (SMGs; S(850 μm) > 5 mJy; e.g., Blain et al. 2002) opened up a window onto the formation of the most intense episodes of star formation in massive distant galaxies. It is now widely accepted that understanding the nature of SMGs is a prerequisite to understanding the formation of local massive ellipticals and brightest cluster galaxies (BCGs) at the center of the most massive virialized structures. Defining starbursts as galaxies above the star formation main sequence (MS, star formation rate vs. stellar mass, SFR-M_⋆ correlation; Noeske et al. 2007; Elbaz et al. 2007, 2011; Daddi et al. 2007; Rodighiero et al. 2011), it has been found that some SMGs fall into the starburst category while some are MS galaxies (massive MS galaxies at z = 3–5 exhibit very high star formation rates (SFRs), Schreiber et al. 2015), some exhibit compact star formation, and others exhibit extended star formation (see e.g., Hodge et al. 2019; Rujopakarn et al. 2016, 2019; Elbaz et al. 2018).

Using ALMA to follow up on 63 optically dark galaxies located in the three Southern Hemisphere CANDELS fields (totalizing 600 arcmin² in the COSMOS, EGS, and GOODS-South fields), Wang et al. (2019) showed that these galaxies contribute ten times more than equivalently massive UV-bright galaxies to the cosmic SFR density at z > 3. Despite being bright – although three to ten times fainter than the classical SMG population – in the submillimeter (here 870 μm), these galaxies unseen in the deepest UV to near-infrared surveys were selected as H-dropouts with an Infrared Array Camera (IRAC) counterpart brighter than 24 AB at 4.5 μm. Thanks to the angular resolution of ALMA (06 here), it is indisputable that these galaxies exhibit no counterpart from the UV to the H-band. While in the observations from single dish telescopes such as Submillimeter Common-User Bolometer Array (SCUBA), LABOCA LArge APEX BOlometer (LABOCA), the IRAM 30 m, and Large Millimeter Telescope (LMT), several optical sources can play the role of candidate counterparts due to the large beam size. It took 14 years (Walter et al. 2012) to identify the galaxy responsible for the brightest source in the HUDF (Hughes et al. 1998, HDF850.1). The H-dropout (or optically dark) galaxies followed up by ALMA experience lower, hence more normal, star formation activity than the more classical SMGs, such as those followed up by ALMA as part of the ALESS project (see e.g., da Cunha et al. 2015, Swinbank et al. 2014). The ALMA detected ones exhibit ⟨SFR⟩ = 200 M_⊙ yr⁻¹ typical of MS galaxies at these redshifts although with a large scatter, whereas the undetected ones experience much lower SFRs (Wang et al. 2019).

The clustering properties of these H-dropout galaxies – derived from their cross-correlation with H-detected galaxies with similar redshifts and masses – show that they inhabit massive dark matter halos (M_h ∼ 10^{13 ± 0.3}h⁻¹ M_⊙ at z = 4) and constitute excellent candidate progenitors of the most massive galaxies presently located in the center of massive groups and galaxy clusters (Wang et al. 2019). Recent simulations strengthen this hypothesis by showing that galaxies in (proto)clusters at z > 3 may produce as much as 50% of the cosmic SFR density (Chiang et al. 2017). Numerical simulations can be used to trace back the positions of the galaxies that will end up in a present-day galaxy cluster. Chiang et al. (2013) followed up these galaxies in the stage preceding virialization, which they called the proto-cluster phase. They also defined a useful reference effective radius encompassing 65% of the mass that will end up in the cluster. At a redshift of z ∼ 3, which will be of interest in our analysis, this radius ranges from 5 comoving-Mpc (hereafter cMpc) for a M(z = 0) = 3 × 10¹⁴ M_⊙ bound halo to > 8 cMpc for a Coma-like M(z = 0) = 10¹⁵ M_⊙ halo. This is in line with the enclosed sizes, which contain 90% of the stellar mass in proto-clusters, found in Muldrew et al. (2015). Identifying a proto-cluster from observations is much more complex since it is difficult to determine the future of any given overdensity of galaxies. However, when studying galaxy overdensities the effective radius defined from numerical simulations may be used as a reference size for a potential proto-cluster. Based on the standard hierarchical models of galaxy formation, numerical simulations demonstrate that massive galaxies at high redshifts are useful tracers of overdense regions since they preferentially form within the peaks in the density field (Springel et al. 2005; De Lucia & Blaizot 2007). At the same time, recent ALMA surveys have revealed that ALMA is a powerful tool for detecting massive galaxies at high redshift (Dunlop et al. 2017; Hatsukade et al. 2018; Franco et al. 2018) taking advantage of its high sensitivity and high spatial-resolution. As a result, not only can such a population represent a unique probe of star formation at cosmic epochs when dust attenuation is largely unknown (i.e., at z > 3), but it can also be seen as a potential beacon to identify clusters at the epoch of their formation. Studying the connection of optically dark ALMA galaxies with the formation of large-scale structures is the main goal of the present paper.

The study on how galaxy clusters form and evolve across cosmic time is critical for testing models of large-scale structure formation. The formation of galaxy clusters is related to the collapse of the gravitationally bound overdensities originating from the peaks in the initial density fluctuations since the Big Bang (Kravtsov & Borgani 2012). It is also well known that galaxies evolve differently in different environments at early times: Galaxies in denser regions tend to form earlier and evolve faster, as indicated by their ages and enhanced SFR, than galaxies in less dense regions (Thomas et al. 2005; Elbaz et al. 2007; Gobat et al. 2008). Therefore, by exploring the assembly of the most massive dark matter halos we hope to better understand how massive galaxies formed, in the framework of large-scale structure formation.

Various techniques have been developed to search for high-redshift galaxy clusters. However, a complete search for galaxy clusters at high redshift (≳2) becomes more difficult as the cluster properties evolve dramatically at an early age. Some galaxy clusters are mature systems with evidence of quiescent galaxies dominating the center (Papovich et al. 2012; Newman et al. 2014; Willis et al. 2020), while others still show substantial star formation activities in the center (Brodwin et al. 2013; Webb et al. 2015; Wang et al. 2016). As a result, infrared and submilimeter surveys can be beneficial in revealing the overdensities (Daddi et al. 2009; Capak et al. 2011; Miller et al. 2018; Oteo et al. 2018; Gómez-Guijarro et al. 2019; Castignani et al. 2020), and extensive in-depth deep surveys from Atacama Large Millimeter Array (ALMA) can make significant contributions.

In this paper, we present the physical properties of six optically dark galaxies detected in the GOODS-ALMA field based on the new ALMA data obtained from the spectroscopic follow-up of five of these galaxies. We also discuss the connection between their optically dark nature and the dense environment they reside in.

The GOODS-ALMA survey makes use of the 69 arcmin² ALMA image covering the deepest region of the Great Observatories Origins Deep Survey South (GOODS-South) field. In this large blind survey, massive (> 10^10.5 M_⊙) and dust-rich galaxies at high redshift are efficiently detected due to their low surface number density in the early Universe. For more details of this survey we refer the readers to Franco et al. (2018, hereafter F18).

We use cosmological parameters of H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_M = 0.3, and Ω_Λ = 0.7. A Salpeter (1955) initial mass function (IMF) is adopted to derive stellar masses and SFRs. We convert values obtained by other studies from a Chabrier (2003) IMF to a Salpeter (1955) IMF by multiplying SFR and stellar masses by the same factor, 1.74.

2. Data and observations

2.1. ALMA data and observations

We present in the following the GOODS-ALMA survey data and the ALMA spectroscopic scan data we use for the optically dark galaxies.

2.1.1. GOODS-ALMA survey data

We use data from GOODS-ALMA, a 1.1 mm cosmological survey of 69 arcmin² (69 × 10′) in the GOODS-South field (program 2015.1.00543.S, PI: D.Elbaz). GOODS-ALMA reaches an rms sensitivity of σ ≃ 0.18 mJy beam⁻¹ in the 060 tapered mosaic (F18). The comparison between the redshift distribution of the detections in the GOODS-ALMA survey (F18, F20a) and the one in the smaller pencil beam ALMA survey of the Hubble Ultra Deep Field (4.5 arcmin², Dunlop et al. 2017) shows that, for the same total observing time of about 20 hours, the shallower but 15 times wider GOODS-ALMA survey enables the detection of systematically more distant galaxies (F18, F20a). This is a natural consequence of the fact that high-redshift dusty galaxies are found to be massive. The combination of the low surface density and brightness of these galaxies then favors shallower and wider surveys.

A total of 35 galaxies have been detected above 3.5σ in GOODS-ALMA (Franco et al. 2020a, hereafter F20a). These include 19 galaxies detected above the 4.8σ limit from a blind detection approach (F18) and 16 galaxies within 3.5 and 4.8σ detected using ancillary information, mainly Spitzer-IRAC prior positions (F20a). The median redshift and stellar mass of the S/N  ≥  4.8 sources are z = 2.73 and M_⋆ = 1.0 × 10¹¹ M_⊙, whereas the sources with a 4.8 > S/N  ≥  3.5, are both slightly closer (z = 2.40) and less massive (M_⋆ = 7.2 × 10¹⁰ M_⊙).

2.1.2. ALMA follow-up observations of five optically dark galaxies

Out of the 35 ALMA galaxies presented in Sect. 2.1.1, an ensemble of six galaxies are optically dark, or also called HST-dark since they are not detected in the optical to near-infrared HST images down to a 5σH-band detection limit (HST-WFC3/F160W) of H = 28.16 AB mag (Guo et al. 2013). As a result, about 17% of the sources detected above 3.5σ in the GOODS-ALMA blind survey are optically dark (and 21% above 4.8σ). Out of these six sources initially classified as optically dark, we will show in Sects. 4.1 and 5.3.2 that two are associated with an H-band counterpart. They were missed from the H-band catalog due to blending.

We have obtained ALMA spectroscopic follow-up data for five of these optically dark galaxies (AGS4, 11, 15, 17, 24) that we discuss in Sect. 4. The sixth one, AGS25, was identified after the ALMA call hence no follow-up data have been obtained for this one.

We followed up these five optically dark sources with ALMA band 4 using the spectroscopic scan mode in Cycle 7 (project 2018.1.01079.S PI: M. Franco). The band 4 observations were performed by combining four tunings to cover a total frequency range of 137.2–165.8 GHz that ensures a broad redshift coverage from z = 2 to z = 5, with at least one line detected among the possible CO transitions from J = 4 to J = 7, CI and H₂O. The observations were carried out with the array configuration C43-4, with a synthesized beam of ∼088. Each source was observed for 13 min on-source time in each tuning, reaching a typical rms sensitivity of 0.2 mJy per 400 km s⁻¹. We reduced the raw data using the standard ALMA pipeline with Common Astronomy Software Application package (CASA, McMullin et al. 2007). Then we converted the data into uvfits format and performed analysis with the IRAM GILDAS tool (Gildas Team 2013) using the same method as in Jin et al. (2019).

2.2. Ancillary data

GOODS-ALMA benefits from abundant ancillary data. We summarize in the following Sections the data that we use in the present study. For more details we refer the reader to Franco et al. (2018, Sect. 2.4.).

2.2.1. Optical/near-infrared data

We use the data from the Cosmic Assembly Near-IR Deep Extragalactic Legacy Survey (CANDELS; Grogin et al. 2011) obtained using images from the Wide Field Camera 3 (WFC3) and the Advanced Camera for Surveys (ACS) on board the Hubble Space Telescope (HST) in the filters: Y₁₂₅, J₁₂₅, B₄₃₅, V₆₀₆, i₇₇₅, i₈₁₄ and z₈₅₀. We also use images from the VLT, obtained in the U-band with VIMOS (Nonino et al. 2009), and in the K_s-band with ISAAC (Retzlaff et al. 2010) and HAWK-I (Fontana et al. 2014) and the associated CANDELS multiwavelength catalog (Guo et al. 2013).

In the near-infrared, we also use images and data from the FourStar galaxy evolution survey (ZFOURGE, PI: I. Labbé; Straatman et al. 2016). These include narrow-band filter images in the J (J₁, J₂, J₃) and H (H_s, H_l) bands useful for the photometric redshift determination. The ZFOURGE team produced the deepest K_s-band detection image by combining images from the survey itself together with images from previous surveys and reaching a 5σ limiting depth of 26.2–26.5 AB mag in the whole Chandra Deep Field South (Giacconi et al. 2002, CDFS).

2.2.2. Mid-infrared data

We use Spitzer Infrared Array Camera (IRAC) images at 3.6 and 4.5 μm, reaching a depth of 26.5 AB mag, from the ultradeep mosaics of the IRAC Ultra Deep Field (Labbé et al. 2015, IUDF). These were obtained by combining the observations from the IUDF (PI: Labbé) and IRAC Legacy over GOODS (IGOODS, PI: Oesch) programs, together with archival data from GOODS (PI: Dickinson), S-CANDELS (PI: Fazio), ERS (PI: Fazio), and UDF2 (PI: Bouwens).

2.2.3. Radio data

We use radio images at 3 GHz and 6 GHz from the Karl G. Jansky Very Large Array (VLA; PI: W. Rujopakarn), which cover the entire GOODS-ALMA field. The 6 GHz map has an angular resolution of 031 × 061. The exposure time is in total 177 hours reaching an rms noise at the pointing center of 0.32 μJy beam⁻¹. The 3 GHz map has an angular resolution of 119 × 059 and an average rms noise of 1.03 μJy beam⁻¹. At 3 GHz, AGS4, 17 and 24 are detected with an S/N of 10.4, 9.9, and 5.6 using PyBDSF¹ (Rujopakarn et al., in prep.). We examined the image and noticed that the peak fluxes of AGS11, 15 and 25 are at 3 < S/N < 5. At 6 GHz, only AGS4 is detected with S/N > 5 using PyBDSF¹, AGS17, 24 and 25 show peak fluxes at 3 < S/N < 5 in the image. We list the fluxes in Table 1.

2.2.4. X-ray data

We use the catalog of Luo et al. (2017) obtained from the 7 Msec X-ray data in the CDFS field observed in three bands: 0.5–7.0 keV, 0.5–2.0 keV, and 2–7 keV. This catalog includes a classification of candidate X-ray active galactic nuclei (AGNs) that we discuss in the following of the paper.

2.3. Origin of the redshifts and stellar masses

In total, 11 621 galaxies located in the GOODS-ALMA field have been attributed either a photometric redshift (92% of the sample) or a spectroscopic redshift (894 galaxies). We discuss below the origin of these redshifts, including the newly obtained ones from the VANDELS survey that we use to improve the contrast on the z ∼ 3.5 structure that we will discuss in the following sections.

2.3.1. Photometric redshifts and stellar masses

The photometric redshifts used in the present study come from the catalog generated by the ZFOURGE team (Straatman et al. 2016; Forrest et al. 2017). We use the redshift maximizing the likelihood (minimized χ²), z_peak, as derived by fitting the observed SEDs with the code EAzY (Brammer et al. 2008). These photometric redshifts present an excellent agreement with existing spectroscopic redshifts with σ_z = 0.020 at z > 1.5 (Straatman et al. 2016), where σ_z is the normalized median absolute deviation (NMAD) of Δz/(1 + z), that is, 1.48× the median absolute deviation (MAD) of |z_phot  −  z_spec|/(1 + z_spec).

We adopt the stellar masses of the ZFOURGE catalog. The stellar masses were derived with Bruzual & Charlot (2003) stellar population synthesis models (Straatman et al. 2016) assuming exponentially declining star formation histories (SFH) and a Calzetti et al. (2000) dust attenuation law.

2.3.2. Spectroscopic redshifts

The spectroscopic redshifts are both from the literature, as compiled in the ZFOURGE catalog, and from the latest data release from the VANDELS survey (DR3²). VANDELS targets star-forming galaxies at 5.5 ≥ z ≥ 2.5 and massive passive galaxies at 2.5 ≥ z ≥ 1.0 in the CANDELS CDFS and UDS fields. It is an ultra-deep spectroscopic survey of high-redshift galaxies with exposure times ranging from 20 to 80 hours per source (Pentericci et al. 2018; McLure et al. 2018). In total, we obtain 132 new spectroscopic redshifts, z_spec, in the GOODS-ALMA field from VANDELS and 83 of them fall in the redshift range 4 ≥ z ≥ 3.

2.4. Derived parameters of the optically dark galaxies

2.4.1. Redshifts

To derive the photometric redshifts, we first de-blend³ AGS4, 15, and 24 from their neighbors in the observed optical to mid-infrared images (ancillary data in Sect. 2.2), as to be discussed in detail in Sect. 4. The photometric data are then used to measure the redshifts by fitting the optical to mid-infrared (rest-frame UV to near-infrared) SEDs with EAzY (Brammer et al. 2008). The photometric redshift of AGS25 comes from its K_s-band counterpart in the ZFOURGE catalog, ID_ZFOURGE = 11 353. AGS11 and 17 are either too faint or too close to a bright emitter in the images to be de-blended properly and derive meaningful redshifts, therefore we do not have photometric redshifts for these two galaxies.

The spectroscopic redshifts of AGS4 and AGS17 are obtained based on the detected emission lines in the ALMA spectroscopic follow-up and their SEDs at UV to near-infrared and far-infrared (rest-frame). This will be discussed in detail in Sect. 3.

2.4.2. Stellar masses

The stellar masses of AGS4, 15 and 24 are derived from the observed optical to mid-infrared SED fitted with FAST++⁴ assuming a delayed exponentially declining SFH and a Calzetti et al. (2000) attenuation law. The stellar masses of AGS11 and 17 are calculated based only on the IRAC emission. The stellar mass of AGS25 comes from its K_s-band counterpart in the ZFOURGE catalog.

2.4.3. Infrared luminosities

Infrared luminosities are measured by the SED fitting of the Herschel and ALMA flux densities for AGS4 and AGS17 using CIGALE (Boquien et al. 2019) as in F20b. The other four galaxies are not detected by Herschel, then their infrared luminosity was derived by fitting the ALMA 1.1 mm point with the SED templates in Schreiber et al. (2018a), as described in F20b.

2.4.4. SFRs

The SFRs are measured in Franco et al. (2020b, hereafter, F20b), as SFR_tot = SFR_IR + SFR_UV. We calculate the SFR_IR based on the correlation between infrared luminosity and SFR in Kennicutt (1998) and SFR_UV contributes only 1% to SFR in the GOODS-ALMA detections (F20b), therefore it is negligible.

All six galaxies are detected in the radio at 3 GHz with flux densities ranging from 4 μJy to 40 μJy (Table 1). This allows us to estimate a radio SFR. With such extremely faint flux densities, they are not expected to present a radio excess that could be attributed to an AGN, but we checked this as follows. We calculated the luminosity at 1.4 GHz (restframe) from the flux at 3 GHz assuming a radio spectral index of α = = −0.8, the typical value for galaxies at z > 2 (Delhaize et al. 2017):

(1)

where D_L is the luminosity distance to the object. The results are listed in Table 1. Then we compared the ratio L_1.4 GHz/SFR with the radio excess threshold as in Delvecchio et al. (2017). All the six optically dark galaxies have L_1.4 GHz/SFR well below threshold, with the one of AGS24 being five time lower and the others at least ten times lower. This is fully consistent with star-forming galaxies at z ∼ 3.5. A radio AGN contribution cannot be ruled out, but it must be contributing for a very small fraction to the total radio emission. Nevertheless, due to vigorous star formation dominating the total radio emission at low resolution, only Very Long Baseline Interferometry (VLBI) observations could unambiguously shed light on the possible presence of radio AGN activity on circum-nuclear (< 100 pc) scales.

We also checked in the CDFS 7Ms catalog (Luo et al. 2017) and found that an X-ray source is located equally close to AGS4 and another galaxy, ID_ZFOURGE = 12 333. Therefore, it is not clear whether this X-ray source is associated with our optically dark source, AGS4. A detailed study of the AGN contamination to the SFR of AGS4 is presented in F20b where it is found to be negligible.

3. Results of the ALMA spectroscopic follow-up

In this section, we present the results of the ALMA spectroscopic follow-up. As introduced in Sect. 2.1.2, five of the six optically dark galaxies are followed up, as AGS25 was identified after the ALMA call. One emission line is detected in AGS4 and AGS17, and we constrained the upper limits for AGS11, 15 and 24.

3.1. AGS4

The spectrum of AGS4 shows a clear line detection with an S/N = 8.58 at ν_obs = 151.44 GHz (see Fig. 1-left). One ALMA spectral line alone is not sufficient to obtain a definitive spectroscopic redshift for AGS4. However, out of the various possible lines that may be responsible for the one that is detected, only three lines can be reconciled with the 4000 Å break that is seen in the optical to near-infrared SED of AGS4 (Fig. 2-top-left, see more discussion on the SED in Sect. 4.1), namely the CO(6-5), CO(7-6) and H₂O(2₁₁-2₀₂) lines.

Fig. 1. Left: part of the spectrum obtained from the ALMA band 4 observation of AGS4. The yellow shade highlights the detected emission line. The emission line map is shown in the bottom-left corner. The contours denote 1 to 5σ, with the step of 1σ. Right: same as on the left, but for AGS17. The contours in the map denote 1 to 6σ, with the step of 1σ.

Fig. 2. Top-left: SED of AGS4. Top: photometric redshift PDF derived by EAzY (Brammer et al. 2008). Bottom: SED fittings. The green curve shows the best-fit of AGS4 at zphot = 4.02 using EAzY, and we shift the SED to the spectroscopic redshift (zspec = 3.556) of AGS4 in light blue for comparison. The gray curve shows the SED of the optically bright neighbor, IDZFOURGE = 12 333. Top-right: SED of AGS15. The blue curve shows the best-fit at zphot = 3.96 derived by EAzY and we shift the SED to the redshift peak of the overdensity, at z = 3.472, in light blue for comparison. The gray curve shows the SED of the optically bright neighbor, IDCANDELS = 3818. Bottom-left: SED of AGS24. Same as AGS15. Bottom-right: SED of AGS25 at zphot = 4.70, fitted by EAzY with the photometry retrieved from the ZFOURGE catalog.

Another reason to favor these lines comes from the far-infrared SED (Fig. 22 in F18). Being one of the brightest GOODS-ALMA sources, AGS4 is detected in five Herschel bands (at 100, 160, 250, 350, and 500 μm), which together with the ALMA measurements, allowed F18 to obtain a robust identification of the peak of the far-infrared SED, which falls close to 350 μm. Redshifts below z = 2 (such as the one obtained for the 4-3 CO transition) would lead to dust temperatures (peak around 120 μm) much colder than typical z ∼ 2 galaxies (see Schreiber et al. 2018a).

We exclude H₂O(2₁₁-2₀₂). It is typically less bright or as bright (Yang et al. 2016, flux ratio from 0.4 to 1.1) as the neighboring high-J CO transition lines (7-6 or 6-5). However, we do not find any evidence for a second line in the ALMA spectrum that can be comparable to the 8.7σ detection.

The CO(7-6) line with a rest-frame frequency of ν_RF = 806.7 GHz provides a spectroscopic redshift of = 4.326, which agrees with one of the two peaks of the probability distribution function (PDF) of the optical photometric redshift. We note that the [CI](2-1) (ν_rf = 809.3 GHz) line falls very close to the CO(7-6) line and is not detected here with a 3.5σ limit of 0.7 mJy, which corresponds to a [CI](2-1)/CO(7-6) ratio of 0.40. This flux ratio falls at the lower end of the trend defined by the MS galaxies, starburst galaxies and AGNs in the local universe and at high-z (data compiled and described in Valentino et al. 2018, 2020), as shown in Fig. 3. Chance is small to have AGS4 to be an extreme case of the [CI](2-1)/CO(7-6) flux ratio. The CO(6-5) line with a rest-frame frequency of ν_RF = 806.7 GHz gives a = 3.566, which agrees with the 4000 Å discontinuity in the optical SED. Therefore, we adopt the redshift to be = 3.556 assuming the line to be CO(6-5).

Fig. 3. Upper limit of the [CI](2-1)/CO(7-6) flux ratio of AGS4 (blue arrows), compared with existing observations of local and high-redshift MS galaxies, starburst galaxies and AGNs (data compiled in Valentino et al. 2018, 2020). If the detected emission line is indeed CO(7-6), then the flux ratio would fall at the lower end of the trend shown by the existing observations.

3.2. AGS17

Our spectroscopic scan follow-up of AGS17 shows an emission line at 154.78 GHz with an S/N = 10.23. The far-infrared SED in this region peaks at around 400 μm (F18, Fig. 22). Taking this far-infrared SED into account, the detected line is likely to be CO(5-4) or CO(6-5) or CO(7-6) since the higher-J or lower-J transitions give redshifts leading to unrealistic dust temperatures with dust emission peaking at λ < 67 μm or at λ > 134 μm, restframe. We checked the three L_IR − L_CO relation for the three transitions. We found that the ALMA flux of AGS17 and the corresponding L_IR for the three different redshifts, which are z = 2.723, 3.467, 4.212, will bring the galaxy to the positions that are consistent with the L_IR – L_CO relations of the three transitions (Liu et al. 2015). Therefore, we cannot disentangle these possibilities based on the single line.

However, if we assume the line to correspond to the CO(7-6) transition (ν_rf = 806.7 GHz), a second line, [CI] (ν_rf = 809.3 GHz), is supposed to be detected at similar significance (Valentino et al. 2020) and it is not the case for AGS17. We are then left with two possibilities: a spectroscopic redshift of = 2.723 in the case of the CO(5-4) transition and the other of = 3.467 for the CO(6-5) transition. The second redshift of = 3.467 turns out to fall within a redshift peak in the redshift distribution of the galaxies in the field, which also exhibits a concentration of galaxies right around AGS17, as we will discuss in more details in Sect. 5.3. For that reason, despite the fact that we cannot strictly rule out the possibility that we are here observing the CO(5-4) transition, we decided to favor the CO(6-5) transition from a purely probabilistic point of view. A mild confirmation of this choice comes from the dust temperature that we obtain when assuming one or the other transition and redshift. If the transition is CO(5-4), the spectroscopic redshift = 2.723 implies a dust temperature T_dust ∼ 30 K, using a simple conversion of the peak of the SED at around 400 μm to the dust temperature. This is lower than the typical temperature T_dust ∼ 38 K for the NS galaxies at this redshift (Schreiber et al. 2018a). Instead, if we assume the transition to be CO(6-5), the = 3.467 indicates a dust temperature T_dust ∼ 37 K that is close to the typical dust temperature T_dust ∼ 40 K at z = 3.5 (Schreiber et al. 2018a) and is consistent with the mean dust temperature, T_dust = 37 K, of the 39 optically dark galaxies in Wang et al. (2019). We note that we also tried to obtain an independent confirmation of the redshift of AGS17 from its UV to near-IR SED as we did for AGS4. However, as we will discuss in Sect. 5.3.2 and as shown in Fig. 4-bottom, the optical emission of AGS17 suffers from a strong blending with bright neighbors that prevents us from extracting a meaningful SED that we could use to obtain a robust, or even tentative, model fit and photometric redshift derivation.

Fig. 4. Top: images of AGS11 in the J band, H band, Ks band and IRAC 3.6 μm. The magenta “x” denotes the detection in Ks-band by ZFOURGE, IDZFOURGE = 7589 (zphot = 4.82). Bottom: same as above, but for AGS17. The cyan “x”s denote the optical bright neighbors, IDCANDELS = 4414 (zphot = 0.03, below) and IDCANDELS = 4436 (zphot = 0.95, above). The magenta “x” denotes the detection in Ks-band by ZFOURGE, IDZFOURGE = 6964 (zphot = 1.85).

3.3. Upper limits of AGS11, AGS15, and AGS24

Our spectroscopic scan follow-up of AGS11, AGS15, and AGS24 does not show emission line detections higher than 5σ. Based on their assumed redshifts, as will be discussed in Sect. 5.3, the ALMA spectroscopic scans encompass the position of the CO (6-5) emission line in the three galaxies. At the assumed redshifts, the infrared luminosities can be converted to CO(6-5) luminosities using the relation defined by Liu et al. (2015) and then to CO(6-5) fluxes accordingly. We find that the non-detection of the CO(6-5) line agrees with the dispersion of the relation with values below 3σ for the three galaxies (2.8σ, 2.3σ, 0.7σ for AGS11, AGS15, and AGS24, respectively).

4. GOODS-ALMA optically dark galaxies

As described in Sect. 2.1.1, six out of the 35 GOODS-ALMA detections do not have counterparts in the CANDELS catalog (Guo et al. 2013). Despite their non detection in the H-band images, we are able to measure the residual emission in the optical to near-infrared images after subtracting the contribution of neighboring galaxies and to use the resulting measurements to obtain photometric redshift determinations as discussed in the next sections. After a careful analysis of the H-band images of these galaxies, we realized that in two cases, AGS4 and AGS15, there was a clear association in the H-band image that was matching the ALMA contours. In both cases, the H-band detection was close to a bright H-band neighbor and both objects were interpreted in the CANDELS catalog as a single object, hence neglecting the ALMA source. If we were to extrapolate our small number statistics to the population of optically dark galaxies at these depths, this would imply that about 35% of the optically dark sources are mistaken as such due to confusion. We note, however, that these galaxies are detected here because of the very deep H-band image since their AB magnitudes are H = 25.23 and 27.11 AB for AGS4, and AGS15, respectively.

All four of the optically dark galaxies detected from a blind detection above 4.8σ are detected with IRAC and listed in the catalog of Ashby et al. (2015). These galaxies are AGS4, 11, 15 and 17. The remaining two optically dark galaxies, AGS24 and AGS25, were detected in the ALMA image using IRAC priors in F20a. However, they were not listed in the catalog of Ashby et al. (2015). In both cases, F20a (see also Sect. 5.3.1) noticed that despite being well above the IRAC detection limit, these two galaxies were missed because of the presence of a bright neighbor due to confusion (see Fig. 5). After carefully subtracting the surrounding neighbors, the flux densities measured at 3.5 μm for both galaxies are F_3.6 μm(AGS24) = 4.7 ± 0.5 μJy and F_3.6 μm(AGS25) = 1.5 ± 0.1 μJy.

Fig. 5. Images of AGS24 in the B, V, H and Ks band, IRAC 3.6 and 4.5 μm, JVLA 6 GHz (3.7σ) and 3 GHz (5.7σ). The contours are the same as in Fig. 6 but for AGS24. The three neighbor galaxies are labeled with IDZFOURGE or IDCDFS and their redshift separately.

Five of the GOODS-ALMA optically dark galaxies (AGS4, 11, 15, 17 and 24) were followed up with ALMA band 4 (see Sect. 2.1.2) and all are confirmed with an S/N of 14.7, 16.8, 15.3, 15.9 and 8.7 in the 2 mm continuum. AGS25 was not included in the ALMA proposal because the list of sources below the 4.8σ limit that includes it was still under construction at the time. All the six optically dark galaxies are also detected above 4.5σ at 870 μm continuum (Cowie et al. 2018). Hence there is no doubt that they are all real.

Another confirmation of the robustness of three of the optically dark galaxies – AGS4, 17 and 24– comes from their detection in the radio with the VLA at 3 GHz with an S/N of 10.4, 9.9, and 5.6, respectively. AGS11, 15 and 25 are not detected with S/N higher than 5 (Rujopakarn et al., in prep.). The first two of the detected galaxies, AGS4 and AGS17, are the only optically dark galaxies in the sample that are also detected with Herschel and their far-infrared SED is well-measured with no less than five Herschel data points from 100 to 500 μm in the observed frame. We discuss in the following two sections the properties of AGS4 and AGS25.

4.1. AGS4, an extremely massive galaxy at z = 3.556 and a case of blending in the Hubble H-band image

AGS4 is detected in GOODS-ALMA with an S/N = 9.7 at 1.1 mm (F18, see also red solid contours in Fig. 6). Using the ALMA follow-up observations presented in Sect. 3.1, we now have a confirmation of the robustness of this source through its detection with ALMA at 2 mm with an even higher S/N = 14.7 (green dashed contours in Fig. 6).

Fig. 6. Images of AGS4 from the observed B (435 nm, rest-frame UV) to near-infrared observed Ks (2.2 μm, rest-frame B) bands. The red contours denote the ALMA detection at 1.1 mm at the resolution of 029. The green dashed contours denote the ALMA detection at 2 mm at the resolution of 088. The astrometric corrections between ALMA and HST images are applied throughout this work. The cyan “x” denotes the optical bright neighbor, IDCANDELS = 8923 (zphot = 0.24). The magenta “x” denotes the detection by ZFOURGE based on the Ks-band image, IDZFOURGE = 12 333 (zphot = 3.76). If not stated otherwise, the images in this paper are all oriented with the north up and east to the left.

AGS4 is an illustration of what we may call spectral de-confusion. The two catalogs CANDELS and ZFOURGE have identified a single object associated with this position on the sky. The CANDELS position is close to the center of the bright optical emitter whereas the ZFOURGE position is close to the bright K_s-band near-infrared emitter. The position of the ALMA detection exhibits a clear offset with respect to the bright optical emitter and falls on the northwestblob that is the brightest in the near-infrared. As shown in Fig. 6, AGS4 is at 05 from a bright optical emitter. They are considered as two regions of the same galaxy (ID_CANDELS = 8923, z_phot = 0.24, cyan “x” in Fig. 6) in the CANDELS catalog (Guo et al. 2013). The ZFOURGE position is at 016 (ID_ZFOURGE = 12 333) to AGS4 (magenta “x” in Fig. 6), and the photometic redshift derived by the ZFOURGE team is = 3.76.

As discussed above, AGS4 exhibits an H-band counterpart well above the WFC3 detection limit in GOODS-South of H = 25.23 AB that was missed due to blending with a bright neighbor. Therefore, despite being truly optically dark in the optical bands, AGS4 is not strictly speaking an optically dark galaxy since it is detected in the H-band. After de-blending³ both emitters – none of which corresponds to the centroid of the CANDELS or ZFOURGE positions that are both in between the two emitters – we can clearly see two different SEDs (Fig. 2-top-left). We note that the two ALMA detections at 1.1 mm and 2 mm are well-centered on the northwestemitter. The photometric redshift associated to this object, once deconvolved from the bright optical emitter, is = 4.02. We note that the optical source identified in the CANDELS catalog (ID_CANDELS = 8923) is clearly offset from the ALMA contours, leaving little doubt that these are two separate galaxies. It is only in the low resolution K_s-band image that some confusion takes place. Being near-IR based, the ZFOURGE catalog identified the position of the peak K_s emission that matches the ALMA galaxy and attributed to it a photometric redshift of = 3.76 (ID_ZFOURGE = 12 333). The SED fittings of the photometric points associated to the CANDELS (gray) and ZFOURGE (green) galaxies are both shown in Fig. 2-top-left. The gray SED does not leave any doubt on the fact that this is a very nearby galaxy ( = 0.24), with a clear drop of the emission above 0.8 μm, whereas the green SED exhibits a break around 1.6 μm ( = 3.76). The very strong difference between the SEDs of the CANDELS and ZFOURGE sources clearly show that they cannot be considered as two sides of the same galaxy, one visible in the optical the other one in the near to far infrared. These two redshifts suggest that there are indeed two galaxies in projection here.

Our ALMA spectral-scan follow-up presented in Sect. 3.1 provides complementary evidence for this spectral de-confusion between two super-imposed galaxies with a clear line detection that matches the CO(6-5) transition at a redshift of z_spec = 3.566 (see Fig. 1-left), well within the range of possibilities of our photometric redshift estimate. Adopting the redshift at = 3.556, the stellar mass of AGS4 is M_⋆ = 10^{11.09 ± 0.08} M_⊙, which makes AGS4 one of the most massive galaxies in the early Universe. We note that adopting = 4.326, as for CO(7-6), the stellar mass that we obtain for AGS4 is M_⋆ = 10^{11.45 ± 0.20} M_⊙, which would make AGS4 the most massive galaxies known in the Universe above z = 4.

We investigated the possibility that the stellar mass of AGS4 may be contaminated by the presence of the continuum from an AGN, but we did not find any evidence of a power-law continuum emission in the SED presented in the Fig. 23 of F18 (also in F20b), nor any radio emission associated to it. There is an X-ray detected AGN that is associated with its neighbor (ID_ZFOURGE = 12 333), which could potentially affect the measurement if attributed to the wrong counterpart, but again we do not see any clear signature of such effect. The q_TIR-parameter (see Eq. (2)) measuring the infrared to radio flux ratio parameter is equal to q_TIR = 2.49 ± 0.06, higher than the typical value for star-forming galaxies at z ∼ 3.5 proposed by Delhaize et al. (2017), q_TIR = 2.16  ±  0.06, indicating no radio excess of this galaxy.

(2)

where L_TIR is total infrared luminosity (8–1000 μm) in unit W, and the L_{1.4 GHz} = 2.73 × 10²⁴ W Hz⁻¹ is derived from the 10.4σ detection of AGS4 at 3 GHz assuming a radio spectral index of α = = −0.8, the typical value for galaxies at z > 2 (Delhaize et al. 2017), using Eq. (1).

4.2. AGS25, the most distant optically dark galaxy in GOODS-ALMA

AGS25 is the only optically dark galaxy in the sample that is identified in the ZFOURGE catalog and is not confused with a neighboring galaxy. The counterpart of AGS25 in the K_s-band image is clearly visible in Fig. 7 (middle panel). It is associated with the ZFOURGE galaxy ID_ZFOURGE = 11 353 at = 4.64. We fitted its photometric measurements as listed in the ZFOURGE catalog using the code EAzY for consistency with the other galaxies, and found a PDF of the photometric redshift well peaked at the redshift given in the ZFOURGE catalog (Fig. 2-bottom-right). We find a peak redshift of = 4.71 that encompasses the value from ZFOURGE that we keep for consistency with F20a. This makes AGS25 the most distant optically dark galaxy in GOODS-ALMA. The IRAC emission of the galaxy is polluted by the contribution of a bright neighbor at a distance of 3 arcsec (outside of the postage stamps shown in Fig. 7). This explains why the galaxy is not listed in the IRAC S-CANDELS catalog (Ashby et al. 2015). After de-convolution of the bright neighbor, a clear IRAC detection is obtained (see Fig. 7 right panel) that was used to determine its flux density in F20a. The stellar mass of AGS25 at this redshift is M_⋆ = 10 M_⊙ (F20a). AGS25 is not detected by Herschel, hence its SFR = 839 ±131 M_⊙ yr⁻¹ was derived in F20a by adjusting the SED of Schreiber et al. (2018a) to the 1.1 mm ALMA flux density corresponding to 195 μm in the rest-frame, which is not far from the peak of the far-infrared SED.

Fig. 7. Images of AGS25 in the H-band (1.6 μm), Ks-band (2.2 μm) and IRAC1-band (3.6 μm). The contribution of a 3″distant bright neighbor of AGS25 was subtracted to produce the IRAC image as in F20a. The red contours denote the ALMA detection at 1.1 mm at the resolution of 060. The magenta “x” is the ZFOURGE counterpart IDZFOURGE = 11 353.

5. An overdensity at z ∼ 3.5 in GOODS-ALMA

5.1. Clustering properties of optically dark galaxies

A large statistical sample of IRAC detected optically dark galaxies in the CANDELS fields studied by Wang et al. (2019) shows that the optically dark galaxies are the most massive galaxies in the early Universe and are probably the progenitors of the most massive galaxies we see in the groups or clusters of galaxy in the local Universe.

Here in the 69 arcmin² GOODS-ALMA field, four of the six optically dark galaxies reside in a small region of 5 arcmin² only. We calculate that the probability to have four out six galaxies randomly falling in such a small area is only 0.4% from four million Monte Carlo mock realizations. In addition, all the four optically dark galaxies present a redshift consistent with z = 3.5 where field galaxies demonstrate a clear peak in their redshift distribution.

5.2. A clear peak at z ∼ 3.5 in the redshift distribution

The redshift distribution of the galaxies with either a photometric or a spectroscopic redshift in the 6.9′×10′ GOODS-ALMA field presents a peak at z ∼ 3.5 (Fig. 8). This peak is at a 3.5σ significance among the 1373 galaxies at 3 ≤ z ≤ 4 ( Fig. 8). The redshift bin size is taken to be 0.03 to optimize the peak in the redshift distribution that was previously identified in Straatman et al. (2016, Fig. 23) using the ZFOURGE catalog.

Fig. 8. Redshift distribution of the 11 674 galaxies with a redshift (either photometric or spectroscopic) located within the 6.9′ × 10′ GOODS-ALMA field (gray). The subsample of 395 galaxies located within 50″ of the most massive galaxy AGS24, are shown in blue. The insert shows the distribution of the spectroscopic redshifts of the 121 galaxies at 3 ≤ zspec ≤ 4 over the GOODS-ALMA field. Ten galaxies fall in the redshift bin 3.470 ≤ zspec ≤ 3.480. The mean redshift of the ten galaxies is zspec = 3.472. The redshifts are from the ZFOURGE (Straatman et al. 2016) and VANDELS DR3 (McLure et al. 2018; Pentericci et al. 2018) catalogs.

Here we increase by a factor of three the total number of galaxies with a spectroscopic redshift that fall within the redshift range 3 ≤ z ≤ 4 in GOODS-ALMA. As many as 83 new spectroscopic redshifts have recently been measured and released in the data release 3 (DR3) of the VANDELS survey (McLure et al. 2018; Pentericci et al. 2018) that fall within 3 ≤ z ≤ 4 in GOODS-ALMA. The VANDELS redshifts supplement the already existing 38 spectroscopic redshifts originally listed in the ZFOURGE catalog and used by Straatman et al. (2016).

The 121 galaxies with spectroscopic redshifts at 3 ≤ z_spec ≤ 4 show a clear peak at 3.470 ≤ z_spec ≤ 3.480 containing ten galaxies (Fig. 8). These ten galaxies have a mean redshift at z_spec = 3.472. This is consistent with the peak at z ∼ 3.5 shown by the 1373 galaxies with either spectroscopic or photometric redshifts.

5.3. optically dark galaxies at z ∼ 3.5

We find that four out of the six optically dark galaxies exhibit a redshift that is consistent with being located in the same redshift peak at z ∼ 3.5, which corresponds to an overdensity centered on the most massive of them, AGS24, as to be discussed below. We may even count five out of the six if we include AGS4, which is at z_spec = 3.556 and at a distance of 3′ from AGS24. Accounting for the redshifts of the sources, this corresponds to 70.8 cMpc (comoving) and 15.5 pMpc (proper). The comoving distance is typical of those found in distant super-clusters while the proper distance is consistent with those of proto-clusters, which suggests that AGS4 could potentially be a member of the proto-cluster as well. We start by describing in detail the properties of AGS24, because it is the most massive of all the optically dark galaxies in GOODS-ALMA and even among all GOODS-South galaxies at z > 3 and constitutes an excellent candidate for a future BCG galaxy. The remaining three optically dark galaxies at z ∼ 3.5 are discussed in the next section.

5.3.1. AGS24, the most massive galaxy at z > 3 in GOODS-South

AGS24 has been detected with an S/N of 4.9σ and 8.7σ in the GOODS-ALMA 1.1 mm image and in the 2 mm continuum image of the ALMA spectroscopic scan follow-up (red and green contours in Fig. 5). This galaxy exhibits a lower S/N in the GOODS-ALMA image tapered at 06 (3.9σ) than in the image at the original resolution of 029 (4.9σ), which suggests that the dust emission is very compact. This optically dark galaxy exhibits a clear radio counterpart at 6 GHz (3.7σ) and 3 GHz (5.7σ). Despite being absent from the list of IRAC sources in the field given in the S-CANDELS catalog (Ashby et al. 2015), this source does exhibit clear IRAC emission at 3.6 and 4.5 μm that is partly blended with the three surrounding very nearby galaxies that are clearly detected in the H and K_s bands. The emission of AGS24 in the H and K_s bands is marginal but was estimated after de-blending the contribution from the three neighboring sources and used to determine its photometric redshift.

After de-blending³, we derive a photometric redshift of z_AGS24 = 3.58. Assuming this redshift, the excess emission in the V band may correspond to the Lyα emission line (Fig. 2-bottom-left). However, we analyzed the VLT/MUSE data obtained in this region but did not find any evidence for Lyα emission despite the high sensitivity of MUSE (rms = 4 × 10²⁰ erg s⁻¹ cm⁻² A⁻¹, Herenz et al. 2017). We conclude that the excess V band emission remains uncertain and may be due to noise fluctuation.

The three neighboring galaxies are located at similar distances of ∼15 (equivalent to 11.3 kpc at z ∼ 3.5) from AGS24 (Fig. 5 and Table 2). They all have a redshift of z∼3.5, including two spectroscopic redshifts – ID_ZFOURGE = 10627, z_spec = 3.478 and ID_CDFS = 007493, z_spec = 3.482 (VANDELS DR3) – and one photometric redshift – ID_ZFOURGE = 10 672, z_phot = 3.53 (Straatman et al. 2016). As a result, the group of four galaxies appears to be physically connected and makes a dense and compact association of galaxies typical of what can be found at the center of distant proto-clusters (Springel et al. 2005; De Lucia & Blaizot 2007).

Altogether, the following reasons suggest that AGS24 is located at the redshift of z ∼ 3.472 (Sect. 5.2) of the proto-cluster: (i) the three neighbors of AGS24 have z ∼ 3.472, (ii) the photometric redshift of AGS24 is consistent with this redshift, (iii) AGS24 belongs to the highly improbable association of optically dark galaxies that fall in the region of the z ∼ 3.472 proto-cluster.

Assuming this redshift, we estimate the total infrared luminosity of AGS24 to be L_IR = 10^{12.31 ± 0.11} L_⊙ fitting the rest-frame 195, 245, 490 μm emission of the galaxy that corresponds to the observed 870 μm, 1.13 and 2 mm using the typical far-infrared SED of a MS galaxy at this redshift (Schreiber et al. 2015). We also checked that we obtained a consistent total infrared luminosity after fitting these measurements with a Draine et al. (2007) SED, which gives the same result within 10%.

Assuming z = 3.472, the stellar mass of AGS24 is M_⋆ = 10 M_⊙ (adopting the Calzetti et al. 2000 attenuation law and a delayed exponentially declining SFH). It is ∼10 times more massive than the galaxy presented in Ginolfi et al. (2017), which used to be the most massive galaxy at z = 3–4 in this field. We note that the stellar mass was derived by fitting the near UV to near-infrared photometry, similar to the method we used for our galaxies. In order to test the robustness of the value of the stellar mass of AGS24, we recompute it assuming various star formation histories (SFH) to see whether this mass could fall down by a large amount. We use the SED fitting code FAST++⁴ to test a delayed exponentially declining and an exponentially declining SFH as well as a truncated SFH, with the possibility to add a second burst. We also test the impact of varying the attenuation law using the code CIGALE (Boquien et al. 2019) with the attenuation laws from Calzetti et al. (2000), Charlot & Fall (2000), and the modified Charlot & Fall (2000) law in MAGPHYS (da Cunha et al. 2008). The results are listed in Table 3. They show that despite these wide variety of SFH and attenuation laws, the stellar mass of AGS24 remains within 10^11.31 and 10^11.52 M_⊙.

With a redshift of z ∼ 3.5, AGS24 is the most massive galaxy in the GOODS-ALMA field above z = 3. The other similarly massive candidate galaxies in this redshift range either have a strong AGN potentially contaminating the estimate of the stellar mass or highly unreliable redshifts. Two galaxies without an AGN at z > 3 have similarly large stellar masses in the GOODS-ALMA field (ID_ZFOURGE = 11 505 and 16 410, marginally detected in the K_s-band with S/N = 5.0, 5.5, at z_phot = 5.36, 4.73). They both have very high χ² values of 126.6 and 43.0 associated with their photometric redshifts. Only 3% and 19% of the ZFOURGE galaxies have photometric redshifts with such high χ² values. Therefore, we consider the redshifts of these two galaxies not reliable.

We inspect the possible contamination of an AGN on the near- and mid-infrared emission, which may result in an overestimation of the stellar mass. We calculate the infrared to radio flux ratio, q_TIR = 2.05 ± 0.18 (see Eq. (2)), where the 1.4 GHz luminosity, L_{1.4 GHz} = 1.87 × 10²⁴ W Hz⁻¹, is converted from the JVLA 3 GHz emission (5.7σ) of AGS24 assuming a radio spectral index of α = −0.8, as we did for AGS4. Within the rms, the infrared to radio flux ratio agrees with the value q_TIR = 2.16  ±  0.06 for star-forming galaxies at redshift z ∼ 3.5 by Delhaize et al. (2017), indicating no radio excess of this galaxy. We also checked that AGS24 does not have a counterpart in the deepest Chandra 7Ms source catalog (Luo et al. 2017). We searched for extended diffuse emission in the Chandra 7Ms image but found no robust evidence for such emission that could be associated with an intra-cluster medium. Therefore, we conclude that AGS24 is the most massive galaxy at z > 3 without an AGN in the GOODS-ALMA field. Forrest et al. (2020) recently spectroscopically confirmed (z_spec = 3.493) a similarly extremely massive quiescent galaxy in the VIDEO XMM-Newton field (Annunziatella et al., in prep.) that is spectroscopically confirmed (z_spec = 3.493) and interpreted as experiencing a rapid downfall of star formation.

The predicted number of galaxies that are as massive or more massive than AGS24 in the GOODS-ALMA field at 3 < z < 4.5 is 1.15, based on the galaxy stellar mass function (SMF) derived from the observations in the COSMOS field (Davidzon et al. 2017). Pillepich et al. (2018) compared the SMF of Davidzon et al. (2017) to the one from the IllustrisTNG simulation and found them consistent. Thus the existence of one massive galaxy, AGS24, is consistent with the existing observation and simulation.

5.3.2. Physical properties of the optically dark galaxies AGS11, AGS15 and AGS17 at z ∼ 3.5

AGS11. AGS11 is detected at 1.1 mm with an S/N = 5.71 in F18 (also see red solid contours in Fig. 4-top) and 2 mm with an S/N = 16.8 (green dashed contours in Fig. 4-top). A neighboring galaxy is identified in ZFOURGE that matches the K_s and IRAC emission (magenta “x” in Fig. 4-top). This object is located at 45 from the ALMA detection. It has been attributed with a photometric redshift of = 4.82 in the ZFOURGE catalog. However, this redshift is highly uncertain with flag use = 0, hence highly unreliable. Data at hand do not allow us to determine whether this object is the counterpart of the ALMA detection, AGS11, or like several other optically dark sources studied here, an independent galaxy seen in projection. The optical emission at the ALMA detection position is within the noise fluctuation. Therefore, we do not attributed the ZFOURGE ID and the associated redshift to AGS11. Instead, we consider the association of the optically dark galaxies with the overdensity at z ∼ 3.5 as potential evidence that AGS11 has a redshift matching this overdensity. In Table 1 we decided to attribute AGS11 with the z_peak = 3.472 of the overdensity. Assuming z_AGS11 = 3.472, the infrared luminosity we derived for AGS11 is L_IR = 10^{12.94 ± 0.07} L_⊙.

AGS15. AGS15 is located at 059 to a optically bright galaxy (ID_CANDELS = 3818, cyan “x” in Fig. 9). It is detected at 5.22σ in the 1 mm image (F18, see also red contours in Fig. 9). This detection is further confirmed in the 2 mm continuum at 15.3σ in the follow up ALMA spectroscopic scan (green dashed contours in Fig. 9). In the K_s-band image, the infrared emission from AGS15 and from the optically bright galaxy is blended, and the two galaxies are confused as one (ID_ZFOURGE = 6755, magenta “x” in Fig. 9) in the ZFOURGE catalog. The redshift derived by the ZFOURGE team is = 3.46 and the stellar mass is M_⋆ = 10^9.86 M_⊙. We notice that there is a high uncertainty in the SED fitting by the ZFOURGE team, with χ² = 45.4, which is higher than 84% of all the galaxies in the ZFOURGE catalog. If it is the counterpart of AGS15, then the low stellar mass will make AGS15 an extreme case, that is, ten times less massive compared to the other ALMA detections. After de-blending AGS15 and the optically bright galaxy, as shown in Figs. 9 and 2-top-right, we find that in fact AGS15 has an H band magnitude higher than the detection limit, H = 27.11 AB. We fit the SED with EAzY (Brammer et al. 2008) and FAST++⁴ to derive the photometric redshift and the stellar mass. The z_peak derived for AGS15 is at = 3.96. Hence the galaxy has, within the uncertainties on the photometry, a photometric redshift that encompasses the redshift peak at z_peak = 3.472. Assuming this redshift, we obtain the infrared luminosity of AGS15 to be L_IR = 10^{12.78 ± 0.08} L_⊙, the SFR to be SFR = 1042 M_⊙ yr⁻¹ and the stellar mass to be M_⋆ = 10 M_⊙.

Fig. 9. Same as in Fig. 6, but for AGS15. The cyan “x” denotes the optical bright neighbor, IDCANDELS = 3818 (zphot = 4.27). The magenta “x” denotes the detection in Ks-band by ZFOURGE, IDZFOURGE = 6755 (zphot = 3.47).

AGS17. AGS17 is detected at 1.1 mm with an S/N = 5.01 (F18, also see red solid contours in Fig. 4-bottom) and with an S/N = 15.9 in the 2 mm continuum (green dashed contours in Fig. 4-bottom). It is located between two bright optical emitters, ID_CANDELS = 4414, ( = 0.03, Pannella et al. 2015) at 027, and ID_CANDELS = 4436, ( = 0.95, fitted with EAzY) at 057 (cyan “x” s in Fig. 4-bottom). These two emitters are also clearly seen in the K_s-band image, however, they are not resolved in the ZFOURGE catalog (magenta “x” in Fig. 4-bottom, ID_ZFOURGE = 6964). The redshift derived from this blended emission is = 1.85, which is not consistent with either the redshifts of the optical emitters.

Assuming that the emission line detected in our ALMA spectroscopic follow-up is the CO(6-5) transition as discussed in Sect. 3.2, AGS17 has a spectroscopic redshift = 3.467, which puts it at a distance of 4.30 cMpc from AGS24, assuming z^AGS24 = z_peak = 3.472. Hence it falls inside the proto-cluster as the distance is shorter than the typical size of a proto-cluster (∼14.5 cMpc) at z ∼ 3.5 in simulations (Muldrew et al. 2015). The low probability to have a galaxy fall in projection with the structure that is located at the redshift corresponding to the CO(6-5) transition makes us favor this transition. Then the stellar mass is M_⋆ = 10 M_⊙ based on the IRAC band emission and the infrared luminosity derived from the SED fitting using Herschel and ALMA data is L_IR = 10^{13.08 ± 0.02} L_⊙.

5.4. Spatial distribution of galaxies at z ∼ 3.5 in the GOODS-ALMA field

5.4.1. Local overdensity (δ) of individual galaxies

We calculate the overdensity (δ = ρ / ) of the 364 galaxies at 3.42 ≤ z ≤ 3.57 in GOODS-ALMA. This redshift range covers a distance of 125 cMpc, which is consistent with the size of a large scale structure, a (proto-)supercluster of ∼ 60 × 60 × 150 cMpc at z ∼ 2.45, identified by Cucciati et al. (2018). We derive a dispersion of σ(δ) = 0.20 among the 364 galaxies. We do that by comparing the projected local galaxy density within 2 cMpc (ρ) of a galaxy, to the average number density over the GOODS-ALMA field ( = 1.32 cMpc⁻² or 5.13 arcmin⁻²).

Seven of the detections in the GOODS-ALMA field reside in this narrow redshift range, AGS1, 4, 5, 11, 15, 17, and 24, five (underlined) out of the seven are optically dark. The = 3.442 of AGS1 comes from VANDELS DR3, zflag = 1, meaning 50% probability to be correct, while we note that it is recorded as = 2.309 in F18. AGS5 has a photometric redshift of = 3.46 (Straatman et al. 2016). As mentioned in Sect. 4.1, AGS4 has a redshift of = 3.556.

At the same time, eight massive galaxies fall in this redshift range in the ZFOURGE catalog. They have a stellar mass, log [M_⋆/M_⊙] > 10.5, which is equivalently massive as the ALMA detections (F18). They have a median stellar mass of M_⋆ = 10^10.82 M_⊙. These eight massive galaxies, as well as the seven ALMA detections all show a local overdensity (δ > 1). Out of the 15 massive galaxies, AGS24 shows the most prominent clustering, δ_AGS24 = 2.06, hence a > 5σ excess with respect to the average galaxy surface density. Hence AGS24 is both the most massvie galaxy at z > 3 in the GOODS-ALMA field with no AGN and it is located in the densest surface density peak at z ∼ 3.5. This makes it an excellent candidate for a future BCG as we discuss in Sect. 5.5. The degree of overdensity of the four other HST-dark galaxies are δ_AGS11 = 1.46 (> 2σ excess), δ_AGS15 = 1.64 (> 3σ excess), δ_AGS17 = 1.76 (∼4σ excess), δ_AGS4 = 1.39 (∼2σ excess), and for the other two ALMA detection with optical counterparts, it is δ_AGS1 = 1.21, δ_AGS5 = 1.03.

Regardless of mass, when focusing on the galaxies in the densest regions, that is to say, δ > 5σ(δ), we found in total 12 galaxies in the ZFOURGE catalog. It turns out that they all lie within 50″of AGS24 (red circle in Fig. 10) and the barycenter of these galaxies falls exactly at AGS24. The most massive of the 12 galaxies is ID_ZFOURGE = 10672, M_⋆ = 10^10.76 M_⊙. It is at 15 from AGS24.

Fig. 10. Projected density of 364 galaxies (40 with zspec) at 3.42 ≤ z ≤ 3.57. The contours show 1, 2, 3, 4σ level of the overdensity. Left panel: field observed in GOODS-South ALMA. Right panel: zoom-in of the extended structure in the bottom right corner of GOODS-ALMA field that encompasses 126(18) galaxies over 15 arcmin2. The red star indicates the position of AGS24. The red circle indicates a 50″radius region (2 arcmin2) including 26(3) galaxies. The yellow stars refer to the HST-galaxies AGS 11, 15 and 17 located also in this larger structure. The green circles represent the AGS 1 and 5, which are also in this redshift range but not optically dark. The white circles show the 40 galaxies with spectroscopic redshifts. We note that the two white circles falling on the red star of AGS24 are two galaxies dissociated from AGS24. The magenta square denotes the center of the proto-cluster CCPC-z34-002 identified by Franck & McGaugh (2016). The cyan square denotes Candels5001, which used to be the most massive galaxy at 3 <  z < 4 in this field (Ginolfi et al. 2017). The 50″scale in the figure is equivalent to a proper distance of 0.37 Mpc at z = 3.5.

In Fig. 11, we show the overdensity as a function of stellar mass of the 364 galaxies at 3.42 < z < 3.57 in the GOODS-ALMA field and the ALMA detected galaxies. We can see that as the stellar mass increases, the galaxies tend to fall in denser regions, which is in line with what Muldrew et al. (2015) found in Millennium Simulation (Springel et al. 2005) on galaxies at z = 2. They conclude that the most massive galaxies at high redshift reside in proto-clusters given the little difference in the environments of only proto-cluster galaxies and those of all the massive galaxies. Besides, the ALMA detections are the most massive ones in the field and reside in the densest regions. AGS24, as the most massive one in the field, resides in the densest region in the field. Therefore, these optically dark galaxies revealed by ALMA make up for what has been missed from optical observations. The relation between their dusty and massive nature and the dense environment they reside in provides evidence of environmental effects on galaxy evolution.

Fig. 11. Local overdensity (δ) of individual galaxies in the GOODS-ALMA field at 3.42 < z < 3.57. The black dots show the sliding medians of δ at different stellar masses. The filled circles indicate all the ALMA detections. The purple circles are the HST dark galaxies and the red circle indicates AGS24.

5.4.2. Projected density map over the GOODS-ALMA field

We construct the projected galaxy density in Fig. 10. Similar to the method used in Forrest et al. (2017), we create a grid with a cell size of 09 × 09 in the GOODS-ALMA field and calculate the overdensity (δ = ρ / ) for each cell. 364 galaxies (40 with z_spec) at 3.42 ≤ z ≤ 3.57 are used to produce the map. The galaxies with spectroscopic redshifts are marked as white dots in Fig. 10. An overdense region comes out in the southeast region as indicated by the zoom-in box in Fig. 10. The 4σ peak of this overdensity happens to be very close to AGS24, the galaxy with the strongest overdensity (δ_AGS24 = 2.06) this is consistent with the fact that it is in a highly clustered environment. Out of the 364/40 (redshift/z_spec) galaxies in the 69 arcmin² field, 126/18 reside in this 15 arcmin² box, and 26/3 galaxies in this redshift slice fall within 50″ (2 arcmin²) of AGS24.

We notice that the overdensity shown in the projected density map encompasses four out of the six HST-dark galaxies detected by ALMA (AGS 11, 15, 17, 24, stars in Fig. 10). The overlap of an overdensity of ALMA detected optically dark galaxies and an overdensity of optically detected galaxies shows that these ALMA detected galaxies are good tracers of structures in the early Universe.

This overdensity has been classified as a proto-cluster in previous publications. Franck & McGaugh (2016) identified it in a list of proto-clusters at 2.74 < z < 3.71 built out of ∼14 000 spectroscopic redshifts (with 604 of them finally used). They defined proto-clusters as galaxy densities greater than the field density by a factor seven within a radius of 2′ and a redshift range of δz < 0.03. One of their proto-clusters, CCPC-z34-002, has a redshift of z = 3.476, hence consistent with the redshift peak of z = 3.472 studied here. We note that the CCPC-z34-002 density peak, roughly estimated from the distribution of only 23 galaxies within a radius of 20 cMpc, does not correspond to the location of AGS24 but falls only 10″ from AGS17, and the typical diameter of the proto-clusters studied by Franck & McGaugh (2016) of 4′ matches the zoom shown in Fig. 10-right. In the following, we will consider the galaxies at the redshift peak that fall in this area as a member of a candidate proto-cluster. This structure encompasses four out of six, hence 67%, of the optically dark galaxies identified in the GOODS-ALMA field.

This overdensity was also identified by Forrest et al. (2017) who used a population of extreme emission line galaxies, with low stellar masses and by the VANDELS survey with a three-dimensional algorithm (Guaita et al. 2020). Moreover, only nine out of 131 the Lyα-emitting galaxies selected from the VANDELS survey reside in their 13 detected overdensities at z > 2, indicating that Lyα-emitting galaxies are not ideal tracers of overdensities. Finally, we note that Ginolfi et al. (2017) observed the molecular gas content of the galaxy CANDELS-5001 (blue square in Fig. 10, M_⋆ = 10 M_⊙, z_spec = 3.473) and found that this galaxy was surrounded by an extended emission of molecular gas spanning 40 kpc traced by the CO(4-3) transition. Furthermore, they detect nine additional CO systems within a radius of 250 kpc from the massive galaxy and mostly distributed in the same direction as the CO elongated structure found in the central 40 kpc, which they interpret as evidence for large-scale gas accretion on the galaxy. This large-scale accretion might be enhanced by the larger scale dark matter halo of the structure (see Rosdahl & Blaizot 2012).

We calculate the distance between AGS1, 4, 5 and AGS24 to be 40.5 cMpc, 70.8 cMpc, and 25.5 cMpc, respectively. Despite their larger distance, AGS1 and 5 might be involved in a wider structure than the proto-cluster that includes AGS11, 15, 17, 24 and potentially also AGS4. This wider structure may be compared to the superclusters at high redshift discussed in Cucciati et al. (2018) and Toshikawa et al. (2020).

5.4.3. Two point correlation function (2PCF)

We calculated the auto-correlation function of the 364 galaxies (40 with z_spec) at 3.42 ≤ z ≤ 3.57 using a uniformly generated random sample in the GOODS-ALMA field, independent of mass (red dots in Fig. 12), the cross-correlation function between these 364 galaxies and the ALMA detections (black triangles in Fig. 12), and the cross-correlation function between these 364 galaxies and the massive galaxies that are equivalently massive as the ALMA detections (M_⋆ > 10^10.5 M_⊙) out of the 364, (green triangles in Fig. 12). A Landy & Szalay (1993) estimator is adopted. Error bars are estimated from 1000 bootstrap samples. Despite of the large error bars, the two cross-correlations both show systematically higher excess possibility compared to the auto-correlation of the 364 galaxies in this redshift range. This indicates a clustering feature of the ALMA detected optically dark galaxies, as well as the equivalently massive galaxies.

Fig. 12. Comparison of the auto-correlation function (ACF) of the 364 galaxies at 3.42 ≤ z ≤ 3.57 in the ZFOURGE catalog (red dots), the cross-correlation function (CCF) between these galaxies and the ALMA detections (black triangles), and the CCF between these galaxies and the galaxies that are equivalently massive (M⋆ > 1010.5 M⊙) as the ALMA detections (green triangles).

5.5. Dynamical state of the proto-cluster at z ∼ 3.5

As discussed in Sect. 5.3.1, AGS24 is the most massive galaxy in the field at z > 3 (M_⋆ = 10 M_⊙) and in particular at z ∼ 3.5, making it a good candidate for the future BCG of the overdensity of galaxies if it becomes a galaxy cluster. The fact that it is surrounded by three galaxies with the same redshift of z ∼ 3.5 within 15 from AGS24 (11.3 kpc at z ∼ 3.5 ) strengthens this hypothesis (Fig. 5 and Table 2).

To estimate the dark matter mass of this structure, we first calculated the total stellar mass of all the galaxies within 50″of AGS24 and then determined an upper limit of M_h = 10^15.0 M_⊙ to the dark matter mass based on the relation in Behroozi et al. (2013), assuming that this structure is already virialized. We also determined a lower limit of M_h = 10^13.3 M_⊙ by summing up the individual halo masses of the each galaxy within 50″of AGS24.

We note that the projected number density profile of the galaxies in the structure do not follow the slope expected for the NFW profile (Navarro et al. 1996) of a gravitationally bound cluster, but instead exhibit a shallow slope (Fig. 13). We show for comparison the projected number density profile of the most distant galaxy cluster to our knowledge, CL J1001+0220 at z = 2.506 (Wang et al. 2016).

Fig. 13. Projected numbers of galaxies in the overdensity as a function of their distance to AGS24 (black circles) and the black dashed line is the linear fitting to the profile. The red line shows the best-fit projected NFW profile (Navarro et al. 1996) of a galaxy cluster at z = 2.506, CL J1001+0220 (Wang et al. 2016).

However, the structure shown in the density map (Fig. 10), as well as the fact that AGS24 is the most massive galaxy that falls at the density peak and has a reduced SFR compared to MS galaxies, suggests that this structure is in the process of viralization and that such an environment is influencing the evolution of AGS24, diminishing its star formation activity. Galaxies like AGS24 have been missed by previous studies given their dusty nature, despite being the most massive candidate galaxies that can be used as a beacon to trace the barycenter of the proto-clusters in the process of virialization.

6. Conclusions

We have analyzed the properties of an ensemble of six optically dark galaxies, AGS4, 11, 15, 17, 24 and 25 identified in the GOODS-ALMA survey. They do not have optical counterparts in the deepest H-band based catalog down to H = 28.16 AB. Out of the six sources initially classified as optically dark, a careful analysis revealed that two were associated with an H-band counterpart that was not identified in the H-band catalog missed due to blending. AGS4 and AGS15, that were identified as optically dark in F18, are extremely close to a bright optical source at 050 and 027, respectively, and as such they were wrongly considered as parts of the bright neighbor. However, we showed that their SED can be used to perform a spectral de-confusion of the sources at very different redshifts. We note that the H-band magnitude of these galaxies is extremely faint since it is 25.2 AB and 27.1 AB for AGS4 and AGS15, respectively. The second one is below the detection limit of most H-band images. The definition of HST-dark or optically dark is obviously a function of the depth of the HST images. All of six sources except AGS24 do exhibit a counterpart in the K_s-band image but they were mistakenly associated with a bright neighbor, hence none was listed in the ZFOURGE catalog.

We performed a spectroscopic scan follow-up of five of them with ALMA (AGS4, 11, 15, 17 and 24). All five exhibit a clear 2 mm counterpart in the continuum with S/N = 14.7, 16.8, 15.3, 15.9, and 8.7, respectively, which reinforces their robustness. We detected one emission line for AGS4 (ν_obs = 151.44 GHz with an S/N = 8.58) and for AGS17 (ν_obs = 154.78 GHz with an S/N = 10.23). Taking into consideration the PDF of the photometric redshifts derived from the optical to near-infrared SEDs and the far-infrared SEDs as well as the infrared luminosities, we conclude that the spectroscopic redshifts of the two galaxies match the CO(6-5) transition and are = 3.556 and = 3.467. Existing spectroscopic confirmation of optically dark sources have been so far limited to a few sources, such as, z = 5.183: Walter et al. (2012); z = 3.717: Schreiber et al. (2018b); z = 3.097 & 5.113: Wang et al. (2019).

Close to 70% (4/6) of the optically dark galaxies reside in an overdensity of galaxies at z ∼ 3.5 (AGS11, 15, 17 and 24, Fig. 10). In Wang et al. (2019), the HST-dark galaxies exhibit a highly heterogeneous spatial distribution illustrated by their strong clustering measured by the cross-correlation functions. This is consistent with the present finding of a strong association of most of the blindly detected optically dark galaxies in GOODS-ALMA of four out of six galaxies located in the same proto-cluster. In addition, the redshifts of these galaxies are all consistent with this overdensity (Table 1), indicating that these galaxies are tracing a galaxy cluster in formation.

We also notice that AGS24 is not only in the barycenter of the structure at z ∼ 3.5 but also the most massive galaxy in this structure (M_⋆ = 10 M_⊙), suggesting that AGS24 is a candidate BCG in formation. In fact, AGS24 is the most massive galaxy without an AGN at z > 3 in the GOODS-ALMA field. Compared to the other five optically dark galaxies with high starburstness, AGS24 is the only MS galaxy (Table 1) indicating that it is influenced by the environment it resides⁵ in. The fact that AGS24 is a is a candidate BCG also indicates that the proto-cluster⁶ is in the process of virialization.

Acknowledgments

We thank the anonymous referee for constructive comments that helped improve the overall quality and consistency of the paper. We thank Wiphu Rujopakarn, Kristina Nyland and Preshanth Jagannathan for sharing the HUDF 3 GHz image. We thank Francesco Valentino for providing information on the CO(7-6) and [CI](2-1) emission of existing observations. We thank Ivan Delvecchio for helpful discussions on radio galaxies. L. Z. and Y. S. acknowledge the support from the National Key R&D Program of China (No. 2017YFA0404502, No. 2018YFA04027040404502) and the National Natural Science Foundation of China (NSFC grants 11825302, 11733002 and 11773013). L. Z. also acknowledges China Scholarship Council (CSC). This work was supported by the Programme National Cosmology et Galaxies (PNCG) of CNRS/INSU with INP and IN2P3, co-funded by CEA and CNES. M. F. acknowledges support from the UK Science and Technology Facilities Council (STFC) (grant number ST/R000905/1). R. D. gratefully acknowledges support from the Chilean Centro de Excelencia en Astrofísica y Tecnologías Afines (CATA) BASAL grant AFB-170002. H. I. acknowledges support from JSPS KAKENHI Grant Number JP19K23462. M. P. is supported by the ERC-StG “ClustersXCosmo”, Grant agreement 71676 VANDELS survey is based on data products created from observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programme 194.A-2003(E-T). This paper makes use of the following ALMA data: ADS/JAO.ALMA #2015.1.00543.S and #2018.1.01079.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. This paper employed Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018); APLpy, an open-source plotting package for Python (Robitaille & Bressert 2012; Robitaille 2019); Matplotlib (Hunter 2007); Numpy (Oliphant 2006); SciPy (Virtanen et al. 2020); CASA (McMullin et al. 2007); GILDAS (Gildas Team 2013).

References

All Tables

All Figures

	Fig. 1. Left: part of the spectrum obtained from the ALMA band 4 observation of AGS4. The yellow shade highlights the detected emission line. The emission line map is shown in the bottom-left corner. The contours denote 1 to 5σ, with the step of 1σ. Right: same as on the left, but for AGS17. The contours in the map denote 1 to 6σ, with the step of 1σ.
In the text

Fig. 2. Top-left: SED of AGS4. Top: photometric redshift PDF derived by EAzY (Brammer et al. 2008). Bottom: SED fittings. The green curve shows the best-fit of AGS4 at zphot = 4.02 using EAzY, and we shift the SED to the spectroscopic redshift (zspec = 3.556) of AGS4 in light blue for comparison. The gray curve shows the SED of the optically bright neighbor, IDZFOURGE = 12 333. Top-right: SED of AGS15. The blue curve shows the best-fit at zphot = 3.96 derived by EAzY and we shift the SED to the redshift peak of the overdensity, at z = 3.472, in light blue for comparison. The gray curve shows the SED of the optically bright neighbor, IDCANDELS = 3818. Bottom-left: SED of AGS24. Same as AGS15. Bottom-right: SED of AGS25 at zphot = 4.70, fitted by EAzY with the photometry retrieved from the ZFOURGE catalog.

In the text

	Fig. 3. Upper limit of the [CI](2-1)/CO(7-6) flux ratio of AGS4 (blue arrows), compared with existing observations of local and high-redshift MS galaxies, starburst galaxies and AGNs (data compiled in Valentino et al. 2018, 2020). If the detected emission line is indeed CO(7-6), then the flux ratio would fall at the lower end of the trend shown by the existing observations.
In the text

Fig. 4. Top: images of AGS11 in the J band, H band, Ks band and IRAC 3.6 μm. The magenta “x” denotes the detection in Ks-band by ZFOURGE, IDZFOURGE = 7589 (zphot = 4.82). Bottom: same as above, but for AGS17. The cyan “x”s denote the optical bright neighbors, IDCANDELS = 4414 (zphot = 0.03, below) and IDCANDELS = 4436 (zphot = 0.95, above). The magenta “x” denotes the detection in Ks-band by ZFOURGE, IDZFOURGE = 6964 (zphot = 1.85).

In the text

	Fig. 5. Images of AGS24 in the B, V, H and Ks band, IRAC 3.6 and 4.5 μm, JVLA 6 GHz (3.7σ) and 3 GHz (5.7σ). The contours are the same as in Fig. 6 but for AGS24. The three neighbor galaxies are labeled with IDZFOURGE or IDCDFS and their redshift separately.
In the text

Fig. 6. Images of AGS4 from the observed B (435 nm, rest-frame UV) to near-infrared observed Ks (2.2 μm, rest-frame B) bands. The red contours denote the ALMA detection at 1.1 mm at the resolution of 029. The green dashed contours denote the ALMA detection at 2 mm at the resolution of 088. The astrometric corrections between ALMA and HST images are applied throughout this work. The cyan “x” denotes the optical bright neighbor, IDCANDELS = 8923 (zphot = 0.24). The magenta “x” denotes the detection by ZFOURGE based on the Ks-band image, IDZFOURGE = 12 333 (zphot = 3.76). If not stated otherwise, the images in this paper are all oriented with the north up and east to the left.

In the text

	Fig. 7. Images of AGS25 in the H-band (1.6 μm), Ks-band (2.2 μm) and IRAC1-band (3.6 μm). The contribution of a 3″distant bright neighbor of AGS25 was subtracted to produce the IRAC image as in F20a. The red contours denote the ALMA detection at 1.1 mm at the resolution of 060. The magenta “x” is the ZFOURGE counterpart IDZFOURGE = 11 353.
In the text

Fig. 8. Redshift distribution of the 11 674 galaxies with a redshift (either photometric or spectroscopic) located within the 6.9′ × 10′ GOODS-ALMA field (gray). The subsample of 395 galaxies located within 50″ of the most massive galaxy AGS24, are shown in blue. The insert shows the distribution of the spectroscopic redshifts of the 121 galaxies at 3 ≤ zspec ≤ 4 over the GOODS-ALMA field. Ten galaxies fall in the redshift bin 3.470 ≤ zspec ≤ 3.480. The mean redshift of the ten galaxies is zspec = 3.472. The redshifts are from the ZFOURGE (Straatman et al. 2016) and VANDELS DR3 (McLure et al. 2018; Pentericci et al. 2018) catalogs.

In the text

	Fig. 9. Same as in Fig. 6, but for AGS15. The cyan “x” denotes the optical bright neighbor, IDCANDELS = 3818 (zphot = 4.27). The magenta “x” denotes the detection in Ks-band by ZFOURGE, IDZFOURGE = 6755 (zphot = 3.47).
In the text

Fig. 10. Projected density of 364 galaxies (40 with zspec) at 3.42 ≤ z ≤ 3.57. The contours show 1, 2, 3, 4σ level of the overdensity. Left panel: field observed in GOODS-South ALMA. Right panel: zoom-in of the extended structure in the bottom right corner of GOODS-ALMA field that encompasses 126(18) galaxies over 15 arcmin2. The red star indicates the position of AGS24. The red circle indicates a 50″radius region (2 arcmin2) including 26(3) galaxies. The yellow stars refer to the HST-galaxies AGS 11, 15 and 17 located also in this larger structure. The green circles represent the AGS 1 and 5, which are also in this redshift range but not optically dark. The white circles show the 40 galaxies with spectroscopic redshifts. We note that the two white circles falling on the red star of AGS24 are two galaxies dissociated from AGS24. The magenta square denotes the center of the proto-cluster CCPC-z34-002 identified by Franck & McGaugh (2016). The cyan square denotes Candels5001, which used to be the most massive galaxy at 3 <  z < 4 in this field (Ginolfi et al. 2017). The 50″scale in the figure is equivalent to a proper distance of 0.37 Mpc at z = 3.5.

In the text

	Fig. 11. Local overdensity (δ) of individual galaxies in the GOODS-ALMA field at 3.42 < z < 3.57. The black dots show the sliding medians of δ at different stellar masses. The filled circles indicate all the ALMA detections. The purple circles are the HST dark galaxies and the red circle indicates AGS24.
In the text

	Fig. 12. Comparison of the auto-correlation function (ACF) of the 364 galaxies at 3.42 ≤ z ≤ 3.57 in the ZFOURGE catalog (red dots), the cross-correlation function (CCF) between these galaxies and the ALMA detections (black triangles), and the CCF between these galaxies and the galaxies that are equivalently massive (M⋆ > 1010.5 M⊙) as the ALMA detections (green triangles).
In the text

	Fig. 13. Projected numbers of galaxies in the overdensity as a function of their distance to AGS24 (black circles) and the black dashed line is the linear fitting to the profile. The red line shows the best-fit projected NFW profile (Navarro et al. 1996) of a galaxy cluster at z = 2.506, CL J1001+0220 (Wang et al. 2016).
In the text