# GOODS-ALMA: Optically dark ALMA galaxies shed light on a cluster in formation at z = 3.5

L. Zhou, D. Elbaz, M. Franco, B. Magnelli, C. Schreiber, T. Wang, L. Ciesla, E. Daddi, M. Dickinson, N. Nagar, et al.

► To cite this version:

L. Zhou, D. Elbaz, M. Franco, B. Magnelli, C. Schreiber, et al. "GOODS-ALMA: Optically dark ALMA galaxies shed light on a cluster in formation at z = 3.5." Astronomy and Astrophysics - A&A, 2020, 642, pp.A155. 10.1051/0004-6361/202038059. hal-03376249

HAL Id: hal-03376249

https://hal.science/hal-03376249

Submitted on 11 Dec 2021

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# GOODS-ALMA: Optically dark ALMA galaxies shed light on a cluster in formation at z = 3.5

L. Zhou<sup>1,2,3</sup>, D. Elbaz<sup>2</sup>, M. Franco<sup>2,4</sup>, B. Magnelli<sup>5</sup>, C. Schreiber<sup>6</sup>, T. Wang<sup>2,7</sup>, L. Ciesla<sup>2,8</sup>, E. Daddi<sup>2</sup>, M. Dickinson<sup>9</sup>, N. Nagar<sup>10</sup>, G. Magdis<sup>11,12,13,14,15</sup>, D. M. Alexander<sup>16</sup>, M. Béthermin<sup>8</sup>, R. Demarco<sup>10</sup>, J. Mullaney<sup>17</sup>, F. Bournaud<sup>2</sup>, H. Ferguson<sup>18</sup>, S. L. Finkelstein<sup>19</sup>, M. Giavalisco<sup>20</sup>, H. Inami<sup>21</sup>, D. Iono<sup>22,23</sup>, S. Juneau<sup>2,9</sup>, G. Lagache<sup>8</sup>, H. Messias<sup>24,25</sup>, K. Motohara<sup>26</sup>, K. Okumura<sup>2</sup>, M. Pannella<sup>27</sup>, C. Papovich<sup>28,29</sup>, A. Pope<sup>20</sup>, W. Rujopakarn<sup>30,31,32</sup>, Y. Shi<sup>1,3</sup>, X. Shu<sup>33</sup>, and J. Silverman<sup>34</sup>

*(Affiliations can be found after the references)*

Received 31 March 2020 / Accepted 7 August 2020

## ABSTRACT

Thanks to its outstanding angular resolution, the Atacama Large Millimeter/submillimeter Array (ALMA) has recently unambiguously identified a population of optically dark galaxies with redshifts greater than $z=3$, which play an important role in the cosmic star formation in massive galaxies. In this paper we study the properties of the six optically dark galaxies detected in the 69 arcmin² GOODS-ALMA 1.1 mm continuum survey. While none of them are listed in the deepest H-band based CANDELS catalog in the GOODS-South field down to $H = 28.16$ AB, we were able to de-blend two of them from their bright neighbor and measure an H-band flux for them. We present the spectroscopic scan follow-up of five of the six sources with ALMA band 4. All are detected in the 2 mm continuum with signal-to-noise ratios higher than eight. One emission line is detected in AGS4 ($v_{obs} = 151.44$ GHz with an $S/N = 8.58$) and AGS17 ($v_{obs} = 154.78$ GHz with an $S/N = 10.23$), which we interpret in both cases as being due to the CO(6–5) line at $z_{spec}^{AGS4} = 3.556$ and $z_{spec}^{AGS17} = 3.467$, respectively. These redshifts match both the probability distribution of the photometric redshifts derived from the UV to near-infrared spectral energy distributions (SEDs) and the far-infrared SEDs for typical dust temperatures of galaxies at these redshifts. We present evidence that nearly 70% (4/6 of galaxies) of the optically dark galaxies belong to the same overdensity of galaxies at $z \sim 3.5$. overdensity The most massive one, AGS24 ($M_{\star} = 10^{11.32} M_{\odot}$), is the most massive galaxy without an active galactic nucleus at $z > 3$ in the GOODS-ALMA field. It falls in the very center of the peak of the galaxy surface density, which suggests that the surrounding overdensity is a proto-cluster in the process of virialization and that AGS24 is the candidate progenitor of the future brightest cluster galaxy.

**Key words.** galaxies: high-redshift – galaxies: evolution – galaxies: star formation – galaxies: photometry – galaxies: structure – galaxies: submillimeter: galaxies

## 1. Introduction

The discovery of bright submm-selected galaxies (SMGs; $S(850~\mu 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_{\star}$ 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 star-burst 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 (totaling 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~\mu 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~\mu m$. Thanks to the angular resolution of ALMA (0.7' 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 $\langle$SFR$\rangle = 200~M_{\odot}$ yr$^{-1}$ 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 \sim 10^{13\pm0.3} h^{-1} M_\odot$ 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 \sim 3$, which will be of interest in our analysis, this radius ranges from 5 comoving-Mpc (hereafter cMpc) for a $M(z=0) = 3 \times 10^{14} M_\odot$ bound halo to $>8$ cMpc for a Coma-like $M(z=0) = 10^{15} M_\odot$ 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 & Blazot 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 ($z \ge 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 submillimeter 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_\odot$) 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_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$, $\Omega_M = 0.3$, and $\Omega_\Lambda = 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² (6.9 × 10⁷) in the GOODS-South field (program 2015.1.00543.S, PI: D.Elbaz). GOODS-ALMA reaches an rms sensitivity of $\sigma \approx 0.18$ mJy beam⁻¹ in the 0.7'-60 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\sigma$ in GOODS-ALMA (Franco et al. 2020a, hereafter F20a). These include 19 galaxies detected above the $4.8\sigma$ limit from a blind detection approach (F18) and 16 galaxies within 3.5 and $4.8\sigma$ detected using ancillary information, mainly Spitzer-IRAC prior positions (F20a). The median redshift and stellar mass of the $S/N \ge 4.8$ sources are $z=2.73$ and $M_\star = 1.0 \times 10^{11} M_\odot$, whereas the sources with a $4.8 > S/N \ge 3.5$, are both slightly closer ($z=2.40$) and less massive ($M_\star = 7.2 \times 10^{10} M_\odot$).

#### 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 ~0.88. 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⁻<fnref n="1" />. 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*<sub>125</sub>, *J*<sub>125</sub>, *B*<sub>435</sub>, *V*<sub>606</sub>, *i*<sub>775</sub>, *i*<sub>814</sub> and *z*<sub>850</sub>. We also use images from the VLT, obtained in the *U*-band with VIMOS (Nonino et al. 2009), and in the *K*<sub>s</sub>-band with ISAAC (Retzlaff et al. 2010) and HAWK-I (Fontana et al. 2014) and the associated CANDELS multiwave-length 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*<sub>1</sub>, *J*<sub>2</sub>, *J*<sub>3</sub>) and *H* (*H*<sub>s</sub>, *H*<sub>f</sub>) bands useful for the photometric redshift determination. The ZFOURGE team produced the deepest *K*<sub>s</sub>-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 0.31' × 0.61". The exposure time is in total 177 hours reaching an rms noise at the pointing center of 0.32 µJy beam⁻<fnref n="1" />. The 3 GHz map has an angular resolution of 1.19' × 0.59" and an average rms noise of 1.03 µJy beam⁻<fnref n="1" />. At 3 GHz, AGS4, 17 and 24 are detected with an S/N of 10.4, 9.9, and 5.6 using PyBDSF<fnref n="1" /> (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<fnref n="1" />, 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*<sub>peak</sub>, 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 *σ*<sub>*z*</sub> = 0.020 at *z* > 1.5 (Straatman et al. 2016), where *σ*<sub>*z*</sub> is the normalized median absolute deviation (NMAD) of Δ*z*/(1 + *z*), that is, 1.48× the median absolute deviation (MAD) of |*z*<sub>phot</sub> − *z*<sub>spec</sub>|/(1 + *z*<sub>spec</sub>).

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.

<footnote n="1">¹ https://www.astron.nl/citt/pybdsf/</footnote>
<table><caption>**Table 1.** Properties of the optically dark galaxies.</caption><thead><tr><th>ID<br>(1)</th><th>RA<br>[deg]<br>(2)</th><th>Dec<br>[deg]<br>(3)</th><th>z<sub>phot</sub><br>(4)</th><th>z<sub>spec</sub><br>(5)</th><th>assumed z<br>(6)</th><th>log M<sub>*</sub><br>[log M<sub>☉</sub>]<br>(7)</th><th>log L<sub>IR</sub><br>[log L<sub>☉</sub>]<br>(8)</th><th>SFR<br>[M<sub>☉</sub>yr<sup>-1</sup>]<br>(9)</th><th>R<sub>SB</sub><br>(10)</th><th>F<sub>1.1mm</sub><br>[mJy]<br>(11)</th><th>F<sub>2mm</sub><br>[mJy]<br>(12)</th><th>log L<sub>CO</sub><br>log L<sub>☉</sub><br>(13)</th><th>H<br>[mag]<br>(14)</th><th>F<sub>3GHz</sub><br>[μJy]<br>(15)</th><th>F<sub>6GHz</sub><br>[μJy]<br>(16)</th><th>L<sub>1.4GHz</sub><br>[WHz<sup>-1</sup>]<br>(17)</th></tr></thead><tbody><tr><td>AGS4</td><td>53.148839</td><td>-27.821192</td><td>-</td><td>3.556</td><td>3.556</td><td>11.09<sup>+0.06</sup><br>-0.18</td><td>12.93 ± 0.02</td><td>1435<sup>+76</sup><sub>-83</sub></td><td>3.8<sup>+1.9</sup><sub>-0.5</sub></td><td>1.72 ± 0.20</td><td>0.29 ± 0.02</td><td>4.95 ± 0.93</td><td>25.23</td><td>17.28 ± 1.66</td><td>8.64 ± 0.77</td><td>2.73 × 10<sup>24</sup></td></tr><tr><td>AGS11</td><td>53.108818</td><td>-27.869055</td><td>-</td><td>3.472</td><td>3.472</td><td>10.24<sup>+0.75</sup><br>-0.01</td><td>12.94 ± 0.07</td><td>1492<sup>+52</sup><sub>-53</sub></td><td>28.7<sup>+0.1</sup><sub>-0.1</sub></td><td>1.34 ± 0.25</td><td>0.30 ± 0.02</td><td>-</td><td>-</td><td>4.48 ± 1.06</td><td>-</td><td>6.72 × 10<sup>23</sup></td></tr><tr><td>AGS15</td><td>53.074847</td><td>-27.875880</td><td>3.96<sup>+0.80</sup><sub>-0.78</sub></td><td>-</td><td>3.472</td><td>10.56<sup>+0.01</sup><br>-0.01</td><td>12.78 ± 0.08</td><td>1034<sup>+180</sup><sub>-179</sub></td><td>9.5<sup>+1.17</sup><sub>-2.6</sub></td><td>1.21 ± 0.11</td><td>0.32 ± 0.02</td><td>-</td><td>27.11</td><td>5.38 ± 1.06</td><td>-</td><td>8.08 × 10<sup>23</sup></td></tr><tr><td>AGS17</td><td>53.079374</td><td>-27.870770</td><td>-</td><td>3.467</td><td>3.467</td><td>10.52<sup>+0.40</sup><br>-0.06</td><td>13.08 ± 0.02</td><td>2070<sup>+112</sup><sub>-117</sub></td><td>20.9<sup>+3.1</sup><sub>-1.2</sub></td><td>2.30 ± 0.44</td><td>0.35 ± 0.02</td><td>5.22 ± 1.01</td><td>-</td><td>39.00 ± 3.94</td><td>1.16 ± 0.32</td><td>5.86 × 10<sup>24</sup></td></tr><tr><td>AGS24</td><td>53.087178</td><td>-27.840217</td><td>3.58<sup>+0.49</sup><sub>-0.38</sub></td><td>-</td><td>3.472</td><td>11.32<sup>+0.19</sup><br>-0.19</td><td>12.31 ± 0.11</td><td>353<sup>+80</sup><sub>-82</sub></td><td>0.6<sup>+0.1</sup><sub>-0.1</sub></td><td>0.88 ± 0.22</td><td>0.19 ± 0.02</td><td>-</td><td>-</td><td>12.43 ± 2.19</td><td>1.41 ± 0.32</td><td>1.87 × 10<sup>24</sup></td></tr><tr><td>AGS25</td><td>53.183710</td><td>-27.836515</td><td>4.64</td><td>-</td><td>4.64</td><td>10.39</td><td>12.68 ± 0.07</td><td>832<sup>+126</sup><sub>-135</sub></td><td>8.0<sup>+1.3</sup><sub>-1.3</sub></td><td>0.82 ± 0.19</td><td>-</td><td>-</td><td>3.92 ± 1.06</td><td>1.07 ± 0.32</td><td>5.90 × 10<sup>23</sup></td></tr></tbody></table>

Columns: (1) source names; (2) (3) coordinates in the ALMA image (J2000); see F18 and F20a; (4) photometric redshifts; see Sect. 2.4.1; (5) spectroscopic redshifts: the spectroscopic redshifts of AGS4 and AGS17 are derived from the CO(6-5) line measured in our ALMA spectroscopic scan follow-up; see Sect. 3; (6) assumed redshifts: We assign the redshift of AGS11, AGS15 and AGS24 to be the same as the central position of the redshift peak that presents an overdensity centered on AGS24, $z_{peak} = 3.472$, as discussed in Sects. 5.2 and 5.3, which is consistent with their photometric redshifts. For each galaxy, all the properties are measured assuming the redshift listed here: (7) stellar masses; see Sect. 2.4.2; (8) infrared luminosities; see Sect. 2.4.3; (9) SFRs; see Sect. 2.4.4; (10) starburstness; as in F20a, $R_{SB} = SFR/SFR_{MS}$, where $SFR_{MS}$ is the average SFR of MS galaxies as defined in Schreiber et al. (2018a); (11) flux densities at 1.1 mm: Peak fluxes measured using Blobcat as in F18 and the fluxes of AGS15 and AGS17 are updated in F20a; (12) peak fluxes at 2 mm from the ALMA band4 follow-up; (13) CO luminosities: the CO(6-5) line detected in our ALMA spectroscopic scan follow-up; (14) H-band AB magnitudes: After de-blending, AGS4 and AGS15 turn out to have H-band magnitudes higher than the detection limit; (15) flux densities at 3 GHz: the flux densities of AGS4, 17 and 24 are the peak flux derived using PyBDSF¹ (Rujopakarn et al., in prep.), flux densities of AGS11, 15 and 25 are the peak flux densities in the 6 GHz image, AGS11 and AGS15 do not have peak flux density higher than $3\sigma$; (17) luminosity at 1.4 GHz: derived from the flux density at 3 GHz using Eq. (1).

### 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 \ge z \ge 2.5$ and massive passive galaxies at $2.5 \ge z \ge 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 \ge z \ge 3$.

## 2.4. Derived parameters of the optically dark galaxies

### 2.4.1. Redshifts

To derive the photometric redshifts, we first de-blend<fnref n="3" /> 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} = 11353$. 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++<fnref n="4" /> 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.

<footnote n="2">² VANDELS DR3: http://archive.eso.org/cms/eso-archive-news/new-data-release-of-spectra-and-catalog-from-vandels.html; http://vandels.inaf.it/db/query_catalogs.jsp</footnote>

<footnote n="3">³ By applying the de-blending method described in Schreiber et al. (2018b), also see code developed by Corentin Schreiber: https://github.com/cschreib/qdeblend</footnote>

<footnote n="4">⁴ FAST++: https://github.com/cschreib/fastpp, a version of the code FAST Kriek et al. (2009) fully rewritten that can handle much larger parameter grids and offers additional features.</footnote>
<figure><img src="image_2.png" /><figcaption>**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σ.</figcaption></figure>

### 2.4.4. SFRs

The SFRs are measured in Franco et al. (2020b, hereafter, F20b), as SFR<sub>tot</sub> = SFR<sub>IR</sub> + SFR<sub>UV</sub>. We calculate the SFR<sub>IR</sub> based on the correlation between infrared luminosity and SFR in Kennicutt (1998) and SFR<sub>UV</sub> 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 α = α<sub>1.4 GHz</sub><sup>3 GHz</sup> = -0.8, the typical value for galaxies at z > 2 (Delhaize et al. 2017):

$$L_{1.4\,\text{GHz}} (\text{W Hz}^{-1}) = \frac{4\pi D_L^2}{(1+z)^{\alpha+1}} \left(\frac{1.4}{3}\right)^\alpha S_{3\,\text{GHz}}, \quad (1)$$

where D<sub>L</sub> is the luminosity distance to the object. The results are listed in Table 1. Then we compared the ratio L<sub>1.4 GHz</sub>/SFR with the radio excess threshold as in Delvecchio et al. (2017). All the six optically dark galaxies have L<sub>1.4 GHz</sub>/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<sub>ZFOURGE</sub> = 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 v<sub>obs</sub> = 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<sub>2</sub>O(2<sub>11</sub>-2<sub>02</sub>) lines.

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<sub>2</sub>O(2<sub>11</sub>-2<sub>02</sub>). 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 v<sub>RF</sub> = 806.7 GHz provides a spectroscopic redshift of z<sub>spec</sub><sup>AGS4</sup> = 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) (v<sub>rf</sub> = 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 v<sub>RF</sub> = 806.7 GHz gives a z<sub>spec</sub><sup>AGS4</sup> = 3.566, which agrees with the 4000 Å discontinuity in the optical SED. Therefore, we adopt the redshift to be z<sub>spec</sub><sup>AGS4</sup> = 3.556 assuming the line to be CO(6-5).

### 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
<figure><img src="image_3.png" /><figcaption>**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 $z_{phot} = 4.02$ using EAZY, and we shift the SED to the spectroscopic redshift ($z_{spec} = 3.556$) of AGS4 in light blue for comparison. The gray curve shows the SED of the optically bright neighbor, $ID_{ZFOURGE} = 12333$. *Top-right:* SED of AGS15. The blue curve shows the best-fit at $z_{phot} = 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, $ID_{CANDELS} = 3818$. *Bottom-left:* SED of AGS24. Same as AGS15. *Bottom-right:* SED of AGS25 at $z_{phot} = 4.70$, fitted by EAZY with the photometry retrieved from the ZFOURGE catalog.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**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.</figcaption></figure>

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 $\lambda < 67$ µm or at $\lambda > 134$ µm, rest-frame. 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 ($v_{rf} = 806.7$ GHz), a second line, [CI] ($v_{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 $z_{spec}^{AGS17} = 2.723$ in the case of the CO(5-4) transition and the other of $z_{spec}^{AGS17} = 3.467$ for the CO(6-5) transition. The second redshift of $z_{spec}^{AGS17} = 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 $z_{spec}^{AGS17} = 2.723$ implies a dust
<figure><img src="image_5.png" /><figcaption>**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, ID<sub>ZFOURGE</sub> = 7589 (z<sub>phot</sub> = 4.82). Bottom: same as above, but for AGS17. The cyan “x”s denote the optical bright neighbors, ID<sub>CANDELS</sub> = 4414 (z<sub>phot</sub> = 0.03, below) and ID<sub>CANDELS</sub> = 4436 (z<sub>phot</sub> = 0.95, above). The magenta “x” denotes the detection in Ks-band by ZFOURGE, ID<sub>ZFOURGE</sub> = 6964 (z<sub>phot</sub> = 1.85).</figcaption></figure>

temperature $T_{\text{dust}} \sim 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_{\text{dust}} \sim 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 $z_{\text{spec}}^{\text{AGS17}} = 3.467$ indicates a dust temperature $T_{\text{dust}} \sim 37$ K that is close to the typical dust temperature $T_{\text{dust}} \sim 40$ K at $z=3.5$ (Schreiber et al. 2018a) and is consistent with the mean dust temperature, $T_{\text{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.

### 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\sigma$. 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 $L_{\text{IR}} - L'_{\text{CO(6-5)}}$ relation with values below $3\sigma$ for the three galaxies ($2.8\sigma$, $2.3\sigma$, $0.7\sigma$ 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\sigma$ 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\mu m}(\text{AGS24}) = 4.7 \pm 0.5$ µJy and $F_{3.6\mu m}(\text{AGS25}) = 1.5 \pm 0.1$ µJy.

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\sigma$ limit that includes it was still under construction at the time. All the six optically dark galaxies are also detected above $4.5\sigma$ 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
<figure><img src="image_6.png" /><figcaption>**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 ID<sub>ZFOURGE</sub> or ID<sub>CDFS</sub> and their redshift separately.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**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 0.′29. The green dashed contours denote the ALMA detection at 2 mm at the resolution of 0.′88. The astrometric corrections between ALMA and HST images are applied throughout this work. The cyan “x” denotes the optical bright neighbor, ID<sub>CANDELS</sub> = 8923 (z<sub>phot</sub> = 0.24). The magenta “x” denotes the detection by ZFOURGE based on the Ks-band image, ID<sub>ZFOURGE</sub> = 12 333 (z<sub>phot</sub> = 3.76). If not stated otherwise, the images in this paper are all oriented with the north up and east to the left.</figcaption></figure>

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).

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 Ks-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 northwest blob that is the brightest in the near-infrared. As shown in Fig. 6, AGS4 is at 0.′5 from a bright optical emitter. They are considered as two regions of the same galaxy (ID<sub>CANDELS</sub> = 8923, z<sub>phot</sub> = 0.24, cyan “x” in Fig. 6) in the CANDELS catalog (Guo et al. 2013). The ZFOURGE position is at 0.′16 (ID<sub>ZFOURGE</sub> = 12 333) to AGS4 (magenta “x” in
Fig. 6), and the photometric redshift derived by the ZFOURGE team is $z_{\text{phot}}^{12\,333} = 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 northwestern emitter. The photometric redshift associated to this object, once deconvolved from the bright optical emitter, is $z_{\text{phot}}^{\text{AGS4}} = 4.02^{+0.36}_{-0.60}$. 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 $z_{\text{phot}}^{\text{ZFOURGE}} = 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 ($z_{\text{phot}}^{\text{CANDELS}} = 0.24$), with a clear drop of the emission above 0.8 µm, whereas the green SED exhibits a break around 1.6 µm ($z_{\text{phot}}^{\text{ZFOURGE}} = 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_{\text{spec}} = 3.566$ (see Fig. 1-left), well within the range of possibilities of our photometric redshift estimate. Adopting the redshift at $z_{\text{spec}}^{\text{AGS4}} = 3.556$, the stellar mass of AGS4 is $M_{\star} = 10^{11.09\pm0.08} M_{\odot}$, which makes AGS4 one of the most massive galaxies in the early Universe. We note that adopting $z_{\text{spec}}^{\text{AGS4}} = 4.326$, as for CO(7-6), the stellar mass that we obtain for AGS4 is $M_{\star} = 10^{11.45\pm0.20} M_{\odot}$, 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_{\text{TIR}}$-parameter (see Eq. (2)) measuring the infrared to radio flux ratio parameter is equal to $q_{\text{TIR}} = 2.49 \pm 0.06$, higher than the typical value for star-forming galaxies at $z \sim 3.5$ proposed by Delhaize et al. (2017), $q_{\text{TIR}} = 2.16 \pm 0.06$, indicating no radio excess of this galaxy.

$$q_{\text{TIR}} = \log\left(\frac{L_{\text{TIR}}}{3.75 \times 10^{12} \text{ Hz}}\right) - \log\left(\frac{L_{\text{1.4GHz}}}{\text{W Hz}^{-1}}\right), \quad (2)$$

where $L_{\text{TIR}}$ is total infrared luminosity (8–1000 µm) in unit W, and the $L_{1.4\text{GHz}} = 2.73 \times 10^{24} \text{ W Hz}^{-1}$ is derived from the $10.4\sigma$ detection of AGS4 at 3 GHz assuming a radio spectral index of $\alpha = \alpha_{1.4\text{GHz}}^{3\text{GHz}} = -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 $z_{\text{phot}}^{\text{AGS25}} = 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 $z_{\text{phot}}^{\text{AGS25}} = 4.71^{+0.24}_{-0.24}$ 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_{\star} = 10^{10.39\pm0.12} M_{\odot}$ (F20a). AGS25 is not detected by Herschel, hence its SFR = 839 ± 131 $M_{\odot}$ yr$^{-1}$ 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.

# 5. An overdensity at $z \sim 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 \sim 3.5$ in the redshift distribution

The redshift distribution of the galaxies with either a photometric or a spectroscopic redshift in the 6.9'x10' GOODS-ALMA field presents a peak at $z \sim 3.5$ (Fig. 8). This peak is at a $3.5\sigma$ significance among the 1373 galaxies at $3 \le z \le 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.
<figure><img src="image_8.png" /><figcaption>**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 0.60. The magenta "x" is the ZFOURGE counterpart ID<sub>ZFOURGE</sub> = 11353.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Redshift distribution of the 11674 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 \le z_{\text{spec}} \le 4$ over the GOODS-ALMA field. Ten galaxies fall in the redshift bin $3.470 \le z_{\text{spec}} \le 3.480$. The mean redshift of the ten galaxies is $z_{\text{spec}} = 3.472$. The redshifts are from the ZFOURGE (Straatman et al. 2016) and VANDELS DR3 (McLure et al. 2018; Pentericci et al. 2018) catalogs.</figcaption></figure>

Here we increase by a factor of three the total number of galaxies with a spectroscopic redshift that fall within the redshift range $3 \le z \le 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 \le z \le 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 \le z_{\text{spec}} \le 4$ show a clear peak at $3.470 \le z_{\text{spec}} \le 3.480$ containing ten galaxies (Fig. 8). These ten galaxies have a mean redshift at $z_{\text{spec}} = 3.472$. This is consistent with the peak at $z \sim 3.5$ shown by the 1373 galaxies with either spectroscopic or photometric redshifts.

### 5.3. optically dark galaxies at $z \sim 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

<table><caption>**Table 2.** Information of the neighbors of AGS24.</caption><thead><tr><th>ID</th><th>z<sub>spec</sub></th><th>z<sub>phot</sub></th><th>M<sub>*</sub><br>[10<sup>10</sup> × M<sub>☉</sub>]</th><th>SFR<br>[M<sub>☉</sub> yr<sup>-1</sup>]</th></tr></thead><tbody><tr><td>ZFOURGE10672</td><td>—</td><td>3.53</td><td>5.88</td><td>30.2</td></tr><tr><td>ZFOURGE10627</td><td>3.478</td><td>3.48</td><td>2.63</td><td>14.8</td></tr><tr><td>CDFS007493</td><td>3.482</td><td>—</td><td>0.41</td><td>—</td></tr></tbody></table>

peak at $z \sim 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_{\text{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 \sim 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$\sigma$ and 8.7$\sigma$ 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 0.6 (3.9$\sigma$) than in the image at the original resolution of 0.29 (4.9$\sigma$), which suggests that the dust emission is very compact. This optically dark galaxy exhibits a clear radio counterpart at 6 GHz (3.7$\sigma$) and 3 GHz (5.7$\sigma$). 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 Ks bands. The emission of AGS24 in the H and Ks 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_{\text{AGS24}} = 3.58^{+0.49}_{-0.38}$. Assuming this redshift, the excess emission in the V band may correspond to the Ly$\alpha$ 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$\alpha$ emission despite the high sensitivity of MUSE (rms = $4 \times 10^{20}$ erg s$^{-1}$ cm$^{-2}$ A$^{-1}$, 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 ~1.5 (equivalent to 11.3 kpc at $z \sim 3.5$) from AGS24 (Fig. 5 and Table 2). They all have a redshift of $z \sim 3.5$, including two spectroscopic redshifts – ID<sub>ZFOURGE</sub> = 10627, z<sub>spec</sub> = 3.478 and ID<sub>CDFS</sub> = 007493, z<sub>spec</sub> = 3.482 (VANDELS DR3) – and one photometric redshift – ID<sub>ZFOURGE</sub> = 10672, z<sub>phot</sub> = 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 \sim 3.472$ (Sect. 5.2) of the proto-cluster: (i) the three neighbors of AGS24 have $z \sim 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 \sim 3.472$ proto-cluster.

Assuming this redshift, we estimate the total infrared luminosity of AGS24 to be $L_{IR} = 10^{12.31 \pm 0.11} L_{\odot}$ 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_{\star} = 10^{11.32^{+0.02}_{-0.19}} M_{\odot}$ (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_{\odot}$.

With a redshift of $z \sim 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 (IDZFOURGE = 11505 and 16410, 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 $\chi^2$ 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 $\chi^2$ 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 \pm 0.18$ (see Eq. (2)), where the 1.4 GHz luminosity, $L_{1.4GHz} = 1.87 \times 10^{24} \text{ W Hz}^{-1}$, is converted from the JVLA 3 GHz emission ($5.7\sigma$) of AGS24 assuming a radio spectral index of $\alpha = -0.8$, as we did for AGS4. Within the rms, the infrared to radio flux ratio agrees with the value $q_{TIR} = 2.16 \pm 0.06$ for star-forming galaxies at redshift $z \sim 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

<table><caption>**Table 3.** Stellar masses of AGS24 assuming different SFHs and attenuation laws.</caption><thead><tr><th colspan="2">FAST++<sup>4</sup>:<br>Calzetti et al. (2000) attenuation law + vary SFHs</th></tr><tr><th>Models</th><th>log M<sub>★</sub> [log M<sub>☉</sub>]</th></tr></thead><tbody><tr><td>Delayed exponentially declining</td><td>11.37<sup>+0.01</sup><sub>-0.23</sub></td></tr><tr><td>Exponentialy declining</td><td>11.31<sup>+0.11</sup><sub>-0.19</sub></td></tr><tr><td>Truncated</td><td>11.40<sup>+0.00</sup><sub>-0.28</sub></td></tr><tr><td>A second burst</td><td>11.33<sup>+0.04</sup><sub>-0.18</sub></td></tr></tbody></table>

<table><thead><tr><th colspan="2">CIGALE:<br>delayed exponentially declining SFH + vary attenuation laws</th></tr><tr><th>Models</th><th>log M<sub>★</sub> [log M<sub>☉</sub>]</th></tr></thead><tbody><tr><td>Calzetti et al. (2000)</td><td>11.52 ± 0.57</td></tr><tr><td>Charlot & Fall (2000)</td><td>11.41 ± 0.21</td></tr><tr><td>MAGPHYS <sup>(a)</sup></td><td>11.44 ± 0.22</td></tr></tbody></table>

**Notes.** (a) A modified Charlot & Fall (2000) law that is used in MAGPHYS (da Cunha et al. 2008).

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 \sim 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 $0.45$ from the ALMA detection. It has been attributed with a photometric redshift of $z_{phot}^{7589} = 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 attribute the ZFOURGE ID and the associated redshift to AGS11. Instead, we consider the association of the optically dark galaxies with the overdensity at $z \sim 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.
<figure><img src="image_10.png" /><figcaption>Fig. 9. Same as in Fig. 6, but for AGS15. The cyan “x” denotes the optical bright neighbor, $ID_{CANDELS} = 3818$ ($z_{phot} = 4.27$). The magenta “x” denotes the detection in $K_s$-band by ZFOURGE, $ID_{ZFOURGE} = 6755$ ($z_{phot} = 3.47$).</figcaption></figure>

Assuming $z_{AGS11} = 3.472$, the infrared luminosity we derived for AGS11 is $L_{IR} = 10^{12.94\pm0.07} L_{\odot}$.

AGS15. AGS15 is located at 0.′59 to an optically bright galaxy ($ID_{CANDELS} = 3818$, cyan “x” in Fig. 9). It is detected at $5.22\sigma$ 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\sigma$ 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 $z_{phot}^{6755} = 3.46$ and the stellar mass is $M_{\star} = 10^{9.86} M_{\odot}$. We notice that there is a high uncertainty in the SED fitting by the ZFOURGE team, with $\chi^2 = 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++$^4$ to derive the photometric redshift and the stellar mass. The $z_{peak}$ derived for AGS15 is at $z_{phot}^{AGS15} = 3.96^{+0.80}_{-0.78}$. 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\pm0.08} L_{\odot}$, the SFR to be SFR = $1042^{+181}_{-183} M_{\odot} \text{ yr}^{-1}$ and the stellar mass to be $M_{\star} = 10^{10.56^{+0.01}_{-0.41}} M_{\odot}$.

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$, ($z_{phot}^{4414} = 0.03$, Pannella et al. 2015) at 0.′27, and $ID_{CANDELS} = 4436$, ($z_{phot}^{4436} = 0.95$, fitted with EAZY) at 0.′57 (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 $z_{phot}^{6964} = 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 $z_{spec}^{AGS17} = 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 \sim 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_{\star} = 10^{10.59^{+0.25}_{-0.19}} M_{\odot}$ 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\pm0.02} L_{\odot}$.

## 5.4. Spatial distribution of galaxies at $z \sim 3.5$ in the GOODS-ALMA field

### 5.4.1. Local overdensity ($\delta$) of individual galaxies

We calculate the overdensity ($\delta = \rho/\bar{\rho}$) of the 364 galaxies at $3.42 \le z \le 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 \sim 2.45$, identified by Cucciati et al. (2018). We derive a dispersion of $\sigma(\delta) = 0.20$ among the 364 galaxies. We do that by comparing the projected local galaxy density within 2 cMpc ($\rho$) of a galaxy, to the average number density over the GOODS-ALMA field ($\bar{\rho} = 1.32 \text{ cMpc}^{-2}$ or $5.13 \text{ arcmin}^{-2}$).

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 $z_{spec}^{AGS1} = 3.442$ of AGS1 comes from VANDELS DR3, zflag = 1, meaning 50% probability to be correct, while we note that it is recorded as $z_{spec}^{AGS1} = 2.309$ in F18. AGS5 has a photometric redshift of $z_{phot}^{AGS5} = 3.46$ (Straatman et al. 2016). As mentioned in Sect. 4.1, AGS4 has a redshift of $z_{spec}^{AGS4} = 3.556$.
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Projected density of 364 galaxies (40 with $z_{\text{spec}}$) at $3.42 \le z \le 3.57$. The contours show 1, 2, 3, $4\sigma$ 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 arcmin$^2$. The red star indicates the position of AGS24. The red circle indicates a 50"-radius region (2 arcmin$^2$) 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$.</figcaption></figure>

At the same time, eight massive galaxies fall in this redshift range in the ZFOURGE catalog. They have a stellar mass, $\log[M_*/M_{\odot}] > 10.5$, which is equivalently massive as the ALMA detections (F18). They have a median stellar mass of $M_* = 10^{10.82} M_{\odot}$. These eight massive galaxies, as well as the seven ALMA detections all show a local overdensity ($\delta > 1$). Out of the 15 massive galaxies, AGS24 shows the most prominent clustering, $\delta_{\text{AGS24}} = 2.06$, hence a $> 5\sigma$ excess with respect to the average galaxy surface density. Hence AGS24 is both the most massive galaxy at $z > 3$ in the GOODS-ALMA field with no AGN and it is located in the densest surface density peak at $z \sim 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 $\delta_{\text{AGS11}} = 1.46$ ($>2\sigma$ excess), $\delta_{\text{AGS15}} = 1.64$ ($>3\sigma$ excess), $\delta_{\text{AGS17}} = 1.76$ ($\sim 4\sigma$ excess), $\delta_{\text{AGS4}} = 1.39$ ($\sim 2\sigma$ excess), and for the other two ALMA detection with optical counterparts, it is $\delta_{\text{AGS1}} = 1.21$, $\delta_{\text{AGS5}} = 1.03$.

Regardless of mass, when focusing on the galaxies in the densest regions, that is to say, $\delta > 5\sigma(\delta)$, 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 $\text{ID}_{\text{ZFOURGE}} = 10672$, $M_* = 10^{10.76} M_{\odot}$. It is at 1"/5 from AGS24.

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

<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Local overdensity ($\delta$) of individual galaxies in the GOODS-ALMA field at $3.42 < z < 3.57$. The black dots show the sliding medians of $\delta$ 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.</figcaption></figure>

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.
### 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 0.9'×0.9' in the GOODS-ALMA field and calculate the overdensity ($\delta = \rho/\bar{\rho}$) for each cell. 364 galaxies (40 with $z_{\text{spec}}$) at $3.42 \le z \le 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$\sigma$ peak of this overdensity happens to be very close to AGS24, the galaxy with the strongest overdensity ($\delta_{\text{AGS24}} = 2.06$) this is consistent with the fact that it is in a highly clustered environment. Out of the 364/40 (redshift/$z_{\text{spec}}$) galaxies in the 69 arcmin$^2$ field, 126/18 reside in this 15 arcmin$^2$ box, and 26/3 galaxies in this redshift slice fall within 50'' (2 arcmin$^2$) 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 $\delta 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$\alpha$-emitting galaxies selected from the VANDELS survey reside in their 13 detected overdensities at $z>2$, indicating that Ly$\alpha$-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_\star = 10^{10.27} M_{\odot}$, $z_{\text{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,

<figure><img src="image_13.png" /><figcaption>Fig. 12. Comparison of the auto-correlation function (ACF) of the 364 galaxies at $3.42 \le z \le 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_\star > 10^{10.5} M_\odot$) as the ALMA detections (green triangles).</figcaption></figure>

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_{\text{spec}}$) at $3.42 \le z \le 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_\star > 10^{10.5} M_\odot$) 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 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.

### 5.5. Dynamical state of the proto-cluster at $z \sim 3.5$

As discussed in Sect. 5.3.1, AGS24 is the most massive galaxy in the field at $z>3$ ($M_\star = 10^{11.32} M_{\odot}$) and in particular at $z \sim 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 \sim 3.5$ within 1.5'' from AGS24 (11.3 kpc at $z \sim 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_\odot$ 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_\odot$ 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
<figure><img src="image_14.png" /><figcaption>**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).</figcaption></figure>

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).

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 virialization 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 0''50 and 0''27, 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 ($v_{obs} = 151.44$ GHz with an $S/N = 8.58$) and for AGS17 ($v_{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 $z_{\text{spec}}^{\text{AGS4}} = 3.556$ and $z_{\text{spec}}^{\text{AGS17}} = 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 \sim 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 \sim 3.5$ but also the most massive galaxy in this structure ($M_* = 10^{11.32^{+0.02}_{-0.19}} M_{\odot}$), 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<fnref n="5" /> in. The fact that AGS24 is a candidate BCG also indicates that the proto-cluster<fnref n="6" /> is in the process of virialization.

**Acknowledgements.** 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. 2018YFA0402704040502) 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

<footnote n="5">⁵ http://www.astropy.org</footnote>

<footnote n="6">⁶ http://www.iram.fr/IRAMFR/GILDAS</footnote>
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

Ashby, M. L. N., Willner, S. P., Fazio, G. G., et al. 2015, ApJS, 218, 33

Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33

Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123

Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57

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

Boquien, M., Burgarella, D., Roehlly, Y., et al. 2019, A&A, 622, A103

Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, ApJ, 686, 1503

Brodwin, M., Stanford, S. A., Gonzalez, A. H., et al. 2013, ApJ, 779, 138

Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000

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

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

Castignani, G., Combes, F., & Salomé, P. 2020, A&A, 635, L10

Chabrier, G. 2003, PASP, 115, 763

Charlot, S., & Fall, S. M. 2000, ApJ, 539, 718

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

Chiang, Y.-K., Overzier, R. A., Gebhardt, K., & Henriques, B. 2017, ApJ, 844, L23

Cowie, L. L., González-López, J., Barger, A. J., et al. 2018, ApJ, 865, 106

Cucciati, O., Lemaux, B. C., Zamorani, G., et al. 2018, A&A, 619, A49

da Cunha, E., Charlot, S., & Elbaz, D. 2008, MNRAS, 388, 1595

da Cunha, E., Walter, F., Smail, I. R., et al. 2015, ApJ, 806, 110

Daddi, E., Dickinson, M., Morrison, G., et al. 2007, ApJ, 670, 156

Daddi, E., Dannerbauer, H., Stern, D., et al. 2009, ApJ, 694, 1517

Davidzon, I., Ilbert, O., Laigle, C., et al. 2017, A&A, 605, A70

Delhaize, J., Smolčić, V., Delvecchio, I., et al. 2017, A&A, 602, A4

De Lucia, G., & Blaizot, J. 2007, MNRAS, 375, 2

Delvecchio, I., Smolčić, V., Zamorani, G., et al. 2017, A&A, 602, A3

Draine, B. T., Dale, D. A., Bendo, G., et al. 2007, ApJ, 663, 866

Dunlop, J. S., McLure, R. J., Biggs, A. D., et al. 2017, MNRAS, 466, 861

Elbaz, D., Daddi, E., Le Borgne, D., et al. 2007, A&A, 468, 33

Elbaz, D., Dickinson, M., Hwang, H. S., et al. 2011, A&A, 533, A119

Elbaz, D., Leiton, R., Nagar, N., et al. 2018, A&A, 616, A110

Fontana, A., Dunlop, J. S., Paris, D., et al. 2014, A&A, 570, A11

Forrest, B., Tran, K.-V. H., Broussard, A., et al. 2017, ApJ, 838, L12

Forrest, B., Annunziatella, M., Wilson, G., et al. 2020, ApJ, 890, L1

Franck, J. R., & McGaugh, S. S. 2016, ApJ, 817, 158

Franco, M., Elbaz, D., Béthermin, M., et al. 2018, A&A, 620, A152

Franco, M., Elbaz, D., Zhou, L., et al. 2020a, A&A, in press, https://doi.org/10.1051/0004-6361/202038310

Franco, M., Elbaz, D., Zhou, L., et al. 2020b, A&A, in press, https://doi.org/10.1051/0004-6361/202038312

Giacconi, R., Zirm, A., Wang, J., et al. 2002, ApJS, 139, 369

Gildas Team, 2013, GILDAS: Greoble Image and Line Data Analysis Software

Ginolfi, M., Maiolino, R., Nagao, T., et al. 2017, MNRAS, 468, 3468

Gobat, R., Rosati, P., Strazzullo, V., et al. 2008, A&A, 488, 853

Gómez-Guijarro, C., Riechers, D. A., Pavesi, R., et al. 2019, ApJ, 872, 117

Grogin, N. A., Kocevski, D. D., Faber, S. M., et al. 2011, ApJS, 197, 35

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

Guo, Y., Ferguson, H. C., Giavalisco, M., et al. 2013, ApJS, 207, 24

Hatsukade, B., Kohno, K., Yamaguchi, Y., et al. 2018, PASJ, 70, 105

Herenz, E. C., Urrutia, T., Wisotzki, L., et al. 2017, A&A, 606, A12

Hodge, J. A., Smail, I., Walter, F., et al. 2019, ApJ, 876, 130

Hughes, D. H., Serjeant, S., Dunlop, J., et al. 1998, Nature, 394, 241

Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90

Jin, S., Daddi, E., Magdis, G. E., et al. 2019, ApJ, 887, 144

Kennicutt, R. C. Jr. 1998, ApJ, 498, 541

Kravtsov, A. V., & Borgani, S. 2012, ARA&A, 50, 353

Kriek, M., van Dokkum, P. G., Labbé, I., et al. 2009, ApJ, 700, 221

Labbé, I., Oesch, P. A., Illingworth, G. D., et al. 2015, ApJS, 221, 23

Landy, S. D., & Szalay, A. S. 1993, ApJ, 412, 64

Liu, D., Gao, Y., Isaak, K., et al. 2015, ApJ, 810, L14

Luo, B., Brandt, W. N., Xue, Y. Q., et al. 2017, ApJS, 228, 2

McLure, R. J., Pentericci, L., Cimatti, A., et al. 2018, MNRAS, 479, 25

McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in Astronomical Data Analysis Software and Systems XVI, eds. R. A. Shaw, F. Hill, & D. J. Bell, ASP Conf. Ser., 376, 127

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

Muldrew, S. I., Hatch, N. A., & Cooke, E. A. 2015, MNRAS, 452, 2528

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563

Newman, A. B., Ellis, R. S., Andreon, S., et al. 2014, ApJ, 788, 51

Noeske, K. G., Weiner, B. J., Faber, S. M., et al. 2007, ApJ, 660, L43

Nonino, M., Dickinson, M., Rosati, P., et al. 2009, ApJS, 183, 244

Oliphant, T. E. 2006, A Guide to NumPy, 1 (USA: Trelgol Publishing)

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

Pannella, M., Elbaz, D., Daddi, E., et al. 2015, ApJ, 807, 141

Papovich, C., Bassett, R., Lotz, J. M., et al. 2012, ApJ, 750, 93

Pentericci, L., McLure, R. J., Garilli, B., et al. 2018, A&A, 616, A174

Pillepich, A., Nelson, D., Hernquist, L., et al. 2018, MNRAS, 475, 648

Retzlaff, J., Rosati, P., Dickinson, M., et al. 2010, A&A, 511, A50

Robitaille, T. 2019, APLpy v2.0: The Astronomical Plotting Library in Python

Robitaille, T., & Bressert, E. 2012, APLpy: Astronomical Plotting Library in Python. (Astrophysics Source Code Library)

Rodighiero, G., Daddi, E., Baronchelli, I., et al. 2011, ApJ, 739, L40

Rosdahl, J., & Blaizot, J. 2012, MNRAS, 423, 344

Rujopakarn, W., Dunlop, J. S., Rieke, G. H., et al. 2016, ApJ, 833, 12

Rujopakarn, W., Daddi, E., Rieke, G. H., et al. 2019, ApJ, 882, 107

Salpeter, E. E. 1955, ApJ, 121, 161

Schreiber, C., Pannella, M., Elbaz, D., et al. 2015, A&A, 575, A74

Schreiber, C., Elbaz, D., Pannella, M., et al. 2018a, A&A, 609, A30

Schreiber, C., Labbé, I., Glazebrook, K., et al. 2018b, A&A, 611, A22

Springel, V., White, S. D. M., Jenkins, A., et al. 2005, Nature, 435, 629

Straatman, C. M. S., Spitler, L. R., Quadri, R. F., et al. 2016, ApJ, 830, 51

Swinbank, A. M., Simpson, J. M., Smail, I., et al. 2014, MNRAS, 438, 1267

Thomas, D., Maraston, C., Bender, R., & Mendes de Oliveira, C. 2005, ApJ, 621, 673

Toshikawa, J., Malkan, M. A., Kashikawa, N., et al. 2020, ApJ, 888, 89

Valentino, F., Magdis, G. E., Daddi, E., et al. 2018, ApJ, 869, 27

Valentino, F., Magdis, G. E., Daddi, E., et al. 2020, ApJ, 890, 24

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261

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

Wang, T., Elbaz, D., Daddi, E., et al. 2016, ApJ, 828, 56

Wang, T., Schreiber, C., Elbaz, D., et al. 2019, Nature, 572, 211

Webb, T., Noble, A., DeGroot, A., et al. 2015, ApJ, 809, 173

Willis, J. P., Canning, R. E. A., Noordeh, E. S., et al. 2020, Nature, 577, 39

Yang, C., Omont, A., Beelen, A., et al. 2016, A&A, 595, A80

¹ School of Astronomy and Space Science, Nanjing University, Nanjing 210093, PR China

² AIM, CEA, CNRS, Université Paris-Saclay, Université Paris Diderot, Sorbonne Paris Cité, 91191 Gif-sur-Yvette, France e-mail: *luwenjia.zhou@cea.fr*

³ Key Laboratory of Modern Astronomy and Astrophysics (Nanjing University), Ministry of Education, Nanjing 210093, PR China

⁴ Centre for Astrophysics Research, University of Hertfordshire, Hatfield AL10 9AB, UK

⁵ Argelander-Institut für Astronomie, Universität Bonn, Auf dem Hügel 71, 53121 Bonn, Germany

⁶ Department of Physics, University of Oxford, Keble Road, Oxford OX1 3RH, UK

⁷ Institute of Astronomy, University of Tokyo, 2-21-1 Osawa, Mitaka, Tokyo 181-0015, Japan

⁸ Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France

⁹ Community Science and Data Center/NSF's NOIRLab, 950 N. Cherry Ave., Tucson, AZ 85719, USA

¹⁰ Departamento de Astronomía, Facultad de Ciencias Físicas y Matemáticas, Universidad de Concepción, Concepción, Chile

¹¹ Cosmic Dawn Center (DAWN), Copenhagen, Denmark

¹² Niels Bohr Institute, University of Copenhagen, Lyngbyvej 2, 2100 Copenhagen Ø, Denmark

¹³ DTU Space, National Space Institute, Technical University of Denmark, Elektrovej 327, 2800 Kgs. Lyngby, Denmark

¹⁴ Niels Bohr Institute, University of Copenhagen, 2100 Copenhagen, Denmark

¹⁵ Institute for Astronomy, Astrophysics, Space Applications and Remote Sensing, National Observatory of Athens, 15236 Athens, Greece

¹⁶ Centre for Extragalactic Astronomy, Department of Physics, Durham University, Durham DH1 3LE, UK

¹⁷ Department of Physics and Astronomy, The University of Sheffield, Hounsfield Road, Sheffield S3 7RH, UK
18 Space Telescope Science Institute, 3700 San Martin Drive, Balti-
more, MD 21218, USA

19 Department of Astronomy, The University of Texas at Austin, Austin, TX 78712, USA

20 Astronomy Department, University of Massachusetts, Amherst, MA 01003, USA

21 Hiroshima Astrophysical Science Center, Hiroshima University, 1-3-1 Kagamiyama, Higashi-Hiroshima, Hiroshima 739-8526, Japan

22 National Astronomical Observatory of Japan, National Institutes of Natural Sciences, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan

23 SOKENDAI (The Graduate University for Advanced Studies), 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan

24 Universidad de Concepción, Barrio Universitario, Concepción, Chile

25 Instituto de Astrofísica e Ciências Espaciais, Observatório Astronómico de Lisboa, Tapada da Ajuda 1349-018 Lisbon, Portugal

26 Institute of Astronomy, Graduate School of Science, The University of Tokyo, 2-21-1 Osawa, Mitaka, Tokyo 181-0015, Japan

27 Astronomy Unit, Department of Physics, University of Trieste, via Tiepolo 11, 34131 Trieste, Italy

28 Department of Physics and Astronomy, Texas A&M University, College Station, TX 77843-4242, USA

29 George P. and Cynthia Woods Mitchell Institute for Fundamental Physics and Astronomy, Texas A&M University, College Station, TX 77843-4242, USA

30 Department of Physics, Faculty of Science, Chulalongkorn University, 254 Phayathai Road, Pathumwan, Bangkok 10330, Thailand

31 National Astronomical Research Institute of Thailand (Pub-
lic Organization), Donkaew, Maerim, Chiangmai 50180,
Thailand

32 Kavli Institute for the Physics and Mathematics of the Universe (WPI), The University of Tokyo Institutes for Advanced Study, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan

33 Department of Physics, Anhui Normal University, Wuhu, Anhui 241000, PR China

34 Kavli Institute for the Physics and Mathematics of the Universe (WPI), The University of Tokyo Institutes for Advanced Study, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan