© ESO 2018

1. Introduction

In the late 1990s a population of galaxies was discovered at submillimeter wavelengths using the Submillimeter Common-User Bolometer Array (SCUBA; Holland et al. 1999) on the James Clerk Maxwell Telescope (see e.g., Smail et al. 1997; Hughes et al. 1998; Barger et al. 1998; Blain et al. 2002). These “submillimeter galaxies” or SMGs are highly obscured by dust, typically located around z ∼ 2–2.5 (e.g., Chapman et al. 2003; Wardlow et al. 2011; Yun et al. 2012), massive (M_⋆ >  7 × 10¹⁰ M_⊙; e.g., Chapman et al. 2005; Hainline et al. 2011; Simpson et al. 2014), gas-rich (f_gas >  50%; e.g., Daddi et al. 2010), with huge star formation rates (SFR) – often greater than 100 M_⊙ yr⁻¹ (e.g., Magnelli et al. 2012; Swinbank et al. 2014) – making them significant contributors to the cosmic star formation (e.g., Casey et al. 2013), often driven by mergers (e.g., Tacconi et al. 2008; Narayanan et al. 2010) and often host an active galactic nucleus (AGN; e.g., Alexander et al. 2008; Pope et al. 2008; Wang et al. 2013). These SMGs are plausible progenitors of present-day massive early-type galaxies (e.g., Cimatti et al. 2008; Michałowski et al. 2010).

Recently, thanks to the advent of the Atacama Large Millimeter/submillimeter Array (ALMA) and its capabilities to perform both high-resolution and high-sensitivity observations, our view of SMGs has become increasingly refined. The high angular resolution compared to single-dish observations reduces drastically the uncertainties of source confusion and blending, and affords new opportunities for robust galaxy identification and flux measurement. The ALMA sensitivity allows for the detection of sources down to 0.1 mJy (e.g., Carniani et al. 2015), the analysis of populations of dust-poor high-z galaxies (Fujimoto et al. 2016) or main sequence (MS; Noeske et al. 2007; Rodighiero et al. 2011; Elbaz et al. 2011) galaxies (e.g., Papovich et al. 2016; Dunlop et al. 2017; Schreiber et al. 2017), and also demonstrates that the extragalactic background light (EBL) can be resolved partially or totally by faint galaxies (S <  1 mJy; e.g., Hatsukade et al. 2013; Ono et al. 2014; Carniani et al. 2015; Fujimoto et al. 2016). Thanks to this new domain of sensitivity, ALMA is able to unveil less extreme objects, bridging the gap between massive starbursts and more normal galaxies: SMGs no longer stand apart from the general galaxy population.

However, many previous ALMA studies have been based on biased samples, with prior selection (pointing) or a posteriori selection (e.g., based on HST detections) of galaxies, or in a relatively limited region. In this study, we present an unbiased view of a large (69 arcmin²) region of the sky, without prior or a posteriori selection based on already known galaxies, in order to improve our understanding of dust-obscured star formation and investigate the main properties of these objects. We take advantage of one of the most uncertain and potentially transformational outputs of ALMA – its ability to reveal a new class of galaxies through serendipitous detections. This is one of the main reasons for performing blind extragalactic surveys.

Thanks to the availability of very deep, panchromatic photometry at rest-frame UV, optical and NIR in legacy fields such as great observatories origins deep survey-South (GOODS-South), which also includes among the deepest available X-ray and radio maps, precise multiwavelength analysis that include the crucial FIR region is now possible with ALMA. In particular, a population of high redshift (2 <  z <  4) galaxies, too faint to be detected in the deepest HST-WFC3 images of the GOODS-South field has been revealed, thanks to the thermal dust emission seen by ALMA. Sources without an HST counterpart in the H-band, the reddest available (so-called HST-dark) have been previously found by color selection (e.g., Huang et al. 2011; Caputi et al. 2012, 2015; Wang et al. 2016), by serendipitous detection of line emitters (e.g., Ono et al. 2014) or in the continuum (e.g., Fujimoto et al. 2016). We will show that ∼20% of the sources detected in the survey described in this paper are HST-dark, and strong evidence suggests that they are not spurious detections.

The aim of the work presented in this paper is to exploit a 69 arcmin² ALMA image reaching a sensitivity of 0.18 mJy at a resolution of 0″60. We used the leverage of the excellent multiwavelength supporting data in the GOODS-South field: the cosmic assembly near-infrared deep extragalactic legacy survey (Koekemoer et al. 2011; Grogin et al. 2011), the Spitzer extended deep survey (Ashby et al. 2013), the GOODS-Herschel survey (Elbaz et al. 2011), the Chandra deep field-South (Luo et al. 2017) and ultra-deep radio imaging with the VLA (Rujopakarn et al. 2016), to construct a robust catalog and derive physical properties of ALMA-detected galaxies. The region covered by ALMA in this survey corresponds to the region with the deepest HST-WFC3 coverage, and has also been chosen for a guaranteed time observation (GTO) program with the James Webb Space Telescope (JWST).

This paper is organized as follows: in Sect. 2 we describe our ALMA survey, the data reduction, and the multiwavelength ancillary data which support our studies. In Sect. 3, we present the methodology and criteria used to detect sources, we also present the procedures used to compute the completeness and the fidelity of our flux measurements. In Sect. 4 we detail the different steps we conducted to construct a catalog of our detections. In Sect. 5 we estimate the differential and cumulative number counts from our detections. We compare these counts with other (sub)millimeter studies. In Sect. 6 we investigate some properties of our galaxies such as redshift and mass distributions. Other properties will be analyzed in Franco et al. (in prep.) and finally in Sect. 8, we summarize the main results of this study. Throughout this paper, we adopt a spatially flat ΛCDM cosmological model with H₀ = 70 kms⁻¹ Mpc⁻¹, Ω_m = 0.7 and Ω_Λ = 0.3. We assume a Salpeter (Salpeter 1955) initial mass function (IMF). We used the conversion factor of M_⋆ (Salpeter 1955, IMF) = 1.7 × M_⋆ (Chabrier 2003, IMF). All magnitudes are quoted in the AB system (Oke & Gunn 1983).

2. ALMA GOODS-South survey data

2.1. Survey description

Our ALMA coverage extends over an effective area of 69 arcmin² within the GOODS-South field (Fig. 1), centered at α = 3^h32^m30.0^s, δ = −27° 48′00″ (J2000; 2015.1.00543.S; PI: D. Elbaz). To cover this ∼10′×7′ region (comoving scale of 15.1 Mpc × 10.5 Mpc at z = 2), we designed a 846-pointing mosaic, each pointing being separated by 0.8 times the antenna half power beam width (HPBW ∼ 23″3).

Fig. 1. ALMA 1.1 mm image tapered at 0″60. The white circles have a diameter of 4 arcseconds and indicate the positions of the galaxies listed in Table 3. Black contours show the different slices (labeled A–F) used to compose the homogeneous 1.1 mm coverage, with a median rms-noise of 0.18 mJy per beam. Blue lines show the limits of the HST/ACS field and green lines indicate the HST-WFC3 deep region. The cyan contour represents the limit of the Dunlop et al. (2017) survey covering all the Hubble Ultra Deep Field region. All of the ALMA-survey field is encompassed by the Chandra deep field-South.

To accommodate such a large number of pointings within the ALMA Cycle 3 observing mode restrictions, we divided this mosaic into six parallel, slightly overlapping, submosaics of 141 pointing each. To get a homogeneous pattern over the 846 pointings, we computed the offsets between the submosaics so that they connected with each other without breaking the hexagonal pattern of the ALMA mosaics.

Each submosaic (or slice) had a length of 6.8 arcmin, a width of 1.5 arcmin and an inclination (PA) of 70 deg (see Fig. 1). This required three execution blocks (EBs), yielding a total on-source integration time of ∼60 s per pointing (Table 1). We determined that the highest frequencies of the band 6 were the optimal setup for a continuum survey and we thus set the ALMA correlator to Time Division Multiplexing (TDM) mode and optimized the setup for continuum detection at 264.9 GHz (λ = 1.13 mm) using four 1875 MHz-wide spectral windows centered at 255.9 GHz, 257.9 GHz, 271.9 GHz and 273.9 GHz, covering a total bandwidth of 7.5 GHz. The TDM mode has 128 channels per spectral window, providing us with ∼37 km s⁻¹ velocity channels.

Observations were taken between the 1st of August and the 2nd of September 2016, using ∼40 antennae (see Table 1) in configuration C40-5 with a maximum baseline of ∼1500 m. J0334–4008 and J0348–2749 (VLBA calibrator and hence has a highly precise position) were systematically used as flux and phase calibrators, respectively. In 14 EBs, J0522–3627 was used as bandpass calibrator, while in the remaining 4 EBs J0238+1636 was used. Observations were taken under nominal weather conditions with a typical precipitable water vapor of ∼1 mm.

2.2. Data reduction

All EBs were calibrated with CASA (McMullin et al. 2007) using the scripts provided by the ALMA project. Calibrated visibilities were systematically inspected and few additional flaggings were added to the original calibration scripts. Flux calibrations were validated by verifying the accuracy of our phase and bandpass calibrator flux density estimations. Finally, to reduce computational time for the forthcoming continuum imaging, we time- and frequency-averaged our calibrated EBs over 120 s and 8 channels, respectively.

Imaging was done in CASA using the multifrequency synthesis algorithm implemented within the task CLEAN. Submosaics were produced separately and combined subsequently using a weighted mean based on their noise maps. As each submosaic was observed at different epochs and under different weather conditions, they exhibit different synthesized beams and sensitivities (Table 1). Submosaics were produced and primary beam corrected separately, to finally be combined using a weighted mean based on their noise maps. To obtain a relatively homogeneous and circular synthesized beam across our final mosaic, we applied different u, v tapers to each submosaic. The best balance between spatial resolution and sensitivity was found with a homogeneous and circular synthesized beam of 0″29 full width half maximum (FWHM; hereafter 0″29-mosaic; Table 1). This resolution corresponds to the highest resolution for which a circular beam can be synthesized for the full mosaic. We also applied this tapering method to create a second mosaic with a homogeneous and circular synthesized beam of 0″60 FWHM (hereafter 0″60-mosaic; Table 1), in other words, optimized for the detection of extended sources. Mosaics with even coarser spatial resolution could not be created because of drastic sensitivity and synthesized beam shape degradations.

Due to the good coverage in the uv-plane (see Fig. 2) and the absence of very bright sources (the sources present in our image do not cover a large dynamic range in flux densities; see Sect. 4), we decided to work with the dirty map. This prevents introducing potential biases during the CLEAN process and we noticed that the noise in the clean map is not significantly different (< 1%).

Fig. 2. uv-coverage of one of the 846 ALMA pointings constituting this survey. This uv-coverage allows us to perform the source detection in the dirty map.

2.3. Building of the noise map

We built the rms-map of the ALMA survey by a k-σ clipping method. In steps of four pixels on the image map, the standard deviation was computed in a square of 100 × 100 pixels around the central pixel. The pixels, inside this box, with values greater than three times the standard deviation (σ) from the median value were masked. This procedure was repeated three times. Finally, we assigned the value of the standard deviation of the non-masked pixels to the central pixel. This box size corresponds to the smallest size for which the value of the median pixel of the rms map converges to the typical value of the noise in the ALMA map while taking into account the local variation of noise. The step of four pixels corresponds to a subsampling of the beam so, the noise should not vary significantly on this scale. The median value of the standard deviation is 0.176 mJy beam⁻¹. In comparison, the Gaussian fit of the unclipped map gave a standard deviation of 0.182 mJy beam⁻¹. We adopted a general value of rms sensitivity σ = 0.18 mJy beam⁻¹. The average values for the 0″29-mosaic and the untapered mosaic are given in Table 1.

2.4. Ancillary data

The area covered by this survey is ideally located, in that it profits from ancillary data from some of the deepest sky surveys at infrared (IR), optical and X-ray wavelengths. In this section, we describe all of the data that were used in the analysis of the ALMA detected sources in this paper.

2.4.1. Optical and near-infrared imaging

We have supporting data from the Cosmic Assembly Near-IR Deep Extragalactic Legacy Survey (CANDELS; Grogin et al. 2011) with images obtained with the wide field camera 3/Infrared Channel (WFC3/IR) and UVIS channel, along with the advanced camera for surveys (ACS; Koekemoer et al. 2011. The area covered by this survey lies in the deep region of the CANDELS program (central one-third of the field). The 5-σ detection depth for a point-source reaches a magnitude of 28.16 for the H₁₆₀ filter (measured within a fixed aperture of 0.17″ Guo et al. 2013). The CANDELS/Deep program also provides images in seven other bands: the Y₁₂₅, J₁₂₅, B₄₃₅, V₆₀₆, i₇₇₅, i₈₁₄ and z₈₅₀ filters, reaching 5-σ detection depths of 28.45, 28.35, 28.95, 29.35, 28.55, 28.84, and 28.77 mag respectively. The Guo et al. (2013) catalog also includes galaxy magnitudes from the VLT, taken in the U-band with VIMOS (Nonino et al. 2009), and in the K_s-band with ISAAC (Retzlaff et al. 2010) and HAWK-I (Fontana et al. 2014).

In addition, we used data coming from the FourStar galaxy evolution survey (ZFOURGE, PI: I. Labbé) on the 6.5 m Magellan Baade telescope. The FourStar instrument (Persson et al. 2013) observed the CDFS (encompassing the GOODS-South field) through five near-IR medium-bandwidth filters (J₁, J₂, J₃, H_s, H_l) as well as broad-band K_s. By combination of the FourStar observations in the K_s-band and previous deep and ultra-deep surveys in the K-band, VLT/ISAAC/K (v2.0) from GOODS (Retzlaff et al. 2010), VLT/HAWK-I/K from HUGS (Fontana et al. 2014), CFHST/WIRCAM/K from TENIS (Hsieh et al. 2012) and Magellan/PANIC/K in HUDF (PI: I. Labbé), a super-deep detection image has been produced. The ZFOURGE catalog reaches a completeness greater than 80% to K_s <  25.3–25.9 (Straatman et al. 2016).

We used the stellar masses and redshifts from the ZFOURGE catalog, except when spectroscopic redshifts were available. Stellar masses have been derived from Bruzual & Charlot (2003) models (Straatman et al. 2016) assuming exponentially declining star formation histories and a dust attenuation law as described by Calzetti et al. (2000).

2.4.2. Mid/far-infrared imaging

Data in the mid and far-IR are provided by the infrared array camera (IRAC; Fazio et al. 2004) at 3.6, 4.5, 5.8, and 8 μm, Spitzer multiband imaging photometer (MIPS; Rieke et al. 2004) at 24 μm, Herschel photodetector array camera and spectrometer (PACS; Poglitsch et al. 2010) at 70, 100 and 160 μm, and Herschel spectral and photometric imaging receiver (SPIRE; Griffin et al. 2010) at 250, 350, and 500 μm.

The IRAC observations in the GOODS-South field were taken in February 2004 and August 2004 by the GOODS Spitzer legacy project (PI: M. Dickinson). These data have been supplemented by the Spitzer extended deep survey (SEDS; PI: G. Fazio) at 3.6 and 4.5 μm (Ashby et al. 2013) as well as the Spitzer-cosmic assembly near-infrared deep extragalactic survey (S-CANDELS; Ashby et al. 2015) and recently by the ultradeep IRAC imaging at 3.6 and 4.5 μm (Labbé et al. 2015).

The flux extraction and deblending in 24 μm imaging have been provided by Magnelli et al. (2009) to reach a depth of S₂₄ ∼ 30 μJy. Herschel images come from a 206.3 h GOODS-South observational program (Elbaz et al. 2011) and combined by Magnelli et al. (2013) with the PACS evolutionary probe (PEP) observations (Lutz et al. 2011). Because the SPIRE confusion limit is very high, we used the catalog of Wang et al. (in prep.), which is built with a state-of-the-art de-blending method using optimal prior sources positions from 24 μm and Herschel PACS detections.

2.4.3. Complementary ALMA data

As the GOODS-South field encompasses the Hubble Ultra Deep Field (HUDF), we took advantage of deep 1.3-mm ALMA data of the HUDF. The ALMA image of the full HUDF reaches a σ_1.3 mm = 35 μJy (Dunlop et al. 2017), over an area of 4.5 arcmin² that was observed using a 45-pointing mosaic at a tapered resolution of 0.7″. These observations were taken in two separate periods from July to September 2014. In this region, 16 galaxies were detected by Dunlop et al. (2017), three of them with a high S/N (S/N > 14), the other 13 with lower S/Ns (3.51 <  S/N <  6.63).

2.4.4. Radio imaging

We also used radio imaging at 5 cm from the Karl G. Jansky Very Large Array (VLA). These data were observed during 2014 March–2015 September for a total of 177 h in the A, B, and C configurations (PI: W. Rujopakarn). The images have a 0″31  ×  0″61 synthesized beam and an rms noise at the pointing center of 0.32 μJy beam⁻¹ (Rujopakarn et al. 2016). Here, 179 galaxies were detected with a significance greater than 3σ over an area of 61 arcmin² around the HUDF field, with a rms sensitivity better than 1 μJy beam⁻¹. However, this radio survey does not cover the entire ALMA area presented in this paper.

2.4.5. X-ray

The Chandra deep field-South (CDF-S) was observed for 7 Msec between 2014 June and 2016 March. These observations cover a total area of 484.2 arcmin², offset by just 32″ from the center of our survey, in three X-ray bands: 0.5–7.0 keV, 0.5–2.0 keV, and 2–7 keV (Luo et al. 2017). The average flux limits over the central region are 1.9 × 10⁻¹⁷, 6.4 × 10⁻¹⁸, and 2.7 × 10⁻¹⁷ erg cm⁻² s⁻¹ respectively. This survey enhances the previous X-ray catalogs in this field, the 4 Msec Chandra exposure (Xue et al. 2011) and the 3 Msec XMM-Newton exposure (Ranalli et al. 2013). We will use this X-ray catalog to identify candidate X-ray active galactic nuclei (AGN) among our ALMA detections.

3. Source detection

The search for faint sources in high-resolution images with moderate source densities faces a major limitation. At the native resolution (0″25  ×  0″23), the untapered ALMA mosaic encompasses almost four million independent beams, where the beam area is A_beam = π × FWHM²/(4ln(2)). It results that a search for sources above a detection threshold of 4-σ would include as many as 130 spurious sources assuming a Gaussian statistics. Identifying the real sources from such catalog is not possible. In order to increase the detection quality to a level that ensures a purity greater than 80% – in other words, the excess of sources in the original mosaic needs to be five times greater than the number of detections in the mosaic multiplied by (−1) – we have decided to use a tapered image and adapt the detection threshold accordingly.

By reducing the weight of the signal originating from the most peripheral ALMA antennae, the tapering reduces the angular resolution hence the number of independent beams at the expense of collected light. The lower angular resolution presents the advantage of optimizing the sensitivity to point sources – we recall that 0″24 corresponds to a proper size of only 2 kpc at z ∼ 1–3 – and therefore will result in an enhancement of the signal-to-noise ratio (S/N) for the sources larger than the resolution.

We chose to taper the image with a homogeneous and circular synthesized beam of 0″60 FWHM – corresponding to a proper size of 5 kpc at z ∼ 1–3 – having tested various kernels and finding that this beam was optimized for our mosaic, avoiding both a beam degradation and a too heavy loss of sensitivity. This tapering reduces by nearly an order of magnitude the number of spurious sources expected at a 4-σ level down to about 19 out of 600 000 independent beams. However, we will check in a second step whether we may have missed in the process some compact sources by also analyzing the 0″29 tapered map.

We also excluded the edges of the mosaic, where the standard deviation is larger than 0.30 mJy beam⁻¹ in the 0″60-mosaic. The effective area was thus reduced by 4.9% as compared to the full mosaic (69.46 arcmin² out of 72.83 arcmin²).

To identify the galaxies present on the image, we used BLOBCAT (Hales et al. 2012). BLOBCAT is a source extraction software using a “flood fill” algorithm to detect and catalog blobs (see Hales et al. 2012). A blob is defined by two criteria:

• at least one pixel has to be above a threshold (σ_p)

• all the adjacent surrounding pixels must be above a floodclip threshold (σ_f)

where σ_p and σ_f are defined in number of σ, the local rms of the mosaic.

A first guess to determine the detection threshold σ_p is provided by the examination of the pixel distribution of the S/N-map. The S/N-map has been created by dividing the 0″60 tapered map by the noise map. Figure 3 shows that the S/N-map follows an almost perfect Gaussian below S/N = 4.2. Above this threshold, a significant difference can be observed that is characteristic of the excess of positive signal expected in the presence of real sources in the image. However, this histogram alone cannot be used to estimate a number of sources because the pixels inside one beam are not independent of one another. Hence although the non-Gaussian behavior appears around S/N = 4.2 we performed simulations to determine the optimal values of σ_p and σ_f.

Fig. 3. Histogram of pixels of the S/N map, where pixels with noise >0.3 mJy beam−1 have been removed. The red dashed line is the best Gaussian fit. The green dashed line is indicative and shows where the pixel brightness distribution moves away from the Gaussian fit. This is also the 4.2σ level corresponding to a peak flux of 0.76 mJy for a typical noise per beam of 0.18 mJy. The solid black line corresponds to our peak threshold of 4.8σ (0.86 mJy).

We first conducted positive and negative – on the continuum map multiplied by (−1) – detection analysis for a range of σ_p and σ_f values ranging from σ_p = 4 to 6 and σ_f = 2.5 to 4 with intervals of 0.05 and imposing each time σ_p ≥ σ_f. The difference between positive and negative detections for each pair of (σ_p, σ_f) values provides the expected number of real sources.

We then searched for the pair of threshold parameters to find the best compromise between (i) providing the maximum number of detections, and (ii) minimizing the number of spurious sources. The later purity criterion, p_c, is defined as:

(1)

where N_p and N_n are the numbers of positive and negative detections respectively. To ensure a purity of 80% as discussed above, we enforced p_c ≥ 0.8. This led to σ_p = 4.8σ when fixing the value of σ_f = 2.7σ (see Fig. 4-left). Below σ_p = 4.8σ, the purity criterion rapidly drops below 80% whereas above this value it only mildly rises. Fixing σ_p = 4.8σ, the purity remains roughly constant at ∼80 ± 5% when varying σ_f. We did see an increase in the difference between the number of positive and negative detections with increasing σ_f. However, the size of the sources above σ_f = 2.7σ drops below the 0″60 FWHM and tends to become pixel-like, hence non physical. This is because an increase of σ_f results in a reduction of the number of pixels above the floodclip threshold (σ_f) that will be associated with a given source. This parameter can be seen as a percolation criterion that sets the size of the sources in a number of pixels. Reversely reducing σ_f below 2.7σ results in adding more noise than signal and a reduction in the number of detections. We therefore decided to set σ_f to 2.7σ.

Fig. 4. Cumulative number of positive (red histogram) and negative (blue histogram) detections as a function of the σp (at a fixed σf, left panel) and σf (at a fixed σp, right panel) in units of σ. Solid black line represents the purity criterion pc define by Eq. (1), green dashed-line represents the percentage of positive detection with HST-WFC3 counterpart pHST and magenta dashed-line represents the percentage of positive detection with ZFOURGE counterpart pZFOURGE. Gray dashed-lines show the thresholds σp = 4.8σ and σf = 2.7σ and the 80% purity limit.

While we did not wish to impose a criterion on the existence of optical counterparts to define our ALMA catalog, we found that high values of σ_f not only generate the problem discussed above, but also generate a rapid drop of the fraction of ALMA detections with an HST counterpart in the Guo et al. (2013) catalog, p_HST = N_HST/N_p. N_HST is the number of ALMA sources with an HST counterpart within 0″60 (corresponding to the size of the beam). The fraction falls rapidly from around ∼80% to ∼60%, which we interpreted as being due to a rise of the proportion of spurious sources, since the faintest optical sources, for example, detected by HST-WFC3, are not necessarily associated with the faintest ALMA sources due to the negative K-correction at 1.1 mm. This rapid drop can be seen in the dashed green and dotted pink lines of Fig. 4-right. This confirms that the sources that are added to our catalog with a floodclip threshold greater than 2.7σ are most probably spurious. Similarly, we can see in Fig. 4-left that increasing the number of ALMA detections to fainter flux densities by reducing σ_p below 4.8σ leads to a rapid drop of the fraction of ALMA detections with an HST counterpart. Again there is no well-established physical reason to expect the number of ALMA detections with an optical counterpart to decrease with decreasing S/N ratio in the ALMA catalog.

As a result, we decided to set σ_p = 4.8σ and σ_f = 2.7σ to produce our catalog of ALMA detections. We note that we only discussed the existence of HST counterparts as a complementary test on the definition of the detection thresholds but our approach is not set to limit in any way our ALMA detections to galaxies with HST counterparts.

Indeed, evidence for the existence of ALMA detections with no HST-WFC3 counterparts already exist in the literature. Wang et al. (2016) identified H-dropouts galaxies, that is galaxies detected above the H-band with Spitzer-IRAC at 4.5 μm but undetected in the H-band and in the optical. The median flux density of these galaxies is F_870 μm ≃ 1.6 mJy (Wang et al., in prep.). By scaling this median value to our wavelength of 1.1 mm (the details of this computation are given in Sect. 5.4), we obtain a flux density of 0.9 mJy, close to the typical flux of our detections (median flux ∼1 mJy, see Table 3).

4. Catalog

4.1. Creation of the catalog

Using the optimal parameters of σ_p = 4.8σ and σ_f = 2.7σ described in Sect. 3, we obtained a total of 20 detections down to a flux density limit of S_1.1 mm ≈ 880 μJy that constitute our main catalog. These detections can be seen ranked by their S/N in Fig. 1. The comparison of negative and positive detections suggests the presence of 4 ± 2 (assuming a Poissonian uncertainty on the difference between the number of positive and negative detections) spurious sources in this sample.

In the following, we assume that the galaxies detected in the 0″60-mosaic are point-like. This hypothesis will later be discussed and justified in Sect. 4.5. In order to check the robustness of our flux density measurements, we compared different flux extraction methods and softwares: PyBDSM (Mohan & Rafferty 2015); Galfit (Peng et al. 2010); Blobcat (Hales et al. 2012). The peak flux value determined by Blobcat refers to the peak of the surface brightness corrected for peak bias (see Hales et al. 2012). The different results were consistent, with a median ratio of and . The fluxes measured using psf-fitting (Galfit) and peak flux measurement (Blobcat) for each galaxy are listed in Table 3. We also ran CASA fitsky and a simple aperture photometry corrected for the ALMA PSF and also found consistent results. The psf-fitting with Galfit was performed inside a box of 5 × 5″ centered on the source.

The main characteristics of these detections (redshift, flux, S/N, stellar mass, counterpart) are given in Table 3. We used redshifts and stellar masses from the ZFOURGE catalog (see Sect. 2.4.1).

We compared the presence of galaxies between the 0″60-mosaic and the 0″29-mosaic. Of the 20 detections found in the 0″60 map, 14 of them are also detected in the 0″29 map. The presence of a detection in both maps reinforces the plausibility of a detection. However, a detection in only one of these two maps may be a consequence of the intrinsic source size. An extended source is more likely to be detected with a larger beam, whereas a more compact source is more likely to be missed in the maps with larger tapered sizes and reduced point source sensitivity.

A first method to identify potential false detections was to compare our results with a deeper survey overlapping with our area of the sky. We compared the positions of our catalog sources with the positions of sources found by Dunlop et al. (2017) in the HUDF. This 1.3-mm image is deeper than our survey and reaches a σ ≃ 35 μJy (corresponding to σ = 52 μJy at 1.1 mm) but overlaps with only ∼6.5% of our survey area. The final sample of Dunlop et al. (2017) was compiled by selecting sources with S_1.3 > 120 μJy to avoid including spurious sources due to the large number of beams in the mosaics and due to their choice of including only ALMA detections with optical counterparts seen with HST.

With our flux density limit of S_1.1 mm ≈ 880 μJy any non-spurious detection should be associated with a source seen at 1.3 mm in the HUDF 1.3 mm survey, the impact of the wavelength difference being much smaller than this ratio. We detected three galaxies that were also detected by Dunlop et al. (2017), UDF1, UDF2 and UDF3, all of which have S_1.3 mm >  0.8 mJy. The other galaxies detected by Dunlop et al. (2017) have a flux density at 1.3 mm lower than 320 μJy, which makes them undetectable with our sensitivity.

We note however that we did not impose as a strict criterion the existence of an optical counterpart to our detections, whereas Dunlop et al. (2017) did. Hence if we had detected a source with no optical counterpart within the HUDF, this source may not be included in the Dunlop et al. (2017) catalog. However, as we will see, the projected density of such sources is small and none of our candidate optically dark sources fall within the limited area of the HUDF. We also note that the presence of an HST-WFC3 source within a radius of 0″6 does not necessarily imply that is the correct counterpart. As we will discuss in detail in Sect. 4.4, due to the depth of the HST-WFC3 observations and the large number of galaxies listed in the CANDELS catalog, a match between the HST and ALMA positions may be possible by chance alignment alone (see Sect. 4.4).

4.2. Supplementary catalog

After the completion of the main catalog, three sources that did not satisfy the criteria of the main catalog presented strong evidence of being robust detections. We therefore enlarged our catalog, in order to incorporate these sources into a supplementary catalog.

These three sources are each detected using a combination of σ_p and σ_f giving a purity factor greater than 80%, whilst also ensuring the existence of an HST counterpart.

The galaxy AGS21 has an S/N = 5.83 in the 0″29 tapered map, but is not detected in the 0″60 tapered map. The non-detection of this source is most likely caused by its size. Due to its dilution in the 0″60-mosaic, a very compact galaxy detected at 5σ in the 0″29-mosaic map could be below the detection limit in the 0″60-mosaic. The ratio of the mean rms of the two tapered maps is 1.56, meaning that for a point source of certain flux, a 5.83σ measurement in the 0″29-mosaic becomes 3.74σ in the 0″60-map.

The galaxy AGS22 has been detected with an S/N = 4.9 in the 0″60 tapered map (σ_p = 4.9 and σ_f = 3.1). With σ_p and σ_f values more stringent than the thresholds chosen for the main catalog, it may seem paradoxical that this source does not appear in the main catalog. With a floodclip criterion of 2.7σ, this source would have an S/N just below 4.8, excluding it from the main catalog. This source is associated with a faint galaxy that has been detected by HST-WFC3 (ID_CANDELS = 28 952) at 1.6 μm (6.6σ) at a position close to the ALMA detection (0″28). Significant flux has also been measured at 1.25 μm (3.6 σ) for this galaxy. In all of the other filters, the flux measurement is not significant (< 3σ). Due to this lack of information, it has not been possible to compute its redshift. AGS22 is not detected in the 0″29-mosaic map with p_c >  0.8. The optical counterpart of this source has a low H-band magnitude (26.8 ± 0.2 AB), which corresponds to a range for which the Guo et al. (2013) catalog is no longer complete. This is the only galaxy (except the three galaxies most likely to be spurious: AGS14, AGS16 and AGS19) that has not been detected by IRAC (which could possibly be explained by a low stellar mass). The probability of the ALMA detection being spurious, within the association radius 0″6 of a H-band source of this magnitude or brighter, is 5.5%. For these reasons, we did not consider it as spurious.

The galaxy AGS23 was detected in the 0″60 map just below our threshold at 4.8σ, with a combination σ_p = 4.6 and σ_f = 2.9 giving a purity criterion greater than 0.9. This detection is associated with an HST-WFC3 counterpart. It is for these two reasons that we include this galaxy in the supplementary catalog. The photometric redshift (z = 2.36) and stellar mass (10^11.26 M_⊙) both reinforce the plausibility of this detection.

4.3. Astrometric correction

The comparison of our ALMA detections with HST (Sect. 4.1) in the previous section was carried out after correcting for an astrometric offset, which we outline here. In order to perform the most rigorous counterpart identification and take advantage of the accuracy of ALMA, we carefully investigated the astrometry of our images. Before correction, the galaxy positions viewed by HST were systematically offset from the ALMA positions. This offset has already been identified in previous studies (e.g., Maiolino et al. 2015; Rujopakarn et al. 2016; Dunlop et al. 2017).

In order to quantify this effect, we compared the HST source positions with detections from the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS). This survey has the double advantage to cover a large portion of the sky, notably the GOODS-South field, and to observe the sky at a wavelength similar to HST-WFC3. We used the Pan-STARRS DR1 catalog provided by Flewelling et al. (2016) and also included the corresponding regions issued from the Gaia DR1 (Gaia Collaboration 2016).

Cross-matching was done within a radius of 0″5. In order to minimize the number of false identifications, we subtracted the median offset between the two catalogs from the Guo et al. (2013) catalog positions, after the first round of matching. We iterated this process three times. In this way, 629 pairs were found over the GOODS-South field.

To correct for the median offset between the HST and ALMA images, the HST image coordinates needed to be corrected by −96 ± 113 mas in right ascension, α, and 261 ± 125 mas in declination, δ, where the uncertainties correspond to the standard deviation of the 629 offset measurements. This offset is consistent with that found by Rujopakarn et al. (2016) of Δα = −80 ± 110 mas and Δδ =  260 ± 130 mas. The latter offsets were calculated by comparing the HST source positions with 2MASS and VLA positions. In all cases, it is the HST image that presents an offset, whereas ALMA, Pan-STARRS, Gaia, 2MASS and VLA are all in agreement. We therefore deduced that it is the astrometric solution used to build the HST mosaic that introduced this offset. As discussed in Dickinson et al. (in prep.), the process of building the HST mosaic also introduced less significant local offsets, that can be considered equivalent to a distortion of the HST image. These local offsets are larger in the periphery of GOODS-South than in the center, and close to zero in the HUDF field. The local offsets can be considered as a distortion effect. The offsets listed in Table 2 include both effects, the global and local offsets. The separation between HST and ALMA detections before and after offset correction, and the individual offsets applied for each of the galaxies are indicated in Table 2 and can be visualized in Fig. 5. We applied the same offset corrections to the galaxies listed in the ZFOURGE catalog.

Fig. 5. Positional offset (RAHST–RAALMA, DecHST–DecALMA) between HST and ALMA before (red crosses) and after (blue crosses) the correction of both a global systematic offset and a local offset. The black dashed circle corresponds to the cross-matching limit radius of 0″6. The gray dashed circles show a positional offset of 0″2 and 0″4 respectively. The magenta lines indicate the HST galaxies previously falsely associated with ALMA detections.

This accurate subtraction of the global systematic offset, as well as the local offset, does not however guarantee a perfect overlap between ALMA and HST emission. The location of the dust emission may not align perfectly with the starlight from a galaxy, due to the difference in ALMA and HST resolutions, as well as the physical offsets between dust and stellar emission that may exist. In Fig. 6, we show the ALMA contours (4–10σ) overlaid on the F160W HST-WFC3 images after astrometric correction. In some cases (AGS1, AGS3, AGS6, AGS13, AGS21 for example), the position of the dust radiation matches that of the stellar emission; in other cases, (AGS4, AGS17 for example), a displacement appears between both two wavelengths. Finally, in some cases (AGS11, AGS14, AGS16, and AGS19) there are no optical counterparts. We will discuss the possible explanations for this in Sect. 7.

Fig. 6. Postage stamps of 1.8 × 1.8 arcsec. ALMA contours (4, 4.5 then 5–10-σ with a step of 1-σ) at 1.1 mm (white lines) are overlaid on F160W HST/WFC3 images. The images are centered on the ALMA detections. The shape of the synthesized beam is given in the bottom left corner. Astrometry corrections described in Sect. 4.3 have been applied to the HST images. In some cases (AGS1, AGS3, AGS6, AGS13, AGS21 for example), the position of the dust radiation matches that of the stellar emission; in other cases, (AGS4, AGS17 for example), a displacement appears between both two wavelengths. Finally, in some cases (AGS11, AGS14, AGS16, and AGS19) there are no optical counterparts. We will discuss the possible explanations for this in Sect. 7.

4.4. Identification of counterparts

We searched for optical counterparts in the CANDELS/GOODS-South catalog, within a radius of 0″6 from the millimeter position after applying the astrometric corrections to the source positions described in Sect. 4.3. The radius of the cross-matching has been chosen to correspond to the synthesized beam (0″60) of the tapered ALMA map used for galaxy detection. Following Condon (1997), the maximal positional accuracy of the detection in the 1.1 mm map is given by θ_beam/(2 × S/N). In the 0″60-mosaic, the positional accuracy therefore ranges between 26.5 mas and 62.5 mas for our range of S/N (4.8–11.3), corresponding to physical sizes between 200 and 480 pc at z = 3.

Despite the high angular resolution of ALMA, the chance of an ALMA-HST coincidence is not negligible, because of the large projected source density of the CANDELS/GOODS-South catalog. Figure 7 shows a Monte Carlo simulation performed to estimate this probability. We separate here the deeper Hubble Ultra Deep Field (blue histogram) from the rest of the CANDELS-deep area (orange histogram). We randomly defined a position within GOODS-South and then measured the distance to its closest HST neighbor using the source positions listed in Guo et al. (2013). We repeated this procedure 100 000 times inside and outside the HUDF. The probability for a position randomly selected in the GOODS-South field to fall within 0.6 arcsec of an HST source is 9.2% outside the HUDF, and 15.8% inside the HUDF. We repeated this exercise to test the presence of an IRAC counterpart with the Ashby et al. (2015) catalog (green histogram). The probability to randomly fall on an IRAC source is only 2.1%.

Fig. 7. Probability of a randomly selected position in the area defined by this survey to have at least one HST (blue, orange) or IRAC (green) neighbor as a function of distance. We computed this probability by Monte Carlo simulation using the distribution of galaxies listed in the CANDELS/GOODS-South catalog. Due to the presence of the HUDF within the GOODS-South field, we cannot consider that the density of HST galaxies is uniform, and we consider these two fields separately (blue inside and orange outside the HUDF).

With the detection threshold determined in Sect. 3, 80% of the millimeter galaxies detected have an HST-WFC3 counterpart, and four galaxies remain without an optical counterpart. We cross-matched our detections with the ZFOURGE catalog.

Figure 8 shows 3″5 × 3″5 postage stamps of the ALMA-detected galaxies, overlaid with the positions of galaxies from the CANDELS/GOODS-South catalog (magenta double crosses), ZFOURGE catalog (white circles) or both catalogs (sources with an angular separation lower than 0″4, blue circles). These are all shown after astrometric correction. Based on the ZFOURGE catalog, we found optical counterparts for one galaxy that did not have an HST counterpart: AGS11, a photometric redshift has been computed in the ZFOURGE catalog for this galaxy.

Fig. 8. ALMA 1.1 mm continuum maps for the 23 detections tapered at 0.60 arcsec. Each 3″5 × 3″5 image is centered on the position of the ALMA detection. Cyan double crosses show sources from the GOODS-S CANDELS catalog. White circles show sources from the ZFOURGE catalog. Blue circles show common sources from both optical catalogs (sources with an angular separation lower than 0″4). The shape of the synthesized beam is given in the bottom left corner.

The redshifts of AGS4 and AGS17 as given in the CANDELS catalog are unexpectedly low (z = 0.24 and z = 0.03, respectively), but the redshifts for these galaxies given in the ZFOURGE catalog (z = 3.76 and z = 1.85, respectively) are more compatible with the expected redshifts for galaxies detected with ALMA. These galaxies, missed by the HST or incorrectly listed as local galaxies are particularly interesting galaxies (see Sect. 7). AGS6 is not listed in the ZFOURGE catalog, most likely because it is close (<0″7) to another bright galaxy (ID_CANDELS = 15 768). These galaxies are blended in the ZFOURGE ground-based K_s-band images. AGS6 has previously been detected at 1.3 mm in the HUDF, so we adopt the redshift and stellar mass found by Dunlop et al. (2017). The consensus CANDELS zphot from Santini et al. (2015) is z = 3.06 (95% confidence: 2.92 <  z <  3.40), consistent with the value in Dunlop et al. (2017).

4.5. Galaxy sizes

Correctly estimating the size of a source is an essential ingredient for measuring its flux. As a first step, it is imperative to know if the detections are resolved or unresolved. In this section, we discuss our considerations regarding the sizes of our galaxies. The low number of galaxies with measured ALMA sizes in the literature makes it difficult to constrain the size distribution of dust emission in galaxies. Recent studies (e.g., Barro et al. 2016; Rujopakarn et al. 2016; Elbaz et al. 2018; Ikarashi et al. 2017; Fujimoto et al. 2017) with sufficient resolution to measure ALMA sizes of galaxies suggest that dust emission takes place within compact regions of the galaxy.

Two of our galaxies (AGS1 and AGS3) have been observed in individual pointings (ALMA Cycle 1; P.I. R. Leiton, presented in Elbaz et al. 2018) at 870 μm with a long integration time (40–50 min on source). These deeper observations give more information on the nature of the galaxies, in particular on their morphology. Due to their high S/N (∼100) the sizes of the dust emission could be measured accurately: R_1/2maj = 120 ± 4 and 139 ± 6 mas for AGS1 and AGS3 respectively, revealing extremely compact star-forming regions corresponding to circularized effective radii of ∼1 kpc at redshift z ∼ 2. The Sersic indices are 1.27 ± 0.22 and 1.15 ± 0.22 for AGS1 and AGS3 respectively: the dusty star-forming regions therefore seem to be disk-like. Based on their sizes, their stellar masses (> 10¹¹ M_⊙), their SFRs (> 10³ M_⊙ yr⁻¹) and their redshifts (z ∼ 2), these very compact galaxies are ideal candidate progenitors of compact quiescent galaxies at z ∼ 2 (Barro et al. 2013; Williams et al. 2014; van der Wel et al. 2014; Kocevski et al. 2017, see also Elbaz et al. 2018).

Size measurements of galaxies at (sub)millimeter wavelengths have previously been made as part of several different studies. Ikarashi et al. (2015) measured sizes for 13 AzTEC-selected SMGs. The Gaussian FWHM range between 0″10 and 0″38 with a median of at 1.1 mm. Simpson et al. (2015a) derived a median intrinsic angular size of FWHM = 0″30 ± 0″04 for their 23 detections with a S/N >  10 in the Ultra Deep Survey (UDS) for a resolution of 0″3 at 870 μm. Tadaki et al. (2017) found a median FWHM of 0″11  ±  0.02 for 12 sources in a 0″2-resolution survey at 870 μm. Barro et al. (2016) use a high spatial resolution (FWHM ∼ 0″14) to measure a median Gaussian FWHM of 0″12 at 870 μm, with an average Sersic index of 1.28. For Hodge et al. (2016), the median major axis size of the Gaussian fit is FWHM = 0″42 ± 0″04 with a median axis ratio b/a = 0.53 ± 0.03 for 16 luminous ALESS SMGs, using high-resolution (∼0″16) data at 870 μm. Rujopakarn et al. (2016) found a median circular FWHM at 1.3 mm of 0″46 from the ALMA image of the HUDF (Dunlop et al. 2017). González-López et al. (2017) studied 12 galaxies at S/N ≥ 5, using 3 different beam sizes (0″63 × 0″49), (1″52 × 0″85) and (1″22 × 1″08). They found effective radii spanning < 0″05 to 0″37 ± 0″21 in the ALMA Frontier Fields survey at 1.1 mm. Ikarashi et al. (2017) obtained ALMA millimeter-sizes of 0″08–0″68 (FWHM) for 69 ALMA-identified AzTEC SMGs with an S/N greater than 10. These galaxies have a median size of 0″31. These studies are all broadly in agreement, revealing compact galaxy sizes in the sub(millimeter) regime of typically 0″3 ± 0″1.

Size measurements require a high S/N detection to ensure a reliable result. The S/N range of our detections is 4.8–11.3. Following Martí-Vidal et al. (2012), the reliable size measurement limit for an interferometer is:

(2)

where λ_c is the value of the log-likelihood, corresponding to the cutoff of a Gaussian distribution to have a false detection and β is a coefficient related to the intensity profile of the source model and the density of the visibilities in Fourier space. This coefficient usually takes values in the range 0.5–1. We assumed λ_c = 3.84 corresponding to a 2σ cut-off, and β = 0.75. For θ_beam = 0″60 and a range of S/N between 4.8 and 11.3, the minimum detectable size (FWHM) therefore varies between 0″16 and 0″24. Using the 0″60-mosaic map, the sizes of a large number of detections found in previous studies could therefore not be reliably measured.

To quantitatively test if the millimeter galaxies are resolved in our survey we performed several tests. The first test was to stack the 23 ALMA-detections and compare the obtained flux profile with the profile of the PSF. However, in the mosaic map, each slice has its own PSF. We therefore also needed to stack the PSFs at these 23 positions in order to obtain a global PSF for comparison. Figure 9 shows the different PSFs used in this survey in the 0″60-mosaic. The FWHM of each PSF is identical, the differences are only in the wings. The stack of the 23 PSFs for the 23 detections and the result of the source stacking in the 0″60-mosaic is shown in Fig. 10. The flux of each detection is normalized so that all sources have the same weight, and the stacking is not skewed by the brightest sources.

Fig. 9. East–west profile of the PSFs corresponding to the six different parallel slices composing the ALMA image in the 0″60-mosaic (see Fig. 1).

Fig. 10. Comparison between the stacked PSF (black solid line) and the stack of the 23 ALMA-detections (black dashed line) in the 0″60-mosaic. As each slice has a specific PSF, we stack the PSF corresponding to the position of each detection. The fluxes of each detection have been normalized, so that the brightest sources do not skew the results. Fluxes of the PSF and ALMA detections are normalized to 1. Flux profiles are taken across the East–west direction. The result is consistent with unresolved or marginally resolved sources at this resolution. The insert in the top-right corner shows the same procedure for the 15 sources detected in the 0″29-mosaic (see Table 3).

Size stacking to measure the structural parameters of galaxies is at present a relatively unexplored area. This measurement could suffer from several sources of bias. The uncertainties on the individual ALMA peak positions could increase the measured size in the stacked image, for example. On the other hand, due to the different inclination of each galaxy, the stacked galaxy could appear more compact than the individual galaxies (e.g. Hao et al. 2006; Padilla & Strauss 2008; Li et al. 2016). Alternatively, some studies (e.g. van Dokkum et al. 2010) indicate that size stacking gives reasonably accurate mean galaxy radii. In our case, the result of the size stacking is consistent with unresolved sources or marginally resolved at this resolution which corresponds to a physical diameter of 4.6 kpc at z = 3.

The second test was to extract the flux for each galaxy using PSF-fitting. We used Galfit (Peng et al. 2010) on the 0″60-mosaic. The residuals of this PSF-extraction are shown for the 6 brightest galaxies in Fig. 11. The residuals of 21/23 detections do not have a peak greater than 3σ in a radius of 1″ around the source. Only sources AGS10 and AGS21 present a maximum in the residual map at ∼3.1σ.

Fig. 11. 10″ × 10″ postage stamps, centered on the galaxy detections. Left to right: source in the 0″60-mosaic map, and residuals obtained after PSF, Gaussian and Sersic flux fitting. The residuals are very similar between the three different extraction methods. Only the 6 brightest galaxies are shown.

We compared the PSF flux extraction method with Gaussian and Sersic shapes. As our sources are not detected with a particularly high S/N, and in order to limit the number of degrees of freedom, the Sersic index was frozen to n = 1 (exponential disk profile, in good agreement with Hodge et al. 2016 and Elbaz et al. 2018 for example), assuming that the dust emission is disk-like. Figure 11 shows the residuals for the three different extraction profiles. The residuals are very similar between the point source, Gaussian and Sersic profiles, suggesting that the approximation that the sources are not resolved is appropriate, and does not result in significant flux loss. We also note that, for several galaxies, due to large size uncertainties, the Gaussian and Sersic fits give worse residuals than the PSF fit (AGS4 for example).

For the third test, we took advantage of the different tapered maps. We compared the peak flux for each detection between the 0″60-mosaic map and the 0″29-mosaic map. The median ratio is . This small decrease, of only 10% in the peak flux density between the two tapered maps suggests that the flux of the galaxies is only slightly more resolved in the 0″29-mosaic map.

In order to test the impact of our hypothesis that the sources can be considered as point-like in the mosaic tapered at 0″60, we fit their light profiles with a circular Gaussian in the uv-plane using uvmodelfit in CASA (we also tested the use of an asymmetric Gaussian but the results remained similar although with a lower precision due to the larger number of free parameters in the fit). The sizes that we obtained confirmed our hypothesis that our galaxies are particularly compact since 85% of the sources (17 out of 20 robust detections) exhibit a FWHM smaller than 0″25 (in other words the half-light radius is twice smaller than this value). The median size of our sample of 20 galaxies is 0″18 (see the distribution of sizes in Fig. 12). This analysis shows that two sources are outliers with sizes of 0″41 ± 0″03 and 0″50 ± 0″08, for AGS17 and AGS18 respectively. For these two sources, the assumption of point-like sources is not valid and leads to an underestimate of the actual flux densities by a factor of 2.3 and 1.7 respectively. This correction has been applied to the list of peak flux densities provided in Table 3.

Fig. 12. Size distribution histogram for the 20 robust detections. These sizes are computed by fitting the ALMA detections with a circular Gaussian in the uv-plane using uvmodelfit in CASA. 85% of the sources exhibit a FWHM smaller than 0″25.

Having performed these tests, we concluded that for all of the detections, except AGS17 and AGS18, the approximation that these sources appear point-like in the 0″60-mosaic map is justified. For the two remaining sources, we applied a correction given above. Our photometry was therefore performed under this assumption.

5. Number counts

5.1. Completeness

We assessed the accuracy of our catalog by performing completeness tests. The completeness is the probability for a source to be detected in the map given factors such as the depth of the observations. We computed the completeness of our observations using Monte Carlo simulations performed on the 0″60-mosaic map. We injected 50 artificial sources in each slice. Each source was convolved with the PSF and randomly injected on the dirty map tapered at 0″60. In total, for each simulation run, 300 sources with the same flux were injected into the total map. In view of the size of the map, the number of independent beams and the few number of sources detected in our survey, we can consider, to first order, that our dirty map can be used as a blank map containing only noise, and that the probability to inject a source exactly at the same place as a detected galaxy is negligible. The probability that at least two point sources, randomly injected, are located within the same beam (p_b) is:

(3)

where N_b is the number of beams and n is the number of injected sources. For each one of the six slices of the survey, we count ∼100 000 independent beams. The probability of having source blending for 50 simulated sources in one map is ∼1%.

We then counted the number of injected sources detected with σ_p = 4.8σ and σ_f = 2.7σ, corresponding to the thresholds of our main catalog. We injected 300 artificial sources of a given flux, and repeated this procedure 100 times for each flux density. Our simulations cover the range S _ν = 0.5–2.4 mJy in steps of 0.1 mJy. Considering the resolution of the survey, it would be reasonable to expect that a non-negligible number of galaxies are not seen as point sources but extended sources (see Sect. 4.5). We simulated different sizes of galaxies with Gaussian FWHM between 0″2 and 0″9 in steps of 0″1, as well as point-source galaxies, to better understand the importance of the galaxy size in the detectability process. We matched the recovered source with the input position within a radius of 0″6.

Figure 13 shows the resulting completeness as a function of input flux, for different FWHM Gaussian sizes convolved by the PSF and injected into the map.

Fig. 13. Median source detection completeness for simulated point-like and Gaussian galaxies as a function of integrated flux, for different FWHMs (see figure’s legend). The shaded regions correspond to the standard deviation of 100 runs, each containing 300 simulated sources.

As a result of our simulations, we determined that at 1.2 mJy, our sample is 94 ± 1% complete for point sources. This percentage drastically decreases for larger galaxy sizes. For the same flux density, the median detection rate drops to 61 ± 3% for a galaxy with FWHM ∼ 0″3, and to 9 ± 1% for a FWHM ∼ 0″6 galaxy. This means, that for a galaxy with an intrinsic flux density of 1.2 mJy, we are more than ten times more likely to detect a point source galaxy than a galaxy with FWHM ∼ 0″6.

The size of the millimeter emission area plays an essential role in the flux measurement and completeness evaluation. We took the hypothesis that ALMA sizes are 1.4 times smaller than the size measured in HST H-band (as derived by Fujimoto et al. 2017 using 1034 ALMA galaxies). We are aware that this size ratio is poorly constrained at the present time, but such relations have been observed in several studies (see Sect. 4.5). For example, of the 12 galaxies presented by Laporte et al. (2017), with fluxes measured using ALMA at 1.1 mm (González-López et al. 2017), seven of them have a size measured by HST F140W/WFC3 similar to the size measured in the ALMA map. On the other hand, for the remaining five galaxies, their sizes are approximately two times more compact at millimeter wavelengths than at optical wavelengths. This illustrates the dispersion of this ratio.

5.2. Effective area

As the sensitivity of our 1.1 mm ALMA map is not uniform, we defined an effective area where a source with a given flux can be detected with an S/N >  4.8σ, as shown in Fig. 14. Our map is composed of six different slices – one of them, slice B, presents a noise 30% greater than the mean of the other five, whose noise levels are comparable. The total survey area is 69.46 arcmin², with 90% of the survey area reaching a sensitivity of at least 1.06 mJy beam⁻¹. We considered the relevant effective area for each flux density in order to compute the number counts. We considered the total effective area over all slices in the number counts computation.

Fig. 14. Effective area as a function of flux density, where a source with a given flux can be detected with an S/N >  4.8σ. Ninety percent of the survey area reaches a sensitivity of at least 1.06 mJy beam−1.

5.3. Flux boosting and Eddington bias

In this section, we evaluate the effect of flux boosting. Galaxies detected with a relatively low S/N tend to be boosted by noise fluctuations (see Hogg & Turner 1998; Coppin et al. 2005; Scott et al. 2002). To estimate the effect of flux boosting, we used the same set of simulations that we used for completeness estimations.

The results of our simulations are shown in Fig. 15. The boosting effect is shown as the ratio between the input and output flux densities as a function of the measured S/N. For point sources, we observed the well-known flux boosting effect for the lowest S/Ns. This effect is not negligible for the faintest sources in our survey. At 4.8σ, the flux boosting is ∼15%, and drops below 10% for an S/N greater than 5.2. We estimated the de-boosted flux by dividing the measured flux by the median value of the boosting effect as a function of S/N (red line in Fig. 15).

Fig. 15. Flux boosting as a function of measured S/N estimated from simulations. The median of the boosting is shown by a solid red line. The 1σ confidence intervals (dashed red lines) are overplotted. The solid black horizontal line corresponds to Fout = Fin (see text for details). We used the same set of simulations that we used for the completeness analysis.

We also corrected for the effects of the Eddington bias (Eddington 1913). As sources with lower luminosities are more numerous than bright sources, Gaussian distributed noise gives rise to an overestimation of the number counts in the lowest flux bins. We simulated a realistic number of sources (the slope of the number counts were computed using the coefficients given in Table 5) and added Gaussian noise to each simulated source. The correction factor for each flux bin was therefore the ratio between the flux distribution before and after adding the noise.

5.4. Cumulative and differential number counts

We used sources with a S/N greater than 4.8 from the main catalog to create cumulative and differential number counts. We needed to take into account the contamination by spurious sources, completeness effects, and flux boosting in order to compute these number counts.

The contribution of a source with flux density S_i ± dS_i to the cumulative number count is given by:

(4)

where p_c(S_i) is the purity criterion as defined in Eq. (1) at the flux density S_i, A_eff(S_i) and are the effective area and the completeness for the flux interval dS_i, as shown in Figs. 13 and 14. The completeness is strongly correlated with the sizes of the galaxies. To estimate the completeness, galaxies that do not have measured sizes in the H-band (van der Wel et al. 2012) were considered as point sources, otherwise we used (see Sect. 5.1).

The cumulative number counts are given by the sum over all of the galaxies with a flux density higher than S:

(5)

Errors are computed by Monte-Carlo simulations, added in quadrature to the Poisson uncertainties. The derived number counts are provided in Table 4. AGS19 is located at a position where the noise is artificially low, and has therefore not been taken into account.

In Fig. 16, we compare our results with previous studies (Lindner et al. 2011; Scott et al. 2012; Karim et al. 2013; Hatsukade et al. 2013, 2016; Simpson et al. 2015b; Oteo et al. 2016; Aravena et al. 2016; Fujimoto et al. 2016; Umehata et al. 2017; Geach et al. 2017; Dunlop et al. 2017). To standardize these previous studies, the different flux densities are scaled to 1.1 mm using a Modified Black Body (MBB) model, assuming a dust emissivity index β = 1.5 (e.g., Gordon et al. 2010), a dust temperature T_d = 35 K (e.g. Chapman et al. 2005; Kovács et al. 2006; Coppin et al. 2008), and a redshift of z = 2.5 (e.g., Wardlow et al. 2011; Yun et al. 2012). These values have also been chosen to be consistent with Hatsukade et al. (2016). The different fluxes were therefore scaled to 1.1 mm using the relations S_1.1 mm/S_1.2 mm = 1.29, S_1.1 mm/S_1.3 mm = 1.48 and S_1.1 mm/S_870 μm = 0.56. It is a real challenge to standardize these previous studies because instruments, observational techniques or resolution often vary between studies. Some of these counts have been computed from individual pointings, by brightness selection, or by serendipitous detections. Observations with a single dish or a low resolution can also overestimate the number counts for the brightest galaxies, because of blending effects (see Ono et al. 2014). Another non-negligible source of error can come from an inhomogeneous distribution of bright galaxies. An underdensity by a factor of two of submillimeter galaxies with far infrared luminosities greater than 2 × 10¹² L_⊙ in the extended Chandra deep field south (ECDFS) compared to other deep fields has been revealed by Weiß et al. (2009).

Fig. 16. 1.1-mm cumulative (left panel) and differential (right panel) number counts derived using the corrections described in Sect. 5.4, for the sources detected at > 4.8σ in the main catalog. AGS19 is located at a position where the noise is artificially low, and has therefore not been taken into account. Previous (sub)millimeter cumulative number counts are also shown (Lindner et al. 2011; Scott et al. 2012; Karim et al. 2013; Hatsukade et al. 2013, 2016; Ono et al. 2014; Simpson et al. 2015b; Oteo et al. 2016; Carniani et al. 2015; Aravena et al. 2016; Fujimoto et al. 2016; Umehata et al. 2017; Geach et al. 2017; Dunlop et al. 2017). The different fluxes are scaled to 1.1 mm flux densities using S1.1 mm/S1.2 mm = 1.29, S1.1 mm/S1.3 mm = 1.48, S1.1 mm/S870 μm = 0.56. From the Umehata et al. (2017) study, we used only sources which do not have z = 3.09, (which means we are excluding the protocluster members). Results from single-dish surveys are shown with unfilled pentagon markers and are only indicative, they are not considered for model fitting. The gray curve shows the best-fit Schechter function (with 1-σ) uncertainties, the red curve shows the best-fit DPL function (with 1-σ).

Despite those potential caveats, the results from our ALMA survey in the GOODS-South field are in good agreement with previous studies for flux densities below 1 mJy. For values above this flux density, two different trends coexist as illustrated in Fig. 16: our counts are similar to those found by Karim et al. (2013), but below the trend characterized by Scott et al. (2012). These two previous studies have been realized under different conditions. The effects of blending, induced by the low resolution of a single dish observation, as with Scott et al. (2012), tend to overestimate the number counts at the bright-end (Ono et al. 2014; Karim et al. 2013; Béthermin et al. 2017). We indicate these points on the Fig. 16 on an indicative basis only.

The differences in wavelength between the different surveys, even after applying the scaling corrections above, can also induce scatter in the results, especially for wavelengths far from 1.1 mm. The cumulative source counts from the 20 detections in this study and the results from other multidish blank surveys are fitted with a Double Power Law (DPL) function (e.g., Scott et al. 2002) given by:

(6)

and a modified Schechter (Schechter 1976) function (e.g., Knudsen et al. 2008):

(7)

where N₀ is the normalization, S₀ the characteristic flux density and α is the faint-end slope. β is the bright-end slope of the number of counts in Eq. (6). We used a least squares method with the trust region reflective algorithm for these two fitted-functions. The best-fit parameters are given in Table 5.

One of the advantages of using differential number counts compared to cumulative number counts is the absence of correlation of the counts between the different bins. However, the differential number counts are sensitive to the lower number of detections per flux density bin. Here we used Δlog S_ν = 0.2 dex flux density bins.

We compare our results with an empirical model that predicts the number counts at far-IR and millimeter wavelengths, developed by Béthermin et al. (2017). This simulation, called SIDES (Simulated Infrared Dusty Extragalactic Sky), updates the Béthermin et al. (2012) model. These predictions are based on the redshift evolution of the galaxy properties, using a two star-formation mode galaxy evolution model (see also Sargent et al. 2012).

The Béthermin et al. (2017) prediction is in good agreement with the number counts derived in this study, for the two bins with the lowest fluxes. For the highest-flux bin, the model is slightly above the data (∼1σ above the best Schechter fit for fluxes greater than 1 mJy). However, both the Béthermin et al. (2017) model and our data points are below the single-dish measurements for fluxes greater than 1 mJy. This disagreement between interferometric and single-dish counts is expected, because the boosting of the flux of single-dish sources by their neighbor in the beam (Karim et al. 2013; Hodge et al. 2013; Scudder et al. 2016). Béthermin et al. (2017) derived numbers counts from a simulated single-dish map based on their model and found a nice agreement with single-dish data, while the intrinsic number counts in the simulation are much lower and compatible with our interferometric study.

Cosmic variance was not taken into account in the calculation of the errors. Above z = 1.8 and up to the redshift of the farthest galaxy in our catalog at z = 4.8, the strong negative K-correction at this wavelength ensures that the selection of galaxies is not redshift-biased. The cosmic variance, although significant for massive galaxies in a small solid angle, is counterbalanced by the negative K-correction, which makes the redshift interval of our sources (Δz = 3 in Eq. (12) in Moster et al. 2011) relatively large, spanning a comoving volume of 1400 Gpc³. Based on Moster et al. (2011), the cosmic variance for our sources is ∼15%, which does not significantly affect the calculation of the errors on our number counts.

5.5. Contribution to the cosmic infrared background

The extragalactic background light (EBL) is the integrated intensity of all of the light emitted throughout cosmic time. Radiation re-emitted by dust comprises a significant fraction of the EBL, because this re-emitted radiation, peaking around 100 μm, has an intensity comparable to optical background (Dole et al. 2006). The contribution of our ALMA sources to the EBL is derived by integrating the derived number counts down to a certain flux density limit. Using the 20 (>4.8σ) sources detected, we computed the fraction of the 1.1 mm EBL resolved into discrete sources. The integrated flux density is given by:

(8)

We used the set of parameters given in Table 5 on the differential number counts. We compared our results with observations from the far infrared absolute spectrophotometer (FIRAS) on the cosmic background explorer (COBE), knowing that uncertainties exist on the COBE measurements (e.g., Yamaguchi et al. 2016). We used the equation given in Fixsen et al. (1998) to compute the total energy of the EBL:

(9)

where ν₀ = 100 cm⁻¹, and P_ν is the familiar Planck function with I_ν in erg s⁻¹ cm⁻¹ Hz⁻¹ sr⁻¹. From this equation, we found that at 1.1 mm, the energy of the EBL is 2.87 nW m⁻² sr⁻¹. From Eq. (8) we can estimate the integrated EBL light. Figure 17 shows this total integrated flux density. For our data, the lowest flux density bin for differential counts S_lim is 0.8 mJy, and we extrapolate to lower flux densities. We have resolved only of the EBL into individual galaxies at 0.8 mJy. This result is in good agreement with studies such as Fujimoto et al. (2016). In order to have the majority of the EBL resolved (e.g., Hatsukade et al. 2013; Ono et al. 2014; Carniani et al. 2015; Fujimoto et al. 2016), we would need to detect galaxies down to 0.1 mJy (about 50 % of the EBL is resolved at this value).

Fig. 17. Resolved 1.1 mm EBL computed from the best-fit Schechter function to the differential number counts. The green, red, purple and brown lines are from the number counts estimated by Fujimoto et al. (2016), Ono et al. (2014), Hatsukade et al. (2016), and Carniani et al. (2015) respectively. The blue line and shaded region show the results from this work and the associated uncertainty. The solid lines represent the model above the detection limits, and the dashed lines show the extrapolation below these limits. The gray shaded region shows the 1.1 mm cosmic infrared background measured by COBE (Fixsen et al. 1998).

The extrapolation of the integrated flux density below S_lim suggests a flattening of the number counts. The population of galaxies that dominate this background is composed of the galaxies undetected in our survey, with a flux density below our detection limit.

6. Galaxy properties

We now study the physical properties of the ALMA detected sources, taking advantage of the wealth of ancillary data available for the GOODS-South field.

6.1. Redshift distribution

Among the 17 ALMA detected sources for which redshifts have been computed, six have a spectroscopic redshift (AGS1, AGS2, AGS3, AGS9, AGS12, AGS13 and AGS18) determined by Kurk et al. (2013), and recently confirmed by Barro et al. (2017), Momcheva et al. (2016), Vanzella et al. (2008), Mobasher (priv. comm.), Inami et al. (2017), Kriek et al. (2008), and Dunlop et al. (2017) – from a private communication of Brammer – respectively. The redshift distribution of these 17 ALMA sources is presented in Fig. 18, compared to the distributions of four other deep ALMA blind surveys (Dunlop et al. 2017; Aravena et al. 2016; González-López et al. 2017; Ueda et al. 2018). Of the 17 sources, 15 are in the redshift range z = 1.9–3.8. Only two galaxies (AGS4 and AGS11) have a redshift greater than 4 (z_phot = 4.32 and 4.82 respectively). We discuss these galaxies further in Sect. 7. The mean redshift of the sample is z = 3.03 ± 0.17, where the error is computed by bootstrapping. This mean redshift is significantly higher than those found by Dunlop et al. (2017), Aravena et al. (2016), González-López et al. (2017), and Ueda et al. (2018) who find distributions peaking at 2.13, 1.67, 1.99 and 2.28 respectively. The median redshift of our sample is 2.92 ± 0.20, which is a little higher than the value expected from the models of Béthermin et al. (2015), which predict a median redshift of 2.5 at 1.1 mm, considering our flux density limit of ∼874 μJy (4.8σ).

Fig. 18. Redshift distributions (photometric or spectroscopic) for millimeter-selected galaxies. The blue solid line shows the redshift distribution of our ALMA GOODS-South blind survey. The green dashed line shows the Hubble Ultra Deep Field Survey redshift distribution (Dunlop et al. 2017), the black dash-dotted line shows the ASPECS sample (Aravena et al. 2016), the red dotted line shows the ALMA Frontier Fields survey (González-López et al. 2017) and the yellow dotted line shows the ASAGAO survey (Ueda et al. 2018). Short colored lines at the top of the figure indicate the median redshifts for these four studies.

Our limiting sensitivity is shallower than that of previous blind surveys: 0.184 mJy here compared with 13 μJy in Aravena et al. (2016), 35 μJy in Dunlop et al. (2017), (55–71) μJy in González-López et al. (2017) and 89 μJy in Ueda et al. (2018). However, our survey covers a larger region on the sky: 69 arcmin² here, compared to 1 arcmin², 4.5 arcmin², 13.8 arcmin² and 26 arcmin² for these four surveys respectively. The area covered by our survey is therefore a key parameter in the detection of high redshift galaxies due to a tight link between 1.1 mm luminosity and stellar mass as, we will show in the next section. The combination of two effects: a shallower survey allowing us to detect brighter SMGs, which are more biased toward higher redshifts (e.g., Pope et al. 2005), as well as a larger survey allowing us to reach more massive galaxies, enables us to open the parameter space at redshifts greater than 3, as shown in Fig. 18. This redshift space is partly or totally missed in smaller blind surveys.

We emphasize that the two HST-dark galaxies (see Sect. 7) for which the mass and redshift could be determined (AGS4 and AGS11) are the two most distant galaxies in our sample, with redshifts greater than 4.

6.2. Stellar masses

Over half (10/17) of our galaxies have a stellar mass greater than 10¹¹ M_⊙ (median mass of M_⋆ = 1.1 × 10¹¹ M_⊙). The population of massive and compact star-forming galaxies at z ∼ 2 has been documented at length (e.g., Daddi et al. 2005; van Dokkum et al. 2015), but their high redshift progenitors are to-date poorly detected in the UV. Our massive galaxies at redshifts greater than 3 might therefore give us an insight into these progenitors.

Figure 19 shows the stellar mass as a function of redshift for all of the UVJ active galaxies, listed in the ZFOURGE catalog, in our ALMA survey field of view. Star-forming galaxies (SFGs) have been selected by a UVJ color–color criterion as given by Williams et al. (2009) and applied at all redshifts and stellar masses as suggested by Schreiber et al. (2015):

(10)

Fig. 19. Stellar mass versus redshift for the galaxies detected in our ALMA GOODS-South blind survey (red points). For comparison, the distribution of all of the galaxies, listed in the ZFOURGE catalog, in the same field of view is given in blue. Only UVJ active galaxies are shown. The two HST-dark galaxies for which we have redshifts (AGS4 and AGS11) are represented by open circles. The redshift of AGS11 is however uncertain. The green dashed line shows the position that would be occupied by a typical star-forming galaxy – lying on the median of the SFR–M⋆ star-formation main sequence (MS) – that would produce a 1.1 mm flux density equal to our average detection limit of 0.88 mJy (4.8-σ) using the spectral energy distribution (SED) library of Schreiber et al. (2018). The dotted line illustrates the position of galaxies 3 times above the MS using the appropriate SEDs from the same library. Galaxies hosting an AGN that are undetected or detected by ALMA are identified with black dots and yellow stars respectively. Inside the black dashed rectangle, 50% of the galaxies detected by ALMA host an AGN, while only 14% of the UVJ active galaxies undetected by ALMA host an AGN.

All galaxies not fulfilling these color criteria are considered as quiescent galaxies and are excluded from our comparison sample (9.3% of the original sample). The ALMA detected galaxies in our survey are massive compared to typical SFGs detected in deep optical and near-IR surveys like CANDELS, in the same redshift range (2 <  z <  4), as shown in Fig. 19.

The high proportion of massive galaxies among the ALMA detected sources suggests that stellar mass can be a strong driver for a source to be detected by ALMA at high redshift (Dunlop et al. 2017). The strong link between detection and stellar mass is related to the underlying relation between stellar mass and star-formation rate of SFGs (e.g., Noeske et al. 2007; Elbaz et al. 2011). Almost one third (7/24) of the galaxies previously cataloged in the field of view of this study with M_⋆ >  10¹¹ and 2 <  z <  3 are also detected with ALMA. The position of our galaxies along the main sequence of star-formation will be studied in a following paper (Franco et al., in prep.).

We observe a lack of detections at redshift z <  2, driven by both a strong positive K-correction favoring higher redshifts and a decrease in the star formation activity at low redshift. Indeed, the specific star-formation rate (sSFR), defined as the ratio of galaxy SFR to stellar mass, drops quickly at lower redshifts (z <  2), whereas this rate increases continuously at greater redshifts (e.g., Schreiber et al. 2015). In addition, very massive galaxies (stellar mass greater than 10¹¹ M_⊙) are relatively rare objects in the smaller co-moving volumes enclosed by our survey at lower redshift. To detect galaxies with these masses, a survey has to be sufficiently large. The covered area is therefore a critical parameter for blind surveys to find massive high redshift galaxies.

In order to estimate the selection bias relative to the position of our galaxies on the main sequence, we show in Fig. 19 the minimum stellar mass as a function of redshift that our survey can detect, for galaxies on the MS of star-formation (green dashed line), and for those with a SFR three times above the MS (green dotted line).

To determine this limit, we calculated the SFR of a given MS galaxy, based on the galaxy stellar mass and redshift as defined in Schreiber et al. (2015). From this SFR and stellar mass, the galaxy SED can also be calculated using the Schreiber et al. (2018) library. We then integrated the flux of this SED around 1.1 mm.

It can be seen that the stellar mass detection limit corresponding to MS galaxies lies at higher stellar mass than all of the galaxies detected by our ALMA survey (as well as all but one of the other star-forming galaxies present in the same region). This means that our survey is unable to detect star-forming galaxies below the main sequence. We can quantify the offset of a galaxy from the main sequence, the so-called “starburstiness” (Elbaz et al. 2011), by the ratio SFR/SFR_MS, where SFR_MS is the average SFR of “main sequence” galaxies computed from Schreiber et al. (2015). We also indicated our detection limit for galaxies with SFR/SFR_MS = 3. In this case, 7/17 galaxies shown lie above the limit. To have been detected, these galaxies must therefore have SFRs at least larger than the SFR_MS, the other ten galaxies must have a SFR at least three times above the MS. This highlights that our survey is biased toward galaxies with high SFRs.

6.3. AGN

In this Section, we discuss the presence of AGN within the 20 most robust ALMA detections, in other words, rejecting the three spurious detections with no IRAC counterpart (AGS14, AGS16 and AGS19 marked with a star in Table 3) but including three of the supplementary sources (AGS21, AGS22 and AGS23). We found an X-ray counterpart for 65% of them (13/20) in the 7 Msec X-ray survey of GOODS-South with Chandra (Luo et al. 2017). Most of these galaxies were classified as AGN in the catalog of Luo et al. (2017) that identifies as AGN all galaxies with an intrinsic 0.5–7.0 keV luminosity higher than L_X, int = 3 × 10⁴² erg s⁻¹, among other criteria. However, as our ALMA galaxies are biased toward highly star-forming galaxies, we decided to increase the minimum X-ray luminosity to a three times stronger X-ray luminosity threshold to avoid any contamination by star-formation. We also considered as AGN the galaxies exhibiting a hard X-ray spectrum. We therefore adopted here the following criteria to identify AGN: either (i)L_X, int >  10⁴³ erg s⁻¹ (luminous X-ray sources) or (ii)Γ <  1.0 (hard X-ray sources).

As the redshifts adopted by Luo et al. (2017) are not always the same as ours, when necessary we scaled the X-luminosities to our redshifts using Eq. (1) from Alexander et al. (2003), and assuming a photon index of Γ = 2.

Using these conservative criteria, we found that eight ALMA galaxies host an X-ray AGN (marked with a yellow star in Fig. 19). In order to compare the AGN fraction among ALMA detections with galaxies undetected by ALMA with similar masses and redshifts, we restricted our comparison to galaxies with M_⋆ >  3 × 10¹⁰ M_⊙ and 1.8 <  z <  4.5 (rectangle in black dotted lines in Fig. 19). In this area encompassing 16 ALMA detections, we found that 50% of the ALMA sources host an AGN (8/16) as compared to only 14% (23/160) of the star-forming galaxies undetected by ALMA located in this same area (selected using the UVJ criteria recalled in Eq. (10) in the ZFOURGE catalog).

The presence of a high percentage of AGN among the galaxies detected by ALMA may reflect the fact that the ALMA sources are experiencing a starburst (well above the MS marked with a green dashed line in Fig. 19), possibly triggered by a merger that may dramatically reduce the angular momentum of the gas and drive it towards the center of the galaxies (e.g., Rovilos et al. 2012; Gatti et al. 2015; Lamastra et al. 2013) or violent disk instabilities (Bournaud et al. 2012). In addition, the high AGN fraction may be driven by the link between the presence of an AGN and the compactness of their host galaxy. Elbaz et al. (2018), Chang et al. (2017), and Ueda et al. (2018) suggest that the proportion of galaxies hosting an AGN increases with IR luminosity surface density. As discussed in Sect. 5.1, the size, and therefore the compactness of a galaxy, increases the likelihood of an ALMA detection at our angular resolution. Alternatively, ALMA might preferentially detect galaxies with a high gas, hence also dust, content, more prone to efficiently fuel the central black hole and trigger an AGN.

This fraction of galaxies with a high X-ray luminosity (L_X, int >  10⁴³ erg s⁻¹) seems to be significantly higher than that found in some other ALMA surveys, in particular in Dunlop et al. (2017; 2/16) or Ueda et al. (2018; 4/12).

7. HST-dark galaxies

Some galaxies without H-band HST-WFC3 (1.6 μm) counterparts have been discovered. We discuss below the possibility that these detections may be real HST-dark galaxies. Some ALMA detections previously attributed to an HST counterpart seem in fact to be either more distant galaxies, extremely close on the line of sight to another galaxy, hidden by a foreground galaxy, or too faint at optical rest-frame wavelengths to be detected by HST.

It is already known that some of the most luminous millimeter or submillimeter galaxies can be completely missed at optical wavelengths (Wang et al. 2016), even in the deepest optical surveys, due to dust extinction. Some of these galaxies can also be undetected in the NIR (Wang et al. 2009).

Among the sources that do not have detections in the H-band of HST-WFC3, we distinguish the sources not detected by HST but detected by other instruments (we will discuss the importance of the IRAC filters), and sources undetected by HST and all of the other available instruments in the GOODS-South field (described in the Sect. 2.4).

Of the 20 galaxies detected in our main catalog, seven (35%) do not present an obvious HST counterpart. This number is slightly higher than the expected number of spurious sources (4 ± 2), predicted by the statistical analysis of our survey. To be more accurate, for three of these seven galaxies (AGS4, AGS15 and AGS17), an HST galaxy is in fact relatively close in the line of sight, but strong evidence, presented below, suggests that the HST galaxies are not the counterpart of the ALMA detections, and without the resolution of ALMA we would falsely associate the counterpart. For the four other ALMA detections without HST-WFC3 counterparts within a radius of 0″60, one of them (AGS11) has also been detected at other wavelengths. In this section, we will discuss four particularly interesting cases of HST-dark galaxies (AGS4, AGS11, AGS15 and AGS17), and discuss our reasons for classifying the other three as spurious sources.

Our four HST-dark galaxies (AGS4, AGS11, AGS15 and AGS17) have at least one feature in common, the presence of an IRAC detection and the fact that this IRAC detection is closer on the sky than the unrelated HST detection (see Table 6). The IRAC detections come from the Ashby et al. (2015) catalog, except for AGS15 where the position comes from the ZFOURGE catalog, using the Labbé et al. (2015) survey. The offset between the IRAC and HST sources might suggest that they are different sources. Figure 20 shows the IRAC contours at 3.6 μm centered on the ALMA detection, superimposed over the HST H-band image. The presence of IRAC detections at these distances from the ALMA galaxies is a very strong driver for the identification of sources. The probability of random IRAC association is between one and two orders of magnitude less likely than random HST association for this range of distances, as shown in Fig. 7 and Table 6. The selection of ALMA candidates from galaxies detected in IRAC channels 1 and 2 but missed by HST-WFC3 at 1.6 μm has already been experimented successfully by Wang et al. (in prep.), and seems to be a good indicator to detect HST-dark ALMA galaxies.

Fig. 20. IRAC 3.6 μm (red contours, 3 μJy–30 μJy in steps of 3.0 μJy) and ALMA 1.1 mm (white contours, 4, 4.5 then 5 to 10-σ in steps of 1-σ) overlaid on 8″3 × 8″3 HST H-band images. The position of the previously associated HST counterpart is shown by a cyan circle.

As each of our HST-dark galaxies have different features, we will discuss each galaxy individually.

7.1. AGS4

AGS4 is a close neighbor of ID_CANDELS 8923. AGS4 is the fourth brightest detection in our survey with an S/N greater than 9. The center of the ALMA detection is located at only 0″38 from ID_CANDELS 8923, its closest neighboring galaxy. Before astrometric correction, this distance was only 0″21. This is therefore an example where the astrometric correction moves the ALMA galaxy away from the supposed counterpart. In Fig. 21, we can clearly see that the ALMA emission is offset from the observed H-band galaxy shown by the white arrow in Fig. 21. This offset could be explained physically, for example, as a region extremely obscured by dust, within the same galaxy, greatly extinguishing the optical rest-frame emission that is revealed by ALMA. However, for AGS4, a series of clues suggest another explanation for this offset.

Fig. 21. Postage stamps of 10 × 10 arcsec from HST-WFC3 at 0.606 μm to VLA at 5 cm, for the four optically-dark galaxies discussed in Sect. 7. For the two ALMA images at 1.1 mm, those marked by 1correspond to the non-tapered images, those marked by 2correspond to the 0″60-mosaic images. The Ks-band thumbnail comes from the super-deep detection image described in Sect. 2.4.1. All images are centered on the ALMA detection. We indicate with white arrows the position of the previously associated HST counterpart.

The first piece of evidence is the comparison between the IR SED at the position of the ALMA detection (the SEDs of all of the galaxies detected in this paper will be presented in a future publication, Franco et al., in prep.) and the redshift of the optical galaxy. The redshift of the optical galaxy is z = 0.241, whereas the far IR SED peaks around 350 μm (see Fig. 22). If AGS4 was a dusty star-forming region on the outskirts of 8923, this infrared SED would suggest an abnormally cold dust temperature. It is therefore more probable that AGS4 is not part of ID_CANDELS 8923, and is a dusty distant galaxy. The fuzzy emission in the H-band HST image, exactly centered at the position of the ALMA detection (see Fig. 21) has not led to any detection in the CANDELS catalog. In the V-band HST images, only ID_CANDELS 8923 is present, seen to the South-East of the position of the ALMA detection (indicated by a white cross). No emission is visible at the exact position of the ALMA detection. In the z-band, a barely visible detection appears extremely close to the center of the image.

Fig. 22. Spectral energy distributions (SED) of the two optically-dark galaxies AGS4 and AGS17. The flux densities from 100 to 500 μm are from GOODS-HerschelElbaz et al. (2011). AGS4 (top) (z = 4.32, see Sect. 7) is fitted with the model of Schreiber et al. (2018). The SED of AGS17 (bottom panel), which has no known redshift, is simply presented with an interpolation between the observed flux densities to illustrate that it peaks around 400 μm. This peak is inconsistent with the redshifts of the two optical sources with IDCANDELS 4414 (z = 1.85) and 4436 (z = 0.92)

The second clue is the detection of a galaxy with redshift z = 3.76 in the FourStar galaxy evolution survey, 0″16 from the ALMA detection. This redshift is much more consistent with the peak of the IR SED. The ZFOURGE survey is efficient at detecting galaxies with redshifts between 1 and 4 by using a K_s-band detection image (instead of H-band as used for the CANDELS survey), and also due to the high spectral resolution (λ/Δλ ≈ 10) of the medium-bandwidth filters which provide fine sampling of the Balmer/4000 Å spectral break at these redshifts (Tomczak et al. 2016). Furthermore, the stellar mass derived in the ZFOURGE catalog (10^10.50 M_⊙ compared with 10^7.64 M_⊙ in the CANDELS catalog) is more consistent with the expected mass of galaxies detected by ALMA. Indeed as shown in this paper, and as already shown by Dunlop et al. (2017), ALMA tends to reveal the most massive dusty galaxies.

The third piece of evidence is the presence, in the Spitzer-CANDELS catalog (Ashby et al. 2015), of a galaxy detected with the IRAC filers only 0″1 from our ALMA detection. This IRAC galaxy has a magnitude of 22.51 at 3.6 μm, measured within an aperture of 2″4 radius.

We also note that Rujopakarn et al. (2016) detect a radio galaxy at S/N ≈ 17 only 55 mas from the center of the ALMA detection shown in Fig. 21 (the positional accuracy of this VLA image is 40 mas). Additionally, AGS4 is detected in two of the three Chandra bands: 0.5–7.0 keV (full band; FB) and 0.5–2.0 keV (soft band; SB), but not at 2–7 keV (hard band; HB) from the 7 Ms Chandra observations of the GOODS-South field. The integrated X-ray flux is only 6.86 × 10⁴⁰ erg s⁻¹, but this galaxy is classified as an AGN in the 7 Ms catalog.

The detection of a local galaxy at this position has been largely documented (e.g., Hsu et al. 2014; Skelton et al. 2014; Santini et al. 2015). In contrast, some studies present the galaxy located at this location as a distant galaxy. Cardamone et al. (2010) take advantage of the 18-medium-band photometry from the Subaru telescope and the photometric redshift code EAzY (Brammer et al. 2008) to derive a redshift z = 3.60. Wuyts et al. (2008) find a redshift of z = 3.52 also using EAzY. We can also add a redshift determination by Rafferty et al. (2011) using the Zurich extragalactic bayesian redshift analyzer (ZEBRA; Feldmann et al. 2006), at z = 2.92. These determinations of high redshift by independent studies support the existence of a distant galaxy at this position.

Although close, the two sources (ID_CANDELS 8923 and 8923b) were successfully de-blended using two light-profile models, determined by fitting the HST H-band image with Galfit. The two sources were then fit simultaneously using these two models on all of the available images, fixing the profile to that observed in the H band. The SEDs of these two galaxies are shown in Fig. 23, in blue for the HST galaxy and in orange for the ALMA galaxy, together with the photometric redshift probability distribution for AGS4. The redshifts were estimated using EAzY. For the blue HST galaxy we found , in good agreement with that found by Skelton et al. (2014). On the other hand, the redshift found for AGS4 is slightly higher than in ZFOURGE, with . However, we can also see a secondary peak in the redshift probability distribution, at the position of the ZFOURGE redshift. As the Balmer break is well established in the K-band, we consider that the redshift determination () is robust and we adopt this redshift for AGS4. The stellar mass of the ALMA galaxy was then computed with FAST (Kriek et al. 2009), and we found 10^{11.45 ± 0.2} M_⊙ (probably slightly overestimated due to the presence of an AGN, suggested by a flux excess in the IRAC bands). The IR SED of this galaxy is shown in Fig. 22. For the first time, thanks to ALMA, we can argue that there exists, at this position, not one but two galaxies, close to each other on the line of sight.

Fig. 23. Spectral energy distributions of AGS4 and IDCANDELS 8923. Aperture photometry allows the separation between the local galaxy detected by HST (blue, and indicated by a white arrow in Fig. 21, IDCANDELS 8923) and the distant galaxy detected by ALMA (orange). The top panel shows the photometric redshift probability distribution of AGS4. As the Balmer break is well established in the K-band, we consider that this redshift determination is robust, and we adopt the derived redshift and stellar mass (1011.45 ± 0.2 M⊙) values for AGS4.

7.2. AGS11

AGS11 is detected at 1.1 mm with a flux of 1.4 mJy (S/N ∼ 8) without any counterpart in the deep HST image. However, the galaxy is also detected by IRAC, confirming the existence of a galaxy at this position. A galaxy was recently found, for the first time, in the ZFOURGE catalog at 0″18 from the center of the ALMA position. This galaxy was not detected directly in the Magellan image but in a super-deep combined K_s-band image at 4.5σ. From this position, the flux in the IRAC-bands have been extracted with S/Ns of 26, 34, 8 and 8 at 3.6 μm, 4.5 μm, 5.8 μm and 8.0 μm respectively.

This HST-dark galaxy falls in a projected overdensity on the sky, consisting of sources in the redshift range 3.42 ≤ z ≤ 3.56 and brighter than K_s <  24.9 (Forrest et al. 2017). This density has been computed by Forrest et al. (2017) using the 7th nearest-neighbor technique (Papovich et al. 2010). This overdensity, centered at RA = 53.08°, Dec = −27.85°, extends beyond approximately 1.8 Mpc.

The redshift derived in the ZFOURGE catalog is z = 4.82, making it the farthest galaxy detected in this blind survey. However, we remain cautious regarding this redshift, as this entry has been flagged in the ZFOURGE catalog (use = 0) due to the S/N of this galaxy (4.7) being below the limit defining galaxies with good photometry (S/N ≥ 5). This galaxy is the only galaxy in our catalog flagged in the ZFOURGE catalog. For this reason, we represent it with an empty circle Fig. 19. AGS11 has not been detected in the 7 Ms Chandra survey.

The stellar mass, derived in the ZFOURGE catalog, 3.55 × 10¹⁰ M_⊙, is consistent with the masses of all of the other ALMA galaxies found in this survey. What is particularly interesting in the multiwavelength images of this galaxy is that AGS11 is detectable only by ALMA and in the IRAC-bands (in non-stacked images). Outside of these wavelengths, no emission is detectable.

7.3. AGS15

AGS15 is at a distance of 0″59 from its possible HST counterpart (ID_CANDELS 3818) after astrometric correction, corresponding to a physical distance of 4.33 kpc. This is the largest HST-ALMA offset in our entire catalog. The IRAC position, in contrast, matches much more closely with the ALMA position, with an offset of only 0″14. The stellar mass of the optical galaxy, given by the ZFOURGE catalog (7.24 × 10⁹ M_⋆) would have made AGS15 a galaxy lying far from the median stellar mass (1.1 × 10¹¹ M_⊙) of our survey. The redshift of ID_CANDELS 3818 (z = 3.46) is nevertheless consistent with the other redshifts found in this study.

7.4. AGS17

AGS17 is a close neighbor (0″27) of ID_CANDELS 4414 (z = 1.85). AGS17 is one of the three galaxies detected by Hodge et al. (2013) at 870 μm in the ALMA field of view (along with AGS8 and AGS15 previously discussed). The counterpart of AGS17 was attributed to ID_CANDELS 4414 by Wiklind et al. (2014) with an offset between the ALMA detection and the corresponding F160W object of 0″32. Again, there are indications that the identification may be false: the peak of the IR SED is ∼400 μm (see Fig. 22), suggesting a more distant galaxy. To be detected with the flux densities reported in Table 4, a galaxy at z = 1.85 would have an extraordinarily high star formation rate (∼820 ± 240 M_⊙ yr⁻¹), using the IR SEDs of Schreiber et al. (2018). If truly associated with the CANDELS counterpart, this galaxy would be an extreme starburst with an SFR 59 ± 17 times greater than the SFR_MS. Galaxies with these properties cannot be ruled out, as galaxies with much higher star formation rates (and offsets from the main sequence) have already been observed (e.g., Pope et al. 2005; Fu et al. 2013). However, such objects are relatively rare. In addition, the stellar mass of ID_CANDELS 4414 (10¹⁰M_⊙) is inconsistent with the trend of the other detections (more than one order of magnitude below the median stellar mass of our catalog).

Another galaxy (ID_CANDELS 4436) is relatively close (0″57) to the ALMA detection. The position of the ALMA detection, which is between ID_CANDELS 4414 and ID_CANDELS 4436, could be the signature of a major merger occurring between these two galaxies. The emission observed by ALMA could result in this case from the heating of the dust caused by the interaction of these two galaxies, but the redshift determination of by Le PHARE (Arnouts et al. 1999; Ilbert et al. 2006) dismisses this hypothesis.

After subtraction of the 2 galaxies close on the line of sight (ID_CANDELS 4414 and 4436) in the HAWK-I image, a diffuse source is revealed (half-light radius = 1″55 ± 0″12, sersic index = 1.0). Lower resolution ALMA observations would be needed in order to correctly measure the total submm flux of this extended source.

We also note the position of the IRAC source, located only 0″06 from the ALMA detection.

7.5. Discussion

Of the total 23 detections in this survey, seven do not show an HST H-band counterpart. This lack of counterpart could arise from an occultation of the optical counterpart by a foreground galaxy, faint emission at optical wavelengths, or a spurious ALMA detection.

For the four galaxies previously discussed (AGS4, AGS11, AGS15 and AGS17), we observe a signal with IRAC at 3.6 μm and 4.5 μm, despite the limiting sensitivity of IRAC (26 AB mag at 3 σ for both 3.6 and 4.5 μm; Ashby et al. 2015) being lower than HST-WFC3 (28.16 AB mag at 5σ for F160W; Guo et al. 2013) in the respective images. Furthermore, two of the galaxies (AGS15 and AGS17) have already been detected at submillimeter wavelengths (870 μm) by Hodge et al. (2013).

The other three galaxies (AGS14, AGS16 and AGS19) are not detected at any other wavelength hence there is a high probability that they are spurious. This number is in good agreement with the expected number of spurious sources for our sample (4 ± 2).

Figure 5 gives us a glimpse into how sources can be falsely associated with an HST galaxy. When the offset correction is applied, the three galaxies shown with magenta lines move further away from the center position (Δδ = 0, Δα = 0), rather than closer to it. Another source also appears to show this behavior: AGS20, seen in the lower left quadrant of Fig. 5. The ALMA detection of AGS20 seems to be clearly offset from an HST galaxy, similar to AGS4. To ensure that there is not a more distant counterpart for AGS20 obscured by the HST source, we performed the same analysis as described in Sect. 7 and illustrated in Fig. 23. The result of the decomposition suggests that the ALMA and HST sources are either two components of the same galaxy or two galaxies merging at this position. A spectroscopic analysis of AGS20 would allow for a distinction between these two possibilities.

The IRAC detections seem to be particularly useful to confirm the existence of a source. In the main catalog, except for the three galaxies that we consider as spurious, all others are also detected in the IRAC filters.

In conclusion, we have detected 20% HST-dark galaxies (4 out of 20 robust detections) with a counterpart confirmed at least by IRAC. This proportion may depend in a manner that we cannot address here on the depth of the optical and millimeter images. Knowing that these HST-dark galaxies are dust, hence metal, rich they are likely progenitors of the most massive galaxies seen at z = 0, hence potentially hosted by massive groups or clusters of galaxies. Two of these HST-dark galaxies have a tentative redshift of z  =  4.82 and z = 3.76, we therefore expect these galaxies to be located on average within z ∼ 4–5. These two galaxies are already massive (10^10.55 M_⊙ and 10^10.50 M_⊙ respectively), suggesting that this population of galaxies is particularly interesting for understanding massive galaxy formation during the first billion years after the Big Bang. Spectroscopy with the JWST NIRSpec instrument will permit very sensitive spectroscopic detection of Hα emission at z <  6.6, and hence an important new tool to measure redshifts of these HST-dark galaxies. GOODS-South will undoubtedly be a venue for extensive JWST spectroscopy, including Guaranteed Time Observations. Spectral scan observations with ALMA can also be a powerful tool to determine the distances, and hence physical properties, of this intriguing population of HST-dark galaxies.

8. Summary and conclusions

The GOODS-ALMA survey covers an area of 69 arcmin² matching the deepest HST-WFC3 coverage of the GOODS-South field at 1.1 mm and at a native resolution of ∼0″24. We used a 0″60 tapered mosaic due to the large number of independent beams at the native resolution. A comparison of the HST source positions with existing catalogs such as Pan-STARRS allowed us to correct the HST astrometry of the GOODS-South field from both a global and local offset (equivalent to a distortion map, see also Dickinson et al., in prep.). We found a median offset between the HST and ALMA images of −96 ± 113 mas in right ascension, α, and 261 ± 125 mas in declination, δ. The main conclusions from our study are listed below.

1. 20 galaxies brighter than 0.7 mJy at 1.1 mm. We detect in total 20 sources above a detection threshold that guarantees an 80% purity (less than 20% chance to be spurious). Among these 20 galaxies (with an S/N >  4.8), we expect 4 ± 2 spurious galaxies from the analysis of the inverted map and we identify 3 probably spurious detections with no HST nor Spitzer-IRAC counterpart, consistent with the expected number of spurious galaxies. An additional three sources with HST counterparts are detected either at high significance in the higher resolution map, or with different detection-algorithm parameters ensuring a purity greater than 80%. Hence we identify in total 20 robust detections.

2. Pushing further in redshift the blind detection of massive galaxies with ALMA. The sources exhibit flux densities ranging from 0.6 to 2 mJy, have a median redshift (and rms) of z = 2.92 ± 0.20 and stellar mass of M_⋆ = (1.1 ± 0.4)×10¹¹ M_⊙. By comparison with deeper but smaller ALMA extragalactic surveys (Aravena et al. 2016; Dunlop et al. 2017; González-López et al. 2017; Ueda et al. 2018), our redshift distribution is shifted to higher values even though our survey is shallower. This is due to the low surface density of massive, metal hence dust-rich, galaxies at high redshifts. The size of the ALMA survey is therefore a key parameter to detect high redshift galaxies.

3. 20% HST-dark galaxies. The detection criteria of this main catalog allowed us to identify sources with no HST counterparts. Out of the 20 galaxies listed above, and excluding the three candidate spurious detections, we identified four optically-dark or HST-dark galaxies with the request of 80% purity and with a Spitzer-IRAC counterpart at 3.6 μm and 4.5 μm, confirming the existence of a galaxy at the position of the ALMA detection. It is not the first time that such HST-dark sources have been found using e.g., infrared color selections (H-dropouts, see e.g., Wang et al. 2016), but their identification in an unbiased survey at the depth of ALMA in the millimeter range allows us to determine that 20% of the ALMA sources detected at 1.1mm above ∼0.7 mJy are HST-dark (4/20 sources in the main catalog). Two of these sources are detected in the near-infrared in the ZFOURGE catalog, with a photometric redshift of z_phot = 4.32 (derived in this study; AGS4, also detected in the radio with VLA) and 4.82 (AGS11). The other two sources (AGS15 & AGS17) were detected with the LABOCA ECDFS Submillimeter Survey (LESS) at 870 μm and with ALMA after a follow-up at the same wavelength, confirming that they were not the result of source blending (Hodge et al. 2013).

4. Exceptionally high AGN fraction. We found a high proportion of AGNs in our ALMA 1.1 mm sample with 40% (8 out of 20 robust detections) detected in the 7 Msec Chandra X-ray survey of GOODS-South in the 0.5–7.0 keV band with a X-luminosity greater than 10⁴³ erg s⁻¹. Limiting our analysis to the ALMA sources with a redshift and stellar mass determination, we found that 50% of the ALMA sources located in a well-defined stellar mass (M_⋆ >  3 × 10¹⁰ M_⊙) – redshift (z ∼ 1.8–4.5) range host an AGN as compared to only 14% for the galaxies located within the same zone but undetected by ALMA. This excess AGN contribution may be due to the fact that the ALMA galaxies are preferentially in a starburst mode due to our detection limit – hence possibly experiencing a merger – or/and that the high-resolution of ALMA favors unresolved, hence compact, sources knowing that the mechanism that leads to such compact star-formation may also trigger an AGN.

5. Alleviating the degeneracy of the bright end of the ALMA counts. The differential and cumulative number counts of our 20 primary detections allowed us to partly alleviate the degeneracy observed above 1 mJy beam⁻¹ in previous (sub)millimeter studies. We show that ∼15% of the extragalactic background light is resolved into individual sources at 0.75 mJy. By extrapolation, ∼50% of the EBL is resolved at 0.1 mJy.

Acknowledgments

We wish to thank Ivo Labbé for sharing with us the IRAC/GREATS images and Bahram Mobasher for sharing the spectroscopic redshift of AGS4 with us. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2015.1.00543.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. RD and NN gratefully acknowledge the support provided by the BASAL Center for Astrophysics and Associated Technologies (CATA) through grant PFB-06 Etapa II. Support for BM was provided by the DFG priority program 1573 “The physics of the interstellar medium”. DMA thanks the Science and Technology Facilities Council (STFC) for support from grant ST/L00075X/1. WR is supported by JSPS KAKENHI Grant Number JP15K17604, Thailand Research Fund/Office of the Higher Education Commission Grant Number MRG6080294, and Chulalongkorn University’s CUniverse. NN acknowledges funding from Fondecyt 1171506.

References

All Tables

All Figures

Fig. 1. ALMA 1.1 mm image tapered at 0″60. The white circles have a diameter of 4 arcseconds and indicate the positions of the galaxies listed in Table 3. Black contours show the different slices (labeled A–F) used to compose the homogeneous 1.1 mm coverage, with a median rms-noise of 0.18 mJy per beam. Blue lines show the limits of the HST/ACS field and green lines indicate the HST-WFC3 deep region. The cyan contour represents the limit of the Dunlop et al. (2017) survey covering all the Hubble Ultra Deep Field region. All of the ALMA-survey field is encompassed by the Chandra deep field-South.

In the text

	Fig. 2. uv-coverage of one of the 846 ALMA pointings constituting this survey. This uv-coverage allows us to perform the source detection in the dirty map.
In the text

Fig. 3. Histogram of pixels of the S/N map, where pixels with noise >0.3 mJy beam−1 have been removed. The red dashed line is the best Gaussian fit. The green dashed line is indicative and shows where the pixel brightness distribution moves away from the Gaussian fit. This is also the 4.2σ level corresponding to a peak flux of 0.76 mJy for a typical noise per beam of 0.18 mJy. The solid black line corresponds to our peak threshold of 4.8σ (0.86 mJy).

In the text

Fig. 4. Cumulative number of positive (red histogram) and negative (blue histogram) detections as a function of the σp (at a fixed σf, left panel) and σf (at a fixed σp, right panel) in units of σ. Solid black line represents the purity criterion pc define by Eq. (1), green dashed-line represents the percentage of positive detection with HST-WFC3 counterpart pHST and magenta dashed-line represents the percentage of positive detection with ZFOURGE counterpart pZFOURGE. Gray dashed-lines show the thresholds σp = 4.8σ and σf = 2.7σ and the 80% purity limit.

In the text

Fig. 5. Positional offset (RAHST–RAALMA, DecHST–DecALMA) between HST and ALMA before (red crosses) and after (blue crosses) the correction of both a global systematic offset and a local offset. The black dashed circle corresponds to the cross-matching limit radius of 0″6. The gray dashed circles show a positional offset of 0″2 and 0″4 respectively. The magenta lines indicate the HST galaxies previously falsely associated with ALMA detections.

In the text

Fig. 6. Postage stamps of 1.8 × 1.8 arcsec. ALMA contours (4, 4.5 then 5–10-σ with a step of 1-σ) at 1.1 mm (white lines) are overlaid on F160W HST/WFC3 images. The images are centered on the ALMA detections. The shape of the synthesized beam is given in the bottom left corner. Astrometry corrections described in Sect. 4.3 have been applied to the HST images. In some cases (AGS1, AGS3, AGS6, AGS13, AGS21 for example), the position of the dust radiation matches that of the stellar emission; in other cases, (AGS4, AGS17 for example), a displacement appears between both two wavelengths. Finally, in some cases (AGS11, AGS14, AGS16, and AGS19) there are no optical counterparts. We will discuss the possible explanations for this in Sect. 7.

In the text

Fig. 7. Probability of a randomly selected position in the area defined by this survey to have at least one HST (blue, orange) or IRAC (green) neighbor as a function of distance. We computed this probability by Monte Carlo simulation using the distribution of galaxies listed in the CANDELS/GOODS-South catalog. Due to the presence of the HUDF within the GOODS-South field, we cannot consider that the density of HST galaxies is uniform, and we consider these two fields separately (blue inside and orange outside the HUDF).

In the text

Fig. 8. ALMA 1.1 mm continuum maps for the 23 detections tapered at 0.60 arcsec. Each 3″5 × 3″5 image is centered on the position of the ALMA detection. Cyan double crosses show sources from the GOODS-S CANDELS catalog. White circles show sources from the ZFOURGE catalog. Blue circles show common sources from both optical catalogs (sources with an angular separation lower than 0″4). The shape of the synthesized beam is given in the bottom left corner.

In the text

	Fig. 9. East–west profile of the PSFs corresponding to the six different parallel slices composing the ALMA image in the 0″60-mosaic (see Fig. 1).
In the text

Fig. 10. Comparison between the stacked PSF (black solid line) and the stack of the 23 ALMA-detections (black dashed line) in the 0″60-mosaic. As each slice has a specific PSF, we stack the PSF corresponding to the position of each detection. The fluxes of each detection have been normalized, so that the brightest sources do not skew the results. Fluxes of the PSF and ALMA detections are normalized to 1. Flux profiles are taken across the East–west direction. The result is consistent with unresolved or marginally resolved sources at this resolution. The insert in the top-right corner shows the same procedure for the 15 sources detected in the 0″29-mosaic (see Table 3).

In the text

	Fig. 11. 10″ × 10″ postage stamps, centered on the galaxy detections. Left to right: source in the 0″60-mosaic map, and residuals obtained after PSF, Gaussian and Sersic flux fitting. The residuals are very similar between the three different extraction methods. Only the 6 brightest galaxies are shown.
In the text

	Fig. 12. Size distribution histogram for the 20 robust detections. These sizes are computed by fitting the ALMA detections with a circular Gaussian in the uv-plane using uvmodelfit in CASA. 85% of the sources exhibit a FWHM smaller than 0″25.
In the text

	Fig. 13. Median source detection completeness for simulated point-like and Gaussian galaxies as a function of integrated flux, for different FWHMs (see figure’s legend). The shaded regions correspond to the standard deviation of 100 runs, each containing 300 simulated sources.
In the text

	Fig. 14. Effective area as a function of flux density, where a source with a given flux can be detected with an S/N >  4.8σ. Ninety percent of the survey area reaches a sensitivity of at least 1.06 mJy beam−1.
In the text

	Fig. 15. Flux boosting as a function of measured S/N estimated from simulations. The median of the boosting is shown by a solid red line. The 1σ confidence intervals (dashed red lines) are overplotted. The solid black horizontal line corresponds to Fout = Fin (see text for details). We used the same set of simulations that we used for the completeness analysis.
In the text

Fig. 16. 1.1-mm cumulative (left panel) and differential (right panel) number counts derived using the corrections described in Sect. 5.4, for the sources detected at > 4.8σ in the main catalog. AGS19 is located at a position where the noise is artificially low, and has therefore not been taken into account. Previous (sub)millimeter cumulative number counts are also shown (Lindner et al. 2011; Scott et al. 2012; Karim et al. 2013; Hatsukade et al. 2013, 2016; Ono et al. 2014; Simpson et al. 2015b; Oteo et al. 2016; Carniani et al. 2015; Aravena et al. 2016; Fujimoto et al. 2016; Umehata et al. 2017; Geach et al. 2017; Dunlop et al. 2017). The different fluxes are scaled to 1.1 mm flux densities using S1.1 mm/S1.2 mm = 1.29, S1.1 mm/S1.3 mm = 1.48, S1.1 mm/S870 μm = 0.56. From the Umehata et al. (2017) study, we used only sources which do not have z = 3.09, (which means we are excluding the protocluster members). Results from single-dish surveys are shown with unfilled pentagon markers and are only indicative, they are not considered for model fitting. The gray curve shows the best-fit Schechter function (with 1-σ) uncertainties, the red curve shows the best-fit DPL function (with 1-σ).

In the text

Fig. 17. Resolved 1.1 mm EBL computed from the best-fit Schechter function to the differential number counts. The green, red, purple and brown lines are from the number counts estimated by Fujimoto et al. (2016), Ono et al. (2014), Hatsukade et al. (2016), and Carniani et al. (2015) respectively. The blue line and shaded region show the results from this work and the associated uncertainty. The solid lines represent the model above the detection limits, and the dashed lines show the extrapolation below these limits. The gray shaded region shows the 1.1 mm cosmic infrared background measured by COBE (Fixsen et al. 1998).

In the text

Fig. 18. Redshift distributions (photometric or spectroscopic) for millimeter-selected galaxies. The blue solid line shows the redshift distribution of our ALMA GOODS-South blind survey. The green dashed line shows the Hubble Ultra Deep Field Survey redshift distribution (Dunlop et al. 2017), the black dash-dotted line shows the ASPECS sample (Aravena et al. 2016), the red dotted line shows the ALMA Frontier Fields survey (González-López et al. 2017) and the yellow dotted line shows the ASAGAO survey (Ueda et al. 2018). Short colored lines at the top of the figure indicate the median redshifts for these four studies.

In the text

Fig. 19. Stellar mass versus redshift for the galaxies detected in our ALMA GOODS-South blind survey (red points). For comparison, the distribution of all of the galaxies, listed in the ZFOURGE catalog, in the same field of view is given in blue. Only UVJ active galaxies are shown. The two HST-dark galaxies for which we have redshifts (AGS4 and AGS11) are represented by open circles. The redshift of AGS11 is however uncertain. The green dashed line shows the position that would be occupied by a typical star-forming galaxy – lying on the median of the SFR–M⋆ star-formation main sequence (MS) – that would produce a 1.1 mm flux density equal to our average detection limit of 0.88 mJy (4.8-σ) using the spectral energy distribution (SED) library of Schreiber et al. (2018). The dotted line illustrates the position of galaxies 3 times above the MS using the appropriate SEDs from the same library. Galaxies hosting an AGN that are undetected or detected by ALMA are identified with black dots and yellow stars respectively. Inside the black dashed rectangle, 50% of the galaxies detected by ALMA host an AGN, while only 14% of the UVJ active galaxies undetected by ALMA host an AGN.

In the text

	Fig. 20. IRAC 3.6 μm (red contours, 3 μJy–30 μJy in steps of 3.0 μJy) and ALMA 1.1 mm (white contours, 4, 4.5 then 5 to 10-σ in steps of 1-σ) overlaid on 8″3 × 8″3 HST H-band images. The position of the previously associated HST counterpart is shown by a cyan circle.
In the text

Fig. 21. Postage stamps of 10 × 10 arcsec from HST-WFC3 at 0.606 μm to VLA at 5 cm, for the four optically-dark galaxies discussed in Sect. 7. For the two ALMA images at 1.1 mm, those marked by 1correspond to the non-tapered images, those marked by 2correspond to the 0″60-mosaic images. The Ks-band thumbnail comes from the super-deep detection image described in Sect. 2.4.1. All images are centered on the ALMA detection. We indicate with white arrows the position of the previously associated HST counterpart.

In the text

Fig. 22. Spectral energy distributions (SED) of the two optically-dark galaxies AGS4 and AGS17. The flux densities from 100 to 500 μm are from GOODS-HerschelElbaz et al. (2011). AGS4 (top) (z = 4.32, see Sect. 7) is fitted with the model of Schreiber et al. (2018). The SED of AGS17 (bottom panel), which has no known redshift, is simply presented with an interpolation between the observed flux densities to illustrate that it peaks around 400 μm. This peak is inconsistent with the redshifts of the two optical sources with IDCANDELS 4414 (z = 1.85) and 4436 (z = 0.92)

In the text

Fig. 23. Spectral energy distributions of AGS4 and IDCANDELS 8923. Aperture photometry allows the separation between the local galaxy detected by HST (blue, and indicated by a white arrow in Fig. 21, IDCANDELS 8923) and the distant galaxy detected by ALMA (orange). The top panel shows the photometric redshift probability distribution of AGS4. As the Balmer break is well established in the K-band, we consider that this redshift determination is robust, and we adopt the derived redshift and stellar mass (1011.45 ± 0.2 M⊙) values for AGS4.

In the text