A$^3$COSMOS: the infrared luminosity function and dust-obscured star formation rate density at $0.5

Aims: We leverage the largest available Atacama Large Millimetre/submillimetre Array (ALMA) survey from the archive (A$^3$COSMOS) to study to study infrared luminosity function and dust-obscured star formation rate density of sub-millimeter/millimeter (sub-mm/mm) galaxies from $z=0.5\,-\,6$. Methods: The A$^3$COSMOS survey utilizes all publicly available ALMA data in the COSMOS field, therefore having inhomogeneous coverage in terms of observing wavelength and depth. In order to derive the luminosity functions and star formation rate densities, we apply a newly developed method that corrects the statistics of an inhomogeously sampled survey of individual pointings to those representing an unbiased blind survey. Results: We find our sample to mostly consist of massive ($M_{\star} \sim 10^{10} - 10^{12}$ $\rm M_{\odot}$), IR-bright ($L_* \sim 10^{11}-10^{13.5} \rm L_{\odot}$), highly star-forming (SFR $\sim 100-1000$ $\rm M_{\odot}$ $\rm yr^{-1}$) galaxies. We find an evolutionary trend in the typical density ($\Phi^*$) and luminosity ($L^*$) of the galaxy population, which decrease and increase with redshift, respectively. Our IR LF is in agreement with previous literature results and we are able to extend to high redshift ($z>3$) the constraints on the knee and bright-end of the LF, derived by using the Herschel data. Finally, we obtain the SFRD up to $z\sim 6$ by integrating the IR LF, finding a broad peak from $z \sim 1$ to $z \sim 3$ and a decline towards higher redshifts, in agreement with recent IR/mm-based studies, within the uncertainties, thus implying the presence of larger quantities of dust than what is expected by optical/UV studies.


Introduction
Understanding how galaxies formed and evolved is one of the open questions of modern astrophysics.This can be addressed in many different ways, using information across the whole electromagnetic spectrum.One of the best approaches consists in the study of galaxy samples across a wide range of redshift and luminosity, enabling the derivation of statistical properties as a function of the redshift, such as the luminosity function (LF) and the cosmic star formation rate density (cSFRD).These quantities are fundamental to probe the statistical nature of various galaxy populations at different cosmic times, as well as to study the mass assembly process in galaxies at different epochs.
Optical/UV based studies have made it possible to compute the SFRD up to z ∼ 7 − 8 (see e.g., Bouwens et al. 2014;Oesch et al. 2015;Laporte et al. 2016;Oesch et al. 2018) and even z ∼ 10 (e.g., Harikane et al. 2023) with the James Webb Space Telescope (JWST), extending our knowledge on the star formation activity to very early epochs of the universe.However, these studies can be potentially biased by the observing band, i.e., the rest-frame UV, which is highly affected by dust obscuration.Indeed, the dusty contribution is only retrieved from dust correction measured in the UV.These corrections are still very uncertain, as thermalization and re-emission by dust in the IRmm bands has been shown to be significantly contributing to the SFRD even at z > 3 − 4 (see e.g., Magdis et al. 2012;Magnelli et al. 2013;Béthermin et al. 2015;Casey et al. 2018;Zavala et al. 2018;Gruppioni et al. 2020;Algera et al. 2022).Therefore, it is crucial to determine the contribution of galaxies selected in the IR band, where this emission is produced.While attempts have been made to constrain the SFRD IR at z > 3 in the past using single-dish IR-mm surveys, only the advent of ALMA opened up this possibility.The unprecedented sensitivity reached by ALMA, coupled with the assembly of unbiased samples in the mm bands (Hodge et al. 2013;Staguhn et al. 2014;Zavala et al. 2018;Franco et al. 2018a), allows for the study of evolution of these galaxies up to higher redshifts, thus covering the z > 3 range still affected by many uncertainties (poor statistics, bias in the IR luminosity).Recent works using sub-mm/mm samples (e.g., Gruppioni et al. 2020;Algera et al. 2022) support a scenario in which the SFRD shows a plateau rather than a significant decrease at z = 2 − 6.These first studies are, however, limited by statistics and larger samples are required to better constrain the SFRD at higher redshift.
In this perspective, the Automated mining of the ALMA Archive in COSMOS (A 3 COSMOS, Liu et al. 2019a,b), which is a compilation of all ALMA observations in the COSMOS field, represents the largest ALMA survey to date.The fact that the survey is in the COSMOS field (Scoville et al. 2007), where a large wealth of multiwavelength data are available, including the COSMOS2020 catalog (Weaver et al. 2022), makes it ideal to perform statistical studies on the nature and evolution of star forming galaxies over a large range of redshifts and luminosities.However, it is not a purely blind survey, which instead is needed to perform statistical studies.Moreover, the individual pointings composing the survey are at different observing wavelengths and have different sensitivities, making it even more inhomogeneous.For these reasons, we developed, within the A 3 COSMOS collaboration, a new method specifically tailored to turn a targeted survey, composed by an arbitrary number of pointings (isolated or overlapping), each with its own limiting flux (radially varying within the pointing) and observing band, into a "blind-like" (targeted unbiased) one, thus allowing the derivation of the main statistical properties of large galaxy samples over cosmic time.However, it is important to bear in mind that the conversion into a blind survey can be affected by uncertainties related to the assumptions made (see Section 4), mainly on the removal of the pointed target and/or clustered sources as well as on the conversion curves for the RMS.In order to take advantage of the most recent A 3 COSMOS 1 and multiwavelength (COSMOS2020 Weaver et al. 2022) catalogs, we have performed a new catalog match and SED-fitting analysis using the python "Code Investigating GALaxy Emission" (CIGALE; Boquien et al. 2019) SED fitting tool.We then derived the IR LF and we present new estimates for the dust-obscured SFRD from z ∼ 0.5 to z ∼ 6.The paper is organized as follows: in Section 2 we present in detail the A 3 COSMOS survey and the multiwavelength sample, as well as the available photometry; in Section 3 we describe the SED fitting and the results obtained for the physical properties; in Section 4 we illustrate our new method to turn the A 3 COSMOS survey from a pointed into a blind one; in Section 5 we present our estimates of the IR LF; in Section 6 we derive the SFRD; finally, in Section 7 we report a summary of our conclusions.

Catalog descriptions
In this section we present the catalogs used in this work and describe their main features.We provide a brief description of the automated pipeline used to obtain the catalogue.A more detailed description of how the pipeline was developed can be found in Liu et al. (2019b) and Adscheid et al., prep; here we briefly describe it.The survey is built by downloading the ALMA pointings from the archive with automatic pipelines, using the Python package astroquery (Ginsburg et al. 2019).For calibration and creation of continuum images from the raw data, the Common Astronomy Software Application (CASA; McMullin et al. 2007) was used.Sources were extracted blindly and in prior mode, using Python Blob Detector and Source Finder (PyBDSF; Mohan & Rafferty 2015, blind) and GALFIT (Peng et al. 2002(Peng et al. , 2010, prior), prior), utilizing prior source positions from multi-wavelength catalogs covering the COSMOS field.The matching with the priors was done with a 1" radius to reduce multiple associations.Finally, to limit the number of spurious sources, a minimum peak signal-to-noise threshold of 4.35 σ (for the prior) was applied, resulting in a global spurious fraction of less than 12% (see Figure 8 from Liu et al. 2019b).

A 3 COSMOS catalog
The large amount of observation collected by the ALMA interferometer has already been explored in some recent works (see, e.g., da Cunha et al. 2015;Bouwens et al. 2016;Fujimoto et al. 2017;Scoville et al. 2017;Zavala et al. 2018), focusing on the physical properties of the high-z galaxy population.However, in order to investigate galaxy properties (such as the gas and dust content) of statistically significant samples, a systematic mining of the archive is needed.In this context, the A 3 COSMOS project (Liu et al. 2019b,a) aims at building a catalog of galaxies from ALMA archival images, processed homogeneously, in the COSMOS field.In this way, it is possible to retrieve both targeted and serendipitously detected objects in the individual pointings, building a statistically robust catalog of sub-mm detected sources.Within the A 3 COSMOS survey, two different catalogs are available.The first contains sources blindly extracted from the images, while the second one is a prior-based catalog, using optical/near-infrared positional priors (see Section 2.2).In this work, we used the prior version of the catalog (1620 sources), since it allows us to construct the broad-band (from UV up to sub-mm/mm) spectral energy distribution (SED) of our sample.
The most recent release of the COSMOS photometric catalog (i.e.COSMOS2020, Weaver et al. 2022) is characterized by the addition of new data from the Hyper Suprime-Cam (HSC) Subaru Strategic Program (SSP) PDR2 (Aihara et al. 2019), new data from the DR4 (Moneti et al. 2023) of the Visible Infrared Survey Telescope for Astronomy (VISTA) and all the Spitzer IRAC data in the COSMOS field.Moreover, it contains two independently derived photometric datasets.One (the CLASSIC) retrieved with classical aperture photometry on PSFhomogenized images using IRACLEAN (Hsieh et al. 2012), and another one (the FARMER) derived using a PSF-fitting tool (the Tractor; Lang et al. 2016) to extract the photometry.The covered area is ∼ 1.77 deg 2 and the total number of sources in the CLASSIC version is 1,720,700.See Weaver et al. (2022) for a detailed description of the two methods and catalogs.In this paper we use the CLASSIC version of the COSMOS2020 catalog.

Our sample
In this work we use the latest version from the A 3 COSMOS database (Adscheid et al. in prep.).This version combines the already tested process of the automatic mining of the ALMA archive with the new photometry presented in the COS-MOS2020 catalog (the new A 3 COSMOS catalog based on prior extraction from COSMOS2020 catalog is hereafter called A 3 C20 ).A 3 C20 is composed of 3215 individual pointings, coming from 171 different ALMA projects, covering ALMA bands from 3 to 9 In Figure 1, the wavelength distribution of the different pointings is reported.The ALMA bands are shown with different color-shaded areas.It is clear that the vast majority (∼ 80%) of the sample comes from ALMA bands 7 (orange area) and 6 (blue region), with more than 2500 observations available.In Figure 2  the survey, color-coded by the observed-frame wavelength.We also highlight three different pointing configurations, as representative of their complex spatial distribution in the survey: a) a case of partially overlapping pointings in the same band; b) concentric pointings in different bands; c) an extreme case of N>10 overlapping pointings in different bands.For further details see Section 4.
We select galaxies above 4.35σ (with σ being the local RMS at the position of each source, see Adscheid et al., in prep.), and it consists of 1620 sources with flux in at least one ALMA band.For our sample, 25% (441/1620) of the sources have a spectroscopic redshift (spec-z) and for 1069/1620 we have used the photo-z in the COSMOS2020 catalog.The spec-z are compiled from the literature (e.g., Riechers et al. 2013;Capak et al. 2015;Smolcic et al. 2015;Brisbin et al. 2017;Lee et al. 2017), the catalog by M. Salvato et al. (version 2017 September 1; available internally to the COSMOS collaboration) and the ALMA archive (Liu et al. 2019b).The photo-z used here are from Salvato et al. (2011); Davidzon et al. (2017); Delvecchio et al. (2017);Jin et al. (2018) and derived from either the CLASSIC or FARMER version of the COSMOS2020 catalog.Finally, 110/1620 sources do not have any redshift information.For this subsample, we computed their photo-z using CIGALE (Boquien et al. 2019) as described in the next section.

CIGALE SED fitting
In this work, we decided to perform the SED fitting of the A 3 COSMOS galaxies using CIGALE, a python SED fitting code based on the energy balance between the UV/optical emission by stars and the re-emission in the IR/mm by the dust.CIGALE is a highly flexible code, which allows one to choose among different individual templates for each emission components (e.g., stellar optical/UV emission, cold dust emission, AGN), across a broad parameter space.Furthermore, one of the most important features is the availability of active galactic nucleus (AGN) templates, which can be easily included in the fit, allowing a decomposition between star formation-powered and AGN-powered IR emission.In particular, derive the IR luminosity from the SED is crucial to compute the IR LF.In the following sections, we report the available photometry and the individual components used to perform the SED fitting, following recent SED-based studies (Ciesla et al. 2017;Lo Faro et al. 2017;Małek et al. 2018;Pearson et al. 2018;Buat et al. 2019;Donevski et al. 2020).Moreover, when needed, we included in the fit an input grid of redshifts spanning between z = 0 and z = 8 (with a step of ∆z = 0.1) to derive the best photo-z, if missing.

Photometric coverage
The A 3 C20 catalog, being the combination of COSMOS2020 catalog and the archival ALMA observations, takes advantage of a large photometric coverage, from the UV to the FIR/mm.To perform the SED-fitting, we have considered the following filters available in the COSMOS2020 catalog (Weaver et al. 2022): CFHT MegaCam u; Subaru Suprime Cam i, B, V, r, z; Subaru HSC y; VISTA VIRCAM Y, J, H, Ks; Spitzer IRAC channel 1, 2, 3 and 4; Spitzer MIPS 24 µm; Herschel PACS at 100 and 160 µm; Herschel SPIRE at 250, 350, 500 µm; JMCT SCUBA2 at 850 µm; ASTE AzTEC (1 mm) and IRAM MAMBO (1.2 mm).
To deal with all the ALMA frequency setup present in the A 3 COSMOS database, we have built artificial filters to be provided to CIGALE, corresponding to each observing wavelength in the A 3 COSMOS catalog.The filters are centered at that specific wavelength and are box-like, with their width being equal to 16 GHz.With this procedure, we added up to 330 continuous ALMA filters between 446 µm and 3325 µm to the CIGALE database.

CIGALE input templates
We modeled the stellar population with the bc03 stellar population synthesis model (Bruzual & Charlot 2003), to build the SED stellar component, and a delayed star formation history (SFH) with an optional second burst of star formation.We selected the dustatt_modified_CF00 (Charlot & Fall 2000) template to model the attenuation by dust and the dl2014 (Draine et al. 2014) to model the dust emission, both based on the assumption of having two different attenuation/emission sources represented by birth clouds and diffuse ISM.
Finally, we use the fritz2006 module (Fritz et al. 2006;Feltre et al. 2012) to model the AGN component in the SED.The AGN emission is described with a radiative transfer model which takes into account primary emission coming from the central engine (i.e. the accretion disc), scattered emission produced by dust and a thermal component of the dust emission.The individual input parameters are described in detail in Appendix A.

SED fitting results
We performed the SED fitting for the 1620 sources of our sample using CIGALE with the modules described above.For our purposes, we obtained the following output parameters: dust luminosity (i.e., the IR 8-1000 µm luminosity), stellar mass, AGN fraction ( f AGN ) contributing to the 5−40 µm total emission and τ (equatorial optical depth at 9.7µm) and Ψ (angle between equatorial axis and line of sight) parameters of the fritz2006 model and redshift for 110 sources.
Among the 1620 sources for which we performed the fit, we first remove those presenting a high reduced χ 2 (>10).Since our goals are strictly linked to the IR emission of these sources, we also computed the ratio between the ALMA observed flux and the best-fit flux at the same wavelength; sources with a ratio greater than the 5σ of the ratio distribution were removed and classified as bad SED, if we were not able to re-perform an acceptable fit.The 5σ threshold has been selected to be consistent with the "good SED" selection performed by Liu et al. (2019b) and Adscheid et al. (prep) in the catalog construction.We obtained a good fit for 1411/1620 galaxies, with 43/110 sources having no initial associated redshift.The remaining 67/110 galaxies without redshift had a bad fit (i.e., the IR-mm photometry is not matching the optical part of the SED), so we decided to exclude them from the final sample, but we consider them in the incompleteness estimate.
In Figure 3 we show the best-fitting SED of three representative objects: a Type II AGN, a Type I AGN and a galaxy without an AGN component.
In Figure 4 we report the redshift distribution (top left panel) of the sample, as well as the results from the SED fitting for stellar mass (bottom left), dust luminosity (bottom right) and AGN fraction (top right panel), for both the initial (1411) and final (i.e., after the cut, see Section 4, 189 galaxies).The redshift distribution peaks between z ∼ 1.5 and z ∼ 3.
The galaxies in the sample are massive with a peak in the distribution of stellar mass at log(M ⋆ /M ⊙ ) ∼ 11, consistent with them being DSFGs (Chapman et al. 2005;Simpson et al. 2014), although our sample contains sources with masses as low as 10 8 M ⊙ .The A 3 C20 galaxies are on average IR-bright, most of which with IR luminosities spanning from 10 11 L ⊙ to 10 13 L ⊙ .Using the best-fitting templates for each source, we are able to compute the AGN luminosity in a given wavelength range and, thus, the fractional AGN contribution to the total emission in that range.The 5-40 µm range is particularly sensitive to the presence of an AGN, therefore we considered this wavelength interval to derive the fraction of AGN, f AGN , contributing to the luminosity in this range.It can be seen that the distribution is bimodal, with most sources (∼ 65%) having an AGN fraction near 0 (i.e., they are likely not hosting an AGN) and the remainder having a f AGN higher than ∼ 0.2, peaking at ∼ 0.4 with a tail up to 0.8−0.9.The gap between f AGN = 0 and higher values is due to the grid we adopted for the SED fitting.Despite the ∼ 35% of sources is including an AGN component (the mean f AGN is ∼ 0.4) in the best-fit SED, only a small fraction (∼ 15%) is AGN-dominated (i.e., f AGN > 0.5).Also, since the AGN emission is mostly present in the 5 − 40 µm range, the contribution to the 8 − 1000 µm is only a small fraction.For consistency with Gruppioni et al. (2013) and other SFRD estimates we use the total (AGN+galaxy) IR emission, although we stress the AGN is strongly sub-dominant, and a more quantitative analysis is postponed to a future stand-alone work.For these reasons, we do not remove the AGN contribution to the total L IR .
Finally, we derived the dust obscured SFR for each galaxy using the Kennicutt (1998) relation, which assumes the SFR to be proportional to the IR luminosity: We show the SFR against stellar mass distribution (main sequence, MS) in Figure 5.Our sample is characterized by highly star forming galaxies, with SFRs up to ∼ 10 3 M ⊙ yr −1 , typical of starbursting galaxies at the considered redshift.In each panel, as a black solid line, we plot the main sequence from Speagle et al. ( 2014) computed at the mean redshift of the bin and an upper line (in black dashed) corresponding to 4 times the MS, which is indicative of the starbursting regime (Rodighiero et al. 2011).Starburst galaxies tend to cluster on a "sequence" above and parallel to the MS, as also noticed in other works (Caputi et al. 2017;Bisigello et al. 2018).As our SFR is derived directly from the IR luminosity, this bimodality can not be simply explained with the parametrization of the SFH or the presence of dust, but the exact origin of this effect needs further investigation with additional data.However, in our case, PIs selection can affect the MS distribution.In particular, we show in Figure 5 the sources distribution on the MS plane before (1411 sources) and after (189 sources) removing the targeted sources (i.e., the target of the individual pointing) in blue and red.As it can be noticed, the red points occupy mostly the region on the MS.Finally, for 6 out of 189 sources, we were unable to obtain a reasonable fit.The common characteristic of these sources is that they have an SED characterized by optical and ALMA photometry that does not appear to belong to the same object and can only be fitted by assuming extreme, unrealistic templates of dust emission.For this reason, and by inspecting the ALMA and optical cutouts, we have decided to treat these sources statistically in the derivation of the LF.The details of the procedure, as well as some example, are shown in section 5.1.

Turning the A 3 C20 sample into a blind survey
The A 3 COSMOS survey is based on observations with different sensitivities (i.e., limiting fluxes), resolutions and ALMA observed-frame wavelenghts.Resulting in a survey of observations with unique selection functions.Being a collection of several observations, almost each pointing is centered on a targeted source of interest for that observation.In this perspective, the use of the A 3 COSMOS survey for statistical purposes (e.g., LF, SFRD) needs a dedicated analysis for turning it into a blind survey.
In this section, we discuss the method used to turn a generic (F)IR-mm ensemble of pointed observations into a blind survey (see also Adscheid et al. in prep).

Blind surveys
To derive statistical properties of the sample through the corresponding areal coverage, a blind survey is needed.In order to do this, the next step should be followed and features fulfilled:  -In determining the limiting flux of A 3 COSMOS, we scale the root mean square (RMS) of each individual pointing to the corresponding RMS at our reference wavelength of 1.2 mm (Section 4.2); -We determine the total area spanned by our survey, accounting for primary beam attenuation and overlapping pointings (Section 4.3); -Finally, any possible bias due to the presence of a target has be taken into account (Section 4.4).

Root mean square conversion between different λ obs
The first step needed to achieve an unbiased survey from the A 3 COSMOS sample is the conversion of all observing wavelengths to a reference one.In our case, we choose the λ ref = 1200 µm wavelength in the observed frame, which falls in ALMA Band 6, being the most populated.This way, we can rescale all the fluxes at each observing wavelength (λ obs ) to a reference one (λ ref ), using the fluxes at the observed SEDs, in order to infer the observed ratio between λ ref and λ obs for each source.In particular we expect a decrease of the flux ratios when going to higher z and shorter wavelengths, approaching the rest-frame peak of dust emission in the SED.
As reported in Figure 6, we divided the sample into five redshift bins (z<1.5, 1.5-2.5, 2.5-3.5, 3.5-4.5,4.5-6).For visualization purposes, we plot only the most populated bands as grey points.The colored points indicate the median value for a certain λ obs for a specific redshift bin.It can be seen that, for shorter λ obs (close to λ ref ) the ratios are unsurprisingly low and almost redshift independent, while at longer λ obs , the ratios become larger and more redshift-dependent, especially in the highest z-bin.
Figure 7 shows the differences between the ratios in the different z-bins.We note that for z < 4.5, the correction curves are very similar to each other.For this reason we used only one mean curve (up to z = 4.5) and a different one for z > 4.5 in our analysis.

Radial limiting flux
We used the conversion curves obtained as explained in the previous section to convert the root mean square (RMS) in each pointing as if it was observed at 1200 µm.The pointings are now characterized by a uniform (in wavelength) sensitivity across their field of view (FoV).However, observing with the ALMA interferometer leads to a primary beam in which the sensitivity of the observation varies radially from the center to the outer regions.
For this reason, we divided the sky region covered by our pointings into 1" pixels and flagged those inside a pointing with a flux value corresponding to the sensitivity of that pointing (i.e., primary beam correction higher than 0.2).Then, we can  value of 1, decreasing to 0 towards the outer regions: with d being the distance of a pixel from the center of the pointing and σ being the FWHM/2.35.

[µJy / beam]
Fig. 8. Two overlapping pointings in band 6.The black circle represents the adopted primary beam size (out to a pbcor of 0.2), being color-coded by the inferred RMS.Lower RMS values are represented by darker colors and higher values by lighter ones.Outside our primary beam boundary, we have set values to the arbitrary high value of 10 6 µJy.In these type of cases, we adopted the area from the deepest pointing, in the overlapping region.
By dividing pixel by pixel the RMS by the pbcor corresponding to that pixel distance, we obtained a corrected RMS, which is then converted to the limiting flux by multiplying it by 4.35, which is the sensitivity cut of the A 3 COSMOS prior catalog.After this procedure we end up with a pixel map of limiting fluxes inside the pointings.However, for the pointings that are overlapping, the selection of an exclusive area is not straightforward.Indeed, the common area between two overlapping pointings has to be counted only once, in which a limiting flux from one of the two pointings needs to be assigned.
This issue can be dealt with by following the Avni & Bahcall (1980) method which coherently combines multiple samples at different depths.As shown in Figure 8, the limiting flux in the common area covered by two or more pointings is the one corresponding to the most sensitive observation at the reference frequency (i.e., the lowest limiting flux among the pointings).

Areal coverage
We proceed with deriving the cumulative areal coverage of our survey, by combining the effective area accessible by all pointings, at each limiting flux.This is done by counting the number of pixels with a certain limiting flux and deriving their cumulative distribution function (CDF).In this way, it is possible to associate an area to a limiting flux (S LIM ), which is the portion of area having a flux lower than S LIM .
In order to minimize potential biases due to pre-targeted ALMA sources contaminating the selection function, we disregard central targets as follows.We removed pointings in which a target source is not present in the central 1" (2060/3215).In this way, we took into account the possible positional offset of SMGgalaxy follow up with ALMA, which is of about 5" and, therefore, could be wrongly addressed as a serendipitous detection.By removing pointings without a target, we are sure to not bias our sample with possible off-center targeted galaxies.The 1" masking has been chosen through simulation varying the mask radius, increasing its size, and comparing the retrieved masked number counts with an input distribution (Adscheid et al., prep.).It has been shown that by increasing the central mask radius, no benefit is observed in the convergence of the retrieved number counts.In addition, we also removed 104/1620 sources that are potentially clustered (i.e.having a similar redshift to the target) following the criterion from Weaver et al. (2022), for which a source can be considered at the same redshift of another one if the difference in module between redshifts (z s and z t ) is smaller than 0.04(1+z t ).Finally, we removed the target of each pointing.
The final sample, obtained by considering both the target and redshift cuts, consists of 189 sources out of our initial sample of 1620 sources.We show in Figure 9 the areal coverage for the remaining 1155 pointings, i.e, after removing pointings without a central target detection.

A 3 COSMOS luminosity function
Deriving the areal coverage of our survey allows us to properly compute the luminosity function with the 1/V MAX method (Schmidt 1970).In the following sections we describe the method applied to derive the total infrared luminosity function and compare it with previous works.

The method
We derived the LF using the Schmidt (1970) method, that, based on the data, allows us to derive the LF without making any assumption relative to the LF shape.As already mentioned in section 4, the A 3 COSMOS survey is composed of ∼ 1155 individual pointings (after the cut), that can be considered as independent fields.For this reason, by following the Avni & Bahcall (1980) method, we are able to derive the effective areal coverage at each S LIM , across all pointings.The relation between area and limiting fluxes (see Figure 9) is then used to associate an accessible area above a certain flux for each source.
To compute the LF, we divided our sample in eight redshift bins, from z ∼ 0.5 to z ∼ 6 and into luminosity bins of 0.5 dex width, from log(L IR /L ⊙ ) = 10 to log(L IR /L ⊙ ) = 14.For each source, in a z-L IR bin, we measure the contribution to the LF in that bin by applying a redshift step of dz = 0.02 and Kcorrecting its SED from the lower to the upper boundary of the corresponding redshift bin, each time computing the observed flux at 1200 µm.This flux is used to infer the corresponding areal coverage at each dz by interpolating the previously derived areal coverage curve (Figure 9).We finally combine the effective area obtained in this way with the element of volume at each redshift step and obtain a comoving volume over which a given source is accessible: where V zmax and V zmin are the sum of the sub-volume in each dz shell up to the upper and lower limit of the bin, respectively.In particular V zmax can either be the volume at the upper bound of each dz bin, or the maximum volume reachable by considering the S /N limit of the survey (i.e., corresponding to the z at which the area would be 0).Finally, we correct the V MAX by taking into account completeness and spuriousness corrections, derived by Liu et al. (2019b), and we obtain the Φ(L, z) by summing each 1/V MAX in a certain luminosity-redshift bin.

Statistical treatment for potential HST-dark sources
As described in Section 3.3, we were unable to obtain a satisfactory fit for the SED of 6 out of 189 sources.Additionally, we identified these 6 sources as potential HST-dark objects.Consequently, due to having only ALMA flux information (the photo-z is associated with optical photometry), these objects were treated statistically in the LF analysis.Specifically, we assumed that these sources exist at a redshift greater than 3 and followed this procedure: firstly, we computed the cumulative redshift distribution for the rest of the sample; next, for each N-th (N = 100) random extraction during the LF calculation, we assigned a redshift to each of the sources by drawing from the cumulative distribution function (CDF) as a random sampler; once a redshift was drawn, we used the median SED of the sample to "fit" the ALMA flux and, consequently, obtain an infrared luminosity to be incorporated into the LF calculation.

IR LF
We obtained the IR LF by using the total (i.e., including all the SED components) L IR computed in the 8 − 1000 µm range.The z-bins (0.5-1.0; 1.0-1.5;1.5-2.0;2.0-2.5;2.5-3.0;3.0-3.5;3.5-4.5;4.5-6.0)are nearly equally populated, apart from the first bin which contains slightly less sources than the others.The LF is shown in Figure 11 as black points and red boxes representing the Poissonian uncertainties and the values are reported in Table 1.In each bin we overplot the completeness limit as a vertical black dashed line.This threshold is computed by re-scaling all the observed 1200µm fluxes of each SED to the same limiting flux, and then taking the L IR of SED which gives the highest L IR value at that limit.This latter value represents the value of L IR below which our sample is not 100% complete.The points below the completeness limit are reported as empty boxes.We also report existing estimates obtained from other studies.In particular, we compared our results with best fit from previous IR works, either from ALMA (ALPINE, Gruppioni et al. 2020) or Herschel (PEP+GOODS-H and PEP+HerMES, Magnelli et al. 2013;Gruppioni et al. 2013), when available for redshifts that are similar to ours.As it can be seen from Figure 11, we are mostly sampling the bright-end of the luminosity function, with our complete points being above log(L IR /L ⊙ ) ≃ 12 at each redshift bin.This is mainly a consequence of the fact that, at lower redshift, ALMA is sampling down the Rayleigh-Jeans tail and therefore even if deep pointing are available, they are inefficient compare to Herschel which sample the peak of the IR SED.The peak of dust emission can be better probed by ALMA (band 6 and 7 mainly) at high redshifts (z > 3).Despite this, our data points are in good agreement, within the errors, with the black triangles from Herschel (see Figure 12), whose z < 3 LF estimates are characterized by much higher statistics (smaller errorbars) and better sensitivity, being able to probe the IR LF down to the faint-end and the knee region (log(L IR /L ⊙ ) ≃ 11).From this consistency with the LF from Gruppioni et al. (2013) we can assert that, despite the poorer statistics in terms of the number of sources, the method used to turn the A 3 COSMOS into a blind survey is valide and accurate.At z > 3 − 3.5, where Herschel is probing further down the peak of the dust emission while ALMA start probing the peak, our estimates show a systematically higher normalization.Indeed, as stated by Gruppioni et al. (2013), in the 3 < z < 4.2 Herschel redshift bin, most of the sources have a photometric redshift and the PEP selection may be missing a fraction of high redshift galaxies, making this estimate a lower limit.

LF best-fit and evolution
In order to trace the number density of galaxies at different redshifts and infrared luminosities, we quantify the LF using a Schechter function.To this purpose, we modeled our LF with a modified Schechter function (Saunders et al. 1990), described by four free parameters: where α S and σ S are the slope of the faint end and the parameter shaping the bright end slope, respectively, and L * and Φ * represent the luminosity and normalization at the knee, respectively.The modified Schechter function is similar to a power law for L ≪ L * and behaves as a Gaussian for L ≫ L * .To find the L * and Φ * that best reproduce our LF, we performed Monte Carlo Markov Chain (MCMC) analysis using the PYTHON package emcee (Foreman-Mackey et al. 2013).It uses a set of walkers to explore the parameter space simultaneously.We carried out the MCMC analysis using 50 walkers, with 10000 steps (draws) and discarding the first 1000 sampled draws of each walker (burnin).The likelihood has been built in the form: We ran the MCMC using flat prior distributions for the two free parameters and with α S and σ S fixed to the values of Gruppioni et al. (2013) (i.e., α S =1.2 and σ S =0.5), with log(Φ * ) between −5 and −2 and log(L * ) between 10 and 13.The prior distribution is then combined with the likelihood function to obtain the posterior distribution.
Firstly, we fitted the ALMA points alone (Figure 11); in the lower redshift bins (i.e., 0.5 < z < 1.0 and 1.0 < z < 1.5), the individual fit shows a very low normalization with respect to the Herschel best-fit and large error bands, which do not allow us to constrain the fit parameters L * and Φ * in an accurate way.Between redshift 1.5 and 3.0 the best-fit LF is in a good agreement with Herschel, except for a slightly higher, but consistent, normalization in the 1.5 < z < 2.0 redshift bin.Finally, between z = 3 and z = 3.5, our best-fit is consistent with that of the ALPINE survey (Le Fèvre et al. 2020), but has lower Φ * at knee at z > 3.5.
Since the ALMA-only LF (red boxes) is not able to trace the lowest luminosities, we decided to take advantage of the great statistic and coverage in luminosity of the PEP+HerMES survey by Herschel (Gruppioni et al. 2013) up to z ∼ 3, with which our sample is consistent, while extending to higher redshift using our ALMA measurements.By means of this, it is possible to better characterize the faint-end and the knee of the LF at z > 3. Indeed, by using Herschel's LF data at all available redshifts (in Figure 12, only those within the redshift bins from 0.5 to 6 are shown), the MCMC analysis we conducted manages to better characterize the parameters of the Schechter function and subsequently trace its evolution at higher z.
Previous works have already pointed toward a redshift evolution of both typical density and luminosity of the IR galaxy population (see e.g., Caputi et al. 2007;Béthermin et al. 2011;Marsden et al. 2011;Gruppioni et al. 2013), characterized by an increase in luminosity and decrease in density with increasing redshift.In particular both Caputi et al. (2007); Béthermin et al. (2011) and Gruppioni et al. (2013) found a break in the redshift evolution, resulting in a steepening of the density evolution a z > 1 and a flattening for the luminosity at z > 2. For these reasons, we performed a second MCMC fit using the multiredshift information from each redshift bin, assuming an exponential shape and two different z break values (z ρ0 and z l0 ), for the evolution of Φ(z) * and L(z) * , expressed as: where Φ * 0 and L * 0 are the normalization and characteristic luminosity at z = 0, k ρ1 , k ρ2 , k l1 and k l2 are the exponents for values lower and greater than z ρ0 and z l0 for Φ and L, respectively.In this second fit, each point of the LF is associated to a redshift, corresponding to the median value of the galaxy population inside the individual redshift bin.By using these values, we can trace a more accurate evolution of Φ * and L * .From this fit we obtain two evolution curves (L * and Φ * vs z), which take into account the shape of the LF at each z.The priors used in this MCMC are given for the parameters that regulate the evolution of the LF, while α S and σ S are fixed to 1.2 and 0.5 (found by Gruppioni et al. (2013) at z ∼ 0.15); and anchor Φ * and L * to their value at z ∼ 0.7 found with Herschel-only measurements.
In particular, we fixed the L * and Φ * at z ∼ 0.7 to be the same found by Herschel.The best-fit curve at each z-bin is showed in Figure 12 as a dark-red solid line and shaded errors.
The trend with redshift of L * and Φ * is reported in Figure 13 and the results of the MCMC fit are reported in Tables 2 and  3, representing the 16th and 84th percentiles errorbands and the best-fit values at z = 0.The curve from the evolutionary MCMC and uncertainties are reported as dark-red line and shaded areas.The estimates of the individual fit are instead plotted for comparison as pink diamond for the ALMA-only case.From Figure 13 it is possible to notice that the values of L * and Φ * estimated by fitting only the A 3 COSMOS data points are slightly inconsistent with the results obtained instead using A 3 COSMOS plus Herschel (at z > 2).As mentioned before, this is due to the limited ability of ALMA to trace the dust peak emission at z < 2 and to the larger weight of Herschel data (containing more sources) in the combination.The evolution of Φ changes from z < z ρ0 = 0.89 +0.07 −0.14 to lower values, having an evolutionary trend of Φ * ∝ (1 + z) −0.55 for z < z ρ0 and Φ * ∝ (1 + z) −3.41 for z > z ρ0 .Overall, a decreasing trend is observed.This can be interpreted as a decreasing density of the bulk of star-forming galaxies at a given redshift.L * is showing two different trends with redshift below and above the break.In particular, for z < z l0 = 3.03 +0.87  −0.73 , the L * evolves as (1 + z) 3.41 and as L * ∝ (1 + z) 0.59 for z > z l0 , thus becoming flatter.The evolutionary trend of L * can be as-cribed to downsizing (Thomas et al. 2010), i.e., brighter (massive according to the S FR − M relation) galaxies formed earlier than their fainter counterparts.A similar analysis has already been performed by Gruppioni et al. (2013) using the Herschel PEP/HerMES LF (black empty circles in Figure 13), finding a consistently decreasing Φ * and increasing L * .However, their value of z l0 is lower than what we obtain (z l0 ∼ 3), being z l0 ∼ 2. The evolution is in a very good agreement with the Herschel only results (black empty points) and with the ALPINE estimates at z > 4.

Dust-obscured star formation rate density
Finally, using the best-fit luminosity function in each redshift bin, we derive the comoving SFRD.In particular, we first obtain the IR luminosity density (ρ IR ) by integrating the IR LF in each z bin, from log(L IR ) = 8 (to be consistent with others IR-based SFRDs) to log(L IR ) = 14 : Then, we can convert the IR luminosity density into a SFRD by applying the Kennicutt (1998) relation to convert L IR to SFR, assuming a Chabrier (2003) IMF.As described in Section 5.4, in addition to the individual fit, we performed an MCMC fitting taking the information from all the redshift bins together to derive the redshift evolution of L * and Φ * .In this way, we obtained the LF at each redshift from z ∼ 0.5 to z ∼ 6 and we integrated the LF[L * (z),Φ * (z))] in the full redshift range (for clarity we report the results of the fit for the eight bins of our IR LF).The result is an evolution curve of the SFRD with the redshift,     * The second column is the median value, while third and fourth columns show the lower and upper 16th boundaries.
shown in Figure 14.The values of the SFRD and errorbands, are reported in Table 4.The dark-red solid lines, with the shaded area, represent the lower and upper boundaries of the curve, derived as the errors on the integration of the LF at each redshift.
In order to compare our values with those from literature, we overplot previous results on the SFRD from Madau & Dickinson (2014); Béthermin et al. (2017); Lagache (2018) and from the IllustrisTNG simulation (Pillepich et al. 2018).Estimates of the dust-obscured SFRD from other IR-mm works (Gruppioni et al. 2013;Magnelli et al. 2013;Gruppioni et al. 2015;Rowan-Robinson et al. 2016;Liu et al. 2018;Gruppioni et al. 2020) are reported as red empty points with different markers.Blue empty points represent the SFRD from UV works, corrected from dust attenuation, by Dahlen et al. (2007); Reddy & Steidel (2009); Cucciati et al. (2012); Bouwens et al. (2012Bouwens et al. ( , 2015)).Thanks to the usage of Herschel data points, combined with the ALMA ones, in the LF fit, we were able to derive accurate estimates of the dust-obscured SFRD also at z < 3. We found the SFRD computed at z ∼ 0.5 − 1.0 follow the rise described by the Madau & Dickinson (2014) black curve and is consistent with the IR and UV estimates from previous works.Although our points up to z ∼ 2 have a lower normalization with respect to those from Gruppioni et al. (2013) and Magnelli et al. (2013), this discrepancy may be explained by the different type of fit performed which results in a different extrapolation to the faintend.From z ∼ 2 to z ∼ 3.5 the SFRD is following a flat trend, being above the Madau & Dickinson (2014).This flattening of the SFRD is compatible with predictions from models that envisage the early formation of massive spheroids (see Calura & Matteucci 2003).As mentioned before, at z > 3 our estimates are derived using ALMA alone (last three redshift bins).The observed trend is a decrease of the dust-obscured SFRD up to z ∼ 6, even though the limited number of sources does not allow us to strictly constrain it at those redshifts.We found the z > 3 dust-obscured SFRD to be consistent with the dust-corrected UV estimates and with the Madau & Dickinson (2014) one and, in the lower boundary, with the dust-obscured SFRD derived by Zavala et al. (2021), although with a ∼ 0.5 dex difference between the mean values.The upper boundary is consistent with Gruppioni et al. (2020), pointing towards a possibly underestimated SFRD at z > 3.However, this underestimation seems not to be as extreme as suggested by Rowan-Robinson et al. (2016) and Gruppioni et al. (2020).

Summary and Conclusions
Owing to the COSMOS2020 + ALMA photometric coverage, we were able to perform SED fitting using the CIGALE code to constrain the SFR, L IR and M ⋆ of the full ALMA sample.Applying the Avni & Bahcall (1980) method we were able to deal with the overlapping regions between several pointings, to evaluate correctly the depth of those regions.Moreover, we derived an SED-based method to homogenize the pointings at different observing wavelengths converting the RMS noise into that of the reference wavelength.By removing pointings without a detected target and by correcting for clustering, we were able to turn A 3 COSMOS into an untargeted, "blind-like" survey.We thus derived the total infrared luminosity function in the 8 − 1000µm band using the 1/V MAX (Schmidt 1970) method and using the already available data from Herschel.We estimated the LF from z ∼ 0.5 to z ∼ 6 and performed a MCMC simulation to fit the LF data, including the possibility for the characteristic parameters (Φ * and L * ) to evolve with redshift.By integrating the LF we derived the infrared luminosity density and, thus, the dust-obscured star formation rate density up to z ∼ 6.
We summarize the conclusions of this work as follows: -We find the A 3 COSMOS sample to be characterized by galaxies that are massive (log(M ⋆ /M ⊙ ) = 10−12) and bright in the IR (8 − 1000µm) band, with log(L IR /L ⊙ ) = 11 − 13.5.
-We find our LF to be in a good agreement with existing literature (in particular at z > 1), though pushing the SFRD to z ∼ 6 with unprecedented statistics.and the orange shaded area is the dust-obscured SFRD by Zavala et al. (2021), using the MORA survey.SFRD from the UV is plotted as blue empty points with different markers (Dahlen et al. 2007;Reddy & Steidel 2009;Cucciati et al. 2012;Bouwens et al. 2012Bouwens et al. , 2015) ) .
-Our MCMC analysis suggests a joint redshift-decreasing number density and a redshift-increasing IR luminosity for ALMA selected star-forming galaxies.This result is consistent with these galaxies being less frequent and more luminous (i.e., massive) going towards higher redshift.
-Our estimates of the dust-obscured SFRD are consistent with those from other IR surveys and with the UV dust-corrected estimates.Also, we found a broader peak in the SFRD, with a smooth decrease at z > 4, which suggests a significant contribution from the obscured SFRD at high redshift.Furthermore, a contribution to the SFRD, particularly at high redshifts, could originate from HST-dark sources (Franco et al. 2018b;Wang et al. 2019;Talia et al. 2021;Behiri et al. 2023), which are not a priori included in a prior sample, with the exception of the 6 added sources.
In conclusion, our study of the physical (i.e., stellar mass, dust luminosity, star formation rate) and statistical (LF and SFRD) properties of the ALMA sample of galaxies in the COS-MOS field resulted in the detection and characterisation of a starforming and starbursting dominated population, ∼ 40% of which is likely hosting an AGN.The A 3 C20 sample is also found to evolve both in luminosity (L * ) and density (Φ * ), and to significantly contribute to the total (IR+UV) SFRD also at z > 3.In order to improve our knowledge and to put tighter constraints to the evolution and formation of galaxies at higher z, more statistics is however needed.Being constantly updated, the A 3 COSMOS catalog will represent a key survey to reach newer results with a better statistics.The future COSMOS-Web survey (Casey et al. 2022) will also play an important role on covering the COSMOS field with JWST, allowing to obtain a better photometric coverage and, thus, allowing the physical properties of galaxies to be better constrained.

Fig. 1 .
Fig. 1.Number of pointings as a function of observing wavelength in the A 3 COSMOS database.The wavelength ranges of the ALMA bands are plotted as color-shaded regions.The most populated bands are 6 and 7, with ∼ 2000 pointings.

Fig. 2 .
Fig. 2. Archival ALMA observations in the COSMOS field used in this study.Individual pointings are plotted with different colors representing the observed ALMA bands.In the background, the COSMOS field is shown in grey.We show three zoom-in regions representative of possible classes of pointing configurations.The top-left panel (a) shows three overlapping pointings in the same band; the top-right panel (b) presents three concentric pointings in different bands; in the bottom-right panel (c) overlapping and concentric pointings, in different bands, are reported.

Fig. 3 .Fig. 4 .
Fig. 3. Examples of SED-fitting results for three different classes of objects.From left to right: unobscured AGN SED, SED without an AGN contribution and obscured AGN component.The blu filled circles and the red triangles represent data points and upper limits, respectively.The best-fit model is plotted as a black solid line.The stellar attenuated and unattenuated, the dust emission and the AGN emission are reported as yellow solid line, blue dashed line, red and green solid lines, respectively.

Fig. 5 .Fig. 6 .
Fig. 5. Main sequence of galaxies computed in six different redshift bins: 0.5-1.2,1.2-1.5, 1.5-2.0,2.0-3.0,3.0-4.0,4.0-6.0.Data of the full sample are reported as blue circles, while the red diamonds represent the final sample (4.4) used in our analysis.The black solid line and shaded area indicate the main sequence(Speagle et al. 2014) (computed at the mean values of the redshift bins) and 1σ dispersion.Finally, the black dashed lines represent the 4x MS, indicative of the starburst regime.

Fig. 7 .
Fig. 7. Ratio between the flux at 1200µm and in other bands against observed wavelength.Each point is the median value computed as described in Section 4.2.Red, brown, grey and yellow points are at z < 4.5, while the green ones represent the last redshift bin at z > 4.5.

Fig. 9 .
Fig. 9.Total areal coverage of the 1155 pointings after the cut for the lack of target detection within 1 arcsec.

Fig. 10 .
Fig. 10.Example of a source with optical identification (photo-z = 0.74) likely distinct from the ALMA galaxy.Left panel (from top to bottom): ALMA, acs-I, UV-j, and IRAC1 images are displayed with ALMA contours overlaid.It can be inferred that the optical/NIR object near the ALMA position is not centered in the ALMA galaxy.On the right: SED of the same object, showing the impossibility of fitting the optical and ALMA mm flux simultaneously.
Fig.11.A 3 COSMOS luminosity function (black points and red boxes).The individual redshift bin MCMC best-fit is plotted as pink solid lines and shaded error bands.The redshift ranges are reported in the upper-right corner of each subplot, while the luminosity bins are centered at each 0.25 dex, with a width of 0.5 dex (overlapping bins).The black vertical dashed lines represent the completeness limit of the L IR .The orange dashed, purple dotted and green dashed lines are the best fit LF obtained byGruppioni et al. (2013);Magnelli et al. (2013);Gruppioni et al. (2020).
Fig.12.A 3 COSMOS (red boxes and black points) + Herschel (black triangles) luminosity function.The dark-red lines and shaded areas are the best-fit obtained by using all the LF point from different z-bins together.The redshift ranges are reported in the upper-right corner of each subplot, while the luminosity bins are centered at each 0.25 dex, with a width of 0.5 dex (overlapping bins).The black vertical dashed lines represent the completeness limit of the L IR .The orange dashed, purple dotted and green dashed lines are the best fit LF obtained byGruppioni et al. (2013);Magnelli et al. (2013);Gruppioni et al. (2020).

Fig. 14 .
Fig. 14.Comoving star formation rate density evolution with z.The red diamonds are obtained through the integration of the LF best-fit in each z bin, while the dark-red shaded area is obtained by integrating the evolutionary LF from z = 0.5 to 6.Comparison curves from literature are reported as grey dash-dotted (pessimistic and optimistic cases, Lagache et al. 2018), dashed (IllustrisTNG), dotted (Bethermin et al. 2017) and black solid (Madau & Dickinson 2014) lines.The red empty points with different symbols represent different IR-mm estimates of the SFRD from Gruppioni et al. (2013); Magnelli et al. (2013); Gruppioni et al. (2015); Rowan-Robinson et al. (2016); Liu et al. (2018); Gruppioni et al. (2020) and the orange shaded area is the dust-obscured SFRD byZavala et al. (2021), using the MORA survey.SFRD from the UV is plotted as blue empty points with different markers(Dahlen et al. 2007;Reddy & Steidel 2009;Cucciati et al. 2012;Bouwens et al. 2012Bouwens et al. , 2015) ) .

Table 1 .
IR luminosity function as inferred using the A 3 COSMOS database.
* In the first column the luminosity bins are reported, while columns 2-8 report the Φ values in each luminosity and redshift bins.Bold (or italic) values represent independent luminosity bins.Values in round brackets indicate luminosity bins which are below the completeness limit and numbers in square brackets are the number of sources in each L-z bin.

Table 2 .
Best fit parameters at the knee of the IR-LF.

Table 3 .
Values at z = 0 obtained from the MCMC evolutive fit.Traina et al.: The IR LF and SFRD at z ∼ 0.5 − 6 from the A 3 COSMOS survey Evolutionary trends of Φ * (top panel) and L * (bottom panel) with redshift.The solid dark-red lines and shaded areas show the behaviours of Φ * and L * at all the redshifts, while pink diamonds are the estimates obtained fitting each redshift bin individually in the ALMA only case.Finally, black empty circles and yellow points represent estimates for L * and Φ * from Gruppioni et al. (2013) and Gruppioni et al. (2020), reported for comparison.
Gruppioni et al. (2013)ned from the MCMC evolutive fit.α S and σ S are fixed to the values found byGruppioni et al. (2013).Article number, page 12 of 17 A.

Table 4 .
SFRD obtained by integrating the LF best-fit in our eight redshift bins.