The ALMA Spectroscopic Survey in the HUDF: CO luminosity functions and the molecular gas content of galaxies through cosmic history

TL;DR

This study uses ALMA data from the HUDF to measure how the CO luminosity function and molecular gas content of galaxies have evolved over cosmic time, revealing a peak in molecular gas density at z~1.5 and a close link to star formation history.

Contribution

It provides the first comprehensive measurement of the CO luminosity function evolution and cosmic molecular gas density from z>1 using deep ALMA observations.

Findings

CO luminosity function evolves with redshift, decreasing in luminosity by an order of magnitude from z~2 to present.

The majority of cosmic CO luminosity at z=1-3 is recovered in the observations.

Cosmic molecular gas density peaks at z~1.5 and declines to the present, closely matching the star formation rate density evolution.

Abstract

We use the results from the ALMA large program ASPECS, the spectroscopic survey in the Hubble Ultra Deep Field (HUDF), to constrain CO luminosity functions of galaxies and the resulting redshift evolution of $\rho$(H$_2$). The broad frequency range covered enables us to identify CO emission lines of different rotational transitions in the HUDF at $z>1$. We find strong evidence that the CO luminosity function evolves with redshift, with the knee of the CO luminosity function decreasing in luminosity by an order of magnitude from $\sim$2 to the local universe. Based on Schechter fits, we estimate that our observations recover the majority (up to $\sim$90%, depending on the assumptions on the faint end) of the total cosmic CO luminosity at $z$=1.0-3.1. After correcting for CO excitation, and adopting a Galactic CO-to-H$_2$ conversion factor, we constrain the evolution of the cosmic molecular gas density $\rho$(H$_2$): this cosmic gas density peaks at $z\sim1.5$ and drops by factor of $6.5_{-1.4}^{+1.8}$ to the value measured locally. The observed evolution in $\rho$(H$_2$) therefore closely matches the evolution of the cosmic star formation rate density $\rho_{\rm SFR}$. We verify the robustness of our result with respect to assumptions on source inclusion and/or CO excitation. As the cosmic star formation history can be expressed as the product of the star formation efficiency and the cosmic density of molecular gas, the similar evolution of $\rho$(H$_2$) and $\rho_{\rm SFR}$ leaves only little room for a significant evolution of the average star formation efficiency in galaxies since $z\sim 3$ (85% of cosmic history).

Figures

Tables

Table 1. Table 1: CO transitions, redshift bins, cosmic volume, and typical H 2 mass limit [at 5- σ 𝜎 \sigma , assuming a line width of 200 km s -1 , CO excitation as in Daddi et al. 2015 , and a CO-to-H 2 conversion factor α CO subscript 𝛼 CO \alpha_{\rm CO} = 3.6 3.6 3.6 M ⊙ (K km s -1  pc 2 ) -1 ] in ASPECS LP 3 mm.

Line	Redshift	Volume	limit $M_{\mathrm{H2}}$
		[cMpc3]	[${10}^{10}$ M⊙]
(1)	(2)	(3)	(4)
CO(1-0)	$0.003-0.369$	338	0.11
CO(2-1)	$1.006-1.738$	8198	0.68
CO(3-2)	$2.008-3.107$	14931	1.8
CO(4-3)	$3.011-4.475$	18242	2.7

Table 2. Table 2: Results of the Schechter fits of the observed CO LFs, assuming a fixed α = − 0.2 𝛼 0.2 \alpha=-0.2 .

Line	log ${\Phi }_{\ast }$	log $L_{\ast }^{\prime }$
	[Mpc-3 dex-1]	[K km s-1 pc2]
(1)	(2)	(3)
All $L^{\prime }$ bins
CO(2-1)	$-{2.79}_{-0.09}^{+0.09}$	${10.09}_{-0.09}^{+0.10}$
CO(3-2)	$-{3.83}_{-0.12}^{+0.13}$	${10.60}_{-0.15}^{+0.20}$
CO(4-3)	$-{3.43}_{-0.22}^{+0.19}$	${9.98}_{-0.14}^{+0.22}$
Independent $L^{\prime }$ bins
CO(2-1)	$-{2.93}_{-0.12}^{+0.11}$	${10.23}_{-0.11}^{+0.16}$
CO(2-1)	$-{2.90}_{-0.14}^{+0.16}$	${10.22}_{-0.22}^{+0.24}$
CO(2-1)	$-{2.77}_{-0.20}^{+0.21}$	${10.12}_{-0.25}^{+0.35}$
CO(2-1)	$-{2.86}_{-0.14}^{+0.15}$	${10.17}_{-0.17}^{+0.17}$
CO(2-1)	$-{3.14}_{-0.19}^{+0.19}$	${10.32}_{-0.18}^{+0.26}$
CO(3-2)	$-{3.65}_{-0.23}^{+0.25}$	${10.49}_{-0.22}^{+0.26}$
CO(3-2)	$-{3.85}_{-0.20}^{+0.21}$	${10.59}_{-0.20}^{+0.23}$
CO(3-2)	$-{3.63}_{-0.17}^{+0.17}$	${10.36}_{-0.21}^{+0.25}$
CO(3-2)	$-{3.55}_{-0.26}^{+0.28}$	${10.22}_{-0.21}^{+0.19}$
CO(3-2)	$-{3.50}_{-0.21}^{+0.22}$	${10.24}_{-0.15}^{+0.21}$
CO(4-3)	$-{3.53}_{-0.28}^{+0.36}$	${10.01}_{-0.21}^{+0.26}$
CO(4-3)	$-{3.55}_{-0.26}^{+0.23}$	${10.10}_{-0.16}^{+0.18}$
CO(4-3)	$-{3.53}_{-0.19}^{+0.18}$	${10.08}_{-0.15}^{+0.20}$
CO(4-3)	$-{3.38}_{-0.28}^{+0.26}$	${9.98}_{-0.20}^{+0.36}$
CO(4-3)	$-{3.59}_{-0.23}^{+0.25}$	${10.21}_{-0.25}^{+0.40}$

Table 3. Table 3: Luminosity functions of the observed CO transitions. (1, 5) Luminosity bin center; each bin is 0.5 dex wide. (2-4, 6-8) CO luminosity functions, reported as the minimum and maximum values of the confidence levels at 1, 2, and 3- σ 𝜎 \sigma .

log $L^{\prime }$	log $\Phi$, 1-$\sigma$	log $\Phi$, 2-$\sigma$	log $L^{\prime }$	log $\Phi$, 1-$\sigma$	log $\Phi$, 2-$\sigma$
[K km s-1 pc2]	[dex-1 cMpc-3]	[dex-1 cMpc-3]	[K km s-1 pc2]	[dex-1 cMpc-3]	[dex-1 cMpc-3]
(1)	(2)	(3)	(4)	(5)	(6)
	CO(1-0)			CO(2-1)
8.0	-2.86 -1.72	-3.48 -1.48	9.4	-2.63 -2.37	-2.75 -2.27
8.1	-2.86 -1.72	-3.48 -1.48	9.5	-2.54 -2.31	-2.66 -2.22
8.2	-2.64 -1.65	-3.17 -1.43	9.6	-2.58 -2.33	-2.69 -2.24
8.3	-3.88 -1.79	-4.75 -1.51	9.7	-2.55 -2.31	-2.67 -2.22
8.4	-3.88 -1.79	-4.75 -1.51	9.8	-2.69 -2.41	-2.83 -2.30
8.5	-3.16 -1.81	-3.79 -1.55	9.9	-3.01 -2.61	-3.21 -2.47
8.6	-3.15 -1.81	-3.78 -1.55	10.0	-3.40 -2.81	-3.72 -2.64
8.7	-3.65 -1.91	-4.23 -1.61	10.1	-3.27 -2.75	-3.56 -2.59
8.8	-3.65 -1.91	-4.23 -1.61	10.2	-3.70 -2.93	-4.16 -2.73
8.9	-3.65 -1.91	-4.23 -1.61	10.3	-3.76 -2.95	-4.25 -2.75
9.0	-5.52 -1.97	-6.39 -1.65	10.4	-3.76 -2.95	-4.25 -2.75
			10.5	-3.76 -2.95	-4.25 -2.75
	CO(3-2)			CO(4-3)
9.6	-3.95 -3.21	-4.31 -3.01	9.6	-3.56 -3.09	-3.78 -2.93
9.7	-4.03 -3.24	-4.41 -3.03	9.7	-3.51 -3.06	-3.72 -2.91
9.8	-3.82 -3.14	-4.16 -2.96	9.8	-3.46 -3.04	-3.65 -2.89
9.9	-3.82 -3.14	-4.16 -2.96	9.9	-3.68 -3.15	-3.91 -2.99
10.0	-3.82 -3.14	-4.16 -2.96	10.0	-4.04 -3.34	-4.26 -3.13
10.1	-4.51 -3.34	-5.20 -3.11	10.1	-4.04 -3.34	-4.26 -3.13
10.2	-3.74 -3.10	-4.10 -2.92	10.2	-4.17 -3.40	-4.35 -3.18
10.3	-4.02 -3.21	-4.51 -3.01	10.3	-5.19 -3.59	-6.06 -3.31
10.4	-4.02 -3.21	-4.51 -3.01

Table 4. Table 4: Inferred CO(1-0) luminosity functions in various redshift bins. (1, 4) Luminosity bin center; each bin is 0.5 dex wide. (2-3, 5-6) CO luminosity functions, reported as the minimum and maximum values of the confidence levels at 1- and 2- σ 𝜎 \sigma .

log $L^{\prime }$	log $\Phi$, 1-$\sigma$	log $\Phi$, 2-$\sigma$	log $L^{\prime }$	log $\Phi$, 1-$\sigma$	log $\Phi$, 2-$\sigma$
[K km s-1 pc2]	[dex-1 cMpc-3]	[dex-1 cMpc-3]	[K km s-1 pc2]	[dex-1 cMpc-3]	[dex-1 cMpc-3]
(1)	(2)	(3)	(4)	(5)	(6)
8.0	-2.86 -1.72	-3.48 -1.48	9.5	-2.66 -2.40	-2.79 -2.30
8.1	-2.86 -1.72	-3.48 -1.48	9.6	-2.56 -2.32	-2.68 -2.23
8.2	-2.64 -1.65	-3.17 -1.43	9.7	-2.57 -2.33	-2.69 -2.23
8.3	-3.88 -1.79	-4.75 -1.51	9.8	-2.60 -2.35	-2.72 -2.25
8.4	-3.88 -1.79	-4.75 -1.51	9.9	-2.71 -2.42	-2.85 -2.32
8.5	-3.16 -1.81	-3.79 -1.55	10.0	-2.94 -2.57	-3.12 -2.44
8.6	-3.15 -1.81	-3.78 -1.55	10.1	-3.20 -2.72	-3.46 -2.56
8.7	-3.65 -1.91	-4.23 -1.61	10.2	-3.31 -2.77	-3.60 -2.61
8.8	-3.65 -1.91	-4.23 -1.61	10.3	-3.42 -2.82	-3.75 -2.65
8.9	-3.65 -1.91	-4.23 -1.61	10.4	-3.55 -2.88	-3.94 -2.69
9.0	-5.52 -1.97	-6.39 -1.65	10.5	-3.74 -2.94	-4.22 -2.74
			10.6	-3.91 -3.01	-4.43 -2.79
			10.7	-4.44 -3.20	-5.00 -2.93
			10.8	-5.71 -3.34	-6.54 -3.03
			10.9	-6.78 -3.35	-7.65 -3.04
10.0	-3.99 -3.23	-4.36 -3.02	10.1	-3.58 -3.10	-3.79 -2.94
10.1	-3.97 -3.22	-4.33 -3.02	10.2	-3.52 -3.07	-3.73 -2.92
10.2	-3.96 -3.21	-4.32 -3.01	10.3	-3.55 -3.08	-3.76 -2.93
10.3	-3.90 -3.18	-4.26 -2.99	10.4	-3.63 -3.13	-3.85 -2.97
10.4	-4.07 -3.24	-4.50 -3.03	10.5	-3.80 -3.22	-4.03 -3.04
10.5	-4.29 -3.30	-4.81 -3.08	10.6	-4.01 -3.33	-4.23 -3.12
10.6	-4.15 -3.26	-4.65 -3.05	10.7	-4.27 -3.43	-4.49 -3.20
10.7	-4.01 -3.21	-4.46 -3.01	10.8	-4.63 -3.54	-4.85 -3.28
10.8	-3.96 -3.19	-4.40 -3.00	10.9	-5.24 -3.64	-5.50 -3.35
10.9	-4.01 -3.20	-4.48 -3.01	11.0	-5.86 -3.69	-6.11 -3.38

Table 5. Table 5: Constraints on ρ 𝜌 \rho (H 2 ) in various redshift bins. The quoted ranges correspond to the 1- σ 𝜎 \sigma and 2- σ 𝜎 \sigma confidence levels in our analysis. We provide our reference estimates based on the whole sample, and assuming intermediate CO excitation (Daddi et al., 2015 ) (see Fig.  9 ), as well as the estimates for the secure sources only (Fig. B 10 ) and for the whole sample, but using different assumptions for the CO excitation (Fig. B 11 ).

Redshift	log $\rho$(H2), 1-$\sigma$	log $\rho$(H2), 2-$\sigma$
bin	[M⊙ Mpc-3]	[M⊙ Mpc-3]
(1)	(2)	(3)
Reference estimate
0.003–0.369	5.89–6.80	5.40–7.01
1.006–1.738	7.74–7.96	7.63–8.05
2.008–3.107	7.50–7.96	7.26–8.10
3.011–4.475	7.20–7.62	6.97–7.77
Secure sources only
0.003–0.369	5.13–6.41	4.25–6.65
1.006–1.738	7.64–7.90	7.51–7.99
2.008–3.107	7.39–7.96	7.08–8.12
3.011–4.475	6.71–7.35	6.37–7.53
Thermalized CO excitation
0.003–0.369	5.89–6.80	5.40–7.02
1.006–1.738	7.59–7.81	7.47–7.90
2.008–3.107	7.03–7.48	6.79–7.63
3.011–4.475	6.70–7.11	6.48–7.25
Milky Way–like CO excitation
0.003–0.369	5.91–6.82	5.43–7.04
1.006–1.738	7.90–8.12	7.79–8.20
2.008–3.107	7.65–8.09	7.41–8.24
3.011–4.475	7.33–7.81	7.08–7.96

Table 6. Table 6: Example of the line candidates entrying one of the realizations of the CO LFs. (1–2) Sky coordinates of the line candidate; (3) adopted CO–based redshift; (4) signal–to–noise; (5) completeness (see GL19); (6) does the line candidate have a counterpart at optical/NIR wavelengths with matching redshift (see text)? (7) fidelity the line candidate (see GL19); (8) inferred line luminosity; (9) rotational quantum number of the upper energy level of the transition.

RA	Dec	$z_{\mathrm{CO}}$	S/N	Compl.	c/p?	fid.	$L^{\prime }$	$J_{\mathrm{up}}$
$[$deg$]$	$[$deg$]$						[K km s-1 pc2]
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)
53.16063	-27.77626	2.5436	36.18	1.00	Y	1.00	$2.69\times {10}^{10}$	3
53.17664	-27.78551	1.3168	17.50	1.00	Y	1.00	$8.77\times {10}^9$	2
53.17086	-27.77545	2.4534	15.25	1.00	Y	1.00	$1.02\times {10}^{10}$	3
53.14350	-27.78324	1.4144	14.74	1.00	Y	1.00	$1.78\times {10}^{10}$	2
53.16569	-27.76991	1.5502	13.98	1.00	Y	1.00	$1.82\times {10}^{10}$	2
53.16616	-27.78754	1.0952	11.98	1.00	Y	1.00	$6.17\times {10}^9$	2
53.18138	-27.77756	2.6956	9.95	1.00	Y	1.00	$2.71\times {10}^{10}$	3
53.14822	-27.77389	1.3822	9.15	1.00	N	1.00	$3.17\times {10}^9$	2
53.17908	-27.78062	1.0365	8.74	1.00	Y	1.00	$7.85\times {10}^9$	2
53.15085	-27.77440	1.3827	7.44	1.00	N	1.00	$3.64\times {10}^9$	2
53.16583	-27.78157	1.0964	7.31	1.00	Y	1.00	$1.97\times {10}^9$	2
53.14817	-27.78451	3.6013	7.12	0.96	Y	1.00	$5.59\times {10}^9$	4
53.14523	-27.77801	1.0985	7.10	0.85	Y	1.00	$5.08\times {10}^9$	2
53.15199	-27.77552	1.0963	6.78	1.00	Y	1.00	$3.56\times {10}^9$	2
53.16635	-27.76873	1.2942	6.11	1.00	Y	0.96	$4.57\times {10}^9$	2
53.17966	-27.78428	0.1129	5.96	1.00	N	0.95	$1.04\times {10}^8$	1
53.15192	-27.78900	1.4953	5.91	1.00	Y	0.95	$5.14\times {10}^9$	2
53.14544	-27.78721	1.1746	5.78	0.87	Y	0.93	$3.08\times {10}^9$	2
53.16513	-27.76394	1.1769	5.77	1.00	N	0.75	$5.15\times {10}^9$	2
53.15104	-27.78691	1.7009	5.72	0.95	Y	0.87	$4.04\times {10}^9$	2
53.14553	-27.77757	3.6038	5.66	0.92	Y	0.61	$5.25\times {10}^9$	4
53.15495	-27.78709	1.0341	5.63	1.00	N	0.36	$3.39\times {10}^9$	2
53.14659	-27.77822	1.3821	5.56	0.93	N	0.76	$3.25\times {10}^9$	2
53.15702	-27.78166	1.1298	5.52	1.00	N	0.72	$2.56\times {10}^9$	2
53.17554	-27.78809	1.3835	5.46	0.46	N	0.74	$4.05\times {10}^9$	2
53.16848	-27.76772	1.2615	5.42	0.92	Y	0.49	$3.71\times {10}^9$	2
53.16946	-27.79258	0.1428	5.26	1.00	N	0.55	$1.03\times {10}^8$	1
53.16465	-27.79427	1.0122	5.26	1.00	N	0.52	$4.60\times {10}^9$	2
53.14334	-27.78797	3.1259	5.25	1.00	Y	0.38	$5.31\times {10}^9$	4
53.14437	-27.77806	0.1873	5.23	1.00	N	0.54	$2.80\times {10}^8$	1
53.16572	-27.79701	3.2263	5.23	0.76	N	0.55	$3.33\times {10}^9$	4
53.17748	-27.78064	1.1530	5.20	0.96	Y	0.53	$4.32\times {10}^9$	2
53.17554	-27.77674	1.2776	5.20	1.00	N	0.34	$2.14\times {10}^9$	2
53.14444	-27.78346	2.2128	5.12	1.00	N	0.55	$5.83\times {10}^9$	3
53.16161	-27.77591	1.2456	5.09	1.00	N	0.24	$2.24\times {10}^9$	2
53.17350	-27.79211	1.3310	5.05	0.93	N	0.28	$3.61\times {10}^9$	2
53.14514	-27.79452	1.0398	4.97	1.00	N	0.45	$4.33\times {10}^9$	2
53.14967	-27.78415	1.5710	4.93	0.97	Y	0.39	$3.78\times {10}^9$	2
53.14690	-27.78514	3.5460	4.92	1.00	N	0.50	$4.15\times {10}^9$	4
53.17144	-27.76966	2.2103	4.92	0.96	N	0.36	$6.23\times {10}^9$	3
53.14261	-27.78733	1.4269	4.90	1.00	Y	0.45	$5.03\times {10}^9$	2

Equations

$$\frac{L^{\prime }}{\mathrm{K}\mathrm{km}{\mathrm{s}}^{-1}\mathrm{p}{\mathrm{c}}^2}=\frac{3.257\times 10^7}{1+z}\frac{F_{\mathrm{line}}}{\mathrm{Jy}\mathrm{km}{\mathrm{s}}^{-1}}{\left(\frac{{\nu }_0}{\mathrm{GHz}}\right)}^{-2}{\left(\frac{D_{\mathrm{L}}}{\mathrm{Mpc}}\right)}^2$$

$$M_{H2}=\frac{{\alpha }_{\mathrm{CO}}}{r_{J1}}{L}^{\prime }$$

$$\log \Phi (L^{\prime })=\log {\Phi }_{\ast }+\alpha \log \left(\frac{L^{\prime }}{L_{\ast }^{\prime }}\right)\frac{1}{\ln 10}\frac{L^{\prime }}{L_{\ast }^{\prime }}+\log (\ln (10))$$

Full text

galaxies: high-redshift — galaxies: ISM — galaxies: star formation

1 Introduction

The molecular phase of the interstellar medium (ISM) is the birthplace of stars, and therefore it plays a central role in the evolution of galaxies (see reviews in Kennicutt & Evans, 2012; Carilli & Walter, 2013; Bolatto et al., 2013). The cosmic history of star formation (see, e.g., Madau & Dickinson, 2014), i.e., the mass of stars formed per unit time in a cosmological volume (or cosmic star formation rate density, $\rho_{\rm SFR}$) throughout cosmic time, increased from early cosmic epochs up to a peak at $z$=1–3, and then declined by a factor $\sim$8 until the present day. This could be explained by a larger supply of molecular gas (the fuel for star formation) in high–$z$ galaxies; by physical properties of the gas, that could more efficiently form stars; or by a combination of both. The characterization of the content and properties of the molecular ISM in galaxies at different cosmic epochs is therefore fundamental to our understanding of galaxy formation and evolution.

The H2 molecule, the main constituent of molecular gas, is a poor radiator: it lacks rotational transitions, and the energy levels of vibrational lines are populated significantly only at relatively high temperatures ($T_{\rm ex}>500$ K) that are not typical of the cold, star–forming ISM (Omont, 2007). On the other hand, the carbon monoxide molecule, 12CO (hereafter, CO) is the second most abundant molecule in the universe. Thanks to its bright rotational transitions, it has been detected even at the highest redshifts ($z\sim 7$; e.g., Riechers et al., 2013; Venemans et al., 2017a; Strandet et al., 2017; Marrone et al., 2018). Redshifted CO lines are observed in the radio and millimeter (mm) transparent windows of the atmosphere, thus becoming accessible to facilities such as the Jansky Very Large Array (JVLA), the IRAM NOrthern Expanded Millimeter Array (NOEMA), and the Atacama Large Millimeter Array (ALMA). CO is therefore the preferred observational probe of the molecular gas content in galaxies at high redshift.

To date, more than 250 galaxies have been detected in CO at $z>1$, the majority of which are quasar host galaxies or sub-mm galaxies (see, Carilli & Walter 2013 or a review); gravitationally–lensed galaxies (e.g., Riechers et al. 2010, Harris et al. 2012, Dessauges-Zavadsky et al. 2015, ; Aravena et al. 2016c; Dessauges-Zavadsky et al. 2017, González-López et al. 2017); and (proto-)clusters of galaxies (e.g., Aravena et al. 2012, Chapman et al. 2015, Seko et al. 2016, Rudnick et al. 2017, Hayatsu et al. 2017, Lee et al. 2017, Hayashi et al. 2018, Miller et al. 2018, Oteo et al. 2018). The remainder are galaxies selected based on their stellar mass ($M_{\rm*}$), star formation rate (SFR), and/or optical/near-infrared colors (e.g., Daddi et al., 2010a, b; Tacconi et al., 2010, 2013, 2018; Genzel et al., 2010, 2011, 2015). These studies were instrumental in shaping our understanding of the interplay between molecular gas reservoirs and star formation in massive $z>1$ galaxies on and above the ‘main sequence’ of star-forming galaxies (Noeske et al., 2007; Elbaz et al., 2011). E.g., these galaxies are found to have high molecular gas fractions $M_{\rm H2}$/$M_{\rm*}$ compared to galaxies in the local universe. The depletion time, $t_{\rm dep}$=$M_{\rm H2}$/SFR, i.e., the time required to consume the entire molecular gas content of a galaxy at the present rate of star formation, is shorter in starburst galaxies than in galaxies on the main sequence (see, e.g., Silverman et al., 2015, 2018; Schinnerer et al., 2016; Scoville et al., 2017; Tacconi et al., 2018). However, by nature these targeted studies are potentially biased towards specific types of galaxies (e.g., massive, star-forming galaxies), and consequently might fail to capture the full diversity of gas-rich galaxies in the universe.

Spectral line scans provide a complementary approach. These are interferometric observations over wide frequency ranges, targeting ‘blank’ regions of the sky. Gas, traced mainly via CO lines, is searched for at any position and frequency, without pre-selection based on other wavelengths. This provides us with a flux–limited census of the gas content in well–defined cosmological volumes. The first molecular scan reaching sufficient depth to detect MS galaxies targeted a $\sim$1 arcmin2 region in the Hubble Deep Field North (HDF-N; Williams et al., 1996) using the IRAM Plateau de Bure Interferometer (PdBI; see Decarli et al. 2014). The scan resulted in the first redshift measurement for the archetypal sub-mm galaxy HDF 850.1 ($z$=$5.183$, see Walter et al. 2012), and in the discovery of massive ($>10^{10}$ M⊙) gaseous reservoirs associated with galaxies at $z\sim 2$, including one with no obvious optical/NIR counterpart (Decarli et al., 2014). These observations enabled the first, admittedly loose constraints on the CO luminosity functions (LFs) and on the cosmic density of molecular gas in galaxies, $\rho$(H2), as a function of redshift (Walter et al., 2014). The HDF-N was also part of a second large observing campaign using the JVLA, the COLDz project. This effort ($>300$ hr of observations) targeted a $\sim$48 arcmin2 area in the GOODS-North footprint (Giavalisco et al., 2004), and a $\sim 8$ arcmin2 region in COSMOS (Scoville et al., 2007), sampling the frequency range 30–38 GHz (Lentati et al., 2015; Pavesi et al., 2018). This exposed the CO(1-0) emission in galaxies at $z\approx 2.0$–$2.8$ and the CO(2-1) emission at $z\approx 4.9$–6.7. The unprecedentedly large area covered by COLDz resulted in the best constraints on the CO LFs at $z>2$ so far, especially at the bright end (Riechers et al., 2019).

In ALMA Cycle 2, we scanned the 3 mm and 1.2 mm windows (84–115 GHz and 212–272 GHz, respectively) in a $\sim$1 arcmin2 region in the Hubble Ultra Deep Field (HUDF; Beckwith et al., 2006). This pilot program, dubbed the ALMA Spectroscopic Survey in the HUDF (ASPECS; Walter et al., 2016), pushed the constraints on the CO LFs at high redshift towards the expected knee of the CO LFs (Decarli et al., 2016a). By capitalizing on the combination of the 3 mm and 1.2 mm data, and on the unparalleled wealth of ancillary information available in the HUDF, Decarli et al. (2016b) were able to measure CO excitation in some of the observed sources, and to relate the CO emission to other properties of the observed galaxies at various wavelengths. Furthermore, the collapsed 1.2 mm data cube resulted in the deepest dust continuum image ever obtained at these wavelengths ($\sigma$=13 $\mu$Jy beam*-1*), which allowed us to resolve $\sim 80$% of the cosmic infrared background (Aravena et al., 2016a). The 1.2 mm data were also exploited to perform a systematic search for [C ii] emitters at $z$=6–8 (Aravena et al., 2016b), as well as to constrain the IRX–$\beta$ relation at high redshift (Bouwens et al., 2016). Finally, the ASPECS Pilot provided first direct measurements of the impact of foreground CO lines on measurements of the cosmic microwave background fluctuations, which is critical for intensity mapping experiments (Carilli et al., 2016).

The ASPECS Pilot program was limited by the small area surveyed. Here we present results from the ASPECS Large Program (ASPECS LP). The project replicates the survey strategy of the ASPECS Pilot, but on a larger mosaic that covers most of the Hubble eXtremely Deep Field (XDF), the region of the HUDF where the deepest near-infrared data are available (Illingworth et al. 2013; Koekemoer et al. 2013; see Fig. 1). Here we present and focus on the ASPECS LP 3 mm data, which have been collected in ALMA Cycle 4. We discuss the survey strategy and observations, the data reduction, the ancillary dataset, and we use the CO detections from the 3 mm data to measure the CO LFs in various redshift bins, and to infer the cosmic gas density $\rho$(H2) as a function of redshift. In González-López et al. (2019) (hereafter, GL19), we present our search for line and continuum sources, and assess their reliability and completeness. Aravena et al. (2019) place the ASPECS LP 3 mm results in the context of the main sequence narrative. Boogaard et al. (2019) capitalize on the sensitive VLT/MUSE Integral Field Spectroscopy of the field, in order to address how our CO detections relate with the properties of the galaxies as inferred from rest-frame optical/UV wavelengths. Finally, Popping et al. (2019) compare the ASPECS LP 3 mm results to state-of-the-art predictions from cosmological simulations and semi-analytical models.

The structure of this paper is as follows: In Sec. 2, we present the survey strategy, the observations, and the data reduction. In Sec. 3 we summarize the ancillary information available for the galaxies in this field. In Sec. 4 we present the main results of this study, and in Sec. 5 we discuss our findings and compare them with similar works in the literature. Finally, in Sec. 6 we infer our conclusions.

Throughout this paper we adopt a $\Lambda$CDM cosmological model with $H_{0}=70$ km s*-1* Mpc*-1*, $\Omega_{\rm m}=0.3$ and $\Omega_{\Lambda}=0.7$ (consistent with the measurements by the Planck Collaboration, 2015). Magnitudes are reported in the AB photometric system. For consistency with the majority of the literature on this field, in our analysis, we adopt a Chabrier (2003) stellar initial mass function.

2 Observations and data processing

2.1 Survey design and observations

The ASPECS LP survey consists of a 150 hr program in ALMA Cycle 4 (Program ID: 2016.1.00324.L). ASPECS LP comprises two scans, at 3 mm and 1.2 mm. The 3 mm survey presented here took 68 hr of telescope time (including calibrations and overheads), and was executed between December 2–21, 2016 (ALMA Cycle 4).

These observations comprised 17 pointings covering most of the XDF (Illingworth et al. 2013; Koekemoer et al. 2013; see Fig. 1). The pointings were arranged in a hexagonal pattern, distanced by $26.4^{\prime\prime}$ (the half–width of the primary beam of ALMA 12m antennas at the high–frequency end of ALMA band 3), thus ensuring Nyquist sampling and spatially–homogeneous noise in the mosaic. For reference, the central pointing is centered at Right Ascension = 03:32:38.5 and Declination = -27:47:00 (J2000.0). The total area covered at the center of the frequency scan ($\approx 99.5$ GHz) with primary beam attenuation $<0.5$ is 4.6 arcmin2. The observing strategy capitalized on the fast slew of the ALMA antennas in order to fully cover the entire mosaic between each phase calibrator observation. The survey was executed with the array in a relatively compact (C40-3) configuration. Baselines ranged between 15 and 700 m. The quasar J0334-4008 was observed as a flux, bandpass, and pointing calibrator, while the quasar J0342-3007 served as phase calibrator. The observations were performed in 5 different frequency settings, covering the frequency range 84–115 GHz. This enables the observation of one or more CO lines over a wide range of redshifts (see Fig. 2). Lower and upper side bands of the frequency settings partially overlap in the central part of the frequency range (96–103 GHz), thus yielding improved sensitivity at these frequencies (see also Fig. 3 in GL19).

2.2 Data reduction, calibration, and imaging

We processed the data using both the CASA pipeline for ALMA data (v. 4.7.0; McMullin et al., 2007) and our own procedures (see, e.g., Aravena et al., 2016a), which follow the general scheme of the official ALMA pipeline. Our independent inspection for data to be flagged allowed us to improve the depth of our scan in one of the frequency settings by up to 20%. In all the other frequency settings, the final rms appears consistent with the one computed from the cube provided by the ALMA pipeline. As the cube created with our own procedures is at least as good (in terms of low noise) as the one from the pipeline, we will refer to the former in the remainder of the analysis.

We imaged the 3 mm cube with natural weighting using the task tclean. The resulting synthesized beam is $\approx 1.75^{\prime\prime}\times 1.49^{\prime\prime}$ (PA=$91.5^{\circ}$) at the center of the observed frequency range. The lack of very bright sources in our cubes allows us to perform our analysis on the ‘dirty’ cube, thus preserving the intrisic properties of the noise. The resulting cube is used for the line search, and in most of the following analysis. In addition, we image the dataset after applying a u,v taper of $3^{\prime\prime}$ (using the uvtaper parameter in tclean) cleaned to 2-$\sigma$. This yields a reconstructed beam of $\sim 3.5^{\prime\prime}$. This latter cube is used only to extract the spectra of the sources identified in the search: thanks to the lower angular resolution, the spectra extracted in this way encapsulate all the emission of the sources, even in the case of sources that are spatially–resolved in the naturally–weighted imaging (see Aravena et al. 2019 for a discussion on the size of the CO emission in ASPECS LP 3 mm).

We rebin the frequency dimension in channels of 7.813 MHz, i.e., $2\times$ the native spectral resolution of the observations. At 99.5 GHz, this corresponds to $\Delta v\approx 23.5$ km s*-1*. We use ‘nearest’ interpolation scheme in order to maintain the independence of the channels despite the small frequency corrections due to the Earth rotation and revolution during the execution of the observations. We reach a sensitivity of $\sim 0.2$ mJy beam*-1* per 7.813 MHz channel throughout the scanned frequency range. For a line width of 200 km s*-1*, these limits correspond to fiducial 5-$\sigma$ CO line luminosity limits of $(1.4,2.1,2.3)\times 10^{9}$ K km s*-1* pc2, for CO(2-1), CO(3-2), and CO(4-3), respectively. Via the working assumptions discussed in section 4.3, we infer H2 mass limits of $6.8\times 10^{9}$ M⊙, $1.8\times 10^{10}$ M⊙, and $2.7\times 10^{10}$ M⊙ at $1.006<z<1.738$, $2.008<z<3.107$, and $3.011<z<4.475$ respectively. Fig. 3 compares these molecular gas mass limits and volume coverage reached in ASPECS LP 3 mm with those of all the other molecular scans performed so far. Tab. 1 lists the CO redshift coverage, fiducial gas mass limits, and the volume of universe of ASPECS LP 3 mm in various CO transitions.

3 Ancillary data

The HUDF is one of the best studied extragalactic regions in the sky. Our observations thus benefit from a wealth of ancillary data of unparalleled quality in terms of depth, angular resolution, wavelength coverage, and richness of spectroscopic information. When comparing with literature multi-wavelength catalogs, we apply a rigid astrometry offset ($\Delta$RA=$+0.076^{\prime\prime}$, $\Delta$Dec=$-0.279^{\prime\prime}$; see Rujopakarn et al. 2016; Dunlop et al. 2017) to available optical/NIR catalogs, in order to account for the different astrometric solution between the ALMA data and optical/NIR data.

The bulk of optical and NIR photometry comes from the Hubble Space Telescope (HST) Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey (CANDELS; Grogin et al. 2011; Koekemoer et al. 2011). These are based both on archival and new HST images obtained with the Advanced Camera for Surveys (ACS) at optical wavelengths, and with the Wide Field Camera 3 (WFC3) in the near-infrared. We refer to the photometric compilation by Skelton et al. (2014), which also includes ground–based optical and NIR photometry from Nonino et al. (2009), Hildebrandt et al. (2006), Erben et al. (2005), Cardamone et al. (2010), Wuyts et al. (2008), Retzlaff et al. (2010), Hsieh et al. (2012), as well as Spitzer IRAC 3.6 $\mu$m, 4.5 $\mu$m, 5.8 $\mu$m, and 8.0 $\mu$m photometry from Dickinson et al. (2003), Elbaz et al. (2011), and Ashby et al. (2013). We also include the Spitzer MIPS 24 $\mu$m photometric information from Whitaker et al. (2014).

The main optical spectroscopy sample in the ASPECS LP footprint comes from the MUSE Hubble Ultra Deep Survey (Bacon et al., 2017), a mosaic of nine contiguous fields observed with the Multi Unit Spectroscopic Explorer at the ESO Very Large Telescope. The surveyed area encompasses the entire HUDF. MUSE provides integral field spectroscopy of a $1^{\prime}\times 1^{\prime}$ square field over the wavelength range 4750–9300 Å. This yields emission–line redshift coverage in the ranges $z<0.857$, $0.274<z<1.495$, $1.488<z<3.872$, $2.906<z<6.648$ for [Oiii]${}_{5000\,\rm\AA}$, [Oii]${}_{3727\,\rm\AA}$, Ciii]${}_{1909\,\rm\AA}$, and Ly$\alpha$, respectively. The redshift catalog based on the MUSE Hubble Ultra Deep Survey consists of $>1500$ galaxies with spectroscopic redshifts in the HUDF (Inami et al., 2017). We also include any additional spectroscopic information based on various studies at optical and NIR wavelengths, as compiled in Le Fèvre et al. (2005), Coe et al. (2006), Skelton et al. (2014), and Morris et al. (2015).

HST grism spectroscopy is also available in the HUDF. These observations allow for integral field spectroscopy with sub-arcsec angular resolution at relatively modest (\lambda/\Delta\lambda\mathrel{\hbox to0.0pt{\lower 3.0pt\hbox{\sim}\hss}\raise 2.0pt\hbox{<}}1000) spectral resolution. While optical grism spectroscopy of the HUDF has been done (Xu et al., 2007), we take particular advantage of the more recent HST grism spectroscopy campaigns at NIR wavelengths, in particular the 3D-HST survey (Momcheva et al., 2016). This complements the MUSE information, providing spectroscopy of H$\alpha$, H$\beta$ and other rest-frame optical lines at $z=1-3$, together with some additional redshift information in the “redshift desert” at $1.5<z<2.9$ where MUSE is less efficient due to the paucity of bright emission lines that are shifted into the MUSE wavelength range at these redshifts.

We create a master catalog of galaxies in the HUDF by combining the Skelton et al. (2014) catalog with the compilations by Le Fèvre et al. (2005), Coe et al. (2006), Xu et al. (2007), Rhoads et al. (2009), McLure et al. (2013), Schenker et al. (2013), Bouwens et al. (2014, 2015), Morris et al. (2015), Inami et al. (2017). The catalogs are merged with a simple geometrical association, with an angular threshold of $0.5^{\prime\prime}$ ($1.0^{\prime\prime}$) for the photometry (spectroscopy). This selection is also cross-matched with the measurements of morphological parameters (size, ellipticity, light concentration index) from van der Wel et al. (2012). The whole catalog, extending over most of the GOODS–South footprint, consists of $>63000$ entries. In the $2.5^{\prime}\times 2.1^{\prime}$ area of the XDF, the catalog includes photometry in $>$30 broad and medium bands for $\sim 7000$ galaxies, 475 of which have a spectroscopic redshift.

The photometric dataset is modeled with the high–$z$ extension of the Spectral Energy Distribution (SED) fitting code MAGPHYS (da Cunha et al., 2008, 2015), in order to infer physical parameters: stellar mass, sSFR (and thus, SFR), dust extinction, IR luminosity, etc. We use the available photometry between $0.37$ $\mu$m and 8.0 $\mu$m, as well as data from the available 1.2 mm imaging of the field. These results are discussed in detail in Boogaard et al. (2019).

4 Analysis and Results

Our goal is to compute CO luminosity functions and measurements of $\rho$(H2) based on the results from the CO line search in the ASPECS LP 3 mm data. Our workflow, sketched in Fig. 4, is articulated in four main blocks: The search for line candidates in the cube, and their characterization in terms of observed quantites (e.g., line fluxes); the assessment of the reliability of the line candidates and of the completeness of our line search; the identification of the line candidates and the measurement of a CO–based redshift; and the construction of high–level data products (e.g., luminosity functions).

4.1 Line search

We extensively discuss the line search approach in GL19, and summarize the main steps here for completeness. The cube is searched for emission at any spatial position and spectral coordinate, without any prior based on data from other wavelengths, in order to minimize biases in our selection function. Among the compilations presented in GL19, here we refer to the results obtained with findclumps. This catalog of line candidates consists of 613 entries at S/N$>$5.0, 70 at S/N$>$5.5, 21 at S/N$>$6.0, and 15 at S/N$>$6.5.

The fidelity or reliability of a line candidate gauges the impact of false positive detections in our search. The idea is to estimate the probability that a given line candidate may be spurious (i.e., a noise feature). The statistics of negative line candidates is used to model the noise properties of the cube, as a function of the S/N and the width of each line candidate111Since we adopt a matched–filter approach in the line search, the line width affects the reliability of the line candidates in that the narrower the filter kernel, the more independent resolution elements are present in the cube. As a consequence, the probability of finding a high–S/N noise peak increases.. The fidelity is then defined as 1-$P$, where $P$ is the probability of a (positive) line candidate to be due to noise. We limit our analysis to line candidates with fidelity $>$20%. We discuss the impact of fidelity on our results in Section 5.

The completeness of our line search is estimated by ingesting in the cube mock lines spanning a range of values for various parameters (3D position in the cube, flux, width), under the assumption that the lines are well described by Gaussian profiles. The line search is then repeated, and the completeness is simply inferred as the ratio between the number of retrieved and ingested mock lines, as a function of all the input parameters. In the construction of the CO LFs, we only consider line candidates with a parameter set yielding a completeness $>$20%.

4.2 Line identification and redshifts

In order to convert the fluxes of the line candidates into luminosities, we need to identify the observed lines: In principle, the spectral range covered in our 3 mm scan is broad enough to encompass multiple CO transitions at specific redshifts, thus offering a robust direct constraint on the line identification. However, as shown in Fig. 2, this happens only at relatively high redshifts (z\mathrel{\hbox to0.0pt{\lower 3.0pt\hbox{\sim}\hss}\raise 2.0pt\hbox{>}}3, if one considers both CO and [C i]). We therefore need to consider different approaches to pin down the redshift of our line candidates. First, we search for a counterpart at optical/NIR wavelengths. If successful, we use the available redshift of the counterpart to associate line candidates and CO transitions: if the counterpart has a redshift $z_{\rm cat}<0.8$, $0.8<z_{\rm cat}<1.9$, $1.9<z<3.2$, etc, we identify the line candidate as CO(1-0), CO(2-1), CO(3-2), etc, respectively. The majority of the 21 line candidates with S/N$>$6 show very good agreement (|\delta z|=|(z_{\rm cat}-z_{\rm CO})/(1+z_{\rm CO})|\mathrel{\hbox to0.0pt{\lower 3.0pt\hbox{\sim}\hss}\raise 2.0pt\hbox{<}}0.01) between CO–based and catalog redshifts (see Fig. 5). Other line candidates have a CO redshift roughly consistent ($|\delta z|<0.3$) with the catalog photometric redshifts. Two galaxies detected at S/N$>$6 in CO have a spectroscopic catalog redshift that is inconsistent with the CO–based redshift. Our detailed analysis of the MUSE data confirms that these cases are examples of overlapping galaxies at different redshifts; i.e., both the catalog values and the CO–based values are confirmed (Boogaard et al., 2019). Fig. 5 shows the comparison between the CO–based and catalog redshifts.

If the line candidates do not have a counterpart at other wavelengths (about 25% of line candidates at S/N$>$5), the line identification is performed through a bootstrap, where the probability of a line candidate to be CO(1-0), CO(2-1), CO(3-2), and CO(4-3) is proportional to the volume of universe sampled in each of these transitions with ASPECS LP at 3 mm. We do not consider transitions at higher $J$ values, since significant CO excitation would have to be invoked in order to explain bright high–$J$ line emission. In Appendix B, we discuss the impact of these assumptions on our results.

In the construction of CO luminosity functions, we only use CO–based redshifts.

4.3 Line luminosities and corresponding H2 mass

The line fluxes are transformed into luminosities following Carilli & Walter (2013):

[TABLE]

where $F_{\rm line}$ is the integrated line flux, $\nu_{0}$ is the rest-frame frequency of the line, and $D_{\rm L}$ is the luminosity distance. We then infer the corresponding CO(1-0) luminosities by adopting the CO[$J$-($J$-1)]–to–CO(1-0) luminosity ratios, $r_{J1}$, from Daddi et al. (2015): $L^{\prime}$ [CO(1-0)] = $L^{\prime}/r_{J1}$, with $r_{J1}=\{1.00$, $0.76\pm 0.09$, $0.42\pm 0.07$, $0.31\pm 0.07\}$, for $J_{\rm up}$=$\{$1, 2, 3, 4$\}$. These values are based on VLA and PdBI observations of multiple CO transitions in 4 main sequence galaxies at $z\approx 1.5$. These galaxies are less extreme than the typical, high IR luminosity galaxies studied in multiple CO transitions at $z>1$, thus likely more representative of the galaxies studied here. We include a bootstrapped realization of the uncertainties on $r_{J1}$ in the conversion. In Appendix B we discuss the impact of the $r_{J1}$ assumptions on our results.

The cosmic microwave background at high redshift enhances the minimum temperature of the cold ISM, and suppresses the observability of CO lines in galaxies because of the lower constrast against the background (for extended discussions, see e.g., da Cunha et al., 2013; Tunnard & Greve, 2016). The net effect is that the observed CO emission is only a fraction of the intrinsic one, with the suppression being larger for lower J transitions and at higher redshifts. This correction is however typically small at $z=1-3$, and often neglected in the literature (e.g., Tacconi et al., 2018). Indeed, for $T_{\rm kin}\approx T_{\rm dust}$, and following the $T_{\rm dust}$ evolution in Magnelli et al. (2014), we find $T_{kin}>30$ K at $z>1$, thus yielding CO flux corrections of \mathrel{\hbox to0.0pt{\lower 3.0pt\hbox{\sim}\hss}\raise 2.0pt\hbox{<}}15% up to $z$=$4.5$. Because of its minimal impact, the associated uncertainties, and for consistency with the literature, we do not correct our measurements for the cosmic microwave background impact.

The resulting CO(1-0) luminosities are converted into molecular gas masses, $M_{\rm H2}$, via the assumption of a CO–to–H2 conversion factor, $\alpha_{\rm CO}$:

[TABLE]

A widespread assumption in the literature on “normal” high–redshift galaxies (e.g., Daddi et al., 2010a; Magnelli et al., 2012; Carilli & Walter, 2013; Tacconi et al., 2013, 2018; Genzel et al., 2015; Riechers et al., 2019) is a value of $\alpha_{\rm CO}\approx 4$ M⊙ (K km s*-1* pc2)-1, consistent with the Galactic value (see, e.g., Bolatto et al., 2013), once the Helium contribution ($\sim 36$%) is removed. Here we adopt $\alpha_{\rm CO}=3.6$ M⊙ (K km s*-1* pc2)-1 (Daddi et al., 2010a). A different, yet constant choise of $\alpha_{\rm CO}$ would result in a linear scaling of our results involving $M_{\rm H2}$ and $\rho$(H2). This is further discussed in Sec. 5.

4.4 CO luminosity functions

The CO luminosity functions are constructed in a similar way as in Decarli et al. (2016a) via a Monte Carlo approach that allows us to simultaneously account for all the uncertainties in the line flux estimates, in the line identification, in the conversion factors, as well as for the fidelity of the line candidates. For each line candidate, we compute the corresponding values of completeness and fidelity, based on the observed line properties (S/N, line width, flux, etc). If the line has been confirmed by, e.g., a counterpart with a matching spectroscopic redshift, we assume that the fidelity is 1. In all the other cases, we conservatively treat our fidelity estimates as upper limits, and adopt a random value of fidelity that is uniformly distributed between 0 and such upper limit (see GL19; Pavesi et al. 2018; Riechers et al. 2019). We extract a random number for each entry; line candidates are kept in our analysis only if the random value is below the fidelity threshold (so that, the lower the fidelity, the lower the chances that the line candidate is kept in our analysis). Typically, 20–40 line candidates survive this selection in each realization.

We split the list of line candidates by CO transitions and in 0.5 dex wide bins of luminosity. In each bin, we compute the Poissonian uncertainties. We then scale up each entry by the inverse of the completeness. The completeness–corrected entry counts in each bin are then divided by the comoving volume covered in each transition. This is computed by counting the area with sensitivity $>$50% of the peak sensitivity obtained at the center of the mosaic in each channel.

The CO luminosity functions are created 1000 times (both for the observed CO transitions, and for the corresponding J=1$\rightarrow$0 ground transition), each time with a different realization of all the parameters that are left uncertain (the fidelity and its error bars, the identification of lines without counterparts, the $r_{J1}$ ratio, etc). The analysis is then repeated five times after a shift of 0.1 dex of the luminosity bins, which allows us to remove the dependence of the reconstructed CO luminosity functions from the bin definition. The final CO luminosity functions are the averages of all the CO LF realizations. The CO and CO(1-0) LFs are listed in Tables 3 and 4 and plotted in Fig. 6 and 7.

The H2 mass functions in our analysis are simply obtained by scaling the CO(1-0) LFs by the (fixed) $\alpha_{\rm CO}$ factor. We then sum the CO–based completeness–corrected H2 masses of each line candidate passing the fidelity threshold in bins of redshift, and we divide by the comoving volume in order to derive the cosmic gas molecular mass density, $\rho$(H2). By construction, we do not extrapolate towards low CO luminosities / low H2 masses. However, in the following we will show that accounting for the faint end would only very marginally affect our results.

4.5 Analytical fits to the CO LFs

We fit the observed CO luminosity functions with a Schechter function (Schechter, 1976), in the logarithmic form used in Riechers et al. (2019):

[TABLE]

where $\Phi(L^{\prime})\,d(\log L^{\prime})$ is the number of galaxies per comoving volume with a CO line luminosity between log $L^{\prime}$ and log $L^{\prime}+d(\log L^{\prime})$; $\Phi_{*}$ is the scale number of galaxies per unit volume; $L^{\prime}_{*}$ is the scale line luminosity which sets the knee of the luminosity function; $\alpha$ is the slope of the faint end. We fit the observed CO LFs in the three redshift bins at $z>1$ considered in this study; the $z<0.37$ bin is ignored because of the modest luminosity range and sample size in our study. The LFs presented in this work are created in bins of 0.5 dex spaced by 0.1 dex, i.e., consecutive bins are not independent. In order to account for this, and to minimize the impact of our bin assumptions, we first fit the LFs using all the available bins, then we repeat the fits on the five independent contiguous subsets of the luminosity bins.

The slope of the faint end of the LF, $\alpha$, is very sensitive to the corrections we apply for fidelity and completeness (see previous section). We therefore opt to conservatively use a fiducial fixed value of $\alpha$=–0.2 in our analysis. This is consistent with findings at $z\approx 0$ (Saintonge et al., 2017, once we take into account the different definition of $\alpha$), as well as with the typical slope of the stellar mass function of field galaxies at various redshifts (e.g., Ilbert et al., 2013). As for the other two parameters, we assume broad ($\sigma$=$0.5$ dex) log normal distributions as priors in $\Phi_{*}$ and $L^{\prime}_{*}$, centered around $10^{-3}$ Mpc*-3* dex*-1*, and $10^{9.5}$ K km s*-1* pc2, respectively. The best fit value and the 1-$\sigma$ confidence levels of the fitted parameters are derived from the 50%, 14%, and 86% quartiles of the marginalized posterior distributions of each parameter. They are listed in Table 2. Fig. 8 compares the observed CO LFs with the fitted Schechter functions.

The fitted parameters do not show strong dependency on the choice of binning, with the results being typically consistent within 1-$\sigma$ uncertainties. We find an indication of a higher $L^{\prime}_{*}$ at $z$=1–3 (log $L^{\prime}_{*}$ [K km s*-1* pc2] $\approx 10.4$) compared to the $z>3$ bin, and most importantly, with the local universe (log $L^{\prime}_{*}$ [K km s*-1* pc2] $\approx 9.9$, although with a different definition of the Schechter function; Saintonge et al. 2017).

The luminosity–weighted integral of the fitted LFs suggests that the ASPECS LP 3 mm data recover 83%, 91%, and 71% of the total CO luminosity at $\langle z\rangle$ = 1.43, 2.61, and at 3.80. In addition, if we adopt the best fit by Saintonge et al. (2017) for the lowest redshift bin, ASPECS LP 3 mm recovers 59% of the total CO(1-0) luminosity in the local universe, although this last measure is strongly affected by cosmic variance due to the small volume probed by ASPECS LP 3 mm.

5 Discussion

Figs. 6 and 7 show that ASPECS LP 3 mm sampled a factor $\sim 20$ in CO luminosity at $z=1-4$ (see also Tables 3 and 4). We find evidence of an evolution in the CO LFs [and in the corresponding CO(1-0) LFs] as a function of redshift, compared to the local universe (Keres et al., 2003; Boselli et al., 2014; Saintonge et al., 2017), suggesting that the characteristic CO luminosity of galaxies at $z$=1–4 is an order of magnitude higher than in the local universe, once we account for CO excitation. This is in line with the findings from other studies, e.g., other molecular scans (Walter et al., 2014; Decarli et al., 2016a; Riechers et al., 2019); targeted CO observations on large samples of galaxies (e.g., Genzel et al., 2015; Aravena et al., 2016c; Tacconi et al., 2018); and similar works based on dust continuum observations (e.g., Magnelli et al., 2013; Gruppioni et al., 2013; Scoville et al., 2017). The CO LFs show an excess at the bright end compared with the predictions by semi–analytical models (Lagos et al., 2011; Popping et al., 2014), and more compatible with empirical predictions (Sargent et al., 2014; Vallini et al., 2016). Fig. 9 demonstrates that a prominent evolution in $\rho$(H2) occurred between $z\approx 4$ and nowadays, with the molecular gas content in galaxies slowly rising since early cosmic epochs, peaking around $z$=1–3, and dropping by a factor $6.5_{-1.4}^{+1.8}$ down to the present age (see also Table 5). The values of $\rho$(H2) used here only refer to the actual line candidates, i.e., we do not attempt to extrapolate towards undetected faint end of the LFs. However, as discussed in Sec. 4.5, our observations recover close to 90% of the total CO luminosity at $z=1.0$–$3.1$ (under the assumption of a slope of $\alpha$=–$0.2$ for the faint end), i.e., the derived $\rho$(H2) values would shift upwards by small factors ($\sim 10-20$%). In Appendix B, we test the robustness of the CO LFs and $\rho$(H2) evolution with redshift against some of the working assumptions in our analysis. A different choice of $\alpha_{\rm CO}$ would linearly affect our results on $\rho$(H2). In particular, by adopting $\alpha_{\rm CO}\approx 2$ M⊙ (K km s*-1* pc2)-1, as the comparison between dust–based and CO–based gas masses suggests (Aravena et al., 2019), we would infer a milder evolution of $\rho$(H2) at $z>1$ and the local measurements.

In the following, we discuss our results in the context of previous studies.

5.1 CO Luminosity functions

Compared to any previous molecular scan at mm wavelengths (Walter et al., 2014; Decarli et al., 2016a; Riechers et al., 2019), ASPECS LP 3 mm provides superior sample statistics, which enables the more detailed analysis described in this series of papers. As shown in Fig. 3, ASPECS LP 3 mm complements very well COLDz in that it samples a smaller volume but reaching a deeper sensitivity. The large volumes sampled by ASPECS LP 3 mm and COLDz, and the different targeted fields, mitigate the impact of cosmic variance. Overall, the CO LFs observed from the ASPECS LP 3 mm data appear in good agreement with the constraints from the first molecular scan observations (see Fig. 6).

Fig. 6 compares our observed CO LFs with the LF predictions by the semi-analytical models presented in Lagos et al. (2011) and Popping et al. (2014). Semi-analytical models tend to underpredict the bright end of the CO luminosity functions at $z>1$, with a larger discrepancy for the Lagos et al. (2011) models at $z>2$. The tension increases if we compare our inferred CO(1-0) LFs with the predictions from models (Fig. 7). This hints at an intrinsic difference on how widespread large molecular gas reservoirs in high–redshift galaxies are as predicted by models, with respect to that suggested by our observations; the tension is somewhat reduced by the different treatment of the CO excitation (see Appendix B).

The CO(1-0) LF inferred in our study at $2.0<z<3.1$ is in excellent agreement with the one derived from the COLDz survey (Riechers et al., 2019) (see Fig. 7). Because of the different parameter space, the COLDz data show larger uncertainties in the faint end, but provide a better constraint on the bright end compared to ASPECS LP 3 mm. The good match between these two independent observations might be considered as supporting evidence that the impact of cosmic variance is relatively modest (the targeted fields are completely independent); and that our assumption on the CO excitation, used to transform CO(3-2) into its corresponding CO(1-0), works reasonably well. Interestingly, the CO(1-0) LFs from ASPECS LP 3 mm and COLDz derived at $z\sim 2.5$ are in good agreement with the empirical predictions by Vallini et al. (2016), based on the Herschel IR luminosity functions (Gruppioni et al., 2013).

It is also interesting to compare the inferred CO(1-0) LFs with the ones measured in the local universe by Keres et al. (2003), Boselli et al. (2014), and Saintonge et al. (2017). The local measurements differ from each other by up to a factor $\sim 2$. Nevertheless, the ASPECS LP 3 mm data show a very clear evolution in the CO(1-0) LFs, with a shift upward of the knee of the luminosity functions by an order of magnitude, or an excess by several orders of magnitudes in the number density of bright ($L^{\prime}>10^{10}$ K km s*-1* pc2) CO(1-0) emitters at $z>1$ compared to the local universe.

5.2 $\rho$(H2) vs redshift

Fig. 9 compares the observed evolution of $\rho$(H2) as a function of redshift from the available molecular scan efforts. The ASPECS LP 3 mm data confirm the results from the PdBI scan in the HDF-N (Walter et al., 2014) and from the ASPECS Pilot (Decarli et al., 2016a), but with much tighter constraints thanks to the superior statistics. The cosmic density of molecular gas in galaxies appears to increase by a factor $6.5_{-1.4}^{+1.8}$ from the local universe [$\rho$(H2)$\approx$1.1$\times$107 M⊙ Mpc*-3*; Keres et al. 2003, Boselli et al. 2014, Saintonge et al. 2017] to $z\sim 1$ [$\rho$(H2)$\approx$7.1$\times$107 Mpc*-3*]222The ASPECS LP 3 mm constraint on $\rho$(H2) at $z<0.3$ is below the estimates from local studies; this is likely due to the fact that we only sampled a small luminosity range in $L^{\prime}_{\rm CO}$ in this redshift bin, in a tiny cosmic volume; furthermore, the HUDF was originally chosen to be relatively underdense of nearby galaxies., then follows a relatively flat evolution or possibly a mild decline towards higher redshifts. This is in excellent agreement with the constraints on $\rho$(H2) from COLDz (Riechers et al., 2019) at $2.0<z<2.8$, and with the empirical predictions derived by Sargent et al. (2014) based on the “2–star formation mode” framework, where the distributions of various galaxy properities (gas fraction, star formation efficiency, metallicity, etc) are inferred based on empirical relations, with a key role due to the offset of galaxies with respect to the “main sequence”. This analysis results in a similar evolution of $\rho$(H2) with redshift as the one found in ASPECS LP 3 mm.

The observed evolution in $\rho$(H2) is also in qualitative agreement with other observational studies. E.g., most studies searching for CO emission in targeted observations of main sequence galaxies find that $z$=1–3 galaxies typically have 5–10 times more gas than galaxies of similar stellar mass in the local universe (see, e.g., Genzel et al., 2015; Schinnerer et al., 2016; Tacconi et al., 2018). This is in line with the ASPECS LP 3 mm results, which point to a larger molecular gas content in typical galaxies at $z>1$ (see also Aravena et al., 2019). A similar trend is also reported by studies tackling the problem using dust as a probe of the gas content in high–$z$ galaxies (e.g. Magnelli et al., 2013; Gruppioni et al., 2013). E.g., Scoville et al. (2017) put indirect constraints on $\rho$(H2) at various redshifts using dust continuum measurements of Herschel–selected galaxies, and scaling by an internally–calibrated dust–to–gas ratio. The evolution of $\rho$(H2) that Scoville et al. (2017) infer is qualitatively similar, although somewhat shallower than the one observed in ASPECS LP 3 mm, spanning only a factor $\sim 2.5$ in $\rho$(H2) compared to a factor $\sim 6.0$ found in ASPECS LP 3 mm.

6 Conclusions

We presented the ASPECS LP 3 mm survey, an ALMA molecular scan encompassing most of the Hubble XDF over a large fraction of the 3 mm transparent band of the atmosphere. We exploited our data to search for massive molecular gas reservoirs (as traced by CO emission) in galaxies throughout $\sim 90$% of cosmic history, with no prior on counterparts at other wavelengths. We detected 70 line candidates with S/N$>$5.5, $>$75% of which with a photometric counterpart at optical/NIR wavelengths. This search allowed us to put stringent constraints on the CO luminosity functions in various redshift bins, as well as to infer the cosmic density of molecular gas in galaxies, $\rho$(H2). We found that:

• i-

The CO luminosity functions undergo significant evolution compared to the local universe. High redshift galaxies appear brighter in CO than galaxies in the local universe. In particular, at $z$=1–3, the characteristic CO(2-1) and CO(3-2) luminosity is $>3\times$ higher than the characteristic CO(1-0) luminosity observed in the local universe. The evolution is even stronger if we account for CO excitation. Analytical fits of our results suggests that we recovered the majority (up to 90%, depending on assumptions on the faint end) of the total CO luminosity at $z$=1.0–3.1.

• ii-

Similarly, $\rho$(H2) shows a clear evolution with cosmic time: It slowly increased since early cosmic epochs, reached a peak around $z$=1–3, and then decreased by a factor $6.5_{-1.4}^{+1.8}$ to the present day. This factor changes if $\alpha_{\rm CO}$ is allowed to evolve with redshift. In particular, the factor would be $\sim$3 if we adopt $\alpha_{\rm CO}$=2 M⊙ (K km s*-1* pc2)-1 for galaxies at $z>1$.

• iii-

Our results are in agreement with those of other molecular scans which targeted different regions of the sky and sampled different parts of the parameter space (in terms of depth, volume, transitions, etc). Similarly, we generally confirm empirical predictions based on dust continuum observations and SED modeling.

• iv-

Our results are in tension with predictions by semi-analytical models, which struggle to reproduce the bright end of the observed CO LFs. The discrepancy might be mitigated with different assumptions on the CO excitation and $\alpha_{\rm CO}$. Popping et al. (2019) quantitatively address the comparison between models and the ASPECS LP 3 mm observations and the underlying assumptions of both.

• v-

Our results hold valid if we restrict our analysis to the subset of galaxies with counterparts at redshifts that strictly match those inferred from our CO observations. The results are qualitatively robust against different assumptions concerning the CO excitation.

The observed evolution of $\rho$(H2) is in quantitative agreement with the evolution of the cosmic star formation rate density ($\rho_{\rm SFR}$; see, e.g., Madau & Dickinson, 2014), which also shows a mild increase up to $z$=1–3, followed by a drop by a factor $\approx$8 down to present day. Given that the star formation rate can be expressed as the product of the star formation efficiency (= star formation per unit gas mass) and the gas content mass, the similar evolution of $\rho$(H2) and $\rho_{\rm SFR}$ leaves little room for a significant evolution of the star formation efficiency throughout 85% of cosmic history ($z\approx 3$), at least when averaged over the entire galaxy population. The history of cosmic star formation appears dominated by the evolution in the molecular gas content of galaxies.

We thank the anonymous referee for their useful feedback which allowed us to improve the quality of the paper. Este trabajo contó con el apoyo de CONICYT + Programa de Astronomía+ Fondo CHINA-CONICYT. J.G.L. acknowledges partial support from ALMA-CONICYT project 31160033. D.R. acknowledges support from the National Science Foundation under grant number AST-1614213. F.E.B. acknowledges support from CONICYT-Chile Basal AFB-170002 and the Ministry of Economy, Development, and Tourism’s Millennium Science Initiative through grant IC120009, awarded to The Millennium Institute of Astrophysics, MAS. I.R.S. acknowledges support from the ERC Advanced Grant DUSTYGAL (321334) and STFC (ST/P000541/1). T.D-S. acknowledges support from ALMA-CONICYT project 31130005 and FONDECYT project 1151239. J.H. acknowledges support of the VIDI research programme with project number 639.042.611, which is (partly) financed by the Netherlands Organisation for Scientific Research (NWO). ALMA data: 2016.1.00324.L. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), NSC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ.

Appendix A Measured CO Luminosity Functions

For the sake of reproducibility of our results, Table 3 reports the measured CO LFs in ASPECS LP 3 mm. Similarly, Table 4 provides the inferred CO(1-0) LFs from this study. Table 5 lists the estimated values of $\rho$(H2) in various redshift bins and under different working hypothesis (see Appendix B). Finally, Table 6 lists the entries of the line candidates used in the construction of the LFs.

Appendix B Robustness of the CO Luminosity Functions

Here we test the robustness of the CO LFs constraints from ASPECS LP 3 mm by creating different realizations of the CO LFs after altering some of the assumptions discussed in the previous section, in particular concerning the fidelity of line candidates, and the CO excitation. The results of these tests are displayed in Figs. 10 and 11.

B.1 Impact of uncertain redshifts / sources with no counterparts

First, we compare our CO LFs and the constraints on the $\rho$(H2) evolution with redshift against the ones we infer, if we only subselect the galaxies for which a catalog redshift is available, and is consistent with the CO–based redshift within $|\delta z|<0.1$ (see Fig. 5). This automatically removes all the line candidates from the line search that lack a counterpart at other wavelengths, as well as potential misassociations with foreground/background galaxies.

The inferred CO luminosity functions are practically unaltered at their bright end. Small deviations are reported at the faint end, likely due to a combination of two reasons: 1) At the faint end, the impact of false positive candidates is larger. These spurious candidates by definition have counterparts only due to chance alignment, and it is unlikely that such counterparts have matching redshifts. 2) For reasonable ranges of the gas fraction $M_{\rm H2}/M_{*}$, fainter CO lines are typically associated with fainter stellar emission; these optical/NIR–faint galaxies might have relatively large redshift uncertainties, and might get scattered out of the $|\delta z|<0.1$ selection.

The direct consequence of these discrepancies is that $\rho$(H2) estimated only using sources with redshift–matching counterparts shows a faster decline at increasing redshifts at $z>3$, compared to our reference estimate, although the two estimates are well within 1-$\sigma$ uncertainties in both the CO LFs and $\rho$(H2) at any redshift. We thus conclude that our results, and in particular the steep evolution in $\rho$(H2) from present day to z\mathrel{\hbox to0.0pt{\lower 3.0pt\hbox{\sim}\hss}\raise 2.0pt\hbox{>}}1, are not significantly affected by our treatment of sources without clear counterparts or with ambiguous redshift associations.

B.2 Impact of CO excitation

We then examine the impact of the CO excitation assumptions on our estimates of the CO(1-0) LFs and on $\rho$(H2) (the CO LFs of the observed transitions are naturally unaffected by this assumption). We do so by repeating our analysis after assuming two extreme cases: a high excitation case corresponding to thermalized CO up to $J_{\rm up}$=4, and a low excitation scenario where the CO emission is modeled based on the Milky Way disk (see, e.g., Weiß et al., 2007; Carilli & Walter, 2013). A higher (lower) excitation implies fainter (brighter) $L^{\prime}_{\rm CO(1-0)}$ for a given line observed in a $J_{\rm up}>1$ transition, and therefore lower (higher) values of $M_{\rm H2}$. For reference, our fiducial assumption based on Daddi et al. (2015) lies roughly half the way between these two extreme cases for the transitions of interest here.

We find that a thermalized CO scenario would mitigate, but not completely solve, the friction between the ASPECS LP 3 mm CO LFs and the predictions by semi–analytical models. This is further explored in Popping et al. (2019). A low–excitation scenario, on the other hand, would exacerbate the tension. Evidence of a strong evolution in $\rho$(H2) between the local universe and $z>1$ is confirmed irrespective of the assumptions on the CO excitation, but for a low–excitation scenario, $\rho$(H2) appears nearly constant at any $z>1$, while it would drop rapidly at increasing redshifts, if a thermalized CO excitation is assumed.