The ALPINE-ALMA [C II] Survey: [C II] 158 μm Emission Line Luminosity Functions at z ∼ 4–6

The [C II] emission line fluxes are measured from the ALMA data taken by ALPINE. A detailed description of this survey is discussed in Le Fèvre et al. (2020). The reduction of ALPINE data and the association with the extensive multiwavelength ancillary data are presented in Béthermin et al. (2020) and Faisst et al. (2019), respectively. A summary of the salient points of this survey is as follows.

ALPINE is a targeted survey, obtaining the [C II] spectral line and continuum observations with ALMA band 7 (275–373 GHz) for 118 main-sequence galaxies, selected by their UV luminosity at 1500 Å (note: not by stellar mass). These 118 galaxies are chosen to have optical spectroscopic redshifts in two intervals, z ∼ 4.40–4.58 and z ∼ 5.14–5.85, and to have 1500 Å absolute magnitude M_1500Å ≤ −20.2 mag. These targets have stellar masses of 10^8.5–10¹¹M_⊙ with an average value of 10^9.7M_⊙ (Faisst et al. 2019). The redshift slice of 4.6–5.12 was excluded because the redshifted [C II] lines would fall in a low transmission window for ALMA. The ALPINE primary targets are relatively luminous, almost all having L_1500Å > 0.5L*, corresponding to M_1500Å < −20.2. With the updated photometric measurements, three [C II] nondetections now have slightly fainter M₁₅₀₀ magnitudes, as shown by the dashed lines in Figure 1. For comparison, the UV LF at z ≃ 4.5–5.5 has M*_1500Å = −20.96 (Ono et al. 2018). The sample median value of L₁₅₀₀ is 1.5L*.

Zoom In Zoom Out Reset image size

As discussed in Béthermin et al. (2020), the mean [C II] line flux rms is 0.14 Jy km s⁻¹ for this sample. The spectral maps have yielded 75 [C II] emission lines with the peak flux to noise ratio SNR_peak ≥ 3.5σ and 43 nondetections (<3.5σ). In contrast to the high fraction of [C II] detections, the continuum maps have found only 23 detections at ≥3σ. The high fraction of [C II] detection supports the correlation between the integrated [C II] luminosity and star formation rate (SFR). Schaerer et al. (2020) investigates this topic using the ALPINE data.

Another significant result is that a good fraction of the z ∼ 4–6 sources is resolved—both spatially and kinematically—and many of them have multiple components and are clearly interacting. Several ALPINE papers (Béthermin et al. 2020; Le Fèvre et al. 2020; Loiacono et al. 2020) discuss this property. Jones et al. (2020) has focused on one triple merger system, deimos_cosmos_818760, with three [C II] components. Therefore, measuring the total line flux correctly is important for our calculations.

Since a large fraction of our sources are marginally resolved, their total line fluxes are generally larger than their peak fluxes, which are good representations of the total line fluxes only for point sources. Béthermin et al. (2020) has designed three different methods to measure total line fluxes. One is a simple circular aperture photometry within a 15 radius (f_aper). The second method is to sum up all fluxes within a 2σ contour above the background in the zero-moment map (f_clipped). And the last is a two-dimensional elliptical Gaussian fitting method (f_fit). Overall, these three methods give good agreement. For example, the 2D-fit and 2σ-clipped fluxes are consistent, with a small offset of 3% and standard deviation of only 0.17. As expected, aperture fluxes suffer larger errors, especially for faint sources (see Figure 15 in Béthermin et al. 2020). The aperture fluxes tend to underestimate the total fluxes compared with that of f_fit. For example, the mean and median of f_aper/f_fit are 0.85 and 0.93 respectively. For this paper, we computed the LFs using f_fit and f_clipped. As shown in Figure 5, the differences are within the errors, see details in Section 4.

One important and complex question is how to deal with [C II] sources with multiple components. This paper adopts the total [C II] flux measurements for all of the sources provided by Béthermin et al. (2020), without deblending (see gold and red contours in Figure C.4-6 in that paper). The rationale is that the sample has only four systems with multiple components at comparable brightness (within a factor of a few). For four of these (vuds_cosmos_5100822662, deimos_cosmos_818760, deimos_cosmos_873321, and deimos_cosmos_434239), each has two components separated by ≤1'' and with velocity differences ≤200–300 km s⁻¹ (see Figure C.1-6 of Béthermin et al. 2020). Our paper considers these two components are physically connected and belong to one system. It is worth noting that deimos_cosmos_818760 is in fact a triple system, with a third, fainter component ≥2'' away. Our calculation does not include the fluxes from this third system (see Jones et al. 2020). Finally, we comment on vuds_cosmos_5101209780, whose primary target (component a) has a companion (b) at a separation of >2'' with a flux brighter by a factor of 2. Our paper only considers Component a. Loiacono et al. (2020) regards Component b as a part of the serendipitous [C II] detection sample.

Most [C II] systems have only faint tails or slightly extended minor components. Multicomponent decomposition will not produce a more accurate measurement of LF because, as we discuss in the sections below, our derived LF is computed in luminosity bin width of 0.5 dex. Small flux boosting due to faint companions will not change our results.

Our [C II] sample is selected by a combination of UV luminosity, redshift, and [C II] line flux cuts. The redshift cuts produce a volume limited sample. The third selection imposes ALMA [C II] flux limits, which vary from pointing to pointing. Such a unique selection makes the LF calculation in some ways simpler but can also be confusing.

We first describe the basic concept of our method. Figure 1 displays UV absolute magnitude M_UV versus [C II] luminosity log₁₀(L_[CII]) for the 75 [C II] detections. The bottom panel of Figure 1 presents the number of [C II] emitters, N_[CII], per luminosity bin Δlog₁₀(L_[CII]). In the simplest form, the [C II] luminosity function at the ith bin is

$$\begin{eqnarray}&&{\phi }_{i}=\displaystyle \frac{1}{{\rm{\Delta }}(\mathrm{log}{L}_{[\mathrm{CII}]})}\displaystyle \frac{{N}_{[\mathrm{CII}],i}}{{V}_{[\mathrm{CII}],i}}\times {C}_{{UV},i}\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&=\ \displaystyle \frac{1}{{\rm{\Delta }}(\mathrm{log}{L}_{[\mathrm{CII}]})}\displaystyle \frac{{N}_{[\mathrm{CII}],i}}{{V}_{[\mathrm{CII}],i}}\times \displaystyle \frac{{\rm{\Delta }}{M}_{{uv}}{\phi }_{{uv}}{V}_{{uv}}}{{N}_{{uv},i}},\end{eqnarray} \tag{ 2 }$$

where N_[CII],i is the number of [C II] emitters in the ith bin, V_[CII],i is the volume available to these sources, and C_UV,i is the incompleteness factor. This incompleteness comes from the fact that the ALPINE targets are a small subset of all galaxies that meet our UV luminosity cut (M_1500A < −20.2) among the parent galaxy population. As shown in Equation (2), C_UV,i is equal to the ratio of the number of expected UV galaxies from the known UV LF to the number of ALPINE UV galaxies within a volume element V_uv. Note that if the $\mathrm{log}{L}_{[\mathrm{CII}]}$ to M_UV relation were exactly 1:1, the above could be directly computed as written for each luminosity bin. However, Figure 1 shows the two variables have significant scatter. This means that in fact we need effectively a weighted C_UV that is based on the M_UV distribution among all ALPINE sources within the ith [CII] luminosity bin or in other words use ${C}_{{UV},i}=(1/{N}_{[\mathrm{CII}],i}){\sum }_{j=1}^{{N}_{[\mathrm{CII}],i}}{C}_{{UV},j}$ where j is the index across all M_UV bins sampled by the ALPINE sources.

The second issue is that the ALPINE survey is a targeted survey as opposed to a flux-limited survey, which means one would expect the volume elements V_[CII] and V_UV to exactly match and therefore cancel out in Equation (2) (where the two volumes appear, respectively, in the denominator and nominator). Recall that a comoving volume is computed as A × [D_M³(z_lower) − D_M³(z_upper)], where A is the areal coverage and D_M is the comoving distance, a function of redshift. Now, whatever A we choose to adopt in our calculation of V_[CII] has to be the same as the one used to calculate the expected number of UV targets in our survey. Thus, the areas in the denominator and nominator cancel. For the ALPINE survey, the [C II] LF calculation does not rely on the areal coverage value and we only need to know the values of z_lower and z_upper. For the UV galaxies, by selection, the volume is defined by its redshift slices with z_lower and z_upper set as 4.403–4.585 and 5.135–5.85, respectively, for the two slices. However, for the [C II] detections, the volume is not the same as defined by the redshift slices because each ALMA pointing has different sensitivity limits. For all [C II] sources in the z ∼ 4.5 bin, their ALMA maps are deep enough, so that their maximum allowed redshifts (thus volumes) are the same as the upper limits of the redshift slice (z_upper). At the z ∼ 5.5 bin, 8 of 29 sources have slightly shallower ALMA maps; therefore, their maximum allowed (or maximum detectable) redshifts are slightly less than z_upper = 5.85. Therefore, the ALPINE survey does have a flux-limited element. So for the vast majority of sources, the volume elements in Equation (2) exactly cancel. But for eight of our sources this is not the case. Putting these arguments together, we have a final version of our LF calculation per luminosity bin i, given by Equation (3).

$$\begin{eqnarray}&&{\phi }_{i}=\displaystyle \frac{1}{{\rm{\Delta }}(\mathrm{log}{L}_{[\mathrm{CII}]})}\displaystyle \sum _{j=1}^{{N}_{i}}\displaystyle \frac{{C}_{{UV},j}}{{V}_{[\mathrm{CII}],j}},\end{eqnarray} \tag{ 3 }$$

where ϕ is for the ith luminosity bin, and summation goes through all [C II] sources within this bin with j = 1, 2, ... N_i. The reader will note that this equation is just a modified 1/V_max method (Schmidt 1968). However, this is not the true 1/V_max method because the calculation is pegged directly to the published UV LF via the C_UV incompleteness correction factor. Therefore any issues with this method that involve small number statistics, clustering effects, as well as knowledge of precise areal coverage are all canceled out as long as one trusts the accuracy of the state-of-the-art UV luminosity function calculation we use (Ono et al. 2018). For example, if our targets happen to be clustered in our redshift ranges, leading to number density of [CII] emitters enhanced by some factor, the ratio of expected UV sources to known UV sources within the ALPINE target would be slightly lower by the same factor fully canceling it out.

Figure 2 shows our calculated C_UV for both the two redshift bins. For the sake of this calculation and display, we adopted A = 2 square degrees (size of the COSMOS field) but this assumption cancels out in the LF calculation as discussed above. Here C_UV was calculated in bins of UV magnitudes of width 0.5, but these values were interpolated so that we can have C_UV,j corrections corresponding to each individual source's UV magnitudes in Equation (3). In Figure 2, the error on C_UV is derived from the standard deviation of the expected number of galaxies per UV magnitude bin calculated from 1000 Poisson draws per bin. This Poisson uncertainty dominates the uncertainty of the UV LF (Ono et al. 2018).

Zoom In Zoom Out Reset image size

3.1.1. Error Estimates for the [C II] LF

In Equation (1), we see that our estimates of the LF per [C II] luminosity bin are based on the product of N_[CII],i and the UV-incompleteness correction factor C_UV,i. As such the uncertainty in ϕ_i is given by the usual rule of error propagation for products:

$$\begin{eqnarray}&&\sigma ({\phi }_{i})={\phi }_{i}\sqrt{{\left(\displaystyle \frac{{\delta }_{{N}_{[\mathrm{CII}],i}}}{{N}_{[\mathrm{CII}],i}}\right)}^{2}+{\left(\displaystyle \frac{{\delta }_{{C}_{{UV},i}}}{{C}_{{UV},i}}\right)}^{2}}.\end{eqnarray} \tag{ 4 }$$

The first of the error terms, the fractional error on the number of sources per [CII] bin is derived assuming Poisson statistics and based on the lower and upper 1σ uncertainties tabulated in Gehrels (1986). Due to the relatively small number of sources in all of our bins, this turns out to be the dominant error term. The second term uses the C_UV uncertainties as described above and shown in Figure 2. Recall that ${C}_{{UV},i}=(1/{N}_{[\mathrm{CII}],i}){\sum }_{j=1}^{{N}_{[\mathrm{CII}],i}}{C}_{{UV},j}$, which translates to $\displaystyle \frac{{\delta }_{{C}_{{UV},i}}}{{C}_{{UV},i}}$ = $\tfrac{\sqrt{\sum {\delta }_{{C}_{{UV},j}}^{2}}}{\sum {C}_{{UV},j}}$.

In our method above, we already accounted for the incompleteness of the ALPINE sample relative to its parent sample of M_UV < −20.2 sources. However, there is a second type of incompleteness due to [C II] nondetections. Of the 118 ALPINE targets, 43 sources only have [C II] flux limits, with 22 and 21 objects in the z ∼ 4.5 and the z ∼ 5.5 (optical spectroscopic redshifts) bins, respectively. To account for their maximum potential contribution to the [C II] LFs, we take the nondetections' noise estimates, σ, and assign them a maximum luminosity of 3σ. This is a reasonable upper limit for the nondetections, since the distribution of SNR_peak for nondetections peaks at ∼2 (Béthermin et al. 2020) and the threshold for detection is SNR_peak ∼ 3.5.

Figure 3 shows the contribution of nondetections to our two lowest luminosity bins. We consider three different measurements of the noise, aggressive (lowest noise), normal, and secure (highest noise) σ, which are described in Section 3.3 in Béthermin et al. (2020). Figure 3 shows that the distribution of 3σ upper limits measured by the secure method is systematically shifted to higher luminosity roughly by 0.2 dex compared to that by the aggressive method. We consider the maximum potential contribution of nondetections to each luminosity bin to come from whichever noise estimate maximizes the number of nondetections in the particular bin (i.e., aggressive method at the lower luminosity bin, and the secure method in the higher luminosity bin). Because the nondetections could in principle lie below the 3σ limit, these are the upper limits. We remind the reader that galaxies that are treated with such upper limits lie exclusively in the lowest two luminosity bins so in all other bins, as we discuss below, our estimates are in fact strict lower limits. In forthcoming figures we show both the LFs computed from the detected sources alone, and those that include the maximum potential contribution from the nondetections.

Zoom In Zoom Out Reset image size

The last type of incompleteness is that ALPINE is not a blind [C II] emission line survey. Rather it is a survey specifically targeting galaxies selected by 1500 Å UV luminosity. The high fraction of [C II] detections in the ALPINE sample implies that the UV selection should detect [C II] emitters efficiently; however, the UV luminosity and [C II] line luminosity do not have a tight linear correlation (see Figure 1). The ALPINE selection would certainly miss luminous [C II] emitters with low UV luminosity, i.e., M_1500Å > −20.2 magnitude. Thus, the calculated LF based on the UV-selected sample sets the lower bounds to the [C II] LF. Including our nondetections by assuming they are at 3σ pushes our estimates up, in particular, in our two lowest luminosity bins. However, these are still contributions from sources above our UV magnitude limit. To fully account for [C II]-emitters that are potentially UV-faint—i.e., below the ALPINE target selection threshold, we need a blind [C II] survey. The ALPINE surveys allow for the next best thing—i.e., search for serendipitous sources within the same ALMA data as our targets. This search was conducted by Loiacono et al. (2020). In Section 4.2 we compare our results with those of Loiacono et al. (2020) who compute the [C II] LFs using serendipitous [C II] sources in the ALPINE survey. Because of potentially uncertain redshift identifications and strong clustering effects, the results from the serendipitous sources should be considered upper bounds. The comparison of the two estimates of the z ∼ 4–6 [C II] luminosity function therefore helps to overcome some of the systematic uncertainties associated with each method and affirm the overall results on the number density and luminosity function of [C II] emitters at these redshifts as seen from the ALPINE survey.

Figure 4 presents our [C II] LFs at z ∼ 4.5 and z ∼ 5.5, with the number of sources in each bin labeled. The downward arrows mark the upper limits if we include the [C II] nondetections at the two lowest luminosity bins. These are connected by the dashed vertical lines to the points derived using only significant [C II] sources. Table 1 gives our derived LF values at z ∼ 4.5 and z ∼ 5.5 using two different binnings.

Zoom In Zoom Out Reset image size

Table 1.  ALPINE [C II] LFs at z ∼ 4–6 in 10⁻⁴ Mpc⁻³ dex⁻¹

Log10(L[CII])	${\phi }_{[\mathrm{CII}]}^{a}$	${\phi }_{[\mathrm{CII}]}^{{up},b}$	ϕ[CII]	${\phi }_{[\mathrm{CII}]}^{{up},c}$
L⊙	z ∼ 4.5		z ∼ 5.5
7.75	⋯	⋯	${1.58}_{-1.02}^{+2.09}$	${11.1}_{-1.02}^{+2.09}$
8.25	${7.98}_{-2.77}^{+3.94}$	${26.4}_{-2.77}^{+3.94}$	${11.1}_{-3.15}^{+4.20}$	${24.0}_{-3.15}^{+4.20}$
8.75	4.12 ± 0.76	⋯	${2.47}_{-0.70}^{+0.94}$	⋯
9.25	${0.55}_{-0.22}^{+0.33}$	⋯	${0.074}_{-0.04}^{+0.072}$	⋯
9.75	${2.0}_{-1.29}^{+2.63}$	⋯	⋯	⋯
7.95	${0.59}_{-0.38}^{+0.78}$	⋯	${3.59}_{-1.32}^{+1.93}$	⋯
8.45	${9.72}_{-2.30}^{+2.30}$	⋯	${10.8}_{-2.97}^{+3.92}$	⋯
8.95	${2.31}_{-0.49}^{+0.49}$	⋯	${0.76}_{-0.26}^{+0.37}$	⋯
9.45	${2.03}_{-0.97}^{+1.61}$	⋯	⋯	⋯

Note. a: ϕ_UV is the [C II] LF derived based on the ALPINE primary targets, i.e. the UV sample. b and c: ϕ_UV^up refers to the corrected [C II] LFs by taking into account of the [C II] nondetections in the ALPINE sample. These mark the upper limits to the two lower luminosity bins because we assume that these nondetections have 3σ line luminosity.

Download table as:  ASCIITypeset image

Figure 4 shows the following conclusions. First, the lower limits on the [C II] LFs at z ∼ 4–6 derived from the ALPINE UV sample are consistent with the z ∼ 0 LFs at the luminosity range of log₁₀(L_[CII]/L_⊙) ∼ 8.25–9.25 (bin center, 0.5 dex bin width). At L_[CII] ≥ 10^9.5L_⊙, the volume density of high-z [C II] emitters could be higher than that of z ∼ 0, but our measurement is only at a (1–2)σ significance level (68%–95%). Such a large uncertainty requires future confirmation.

Second, we do not measure significant differences between [C II] LFs at z ∼ 4.5 and z ∼ 5.5, within the uncertainties. Considering that LF is essentially the volume density of [C II] emitters, we compute the LF at z ∼ 4–6 simply by using the geometric mean of the values at z ∼ 4.5 and z ∼ 5.5, except at the lowest and highest luminosity bins where only z ∼ 4.5 LF measurements are available. For these two bins, we take only the z ∼ 4.5 LF values, assuming the volume density does not change between z ∼ 4.5 and 5.5. For simplicity, the subsequent figures show the [C II] LF in a single redshift bin of 4 < z < 6.

To check the dependence of our results on the specific binning, in Figure 4 we compare our LF with that derived using a different binning (open symbols). The results are consistent with each other.

To check the dependence of our results on the specific photometry method, Figure 5 compares the LFs using f_fit (which we adopt throughout this paper) and f_clipped. Within the errors, the difference is not significant.

Zoom In Zoom Out Reset image size

To check the dependence of our results on the method for the [C II] luminosity function, we also examined an alternative approach. This method uses the M_UV versus L_[CII] distribution in Figure 1 to construct the probability distribution in each UV magnitude bin for a given L_[CII]. The expected number density of sources from the UV LF is spread among the L_[CII] bins according to these distributions. Summing the contributions from each UV bin gives us the total number density expected in each L_[CII] bin. The results are fully consistent with our method described above.

Finally, we carefully check to see if source blending can be an issue, boosting our highest luminosity bins. We check the [C II] maps (Béthermin et al. 2020) for the four most luminous sources included in the highest-L bins in Figure 4. Three of the four sources are [C II] bright and isolated objects, and one has an optical counterpart with very low UV luminosity. The fourth source, deimos_comos_818760, is a triple system and has been published in Jones et al. (2020). The whole system has components C, E, and W with log₁₀(L_[CII]/L_⊙) = 9.48, 9.21, and 8.7 respectively (Jones et al. 2020). With small spatial and velocity separations, components C and E are considered as a major merger (Jones et al. 2020) and is counted as a single source with log₁₀(L_[CII]/L_⊙) = 9.66 in our calculation. Recall that our LF bin width is 0.5 dex. Component W is included in the Loiacono et al. serendipitous sample. The UV counterpart to component C is the ALPINE primary target; however, the counterpart to component E has high dust extinction, and is not detected in the HST F814W band (Koekemoer et al. 2007, 2011), but is detected in the deep K band, with low UV luminosity, M₁₅₀₀ well below the ALPINE cutoff. Treating these two components separately will lead to two (as opposed to one) [C II]-bright but very UV faint sources with corresponding large C_UV corrections and would boost our LF at the highest luminosity. We choose to treat them as one source as a more conservative approach.

While our calculations are robust (as the above tests show), they are in fact only lower bounds on the true [C II] luminosity function at these redshifts. This is because any UV selected sample such as ALPINE is biased against [C II] emitters with low UV luminosity. As Figure 1 shows, the highest luminosity bin of the [C II] LF indeed has one [C II] bright and UV-faint galaxy. Its large C_UV noticeably boosts the LF value in this luminosity bin. This source (COSMOS_DEIMOS_873756) is a so-called "4.5 μm excess" source, which is off the main sequence with a large stellar mass limit (<10^10.2M_⊙) and very low SFR (${\mathrm{Log}}_{10}(\mathrm{SFR})={0.72}_{-0.2}^{+0.5}{{\rm{M}}}_{\odot }$ yr⁻¹) (Faisst et al. 2020). Note, however, that the stellar mass is an upper limit and the SFR has large error bars. If this source is atypical of the general UV population of the same UV magnitude, then our calculated UV incompleteness correction is not applicable here. Therefore, our estimate of the luminosity function in this highest luminosity bin has this additional systematic uncertainty that needs to be kept in mind.

The LFs presented above are lower limits on the number densities of [C II] emitters at z ∼ 4–6 because they miss potential sources that are bright in [C II] but faint in UV. Fortunately, with the 118 ALMA maps along separate pointings, searching for serendipitous [C II] emitters at 4.40 < z < 4.58 and 5.13 < z < 5.85 can quantify this effect. In a companion paper, Loiacono et al. (2020) carried out a blind search for all line emitters within the 118 data cubes. They identified eight [C II] emitters at z ∼ 5 confirmed by the optical spectroscopic redshifts, and four possible [C II] candidates with only photo-z. Of these 12 sources, 11 are associated with their respective primary targets in both spatial and velocity space, and 1 is truly a field object, for details, see Loiacono et al. (2020). These 11 "clustered" sources are used to compute the [C II] LF shown as green points in Figure 6. The one field serendipitous object is used to constrain the [C II] LF shown as a diamond point in salmon color in Figure 6. Figure 6 compares our derived LF at z ∼ 4–6 with that of Loiacono et al. The LF based on the serendipitous sources is strongly affected by cosmic clustering and by uncertainties in the identification of some line emitters. Thus it sets an upper boundary to the bright end of the [C II] luminosity function at z ∼ 4–6, whereas our results are the lower limits. The true [C II] LF lies in the shaded region, covering the purple squares and green dots and extending out to the upper and lower ends of the error bars in Figure 6.

Zoom In Zoom Out Reset image size

Considering the upper and lower limits derived from the serendipitous and UV-selected ALPINE samples, we conclude that the space densities of [C II] emitters at z ∼ 4–6 are at least as high as that of z ∼ 0. At log₁₀(L_[CII]/L_⊙) > 9.5, there could be more high-z luminous [C II] emitters than locally, but our estimate has large errors. This uncertain but tantalizing result should be validated by future studies. Together, the UV-selected and serendipitous ALPINE samples both point toward significant number densities of luminous [C II] emitters at z ∼ 4–6. At the low luminosity end with log₁₀(L_[CII]/L_⊙) ≤ 8.75 both LFs are consistent with the z ∼ 0 LF, suggesting that the space density of [C II] emitters at z ∼ 5 in this luminosity range is similar to that of z ∼ 0. However, the number densities at log₁₀(L_[CII]/L_⊙) > 9 are far in excess of the z ∼ 0 [C II] luminosity function of Hemmati et al. (2017) suggesting potentially considerable evolution.

In this section, we further test our conclusion above by comparing the direct ALPINE [C II] LFs with indirect estimates of the [C II] luminosity function based far-IR, and CO luminosity functions at z ∼ 4–6 and adopting typical conversion relations. Figure 7 shows this comparison. Below we describe in more detail the assumptions behind this figure.

Zoom In Zoom Out Reset image size

Our first comparison sample is again drawn from serendipitous continuum sources in the ALPINE survey. Here we use the results in Gruppioni et al. (2020) who derive an IR luminosity function based on serendipitously detected continuum sources in the ALPINE data. We also use the 250 μm LF from the SCUBA-2 survey (Koprowski et al. 2017). To convert this to a total IR luminosity function, we adopted L_IR/L_250μm = 38.5 calculated using the infrared spectral energy distribution (SED) template, constructed from data including ALPINE continuum measurements by Béthermin et al. (2020). Finally, to convert both of these into [C II] luminosity functions we adopt the L_IR/L_FIR = L_8–1000μm/L_42–122μm = 1.3 and ${\mathrm{log}}_{10}({{\rm{L}}}_{\mathrm{FIR}}/{{\rm{L}}}_{[{\rm{C}}\mathrm{II}]})=2.69$ recently compiled for high-z galaxies (Zanella et al. 2018).

Next we include the [C II] LFs at z ∼ 5.8, converted from the CO(1–0) LF published by the ALMA SPECtroscopic Survey (ASPECS, pink shaded region) (Decarli et al. 2019) and from the CO(2–1) LF from the CO Luminosity Density at High Redshift survey (COLDz, light shaded green) (Riechers et al. 2019). Here we use CO(2–1)/CO(1–0) ratio of 1 (Riechers et al. 2019). For comparison, the average value for this ratio is 0.85 for submillimeter galaxies. The ASPECS and COLDz CO LFs are published in ${{\rm{L}}}_{\mathrm{CO}(1-0)}^{{\prime} }[{\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{pc}}^{2}]$. We transformed the CO(1–0) LF to the estimated [C II] LF by using

$$\begin{eqnarray}&&\begin{array}{l}{\mathrm{log}}_{10}({L}_{\mathrm{CO}(1-0)}/{L}_{\odot })={\mathrm{log}}_{10}({{\rm{L}}}_{\mathrm{CO}(1-0)}^{{\prime} }\\ \quad \ \times \ [{\rm{K}}\,\mathrm{Km}\,{{\rm{s}}}^{-1}\,{\mathrm{pc}}^{2}])-4.31\end{array}\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{\mathrm{log}}_{10}[{L}_{[\mathrm{CII}]}/{L}_{\mathrm{CO}(1-0)}]=3.6.\end{eqnarray} \tag{ 6 }$$

Combining Equations (5) and (6), we have the ratio of ${\mathrm{log}}_{10}[{L}_{\mathrm{CO}(1-0)}^{{\prime} }/{L}_{[\mathrm{CII}]}]=0.9$, which is used to convert the CO(1–0) LF in ${{\rm{L}}}_{\mathrm{CO}(1-0)}^{{\prime} }$. The first relation is based on the equations defined in Section 2.4 in Carilli & Walter (2013). The ${\mathrm{log}}_{10}[{L}_{[{\rm{C}}\mathrm{II}]}/{L}_{\mathrm{CO}(1-0)}]$ ratio can span a range of 3–3.8 as measured from a sample of normal galaxies at z ∼ 0 and z ∼ 1–2 (Stacey et al. 1991, 2010). Zanella et al. (2018) recently examined a larger sample of galaxies with published L_[CII] and L_CO measurements at z ∼ 0–6. We adopt the value of 3.6 in Equation (6), consistent with the calculation of molecular mass described in Section 4.5. Our adopted value makes a conservative conversion to [C II] luminosity. Extremely metal-poor, low-mass, blue compact dwarf galaxies can have one order of magnitude higher [C II]-to-CO(1–0) ratio (Cormier et al. 2014). But these are unlikely to be good analogs, since the averaged star formation efficiency for the ALPINE sample is an order of magnitude higher than the local sample (Madden et al. 2013; Cormier et al. 2014; Faisst et al. 2019). In addition, our adopted value ignores that at z ∼ 4–6, CMB photons will increase line excitation as well as background for the CO(1–0) line (da Cunha et al. 2013). The CO(1–0) line would be suppressed by ∼20% at z ∼ 6 (Vallini et al. 2018). This would potentially push these estimated [C II] LFs to even higher luminosities, but is still well within the overall uncertainties in this conversion relation.

Figure 7 shows that within the large uncertainties, there is good agreement between all these constraints on the [C II] luminosity function at z ∼ 4–6. Combined, they strengthen even further our conclusion that there has been strong evolution in the number densities of L_[CII] > 10⁹L_⊙ sources between z ∼ 5 and z ∼ 0.

Finally, including the data from our work, Loiacono et al. (2020), the converted values using the FIR LFs (from ALPINE continuum data) from Gruppioni et al. (2020) and Koprowski et al. (2017), we perform a fit to the data assuming a Schechter function. We derive ${{log}}_{10}({L}_{[\mathrm{CII}]}^{* }/{L}_{\odot })=9.5\pm 0.6$, ϕ* = (8.4 ± 8) × 10⁻⁴ Mpc⁻³ dex⁻¹ and α = −1.1 ± 0.3.²² The Schechter fit heavily relies on the FIR LFs and have large errors because both our and Loiacono et al.'s results are only the lower and upper limits. It is expected that this Schechter form does not fit to either our data or those of Loiacono et al. Performing the functional form fit here is for convenience when making comparisons and when computing [C II] line luminosity density. At face value, the z ∼ 4–6 ${L}_{[\mathrm{CII}]}^{* }$ has increased by a factor of 20 compared to L_[CII]* ∼ 2.2 × 10⁸L_⊙ at z ∼ 0 (Hemmati et al. 2017). We caution that systematic errors could be large.

Figure 8 shows the model [C II] LFs from Popping et al. (2016) and Lagache et al. (2018) overlaid on the same data as shown in Figure 7. These models combine semianalytic galaxy evolution models with radiative transfer calculations to make predictions on [C II] line emissions and LFs. As already noted in our z ∼ 0 [C II] LF paper (Hemmati et al. 2017), the Popping et al. (2016) model is definitely ruled out by the observations. Lagache et al.'s (2018) model is much improved, but still fails to match the abundant luminous [C II] emitters at z ∼ 4–6. The Lagache et al. (2018) model is reasonably consistent with the data at log₁₀(L_[CII]/L_⊙) > 10, but the data overall are more consistent with a Schechter functional form than the power-law form of this model.

Zoom In Zoom Out Reset image size

Several cosmological simulations have investigated CO(1–0) and [C II] emission at high-z, most of them focusing on the L([CII])-SFR relation and associated deviations. As illustrated in Ferrara et al. (2019), the physical processes of [C II] emission are complex, and the assumed gas density in PDRs (or size and filling factor), ionization factor, metallicity, and starburstness can all affect [C II] luminosity. Some models can reproduce both the luminous and weak [C II] sources in the L([CII])-SFR plane (Vallini et al. 2015; Lagache et al. 2018; Ferrara et al. 2019). The discrepancy between our LF and the models, however, implies that the current models for synthesizing the high-z [C II] populations likely involves other factors that are still not well understood. More work is clearly needed.

It has been widely recognized that the evolution of the molecular gas mass as a function of redshift holds clues to understanding galaxy formation and evolution more broadly (for a review see Carilli & Walter 2013). In particular, we want to compare with the model for ρ_molgas as a function of z from Liu et al. (2019). Liu et al. (2019) combined their ALMA galaxy catalog with 20 other published data sets and made a comprehensive analysis of scaling relations between stellar mass, gas fraction, and redshift. They expressed the molecular gas mass density, ${\rho }_{{mol}}(z)={\sum }_{{M}_{* }}{{\rm{\Phi }}}_{{SMF}}(z,{M}_{* })$ × M_* × μ_gas(z, M_*, ΔMS), as a function of stellar mass function Φ and gas fraction μ. They integrate down to M_* = 10⁹M_⊙. Using the Liu et al. (2019) μ_molgas = M_molgas/M_* ∼ 2 at z ∼ 5 and ${\alpha }_{[\mathrm{CII}]-{M}_{\mathrm{molgas}}}\sim 30{M}_{\odot }$/L_⊙ from Zanella et al. (2018), this lower limit translates to L_[CII] ∼ 7 × 10⁷L_⊙. Therefore we start by integrating the [CII] luminosity function down to this limit giving us a luminosity density of ${{\mathscr{L}}}_{[\mathrm{CII}]}$ that is 1.2 × 10⁶L_⊙ Mpc⁻³. Zanella et al. (2018) find that the ${\alpha }_{[\mathrm{CII}]-{M}_{\mathrm{molgas}}}$ relation has an 1σ of 0.3 dex, translating into a linear ratio M_mol/L_[CII] = 15–60 M_⊙/L_⊙. We calculate the cosmic volume averaged molecular mass density ρ_molgas = (2–7) × 10⁷M_⊙ Mpc⁻³ at z ∼ 4–6. At a 2σ and 95% confidence level, we estimate ρ_molgas = (1–14) × 10⁷M_⊙ Mpc⁻³.

Note, however, that when using the Zanella et al. (2018) relation on individual ALPINE galaxies to convert their L_[CII] into the molecular gas mass, we see that the mean gas fraction is μ_molgas ≈ 3 (Dessauges-Zavadsky et al. 2020). This is slightly higher than in that of the Liu et al. (2019) sample, likely due to the fact that the ALPINE galaxies have lower stellar masses than the Liu et al. (2019) sample. We checked that changing this constant does not change our results in any substantive way (see the discussion below).

In Figure 9, we overlay our estimates of ρ_molgas on the Liu et al. (2019) model predictions (based on the A3COSMOS sample). We also overlay the 1σ and 2σ constraints from ASPECS (Decarli et al. 2019) and COLDz (Riechers et al. 2019), which surveys measured CO transitions using ALMA and JVLA and set constraints on the mass density of molecular gas at z ∼ 3–6. This cross-check is critical because all these surveys have pros and cons in their attempts to constrain the evolution of the cosmic molecular gas density. The key advantage of the ALPINE survey is its large sample (thanks to the [C II] line being a significantly brighter tracer of molecular gas than any CO line), but due to the uncertainties in the selection function we are left with only lower and upper bounds on the LF. The A3COSMOS-based constraint is also affected by selection biases. Both ASPECS and COLDz, while smaller samples, have the key advantage of being unbiased in terms of cold gas selection, since they use CO, which is the most reliable tracer of cold molecular gas mass currently known. Figure 9 shows that, within the uncertainties, all these studies are consistent. As such, these surveys provide supporting evidence that we are obtaining useful constraints as we use the sizeable ALPINE sample to constrain the cold gas density versus redshift relation.

Zoom In Zoom Out Reset image size

Figure 9 shows that our results are well within the range of the Liu et al. (2019) model, only slightly favoring the upper half. The spread in the model shown indicates a different assumption of the molecular gas fraction as a function of redshift. The upper limit is based on Scoville et al. (2017) analyses of the ALMA dust continuum observations at z < 3.8, and the lower limit is based on the Tacconi et al. (2018) relation, primarily from the PHIBBS CO survey at z < 2.5. In contrast to our results, the results of the CO surveys (COLDz and ASPECS) favor the lower half of the models shown. However, given the large uncertainties in all of these estimates, this slight tension between ours and the CO survey is not too significant. Within the 2σ uncertainties, our results are consistent with those from COLDz and ASPECS.

Dessauges-Zavadsky et al. (2020) also examined the molecular gas content of individual ALPINE galaxies at z ∼ 5, using [C II] detections to constrain the total molecular gas mass at z ∼ 5. Like us, they use the scaling relations from Zanella et al. (2018). They find μ_molgas ≈ 3, which leads to slightly lower total mass density and a better agreement with the Tacconi et al. (2018) model (see Figure 7 of their paper). Indeed, given the large uncertainties on these molecular gas mass densities, we cannot really exclude any of the models. The main difference between our analyses, however, is that our calculation is based on our best estimate of the [CII] luminosity function, including all three ALPINE constraints thereof (the UV-sample, the serendipitous sample, and the continuum-sample) as well as the far-IR LF from Koprowski et al. (2017). In this paper we do not look at the ALPINE constraints on the star formation rate density at z ∼ 5. This question is examined in our companion paper, Loiacono et al. (2020).