© M. Béthermin et al. 2020

1. Introduction

Understanding the early formation of the first massive galaxies is an important goal of modern astrophysics. At z > 4, most of our constraints come from redshifted ultraviolet (UV) light, which probes the unobscured star formation rate (SFR). Except for a few very bright objects (e.g., Walter et al. 2012; Riechers et al. 2013; Watson et al. 2015; Capak et al. 2015; Strandet et al. 2017; Zavala et al. 2018; Jin et al. 2019; Casey et al. 2019), we have much less information about dust-obscured star formation, that is, the UV light absorbed by dust and re-emitted in the far infrared. To accurately measure the star formation history in the Universe, we need to know both the obscured and unobscured parts (e.g., Madau & Dickinson 2014; Maniyar et al. 2018).

With its unprecedented sensitivity, the Atacama Large Millimeter Array (ALMA) is able to detect both the dust continuum and the brightest far-infrared and submillimeter lines in “normal” galaxies at z >  4. However, this remains a difficult task for blind surveys. For instance, current deep field observations detect only a few continuum sources at z >  4 after tens of hours of observations (e.g., Dunlop et al. 2017; Aravena et al. 2016; Franco et al. 2018; Hatsukade et al. 2018). Targeted observations of known sources from optical and near-infrared spectroscopic surveys are usually more efficient. For instance, Capak et al. (2015) detected four objects at z >  5 using a few hours of observations.

The [CII] fine structure line at 158um is mainly emitted by dense photodissociation regions, which are the outer layers of giant molecular clouds (Hollenbach & Tielens 1999; Stacey et al. 2010; Gullberg et al. 2015), although it can also trace the diffuse (cold and warm) neutral medium (Wolfire et al. 2003), and to a lesser degree the ionized medium (e.g., Cormier et al. 2012). It is one of the brightest galaxy lines across the electromagnetic spectrum. In addition, at z >  4, it is conveniently redshifted to the > 850 μm atmospheric windows. This line has a variety of different scientific applications since it can be used to probe the interstellar medium (e.g., Zanella et al. 2018), the SFR (e.g., De Looze et al. 2014; Carniani et al. 2018a), the gas dynamics (e.g., De Breuck et al. 2014; Jones et al. 2020), or outflows (e.g., Maiolino et al. 2012; Gallerani et al. 2018; Ginolfi et al. 2020a). It has now been detected in ∼35 galaxies at z >  4, but most of them are magnified by lensing or starbursts and only one third of them are normal star-forming systems (see compilation in Lagache et al. 2018).

Over the past several years, numerous theoretical studies have focused on the exact contribution of the various gas phases (e.g., Olsen et al. 2017; Pallottini et al. 2019) and the effects of metallicity (Vallini et al. 2015; Lagache et al. 2018), gas dynamics (Kohandel et al. 2019), and star-formation feedback (Katz et al. 2017; Vallini et al. 2017; Ferrara et al. 2019) on the [CII] emission, which nowadays is the most studied long-wavelength line at z >  4.

The rest-frame ∼160 μm dust continuum and the [CII] line can be observed simultaneously by ALMA and are the easiest and the most promising features to help gain an understanding of obscured star formation at z >  4. The ALMA Large Program to INvestigate [CII] at Early times (ALPINE) aims to build the first large sample with a coherent selection process at z >  4, increasing the size of the pioneering Capak et al. (2015) sample by an order of magnitude. Le Fèvre et al. (2020) describe the goals of the survey and Faisst et al. (2020) present the sample selection and the properties of galaxies in the sample, which were measured from ancillary data. In this paper, we present the task of processing the ALPINE data from the raw data to the catalogs and the immediate scientific results such as the basic dust and [CII] properties of the ALPINE targets together with the number counts and redshift distribution of the serendipitous continuum detections.

In Sect. 2, we describe the ALPINE data processing and the main products (maps and cubes). In Sect. 3, we explain how we built the continuum source catalog and characterized the performance of our method (purity, completeness, and photometric accuracy). In Sect. 4, we derive a reliable conversion factor from the 158 μm rest-frame dust continuum to the total infrared luminosity (L_IR, 8−1000 μm) and the infrared SFR (SFR_IR) using the stacking of ancillary single-dish data at the position of photometric samples similar to ALPINE. In Sect. 5, we discuss the continuum properties of ALPINE detections and the statistical properties of nontarget sources found in the fields (redshift distribution, number counts). In Sect. 6, we describe the procedure used to generate and validate the [CII] spectra and catalog. In Sect. 7, we discuss the properties of the [CII] detections (luminosity, width, velocity offset) and we briefly discuss the correlation between SFR and [CII] luminosity. In this paper, we assume Chabrier (2003) initial mass function (IMF) and a ΛCDM cosmology with Ω_Λ = 0.7, Ω_m = 0.3, and H₀ = 70 km s⁻¹ Mpc⁻¹.

2. Data processing

2.1. Observations

The ALPINE-ALMA large program (2017.1.00428.L, PI: Le Fèvre) targeted 122 individual 4.4 <  z_spec <  5.9 and SFR ≳ 10 M_⊙ yr⁻¹ galaxies with known spectroscopic redshifts from optical ground-based observations. The construction and the physical properties of the sample is described in Le Fèvre et al. (2020) and Faisst et al. (2020), respectively. The ALPINE sample contains sources from both the cosmic evolution survey (COSMOS) field and the Chandra deep field south (CDFS).

In this redshift range, the [CII] line falls in the band 7 of ALMA (275−373 GHz). To avoid an atmospheric absorption feature, no source has been included between z = 4.6 and 5.1. In order to minimize the calibration overheads, we created many groups of two sources with similar redshift, which are observed using the same spectral setting. In our sample, the typical optical line width is σ ∼ 100 km s⁻¹ (or FWHM ∼ 235 km s⁻¹). At the targeted frequency, the coarse resolution (Δν_channel = 31.250 MHz) offered by the Time Division Mode (TDM) is sufficient to resolve our lines (Δv_channel = 25−35 km s⁻¹) and results in a total size of our raw data below 3 TB for the whole sample. The [CII] lines of the targeted sources are covered by two contiguous spectral windows (1.875 GHz each), while we placed two remaining spectral windows in the other side band to optimize the bandwidth and thus the continuum sensitivity. To maximize the integrated flux sensitivity, we requested compact array configurations (C43-1 or C43-2) corresponding to a > 0.7 arcsec resolution to avoid diluting the flux of our sources into several synthesized beams.

We aimed for a 1-σ sensitivity on the integrated [CII] luminosity L_[CII] of 0.4 × 10⁸ L_⊙ assuming a line width of 235 km s⁻¹. As shown in Sect. 7.3, this sensitivity was reached on average by our observations. At higher redshift (lower frequency), we need to reach a lower noise in Jy beam⁻¹ to obtain the same luminosity (∼0.2 mJy beam⁻¹ in 235 km s⁻¹ band at z = 5.8 versus ∼0.3 mJy beam⁻¹ in the same band at z = 4.4). In contrast, at low frequency, the noise is lower because of the higher atmospheric transmission and the lower receiver temperature. The two effects compensate each other and the integration times are similar for our entire redshift range (15−25 min on source). Each scheduling block containing the observations of the calibrators and two sources can be observed using a single 50 min−1 h15 min execution. In total, we had 61 scheduling blocks (SBs) for a total of 69.3 h including overheads.

ALPINE was selected in cycle 5 and most of the observations were completed during this period. Between 2018/05/08 and 2018/07/16, 102 of our sources were observed. Observations had to be stopped from mid-July to mid-August because of exceptional snowstorms. Two additional sources were observed after the snow storms (2018/08/20). After that, the configuration was too extended and the 18 last sources were carried over in cycle 6. They were observed between 2019/01/09 and 2019/01/11.

We realized during the data analysis that four ALPINE sources were observed two times with different names: vuds_cosmos_5100822662 and DEIMOS_COSMOS_514583, vuds_cosmos_5101288969 and DEIMOS_COSMOS_679410, vuds_cosmos_510786441 and DEIMOS_COSMOS_455022, vuds_efdcs_530029038 and CANDELS_GOODSS_15. In the rest of the paper, we use the VIMOS ultra deep survey (VUDS) name of these objects. We thus combined the two ALPINE observations of each of these sources to obtain deeper cubes and maps. Our final sample contains 118 objects.

2.2. Pipeline calibration and data quality

The data were initially calibrated at the observatory using the standard ALMA pipeline of the Common Astronomy Software Applications package (CASA) software (McMullin et al. 2007). We checked the automatically-generated calibration reports and identified a few antennae with suspicious behaviors (e.g., phase drifts in the bandpass calibration, unstable phase or gain solutions, anomalously low gains or high system temperatures), which were not flagged by the pipeline. For example, we had to flag the DV19 antenna for all the cycle 6 observations, for which the bandpass phase solution drifted by ∼180 deg GHz⁻¹ in the XX polarization. For half of the observations, no problems were found and we used directly the data calibrated by the observatory pipeline. Most of the other observations were usually good with only 1 or 2 antennae with possible problems. Four SBs have between 3 and 5 potentially problematic antennae. Considering the very low impact of a single antenna on the final sensitivity, we thus decided to be conservative and fully flag these suspicious antennae and subsequently excised them from our analysis.

While the reduction process was generally smooth, we encountered a couple minor issues. The pipeline sometimes flagged the channels of a spectral window overlapping with the noisy edge channels of another spectral window. It was solved by adding the fracspw = 0.03125 option to the hifa_flagdata task before re-running the pipeline script from the observatory. This option flags the edge channels corresponding to 3.125% of the width of the spectral window, while the default is to flag two channels on each side in TDM mode, that is 4/128 = 0.0315. In theory, this command is equivalent to the default routine. In practice, it is not affected by the subtle bug flagging the channels of the other spectral windows when they overlap, which solves our problem. In a few cases, the pipeline used an inconsistent numbering of the spectral windows and we had to manually correct these problematic SBs.

2.3. Flux calibrators variability and calibration uncertainties

The stability of the flux calibration over our entire survey is particularly important to interpret the sample statistically. We thus checked that the quasars used as secondary flux calibrators were reasonably stable across the ALPINE observations. These secondary calibrators are J1058+0133 and J0854+2006 for the targets in the COSMOS field and J0522−3627 for the ones in the CDFS. We downloaded the data from their flux monitoring by the observatory and calibrated using a well-known primary calibrator¹. In Fig. 1, we present the evolution of their band-7 flux density and the spectral index determined using their measured band-7 and band-3 fluxes.

Fig. 1. Upper panel: 345 GHz flux density of the flux calibrators used by the ALPINE survey (the J1058+0133, J0854+2006, and J0522−3627 quasars) as a function of time. The gray areas indicate when the ALPINE targets were observed (see Sect. 2.1). The only quasar used for CDFS targets (J0522−3627) is plotted with a dashed line, while solid lines are used for the calibrators of COSMOS sources. Lower panel: spectral index versus time. This spectral index is estimated using the band-7 and band-3 flux from the calibrator monitoring performed by the observatory (see Sect. 2.3).

The three quasars are reasonably stable between two successive observations and in particular during the ALPINE observations (gray area in Fig. 1). The standard deviation of the relative difference between two successive data points is only 0.059, 0.060 and 0.031 for J0522−3627, J0854+2006, and J1058+0133, respectively. In the figure, the variability of J0522−3627 could seem larger than J0854+2006. However, the actual relative variations between two successive observations are similar. The larger flux of J0522−3627 highlighting small relative variations in our linear-scale plot and the presence of long-term trends at the scale of several months can give this wrong impression. The maximum relative deviation between two successive visits is 0.20 and happened in J0854+2006 in November 2018, when ALPINE observations were not scheduled. Except this outlier, the maximal variation is 0.13. Usually, the last measurements performed by the quasar monitoring survey are used to determine the flux reference to calibrate a science observation. We can thus expect that the calibration uncertainty coming from the variability of the quasars is usually 6% with 13% outliers.

The frequency reference used for this monitoring is 345 GHz. However, for the highest redshift object of our sample, the spectral setup is centered around 283 GHz. The observatory uses the previously-measured spectral index measured using band-7 and band-3 data to derive the expected flux at the observed frequency. If this index varies too much between two monitorings, it could be a problem. The standard deviation of the spectral index between two successive monitorings is 0.075, 0.069, and 0.050 for J0522−3627, J0854+2006, and J1058+0133, respectively. This corresponds to an uncertainty of 1.5%, 1.4%, and 1.0% on the extrapolation of the flux from 345 GHz to 283 GHz. The largest jump (0.22 in J0522−3627) corresponds to 4.5 %. The typical 1-σ uncertainty of the calibration thus is 7.5% for J0522−3627 and J0854+2006 and 4% for J1058+0133 combining linearly the flux and spectral index uncertainties to be conservative and for the source requiring the most uncertain frequency interpolation. Our calibration uncertainty caused by quasar variability is thus slightly smaller than the typical 10% of uncertainty of interferometric calibrations.

2.4. Data cube imaging and production of [CII] moment-0 maps

The datacube were imaged using the tclean CASA routine using 0.15 arcsec pixels to well sample the synthesized beam (6 pixels per beam major axis in the field with the sharpest synthesized beam). The clean algorithm is run down to a flux threshold of 3σ_noise, where σ_noise is the standard deviation measured in a previous nonprimary-beam-corrected cube after masking the sources. The determination of the final clean threshold is thus the result of an iterative process. The noise converges very quickly with negligible variations between the second and the third iteration. In practice, the exact choice of the clean threshold has a very low impact on the final flux measurements, since our pointings mostly contain one or a few sources, which are rarely bright. In addition, the natural weighting produces sidelobes and high signal-to-noise ratio (S/N) sources can produce nonnegligible artifacts in the dirty maps or unproperly cleaned maps. We checked that the amplitude of the largest sidelobes are below 10% of the peak of the main beam. The sidelobe residuals after cleaning down to 3σ should thus be below 0.3σ.

The standard ALPINE products were produced using a natural weighting of the visibilities. This choice maximizes the point-source sensitivity and produces a larger synthesized beam than other weighting schemes, which limits the flux spreading across several beams for slightly extended sources. These cubes are thus optimized to measure integrated properties of ALPINE targets.

We also produced continuum-free cubes. The continuum was subtracted in the uv-plane using the uvcontsub CASA routine. This routine takes as input a user-provided range of channels containing line emission, and masks them before fitting a flat continuum model (order 0) to the visibilities. To identify the channels to mask, we used the line properties determined using the method presented in Sect. 6. We use several iterations of the cube production and the line extraction to obtain the final version of these products. To avoid any line contamination, we chose to be conservative and excluded all the channels up to 3-σ_v from the central frequency of the best Gaussian fit of the line. When a [CII] spectrum exhibits a non-Gaussian excess in the wings, we masked manually an additional ∼0.1−0.2 GHz to produce conservative continuum-free cubes.

Finally, we generated maps of the [CII] integrated intensity by summing all the channels containing the line emission, that is the moment-0 maps defined as $M(x,\mathit{y})=\sum _{k=1}^{N_{\mathrm{channel}}}{S}_{\nu }(x,\mathit{y},k)\mathrm{\Delta }{\mathit{v}}_{\mathrm{channel}}(k)$, where S_ν(x, y, k) is flux density in the channel k at the position (x, y) and Δv_channel(k) is the velocity width of channel k. The integration windows were manually defined using the first extraction of the spectra as shown in Fig. C.1. Contrary to the continuum subtracted cubes, the integration window is not defined in a conservative way (see Sect. 6.1), but designed to avoid adding noise from channels without signal in the moment-0 maps.

2.5. Continuum imaging

We produced continuum maps using the similar method as for the cubes (same clean routine, pixel sizes, and weighting as in Sect. 2.4), except that the continuum maps were produced using multi-frequency synthesis (MFS, Conway et al. 1990) rather than the channel-by-channel method used for the cubes. The MFS technique exploits the fact that various continuum channels probe various positions in the uv plane to better reconstruct 2-dimensional continuum maps. We excluded the same line-contaminated channels as for the uv-plane continuum subtraction used to produce the cubes. Only the lines of the ALPINE target sources were excluded. Some off-center continuum sources with lines were serendipitously detected in the field. A specific method has been used to measure their continuum flux (see Sect. 3.4).

Some sources could be significantly more extended than the synthesized beam. To detect them, in addition to natural-weighted maps, we also produced lower-resolution uv-tapered maps, which are maps imaged assigning a lower weighting to the visibilities corresponding to small scales. We used a Gaussian 1.5 arcsec-diameter tapering. In Sect. 3.1, we discuss the extraction of the sources using the normal and the tapered maps simultaneously.

2.6. Achieved beam sizes and sensitivities

The achieved synthesized beam size varies with the frequency and the exact array configuration, when each source was observed. The average size of minor axis is 0.85 arcsec (minimum of 0.72 arcsec and maximum of 1.04 arcsec), while the average major axis size is 1.13 arcsec (minimum of 0.9 arcsec and maximum of 1.6 arcsec). Our data follow the requirements on the beam size (> 0.7 arcsec). The mean ratio between major and minor axis is 1.3 and the largest value is 1.8.

The [CII] sensitivity was measured on the moment-0 maps. The mean integrated line flux root mean square (rms) sensitivity is 0.14 Jy km s⁻¹. The mean sensitivity is better in the low-frequency range (283−315 GHz, 5.1 <  z <  5.9) with 0.11 ± 0.04 Jy km s⁻¹ than in the high-frequency range (345−356 GHz, 4.3 <  z <  4.6) with 0.17 ± 0.04 Jy km s⁻¹. A difference of sensitivity between fields observed at similar frequency can also be caused by different widths of the velocity window used to integrate the line fluxes. In Fig. 2 (upper panel), we show the sensitivity versus frequency achieved in each field after renormalizing the effect caused by the different bandwidths used to produce the moment-0 maps of our targets. The mean sensitivities at the low and the high frequency (red squares) follow very well the trend expected from a constant [CII] luminosity. This is not surprising, since the survey was designed to have this property, but it is good to actually achieve it with the real data. However, beyond this very smooth overall trend, there is a large scatter around the mean behavior, since sources were observed under different weather conditions and variable number of good antennae.

Fig. 2. Achieved [CII] (upper panel) and continuum (lower panel) rms sensitivities. The blue dots indicate the values measured in individual fields and the red squares the mean values in the two redshift windows. The red error bars on the plots are the standard deviation in each redshift range. The actual uncertainties on the mean values are indeed $\sqrt{N}$ times smaller (central-limit theorem) and are smaller than the size of the squares. Since the [CII] sensitivities were measured using different bandwidths because of the different line widths, we normalized the measurements to a bandwidth of 235 km s−1 by dividing our raw measurements by $\sqrt{\mathrm{\Delta }\mathit{v}/(235\mathrm{km}{\mathrm{s}}^{-1})}$. The solid black lines indicate the trend of the [CII] flux I[CII] and continuum flux S(1 + z)158 μm versus frequency (and thus redshift) at constant [CII] luminosity L[CII] and fixed infrared luminosity LIR, respectively.

The continuum sensitivity also varies with the frequency. For the sources in the 4.3 <  z <  4.6 range (345−356 GHz), the mean sensitivity is 50 μJy beam⁻¹. We obtained a better sensitivity for 5.1 <  z <  5.9 sources (283−315 GHz) with an average value of 28 μJy beam⁻¹. The slope of the continuum sensitivity versus frequency is steeper than the continuum flux density versus redshift at fixed infrared luminosity L_IR (see the solid black line in Fig. 2 upper panel, computed assuming the Bethermin et al. 2017 spectral energy distribution template as discussed in Sect. 4). This means that our L_[CII]-limited survey is paradoxically able to detect to detect galaxies with lower L_IR at higher redshift.

The performances obtained in each pointing are listed in Table A.1.

3. Continuum catalog

3.1. Source extraction method, detection threshold and purity

To extract the continuum sources, we created signal-to-noise-ratio (S/N) maps. We started from the nonprimary-beam-corrected map, which are maps not corrected for the low gain of the antennae far from the pointing center (normally the same as the phase center if the pointing is correct). These maps have the convenient property to have a similar noise level in the center and on the edge of the antennae field of view. We checked this by comparing the noise in the inner and outer regions of the maps (the border was set at 10 arcsec from the center) and found only a 0.7% higher noise in the center on average. This small excess may be caused by faint undetected sources in the central region, where the primary gain of the antennae is higher. The noise is computed using the standard deviation of the maps after excluding the pixels closer than 1 arcsec to the phase center (possibly contaminated by our ALPINE target) and applying a 3-σ clipping to avoid any noise overestimation due to serendipitous bright continuum sources. The final S/N maps are obtained by dividing the nonprimary-beam-corrected map by the estimated noise. The source are then extracted by searching for local maxima using the find_peak of astropy (Astropy Collaboration 2013). To avoid missing extended sources, we apply the same procedure to the tapered continuum maps (Sect. 2.5) and merge the two extracted catalogs. For the sources present in both catalogs, we use the position measured in the nontapered maps, where the synthesized beam is sharper. Practically, very few sources have a higher S/N in the tapered map due to their much higher noise.

The choice of the S/N threshold is crucial. If it is too low, the sample is contaminated by peaks of noise and the purity is very low. If it is too high, the faint sources are missed. We estimated the purity of the extracted sample as a function of the S/N by comparing the number of detections in the positive and the negative maps. The purity is computed using:

$$\begin{array}{c}\hfill \mathrm{purity}=\frac{N_{\mathrm{pos}}-N_{\mathrm{neg}}}{N_{\mathrm{pos}}}=\frac{N_{\mathrm{real}}}{N_{\mathrm{real}}+N_{\mathrm{spurious}}},\end{array}$$(1)

where N_pos is the number of detections in the positive map, which is also the sum of the N_real real and the N_spurious spurious sources. The average expected number of spurious sources in the positive and negative maps should be the same because the noise in our data is symmetrical. This is why we use the same N_spurious notation for both. N_neg is the number of detections in the negative map. Since we do not expect any real source with a negative flux in our data, this number is equal to the number of spurious sources (N_neg = N_spurious). Of course, this is only true on average and Eq. (1) is only valid when N is large. The purity of the sample extracted from all the pointings as a function of the S/N threshold is presented in Fig. 3 (in red for the full field). The uncertainties are computed assuming Poisson statistics. The 95% purity is reached for a S/N of 5.05 and we decided to cut our catalog at the standard 5σ.

Fig. 3. Purity as a function of the S/N threshold. The results obtained around the center of the pointings (1 arcsec radius) are in blue. The results in the full field are in red. The dotted lines show the S/N at which the 95% is reached.

Out of the 67 sources detected above 5σ, only 11 of them are close enough to the phase center to be potentially associated to an ALPINE target. However, when trying to detect a source close to the center of the field, we explore a much smaller number of synthesized beams (lower risk to detect high-S/N serendipitous sources) and a larger fraction of these beams are expected to contain a real source (higher ratio between real and spurious detections). Therefore, the S/N at which we reach 95% completeness should be lower than in the entire field. We thus estimated the purity versus S/N considering only the central region of each pointing. The distribution of the distance of the detections to the phase center has a bump at small distance with a 1-σ width of 0.4 arcsec. Spatial offsets are discussed in Faisst et al. (2020). We thus decided to use a 1 arcsec radius to define the central region, which should contain 98.7% of the ALMA continuum counterparts of our targets. In this small region, we found no S/N >  5 source and only two S/N >  3 sources in the negative map. To reduce the statistical uncertainties on N_neg, we computed the number of sources in the total survey and rescaled by the ratio between the sum of the areas of the 118 central regions and the total imaged area of the survey. The final result is presented in Fig. 3 (blue curve). We reach a purity of 95% for a S/N = 3.5 cut. With this new threshold, we obtain 23 detections in the central regions, doubling the number of detected target sources.

We call target sample the sources extracted in the 1 arcsec central regions and nontarget sample the objects found outside of this area. The cutout images of these sources are shown in Fig. B.1. The position and the S/N of our target and nontarget detections are provided in Tables B.1 and B.3, respectively.

3.2. Photometry

Many methods can be used to measure the flux of compact sources in interferometric data. These methods have various strengths and weaknesses. We thus decided to derive flux density values using four different map-based methods: peak flux, elliptical Gaussian fitting, aperture photometry, and integration of the signal in the 2-σ contours. The first three methods are standard to analyze interferometric data. These four measurements are made automatically to allow us to perform easily Monte Carlo simulations to validate them. In Sect. 3.5, we check the consistency between these methods.

All our measurements have been performed in the cleaned maps. Given that complex artifacts can appear during the cleaning process, as a test we performed the same measurements in the uncleaned (dirty) maps and found an excellent agreement in all the pointings, which do not contain bright sources producing side lobes.

The most basic method is to measure the peak flux of the source. The uncertainty is derived by dividing the noise measured in the nonprimary-beam-corrected map by the gain of the primary beam at the position of the source. While this method is optimal to measure point-source flux densities, it underestimates the flux of extended sources.

A simple way to measure the flux of compact marginally-resolved sources is fitting a two-dimension elliptical Gaussian. We used the astropy fitting tools (Astropy Collaboration 2013, 2018) and chose a 3 arcsec fitting box. The flux density of the source is just the integral of this Gaussian divided by the integral of the synthesized beam normalized to unity at its peak. The sources for which this method does not perform well are the extended clumpy or nonaxisymetrical sources, which are not well fit by an elliptical Gaussian. The uncertainties can be difficult to compute, since the noise in interferometric maps is correlated at the scale of the synthesized beam. We use the formalism of Condon (1997), who proposed a simplified formalism to propagate the uncertainties.

Aperture photometry, that is the integration of flux in a circular aperture, relies on fewer assumptions than the previous method. We used the routine from the astropy photutils package. The aperture radius needs to be chosen carefully. If it is too small, it will miss extended flux emission from the source. If it is too large, the relative contribution from the noise increases, which makes the measurements uncertain. By comparing the mean flux measured for our sample with different apertures, we showed that for most of them the flux converges for apertures around 3 arcsec diameter. Beyond that, we do not gain flux anymore, but the measurements become noisier. We thus chose this aperture for the ALPINE catalog. We estimated the noise σ_aper using the following formula:

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{aper}}=\frac{{\sigma }_{\mathrm{center}}}{G_{\mathrm{pb}}}\sqrt{\frac{\pi D^2}{4{\mathrm{\Omega }}_{\mathrm{beam}}}},\end{array}$$(2)

where σ_center is the rms of the nonprimary beam corrected map, which is also the rms expected at the center of a given pointing. G_pb is the gain of the primary beam at the position of the source, which is unity at the phase center and decreases when the distance from it increases. σ_aper is thus higher on the edge of the field than in the center. In theory, the gain slightly varies across the aperture, but we checked that using the value at the center of the aperture is a good approximation. D is the diameter of the aperture and Ω_beam is the solid angle of the synthesized beam. The normalization of the noise by the square root of the ratio between the aperture area and Ω_beam is equivalent to rescaling the noise by the square root of the number of independent primary beams in the aperture (N_ind). We checked the validity of this approximation by measuring the aperture flux at random empty positions. N_ind varies from 4.2 to 9.2 in the various pointings with a mean value of 6.7. The flux uncertainties are thus on average 2.6 times higher for the aperture photometry than the peak measurement. This is the main weakness of this method.

Finally, we used another slightly less standard approach in millimeter interferometry, for which we define a S/N-based custom region, from which we integrate the source flux. This method has the advantage to produce smaller integration area for compact unblended sources than the large standard aperture described previously. It is similar to an isophotal magnitude measurement performed in optical astronomy, except that the integration area is defined in S/N instead of surface brightness. It is also better suited for sources with complex shapes. However, it does not deblend the close sources in multi-component systems, and tends to define very large areas encompassing the full blended systems (see Appendix D.2). Practically, we define our integration region as the contiguous area around the source where the S/N map is higher than 2. This value has been chosen after performing tests on a small subset of our sample. For point sources close to thnone S/N threshold, this region is smaller than the synthesized beam and the flux would be underestimated. We thus compute the correction to apply by measuring the synthesized beam map produced by CASA using a region with the exact same shape. Similarly to the aperture method, we compute the flux uncertainties by rescaling the noise by the square root of the number of independent synthesized beams in the region. For simplicity, this method will be called 2-σ clipped photometry in this paper.

The flux densities measured for our target and nontarget detections can be found in Tables B.1 and B.3. Four of our continuum detections required a manual measurements of their flux because they are either multi-component or blended with a close bright neighbor. These peculiar systems are discussed in Appendix D.1.

3.3. Upper limits for nondetected target sources

A large fraction of the ALPINE targets are not detected in continuum (80%), since our survey is able to detect only the most star-forming objects of our sample (see Sect. 5.1). To produce 3-σ upper limits, the easiest widely-used approach is to take 3 times the rms of the noise. Since the target sources are at the phase center, it is just 3σ_center in our case (see the column called “aggressive” upper limits in Table B.2). However, these upper limits are a bit too aggressive. If an intrinsic 2.999σ_center signal is present at the position of the source and if we assume a flat prior on the flux distribution of the sources, there is ∼50% probability that the source is actually brighter than 3σ_center. Therefore, we produced more robust upper limits by summing 3σ_center with the highest flux measured 1 arcsec around the phase center (“normal” upper limits in Table B.2). In the extreme case of a significantly extended source, the source could also be missed because its peak flux is a small fraction of the integrated source flux. We produced “secure” upper limits (Table B.2) by applying the previous process to the tapered maps. We recommend to use these “secure” upper limits, except in the case of point sources for which the “normal” ones are appropriate.

3.4. Line contamination of the continuum of nontarget sources

Our continuum maps were produced excluding the channels contaminated by the [CII] line of the target sources only. The [CII] or another line can contaminate the flux density measurements of nontarget sources if it is outside of the excluded frequency range (Sect. 3.2). To identify these problematic cases, we extracted their spectra and after visual inspection found 9 objects withnontarget a possible line contamination. The nature of these objects will be discussed in Loiacono et al. (2020). We generated new continuum maps, where we masked the line-contaminated channels of the nontarget source instead of the ALPINE target ones. We then remeasured the continuum flux using the same method as previously. Table 1 summarizes the impact of this line decontamination. The relative impact of this correction can vary from a 58% decrease of the flux density (SC_1_DEIMOS_COSMOS_787780) to a nonstatistically-significant increase of the flux (SC_2_DEIMOS_COSMOS_818760). It might be surprising that the line-free flux does not decrease significantly in some sources compared to the initial measurements (or even increase by a fraction of σ in the case of SC_2_DEIMOS_COSMOS_818760), but the contaminating line can sometimes overlap with the [CII]-contaminated channels of the target source, which were masked initially.

3.5. Consistency of the various photometric methods

Since the photometry of each source was determined using different methods, the consistency between these methods can be used as a robustness check (see Fig. 4). The 2D-fit, aperture, and 2σ-clipped measurements are overall in excellent agreement with each other (see the two upper panels of Fig. 4). Even if most of the measurements are compatible at 1σ with each other, there is a small proportional offset of −3.4% and +1.4% between the aperture photometry and 2σ-clipped photometry, respectively, and the 2D-fit measurements. This remains negligible compared to the typical 10% absolute calibration uncertainties of interferometric observations (see Sect. 2.3).

Fig. 4. Comparison between our various photometric methods described in Sect. 3.2 for S/N >  5 sources. The blue dots are our measurements and the red line is the one-to-one relation. Upper left, upper right, and lower left panels: comparison between the 2D-fit flux density (x-axis) and the aperture, 2σ-clipped, and peak flux densities, respectively. Lower right panel: ratio between the peak flux and the 2D-fit flux as a function of the ratio between the source area (convolved by the synthesized beam) and the synthesized beam area. The dashed line indicates the expected trend (see Sect. 3.5).

In order to check how consistent are our measurements, we computed the uncertainty-normalized difference between two measurements $S_{\nu }^{\text{method} \text{A}}$ and $S_{\nu }^{\text{method} \text{B}}$:

$$\begin{array}{c}\hfill \frac{(S_{\nu }^{\mathrm{method}\mathrm{A}}-S_{\nu }^{\mathrm{method}\mathrm{B}})}{\sqrt{{\sigma }_{\mathrm{method}\mathrm{A}}^2+{\sigma }_{\mathrm{method}\mathrm{B}}^2}},\end{array}$$(3)

where σ_method X is the uncertainty derived for the method X. If the two measurements would be performed on independent realizations of the noise, the standard deviation of the normalized difference measured for a large sample should be close to unity. We found 0.40 and 0.66 for the comparison between aperture and clipped photometry, respectively, and 2D-fit measurements. It shows that the three methods are overall consistent at better than 1σ. It is not surprising to find a value below unity, since our methods are using the same realization of the noise. We did not expect to find zero either, since each method tends to weight the noise in the various pixels in a different way.

The peak photometry does not agree as well with the other methods and is on average 19% lower than the 2D-fit flux (Fig. 4, lower left panel). This clearly indicates that our sources cannot be considered as point like and that the peak flux is not a good way to measure their integrated flux. In a Gaussian-profile case without noise, the ratio between the peak flux and the integrated flux directly depends on the source size and the synthesized beam size. If we note Ω_beam the beam area defined as the integral of the synthesized beam and Ω_source the integral of the profile of an extended source after normalizing its peak to unity, the peak flux S_peak is:

$$\begin{array}{c}\hfill S_{\mathrm{peak}}=S_{\mathrm{int}}\frac{{\mathrm{\Omega }}_{\mathrm{beam}}}{{\mathrm{\Omega }}_{\mathrm{source}}},\end{array}$$(4)

where S_int is the integrated flux. The S_peak/S_int ratio should thus be inversely proportional to Ω_source/Ω_beam. In the lower right panel of Fig. 4, we show that this is exactly the trend followed by our measurements.

3.6. Comparison between map-based and uv-plane photometry

In millimeter interferometry, we can also measure the flux of a source directly in the uv-plane. This technique is particularly powerful to deblend multiple sources and when the uv-coverage is limited. To perform the uv-fitting, we used the GILDAS² software package MAPPING, which allows us to fit models directly to the uv visibilities. The use of GILDAS required beforehand to export our CASA measurement sets to uvfits tables and then to uvt tables, the GILDAS visibility table format³. We could successfully model nine continuum targets⁴ detected at ≥5σ and without any bright neighbor, using an elliptical Gaussian model for which the analytical Fourier transform could be fit to the merged visibilities of all channels of the 4 spectral windows and the two polarizations (excluding only channels contaminated by the [CII] emission line). We derived uv-based flux measurements for all these targets, marginally resolved in most cases. In Fig. 5, we show the comparison between the map-based 2D-fit method and the uv-plane approach. All our sources are compatible at 1σ with the one-to-one relation. This shows that measuring the flux in map space is sufficient in our case. However, uv-plane modeling is critical for size measurements and is presented in Fujimoto et al. (2020).

Fig. 5. Comparison between the 2D-fit flux densities derived in map space (Sect. 3.2) and the flux determined fitting an elliptical Gaussian model in the uv plane (Sect. 3.6). The blue dots are our measurements and the red line is the one-to-one relation.

3.7. Monte-Carlo source injections

To interpret the statistical properties of nontarget detections, we need to know the completeness in our various pointings as a function of the source flux density, the source size, and the distance to the phase center. We used Monte-Carlo source injections to estimate it, but also to test the reliability of our flux measurements.

We performed injections of sources using a grid of 4 different intrinsic sizes (FWHM = 0, 0.333, 0.666, and 1 arcsec) and 18 different nonprimary beam corrected flux densities ranging from 0.02 mJy to 1 mJy spaced by 0.1 dex. We injected 10 sources in any given pointing, which is sufficiently small to avoid overlap problems and sufficiently large to be efficient at getting a large number of injected sources in a reasonable computing time. We decided to repeat this task 10 times per set of properties (size and flux) in order to have 100 objects per size and flux. Because of our limited computing resources, we limited our study to Gaussian circular sources and we injected sources directly in the image space. For each realization, we extracted the sources and measured their flux using the same exact method as for the real maps. We consider that a source is recovered if it is found less than 1 arcsec from its injected position. We checked that the number of recovered sources are not significantly changed if we had used 0.5 arcsec instead.

3.8. Completeness

Using the Monte Carlo source injections described in Sect. 3.7, we can easily derive the completeness for a given injected flux and size by computing the fraction of recovered sources with this property. In practice, the primary beam gain (G_pb) decreases quickly with the distance from the center and the noise is much larger on the edges of the maps. Consequently, the completeness depends strongly on the distance between the source and the center. However, the local noise can be easily computed by dividing the noise in the center σ_center (estimated in the nonprimary-beam corrected map) by the local primary-beam gain (G_pb). If we inject sources with similar nonprimary-beam-corrected flux (G_pb S_inj), they will have similar S/N whatever their distance to the phase center. The actual flux density, which is corrected by the primary beam gain, of these injected sources will be larger on the edge than in the center. In Fig. 6 (upper panel), we present the completeness as a function of G_pb S_inj. For clarity, we only show the results for point sources. While the completeness tends to zero at low flux and unity at high flux, the flux at which the transition appears varies significantly from pointing to pointing.

Fig. 6. Upper panel: completeness as a function of the continuum nonprimary-beam-corrected flux density (Gpb Sinj) achieved for point sources in various pointings. Middle panel: similar figure after having divided the nonprimary-beam-corrected flux by the noise at the center of each pointing (σcenter). Various colors (blue, green, yellow, and red) corresponds to various injected source sizes (FWHM = 0, 0.333, 0.666, and 1 arcsec, respectively). The solid lines indicate the mean trend of the various pointings, while the dashed lines indicate the 1-σ envelop. Lower panel: same plot after normalizing the injected flux by 1/σcenter and by the source area (Ωsource/Ωbeam). These results are discussed in Sect. 3.8.

When we normalized the injected nonprimary-beam corrected flux by 1/σ_center (middle panel), all the pointings have a very similar completeness curve for point sources (in blue). However, the completeness is not the same for all source sizes. At fixed normalized flux G_pb S_inj/σ_center, the completeness is lower for larger sources. A similar trend was found in the ALMA-GOODS deep field (Franco et al. 2018). We can also remark that a larger scatter from field to field is obtained for larger source size.

We used both the normal and the tapered maps to detect our sources. However, the S/N is usually higher in the normal map. The peak flux density in the normal map thus is a better proxy than the integrated flux to guess if a source will be detected or not by our algorithm searching for S/N peaks. We thus divided our previously-normalized flux densities by Ω_source/Ω_beam (see Sect. 3.5) to obtain a good proxy for the effect of the source size on the detectability. With this last correction, the completeness does not depend significantly on the source size and the scatter between pointings is highly reduced for the extended sources (see Fig. 6 lower panel). We derived the average curve for all sizes and pointings. The median distance to this average relation is only 1.2% with a maximum of 4.7%. We can thus reliably estimate the completeness based on this average relation from the source size, the primary-beam gain at its position, and its flux density.

3.9. Photometric accuracy and flux boosting

We also used our Monte Carlo simulations to test the accuracy of our photometry. In Fig. 7, we show the mean ratio between the recovered and injected flux density for our various photometric methods. For the 2D-fit photometry, the aperture photometry, and the peak flux in the case of point sources only, we observe the classical flux boosting effect at low S/N. Indeed, the sources with an injected flux density corresponding to an intrinsic S/N slightly lower than the detection threshold will be detected only if they are on a peak of noise. Their flux densities will thus be overestimated on average. In contrast, at high S/N, we expect that the output-versus-input flux density ratio will tend to unity, since sources located on both positive and negative fluctuations of the noise are detected. The 2σ-clipped method and the peak photometry of extended sources is more problematic and the results vary significantly with the size. In particular, even close to the S/N threshold, the flux densities are underestimated on average for a source size of 1 arcsec. At high S/N, the 2σ-clipped method converges slowly to unity. As expected, there is no convergence for the peak photometry, since the flux of all extended sources is systematically underestimated even in absence of noise and thus at high S/N.

Fig. 7. Ratio between the injected and recovered flux density as a function of measured S/N (see Sect. 3.9). Upper left, upper right, lower left, and lower right panels: results obtained for the 2D-fit, aperture, 2 σ-clipped, and peak photometry, respectively. The solid lines indicate the median and the shaded areas are the 1-σ contours. Various colors (blue, green, yellow, and red) are used to indicate the various sizes used (FWHM = 0, 0.333, 0.666, and 1 arcsec, respectively). The dashed horizontal line indicate the one-to-one relation.

We used these results to compute the flux boosting correction to apply. We computed the flux boosting correction at the S/N of the source for the immediately lower and higher sizes and used a linear interpolation to derive the correction to adopt for our source size.

To summarize, the peak flux density systematically underestimates the actual flux density of extended sources. Concerning the 2σ-clipped method, the flux boosting converges very slowly at high S/N and the flux boosting is highly size-dependent. Both aperture and 2D-fit photometry provide good results. We decided to use the 2D-fit photometry, because of the very small impact of the size on the deboosting correction to apply. In the following sections of this paper, we use the 2D-fit measurements. The raw and deboosted flux densities obtained using the 2D-fit method are listed in Table B.3.

3.10. Effective survey area associated with nontarget sources

To derive surface density of sources (also called number counts, see Sect. 5.3) or luminosity functions of nontarget sources (Gruppioni et al. 2020), we need to know the effective surface area of our survey as a function of the source properties. Of course, it varies with the flux density, since only the brightest sources can be detected on the edges of the pointing. It also depends on source size, since compact sources have usually a better completeness at fixed flux density (Sect. 3.8 and Fig. 6).

In addition, each pointing is observed at a slightly different frequency. We thus have to take into account that a source detected at a given flux density in a pointing will have a slightly different flux density in another pointing because of the different observed frequency, and consequently a slightly different completeness. For this reason, we apply a frequency-dependent correction factor to convert all the flux densities to 850 μm (353 GHz) assuming the z = 2.5 main-sequence spectral energy distribution (SED) template of the Bethermin et al. (2017) model (see the redshift distribution of nontarget sources in Sect. 5.2 and Fig. 12). Since most of the nontarget sources are at z <  4 and thus observed in the Rayleigh-Jeans part of their spectrum, the continuum slope around 850 μm does not vary significantly with the redshift and it is thus a fair assumption to assume a single template.

The effective surface area Ω_eff as a function of the source flux density S₈₅₀ and the source size θ_source is derived from the completeness C(S₈₅₀, θ_source, x, y) at a position (x, y) (see Sect. 3.8) using:

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\mathrm{eff}}(S_{850},{\theta }_{\mathrm{source}})={\displaystyle \sum _{\mathrm{pointings}}\int \int C(S_{850},{\theta }_{\mathrm{source}},x,y)\mathrm{d}\mathrm{\Omega }.}\end{array}$$(5)

Since the nontarget sources are extracted outside the central 1 arcsec-radius region, we exclude this area from the computation of the integral.

The result is presented in Fig. 8. As expected, the surface area at intermediate flux densities varies significantly with the source size. At bright flux densities (> 10 mJy), the completeness tends to unity and the effective surface area is the total area of all our pointings⁵ (24.92 arcmin²). Our survey is ∼3 times smaller than ALMA-GOODS (Franco et al. 2018) for a similar sensitivity in mJy. However, typical galaxies are fainter by a factor of ∼2 at 1.1 mm. Our band-7 serendipitous survey thus is a valuable complement to the band-6 deep fields (Dunlop et al. 2017; Aravena et al. 2016; González-López et al. 2017; Franco et al. 2018).

Fig. 8. Effective surface area of ALPINE as a function of the 850 μm flux after excluding the central 1 arcsec-radius area where target sources are extracted. The blue, green, gold, and red lines are the results obtained for a source size of 0, 0.3, 0.6, and 1 arcsec, respectively. The method used to compute the surface area is described in Sect. 3.10.

4. From rest-frame 158 μm continuum fluxes to SFR

4.1. Dust spectral energy distribution variation from low-redshift to high-redshift Universe

The obscured star formation is directly related to the bolometric luminosity of the dust (SFR_IR = 1 × 10⁻¹⁰ M_⊙ yr⁻¹/L_⊙ × L_IR, Kennicutt 1998 after converting to Chabrier 2003 IMF). L_IR is usually defined as the total luminosity of a galaxy between 8 and 1000 μm. ALPINE continuum photometry is only probing a narrow range of wavelength around 158 μm rest-frame. Since we have only one photometric point available, we thus have to assume a spectral energy distribution (SED) to derive L_IR. As discussed in Bouwens et al. (2016), Fudamoto et al. (2017), and Faisst et al. (2017), for example, the assumption on the dust temperature of z >  4 galaxies has a significant impact on the relation connecting the dust attenuation to the UV continuum slope β or the stellar mass, which will be discussed in Fudamoto et al. (2020).

While the SEDs of z <  2 galaxies have been well studied thanks to Herschel (e.g., Elbaz et al. 2011; Dunne et al. 2011; Magdis et al. 2012; Berta et al. 2013; Symeonidis et al. 2013; Magnelli et al. 2014), we have fewer constraints on the SEDs at higher redshifts. These z <  2 studies revealed that the temperature of normal, star-forming galaxies tends to increase with redshift, which agrees with the theoretical model predictions (e.g., Cowley et al. 2017a; Imara et al. 2018; Behrens et al. 2018). Because of the confusion noise, Herschel can detect only the brightest galaxies (e.g., Nguyen et al. 2010). However, some interesting constraints up to z ∼ 4 were obtained using stacking analysis of galaxies selected using photometric redshifts (e.g., Béthermin et al. 2015a; Schreiber et al. 2015), Lyman-break selections (e.g., Álvarez-Márquez et al. 2016), and low-redshift analogs of z >  5 galaxies (Faisst et al. 2017). According to these studies, temperature seems to continue to increase up to z ∼ 4. So far, we have very few constraints about what happens at z >  4, which is critical to interpret the ALPINE survey.

In this section, we present a stacking analysis adapted from Béthermin et al. (2015b) to derive an average empirically-based conversion from the 158 μm monochromatic continuum flux density to L_IR and SFR.

4.2. Mean stacked SEDs of ALPINE analogs in the COSMOS field

Béthermin et al. (2015b) used a mean stacking analysis (without source weighting) of Herschel and complementary ground-based measurements in the COSMOS field to derive the mean SEDs of z <  4 galaxies. We used the same Herschel⁶ (Pilbratt et al. 2010) data from the PEP (Lutz et al. 2011) and HerMES (Oliver et al. 2012) surveys and AzTEC/ASTE data of Aretxaga et al. (2011) at 1.1 mm. At 850 μm, we used the SCUBA2 data from Casey et al. (2013) instead of the shallower LABOCA ones used in the 2015 analysis.

The 2015 selection of the stacked targets was performed using a stellar mass cut of > 3 × 10¹⁰ M_⊙ in the photometric Laigle et al. (2016) catalog. There are too few ALPINE sources to obtain a sufficiently high S/N in the stacked Herschel data. We thus used a larger photometric sample with properties similar to ALPINE objects. We chose to select sources with an estimated SFR from an optical and near-infrared SED fitting higher than 10 M_⊙ yr⁻¹, which is approximately equivalent to the ALPINE SFR limit (Le Fèvre et al. 2020; Faisst et al. 2020).

We also use higher redshift bins (4 <  z <  5 and 5 <  z <  6) to match the redshift range probed by ALPINE. Finally, we use the more recent COSMOS catalog of Davidzon et al. (2017) as input sample, since it has been optimized to provide more reliable photometric redshifts and physical parameters at z >  4. Our stacked samples contain respectively 5749 and 1883 sources in the 4 <  z <  5 and 5 <  z <  6 ranges.

Our new stacking analysis was performed using the exact same procedure as in Béthermin et al. (2015b). The uncertainties were derived using a bootstrap technique that takes into account both the photometric noise (instrumental and confusion) and the population variance. The contamination of the stacked flux by clustered neighbors is corrected using the method described in Appendix A of Béthermin et al. (2015b). At z >  4, these corrections are relatively small (< 30%) because of the lower global star formation rate density compared to z = 2. Our results are presented in Fig. 9 and Table 2. The 5 <  z <  6 SED is ∼2 times noisier mainly because of the smaller number of stacked objects.

Fig. 9. Comparison between the Álvarez-Márquez et al. (2016, blue dashed line), Schreiber et al. (2018, orange dot-dashed line), and Bethermin et al. (2017, red solid line) IR SED templates and the observed mean SEDs of SFR > 10 M⊙ yr−1 galaxies measured by stacking (black dots, see Sect. 4.2). The black dotted line is the best fit of the λrest − frame >  40 μm data points by a modified blackbody with β fixed to 1.8 (the temperature in the legend is provided in the rest frame). Upper and lower panels: 4 <  z <  5 and 5 <  z <  6, respectively.

4.3. SED template and conversion factors

The final step to compute the conversion factor from monochromatic luminosity to L_IR is to find an SED model or a parametric description fitting the data. Using an agnostic model as a spline is difficult, since we have few constraints on the mid-infrared (λ_rest <  30 μm). In Fig. 9, the SEDs are represented in flux density units (S_ν = dS/dν), which can give the wrong impression that the contribution at short wavelength is negligible, while it contains ∼15% of the energy⁷. For this reason, although fitting well the Herschel data points, a modified blackbody (ν^β B_ν(ν, T), where B_ν is a blackbody law) tends to underestimate the L_IR because of the very low emission in the mid-infrared. For information, we show in Fig. 9 the best fit of our SEDs by a modified blackbody with a fixed β of 1.8, but a free amplitude and temperature. We excluded rest-frame wavelengths below 40 μm from the fit, since the greybody model does not take into account the warm dust and the polycyclic aromatic hydrocarbon (PAH) features dominating in this wavelength range. The fit is excellent at 4 <  z <  5 (χ² = 0.44 for 3 degrees of freedoms) and acceptable at 5 <  z <  6 (χ² = 3.55 for 2 degrees of freedoms).

We thus chose to use empirical template libraries. We compare our observed SEDs with three different templates. Álvarez-Márquez et al. (2016) template⁸ is based on the stacking of 2.5 <  z <  3.5 Lyman-break galaxies. The SED templates for main sequence galaxies of the Bethermin et al. (2017) model evolves with redshift up to z ∼ 4. Above this redshift, no evolution is assumed (⟨U⟩ = 50). These templates are an update of the Magdis et al. (2012) templates calibrated using the Herschel stacking up to z ∼ 4 (Béthermin et al. 2015b). Finally, Schreiber et al. (2018) also built a template evolving with redshift and calibrated it using another independent Herschel stacking analysis. Contrary to the previous templates, they assume an evolution of the rest-frame dust temperature above z = 4 (4.6 K per unit of redshift). Because of the nature of the ALPINE sample (Faisst et al. 2020), we only consider the templates corresponding to galaxies on the main sequence.

In Fig. 9, we show the comparison between our measured SED and the templates described above. We renormalized the templates to fit the data. This is the only free parameter in our analysis. While the Álvarez-Márquez et al. (2016) template is too cold for both redshift bins, both Schreiber et al. (2018) and Bethermin et al. (2017) templates well fit the data (χ² <  4 with 4 degrees of freedom for both templates in both redshift bins). Since the χ² of Bethermin et al. (2017) is marginally better, we decided to use this template. In Table 3, we provide the ratio between the monochromatic luminosity (νL_ν units) and L_IR computed using this template at wavelengths associated with bright fine-structure lines, which can be targeted by ALMA. In practice, for the ALPINE catalog (Appendix B), we use the exact effective wavelength of the ALMA continuum.

4.4. Caveats

The conversion factors derived previously are based on the best effort, but they are clearly not the final answer about this complex topic. First of all, the selection of the stacked sample is not perfect and based on photometric redshifts and SFRs derived from rest-frame UV to near-IR SED fitting. It is also difficult to estimate how similar this SFR selection is compared to the actual ALPINE sample. Even if it is not likely, we could imagine that a population with very peculiar dust SEDs is missing in one of the two samples.

The stacked SED was obtained by averaging all the galaxies from our stacked sample. The derived conversion factors could thus be largely inaccurate for outliers with extreme dusty SEDs. Finally, even if the same weight is attributed to each source, stacking provides luminosity-weighted mean SEDs, since brighter sources will have a larger relative contribution to the final signal. We could imagine that a population, which represents a significant fraction of the sample in number but contributes little to the luminosity, has an extreme SED. The stacking analysis would miss such objects and their individual L_IR estimates could be incorrect.

In Fig. 10, we present the stacking for a larger SFR cut of 100 M_⊙ yr⁻¹. According to optical and near-infrared SED fitting (Faisst et al. 2020), only 11 out of our 118 sources are following this criterion. This analysis is only possible in the 4 <  z <  5 bin, since there is no detection at higher redshift. For these objects, the dust temperature is warmer (47 K versus 41 K) and the Schreiber et al. (2018) template fits better the data. The consequences of a slightly warmer dust at higher SFR will be discussed in Sect. 7.5.

Fig. 10. Same figure as Fig. 9 (upper panel), but using a SFR > 100 M⊙ yr−1 cut.

5. Continuum source properties

In this section, we discuss the properties of our continuum detections. In Sect. 5.1, we discuss briefly the basic properties of the detected target sources. In the following sections, we focus on the properties of the nontarget detections: redshift distribution (Sect. 5.2), number counts (Sect. 5.3), and contribution to the cosmic infrared background (CIB, Sect. 5.4).

5.1. Properties of the target sources

The redshift distribution of the detected target sources is presented in Fig. 11 (upper left panel). While the detections are distributed across most of the redshift range of the total sample, the detection rate is slightly better in the lower redshift window (26 ± 6%) than in the high redshift window (15 ± 5%). Since the sensitivity at fixed luminosity is better in the z >  5 redshift window (Sect. 2.6), we could have expected the opposite trend. However, this is only a 1.4σ difference and the dust content could be lower at higher redshift. The dust attenuation of our detections will be discussed in Fudamoto et al. (2020).

Fig. 11. Upper left panel: redshift distribution of the ALPINE target sources. The blue and red histograms are the distribution of detected sources only and the full sample, respectively. Upper right panel: distribution of the continuum flux densities. The red histogram indicates the distribution expected from the optical and near-IR SED-derived SFR assuming the long-wavelength SED presented in Sect. 4. Lower left panel: same figure as previously but for LIR. Lower right panel: stellar mass distribution of detected (blue) and all (red) sources.

We can also compare the flux density distribution of our detections and the expected distribution from the ancillary data (Faisst et al. 2020, version including Spitzer photometry in the SED fitting). To produce the expected ALPINE flux densities from ancillary data, we estimated the expected L_IR from the SFR based on optical and near-infrared SED fitting assuming a 1 × 10⁻¹⁰ L_⊙/(M_⊙ yr⁻¹) conversion factor (see Sect. 4.1). By doing so, we assume implicitly that the infrared traces the entire star formation. Finally, we use the long-wavelength SED template presented in Sect. 4.2 to predict the flux density. The results are shown in Fig. 11 (upper right panel). The most extreme predicted flux densities (> 2 mJy) are not found in the real sample. These very high SFR are almost certainly due to overestimated dust-attenuation corrections. In contrast, all the detected objects are above the mode of the predicted distribution. This shows that we are sensitive only to the highest SFRs. However, this is not a sharp cutoff. This demonstrates that the measured distribution could not have been predicted from the ancillary data and that submillimeter data are important to derive reliable SFRs. A similar trend is found for the infrared luminosity L_IR (see lower left panel).

Finally, we compared the stellar mass distribution of the full sample and of the detections only (lower right panel). The mass distributions of the full sample and of the detections are significantly different according to the Kolmogorov-Smirnov test (p-value = 3.8 × 10⁻⁵) and only two detections are below the median stellar mass of our full sample. This is an expected consequence of the correlation between the stellar mass and the star formation rate often called main sequence (e.g. Schreiber et al. 2015; Tasca et al. 2015; Khusanova et al. 2020).

5.2. Redshift distribution of the nontarget continuum detections

Contrary to the target sample, determining the redshift of nontarget sources is not trivial. We have to identify the optical/near-infrared counterparts and use photometric redshifts when spectroscopic redshift are not available. Fortunately, this sample lies in survey areas with rich ancillary data (fully described in Faisst et al. 2020) drawn primarily from COSMOS (Scoville et al. 2007), GOODS (Giavalisco et al. 2004) and CANDELS (Grogin et al. 2011; Koekemoer et al. 2011). The counterparts of 42 of our 57 nontarget continuum detections were identified in the Laigle et al. (2016, COSMOS) or the Momcheva et al. (2016, 3DHST) catalogs. The detailed identification of each source and the sources without counterpart in the previously cited catalogs will be discussed in details in Gruppioni et al. (2020).

The redshift distribution of our nontarget sources is presented in Fig. 12. The mean redshift of our sample is z = 2.5 ± 0.2 (median = 2.3 ± 0.3) with a tail up to z = 6. Our uncertainties are computed using a bootstrap technique and thus include sample variance. This is 1-σ lower than the median redshift (z = 2.65 ± 0.13) found by Simpson et al. (2017) following up > 1 mJy sources selected in a single-dish survey using ALMA at the same wavelength. As shown in Béthermin et al. (2015a, see also Hodge & da Cunha 2020), fainter submillimeter sources are paradoxically expected to have a lower mean redshift. This small difference is thus not surprising. To test if this trend is also found inside our own sample, we split it into two equally populated subsamples containing the faint (< 1.47 mJy) and the bright sources (> 1.47 mJy). The faint sources have a mean redshift of z = 2.6 ± 0.3, while we found z = 2.3 ± 0.2 for the bright ones. It agrees with the trend predicted by Béthermin et al. (2015a), but our sample is too small to provide a statistically significant result. We also expect that longer wavelengths probe higher redshifts. As expected, our median redshift is smaller than what is found at 1.1 mm by Franco et al. (2018, z_med = 2.9) and Brisbin et al. (2017, z_med = 2.48 ± 0.05) or 1.4 mm by Strandet et al. (2016, z_med = 3.9). This last sample is lensed and it might push the median redshift to higher values. In contrast, our median redshift is higher than the very faint 1.1 mm sample of Aravena et al. (2016, down to 0.05 mJy, z_med = 1.9 ± 0.4). As shown in Fig. 3 of Béthermin et al. (2015a), it is expected that < 0.1 mJy 1.1 mm sources are at lower redshift than ∼1 mJy 850 μm sources.

Fig. 12. Photometric redshift distribution of the nontarget ALPINE sources (blue filled histogram). The dashed, solid, and dotted histograms are the predictions from the Bethermin et al. (2017) SIDES simulation for a flux cut corresponding to the first-quartile, the median, and the third-quartile of the observed sample, respectively.

Finally, we compare our measured distribution with the predictions of the simulated infrared dusty infrared sky (SIDES) simulation⁹ (Bethermin et al. 2017). Since the depth of our various pointings is not homogenous, our sample is not flux limited. We thus computed the redshift distribution for different flux cuts corresponding to the first quartile, median, and third quartile of the observed sample. The three predicted distributions are compatible with our measurements at 1-σ. Because of galaxy clustering, we could have expected an excess of sources at the same redshift as the ALPINE targets, but we observe only a 1-σ excess between z = 5 and z = 6. We can thus assume that the sample of nontarget sources with optical counterparts is statistically similar to a sample, which would have been obtained using random pointings.

5.3. Number counts

The area probed by the ALPINE survey is sufficiently large to produce new meaningful constraints on the faint galaxy number counts at 850 μm. While the ALMA band 6 (1.1–1.4 mm) was extensively used for deep surveys, the band 7 (∼850 μm) has been much less explored, especially below 3 mJy. Oteo et al. (2016) provided constraints based on the ALMA calibrator survey, but with large uncertainties. In contrast, this wavelength has been widely explored with single-dish instruments (e.g., Coppin et al. 2006; Casey et al. 2013; Chen et al. 2013; Hsu et al. 2016; Geach et al. 2017). However, because of their limited spatial resolution, several galaxies can be blended in the same beam, biasing the bright number counts toward higher values (Hayward et al. 2013; Karim et al. 2013; Bussmann et al. 2015; Cowley et al. 2017b; Scudder et al. 2016; Bethermin et al. 2017). Above 3 mJy, interferometric observing campaigns had followed up single-dish sources to correct for this effect (Karim et al. 2013; Simpson et al. 2015; Stach et al. 2018).

We derived the integral number counts dN/dΩ, which is the surface density of sources above a certain flux cut, by summing the inverse of the effective area Ω_eff for each nontarget source¹⁰:

$$\begin{array}{c}\hfill \frac{\mathrm{d}N}{\mathrm{d}\mathrm{\Omega }}(>S)={\displaystyle \sum _{i=1}^{N_{\mathrm{source}}}\frac{1}{{\mathrm{\Omega }}_{\mathrm{eff}}(S_i,{\theta }_i)},}\end{array}$$(6)

where S_i and θ_i are the deboosted flux density (see Sect. 3.9) and size of the ith source. All the flux densities have been converted to 850 μm at which the effective area (see Sect. 3.10) was computed. As shown in Sect. 3.9 and Fig. 7, the deboosting factor, which is necessary to apply here, can have a 30% uncertainty. We thus computed the difference between the number counts derived using the 1-σ lower and upper envelopes of the flux boosting curve to estimate the associated uncertainties. These uncertainties are combined with the Poissonian error bars. Finally, SC_1_DEIMOS_COSMOS_787780, SC_1 _vuds_cosmos_5101210235, and SC_2_ DEIMOS_ COSMOS_ 773957 are detected by our algorithm only because their flux density is boosted by a line. Else, their S/N without line contimination falls below our threshold of 5. These sources are thus excluded from our continuum number count computation.

A similar method was used to derive the differential number counts. We summed the inverse of the effective area of all the sources in a given bin and divided by the bin size. To reduce the dynamical range on the figures, we normalized the differential counts by S^2.5. With this normalization, the number counts in an Euclidian nonevolving Universe are flat. This is usually the case for very bright fluxes (> 100 mJy) in the submillimeter domain (Planck Collaboration VII 2013), where the detected sources are mainly local. The deviations from this trend at fainter flux densities provide important constraints for galaxy evolution models. It has also the convenient property to reduce the dynamical scale of the plot and help the visual comparison between the models and the data.

We estimated the integral number counts for various thresholds spaced by 0.2 dex. We chose to use 0.35 mJy for the lowest threshold, which corresponds to the deboosted 850 μm-converted flux density of the second faintest object. Concerning the differential number counts, we used the intervals delimited by this list of thresholds. We do not use fainter bins, since the faintest object (0.30 mJy after conversion to 850 μm, SC_2_DEIMOS_COSMOS_773957) is associated to a completeness of 13% and the correction to apply is thus very large. The mean completeness in the faintest bin (0.35–0.56 mJy) is 50%. In all the other bins, the mean completeness is above 80%.

Our measurements are summarized in Tables 4 and 5 and shown in Fig. 13. Even if the redshift distribution of the nontarget sources provides no firm evidence for it (see Sect. 5.2), we cannot formally exclude a small overdensity of sources at the same redshift as ALPINE targets compared to a random position in the sky. We thus estimated the number counts using both the full sample and a secure z <  4 sample, where only the sources with identified optical or near-IR counterparts below the ALPINE redshift range are kept. These two samples provide respectively an upper and a lower limit on the number counts, which would be derived at a random position in the sky. The values derived using these two samples agree at a 1σ level.

Fig. 13. Left panel: differential number counts at 850 μm. The blue circles are the counts derived from the full ALPINE nontarget sample. The red squares are computed from the ALPINE nontarget sources with a confirmed optical or near-IR counterpart below z = 4 (secure z <  4 sample). It is thus a secure lower limit, since sources clustered with ALPINE sources are excluded. The green crosses and pentagons are the interferometric measurements of Karim et al. (2013) and Stach et al. (2018), respectively. The gray left-facing triangles, right-facing triangles, diamonds, and down-facing triangles are the measurements performed using single-dish data with a lower angular resolution by Geach et al. (2017), Hsu et al. (2016), Casey et al. (2013), and Chen et al. (2013), respectively. The solid, dashed, dotted, and dash-dot lines are the models of Bethermin et al. (2017, SIDES simulation), Gruppioni et al. (2011), Casey et al. (2018, low dust), and Casey et al. (2018, high dust). Right panel: integral number counts at 850 μm. The symbols are the same as in the left panel. The interferometric measurements of Oteo et al. (2016) and Simpson et al. (2015), respectively, are represented by green stars and triangles.

Our new measurements exhibit a shallower slope below ∼3 mJy than what was measured by previous surveys at higher flux density. This is the first time that we probe so well the regime below this slope break, since Oteo et al. (2016) data points suffer from an order of magnitude of uncertainties and the Hsu et al. (2016) analysis is affected by their low angular resolution and rely on complex methods to invert the lensing to recover the intrinsic counts. Our results are compatible at 1σ with Oteo et al. (2016, only shown for integral number counts, i.e. the right panel of Fig. 13) and Hsu et al. (2016).

We can also compare our new measurements with various models of galaxy evolution. The Bethermin et al. (2017) model is an update of the Béthermin et al. (2012a) models, which decomposes star-forming galaxies into main-sequence and starbursts galaxies. Each type has a different SED evolving with redshift. It starts from the observed stellar mass function and the measured evolution of the main sequence of star-forming galaxies to predict infrared and (sub-)millimeter observables. This model includes clustering and can reproduce the single-dish low-resolution data and interferometric data simultaneously. However, there were few constraints below 3 mJy at 850 μm, when it was published. The Gruppioni et al. (2011) model updated in Gruppioni & Pozzi (2019) is based on the observed luminosity functions (Gruppioni et al. 2013). This model contains five different populations, whose luminosity functions evolve with redshift (spirals, starbursts, low-luminosity AGN, type-1 AGN, and type-2 AGN). The Casey et al. (2018) model assumes an evolution of both the infrared luminosity function and the SEDs. Two versions were proposed. The first one has a low number of dusty objects at high redshift and the other assumes a large volume density of bright obscured galaxies.

At the spatial resolution of the ALPINE data (∼1 arcsec), we can directly compare the source counts with the galaxy counts in the model, since it is unlikely to have two physically-separated galaxies in the same beam. Overall, all models agree with our data, since they are between the 1σ lower limit of the secure z <  4 and the 1σ upper limit of the full sample. However, the high-dust model of Casey et al. (2018) is higher than the Stach et al. (2018) measurements at high flux density. Below the ∼3 mJy slope break, the predictions of the Gruppioni et al. and the low-dust Casey et al. models are a factor of 2 lower than the high-dust Casey et al. and the Béthermin et al. models. Unfortunately, our uncertainties are still too large to identify the correct scenario.

5.4. Cosmic infrared background

The cosmic infrared background (CIB) is the relic of all dust emission by galaxies across cosmic times (e.g., Dole et al. 2006). Its absolute brightness was measured in the nineties by the COBE/FIRAS instrument (Puget et al. 1996; Hauser et al. 1998; Fixsen et al. 1998; Gispert et al. 2000; Lagache et al. 2000). By combining FIRAS and Planck, Odegard et al. (2019) estimated the absolute CIB level with a better precision than historical analyses. Their values are used as the reference CIB level in this paper. The CIB SED provides a budget of far-infrared photons that galaxies emit during their evolution. While Herschel identified the galaxy populations (luminosity, redshift) emitting the CIB below 500 μm, there were fewer constraints at longer wavelength before ALMA. However, new number counts extracted from band-6 surveys (∼1.3 mm) can explain 50−100% of the CIB absolute level (Carniani et al. 2015; Aravena et al. 2016; Fujimoto et al. 2016). At 850 μm, using single-dish SCUBA2 data of cluster fields, Hsu et al. (2016) found that the full CIB can be explained by galaxies brighter than 0.1 mJy.

The CIB brightness produced by all the sources above a certain flux density threshold is:

$$\begin{array}{c}\hfill B_{\nu }={\displaystyle {\int }_{S_{\lim }}^{\infty }{S}_{\nu }\frac{\mathrm{d}N}{\mathrm{d}{S}_{\nu }}\mathrm{d}{S}_{\nu },}\end{array}$$(7)

where B_ν is the surface brightness density of the CIB produced by sources above the flux density limit S_lim. To compute this integral, we assume a power-law to connect our data points. We extrapolated the contribution of sources fainter than 0.35 mJy by fitting a power-law to the five faintest data points. The slope is poorly constrained and is responsible for large uncertainties on this extrapolation: $\mathrm{d}N/\mathrm{d}{S}_{\nu }\propto S_{\nu }^{-2.2\pm 0.3}$ for the full sample and $\mathrm{d}N/\mathrm{d}{S}_{\nu }\propto S_{\nu }^{-1.8\pm 0.4}$ for the secure sample. Above our brightest data point, we use the Euclidian plateau level measured by Planck Collaboration VII (2013, dN/dSS^2.5 = 15 Jy^1.5 sr⁻¹). The uncertainties are determined by recomputing 100 000 times the integral of the number counts after randomly offsetting the data points according to their error bars.

In Fig. 14, we present the contribution to the CIB as a function of the flux density limit. The contribution of sources brighter than the limit of our number counts (0.35 mJy) is 0.093 ± 0.013 MJy sr⁻¹ for the full sample and 0.054 ± 0.009 MJy sr⁻¹ for the secure sample. This is 69% and 40% of the full CIB measured by Odegard et al. (2019, 0.135 MJy sr⁻¹ assuming a CIB spectrum for the color correction). Below 2 mJy, our measurements agree with Hsu et al. (2016), since their data are between the values determined from our secure z <  4 measurements (lower limit) and from our full sample (upper limit).

Fig. 14. Contribution of galaxies to the cosmic infrared background at 850 μm as a function of the flux density limit. The blue and red solid lines are our results based on the full and the secure z <  4 samples, respectively. The shaded area is the 1-σ confidence region. The dashed green line is the measurement of Hsu et al. (2016) using SCUBA2 data of cluster fields. The black line is the absolute measurement of the total CIB combining COBE/FIRAS and Planck data (Odegard et al. 2019).

6. [CII] catalog

In this paper, we focus only on the [CII] line detections of the ALPINE targets. The serendipitous detections and [CII]-emitting close companions of the targets will be discussed in Loiacono et al. (2020).

6.1. Extraction of the candidates

Extracting lines from a cube can be significantly more difficult than the continuum when the redshift of the source is not well known, since we also have to explore the spectral dimension. This is the case for ALPINE, since we found significant offsets between the [CII] lines and our reference redshifts derived from optical spectroscopy (see Sect. 7.2). To perform this task, we used the semi-automatic procedure described below.

We first ran a customized line finder algorithm (see comparison with the findclumps algorithm in Sect. 6.2) to search for [CII] emission in a 3D area defined as the cylinder with a 1″ aperture around the phase center extending over the full bandwidth, to allow for spatial and spectral offsets. In practice, we produced running averaged cubes using various numbers of channels. The minimum allowed number of collapsed channels is two (50 km s⁻¹) to avoid single-channel spikes of noise, and the maximum is 40 channels (1000 km s⁻¹). In other words, we produce running moment-0 cubes across the full bandwidth. We then computed rms of these nonprimary-beam-corrected maps using all pixels at > 1″ from the phase-center and searched for S/N >  3 peaks in the maps. Finally, we saved the positions, frequencies and S/Ns of each of these [CII]-emitter candidates.

For each [CII]-emitter candidate, we extracted a spectrum from the continuum subtracted cubes at the position of the brightest pixel of the moment-0 map and had a further quality assessment based on visual inspection. For the good candidates, we fit a Gaussian profile of the spectrum and computed central frequency (f_cen) and FWHM. We then produced a new moment-0 map collapsing the channels included in the [f_cen − FWHM, f_cen + FWHM] spectral range. The choice of such frequency interval is an optimal compromise between including most of the line profile and excluding any noise at the edge of the line. This step is needed because the moment-0 map computed during the first step is the one that maximizes the S/N, and therefore, does not necessarily cover the full line profile. We then re-extracted a spectrum inside the 2-σ region measured in this moment-0 map. This new spectrum is more informative than the first one in the case of spatially-resolved sources, since it is extracted from the entire region of emission and not just at the peak. These steps were repeated several times until the shape of the 2-σ contour, in which the spectrum is extracted, is stable. Few iterations were necessary, since the flux measured in both the spectra and in the moment-0 maps converged within the uncertainties after only one iteration. The final spectra are presented in Fig. C.1.

6.2. Detection threshold and reliability

To assess the reliability of the line emitter candidates, we use the same method based on the negative maps as for the continuum (Sect. 3.1), but using the final moment-0 maps produced after optimizing the integration frequency window (see Sect. 6.1). Contrary to continuum maps, we have to take into account that we searched the line over a large velocity range, which could contain several times the line width. We thus normalized the number of detections in the negative map by the ratio between the velocity range in which we expect a line (∼1500 km s⁻¹, see Sect. 7.2) and the actual width of the velocity window used to produce the moment-0 map. We estimated that the 95% purity in the central 1″ region is reached for a S/N >  3.56 threshold. The purity as a function of the S/N is shown in Fig. 15 (blue solid line).

Fig. 15. Purity of [CII] detections as a function of the S/N threshold in the moment-0 map. The purity is provided only for the 1 arcsec central region and was estimated using the methods described in Sect. 6.2. The blue solid line is derived using the method based on the moment-0 maps and the red dashed line using an extraction in the negative cubes (see Sect. 6.2). The dotted vertical lines are the S/N at which the 95% purity is reached.

Since this normalization is an approximation, we cross-checked our number with an independent approach. We extracted the sources in the negative and the positive cubes using the findclumps (Decarli et al. 2016; Walter et al. 2016) routine. The S/N determined by findclump is on average 15% higher than in our moment-0 map. This is expected since findclump chooses an integration window to maximize the S/N. It is usually narrower and only includes the high-S/N central channels of the line, while our moment-0 maps include also the noisier tails. A similar difference of S/N estimates was also found by the ASPECS survey team (González-López et al. 2019). After correcting for this 15% systematic difference to agree with our moment-0 map S/Ns, we find a 95% purity limit for S/N = 3.45 in the 1 arcsec central region. The two methods are thus in excellent agreement. We also compared the FWHM of the detections in the positive cube and the spurious sources found in the negative cubes close from the detection threshold (3.5 <  S/N <  4.5). We found a mean value of 290 ± 40 km s⁻¹ and 216 ± 1 km s⁻¹, respectively¹¹. This 2-σ difference is a reassuring hint that the low S/N sources are not spurious.

We chose to cut our catalog at S/N >  3.56 to be conservative. The purity estimated using this method is presented in Fig. 15. Since the number density of line emitters in the rest of the field is much lower (Sect. 3.1), a much higher S/N threshold is necessary in the rest of the field. This will be discussed in Loiacono et al. (2020). Finally, we checked that the two extraction methods provide compatible source lists and found that this is true except for a couple of objects close to the S/N detection threshold.

It could seem counterintuitive that the S/N threshold to reach 95% purity for target sources is the same for the continuum and lines. Indeed, this is a coincidence. If the surface number density of continuum and line sources were the same at fixed S/N, the purity would be much lower for the lines, since line spurious sources can come from several velocity channels and are thus more numerous at fixed S/N. However, the number of S/N >  3.56 [CII] sources is three times higher than the number of continuum emitters. The higher number of real sources thus compensates the higher number of spurious sources in Eq. (1).

We obtained 75 [CII] detections out of our 118 targets. The moment-0 cutouts of our detected [CII] targets is presented in Fig. C.2. An extensive discussion of their morphology is presented in Le Fèvre et al. (2020).

6.3. [CII] flux measurements and consistency

The line flux of the detected [CII] sources was measured using the same methods as for the continuum (see Sect. 3.2), but using the moment-0 maps instead of continuum maps. This allows us to reliably measure the line flux of spatially-resolved sources. As discussed in Sect. 6.1, the conservative velocity windows used to build the moment-0 maps minimize the flux loss from the high-velocity tails.

In Fig. 16, we compare the results obtained by our various photometric methods. Except the peak flux method, which systematically underestimates the flux of extended sources mainly located at the bright end of the sample (see Sect. 3.2), we find a good overall agreement between the 2D-fit, aperture, and 2-σ-clipped fluxes. However, we identified a small systematic offset of 7% and 3% compared to the 2D-fit flux for the aperture and 2-σ-clipped flux, respectively. The uncertainty-normalized difference defined in Eq. (3) is 0.57 and 0.56, respectively. The systematic uncertainty between the methods is thus smaller than the typical 1-σ flux uncertainties. Finally, we identified a strong outlier for which our various methods disagree with each others. It is a source with a close bright neighbor, which requires manual deblending (see Appendix D.2 and Ginolfi et al. 2020b).

Fig. 16. Comparison between our various [CII] line flux measurement methods described in Sect. 6.3. The blue dots are our measurements and the red line is the one-to-one relation. Upper left, upper right, and lower left panels: comparison between the 2D-fit flux (x-axis) and the aperture, 2-σ-clipped, and peak flux densities, respectively (y-axis). The vuds_cosmos_5101209780 object (in green) has a complex geometry and is discussed in Appendix D.2. The lower right panel shows the comparison between uv-plane fluxes and 2D-fit fluxes for single-component systems.

We also compared our flux measurements performed in the map space with uv-derived values. These measurements were performed using the uvmodelfit CASA routine. The full measurement process and the comparison of the results obtained assuming various profile models will be presented in Fujimoto et al. (in prep.). In this paper, we compare our measurements with the results obtained in the uv plane using an elliptical Gaussian model. In Fig. 16 (lower right panel), we compare the fluxes from the map space and the uv plane for the sources properly fit in the uv space by a single component. The two methods agree within an offset of 7%.

The results are presented in Table C.1. The fluxes in the table are estimated using the 2D-fit method. For the non detections, we derived upper limits using the same methods as for the continuum (Table B.2) but derived from the moment-0 maps. Since the line is not detected, we have to decide which channels to use to produce a moment-0 map. We chose to use a 300 km s⁻¹ window centered on the optical redshift, which is slightly above the median FWHM of the detected sample (see Sect. 7.3). Of course, these upper limits do not apply if the line is particularly broad or if there is a catastrophic error on the reference optical redshift. The results are provided in Table C.2.

7. [CII] target properties

In this section, we describe the basic [CII] properties of the target sources. In Sect. 7.1, we present the redshift distribution of the detections and discuss the detection rate. The velocity offsets between [CII] and optical redshifts are presented in Sect. 7.2. We discuss the line width and the luminosities in Sects. 7.3 and 7.4, respectively. Finally, we use the ALPINE sample to constrain the relation between the SFR and the [CII] luminosity (Sect. 7.5).

7.1. Redshift distribution and detection rate

In Fig. 17 (upper left panel), we show the distribution of the reference optical redshifts of the full ALPINE sample and the redshift distribution of the subsample detected in [CII]. Overall, the distribution of the detections follows the full sample distribution. The detection rate below and above z = 5 is 66% and 58%, respectively. This corresponds to a 1-σ statistical fluctuation and cannot be considered as significant. This is not surprising, since our survey was built to be at constant [CII]-luminosity sensitivity and no variation with redshift of the average sensitivity in luminosity was found in the real data (Sect. 2.6).

Fig. 17. Upper left panel: redshift distribution of the full sample (red open histogram) and the [CII] detections (blue filled histogram). Upper right panel: distribution of the velocity offset between the [CII] and the reference optical redshift. Lower left panel: distribution of the [CII] line FWHM of the detected ALPINE sources (blue filled histogram) and comparison with the lensed galaxy sample of Gullberg et al. (2015, red) and the quasar sample of Decarli et al. (2018, green). Lower right panel: comparison between the [CII] luminosity distribution of the detections (blue filled histogram) and the expected values based on the SFRs determined by optical and near-infrared SED fitting and the De Looze et al. (2014, red open histogram) and Lagache et al. (2018, black open histogram) relations.

7.2. Offsets between optical and [CII] redshifts

For the detected [CII] sources, we also compared our optical reference redshift and our new [CII] redshift (see Fig. 17, upper right panel). We found non negligible offsets up to 1000 km s⁻¹. A posteriori, this justifies our choice to search for the line not only at the optical velocity, but in the entire side band (see Sect. 6.1). These offsets cannot be explained by the uncertainties on the [CII] redshift, since they are usually smaller than 100 km s⁻¹. The optical redshifts could suffer from several effects leading to a small inaccuracy. First, for sources without bright emission lines, redshifts are determined using the continuum and/or weak absorption features and could have a lack of precision (Le Fèvre et al. 2015; Hasinger et al. 2018). Other sources can have a very bright Lyman α line, which has a lot of weight in the determination of the redshift. While [CII] traces the gas in the galaxy and is very close to the systemic velocity, Lyman α radiative transfer is complex and can produce significant offsets (e.g., Steidel et al. 2010; Faisst et al. 2016; Matthee et al. 2017, 2019; Carniani et al. 2018b; Verhamme et al. 2018; Pahl et al. 2020; Behrens et al. 2019). A complete analysis of the origins of these offsets is presented in Faisst et al. (2020). The physics of the velocity offsets between [CII] and Lyman α is discussed in Cassata et al. (2020).

7.3. [CII] line width

The [CII] line width is also an interesting physical constraint, since it is linked to the dynamical mass and the size. The constraints on the dynamical masses will be discussed in more details in Dessauges-Zavadsky et al. (2020). Before ALPINE, studies of the [CII] line width have been mainly performed using lensed galaxy samples (e.g., Gullberg et al. 2015). However, these samples are biased toward more star-forming systems, which could also be more massive, and possibly more compact systems. ALPINE is probing less extreme galaxies and we thus expect a narrower line width on average.

In Fig. 17 (lower left panel), we compare our [CII] FWHM distribution with the SPT SMG lensed sample of Gullberg et al. (2015). As expected, our sources have much narrower lines with a median FWHM of 252 ± 13 km s⁻¹ versus 541 ± 110 km s⁻¹ for SPT SMGs. At fixed integrated flux, we could be slightly biased against broader lines. For instance, Kohandel et al. (2019) showed using numerical simulations that edge-on systems tend to be more difficult to detect because of their broader lines, but it is unlikely to be the cause of the large difference between the ALPINE sources and quasars or SMGs. However, the continuum stacking of non detections tends to indicate that they are mainly lower SFR systems (see Sect. 7.5) and thus have probably a low mass and low FWHM. The average FHWM of our sample is also smaller than the 355 ± 18 km s⁻¹ measured by Decarli et al. (2018) in z >  5.94 in quasar hosts, which are also expected to be particularly massive systems.

We detected two sources with particularly narrow lines: CANDELS_GOODSS_42 with a FWHM of 63 km s⁻¹ and vuds_cosmos_510596653 with 62 km s⁻¹. CANDELS_GOODSS_42 is a low-mass object (log(M_⋆/M_⊙) = 9.3 ± 0.3, Faisst et al. 2020). The line is close to the edge of the spectral window, but the signal seems to go back to the baseline level before the very last channel. However, we cannot firmly exclude an edge effect. vuds_cosmos_510596653 is too faint at short wavelength to reliably estimate a stellar mass from ancillary data, which is compatible with a low-mass and thus low-dispersion system.

7.4. [CII] luminosity distribution

The [CII] luminosity L_[CII] of our targets was computed using the formula provided in Carilli & Walter (2013):

$$\begin{array}{c}\hfill L_{[\mathrm{CII}]}=1.04\times {10}^{-3}\times I_{[\mathrm{CII}]}{D}_{\mathrm{L}}^2{\nu }_{\mathrm{obs}}\frac{L_{\odot }}{\mathrm{Jy}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^2\mathrm{GHz}},\end{array}$$(8)

where I_[CII] is the [CII] flux, D_L is the luminosity distance, and ν_obs is the observed frequency. The CMB effect on the [CII] line measurements is expected to be negligible in the ALPINE redshift range (Vallini et al. 2015; Lagache et al. 2018). We found a median luminosity of (4.8 ± 0.4)×10⁸ L_⊙, which is very close to 3.6 × 10⁸ L_⊙ found for the pilot sample of 10 sources of Capak et al. (2015). We can also compare the results with the expected luminosity distributions based on empirical SFR-[CII] relations. We used the relation of De Looze et al. (2014, fit of HII/starburst galaxies) measured in the local Universe and the relation of Lagache et al. (2018) calibrated combining the pre-ALPINE [CII] detections and semi-analytical modeling. The SFRs were taken from the UV to near-IR SED fitting of the ancillary data (Faisst et al. 2020). We did not take into account any scatter in the relations nor the formal uncertainties on SFR from SED fitting for this simple comparison. The results are presented in Fig. 17 (lower right panel). The bulk of our detections are more luminous than the median expected value of the sample. As expected, we thus tend to miss mainly the faintest systems. We observe a significant excess of luminous [CII] objects compared to the distribution expected from the combination of the Lagache et al. (2018) relation and the SFRs based on optical and near-infrared SED fitting. This excess decreases to 1.2σ if we use the De Looze et al. (2014) relation instead.

7.5. Obscured SFR as a function of [CII] luminosity

To understand the excess of luminous [CII] objects in our sample, we derived the average SFR in various [CII] luminosity and redshift bins. These mean SFRs are determined by combining the UV rest-frame data (Faisst et al. 2020) and the mean obscured SFRs determined using a stacking of ALPINE continuum data. This analysis is based only on stacked ALMA data and does not explore the scatter. A more comprehensive study combining ALPINE data and ancillary data is presented in Schaerer et al. (2020).

We first define four [CII] luminosity bins between 10⁸ and 10^9.33 L_⊙ with the same logarithmic width. We chose to stack by bins of [CII] luminosity and not by bins of continuum flux density, since we have many more [CII] detections. In addition, the continuum is not affected by the velocity offsets presented in Sect. 7.2 and can be stacked knowing a priori the exact redshift and line width. The stacking of the undetected [CII] lines would have been more difficult, since we do not know in which frequency range it is located and using a very broad frequency window would lead to poor S/Ns. A few outliers are below or above the chosen [CII] luminosity bins, but the stacking of less than three objects would not provide meaningful results. We also split the sample into two redshift subsamples, below and above z = 5. We then stacked the continuum of all the sources of a given luminosity bin in image space after masking the continuum nontarget detections. We also stacked the ALMA continuum of the [CII] non detections to estimate their mean SFR. Since the beam size and the source size can vary significantly from one source to the other, the resulting stacked profile cannot be well fit by a single elliptical Gaussian. We thus decided to use the aperture method to derive the flux density (see Sect. 3.2). To be consistent, we used the [CII] luminosities derived using the same technique. Finally, we derived the uncertainties using a bootstrap technique, which takes into account both the noise and the intrinsic population variance of the sources (Béthermin et al. 2012b). The full description of the stacking technique and its validation will be described in Khusanova et al. (2020).

To derive the mean SFR of each subsample, we convert the stacked continuum fluxes into L_IR using the method described in Sect. 4 assuming the mean redshift of the subsample. The far-UV luminosity, L_FUV, is provided by the ancillary ALPINE catalog presented in Faisst et al. (2020). The full SFR is derived combining the UV and infrared SFRs (Madau & Dickinson 2014):

$$\begin{array}{c}\hfill {\mathrm{SFR}}_{\mathrm{tot}}={\mathrm{SFR}}_{\mathrm{IR}}+{\mathrm{SFR}}_{\mathrm{UV}}={\kappa }_{\mathrm{IR}}{L}_{\mathrm{IR}}+{\kappa }_{\mathrm{FUV}}{L}_{\mathrm{FUV}},\end{array}$$(9)

where SFR_tot is the total SFR, that is the sum of the dust-obscured component SFR_IR and the unobscured one SFR_UV. κ_IR and κ_FUV are conversion factors from luminosity to SFR and their values are 1.02 × 10⁻¹⁰ M_⊙ yr⁻¹$L_{\odot }^{-1}$ and 1.47 × 10⁻¹⁰ M_⊙ yr⁻¹$L_{\odot }^{-1}$, respectively¹².

Our results are presented in Fig. 18 and in Table 6. In addition to detections, we performed the same stacking procedure on the non detections and computed the mean of the secure upper limits on their luminosity. All our results are within the scatter of the De Looze et al. (2014) and Lagache et al. (2018) relations¹³. At low luminosity (log(L_CII) < 8.66), our mean data points nicely agree with the center of their relation. However, at high [CII] fluxes (log(L_CII) > 8.66), our mean stacked z <  5 data point is on the 1-σ upper envelop of these relations. There is thus a mild average [CII] excess in the most luminous sources of our sample. Even if it is rare, [CII]-excesses have been observed in the local Universe (Smirnova-Pinchukova et al. 2019), but we usually observe a deficit in bright sources (Díaz-Santos et al. 2013). This apparent [CII] excess could also be an SED effect. Indeed, if brighter sources have warmer dust than the average stacked value, their obscured SFR would be underestimated. This explanation is compatible with the warmer dust SED measured by stacking for the SFR > 100 M_⊙ yr⁻¹ sources at the same redshift in the COSMOS field (see Sect. 4.4). If we use the Schreiber et al. (2018) SED template, which fits better the high-SFR stacked SED, instead of the Bethermin et al. (2017) one, we find that SFR_IR is a factor of 1.4 higher. In Fig. 18, we illustrate the impact of a warmer dust using a black arrow, showing that it can explain this small [CII] excess. Finally, we tested if this excess could be caused by merging systems by excluding objects classified as merger (type = 2) in Le Fèvre et al. (2020) and redoing the full procedure, but the excess remained. These results and a detailed analysis of the [CII]–SFR relation at high redshift are discussed in depth in Schaerer et al. (2020).

Fig. 18. Mean [CII] luminosity as a function of SFR. The obscured SFR is measured by stacking the ALPINE continuum maps for various subsamples selected in redshift and [CII] luminosity. It is combined with the mean uncorrected UV SFR from ancillary data (Faisst et al. 2020). The method is described in Sect. 7.5. The blue squares are our measurements at z <  5 and the red diamonds are for z >  5. The y-axis upper limits correspond to the mean total SFR of nondetected [CII] samples (the y position corresponds to mean of the secure [CII] upper limits). The blue and red dashed lines indicate the relation of Lagache et al. (2018) at z = 4.5 and z = 5.5, respectively. The gray shaded area is the 1-σ region of the relation of De Looze et al. (2014). The black arrow represents the shift in SFR if we use the warmer Schreiber et al. (2018) SED templates instead of the Bethermin et al. (2017) ones to determine the infrared luminosity.

8. Conclusion

In this paper, we presented the data processing and the catalog construction of the ALPINE ALMA large program. The performance of our survey is fully compatible with our initial goal. Our main technical results are:

• After flagging a few badly calibrated antennae, we produced continuum and [CII] moment-0 maps and reached a typical sensitivity of 30 μJy rms in continuum and 0.14 Jy km s⁻¹ for [CII] in a 235 km s⁻¹ bandwidth. The average beam size is 1.13″ × 0.85″.

• We investigated the stability of quasars used as flux calibrators during the observation campaign and found that the fluctuations between two successive calibrator survey observations are lower than the standard 10% calibration uncertainty.

• We detected 23 of our 118 targets in continuum above 3.5σ threshold corresponding to a 95% purity. In the rest of the field, we had to extract sources above 5σ because of the lower number density of emitters and detected 57 nontarget additional continuum sources.

• We measured the flux densities of our detections using five different methods. They well agree with each other, except the one based on peak flux density because of the large fraction of marginally resolved objects in our sample.

• We performed Monte Carlo simulations to estimate the completeness and the flux boosting for the nontarget continuum sources and proposed a simple way to derive them as a function of the source flux density, size, and the local noise.

• After adjusting our extraction algorithm to reach 95% purity, we detected 75 of our 118 targets in [CII]. Similarly to the continuum, we checked the robustness of the photometry by comparing five different methods.

In addition to these technical results, we obtained promising first scientific results. Our main findings are:

• To determine the conversion factor from the 158 μm continuum luminosity to the total infrared luminosity (L_158 μm = 0.093 L_IR), we measured average dust SEDs by stacking single-dish data at the position of COSMOS photometric sources similar to ALPINE ones.

• The target sources detected in continuum have a median flux density of 0.26 mJy, a median L_IR of 4.4 × 10¹¹ L_⊙, and a median stellar mass of 1.1 × 10¹⁰ M_⊙. As expected, the detections are among the most massive and star-forming systems of our sample.

• The nontarget continuum detections have a mean redshift of 2.5 ± 0.2 (median = 2.3) and are mainly lower-redshift sources without a physical connection with the ALPINE targets.

• We derived number counts probing the 0.35−5.6 mJy range. We identified a slope break in the counts around 3 mJy and estimated that the contribution of > 0.35 mJy sources to the CIB is 40–69%.

• The detected [CII] targets have a median [CII] luminosity of 4.8 × 10⁸ L_⊙ and a median FWHM of 252 km s⁻¹. We also observed significant offsets between the optical and [CII] redshifts (up to 1000 km s⁻¹).

• We measured the mean obscured SFR in various [CII] luminosity bins by stacking ALPINE continuum data and combined it with the unobscured SFR from ancillary rest-frame UV data. Our data agrees with the local and predicted SFR–L_[CII] relations of De Looze et al. (2014) and Lagache et al. (2018).

A series of companion papers discuss the other scientific results of ALPINE (e.g., Le Fèvre et al. 2020; Jones et al. 2020; Ginolfi et al. 2020b,a; Faisst et al. 2019; Cassata et al. 2020; Schaerer et al. 2020; Fudamoto et al. 2020; Fujimoto et al. 2020; Dessauges-Zavadsky et al. 2020; Yan et al. 2020; Loiacono et al. 2020; Gruppioni et al. 2020; Romano et al. 2020)¹⁴.

Acknowledgments

We thanks the ALMA observatory staff for their support and in particular Edwige Chapillon, our ALMA contact scientist, for her very useful advice. MB acknowledges Mladan Novak for his excellent advice about the photometry in interferometric surveys. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2017.1.00428.L. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. This research made use of Astropy (http://www.astropy.org), a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018). This program receives funding from the CNRS national program Cosmology and Galaxies. A.C., C.G., F.L., F.P., and M.T. acknowledge the support from grant PRINMIUR 2017 20173ML3WW_001. DR acknowledges support from the National Science Foundation under grant numbers AST-1614213 and AST-1910107 and from the Alexander von Humboldt Foundation through a Humboldt Research Fellowship for Experienced Researchers. G.C.J. acknowledges ERC Advanced Grant 695671 “QUENCH” and support by the Science and Technology Facilities Council (STFC). E.I. acknowledges partial support from FONDECYT through grant N° 1171710. GL acknowledges support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (project CONCERTO, grant agreement No 788212) and from the Excellence Initiative of Aix-Marseille University-A*Midex, a French “Investissements d’Avenir” programme. This paper is dedicated to the memory of Olivier Le Fèvre, PI of the ALPINE survey.

References

Appendix A: Pointing list, beam size, and depth

Appendix B: Continuum cutouts and catalogs

Fig. B.1. Cutout images of the continuum-detected target sources. The red and the gold contours indicate the 3.5σ and 5σ levels, respectively. The white crosshair indicate the phase center.

Fig. B.2. Cutout images of the continuum-detected nontarget sources. The red and the gold contours indicate the 3.5σ and 5σ levels, respectively. The white crosshair indicate the phase center. The sources are ordered by decreasing S/N.

Fig. B.2. continued.

Appendix C: [CII] spectra, moment-0 maps, and catalogs

Fig. C.1. Spectra of the detected ALPINE sources (S/N >  3.5). The blue area is the frequency range used to produce the moment-0 map. The red solid line is the best-fit by a Gaussian. The process used to produce these spectra is described in Sect. 6.1.

Fig. C.1. continued.

Fig. C.1. continued.

Fig. C.2. Cutout images of the moment-0 maps of our detected [CII] targets. The red and the gold contours indicate the 3.5σ and 5σ levels, respectively. The white crosshair indicate the phase center. The sources are ordered by decreasing S/N.

Fig. C.2. continued.

Fig. C.2. continued.

Appendix D: Sources with complex photometry

D.1. Continuum complex sources

For most of our sources, the automatic photometric measurements described in Sect. 3.2 are sufficient. However, a manual deblending is required for some sources with one or several neighbors. The most complex object is DEIMOS_COSMOS_881725, which was fit using three components simultaneously. These three components are shown in Fig. D.1. Three other continuum sources have a close neighbor, which disturbs our automatic photometric procedure. They were also measured manually. All the flux densities are listed in Table D.1. As pointed out in Le Fèvre et al. (2020), about 40% of ALPINE targets are classified as merging systems (see the criteria in Le Fèvre et al.). For these systems, some deblending will be needed to properly distribute the observed line or continuum emission between the different components of the merger. This deblending of components close to the scale of the synthesized beam based on optical priors will be presented in a future paper.

Fig. D.1. Cutout of the continuum image of DEIMOS_COSMOS_881725. The three components are indicated and the figure and their flux densities are listed in Table D.1.

D.2. [CII] complex source

We had less problems with the deblending of [CII] target sources. Indeed, since the flux is measured in the moment-0 maps, it is very unlikely to find a CO interloper at the same exact frequency. The only possible contaminant is another physically related object at the same velocity as the ALPINE target. For instance, DEIMOS_COSMOS_842313 has an extremely bright neighbor known as J1000+0234 (Schinnerer et al. 2008) or AzTEC/C17 (Brisbin et al. 2017). This nearby submillimeter galaxy is identified as SC_2_DEIMOS_COSMOS_842313 in our nontarget catalog, but there is no signal at this position on the moment-0 map, since the [CII] of this source is at a different velocity.

In the case of close mergers, our three photometric methods are compatible and tend to measure the total flux of the system. However, vuds_cosmos_5101209780 is a bit more problematic. The ALPINE target is faint, but it has a bright southern companion also detected in the moment-0 map. Some photometric methods include this component while other exclude it. We thus decided to perform a manual extraction after having defined manually two regions corresponding to each component (see Fig. D.2). The fluxes are provided in Table D.2.

Fig. D.2. Cutout of the moment-0 map of vuds_cosmos_5101209780 indicating how we defined the two components.

All Tables

All Figures

Fig. 1. Upper panel: 345 GHz flux density of the flux calibrators used by the ALPINE survey (the J1058+0133, J0854+2006, and J0522−3627 quasars) as a function of time. The gray areas indicate when the ALPINE targets were observed (see Sect. 2.1). The only quasar used for CDFS targets (J0522−3627) is plotted with a dashed line, while solid lines are used for the calibrators of COSMOS sources. Lower panel: spectral index versus time. This spectral index is estimated using the band-7 and band-3 flux from the calibrator monitoring performed by the observatory (see Sect. 2.3).

In the text

Fig. 2. Achieved [CII] (upper panel) and continuum (lower panel) rms sensitivities. The blue dots indicate the values measured in individual fields and the red squares the mean values in the two redshift windows. The red error bars on the plots are the standard deviation in each redshift range. The actual uncertainties on the mean values are indeed $\sqrt{N}$ times smaller (central-limit theorem) and are smaller than the size of the squares. Since the [CII] sensitivities were measured using different bandwidths because of the different line widths, we normalized the measurements to a bandwidth of 235 km s−1 by dividing our raw measurements by $\sqrt{\mathrm{\Delta }\mathit{v}/(235\mathrm{km}{\mathrm{s}}^{-1})}$. The solid black lines indicate the trend of the [CII] flux I[CII] and continuum flux S(1 + z)158 μm versus frequency (and thus redshift) at constant [CII] luminosity L[CII] and fixed infrared luminosity LIR, respectively.

In the text

	Fig. 3. Purity as a function of the S/N threshold. The results obtained around the center of the pointings (1 arcsec radius) are in blue. The results in the full field are in red. The dotted lines show the S/N at which the 95% is reached.
In the text

Fig. 4. Comparison between our various photometric methods described in Sect. 3.2 for S/N >  5 sources. The blue dots are our measurements and the red line is the one-to-one relation. Upper left, upper right, and lower left panels: comparison between the 2D-fit flux density (x-axis) and the aperture, 2σ-clipped, and peak flux densities, respectively. Lower right panel: ratio between the peak flux and the 2D-fit flux as a function of the ratio between the source area (convolved by the synthesized beam) and the synthesized beam area. The dashed line indicates the expected trend (see Sect. 3.5).

In the text

	Fig. 5. Comparison between the 2D-fit flux densities derived in map space (Sect. 3.2) and the flux determined fitting an elliptical Gaussian model in the uv plane (Sect. 3.6). The blue dots are our measurements and the red line is the one-to-one relation.
In the text

Fig. 6. Upper panel: completeness as a function of the continuum nonprimary-beam-corrected flux density (Gpb Sinj) achieved for point sources in various pointings. Middle panel: similar figure after having divided the nonprimary-beam-corrected flux by the noise at the center of each pointing (σcenter). Various colors (blue, green, yellow, and red) corresponds to various injected source sizes (FWHM = 0, 0.333, 0.666, and 1 arcsec, respectively). The solid lines indicate the mean trend of the various pointings, while the dashed lines indicate the 1-σ envelop. Lower panel: same plot after normalizing the injected flux by 1/σcenter and by the source area (Ωsource/Ωbeam). These results are discussed in Sect. 3.8.

In the text

Fig. 7. Ratio between the injected and recovered flux density as a function of measured S/N (see Sect. 3.9). Upper left, upper right, lower left, and lower right panels: results obtained for the 2D-fit, aperture, 2 σ-clipped, and peak photometry, respectively. The solid lines indicate the median and the shaded areas are the 1-σ contours. Various colors (blue, green, yellow, and red) are used to indicate the various sizes used (FWHM = 0, 0.333, 0.666, and 1 arcsec, respectively). The dashed horizontal line indicate the one-to-one relation.

In the text

	Fig. 8. Effective surface area of ALPINE as a function of the 850 μm flux after excluding the central 1 arcsec-radius area where target sources are extracted. The blue, green, gold, and red lines are the results obtained for a source size of 0, 0.3, 0.6, and 1 arcsec, respectively. The method used to compute the surface area is described in Sect. 3.10.
In the text

Fig. 9. Comparison between the Álvarez-Márquez et al. (2016, blue dashed line), Schreiber et al. (2018, orange dot-dashed line), and Bethermin et al. (2017, red solid line) IR SED templates and the observed mean SEDs of SFR > 10 M⊙ yr−1 galaxies measured by stacking (black dots, see Sect. 4.2). The black dotted line is the best fit of the λrest − frame >  40 μm data points by a modified blackbody with β fixed to 1.8 (the temperature in the legend is provided in the rest frame). Upper and lower panels: 4 <  z <  5 and 5 <  z <  6, respectively.

In the text

	Fig. 10. Same figure as Fig. 9 (upper panel), but using a SFR > 100 M⊙ yr−1 cut.
In the text

Fig. 11. Upper left panel: redshift distribution of the ALPINE target sources. The blue and red histograms are the distribution of detected sources only and the full sample, respectively. Upper right panel: distribution of the continuum flux densities. The red histogram indicates the distribution expected from the optical and near-IR SED-derived SFR assuming the long-wavelength SED presented in Sect. 4. Lower left panel: same figure as previously but for LIR. Lower right panel: stellar mass distribution of detected (blue) and all (red) sources.

In the text

	Fig. 12. Photometric redshift distribution of the nontarget ALPINE sources (blue filled histogram). The dashed, solid, and dotted histograms are the predictions from the Bethermin et al. (2017) SIDES simulation for a flux cut corresponding to the first-quartile, the median, and the third-quartile of the observed sample, respectively.
In the text

Fig. 13. Left panel: differential number counts at 850 μm. The blue circles are the counts derived from the full ALPINE nontarget sample. The red squares are computed from the ALPINE nontarget sources with a confirmed optical or near-IR counterpart below z = 4 (secure z <  4 sample). It is thus a secure lower limit, since sources clustered with ALPINE sources are excluded. The green crosses and pentagons are the interferometric measurements of Karim et al. (2013) and Stach et al. (2018), respectively. The gray left-facing triangles, right-facing triangles, diamonds, and down-facing triangles are the measurements performed using single-dish data with a lower angular resolution by Geach et al. (2017), Hsu et al. (2016), Casey et al. (2013), and Chen et al. (2013), respectively. The solid, dashed, dotted, and dash-dot lines are the models of Bethermin et al. (2017, SIDES simulation), Gruppioni et al. (2011), Casey et al. (2018, low dust), and Casey et al. (2018, high dust). Right panel: integral number counts at 850 μm. The symbols are the same as in the left panel. The interferometric measurements of Oteo et al. (2016) and Simpson et al. (2015), respectively, are represented by green stars and triangles.

In the text

Fig. 14. Contribution of galaxies to the cosmic infrared background at 850 μm as a function of the flux density limit. The blue and red solid lines are our results based on the full and the secure z <  4 samples, respectively. The shaded area is the 1-σ confidence region. The dashed green line is the measurement of Hsu et al. (2016) using SCUBA2 data of cluster fields. The black line is the absolute measurement of the total CIB combining COBE/FIRAS and Planck data (Odegard et al. 2019).

In the text

Fig. 15. Purity of [CII] detections as a function of the S/N threshold in the moment-0 map. The purity is provided only for the 1 arcsec central region and was estimated using the methods described in Sect. 6.2. The blue solid line is derived using the method based on the moment-0 maps and the red dashed line using an extraction in the negative cubes (see Sect. 6.2). The dotted vertical lines are the S/N at which the 95% purity is reached.

In the text

Fig. 16. Comparison between our various [CII] line flux measurement methods described in Sect. 6.3. The blue dots are our measurements and the red line is the one-to-one relation. Upper left, upper right, and lower left panels: comparison between the 2D-fit flux (x-axis) and the aperture, 2-σ-clipped, and peak flux densities, respectively (y-axis). The vuds_cosmos_5101209780 object (in green) has a complex geometry and is discussed in Appendix D.2. The lower right panel shows the comparison between uv-plane fluxes and 2D-fit fluxes for single-component systems.

In the text

Fig. 17. Upper left panel: redshift distribution of the full sample (red open histogram) and the [CII] detections (blue filled histogram). Upper right panel: distribution of the velocity offset between the [CII] and the reference optical redshift. Lower left panel: distribution of the [CII] line FWHM of the detected ALPINE sources (blue filled histogram) and comparison with the lensed galaxy sample of Gullberg et al. (2015, red) and the quasar sample of Decarli et al. (2018, green). Lower right panel: comparison between the [CII] luminosity distribution of the detections (blue filled histogram) and the expected values based on the SFRs determined by optical and near-infrared SED fitting and the De Looze et al. (2014, red open histogram) and Lagache et al. (2018, black open histogram) relations.

In the text

Fig. 18. Mean [CII] luminosity as a function of SFR. The obscured SFR is measured by stacking the ALPINE continuum maps for various subsamples selected in redshift and [CII] luminosity. It is combined with the mean uncorrected UV SFR from ancillary data (Faisst et al. 2020). The method is described in Sect. 7.5. The blue squares are our measurements at z <  5 and the red diamonds are for z >  5. The y-axis upper limits correspond to the mean total SFR of nondetected [CII] samples (the y position corresponds to mean of the secure [CII] upper limits). The blue and red dashed lines indicate the relation of Lagache et al. (2018) at z = 4.5 and z = 5.5, respectively. The gray shaded area is the 1-σ region of the relation of De Looze et al. (2014). The black arrow represents the shift in SFR if we use the warmer Schreiber et al. (2018) SED templates instead of the Bethermin et al. (2017) ones to determine the infrared luminosity.

In the text

	Fig. B.1. Cutout images of the continuum-detected target sources. The red and the gold contours indicate the 3.5σ and 5σ levels, respectively. The white crosshair indicate the phase center.
In the text

	Fig. B.2. Cutout images of the continuum-detected nontarget sources. The red and the gold contours indicate the 3.5σ and 5σ levels, respectively. The white crosshair indicate the phase center. The sources are ordered by decreasing S/N.
In the text

	Fig. B.2. continued.
In the text

	Fig. C.1. Spectra of the detected ALPINE sources (S/N >  3.5). The blue area is the frequency range used to produce the moment-0 map. The red solid line is the best-fit by a Gaussian. The process used to produce these spectra is described in Sect. 6.1.
In the text

	Fig. C.1. continued.
In the text

	Fig. C.1. continued.
In the text

	Fig. C.2. Cutout images of the moment-0 maps of our detected [CII] targets. The red and the gold contours indicate the 3.5σ and 5σ levels, respectively. The white crosshair indicate the phase center. The sources are ordered by decreasing S/N.
In the text

	Fig. C.2. continued.
In the text

	Fig. C.2. continued.
In the text

	Fig. D.1. Cutout of the continuum image of DEIMOS_COSMOS_881725. The three components are indicated and the figure and their flux densities are listed in Table D.1.
In the text

	Fig. D.2. Cutout of the moment-0 map of vuds_cosmos_5101209780 indicating how we defined the two components.
In the text