The ALPINE-ALMA [CII] Survey: data processing, catalogs, and statistical source properties

The ALPINE-ALMA large program targets the [CII] 158 $\mu$m line and the far-infrared continuum in 118 spectroscopically confirmed star-forming galaxies between z=4.4 and z=5.9. It represents the first large [CII] statistical sample built in this redshift range. We present details of the data processing and the construction of the catalogs. We detected 23 of our targets in the continuum. To derive accurate infrared luminosities and obscured star formation rates, we measured the conversion factor from the ALMA 158 $\mu$m rest-frame dust continuum luminosity to the total infrared luminosity (L$_{\rm IR}$) after constraining the dust spectral energy distribution by stacking a photometric sample similar to ALPINE in ancillary single-dish far-infrared data. We found that our continuum detections have a median L$_{\rm IR}$ of 4.4$\times 10^{11}$ L$_\odot$. We also detected 57 additional continuum sources in our ALMA pointings. They are at lower redshift than the ALPINE targets, with a mean photometric redshift of 2.5$\pm$0.2. We measured the 850 $\mu$m number counts between 0.35 and 3.5 mJy, improving the current interferometric constraints in this flux density range. We found a slope break in the number counts around 3 mJy with a shallower slope below this value. More than 40 % of the cosmic infrared background is emitted by sources brighter than 0.35 mJy. Finally, we detected the [CII] line in 75 of our targets. Their median [CII] luminosity is 4.8$\times$10$^8$ L$_\odot$ and their median full width at half maximum is 252 km/s. After measuring the mean obscured SFR in various [CII] luminosity bins by stacking ALPINE continuum data, we find a good agreement between our data and the local and predicted SFR-L$_{\rm [CII]}$ relations of De Looze et al. (2014) and Lagache et al. (2018).


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).
The catalogs are also available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsarc. u-strasbg.fr/viz-bin/cat/J/A+A/643/A2 The ALPINE products are publicly available at https://cesam. lam.fr/a2c2s/ Targeted observations of known sources from optical and nearinfrared 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)

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 −1 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 −1 (or FWHM ∼ 235 km s −1 ). 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 −1 ) 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 8 L assuming a line width of 235 km s −1 . 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 −1 to obtain the same luminosity (∼0.2 mJy beam −1 in 235 km s −1 band at z = 5.8 versus ∼0.3 mJy beam −1 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.

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 −1 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 A2, page 2 of 43 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.

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 1 . 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.
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 1 https://almascience.eso.org/sc/ 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.

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 A2, page 3 of 43  √ 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 ∆v/(235 km 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 L IR , respectively.
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, y) = N channel k=1 S ν (x, y, k) ∆v 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.

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 linecontaminated 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 naturalweighted 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.

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 −1 . 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 −1 than in the high-frequency range (345−356 GHz, 4.3 < z < 4.6) with 0.17 ± 0.04 Jy km s −1 . 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 A2, page 4 of 43 [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.
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 −1 . We obtained a better sensitivity for 5.1 < z < 5.9 sources (283−315 GHz) with an average value of 28 µJy beam −1 . 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

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 nonprimarybeam-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: 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σ. 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

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 marginallyresolved 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: 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.

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.

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 A2, page 6 of 43 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.

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).
In order to check how consistent are our measurements, we computed the uncertainty-normalized difference between two measurements S method A ν and S method B ν : 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: 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.

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 2 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 3 . We could successfully model nine continuum targets 4 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).

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  Map versus uv-plane continuum photometry 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.

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

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

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 , 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 850 and the source size θ source is derived from the completeness C(S 850 , θ source , x, y) at a position (x, y) (see Sect. 3.8) using: 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 5 (24.92 arcmin 2 ). 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 Aravena et al. 2016;González-López et al. 2017;Franco et al. 2018).

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 −10 M yr −1 /L × L IR , Kennicutt 1998after converting to Chabrier 2003. 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), Lymanbreak selections (e.g., Álvarez-Márquez et al. 2016), and lowredshift 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.

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 groundbased measurements in the COSMOS field to derive the mean SEDs of z < 4 galaxies. We used the same The 2015 selection of the stacked targets was performed using a stellar mass cut of >3 × 10 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 −1 , which is approximately equivalent to the ALPINE SFR limit 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.

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 7 . 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 (χ 2 = 0.44 for 3 degrees of freedoms) and acceptable at 5 < z < 6 (χ 2 = 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 8 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 10 3 Observed wavelength [ m] . 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.
per unit of redshift). Because of the nature of the ALPINE sample , 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 (χ 2 < 4 with 4 degrees of freedom for both templates in both redshift bins). Since the χ 2 of Bethermin et al. (2017) is marginally better, we 0.08 ± 0.04 0.10 ± 0.07 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.

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 −1 . According to optical and near-infrared SED fitting , 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.

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). Table 3. Ratio (without unit) between the monochromatic continuum luminosity νL ν and the total infrared luminosity L IR at different rest-frame wavelengths associated with important far-IR lines.

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). 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 −10 L /(M yr −1 ) 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 −5 ) 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).

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 COS-MOS (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  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. , it is expected that <0.1 mJy 1.1 mm sources are at lower redshift than ∼1 mJy 850 µm sources. Finally, we compare our measured distribution with the predictions of the simulated infrared dusty infrared sky (SIDES) simulation 9 (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.

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 Table 4. Integral number counts at 850 µm derived from ALPINE for the full nontarget sample (all) and the nontarget sources with an optical counterpart at z < 4 (secure, z < 4).  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 10 : 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 µmconverted 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.
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 lowresolution 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)

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  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: 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: dN/dS ν ∝ S −2.2±0.3 ν for the full sample and dN/dS ν ∝ S −1.8±0.4 ν for the secure sample. Above our brightest data point, we use the Euclidian plateau level measured by Planck Collaboration VII (2013, dN/dS S 2.5 = 15 Jy 1.5 sr −1 ). 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 −1 for the full sample and 0.054 ± 0.009 MJy sr −1 for the secure sample. This is 69% and 40% of the full CIB measured by Odegard et al. (2019, 0.135 MJy sr −1 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).

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

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 −1 ) to avoid single-channel spikes of noise, and the maximum is 40 channels (1000 km s −1 ). 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.

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 −1 , 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). Since this normalization is an approximation, we crosschecked our number with an independent approach. We extracted the sources in the negative and the positive cubes using the findclumps 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 −1 and 216 ± 1 km s −1 , respectively 11 . 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 11 The error bars on the mean FWHM are the uncertainty on the mean of the sample and not the scatter (see caption of Fig. 2). The spurious sources are extracted from the full cubes and not the 1 arcsec-radius central regions. They are thus much more numerous. The uncertainty on their mean FWHM is thus smaller than the detections in the central region. 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).

[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 highvelocity 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).
We also compared our flux measurements performed in the map space with uv-derived values. These measurements were performed using the uvmodel f it 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 −1 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.

[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 A2, page 18 of 43 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).

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

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 −1 . 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 −1 . The optical redshifts could suffer from several effects leading to a small inaccuracy.

[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 −1 versus 541 ± 110 km s −1 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 −1 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 −1 and vuds_cosmos_510596653 with 62 km s −1 . CAN-DELS_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.

[CII] luminosity distribution
The [CII] luminosity L [CII] of our targets was computed using the formula provided in Carilli & Walter (2013): 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 8 L , which is very close to 3.6 × 10 8 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- The SFRs were taken from the UV to near-IR SED fitting of the ancillary data ). 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 A2, page 20 of 43  ). 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.

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 ) 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 8 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): where SFR tot is the total SFR, that is the sum of the dustobscured 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 −10 M yr −1 L −1 and 1.47 × 10 −10 M yr −1 L −1 , respectively 12 .
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.  (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 −1 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  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).

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: 12 The κ coefficients from Madau & Dickinson (2014) have been converted from Salpeter to Chabrier IMF by dividing them by a factor of 1.7. 13 The scatter on the Lagache et al. (2018) is not shown on the figure for clarity, but it is larger than De Looze et al. (2014) (∼0.55 dex). -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 −1 for [CII] in a 235 km s −1 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 11 L , and a median stellar mass of 1.1×10 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%.  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) 14 .     Notes. The frequency is the mean of the channels used to derive the continuum. The flux S ν was measured using the 2D-fit approached described in Sect. 3.1. The monochromatic luminosity νL ν is derived at (1 + z) ν obs and does not need any assumption on the galaxy SEDs. The total infrared luminosity L IR (8−1000 µm) is derived assuming Bethermin et al. (2017) SED (see discussion in Sect. 4). The measurement uncertainties do not include the calibration uncertainties (∼4.5% on average and up to 10%, see Sect. 2.3). Notes. Three methods are used to compute them (aggressive, normal, and secure) and are described in Sect. 3.3. In the absence of any external information, we recommend to use the secure upper limits. If the source is known to be point-like, the normal upper limits can be used. The luminosities are derived using the same method as in Table B.1.

Appendix A: Pointing list, beam size, and depth
A2, page 29 of 43 A&A 643, A2 (2020)    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.      Notes. The three methods used to compute them (aggressive, normal, and secure) are the same as in Table B.2 (method described in Sect. 3.3), but applied to the moment-0 maps. In the absence of any external information, we recommend to use the secure upper limits. If the source is known to be point-like, the normal upper limits can be used.  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  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.

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.