<figure><img src="image_1.png" /></figure>

# The ALMA Spectroscopic Survey in the HUDF: The Cosmic Dust and Gas Mass Densities in Galaxies up to $z \sim 3$

Benjamin Magnelli¹, Leindert Boogaard², Roberto Decarli³, Jorge Gónzalez-López⁴,⁵, Mladen Novak⁶, Gergö Popping⁶,⁷, Ian Smail⁸, Fabian Walter⁶,⁹, Manuel Aravena¹⁰, Roberto J. Assef¹⁰, Franz Erik Bauer¹¹,¹²,¹³, Frank Bertoldi¹, Chris Carilli⁹,¹⁴, Paulo C. Cortes¹⁵,¹⁶, Elisabete da Cunha¹⁷, Emanuele Daddi¹⁸, Tanio Díaz-Santos¹⁰, Hanae Inami¹⁹, Robert J. Ivison⁷,²⁰, Olivier Le Fèvre²¹, Pascal Oesch²²,²³, Dominik Riechers⁶,²⁴,<fnref n="28" />, Hans-Walter Rix⁶, Mark T. Sargent²⁵, Paul van der Werf², Jeff Wagg²⁶, and Axel Weiss²⁷

¹Argelander Institut für Astronomie, Universität Bonn, Auf dem Hügel 71, Bonn, D-53121, Germany; magnelli@astro.uni-bonn.de

²Leiden Observatory, Leiden University, P.O. Box 9513, NL-2300 RA Leiden, The Netherlands

³INAF-Osservatorio di Astrofisica e Scienza dello Spazio, via Gobetti 93/3, I-40129, Bologna, Italy

⁴Las Campanas Observatory, Carnegie Institution of Washington, Casilla 601, La Serena, Chile

⁵Núcleo de Astronomía de la Facultad de Ingeniería y Ciencias, Universidad Diego Portales, Avenida Ejército Libertador 441, Santiago, Chile

⁶Max Planck Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany

⁷European Southern Observatory, Karl Schwarzschild Strasse 2, D-85748 Garching, Germany

⁸Centre for Extragalactic Astronomy, Department of Physics, Durham University, South Road, Durham, DH1 3LE, UK

⁹National Radio Astronomy Observatory, Pete V. Domenici Array Science Center, P.O. Box O, Socorro, NM 87801, USA

¹⁰Núcleo de Astronomía, Facultad de Ingeniería y Ciencias, Universidad Diego Portales, Av. Ejército 441, Santiago, Chile

¹¹Instituto de Astrofísica, Facultad de Física, Pontificia Universidad Católica de Chile Av. Vicuña Mackenna 4860, 782-0436 Macul, Santiago, Chile

¹²Millennium Institute of Astrophysics (MAS), Nuncio Monseñor Sótero Sanz 100, Providencia, Santiago, Chile

¹³Space Science Institute, 4750 Walnut Street, Suite 205, Boulder, CO 80301, USA

¹⁴Battcock Centre for Experimental Astrophysics, Cavendish Laboratory, Cambridge CB3 0HE, UK

¹⁵Joint ALMA Observatory—ESO, Av. Alonso de Córdova, 3104, Santiago, Chile

¹⁶National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA, 22903, USA

¹⁷International Centre for Radio Astronomy Research, The University of Western Australia, 35 Stirling Highway, Crawley, WA 6009, Australia

¹⁸Laboratoire AIM, CEA/DSM-CNRS-Université Paris Diderot, Irfu/Service d'Astrophysique, CEA Saclay, Orme des Merisiers, F-91191 Gif-sur-Yvette cedex, France

¹⁹Hiroshima Astrophysical Science Center, Hiroshima University, 1-3-1 Kagamiyama, Higashi-Hiroshima, Hiroshima, 739-8526, Japan

²⁰Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK

²¹Aix Marseille Université, CNRS, LAM (Laboratoire d'Astrophysique de Marseille) UMR 7326, F-13388, Marseille, France

²²Department of Astronomy, University of Geneva, 51 Chemin des Maillettes, 1290 Versoix, Switzerland

²³International Associate, Cosmic Dawn Center (DAWN) at the Niels Bohr Institute, University of Copenhagen and DTU-Space, Technical University of Denmark, Copenhagen, Denmark

²⁴Department of Astronomy, Cornell University, Space Sciences Building, Ithaca, NY 14853, USA

²⁵Astronomy Centre, Department of Physics and Astronomy, University of Sussex, Brighton, BN1 9QH, UK

²⁶SKA Organization, Lower Withington Macclesfield, Cheshire SK11 9DL, UK

²⁷Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany

Received 2019 September 27; revised 2019 December 18; accepted 2020 February 19; published 2020 March 30

## Abstract

Using the deepest 1.2 mm continuum map to date in the Hubble Ultra Deep Field, which was obtained as part of the ALMA Spectroscopic Survey (ASPECS) large program, we measure the cosmic density of dust and implied gas ($H_2+H I$) mass in galaxies as a function of look-back time. We do so by stacking the contribution from all H-band selected galaxies above a given stellar mass in distinct redshift bins, $\rho_{dust}(M_* > M, z)$ and $\rho_{gas}(M_* > M, z)$. At all redshifts, $\rho_{dust}(M_* > M, z)$ and $\rho_{gas}(M_* > M, z)$ grow rapidly as $M$ decreases down to $10^{10} M_\odot$, but this growth slows down toward lower stellar masses. This flattening implies that at our stellar mass-completeness limits ($10^8 M_\odot$ and $10^{8.9} M_\odot$ at $z \sim 0.4$ and $z \sim 3$), both quantities converge toward the total cosmic dust and gas mass densities in galaxies. The cosmic dust and gas mass densities increase at early cosmic time, peak around $z \sim 2$, and decrease by a factor $\sim 4$ and 7, when compared to the density of dust and molecular gas in the local universe, respectively. The contribution of quiescent galaxies (i.e., with little on-going star formation) to the cosmic dust and gas mass densities is minor ($\lesssim 10\%$). The redshift evolution of the cosmic gas density resembles that of the SFR density, as previously found by CO-based measurements. This confirms that galaxies have relatively constant star formation efficiencies (within a factor $\sim 2$) across cosmic time. Our results also imply that at $z \sim 0$, a large fraction ($\sim 90\%$) of dust formed in galaxies across cosmic time has either been destroyed or ejected to the intergalactic medium.

*Unified Astronomy Thesaurus concepts: High-redshift galaxies (734); Galaxy evolution (594); Galaxy formation (595)*

## 1. Introduction

The cosmic star formation rate density (SFRD) of the universe (i.e., mass of stars formed per unit time and comoving

volume; $\rho_{SFR}$) evolves significantly with redshift (see Madau & Dickinson 2014, for a review). It increased from early cosmic epochs, peaked at $z = 1-3$, and then decreased steadily until the present day. To understand this evolution, and therefore how galaxies formed and evolved throughout cosmic time, it is

<footnote n="28">²⁸ Humboldt Research Fellow.</footnote>
necessary to study their molecular gas reservoirs (i.e., the phase out of which stars form) and measure the evolution of the cosmic molecular gas mass density. There are many different approaches to measuring these gas reservoirs, the fundamental problem being that molecular hydrogen (H₂, the main constituent of the molecular gas) cannot be observed easily at the mass-weighted temperatures and density of the cold star-forming interstellar medium (ISM).

Traditionally, the emission lines of the different rotational states of the carbon monoxide molecule (¹²CO; hereafter, CO) have been used as tracers of the molecular gas (see Bolatto et al. 2013, for a review), but there are also other tracers. Most notably, the continuum emission from dust is frequently used as an alternative tracer of the gas, though including both the molecular (H₂) and atomic (H I) phases. With the advent of the Herschel Space Observatory, such dust-based gas mass estimates have been used for high-redshift galaxies (e.g., Magdis et al. 2012; Magnelli et al. 2012a; Santini et al. 2014; Genzel et al. 2015), by fitting their far-infrared-to-submillimeter emission with dust models (e.g., Draine & Li 2007) and using the local gas-to-dust mass ratio relation (e.g., Leroy et al. 2011). Most recently, Scoville et al. (2014, 2016, 2017) advocated that accurate dust-based gas mass estimates could be inferred using a single dust emission measurement in the Rayleigh-Jeans tail. This method relies on the assumption that the mass-dominant dust component of the ISM of most galaxies is at around 25 K, that this component accounts for the bulk of their Rayleigh-Jeans emission, and that the emission is optically thin (see Section 4.1). With this approach, Scoville et al. (2016) calibrated a single conversion factor ($\alpha_{850~\mu m}$) from the Rayleigh-Jeans dust emission to the (molecular) gas mass of massive galaxies (>10¹⁰ M☉). Interestingly, this conversion factor is consistent within a few percent of that inferred from Draine & Li (2007)'s dust mass absorption cross section and the typical gas-to-dust mass ratio of 100 for massive galaxies at $z \sim 0$ (Leroy et al. 2011). Dust-based gas mass estimates from far-infrared-to-submillimeter fits using Draine & Li (2007)'s model and from Scoville et al. (2016)'s Rayleigh-Jeans method are thus consistent for massive galaxies with a typical gas-to-dust mass ratio of 100 (e.g., Magnelli et al. 2019).

In recent years, using both methods, many studies have reached a common conclusion: at high-redshift ($z \gtrsim 0.2$) dust-based gas mass estimates are consistent within ~0.2 dex with those inferred from CO line emission (Genzel et al. 2015; Scoville et al. 2016, 2017; Tacconi et al. 2018; Kaasinen et al. 2019). This demonstrates the reliability of dust-based gas mass measurements for massive (>10¹⁰ M☉) galaxies and suggests that the gas in these galaxies is mostly dominated by its molecular phase. It has thereby facilitated the study of the gas content in high-redshift massive galaxies, which can be measured in just a few minutes of observing time with the Atacama Large Millimeter Array (ALMA).

From these dust-based gas mass estimates and the ever growing sample of CO measurements, much has been learned in recent years about the gas reservoirs of massive high-redshift galaxies (e.g., Magdis et al. 2012; Magnelli et al. 2012a; Bothwell et al. 2013; Saintonge et al. 2013; Santini et al. 2014; Genzel et al. 2015; Aravena et al. 2016; Decarli et al. 2016; Schinnerer et al. 2016; Scoville et al. 2016, 2017; Tacconi et al. 2018; Kaasinen et al. 2019). It is now robustly established that the gas fraction, $f_{gas} = M_{gas}/M_*$, in massive galaxies steadily decreases between $z \sim 4$ and $z \sim 0$, while their star formation

efficiency (i.e., SFR/M<sub>gas</sub>) only slightly decreases within this redshift range. Therefore, larger gas supply rather than enhanced star formation efficiency seems to explain the elevated specific star formation rate (SSFR; SFR/M<sub>*</sub>) of massive high-redshift galaxies compared to galaxies in the local universe (e.g., Schreiber et al. 2015).

Although these results are of utmost importance for galaxy evolution models, they cannot easily be extrapolated to infer the redshift evolution of the cosmic gas mass density in galaxies. Indeed, these targeted studies are biased toward massive, star-forming galaxies and could thus miss a significant fraction of gas-rich galaxies in the universe. Blind spectroscopic surveys at millimeter and radio wavelengths provide a complementary approach. The ALMA Spectroscopic Survey pilot and large program (ASPECS pilot and ASPECS LP, respectively; Decarli et al. 2016, 2019; Walter et al. 2016), and also the Jansky Very Large Array COLDz survey (Riechers et al. 2019) have been used to constrain the CO luminosity function and thereby the cosmic molecular gas mass density from $z \sim 4$ to $z \sim 0.3$. These studies have revealed that the cosmic molecular gas mass density closely matches the evolution of the cosmic SFRD, which implies that the average star formation efficiency in galaxies did not significantly evolve with redshift. Naturally, these studies also suffer from a number of limitations, particularly their dependencies on the assumed CO excitation and CO-to-H₂ conversion factors. To alleviate some of these uncertainties, independent constraints on the cosmic gas mass density using dust-based gas mass estimates are needed. These studies would simultaneously measure the redshift evolution of the cosmic dust mass density in galaxies, which to date remains only sparsely constrained (e.g., Dunne et al. 2003, 2011; Driver et al. 2018; Pozzi et al. 2019). This latter measurement would be instrumental for the growing number of galaxy evolution models that self-consistently track the production and destruction of dust (e.g., Popping et al. 2017; Aoyama et al. 2018; Davé et al. 2019; Vijayan et al. 2019).

As part of the ASPECS LP, we obtained the deepest ALMA 1.2 mm continuum map of the Hubble Ultra Deep Field (HUDF; $1\sigma$ sensitivity of 9.5 $\mu$Jy beam⁻¹; area of 4.2 arcmin²; Walter et al. 2016; González-López et al. 2020). The wavelength of this survey probes the Rayleigh-Jeans dust emission of galaxies up to $z \sim 4$, and is thus ideal to measure the dust and implied gas mass of high-redshift galaxies using Scoville et al. (2016)'s method. In addition, because this is a blind survey, it provides an unbiased view of the observed-frame 1.2 mm emission from all galaxies<fnref n="29" /> within the comoving volume probed by our map. Although a fraction of this 1.2 mm emission is included in individual detections, a large portion could reside below our detection threshold, even in the case of this deep 1.2 mm map. Fortunately, the HUDF is one of the best studied extragalactic regions in the sky. It thus benefits from a remarkable wealth of ancillary data, providing a unique opportunity to build a stellar mass-complete sample of galaxies down to, e.g., ~10⁸ M☉ and ~10⁸.⁹ M☉ at $z \sim 0.4$ and 3, respectively (e.g., Mortlock et al. 2015). Knowing a priori the positions of this stellar mass-complete sample of galaxies, we can thus sum up their 1.2 mm emission within a

<footnote n="29">²⁹ The only exception being galaxies with very extended emission ($\gtrsim 10''$), which could be missed by our observations due to the lack of very short baselines. However, such extended emission would only be associated to low-redshift galaxies ($z \lesssim 0.1$) which are not studied in the present paper.</footnote>
given comoving volume, irrespective of their individual detectability in our 1.2 mm map. Converting this 1.2 mm emission per unit comoving volume into dust and gas masses, assuming $T_{\text{dust}} = 25$ K (Scoville et al. 2016) and the local gas-to-dust mass ratio relation (Leroy et al. 2011), we can measure the cosmic dust and gas mass densities of all the known galaxies in the HUDF above a given stellar mass and as a function of look-back time; i.e., $\rho_{\text{dust}}(M_{*} > M, z)$ and $\rho_{\text{gas}}(M_{*} > M, z)$. These measurements provide constraints for galaxy evolution models and complement those inferred from the ASPECS-CO survey (Decarli et al. 2016, 2019).

The structure of this paper is as follows. In Section 2, we present the ASPECS LP 1.2 mm continuum map. In Section 3, we summarize the ancillary data used in this study and the build-up of our stellar mass-complete sample of galaxies. In Section 4, we present the method that we used to measure the cosmic dust and gas mass densities through stacking of the ASPECS LP 1.2 mm map. In Sections 5 and 6, we present the main results of this study; i.e., the evolution of cosmic dust and gas mass densities as a function of stellar mass and look-back time. In Section 7, we compare these results with outputs from simulations and discuss some of the implications for galaxy evolution models. Finally, we present our conclusions in Section 8.

Throughout this paper, we assume a $\Lambda$CDM cosmology, adopting $H_0 = 67.8$ (km s$^{-1}$) Mpc$^{-1}$, $\Omega_M = 0.308$ and $\Omega_\Lambda = 0.692$ (Planck Collaboration et al. 2016). At $z = 1$, 1" corresponds to 8.229 kpc. A Chabrier (2003) initial mass function (IMF) is used for all stellar masses quoted in this article.

## 2. Data

The ASPECS LP 1.2 mm survey covers 4.2 arcmin$^2$ in the HUDF, centered at $\alpha = 3^h\ 32^m\ 38^s5$, $\delta = -27^\circ\ 47'\ 00''$ (J2000; 2016.1.00324.L). The survey strategy, and also the data calibration and imaging are described in detail by González-López et al. (2020). Here, we only summarize the most important information.

The ASPECS LP 1.2 mm continuum map was obtained by combining eight spectral tunings that cover most of the ALMA band 6. The mosaic consists of 85 pointings and is Nyquist-sampled at all wavelengths. The data was calibrated with the Common Astronomy Software Applications (CASA; McMullin et al. 2007) calibration pipeline using the script provided by the Joint ALMA Observatory (JAO). Imaging was also done in CASA using the multi-frequency synthesis algorithm implemented within the task TCLEAN, which combines all pointings together down to a primary beam (PB) gain of 0.1. We used natural weighting and “cleaned” down to $20~\mu\text{Jy beam}^{-1}$ all sources with a signal-to-noise ratio (S/N) greater than 5. The synthesized beam has a full width at half maximum (FWHM) of $1.53 \times 1.08$ at a position angle of $-84^\circ$. The mosaic covers 2.9 and 4.2 arcmin$^2$ of the HUDF, within a combined PB coverage$^{30}$ of 50% and 10%, respectively. The deepest region in the map (i.e., with a combined PB coverage $\ge 90\%$) has an rms of $9.5~\mu\text{Jy beam}^{-1}$. Where the combined PB coverage is better than 75% (i.e., the region of interest of our study; see Section 4.2), we detected 22 galaxies with a S/N greater than 3 and a “Fidelity” factor greater than 0.5.

## 3. Sample

The ASPECS LP 1.2 mm survey covers most of the *Hubble* eXtremely Deep Field (XDF; Koekemoer et al. 2013; Illingworth et al. 2013), which is itself located within the HUDF (Beckwith et al. 2006). This is one of the best studied extragalactic regions in the sky, and thereby benefits from a remarkable wealth of ancillary data. The compilation of our master catalog of galaxies and the modeling of their spectral energy distribution (SED) are described by Decarli et al. (2019) and Boogaard et al. (2019), respectively. Here, we summarize the most important information.

In the XDF, the bulk of the optical and near-infrared observations comes from the *Hubble Space Telescope (HST)* as part of the Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey (CANDELS; Grogin et al. 2011; Koekemoer et al. 2011) and the HUDF09 (e.g., Bouwens et al. 2011) and HUDF12 (Koekemoer et al. 2013) surveys. These observations were obtained with the Advanced Camera for Surveys (ACS) at optical wavelengths and with the Wide Field Camera 3 (WFC3) in the near-infrared (NIR). By combining all of the *HST* observations, Skelton et al. (2014) performed a multi-wavelength photometric analysis, which also included publicly available ground-based optical/NIR (see Skelton et al. 2014, and reference therein) and also *Spitzer*-IRAC images (Labbé et al. 2015). Complemented with *Spitzer* MIPS-24 µm photometry from Whitaker et al. (2014), this constitutes our master photometric catalog. This provides measurements in $>30$ broad and medium bands for 1481 sources in the region of interest in the XDF; i.e., where the combined PB coverage of the ASPECS LP 1.2 mm survey is better than 75% (see Section 4.2).

The spectroscopic redshifts of 443 of these galaxies were obtained from a variety of studies: the MUSE *Hubble* Ultra Deep Survey (Bacon et al. 2017; Inami et al. 2017); *HST* grism spectroscopy in the optical (Xu et al. 2007) and in the NIR (3D-*HST* survey; Momcheva et al. 2016); and spectroscopic compilations from Le Fèvre et al. (2005), Coe et al. (2006), Skelton et al. (2014), and Morris et al. (2015). For galaxies with no available spectroscopic redshift, we use photometric redshifts determined in Skelton et al. (2014) by means of the *EAZY* code, with a typical uncertainty, $\sigma[|z_{\text{phot}} - z_{\text{spec}}|(1+z_{\text{spec}})]$, of 0.010, and only 5.4% of objects with $|z_{\text{phot}} - z_{\text{spec}}|(1+z_{\text{spec}}) > 0.1$.

The stellar mass of each galaxy was obtained by modeling their SED, using the high-redshift extension of the *MAGPHYS* code (da Cunha et al. 2008, 2015). We used the galaxy’s photometry between 0.37 and 8.0 µm, as well as their 1.2 mm flux density (or upper limit). These stellar masses are on average ~0.2 dex larger than those measured by Skelton et al. (2014) using the *FAST* code. We verified that our results on the cosmic dust and gas mass densities remain unchanged—they simply shifted toward lower stellar masses, while using the stellar masses of Skelton et al. (2014).

From this master catalog, we kept only the 555 galaxies with an observed total$^{31}$ *H*-band magnitude brighter than 27, of which 281 had a spectroscopic redshift. This NIR selection was chosen because the observed-frame *H*-band luminosity of a galaxy is known to correlate with its stellar mass up to $z \sim 3$ (Figure 1). Furthermore, in our master catalog, the *H*-band magnitude distribution (i.e., number of galaxies per bin of

<footnote>$^{30}$ This corresponds to the “.pb” array output by the task TCLEAN when used in “mosaic” mode.</footnote>

<footnote>$^{31}$ We refer the reader to Skelton et al. (2014) for details on how these “total” *H*-band magnitudes were measured.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Correlation in different redshift bins between the stellar mass and observed-frame *H*-band luminosity of our *H*-band-selected galaxies (i.e., *H* ≲ 27; filled circles; 555 galaxies; green filled circles represent galaxies individually detected at 1.2 mm—21 galaxies, while the red filled circles are galaxies undetected at 1.2 mm—534 galaxies). Open circles present galaxies within the XDF but with *H* > 27, blue squares show galaxies detected by IRAC at 5.8 µm. For galaxies detected by IRAC at 5.8 µm but not detected in the *H*-band (2 galaxies), we artificially set their observed-frame *H*-band luminosity to 10⁷ *L*☉. The shaded areas define the range of our *H* = 27 selection limit at the lowest and highest redshift of each bin. The black lines are the best-fit relations, while the dashed lines show their 1σ dispersion. The thick-black vertical lines correspond to the stellar mass-completeness limits of our *H*-band-selected sample, defined at the intersection between our *H*-band selection limits at the highest redshift of the bin (upper horizontal boundary of the gray-shaded area) and the lower 1σ boundary of the best-fit relation. The thin dotted lines show the one-to-one relation. At *z* ≳ 3.2, the observed-frame *H*-band luminosities of galaxies do not strongly correlate with their stellar masses, as expected (see text). Thus, at these redshifts, our *H*-band-selected sample cannot be considered as stellar mass-selected and therefore stellar mass-complete.</figcaption></figure>

magnitude) peaks at *H* = 27 and rapidly decreases toward fainter magnitudes. This rapid decrease at *H* > 27 suggests that at such faint magnitudes, our catalog is affected by large photometric incompleteness. We thus restricted our analysis to sources with *H* ≲ 27. This *H* ≲ 27 selection also ensured that the number of broad and medium bands available for each source was high enough (15<sup>+8</sup><sub>-7</sub>, median and 16th and 84th percentiles) for accurate SED modeling and thus stellar mass estimates.

In the last step, we determined down to which stellar masses our *H*-band-selected galaxy sample can be considered stellar mass-complete. To this end, we used the empirical approach described in Schreiber et al. (2015). For each redshift bin of interest, we fitted the correlation between the observed-frame *H*-band (1.6 µm) luminosities of galaxies and their stellar masses (derived as described above) with a simple power law; i.e., $M_* = CL^\alpha$ (Figure 1). Then, we estimated the scatter around this correlation, which is caused by differences of age, attenuation and *k*-correction between these galaxies. Finally, for a given redshift bin, the stellar mass-completeness limit was set by the stellar mass corresponding to the *H*-band luminosity cut plus the 1σ dispersion of the $L_{1.6\mu m/(1+z)} - M_*$ relation. At this stellar mass-completeness limit, only 16% of galaxies are expected to be missed because of our *H* ≲ 27 selection

criterion, while this percentage drops rapidly to 0% toward higher stellar masses (see open symbols in Figure 1). At *z* ≳ 3.2, because the *H*-band probes rest wavelengths that are shorter than 4000 Å (i.e., Balmer break), the observed-frame *H*-band luminosity of a galaxy no longer correlates strongly with its stellar mass (Figure 1). At these redshifts, our *H*-band-selected sample cannot be considered as stellar mass-selected and therefore stellar mass-complete.

We note that while the observed-frame IRAC luminosities of galaxies correlate in principle better with their stellar masses than the observed-frame *H*-band luminosities, IRAC observations in the XDF do not provide stellar mass-complete samples as deep as those provided by the *H*-band. To illustrate this point, we plot in Figure 1 the locus of all galaxies detected by IRAC at 5.8 µm within our region of interest in the XDF, including those not detected in the *H*-band (2 galaxies) and to which we artificially attributed an observed-frame *H*-band luminosity of 10⁷ *L*☉. IRAC 5.8 µm selected galaxies are associated with the brightest and most massive galaxies in our *H*-band-selected catalog. At *z* < 3, only one IRAC 5.8 µm selected galaxies is missed by our *H*-band selection, which has a stellar mass below our stellar mass completeness limit. This supports the assumption that the observed-frame *H*-band luminosity of a galaxy is a good proxy of its stellar mass up
<table><caption>**Table 1**

Cosmic Dust Mass Density in Galaxies in Different Redshift Bins and Above a Given Stellar Mass; i.e., $\rho_{\text{dust}}(M_{*} > M, z)$</caption><thead><tr><th rowspan="2">Redshift</th><th colspan="7">ρ<sub>dust</sub>(M<sub>*</sub> > M, z)</th></tr><tr><th colspan="7">[×10<sup>5</sup> M<sub>☉</sub> Mpc<sup>-3</sup>]</th></tr></thead><tbody><tr><td>Bin</td><td>&gt;10<sup>8</sup> M<sub>☉</sub></td><td>&gt;10<sup>8.5</sup> M<sub>☉</sub></td><td>&gt;10<sup>9</sup> M<sub>☉</sub></td><td>&gt;10<sup>9.5</sup> M<sub>☉</sub></td><td>&gt;10<sup>10</sup> M<sub>☉</sub></td><td>&gt;10<sup>10.5</sup> M<sub>☉</sub></td><td>&gt;10<sup>11</sup> M<sub>☉</sub></td></tr><tr><td>0.3 ≤ z &lt; 0.6</td><td>Nb = 8<br/>1.1 ± 0.8 (0.5)<br/>[1.2<sup>+0.6</sup><sub>-0.5</sub>]</td><td>Nb = 3<br/>1.0 ± 0.7 (0.3)</td><td>Nb = 2<br/>0.8 ± 0.6 (0.2)<br/>[0.5<sup>+0.3</sup><sub>-0.2</sub>]</td><td>Nb = 1<br/>0.3 ± 0.3 (0.2)</td><td>...</td><td>...</td><td>...</td></tr><tr><td>0.6 ≤ z &lt; 1.0</td><td>...</td><td>Nb = 50<br/>3.0 ± 0.9 (0.7)</td><td>Nb = 30<br/>2.7 ± 0.8 (0.6)<br/>[2.6<sup>+0.7</sup><sub>-0.6</sub>]</td><td>Nb = 18<br/>1.9 ± 0.7 (0.4)</td><td>Nb = 12<br/>1.4 ± 0.6 (0.3)<br/>[1.4<sup>+0.4</sup><sub>-0.6</sub>]</td><td>Nb = 8<br/>0.6 ± 0.3 (0.3)</td><td>Nb = 2<br/>0.2 ± 0.2 (0.1)</td></tr><tr><td>1.0 ≤ z &lt; 1.6</td><td>...</td><td>Nb = 100<br/>4.7 ± 1.4 (0.5)<br/>[4.3<sup>+1.8</sup><sub>-1.4</sub>]</td><td>Nb = 48<br/>4.3 ± 1.3 (0.3)<br/>[4.2<sup>+1.1</sup><sub>-1.3</sub>]</td><td>Nb = 27<br/>4.2 ± 1.3 (0.3)</td><td>Nb = 13<br/>3.1 ± 1.3 (0.2)<br/>[3.6<sup>+1.1</sup><sub>-1.3</sub>]</td><td>Nb = 6<br/>2.4 ± 1.2 (0.1)</td><td>Nb = 2<br/>1.6 ± 1.2 (0.07)</td></tr><tr><td>1.6 ≤ z &lt; 2.3</td><td>...</td><td>...</td><td>Nb = 27<br/>3.5 ± 1.1 (0.2)<br/>[3.6<sup>+0.9</sup><sub>-1.0</sub>]</td><td>Nb = 15<br/>3.3 ± 1.1 (0.1)</td><td>Nb = 10<br/>3.0 ± 1.1 (0.1)<br/>[2.6<sup>+0.9</sup><sub>-0.1</sub>]</td><td>Nb = 4<br/>2.1 ± 1.0 (0.07)</td><td>...</td></tr><tr><td>2.3 ≤ z &lt; 3.2</td><td>...</td><td>...</td><td>Nb = 27<br/>3.1 ± 1.1 (0.2)<br/>[3.2<sup>+0.9</sup><sub>-0.7</sub>]</td><td>Nb = 15<br/>3.0 ± 1.1 (0.1)</td><td>Nb = 8<br/>2.9 ± 1.1 (0.1)<br/>[2.1<sup>+0.6</sup><sub>-0.4</sub>]</td><td>Nb = 2<br/>1.1 ± 0.8 (0.04)</td><td>...</td></tr><tr><td>3.2 ≤ z &lt; 4.5</td><td colspan="7">Nb = 44 0.6 ± 0.2 (0.1)</td></tr><tr><td>4.5 ≤ z &lt; 5.5</td><td colspan="7">Nb = 9 0.05 ± 0.07 (0.07)</td></tr></tbody></table>

Note. The total uncertainties correspond to the quadratic sum of the measurement uncertainties ($\sigma_{\rho}^{m}$) and the Poissonian uncertainties ($\sigma_{\rho}^{\text{Poisson}}$). The measurement uncertainties are provided in parentheses and should be used to assess the detection significance in the stacked stamps. Nb gives the number of stacked galaxies. In squared brackets, we provide the cosmic dust mass densities inferred by fitting the variation of $\rho_{\text{dust}}(M_{*} > M, z)$ using the SMF of Mortlock et al. (2015) and solving for $f_{\text{dust}}$ assuming $f_{\text{dust}} = M_{\text{dust}}/M_{*} = A \times (M / 10^{10.7} M_{\odot})^{B}$ (see Section 5.2). At $3.2 \le z < 4.5$ and $4.5 \le z < 5.5$, where the sample cannot be considered as stellar mass-selected, we summed the contribution of all galaxies with $H \le 27$, irrespective of their stellar masses. Differential measurements can be inferred by subtracting accordingly the values of interest.

to $z < 3$, and that in the XDF, a H-band-selected catalog has a much lower stellar mass-completeness limits than an IRAC 5.8 µm selected catalog. Repeating this analysis with IRAC 3.6 µm leads to the same conclusions, though with stellar mass-completeness limits getting closer to that of our H-band-selected catalog.

To verify that up to $z \sim 3$ our H-band-selected catalog was indeed “complete” down to our stellar mass-completeness limits, we measured the stellar mass functions by counting the number of galaxies in bins of redshifts and stellar masses, normalized by the volume probed by the XDF. Down to our stellar mass-completeness limits, these stellar mass functions agree, within the uncertainties, with the fits inferred by Mortlock et al. (2015) and Davidzon et al. (2017) in the CANDELS and COSMOS fields, respectively.

Finally, we verified that our catalog did not miss any obvious dust emitters; i.e., galaxies already detected by the ASPECS LP 1.2 mm survey but not in our H-band selected sample. There are 22 galaxies detected by González-López et al. (2020) in the ASPECS LP 1.2 mm continuum map where the combined PB coverage is better than 75% (i.e., the region of interest of our study; see Section 4.2). Among these sources, 21 have a counterpart in our H-band selected catalog within the synthesized beam of the ASPECS observations (i.e., a radius of $0''/6$). The remaining source is one of the faintest sources detected by González-López et al. (2020), who also reported no clear NIR counterpart. Assuming that this source is real and at a redshift of $z \sim 0.45$, $z \sim 0.80$, $z \sim 1.30$, $z \sim 1.95$, $z \sim 2.75$, $z \sim 3.9$ or $z \sim 5$, then it would increase the cosmic dust mass

densities inferred here by ~50%, 10%, 3%, 3%, 3%, 10%, 15%, respectively; i.e., well within our total uncertainties (see Table 1). For the cosmic gas mass densities, the impact of this source would depend on its stellar mass: assuming $M_{*} = 10^{9.5} M_{\odot}$, this implies upwards corrections by ~40%, 10%, 6%, 6%, 6%, 10%, 40%, respectively; while assuming $M_{*} = 10^{10.5} M_{\odot}$, this implies upwards corrections by ~30%, 7%, 4%, 4%, 4%, 5%, 16%, respectively. However, because this source is very faint and on the lower end of the S/N (=4.1) and “Fidelity” (=0.78) selection criteria of González-López et al. (2020), it could well be a spurious source (i.e., positive noise peak in the ASPECS LP 1.2 mm survey).

## 4. Method

By combining our stellar mass-complete galaxy sample with our deep ASPECS LP 1.2 mm survey, we can measure in several redshift bins the total dust and gas mass contained in these galaxies and from that infer their cosmic dust and gas mass densities, knowing the comoving volume probed by our survey. In this section, we first summarize the method used to convert observed-frame 1.2 mm flux densities into dust and gas masses. We follow by describing the method used to infer the cosmic dust and gas mass densities from these galaxies through stacking.

### 4.1. Measuring $M_{\text{dust}}$ and $M_{\text{gas}}$ from $S_{1.2\text{mm}}$

In the optically thin approximation, which is almost always valid at the long wavelengths probed by our observations, the dust mass ($M_{\text{dust}}$ in $M_{\odot}$) of a galaxy can be inferred using its
(sub)millimeter flux density ($S_{ν_0}$ in Jy) at the observed-frame frequency, $ν_{obs} = ν_{rest}/(1+z)$, following, for example, Kovács et al. (2010),

$$M_{\text{dust}} = \frac{5.03 \times 10^{-31} \cdot S_{\nu_{\text{obs}}} \cdot D_{\text{L}}^2}{(1+z)^4 \cdot B_{\nu_{\text{obs}}}(T_{\text{obs}}) \cdot \kappa_{\nu_0}} \cdot \left(\frac{\nu_0}{\nu_{\text{rest}}}\right)^{\beta}, \quad (1)$$

where $B_{ν_0}(T_{obs})$ is Planck's blackbody function in Jy sr$^{-1}$ at the observed-frame temperature $T_{obs}$ in Kelvin, which relates to the rest-frame temperature $T$ as $T_{obs} = T/(1+z)$, $D_L$ is the luminosity distance in meter, $β$ is the dust emissivity spectral index and $κ_{ν_0}$ is the photon cross-section to mass ratio of dust (in m$^2$ kg$^{-1}$) at rest-frequency $ν_0$.

As advocated by Scoville et al. (2016), we used a mass-weighted mean dust temperature of $\langle T \rangle_M = T = 25$ K. This is the mass-dominant dust component of the ISM of galaxies and accounts for the bulk of their Rayleigh-Jeans dust emission (Scoville et al. 2016). Part of the dust in the ISM can (and will) be at higher temperatures but only in localized regions with a negligible contribution to the global dust mass and Rayleigh-Jeans emission. A value of 25 K is supported by Herschel-based studies of local and high-redshift galaxies, which find a range of $T \sim 15-35$ K (Dunne et al. 2011; Dale et al. 2012; Auld et al. 2013; Magnelli et al. 2014). It is further supported at high redshift by the recent work of Kaasinen et al. (2019), which compared CO-based and dust-based gas measurements at $z \sim 2$. Note that in the Rayleigh-Jeans tail probed by our observations, dust mass estimates vary as $T^{-1}$. Thus, a dust temperature range of $T \sim 15-35$ K implies a systematic uncertainty for our dust masses of 25%-50%.

As suggested by Scoville et al. (2016), we used $β = 1.8$, which corresponds to the Galactic measurement made by the Planck Collaboration et al. (2011) and is well within the range observed in high-redshift star-forming galaxies (e.g., Chapin et al. 2009; Magnelli et al. 2012a). Varying $β$ within the range predicted by most theoretical models—i.e., $β = 1.5-2.0$ (Draine 2011)—does not impact our results significantly. Indeed, our dust mass estimates would simply be multiplied by 0.99 (1.00), 1.06 (0.96), 1.14 (0.91), 1.23 (0.87) and 1.32 (0.83), at $z \sim 0.45$, $z \sim 0.80$, $z \sim 1.30$, $z \sim 1.95$, and $z \sim 2.75$, while using $β = 1.5$ (2.0) instead of $β = 1.8$, respectively. Finally, we adopted $κ_{ν_0} = 0.0431$ m$^2$ kg$^{-1}$ with $ν_0 = 352.6$ GHz (i.e., 850 µm; Li & Draine 2001). Interestingly, assuming a typical gas-to-dust mass ratio of 100, this dust mass absorption cross section is within a few percent of the “ISM” mass absorption cross section calibrated by Scoville et al. (2016; i.e., their $α_{850}$ µm).

Finally, following da Cunha et al. (2013), we corrected our dust mass measurements for the effect of the cosmic microwave background (CMB), which temperature increases as $T_{CMB}(z) = 2.73 \times (1+z)$. First, the CMB acts as an additional source of heating of the mass-dominant dust component of galaxies, increasing its temperature from 25 K at $z = 0$ to 25.3 K at $z = 5$ (Equation (12) of da Cunha et al. 2013). Second, the CMB acts as a background against which we make our measurements, implying an underestimation of the intrinsic flux densities of galaxies (Equation (18) of da Cunha et al. 2013). Although taken into account, these effects have a relatively minor impact on our results because they yield upward corrections of our measurements by only 1.01, 1.02, 1.03, 1.04, and 1.07, at $z \sim 0.45$, $z \sim 0.80$, $z \sim 1.30$, $z \sim 1.95$, and $z \sim 2.75$, respectively. At $z = 3.9$ and

$z = 5.0$, where our sample cannot be considered as stellar mass complete, these CMB corrections have a value of 1.12 and 1.20, respectively.

Dust masses can be converted into gas masses, assuming a gas-to-dust mass ratio, which can be a function of metallicity (e.g., Leroy et al. 2011; Eales et al. 2012; Magdis et al. 2012; Magnelli et al. 2012b; Santini et al. 2014; Genzel et al. 2015; Scoville et al. 2016; Tacconi et al. 2018). As in Tacconi et al. (2018), we used a gas-to-dust mass ratio ($δ_{GDR}$) that correlates nearly linearly with metallicity and is equal to 100 at solar metallicity, $[12 + \log(O/H)]_⊙ = 8.67$; i.e.,

$$\delta_{\text{GDR}} = \frac{M_{\text{gas}}}{M_{\text{dust}}} = 10^{(+2-0.85 \cdot (12+\log(O/H)_{\text{PP04}}-8.67))}, \quad (2)$$

where $12 + \log(O/H)_{PP04}$ is the gas phase metallicity adopting the Pettini & Pagel (2004) scale. As already mentioned, at solar metallicity (i.e., $δ_{GDR} = 100$), the combination of Equations (1) and (2) is equivalent to the method advocated by Scoville et al. (2016) for massive galaxies. In Equation (2), $M_{gas}$ includes both the molecular ($H_2$) and atomic (H I) phases and a standard 36% mass fraction correction to account for helium. Equation (2) also accounts for the molecular phase that is fully molecular in $H_2$ but dissociated in CO (the so-called CO-dark phase; see Leroy et al. 2011).

The metallicity of our galaxies is inferred using the stellar mass-metallicity relation following Tacconi et al. (2018),

$$\begin{align}
12 + \log(O/H)_{\text{PP04}} &= a - 0.087 \times (\log(M_*) - b)^2, \quad \text{with} \nonumber \\
&\quad a = 8.74, \text{ and} \nonumber \\
&\quad b = 10.4 + 4.46 \times \log(1+z) - 1.78 \times (\log(1+z))^2. \tag{3}
\end{align}$$

We note that while the $δ_{GDR}$-metallicity relation is believed to be nearly linear at $z \sim 0$ down to $12 + \log(O/H) \sim 7.9$ (Leroy et al. 2011; Rémy-Ruyer et al. 2014), at lower metallicities observations suggest that this relation might follow a steeper power-law (~3; Rémy-Ruyer et al. 2014, see also Coogan et al. 2019), yielding larger gas masses per unit dust mass. In our sample, galaxies with stellar masses close to our stellar mass-completeness limits have metallicities in the range 7.7-7.9. Consequently, their gas mass could be underestimated by a factor of a few when using Equation (2). However, the metallicity at which this transition takes place is very uncertain ($7.94 \pm 0.47$; Rémy-Ruyer et al. 2014) and could be irrelevant for our galaxies. The effect of such steep $δ_{GDR}$-metallicity relation at very low metallicities on the inferred cosmic gas mass densities is further discussed in Section 6.1.

## 4.2. Stacking Analysis

We can measure the total observed-frame 1.2 mm flux density of a galaxy population knowing their positions within our survey. To this end, we can start by creating stamps of the ASPECS LP 1.2 mm survey centered around each of these galaxies. Then, for a given redshift bin and stellar mass range of interest, we can stack these stamps together, obtaining thereby at the center of this stacked stamp the total emission of a given galaxy population at 1.2 mm. Finally, using the comoving volume probed by our survey in this redshift bin, this total emission can be converted into comoving dust and gas
mass densities by applying the equations given in Section 4.1. However, to obtain robust measurements out of this relatively straightforward methodology, precautions must be taken, which are summarized below.

As a first step, to account for the different astrometric solution between the ALMA data and optical/NIR data in the XDF, we applied an astrometry offset (ΔR.A. = +0.076", Δdecl. = -0.279") to all of the positions in our master catalog (Rujopakarn et al. 2016; Dunlop et al. 2017; Cibinel et al. 2017; Franco et al. 2018; Decarli et al. 2019).

We then restricted our analysis to the region of the XDF where the combined PB coverage of the ASPECS LP 1.2 mm survey is better than 75%. This ensures that our stacking analysis does not include the relatively noisy edges of the survey. The surface area of the ASPECS LP 1.2 mm survey with a combined PB coverage ≥75% is 2.27 arcmin². This is the surface area that we used to compute the comoving volume probed by the survey in a given redshift bin.

For galaxies that are individually detected (see Section 2 and González-López et al. 2020) by the ASPECS LP 1.2 mm survey, stamps were cut out from the “clean” band 6 continuum mosaic. For galaxies undetected by the survey, stamps were cut out from the “residual” mosaic; i.e., the “clean” mosaic where detected sources were removed using the “clean” synthesized beam (a 2D Gaussian approximation of the synthesized beam created by the task TCLEAN). This ensures that we do not count the flux density of detected galaxies several times in the stacking analysis (for stacked sources within 1"-2" from a detected source) and that the background of our stacked stamps is not dominated by individually detected sources (for stack sources within 2"-10" from a detected source). However, we verified that consistent results are found when stamps for undetected sources were instead cut from the “clean” mosaic. The size of each stamp is 20" × 20", which allows for an ample number of independent beams in the stacked stamps, crucial for accurate noise estimates. Note that even at the depth of the ASPECS LP 1.2 mm survey, confusion noise is negligible (González-López et al. 2020).

The conversion of the observed-frame 1.2 mm flux density of a galaxy into its dust and gas masses (see Equations (1)-(3)) depends on its redshift and metallicity (and thus stellar mass). Therefore, before proceeding with the stacking, we converted each galaxy stamp from Jy beam⁻¹ to comoving M_dust or M_gas density units. To this end, we simply applied the conversion factors to each pixel, normalized by the comoving volume probed by the survey in the redshift bin of interest. The final stacked stamps are thus directly in units of M_☉ Mpc⁻³ beam⁻¹.

The signal in the stacked stamps was measured as follows. First, we measured at their center the signal within an aperture with a 0.6 × FWHM beam radius, optimized for point-source detection (i.e., equivalent to a matched filter). This signal (S_aper) was then compared to the noise (N_aper), defined as the standard deviation of the signal distribution within 200 apertures randomly positioned in the stacked stamp. If S_aper/N_aper ≥ 5, then we fitted the stacked stamp with a 2D Gaussian function. The amplitude (A_fit), as well as the minor (θ_min) and major (θ_maj) FWHM of this 2D Gaussian function centered on the stacked stamp were left as free parameters, accommodating thereby any possible astrometric mismatch between the near-infrared and millimeter centers of these galaxies. The comoving mass density (i.e., ρ_dust or ρ_gas, depending on which conversion was applied beforehand to the 1.2 mm stacked stamps) and its

associated measurement uncertainty (σ_ρ) were then derived using the standard formulae,

$$ρ_{\text{dust}} \text{ or } ρ_{\text{gas}} [M_{\odot} \text{ Mpc}^{-3}] = A_{\text{fit}} \cdot \frac{\theta_{\text{min}}^{\text{fit}} \cdot \theta_{\text{maj}}^{\text{fit}}}{\theta_{\text{min}}^{\text{beam}} \cdot \theta_{\text{maj}}^{\text{beam}}}, \quad (4)$$

$$\sigma_{\rho}^{m} [M_{\odot} \text{ Mpc}^{-3}] = \sigma_{\text{pixel}} \cdot \sqrt{\frac{\theta_{\text{min}}^{\text{fit}} \cdot \theta_{\text{maj}}^{\text{fit}}}{\theta_{\text{min}}^{\text{beam}} \cdot \theta_{\text{maj}}^{\text{beam}}}}, \quad (5)$$

where θ_maj and θ_min are the minor and major FWHM of the synthesized beam; and σ pixel is the standard deviation of the stacked stamp pixel distribution.

In stacked stamps where S_aper/N_aper < 5, the comoving mass density (ρ_dust or ρ_gas) was instead given by the total signal in an aperture with a 1.0 × FWHM beam radius, divided by the area of the synthesized beam in units of pixel; and σ_ρ was given by the standard deviation of the signal distribution within 200 apertures randomly positioned in the stacked stamp, divided by the area of the synthesized beam in units of pixel. This particular radius was chosen to be large enough to encompass the total signal in the stacked stamp and to provide consistent signal measurements with respect to our 2D Gaussian fits when applied on clear detections (i.e., when S_aper/N_aper > 5).

The total uncertainty (σtot) was defined as the quadratic sum of the measurement uncertainty (σ_ρ) and the Poissonian uncertainty (σ_Poisson). The former accounts for the noise in the stacked stamp (i.e., detection significance), while the latter accounts for the uncertainty due to the low number of galaxies stacked in each bin (i.e., ∝ √1/Nb).

We performed the stacking analysis in five redshift bins in which our galaxies can be considered as stellar mass-selected; i.e., 0.3 ≤ z < 0.6, 0.6 ≤ z < 1.0, 1.0 ≤ z < 1.6, 1.6 ≤ z < 2.3, and 2.3 ≤ z < 3.2 (Section 3 and Figure 1). These redshift bins were chosen to sample the cosmic history with a reasonably large number of sources (>8) per < 2 Gyr look-back time intervals. In these redshift bins, we stacked all galaxies with a stellar mass greater than a given threshold, starting at 10¹¹ M_☉ and decreasing it in steps of 0.5 dex down to our stellar mass-completeness limit.

The stacked stamps obtained for the dust are shown in Figure 2, while the measured cosmic dust and gas mass densities—i.e., ρ_dust(M_☉ > M, z) and ρ_gas(M_☉ > M, z)—and associated uncertainties are given in Tables 1 and 2.

From right to left in Figure 2, the stack stamps include ever less massive galaxies. Because these are cumulative stacked stamps, individual panels are not independent. At z > 0.6, the stacking analysis yields significant detections (≥5σ_ρ) for most stellar mass thresholds down to our stellar mass-completeness limits. At 0.3 ≤ z < 0.6, the stacking analysis mostly yields tentative detections (<3σ_ρ) but these measurements provide meaningful constraints on the cosmic dust mass density at this redshift and above these stellar mass thresholds (see Section 5.1). Note that the counter-intuitive increase in detection significance from 0.3 ≤ z < 0.6 (top row in Figure 2) to 2.3 ≤ z < 3.2 (bottom row) is due to (i) the intrinsically lower dust mass content in low-redshift galaxies (see Section 5.1); (ii) the smaller cosmic volume and thus the fewer number of galaxies in the lower redshift bins (see Figure 1); and (iii) to the fact that on the Rayleigh-Jeans tail and at a given observed wavelength, low-redshift galaxies are not significantly brighter per unit dust mass
<figure><img src="image_3.png" /><figcaption>**Figure 2.** 16" × 16" zoom-in on the cumulative stacked stamps corresponding to the comoving dust mass density in galaxies in different redshift bins and above a given stellar mass; i.e., $\rho_{\text{dust}}(M_{*} > M, z)$. Because these are cumulative stacked stamps, individual panels are not independent. To ease the assessment of the detection significance, the color-scale of each stamp is set to vary as $S_{\text{pix}}/N_{\text{pix}}$, where $S_{\text{pix}}$ is the pixel signal and $N_{\text{pix}}$ is the standard deviation of the pixel distribution (both in units of $M_{\odot}$ Mpc$^{-3}$ beam$^{-1}$). The number of sources stacked (Nb) is indicated in each stamp, while in parenthesis is the number of sources among them which are individually detected at 1.2 mm. At $z < 3.2$, we only show stellar masses not affected by incompleteness (see Figure 1). At $z > 3.2$, where our H-band-selected sample cannot be considered as stellar mass-selected, we stacked all $z > 3.2$ galaxies together, irrespective of their stellar masses.</figcaption></figure>

than high-redshift galaxies because of the so-called negative K-correction (see, e.g., Figure 2 in Scoville et al. 2016).

At $z > 3.2$, where the sample cannot be considered as stellar mass-complete, we divided it into two redshift bins, $3.2 \le z < 4.5$ and $4.5 \le z < 5.5$, and stacked all galaxies with $H \le 27$ (see Figure 2). Due to stellar mass-incompleteness, these measurements only provide upper limits on the cosmic dust and gas mass densities in galaxies at these redshifts. However, because our 1.2 mm to gas mass conversion requires accurate metallicity and consequently stellar mass measurements, the inferred upper limits on the cosmic gas mass density at $z > 3.2$ are affected by systematics not included in our total uncertainties.

Figure 3 shows the differential stacked stamps for the dust; i.e., the cosmic dust mass density in galaxies in a given redshift and stellar mass bins, $\rho_{\text{dust}}(M_{*} \in [M \pm 0.25\text{dex}], z)$. In most redshift bins, we obtain significant detections for our massive stellar mass bins while low-mass galaxies are mostly undetected. This implies that, at our stellar mass-completeness limits, the significant detections observed in our cumulative stacked stamps are dominated by the contribution of massive galaxies ($\gtrsim 10^{9.5-10} M_{\odot}$). The fact that $\lesssim 10^{9.5} M_{\odot}$ galaxies only mildly contribute to the cosmic dust (and gas) mass density in galaxies is the main result of the paper (see also Dunlop et al. 2017) and is further illustrated and discussed in Sections 5 and 6. We note that in fact at $0.3 \le z < 0.6$, $0.6 \le z < 1.0$,
<table><caption>**Table 2**
Same as Table 1 But for the Cosmic Gas Mass Density in Galaxies</caption><thead><tr><th rowspan="2">Redshift</th><th colspan="7">ρ<sub>gas</sub>(M<sub>*</sub> > M, z)</th></tr><tr><th colspan="7">[×10<sup>7</sup> M<sub>☉</sub> Mpc<sup>-3</sup>]</th></tr><tr><th>Bin</th><th>&gt;10<sup>8</sup> M<sub>☉</sub></th><th>&gt;10<sup>8.5</sup> M<sub>☉</sub></th><th>&gt;10<sup>9</sup> M<sub>☉</sub></th><th>&gt;10<sup>9.5</sup> M<sub>☉</sub></th><th>&gt;10<sup>10</sup> M<sub>☉</sub></th><th>&gt;10<sup>10.5</sup> M<sub>☉</sub></th><th>&gt;10<sup>11</sup> M<sub>☉</sub></th></tr></thead><tbody><tr><td rowspan="3">0.3 &le; z &lt; 0.6</td><td>Nb = 8</td><td>Nb = 3</td><td>Nb = 2</td><td>Nb = 1</td><td>...</td><td>...</td><td>...</td></tr><tr><td>1.8 &plusmn; 1.6 (1.3)</td><td>1.4 &plusmn; 1.0 (0.4)</td><td>1.0 &plusmn; 0.8 (0.3)</td><td>0.3 &plusmn; 0.3 (0.2)</td><td>[0.2<sup>+0.2</sup><sub>-0.2</sub>]</td><td></td><td></td></tr><tr><td>[2.0<sup>+1.1</sup><sub>-1.1</sub>]</td><td></td><td>[0.8<sup>+0.3</sup><sub>-0.4</sub>]</td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="3">0.6 &le; z &lt; 1.0</td><td>...</td><td>Nb = 50</td><td>Nb = 30</td><td>Nb = 18</td><td>Nb = 12</td><td>Nb = 8</td><td>Nb = 2</td></tr><tr><td>...</td><td>3.8 &plusmn; 1.7 (1.5)</td><td>3.3 &plusmn; 1.2 (1.0)</td><td>2.2 &plusmn; 0.8 (0.5)</td><td>1.4 &plusmn; 0.6 (0.3)</td><td>0.6 &plusmn; 0.3 (0.2)</td><td>0.2 &plusmn; 0.2 (0.1)</td></tr><tr><td>[5.0<sup>+2.1</sup><sub>-2.6</sub>]</td><td></td><td>[3.3<sup>+0.8</sup><sub>-1.3</sub>]</td><td></td><td>[1.4<sup>+0.5</sup><sub>-0.5</sub>]</td><td></td><td></td></tr><tr><td rowspan="3">1.0 &le; z &lt; 1.6</td><td>...</td><td>Nb = 100</td><td>Nb = 48</td><td>Nb = 27</td><td>Nb = 13</td><td>Nb = 6</td><td>Nb = 2</td></tr><tr><td>...</td><td>6.4 &plusmn; 2.0 (1.5)</td><td>5.3 &plusmn; 1.5 (0.6)</td><td>5.1 &plusmn; 1.4 (0.4)</td><td>3.3 &plusmn; 1.2 (0.2)</td><td>2.3 &plusmn; 1.1 (0.1)</td><td>1.5 &plusmn; 1.1 (0.07)</td></tr><tr><td>[5.6<sup>+3.3</sup><sub>-1.7</sub>]</td><td></td><td>[5.4<sup>+1.7</sup><sub>-1.4</sub>]</td><td></td><td>[4.1<sup>+0.9</sup><sub>-1.7</sub>]</td><td></td><td></td></tr><tr><td rowspan="3">1.6 &le; z &lt; 2.3</td><td>...</td><td>...</td><td>Nb = 27</td><td>Nb = 15</td><td>Nb = 10</td><td>Nb = 4</td><td>...</td></tr><tr><td>...</td><td>...</td><td>5.2 &plusmn; 1.5 (0.5)</td><td>4.5 &plusmn; 1.4 (0.3)</td><td>3.9 &plusmn; 1.3 (0.2)</td><td>2.4 &plusmn; 1.2 (0.1)</td><td>...</td></tr><tr><td>[5.6<sup>+1.8</sup><sub>-1.2</sub>]</td><td></td><td>[5.0<sup>+0.9</sup><sub>-0.7</sub>]</td><td></td><td>[3.6<sup>+0.8</sup><sub>-0.6</sub>]</td><td></td><td></td></tr><tr><td rowspan="3">2.3 &le; z &lt; 3.2</td><td>...</td><td>...</td><td>Nb = 27</td><td>Nb = 15</td><td>Nb = 8</td><td>Nb = 2</td><td>...</td></tr><tr><td>...</td><td>...</td><td>5.6 &plusmn; 1.9 (0.5)</td><td>5.3 &plusmn; 1.8 (0.2)</td><td>5.1 &plusmn; 1.8 (0.1)</td><td>1.4 &plusmn; 1.0 (0.04)</td><td>...</td></tr><tr><td>[9.2<sup>+2.9</sup><sub>-3.9</sub>]</td><td></td><td>[6.6<sup>+0.1</sup><sub>-1.8</sub>]</td><td></td><td>[3.5<sup>+0.8</sup><sub>-0.7</sub>]</td><td></td><td></td></tr><tr><td>3.2 &le; z &lt; 4.5</td><td colspan="3">Nb = 44</td><td colspan="3">2.4 &plusmn; 1.6 (1.3)</td></tr><tr><td>4.5 &le; z &lt; 5.5</td><td colspan="3">Nb = 9</td><td colspan="3">1.3 &plusmn; 1.8 (1.7)</td></tr></tbody></table>

1.0 $\le z < 1.6$, $1.6 \le z < 2.3$, and $2.3 \le z < 3.2$, the sources individually detected in the ASPECS LP 1.2 mm map contribute for about 30%, 20%, 70%, 80%, and 80% of the cumulative measurements at our stellar mass-completeness limits, respectively.

Finally, to check the robustness of these results, we repeated these measurements using a bootstrap analysis. For this, we made 200 realizations of the stacking analysis, using for each realization a different sample, drawn from the original one, with the same number of sources but allowing for replacement (i.e., a galaxy can be picked several times). We found that the mean values and standard deviations over these 200 realizations were fully consistent with the values and total uncertainties quoted in Tables 1 and 2. Using the same approach, we also test our method against positional offset between the H-band and 1.2 mm emission. In each bootstrap iteration, we randomly draw the position of each sources within 2D Gaussian functions centered on their original position and with a standard deviation of 0.6 (typical H-band to ALMA offset for $z \sim 2$ star-forming galaxies; Elbaz et al. 2018). Again, the mean values and standard deviations inferred over these 200 realizations were fully consistent with those quoted in Tables 1 and 2. As a last sanity check, we randomized the position of the galaxies in the sample and repeated the stacking analysis. We only obtained non-detections.

## 5. The Cosmic Dust Mass Density in Galaxies

### 5.1. $\rho_{\text{dust}}(M_{*} > M, z)$ versus $M_{*}$

The evolution of the comoving dust mass density in galaxies with stellar masses $>M$, i.e., $\rho_{\text{dust}}(M_{*} > M, z)$, as derived from the stacking above, is shown in Figure 4 and tabulated in Table 1. At $z > 0.6$, $\rho_{\text{dust}}(M_{*} > M, z)$ grows rapidly as $M$ decreases down to $\sim 10^{10} M_{\odot}$, but this growth significantly slows down as $M$ decreases to even lower stellar

masses. At $z < 0.6$, the measurements are unfortunately too uncertain to fully confirm the existence of such a trend. The flattening of $\rho_{\text{dust}}(M_{*} > M, z)$ at low $M$, which happens well above our stellar mass-completeness limits, implies that: (i) at our stellar mass-completeness limits our analysis already accounts for most of the dust in the universe locked in galaxies; and (ii) the contribution of low-mass galaxies ($\lesssim 10^{9.5-10} M_{\odot}$) to the total cosmic dust mass density in galaxies becomes rapidly negligible. This latter finding is clearly illustrated by the differential measurements (i.e., $\rho_{\text{dust}}(M_{*} \in [M \pm 0.25\text{dex}], z)$), which peaks around $10^{10-10.5} M_{\odot}$ in most redshift bins. This characteristic stellar mass of $10^{10-10.5} M_{\odot}$ where most of the dust in galaxies is locked, is consistent with the characteristic stellar mass of star-forming galaxies where most of new stars were formed out to $z \sim 3$ ($10^{10.6 \pm 0.4} M_{\odot}$; Karim et al. 2011). From a more observational point of view, we note that our results are also consistent with the flattening of the cumulative 1.2 mm number counts found by González-López et al. (2020) at the unparalleled depth of the ASPECS LP 1.2 mm survey.

In the sample, about 10%, 20%, 5%, 5%, and 6% of the galaxies at $0.3 \le z < 0.6$, $0.6 \le z < 1.0$, $1.0 \le z < 1.6$, $1.6 \le z < 2.3$, and $2.3 \le z < 3.2$, are classified as quiescent using a standard UVJ selection method (e.g., Mortlock et al. 2015), respectively. Excluding these galaxies from the stacking analysis does not significantly change our $\rho_{\text{dust}}(M_{*} > M, z)$ estimates. At our stellar mass-completeness limits, $\rho_{\text{dust}}(M_{*} > M, z)$ decreases by $< 5 \%$ in the first three redshift bins and by $\sim 13 \%$ in the highest two. Considering that part of this decrease can actually be due to dusty star-forming galaxies contaminating the quiescent UVJ selection (Mortlock et al. 2015; Schreiber et al. 2015), we conclude that the bulk of the dust in galaxies resides in star-forming galaxies.

Given that ours is the first study, to our knowledge, to constrain the evolution of $\rho_{\text{dust}}(M_{*} > M, z)$ with stellar mass,
<figure><img src="image_4.png" /><figcaption>Figure 3. Same as Figure 2 but for the differential stacked stamps; i.e., the cosmic dust mass density in galaxies in a given redshift and stellar mass bin, $\rho_{\text{dust}}(M_{*} \in [M \pm 0.25 \text{ dex}], z)$. In the first row, the second and third panels look very similar. It corresponds to a close pair of galaxies, both at $z = 0.458$ and separated by only $1''7$, and falling into different stellar mass bins.</figcaption></figure>

we can thus only compare our results to the redshift evolution of the total $\rho_{\text{dust}}$ in galaxies<fnref n="32" /> derived by Dunne et al. (2003, 2011) and Driver et al. (2018). For clarity, in Figure 4 we displayed those $\rho_{\text{dust}}$ measurements at a stellar mass of $10^{7.6} M_{\odot}$. Driver et al. (2018) measured the $\rho_{\text{dust}}$ by fitting the optical-to-far-infrared photometry of 200,000 galaxies using the energy-balance code MAGPHYS. Although these estimates relied mostly on dust masses extrapolated from optical dust extinction due to the relatively limited depth of the Herschel observations used by Driver et al., they are in very good agreement with our measurements. Similarly, we find good agreement with Dunne et al. (2003, 2011), who measured $\rho_{\text{dust}}$ by integrating dust mass functions constrained from ground-based single-dish submillimeter observations and assuming $T = 25$ K and $\beta = 2$. This agreement is somewhat surprising, considering that the faint-end slopes of these dust mass functions at high-redshifts were only loosely constrained by those observations and thus fixed to that observed at $z < 0.1$.

The flattening of $\rho_{\text{dust}}(M_{*} > M, z)$ toward low stellar masses, together with the agreements with the total $\rho_{\text{dust}}$ found

by Dunne et al. (2003, 2011) and Driver et al. (2018), suggest that at the stellar mass-completeness limits of our study, we have already accounted for most of the dust in the universe locked in galaxies.

## 5.2. The Dust-to-stellar Mass Ratio of Star-forming Galaxies

At a given redshift, the variations of $\rho_{\text{dust}}(M_{*} > M, z)$ with $M$ can be modeled using the stellar mass function (i.e., SMF($M, z$)) and the mean dust-to-stellar mass ratio of galaxies (i.e., $f_{\text{dust}}(M, z)$),

$$ \rho_{\text{dust}}(M_{*} > M, z) = \int_{M}^{\infty} f_{\text{dust}}(M, z) \times \text{SMF}(M, z) dM. \quad (6) $$

Because the SMF of galaxies is well known up to $z \sim 3$ (e.g., Mortlock et al. 2015), one can solve for $f_{\text{dust}}(M, z)$ by fitting $\rho_{\text{dust}}(M_{*} > M, z)$. We excluded from the fits the measurements below our stellar mass-completeness limits and used the SMF of star-forming galaxies inferred by Mortlock et al. (2015);<fnref n="33" />

<footnote n="32">³² As opposed to literature measurements that could include a significant contribution from the dust in the circumgalactic medium (CGM), such as De Bernardis & Cooray 2012, Ménard and Fukugita (2012), Thacker et al. (2013).</footnote>

<footnote n="33">³³ The SMF inferred from our H-band-selected sample is consistent with that inferred by Mortlock et al. (2015). However, this H-band-selected sample is too small to robustly re-derive the SMF up to $z \sim 3$. Thus, we instead used that from Mortlock et al. (2015).</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Comoving dust mass density in galaxies for different redshift bins. Red and dark gray circles correspond to the cosmic dust mass densities in galaxies with stellar masses >M (i.e., $\rho_{\text{dust}}(M_{*} > M, z)$; cumulative stacking), above and below our stellar mass-completeness limits (i.e., $M_{\text{limit}}$; see Section 3), respectively. Yellow and light-gray circles show the cosmic dust mass densities in galaxies with stellar masses $\in [M \pm 0.25 \text{ dex}]$ (i.e., $\rho_{\text{dust}}(M_{*} \in [M \pm 0.25 \text{ dex}], z)$; differential stacking), above and below our stellar mass-completeness limits, respectively. Some of these differential data points have S/N < 3 (see Figure 3). The black-solid lines show the best-fits of $\rho_{\text{dust}}(M_{*} > M, z)$ using the SMF of Mortlock et al. (2015) and solving for $f_{\text{dust}}$ assuming $f_{\text{dust}} = M_{\text{dust}}/M_{*} = A \times (M/10^{10.7} M_{\odot})^{B}$ (see Section 5.2). The black dashed lines show the best-fits of $\rho_{\text{dust}}(M_{*} > M, z)$ when all redshift bins are fitted simultaneously solving for $f_{\text{dust}}$ assuming $\log(f_{\text{dust}}) = C + D \times \log(1+z) + B \times (\log(M_{*}/M_{z}) - 10.7)$. The light-gray solid and dashed lines show the exact same best-fits but displayed in differential form; i.e., within stellar mass bins which are 0.5 dex wide. Blue diamonds show predictions for galaxies with $>10^{10} M_{\odot}$ using $f_{\text{gas}}(M, z)$ for main-sequence galaxies from Tacconi et al. (2018), the SMF of Mortlock et al. (2015), and assuming a gas-to-dust ratio of 100. Green stars present the total comoving dust mass density in galaxies measured by Driver et al. (2018), applying the energy-balance code MAGPHYS to hundreds of thousands of galaxies in the GAMA/G10-COSMOS/3D-HST surveys. The black square show the total comoving dust mass density of galaxies measured by Dunne et al. (2003) using single-dish (sub)millimeter-selected galaxies.</figcaption></figure>

i.e.,

$$ \text{SMF}(M, z) = \phi^* \cdot \ln(10) \cdot \left( \frac{M}{M^*} \right)^{1+\alpha} \cdot e^{-M/M^*}, \quad (7) $$

where $\phi^*$ is the normalization of the Schechter function, $M^*$ is its turnover mass, and $\alpha$ is its low-mass end slope (Table 3). We did not use the SMF that includes quiescent galaxies because their contribution to $\rho_{\text{dust}}(M_{*} > M, z)$ has been shown to be negligible (see Section 5.1). $f_{\text{dust}}(M, z)$ is thus the mean stellar-to-dust mass ratio of star-forming galaxies.

First, we fitted each redshift bin independently, assuming that $f_{\text{dust}}(M, z)$ follows a simple power-law,

$$ f_{\text{dust}}(M, z) = A \times \left( \frac{M}{10^{10.7} M_{\odot}} \right)^{B}, \quad (8) $$

where the choice of a $10^{10.7} M_{\odot}$ normalization allows direct comparisons with Tacconi et al. (2018). With this parameterization, A has an immediate physical interpretation; i.e., it is the typical dust-to-stellar mass ratios of star-forming galaxies with a stellar mass of $5 \times 10^{10} M_{\odot}$. The results of these fits are

<table><caption>**Table 3**
The Single Schechter Parameters for the Star-forming Galaxy SMFs, as Found in Mortlock et al. (2015)</caption><thead><tr><th>Redshift Bin</th><th>log M</th><th>log φ*</th><th>α</th></tr></thead><tbody><tr><td>0.3 &lt; z &lt; 0.5</td><td>10.83</td><td>-3.31</td><td>-1.41</td></tr><tr><td>0.5 &lt; z &lt; 1.0</td><td>10.77</td><td>-3.28</td><td>-1.45</td></tr><tr><td>1.0 &lt; z &lt; 1.5</td><td>10.64</td><td>-3.14</td><td>-1.37</td></tr><tr><td>1.5 &lt; z &lt; 2.0</td><td>11.01</td><td>-4.05</td><td>-1.74</td></tr><tr><td>2.0 &lt; z &lt; 2.5</td><td>10.93</td><td>-3.93</td><td>-1.77</td></tr><tr><td>2.5 &lt; z &lt; 3.0</td><td>11.08</td><td>-4.41</td><td>-1.92</td></tr></tbody></table>

shown by the black and light gray lines in Figure 4, while the redshift evolution of A and B are shown by the red circles in Figure 5. In each redshift bin, our fit provides an accurate description of our observations, with reduced $\chi^2$ in the range 0.5–1.0. The exponent of the dust-to-stellar mass ratio is found to be negative ($B < 0$) in the first two redshift bins but becomes positive at higher redshifts. However, the uncertainties associated with these exponents render all of them consistent with zero at all redshifts.
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Redshift evolution of the normalization (top panel) and exponent (bottom panel) of the dust-to-stellar mass ratio of star-forming galaxies modeled as a simple power-law of the stellar mass (see inset equation). A is the typical dust-to-stellar mass ratio of star-forming galaxies with a stellar mass of $5 \times 10^{10} M_{\odot}$. Red circles show our constraints while fitting $\rho_{\text{dust}}(M_{*} > M, z)$ in each redshift bin independently using Equation (6). Dashed lines show our constraints while fitting all redshift bins simultaneously and modeling the dust-to-stellar mass ratio as a simple power-law of the stellar mass, with a redshift-independent exponent and a redshift-dependent normalization (see Equation (9)). Gray regions present the range of fits compatible within $1\sigma$ with our observations.</figcaption></figure>

We then solved for $f_{\text{dust}}(M, z)$ by fitting $\rho_{\text{dust}}(M_{*} > M, z)$ but fitting all redshift bins simultaneously and assuming that the exponent does not vary with redshift. For the redshift-dependent normalization, we used a parameterization suggested by Tacconi et al. (2018); i.e.,

$$ \log(f_{\text{dust}}(M, z)) = C + D \times \log(1 + z) + B \times \log(M/(10^{10.7} M_{\odot})). \quad (9) $$

The best-fit is shown by the black dashed lines in Figure 4, while the corresponding normalization (i.e., $\log(A) = C + D \times \log(1 + z)$; with $C = -3.0_{-0.5}^{+0.3}$, and $D = 2.6_{-1.1}^{+1.1}$) and exponent (i.e., $B = 0.1_{-0.7}^{+1.0}$) of the dust-to-stellar mass ratio are shown in Figure 5. This fit accurately described our observations, with a reduced $\chi^2 = 1.0$. The mean dust-to-stellar mass ratio in star-forming galaxies is found to significantly vary with redshift. Taking into account all of the compatible models, within $1\sigma$ with our observations, we found that the mean dust-to-stellar mass ratio of $10^{10.7} M_{\odot}$ star-forming galaxy increases by at least a factor of 2 and at most a factor of 60 between $z = 0.45$ and 2.75. Our best-fit implies an increase by a factor 10, which should be compared to the factor 12 increase of the gas-to-stellar mass ratio found by Tacconi et al. (2018) within this redshift range.

Finally, from the gas-to-stellar mass ratio of star-forming galaxies derived in Tacconi et al. (2018), we predicted $\rho_{\text{dust}}(M_{*} > M, z)$ at $M = 10^{10} M_{\odot}$, assuming a gas-to-dust mass ratio of 100, which is typical for massive galaxies at $z \sim 0$ (Leroy et al. 2011). These predictions are in good

agreement, within the total uncertainties, with our observations, but in our lowest redshift bin (see blue diamonds in Figure 4). Thus, even though at high stellar masses our analysis is affected by the low number of galaxies available within our pencil-beam survey (i.e., cosmic variance), our measurements agree with those inferred using larger, more representative samples of massive star-forming galaxies.

### 5.3. $\rho_{\text{dust}}$ versus Redshift

The redshift evolution of $\rho_{\text{dust}}(M_{*} > M, z)$ at our stellar mass-completeness limits is shown in Figure 6; i.e., $\rho_{\text{dust}}(M_{*} > M_{\text{limit}}, z)$. Along with these measurements, we also included estimates for $M_{*} > 10^8 M_{\odot}$, $M_{*} > 10^9 M_{\odot}$, and $M_{*} > 10^{10} M_{\odot}$, inferred by fitting $\rho_{\text{dust}}(M_{*} > M, z)$ in each redshift bin independently using the method described in Section 5.2. As advocated in Section 5.1, the total cosmic dust mass density in galaxies (i.e., $\rho_{\text{dust}}$) should be well approximated by these $M_{*} > 10^8 M_{\odot}$ measurements.

Our analysis suggests that $\rho_{\text{dust}}$ did not evolve much from $z = 2.75$ to 1.3 but decreased by a factor $\sim 3.6(\pm 2.0)$ from $z = 1.3$ to 0.45. As noted in Section 5.1 and Figure 4, this redshift evolution is consistent with that inferred by Dunne et al. (2003, 2011) and Driver et al. (2018). It also broadly agrees with recent measurements by Pozzi et al. (2019), which were obtained using Herschel observations in the COSMOS field. We only notice a significant disagreement with this later study at $z > 2$; i.e., a redshift range where their observations mostly constrain the bright-end slope of the dust mass function.

From $z = 0.45$ to the present time, $\rho_{\text{dust}}$ did not evolve significantly, because our $z = 0.45$ measurement is already equal to the local ($z = 0.05$) cosmic dust mass density in galaxies constrained by Dunne et al. (2011). We note, however, that these $z < 0.45$ measurements rely mostly on Herschel-250 µm observations that probe the dust peak emission of galaxies, while our measurements rely on their dust Rayleigh-Jeans emission. These two approaches might thus be affected by different systematics, which renders their combination difficult.

From early cosmic time to $z \sim 3$, simulations predict a drastic increase of $\rho_{\text{dust}}$ (see Section 7.1) but further investigations are needed to confirm this trend observationally, such as by using even longer wavelength selected samples (i.e., $\lambda_{\text{obs}} > 2$ µm) that are still sensitive to stellar masses at these redshifts.

The redshift evolution of $\rho_{\text{dust}}$ resembles that of the SFRD ($\rho_{\text{SFR}}$) of the universe (e.g., Madau & Dickinson 2014). To investigate this further, we show in Figure 6 the redshift evolution of $\rho_{\text{SFR}}$ (Madau & Dickinson 2014) scaled assuming a (molecular) gas depletion time ($t_{\text{depl}} = M_{\text{gas}}/\text{SFR}$) of 570 Myr and a gas-to-dust mass ratio ($\delta_{\text{GDR}}$) of 150, i.e.,

$$ \rho_{\text{dust}} = \rho_{\text{SFR}} \times t_{\text{depl}} \times \delta_{\text{GDR}}^{-1}. \quad (10) $$

These particular values of $t_{\text{depl}}$ and $\delta_{\text{GDR}}$ were chosen because they are typical for star-forming galaxies at $z \sim 2$ with a stellar mass of $10^{10.3} M_{\odot}$ (for $t_{\text{depl}}$, Tacconi et al. 2018, for $\delta_{\text{GDR}}$, Equations (2) and (3)); i.e., the characteristic stellar mass where most of the dust in galaxies is locked (Section 5.1). While these predictions describe reasonably well our measurements at $z > 0.45$, they underestimate the observations at low redshifts. A flatter evolution of $\rho_{\text{dust}}$ is predicted and thus a better match to the low-redshift measurements is obtained, when assuming a
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Evolution with look-back time (lower x-axis; or redshift, upper x-axis) of the comoving dust mass density in galaxies derived here using the ASPECS LP 1.2 mm continuum map. Red circles are for galaxies with $M_{*} > M_{\text{limit}}$, as inferred from Figure 4. Red, yellow, and blue hashed regions are inferred from the best-fits of $\rho_{\text{dust}}(M_{*} > M, z)$ (see Figure 4) and correspond to the cosmic dust mass densities in galaxies with $M_{*} > 10^{8} M_{\odot}$, $M_{*} > 10^{9} M_{\odot}$, and $M_{*} > 10^{10} M_{\odot}$, respectively. These best-fits were performed independently for each redshift bin. At $z > 3$, we derive lower limits by stacking all galaxies in our H-band-selected sample, irrespective of their stellar masses. At a similar redshift, we show the lower limit inferred by Magnelli et al. (2019) using a 2 mm selected galaxy sample. Black squares, green stars, and dark blue diamonds show the total comoving dust mass densities in galaxies inferred by Dunne et al. (2003, 2011), Driver et al. (2018), and Pozzi et al. (2019), respectively. The solid line shows predictions from scaling the cosmic SFRD (Madau & Dickinson 2014) assuming a redshift-independent gas depletion timescale ($t_{\text{depl}} = M_{\text{gas}}/\text{SFR}$) of 570 Myr and a gas-to-dust mass ratio of 150; both values being typical for $10^{10.3} M_{\odot}$ main-sequence star-forming galaxies at $z \sim 2$ (Tacconi et al. 2018). The dotted line shows predictions assuming a gas-to-dust mass ratio of 150 but a redshift-dependent gas depletion time as parametrized by Tacconi et al. (2018) for $10^{10.3} M_{\odot}$ main-sequence star-forming galaxies (i.e., $t_{\text{depl}} = 570 \text{ Myr}$ at $z = 2$ and 860 Myr at $z = 0$). Finally, the dashed line shows predictions assuming both a redshift-dependent depletion time and gas-to-dust mass ratio. This latter is derived from Equation (2) using the metallicity of $10^{10.3} M_{\odot}$ galaxies at a given redshift (Equation (3)). This yields $\delta_{\text{GDR}} = 150$ at $z = 2$ and $\delta_{\text{GDR}} = 90$ at $z = 0$. The typical fractional cosmic variance uncertainty of $\sim 45\%$ affecting our measurements and inferred using the methodology advocated by Driver & Robotham (2010) is illustrated in the upper right-hand corner.</figcaption></figure>

more realistic redshift-dependent gas depletion time (dotted line in Figure 6), parametrized using the results from Tacconi et al. (2018) for $10^{10.3} M_{\odot}$ star-forming galaxies (i.e., $t_{\text{depl}} = 570 \text{ Myr}$ at $z = 2$ and 860 Myr at $z = 0$). Finally, assuming both a redshift-dependent depletion time and gas-to-dust mass ratio (dashed line in Figure 6), we predict an even flatter evolution of $\rho_{\text{dust}}$, which matches reasonably well with our measurements. This flatter evolution illustrates the fact that at fix stellar mass (here $10^{10.3} M_{\odot}$), the mean metallicity of galaxies increases from $z \sim 2$ to $z \sim 0$ (see Equation (3)), which implies that their mean gas-to-dust mass ratio decreases over this redshift range ($\delta_{\text{GDR}} = 150$ at $z = 2$ and $\delta_{\text{GDR}} = 90$ at $z = 0$).

The flat evolution of $\rho_{\text{dust}}$ at $z < 0.45$ could also in part be due to the increasing contribution of the atomic phase of the ISM as we approach $z = 0$. This increasing contribution is not taken into account by our toy model because the values for $t_{\text{depl}}$ in Equation (10) were taken from Tacconi et al. (2018) and only include the molecular gas phase. At low redshifts, a significant contribution of the dust in the atomic phase to the total observed dust mass was actually reported by Tacconi et al. (2018), when comparing CO-based and dust-based gas mass estimates.

As a final remark, we note that our estimates of $\rho_{\text{dust}}$ are affected by cosmic variance, which is only partly included in the Poissonian uncertainties. Using the method described by Driver & Robotham (2010), we estimate for our $2.27 \text{ arcmin}^2$ survey a fractional cosmic variance uncertainty of 58%, 49%, 42%, 43% and 42% at $0.3 \le z < 0.6$, $0.6 \le z < 1.0$, $1.0 \le z < 1.6$, $1.6 \le z < 2.3$, and $2.3 \le z < 3.2$, respectively. These uncertainties are not significantly larger than the total uncertainties quoted in Table 1 and displayed in Figure 6.

## 6. The Cosmic Gas Mass Density in Galaxies

### 6.1. $\rho_{\text{gas}}(M_{*} > M, z)$ versus $M_{*}$

The evolution of the comoving gas (H$_2$+H I) mass density in galaxies with stellar masses $> M$—i.e., $\rho_{\text{gas}}(M_{*} > M, z)$—is shown in Figure 7 and tabulated in Table 2. These measurements were inferred by stacking the 1.2 mm emission of these galaxies (Section 4.2) and assuming a metallicity-dependent gas-to-dust mass ratio relation (Section 4.1).

As for the cosmic dust mass density, at $z > 0.6$, $\rho_{\text{gas}}(M_{*} > M, z)$ grows rapidly as $M$ decreases to $\sim 10^{10} M_{\odot}$, and this growth slows down as $M$ decreases to even lower stellar masses. This flattening at low stellar masses suggests once again that at the stellar mass-completeness limits of our
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Same as Figure 4, but for the comoving gas mass density in galaxies with stellar masses >M. Green stars present the *total* comoving molecular gas mass density in galaxies measured by the ASPECS-CO pilot and LP surveys (Decarli et al. 2016, 2019). The black square shows the *total* comoving molecular gas mass density in galaxies measured by the COLDz survey (Riechers et al. 2019). Blue triple-dotted–dashed lines show predictions using $f_{\text{gas}}(M,z)$ for main-sequence galaxies from Tacconi et al. (2018) and the SMF of Mortlock et al. (2015).</figcaption></figure>

sample, the gas mass density measured here already accounts for most of the gas content locked in galaxies and converges toward the *total* cosmic gas mass density in galaxies. We note, however, that the slope of $\rho_{\text{gas}}(M_{*} > M, z)$ at low stellar masses is slightly steeper than that inferred for $\rho_{\text{dust}}(M_{*} > M, z)$. This difference is explained by the decrease of the metallicity and therefore the increase of the gas-to-dust mass ratio at low stellar masses (Equation (2) and (3)). Consequently, the peak of our differential measurements ($\rho_{\text{gas}}(M_{*} \in [M \pm 0.25\text{dex}], z)$) is much broader or somewhat washed out at $z < 1.0$.

Excluding quiescent galaxies from the stacking analysis barely affects our results, decreasing the value of $\rho_{\text{gas}}(M_{*} > M, z)$ at our stellar mass-completeness limits by at most 10%. As in the case of dust, the bulk of the gas in galaxies appears to reside in star-forming galaxies (see also Sargent et al. 2015; Gobat et al. 2018).

Decarli et al. (2016, 2019) and Riechers et al. (2019) measured the *total* molecular gas mass density in galaxies (i.e., $\rho_{\text{H}_2}$) by constraining the CO luminosity function using the ASPECS-CO pilot/LP and COLDz surveys. In all redshift bins, we find a good agreement between our measurements and their estimates. This suggests that even though our gas mass measurements include in principle both the molecular and atomic gas phases, they are mostly dominated by the molecular phase. We note, however, that our method implicitly assumes that the dust emissivity in the atomic and molecular gas phase is the same, while there exists observational evidence of an

enhanced dust emissivity in the dense ISM (see Leroy et al. 2011, and reference therein). This implies that dust-based gas mass estimates might be biased against dust in the atomic phase.

From the gas-to-stellar mass ratio of star-forming galaxies derived in Tacconi et al. (2018), we predicted $\rho_{\text{gas}}(M_{*} > M, z)$ using the SMF of Mortlock et al. (2015). Up to $z \sim 2.0$, these predictions agree, within the uncertainties, with our observations, but in the highest redshift bin they significantly underestimate our measurements (see triple-dotted–dashed lines in Figure 7). At this redshift, the steep dependency with stellar mass of the gas-to-stellar mass ratio found in Tacconi et al. (2018, $B = -0.36$) seems to overestimate the contribution of low stellar mass galaxies to the cosmic gas mass density.

Finally, we note that using a $\delta_{\text{GDR}}$–metallicity relation with a much steeper power-law at metallicity $<7.9$ (Rémy-Ruyer et al. 2014) only affects our estimates below our stellar mass-completeness limits (see Appendix A). This steep power law, which implies much larger gas mass per unit dust mass at low metallicity, increases our measurements by a factor of 2–5 in our lowest stellar mass bins. This leads to a very discontinuous evolution of $\rho_{\text{gas}}(M_{*} > M, z)$ with M, large disagreements with the ASPECS-CO survey's results, and could suggest an increasing contribution of the atomic phase to $\rho_{\text{gas}}(M_{*} > M, z)$ at low metallicities. However, because the metallicity below which this steep power law starts is still very
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Same as Figure 5 but for the gas-to-stellar mass ratio of star-forming galaxies modeled as a simple power-law of the stellar mass (see inset equation). A is the typical gas-to-stellar mass ratio of star-forming galaxies with a stellar mass of $5 \times 10^{10} M_{\odot}$.</figcaption></figure>

uncertain (Rémy-Ruyer et al. 2014), we do not discuss this effect further.

## 6.2. The Gas-to-stellar Mass Ratio of Star-forming Galaxies

Using the method described in Section 5.2, we model the gas-to-stellar mass ratio of star-forming galaxies as a simple power-law function of their stellar mass. First, we independently solve for the normalization and exponent of this power-law function in each redshift bin (thick-black line in Figure 7 and red dots in Figure 8). Then, we fit all our redshift bins simultaneously by assuming a redshift-independent exponent (thick-dashed line in Figure 7 and gray regions in Figure 8; $B = -0.1_{-0.5}^{+0.8}$, $C = -1.1_{-0.5}^{+0.3}$, and $D = 3.0_{-1.0}^{+1.1}$).

The best-fit model obtained from fitting all redshift bins simultaneously yields an exponent of -0.1. This tentative trend, in which massive galaxies have lower gas mass content per unit stellar mass than lower mass galaxies, is, nevertheless, not as steep as that found in massive high-redshift galaxies ($B \sim -0.36$; e.g., Magdis et al. 2012; Genzel et al. 2015; Tacconi et al. 2018). However, in the local universe and considering only the molecular gas phase, Saintonge et al. (2017) found a flatter evolution with stellar mass of the gas-to-stellar mass ratio at $<10^{10.2} M_{\odot}$. Finally, even though we independently constrained the $f_{\text{dust}}-M_{*}$ (Section 5.2) and $f_{\text{gas}}-M_{*}$ relations, these are linked via the gas-to-dust mass ratio (Equation (2)) and stellar mass-metallicity (Equation (3)) relations. By combining Equations (2) and (3), one can predict that the exponent of the $f_{\text{dust}}-M_{*}$ relation should be higher by 0.15-0.3 to that of the $f_{\text{gas}}-M_{*}$ relation.

The mean gas-to-stellar mass ratio in $10^{10.7} M_{\odot}$ star-forming galaxies increases by at least a factor of 2, at most a factor of 70, and for our best-fit a factor of 17 between $z = 2.75$ and 0.45. These values should be compared to the factor 12 increase found by Tacconi et al. (2018) within this redshift range.

## 6.3. $\rho_{\text{gas}}$ versus Redshift

Figure 9 presents the redshift evolution of $\rho_{\text{gas}}$ at our stellar mass-completeness limits, as well as our extrapolations for galaxies with $M_{*} > 10^{8} M_{\odot}$, $M_{*} > 10^{9} M_{\odot}$, and $M_{*} > 10^{10} M_{\odot}$, obtained by fitting $\rho_{\text{gas}}(M_{*} > M, z)$ in each redshift bin independently using the method described in Section 6.2. Because $\rho_{\text{gas}}(M_{*} > M, z)$ clearly flattens toward low stellar masses, extrapolations for $M_{*} > 10^{8} M_{\odot}$ galaxies, should provide a good measurement of the total cosmic gas mass density in galaxies.

The cosmic gas mass density in galaxies decreases by a factor $\sim 1.6(\pm 0.7)$ from $z = 2.75$ and 1.3 and then decreases by a factor $\sim 2.8(\pm 1.7)$ between $z = 1.3$ and 0.45. This redshift evolution is consistent with that inferred using the ASPECS-CO measurements (Decarli et al. 2016, 2019) and the COLDz survey (Riechers et al. 2019). At $z < 0.45$, the decrease of the cosmic gas mass density seems to continue, when considering the molecular gas mass density measured at $z \sim 0$ by Saintonge et al. (2017). However, considering instead the atomic gas mass density measured at $z \sim 0$ by Martin et al. (2010) yields an opposite trend, which illustrates the increased contribution of the atomic phase in the ISM of galaxies. Finally, combined with the ASPECS-CO constraints at $z > 3$, our measurements suggest a rapid increase of $\rho_{\text{gas}}$ from $z = 4$ to 2.75.

The evolution of the (molecular) gas mass densities of galaxies from $z \sim 4$ to $z \sim 0$ resembles that of $\rho_{\text{SFR}}$. As in Section 5.3, to study this further we plot in Figure 9 the redshift evolution of $\rho_{\text{SFR}}$ scaled assuming: (i) a redshift-independent gas depletion time of 570 Myr (solid line), and (ii) a more realistic redshift-dependent gas depletion time for $10^{10.3} M_{\odot}$ star-forming galaxies (dotted line; i.e., $t_{\text{depl}} = 570$ Myr at $z = 2$ and 860 Myr at $z = 0$; Tacconi et al. 2018). Both predictions match strikingly well with the observations from $z \sim 4$ to $z \sim 0$. This finding strongly suggests that: (i) at $z > 0.45$, dust-based gas mass measurements are mostly dominated by the dust in the molecular phase, and (ii) the redshift evolution of the cosmic SFRD is mostly explained by the evolution of the molecular gas reservoir of galaxies (solid line), with variations of their star formation efficiency playing a secondary role (dotted line).

The evolutions of $\rho_{\text{dust}}$ and $\rho_{\text{gas}}$ measured in our study are linked via the gas-to-dust mass ratio and stellar mass-metallicity relations. At a given redshift, the $\rho_{\text{gas}}$-to-$\rho_{\text{dust}}$ ratio correspond thus to the mass-weighted average gas-to-dust mass ratio of star-forming galaxies. At $z \sim 0.45$, $z \sim 0.80$, $z \sim 1.30$, $z \sim 1.95$, and $z \sim 2.75$, we found $166 \pm 115$, $155 \pm 104$, $130 \pm 73$, $140 \pm 63$, and $240 \pm 140$, respectively. We note that in Appendix A, we infer instead the redshift evolution of $\rho_{\text{gas}}$ scaling $\rho_{\text{dust}}$ by a constant gas-to-dust mass ratio of 100, typical for massive galaxies at $z \sim 0$. These measurements underestimate at the $1-2\sigma$ level that from the ASPECS-CO LP at $z = 1-2$.

# 7. Discussion

Our results reveal that up to $z \sim 3$, most of the dust and gas content in galaxies resides in star-forming galaxies with stellar masses $\gtrsim 10^{9.5} M_{\odot}$. Our cosmic gas mass density estimates agree with those inferred from CO observations, which suggests that dust-based measurements are dominated by the dust in the molecular phase. The total dust and gas mass densities in galaxies increase at early cosmic time, peak around
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Same as Figure 6 but for the comoving gas mass density in galaxies. Shaded regions present the $1\sigma$ confidence measurement for the molecular gas from the ASPECS-CO pilot and LP surveys (Decarli et al. 2016, 2019), as well as the 5th-to-95th percentile confidence interval from the COLDz survey (Riechers et al. 2019). At $z \sim 0$, the black triangle and diamond show the cosmic molecular and atomic gas mass densities in galaxies inferred by Saintonge et al. (2017) and Martin et al. (2010), respectively. The solid line shows predictions from scaling the cosmic SFRD (Madau & Dickinson 2014) assuming a redshift-independent gas depletion time ($t_{\text{depl}} = M_{\text{gas}}/\text{SFR}$) of 570 Myr, typical for $10^{10.3} M_{\odot}$ main-sequence star-forming galaxies at $z \sim 2$ (Tacconi et al. 2018). The dotted line shows predictions assuming a more realistic redshift-dependent gas depletion time for $10^{10.3} M_{\odot}$ main-sequence star-forming galaxies (i.e., $t_{\text{depl}} = 570$ Myr at $z = 2$ and 860 Myr at $z = 0$; Tacconi et al. 2018).</figcaption></figure>

$z \sim 1-3$, and then decrease until the present time. However, the dust and gas mass densities decrease at different rates: the former declines by a factor $\sim 4$ between $z \sim 2.5$ and 0 (combining our results with the Herschel-based local measurement of Dunne et al. 2011), while the later declines by a factor $\sim 7$ when only considering the molecular gas phase and using the local CO measurement of Saintonge et al. (2017). The redshift evolution of the cosmic dust and gas mass densities can be modeled by the redshift evolution of the stellar mass function of star-forming galaxies, and that of their dust-to-stellar mass ratio and gas-to-stellar mass ratio, respectively. These models show that the dust and gas content of star-forming galaxies per unit stellar mass continuously decrease from $z=3$ to $z=0$, while they have a mild dependency on stellar masses—with a best-fit power law stellar mass-dependent exponent of 0.1 and -0.1 for the dust-to-stellar mass ratio and gas-to-stellar mass ratio relations, respectively.

In the following subsections, we first compare these new results to the outputs of simulations and we will then put them into the context of galaxy evolution scenarios.

## 7.1. Comparison to Simulations

In the past years, there has been a growing interest in including self-consistent tracking of the production and destruction of dust in cosmological models of galaxy formation (e.g., McKinnon et al. 2017; Popping et al. 2017; Aoyama et al. 2018; Davé et al. 2019; Vijayan et al. 2019). This includes the condensation of dust in the ejecta of asymptotic giant branch

(AGB) stars and supernovae (SNe), the growth of dust in the ISM, as well as the destruction of dust by supernova-induced shocks, star formation (the so-called astraction), reheating, and outflows. Our analysis provides to these models some of the first robust and homogeneously constrained evolution of the comoving dust and gas mass densities in galaxies from $z \sim 3$ to $z \sim 0.4$.

In Figure 10, we compare the cosmic dust and gas mass densities predicted by Popping et al. (2017)'s semi-analytical model for different stellar mass thresholds. These predictions correspond to the dust and gas in the ISM of these simulated galaxies, as opposed to the dust and gas in their CGM and the intergalactic medium (IGM; see Popping et al. 2017, for details). In agreement with our observations, this model predicts that at any redshifts, the bulk of the dust in galaxies is locked in those with a stellar mass $\gtrsim 10^{10} M_{\odot}$. However, while these predictions match the observations from $z \sim 3$ to $z \sim 1$, they overestimate the cosmic dust mass density in galaxies at $z < 1$. In this model, from $z \sim 1$ to $z \sim 0$, the formation and destruction (or ejection) rates of the dust are in quasi-equilibrium, yielding a nearly constant cosmic dust mass density across this redshift range. Instead, observations suggest that this quasi-equilibrium is reached at a later time; i.e., $z < 0.5$. We note, however, that to confirm this trend, one would need to perform a more homogeneous observational analysis across this redshift range. Indeed, low-redshift measurements rely mostly on 250 µm Herschel observations that probe the dust peak emission of galaxies, while our measurements rely on their dust Rayleigh-Jeans emission.
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Evolution with look-back time of the observed comoving dust (top panel) and gas (bottom panel) mass densities in galaxies, compared to predictions from Popping et al. (2017)'s semi-analytical model. For the observed comoving dust and gas mass densities, the symbols are the same as in Figures 6 and 9, respectively. (top panel) Blue, yellow, red and black-solid lines show predictions for galaxies with $M_{*} > 10^{10} M_{\odot}$, $>10^{9} M_{\odot}$, $>10^{8} M_{\odot}$, and $>10^{7} M_{\odot}$, respectively. The dotted–dashed and triple-dotted–dashed lines show predictions from the cosmic dust mass densities in galaxies from the cosmological hydrodynamic simulations of Li et al. (2019; see also Davé et al. 2019) and Aoyama et al. (2018). The gray-shaded area shows the total amount of dust formed in galaxies, assuming a dust yield of 0.004–0.0065 dust masses for every unit of stellar mass formed (Driver et al. 2018) and using Madau & Dickinson (2014)'s cosmic SFRD history. (bottom panel) The dashed and solid lines show predictions from Popping et al. (2017) for the gas (H$_{2}$+H I) and molecular gas mass densities in galaxies, respectively. Lines are color-coded as in the top panel.</figcaption></figure>

These two approaches might thus be affected by different systematics, which renders their combination difficult.

In the top panel of Figure 10, we also overlaid the total cosmic dust mass densities in galaxies predicted by the cosmological hydrodynamic simulations of Aoyama et al. (2018) and Li et al. (2019, see also Davé et al. 2019). However, they are relatively inconsistent with our measurements, over-estimating at most redshifts the cosmic dust mass densities in galaxies.

In the bottom panel of Figure 10, we finally compare the cosmic gas (H$_{2}$+H I) and molecular gas mass densities predicted by Popping et al. (2017) for different stellar mass thresholds. Predictions for the molecular gas phase match reasonably well with the observations from $z \sim 3$ up to $z \sim 0$. In addition, they also correctly predict that most of the molecular gas mass is locked in galaxies with $\gtrsim 10^{9} M_{\odot}$. However, our dust-based measurements should in principle be compared to predictions including both the molecular and
atomic gas phases. In this case, Popping et al. (2017)'s predictions overestimate our measurements by a factor of 2 for galaxies with $M_* > 10^{10} M_\odot$ and a factor of 4 for galaxies with $M_* > 10^8 M_\odot$, illustrating the rising contribution in this model of the atomic gas phase in the ISM of low stellar mass galaxies. These over-estimations suggest that either the model over-estimates the atomic gas content locked in low stellar mass galaxies or that our measurements are biased against this atomic gas phase because of a significant enhancement of the dust emissivity in the dense/molecular ISM.

In Appendix B, we compare our measurements to predictions from the semi-analytical model of Vijayan et al. (2019) for $M_* > 10^8 M_\odot$ and $M_* > 10^9 M_\odot$. For the cosmic dust mass density in galaxies, these predictions are very similar to those of Popping et al. (2017); i.e., this model successfully predicts little contribution to $\rho_{\text{dust}}$ from $M_* < 10^9 M_\odot$ galaxies, but overestimates $\rho_{\text{dust}}$ at $z < 1$. For the cosmic gas mass density in galaxies, the predictions from Vijayan et al. (2019) differ from those of Popping et al. (2017). They overestimate $\rho_{\text{gas}}$ at $z \lesssim 1.5$ but successfully suggest that at high redshift most of the gas in galaxies is in the molecular phase.

## 7.2. Implications for Galaxy Evolution Scenarios

The redshift evolution of $\rho_{\text{dust}}$ and $\rho_{\text{gas}}$ resembles that of the cosmic SFRD, which also peaks at $z \sim 2$ (e.g., Madau & Dickinson 2014). This implies that there is a direct link between star formation, and the dust and gas content of galaxies. For the gas, this link is expected—gas fueling star formation—via the so-called Kennicutt–Schmidt (KS) relation, which connects the gas and star formation rate surface densities of galaxies. The relatively redshift-independent SFRD-to-$\rho_{\text{gas}}$ ratio observed here from $z \sim 3$ to $z \sim 0$ (see solid black line in Figure 9) suggests that the global star formation efficiency of galaxies (i.e., SFR/$M_{\text{gas}}$) does not evolve significantly ($\lesssim \times 2$) across this redshift range and that the main driver of star formation is gas content. At all redshifts, the time needed to deplete the global gas reservoir of the star-forming galaxy population ($\langle t_{\text{depl}} \rangle_{\text{V}} = \rho_{\text{gas}}$/SFRD) is found to be of the order of 600–900 Myr, which is consistent with results inferred from individual galaxies (Genzel et al. 2015; Tacconi et al. 2018). Without a constant replenishment of these gas reservoirs, the star-forming galaxy population observed at, for example, $z \sim 2$, would thus have fully disappeared by $z = 1.5$.

These findings strongly support gas regulator models (e.g., Bouché et al. 2010; Davé et al. 2012; Lilly et al. 2013; Peng & Maiolino 2014; Rathaus & Sternberg 2016), in which galaxy growth is mostly driven by a continuous supply of fresh gas from the cosmic web (Dekel et al. 2009). In these models, gas accretion on halos and subsequently galaxies are controlled by the expansion of the universe. It decreases as $(1+z)^{2-3}$ and scales nearly linearly with halo masses (e.g., Neistein & Dekel 2008). This redshift evolution agrees with the continuous decline of the gas-to-stellar mass ratio of galaxies observed here from $z = 3$ to 0.45. In addition, gas regulator models generally invoke feedback processes—such as stellar winds—to suppress or slow down gas accretion on low-mass halos ($< 10^{11} M_\odot$) (e.g., Bouché et al. 2010; Davé et al. 2012). This introduces a mass-dependency of the gas-to-stellar mass ratio that is consistent with our observations (best-fit exponent of -0.1, though with large uncertainties) and delays the cosmic SFRD peak to $z \sim 2$ because of the relatively low number of halos that are massive enough at early cosmic time

(Bouché et al. 2010). To first order, the rise of the SFRD and $\rho_{\text{gas}}$ from cosmic dawn to $z \sim 2$ would thus be due to the increased number of halos experiencing efficient gas accretion, which is observationally consistent with the increase of the stellar mass function of star-forming galaxies at these epochs. The decrease of the SFRD and $\rho_{\text{gas}}$ at $z < 2$ would then be mostly controlled by the decrease of the gas accretion from the universe's expansion or from shock heating preventing accretion on the most massive halo, which is observationally consistent with the decrease of the gas-to-stellar mass ratio in star-forming galaxies.

A link between the star formation and dust content of galaxies is also expected. Stars produce the metals needed for the formation of dust and at the end of their life-cycle are the locus of significant dust formation, either via an AGB phase ($M \lesssim 8 M_\odot$) or SNe ($M \gtrsim 8 M_\odot$; e.g., Gall et al. 2011, 2014). The coincident peak of the cosmic dust mass density and SFRD suggests that SNe and AGB stars make a very important contribution to dust formation because these formation pathways are linked to star formation on a timescale of 1–2 Gyr (e.g., Dwek et al. 2007; Valiante et al. 2009). However, ISM dust grain growth is supposedly also enhanced in high-redshift star-forming galaxies with high gas densities, and can thus contribute significantly, as well as on a short timescale, to the dust production in these galaxies (Popping et al. 2017). A significant contribution of this latter mechanism to the global dust content of galaxies cannot be ruled out from our observations.

As already noted by Driver et al. (2018), the decrease of $\rho_{\text{dust}}$ at $z < 2$ is at odds with a close-box scenario in which the dust continuously accumulates in galaxies in the absence of destruction or expulsion mechanisms. Our observations suggest instead that at $z < 2$, the dust is destroyed (or expelled) more rapidly than it is formed. Predicting the total amount of dust formed assuming a dust yield of 0.004–0.0065 dust masses for every unit of stellar mass formed as in Driver et al. (2018); the gray-shaded area in Figure 10; i.e., implicitly assuming an redshift-independent initial mass function), we infer that at $z \sim 0$ about 90% of the dust that has been formed in the universe must be destroyed (e.g., astraction, supernovae shocks) or ejected in the IGM (e.g., stellar winds, radiation pressure).

Finally, the decrease of $\rho_{\text{dust}}$ at $z < 2$ is found to be not as pronounced as that of $\rho_{\text{gas}}$ when considering only the molecular phase (i.e., $\rho_{\text{H}_2}(z \sim 0)$ from Saintonge et al. 2017). This flatter evolution of $\rho_{\text{dust}}$ can be explained by the increased contribution of the atomic phase of the ISM and the increased mean metallicity of star-forming galaxies as we approach $z = 0$. Indeed, even though the gas reservoirs of galaxies are replenished by pristine gas from the IGM, their overall metallicity increases from $z \sim 2$ to $z \sim 0$. This is confirmed by the increased zero-point of the stellar mass–metallicity relation, which implies a lower global gas-to-dust mass ratio at $z = 0$ than at $z = 2$.

## 8. Conclusions

We used the deepest ALMA 1.2 mm continuum map to date (rms: 9.5 µJy beam⁻¹) in the HUDF obtained as part of the ASPECS large program to measure the cosmic dust and implied gas ($H_2+H\text{I}$) mass densities as a function of look-back time. We do this by summing (i.e., stacking) the contribution of all the known galaxies in the HUDF above a given stellar
mass in distinct redshift bins; i.e., $\rho_{\text{dust}}(M_{*} > M, z)$ and $\rho_{\text{gas}}(M_{*} > M, z)$. Our galaxy sample is H-band-selected from the available HUDF multi-wavelength catalog, and, up to $z \sim 3$, can be considered as stellar mass-complete down to $\sim 10^{8.9} M_{\odot}$. Dust masses are measured from the 1.2 mm emission of these galaxies assuming a mass-weighted mean dust temperature of $\langle T \rangle_{\text{M}} = 25$ K and a dust emissivity of $\beta = 1.8$ (Scoville et al. 2016). Gas masses are inferred from these dust masses using the local metallicity-dependent gas-to-dust mass ratio and the redshift-dependent stellar mass-metallicity relations (as in Tacconi et al. 2018). With this unique dataset and approach, we find the following:

1. $\rho_{\text{dust}}(M_{*} > M, z)$ and $\rho_{\text{gas}}(M_{*} > M, z)$ grow rapidly as $M$ decreases down to $10^{10} M_{\odot}$, but this growth slows down as $M$ decreases to even lower stellar masses. This flattening implies that at our stellar mass-completeness limits, $\rho_{\text{dust}}(M_{*} > M, z)$ and $\rho_{\text{gas}}(M_{*} > M, z)$ converge already toward the total cosmic dust ($\rho_{\text{dust}}$) and gas ($\rho_{\text{gas}}$) mass densities in galaxies; i.e., with only a minor contribution by galaxies below our stellar mass-completeness limits.

2. The contribution of quiescent galaxies (i.e., galaxies with little on-going star formation and selected using the standard UVJ criterion) to $\rho_{\text{dust}}$ and $\rho_{\text{gas}}$ is negligible ($\lesssim 10\%$). The bulk of the dust and gas in galaxies appears to be locked in star-forming galaxies with $M_{*} \gtrsim 10^{9.5} M_{\odot}$.

3. The gas ($H_2+H\text{I}$) mass densities measured here agree with the molecular gas mass densities inferred from the CO observations of the ASPECS (Decarli et al. 2016, 2019) and COLDz (Riechers et al. 2019) surveys. In the redshift range probed here ($z = 0.45-3.0$), dust-based measurements are thus dominated by the dust in the molecular phase. This suggests that either the bulk of the gas in galaxies is in a molecular phase or that there is a significant enhancement of the dust emissivity in the dense/molecular ISM with respect to that in the more diffuse/atomic ISM.

4. The cosmic dust (gas) mass density increases at early cosmic time, peaks around $z = 1-3$, and decreases by a factor $\sim 4$ (7) until the present time (combining our results with low-redshift measurements; Dunne et al. 2011; Decarli et al. 2016, 2019; Saintonge et al. 2017). The redshift evolution of $\rho_{\text{gas}}$ matches that of the cosmic SFRD, while the decline of $\rho_{\text{dust}}$ at $z < 2$ is less pronounced than that observed for the SFRD.

5. The evolution of $\rho_{\text{dust}}(M_{*} > M, z)$ [$\rho_{\text{gas}}(M_{*} > M, z)$] with stellar masses and redshifts is accurately modeled using the stellar mass function of star-forming galaxies and their average dust[gas]-to-stellar mass ratio. The dust [gas] content of galaxies per unit stellar mass continuously decreases from $z = 3$ to $z = 0$, while having a mild dependency on stellar masses—with a best-fit power-law stellar mass-dependent exponent of 0.1 [$-0.1$].

Our results suggest that galaxies have a relatively constant star formation efficiency (SFR/$M_{\text{gas}}$) across cosmic time

(within a factor $\sim 2$; solid line in Figure 9). Their star formation seems to be mainly controlled by the supply of fresh gas from the cosmic web (Dekel et al. 2009), as advocated by the gas regulator models (e.g., Bouché et al. 2010; Davé et al. 2012; Lilly et al. 2013; Peng & Maiolino 2014; Rathaus & Sternberg 2016). This supply, which varies with cosmic time as $(1+z)^{2.3}$ following the universe's expansion, is in turn the main driver of the continuous decrease of the SFRD at $z < 2$. The decrease of $\rho_{\text{dust}}$ at $z < 2$ implies that a large fraction ($\sim 90\%$; see Driver et al. 2018) of the dust formed in galaxies is, within few Gyr, destroyed (shock, astraction) or ejected to the IGM (wind, radiation pressure).

We would like to thank the referee for their comments that have helped to improve our paper. This research was carried out within the Collaborative Research Centre 956, sub-project A1, funded by the Deutsche Forschungsgemeinschaft (DFG)—project ID 184018867. B.M. thanks Aswin P. Vijayan for providing predictions from the L-GALAXIES semi-analytic model. M.N. and F.W. acknowledge support from ERC Advanced Grant 740246 (Cosmic_Gas). D.R. acknowledges support from the National Science Foundation under grant No. AST-1614213 and from the Alexander von Humboldt Foundation through a Humboldt Research Fellowship for Experienced Researchers. I.R.S. acknowledges support from STFC (ST/P000541/1). R.J.A. was supported by FONDECYT grant No. 1191124. Este trabajo contó con el apoyo de CONICYT + PCI + INSTITUTO MAX PLANCK DE ASTRONOMIA MPG190030. Este trabajo contó con el apoyo de CONICYT + Programa de Astronomía+ Fondo CHINA-CONICYT CAS16026. This paper makes use of the ALMA data ADS/JAO.ALMA#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. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

Facility: ALMA.

## Appendix A
### Different Gas-to-dust Mass Conversions

In Figure 11, we use a $\Delta_{\text{GDR}}$-metallicity relation with a much steeper power-law ($\sim 3$; Rémy-Ruyer et al. 2014) at metallicity $M, z$) with $M$, large disagreements with the ASPECS-CO survey's results, that could suggest an increasing contribution of the atomic phase to $\rho_{\text{gas}}(M_{*} > M, z)$ at low metallicities.

In Figure 12, we infer the redshift evolution of $\rho_{\text{gas}}$ scaling $\rho_{\text{dust}}$ by a constant gas-to-dust mass ratio of 100, typical for massive galaxies at $z \sim 0$. These $\rho_{\text{gas}}$ values are lower than our original measurements by 30% and underestimate at the $1-2\sigma$ level those from the ASPECS-CO LP at $z = 1-2$.
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Same as Figure 7 but here the gas mass densities are inferred using a gas-to-dust mass ratio–metallicity relation with a much steeper power-law at metallicity <7.9 (~3; Rémy-Ruyer et al. 2014).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 12.** Same as Figure 9 but here the cosmic gas mass densities are inferred applying a constant gas-to-dust mass ratio of 100 (typical for massive galaxies at $z \sim 0$) to the cosmic dust mass densities presented in Figure 6.</figcaption></figure>
## Appendix B
### Comparison to Other Simulations

In this appendix, we compare our measurements to predictions from the semi-analytical model of Figure 13 (Vijayan et al. 2019).

<figure><img src="image_14.png" /><figcaption>**Figure 13.** Same as Figure 10 but here we compare our measurements to predictions from the semi-analytical model of Vijayan et al. (2019). For $M_{*} > 10^{9} M_{\odot}$, we used their “Millennium” predictions, while for $M_{*} > 10^{8} M_{\odot}$ we use their “Millennium-II” predictions.</figcaption></figure>
ORCID iDs

Benjamin Magnelli https://orcid.org/0000-0002-6777-6490

Leindert Boogaard https://orcid.org/0000-0002-3952-8588

Roberto Decarli https://orcid.org/0000-0002-2662-8803

Jorge Gónzalez-López https://orcid.org/0000-0003-3926-1411

Mladen Novak https://orcid.org/0000-0001-8695-825X

Gergö Popping https://orcid.org/0000-0003-1151-4659

Ian Smail https://orcid.org/0000-0003-3037-257X

Fabian Walter https://orcid.org/0000-0003-4793-7880

Manuel Aravena https://orcid.org/0000-0002-6290-3198

Roberto J. Assef https://orcid.org/0000-0002-9508-3667

Franz Erik Bauer https://orcid.org/0000-0002-8686-8737

Frank Bertoldi https://orcid.org/0000-0002-1707-1775

Chris Carilli https://orcid.org/0000-0001-6647-3861

Paulo C. Cortes https://orcid.org/0000-0002-3583-780X

Elisabete da Cunha https://orcid.org/0000-0001-9759-4797

Emanuele Daddi https://orcid.org/0000-0002-3331-9590

Tanio Díaz-Santos https://orcid.org/0000-0003-0699-6083

Robert J. Ivison https://orcid.org/0000-0001-5118-1313

Olivier Le Fèvre https://orcid.org/0000-0001-5891-2596

Pascal Oesch https://orcid.org/0000-0001-5851-6649

Dominik Riechers https://orcid.org/0000-0001-9585-1462

Hans-Walter Rix https://orcid.org/0000-0003-4996-9069

Mark T. Sargent https://orcid.org/0000-0003-1033-9684

Paul van der Werf https://orcid.org/0000-0001-5434-5942

Axel Weiss https://orcid.org/0000-0003-4678-3939

References

Aoyama, S., Hou, K.-C., Hirashita, H., Nagamine, K., & Shimizu, I. 2018, MNRAS, 478, 4905

Aravena, M., Decarli, R., Walter, F., et al. 2016, ApJ, 833, 68

Auld, R., Bianchi, S., Smith, M. W. L., et al. 2013, MNRAS, 428, 1880

Bacon, R., Conseil, S., Mary, D., et al. 2017, A&A, 608, A1

Beckwith, S. V. W., Stiavelli, M., Koekemoer, A. M., et al. 2006, AJ, 132, 1729

Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013, ARA&A, 51, 207

Boogaard, L. A., Decarli, R., González-López, J., et al. 2019, ApJ, 882, 140

Bothwell, M. S., Smail, I., Chapman, S. C., et al. 2013, MNRAS, 429, 3047

Bouché, N., Dekel, A., Genzel, R., et al. 2010, ApJ, 718, 1001

Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2011, ApJ, 737, 90

Chabrier, G. 2003, PASP, 115, 763

Chapin, E. L., Pope, A., Scott, D., et al. 2009, MNRAS, 398, 1793

Cibinel, A., Daddi, E., Bournaud, F., et al. 2017, MNRAS, 469, 4683

Coe, D., Benítez, N., Sánchez, S. F., et al. 2006, AJ, 132, 926

Coogan, R. T., Sargent, M. T., Daddi, E., et al. 2019, MNRAS, 485, 2092

da Cunha, E., Charlot, S., & Elbaz, D. 2008, MNRAS, 388, 1595

da Cunha, E., Groves, B., Walter, F., et al. 2013, ApJ, 766, 13

da Cunha, E., Walter, F., Smail, I. R., et al. 2015, ApJ, 806, 110

Dale, D. A., Aniano, G., Engelbracht, C. W., et al. 2012, ApJ, 745, 95

Davé, R., Anglés-Alcázar, D., Narayanan, D., et al. 2019, MNRAS, 486, 2827

Davé, R., Finlator, K., & Oppenheimer, B. D. 2012, MNRAS, 421, 98

Davidzon, I., Ilbert, O., Laigle, C., et al. 2017, A&A, 605, A70

De Bernardis, F., & Cooray, A. 2012, ApJ, 760, 14

Decarli, R., Walter, F., Aravena, M., et al. 2016, ApJ, 833, 69

Decarli, R., Walter, F., Aravena, M., et al. 2016, ApJ, 833, 70

Decarli, R., Walter, F., González-López, J., et al. 2019, ApJ, 882, 138

Dekel, A., Birnboim, Y., Engel, G., et al. 2009, Natur, 457, 451

Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium (Princeton, NJ: Princeton Univ. Press)

Draine, B. T., & Li, A. 2007, ApJ, 657, 810

Driver, S. P., Andrews, S. K., da Cunha, E., et al. 2018, MNRAS, 475, 2891

Driver, S. P., & Robotham, A. S. G. 2010, MNRAS, 407, 2131

Dunlop, J. S., McLure, R. J., Biggs, A. D., et al. 2017, MNRAS, 466, 861

Dunne, L., Eales, S. A., & Edmunds, M. G. 2003, MNRAS, 341, 589

Dunne, L., Gomez, H. L., da Cunha, E., et al. 2011, MNRAS, 417, 1510

Dwek, E., Galliano, F., & Jones, A. P. 2007, ApJ, 662, 927

Eales, S., Smith, M. W. L., Auld, R., et al. 2012, ApJ, 761, 168

Elbaz, D., Leiton, R., Nagar, N., et al. 2018, A&A, 616, A110

Franco, M., Elbaz, D., Bethermin, M., et al. 2018, A&A, 620, A152

Gall, C., Hjorth, J., & Andersen, A. C. 2011, A&ARV, 19, 43

Gall, C., Hjorth, J., Watson, D., et al. 2014, Natur, 511, 326

Genzel, R., Tacconi, L. J., Lutz, D., et al. 2015, ApJ, 800, 20

Gobat, R., Daddi, E., Magdis, G., et al. 2018, NatAs, 2, 239

González-López, J., Novak, M., Decarli, R., et al. 2020, ApJ, in press (arXiv:2002.07199)

Grogin, N. A., Kocevski, D. D., Faber, S. M., et al. 2011, ApJS, 197, 35

Illingworth, G. D., Magee, D., Oesch, P. A., et al. 2013, ApJS, 209, 6

Inami, H., Bacon, R., Brinchmann, J., et al. 2017, A&A, 608, A2

Kaasinen, M., Scoville, N. Z., Walter, F., et al. 2019, ApJ, 880, 15

Karim, A., Schinnerer, E., Martinez-Sansigre, A., et al. 2011, ApJ, 730, 61

Koekemoer, A. M., Ellis, R. S., McLure, R. J., et al. 2013, ApJS, 209, 3

Koekemoer, A. M., Faber, S. M., Ferguson, H. C., et al. 2011, ApJS, 197, 36

Kovács, A., Omont, A., Beelen, A., et al. 2010, ApJ, 717, 29

Labbé, I., Oesch, P. A., Illingworth, G. D., et al. 2015, ApJS, 221, 23

Le Fèvre, O., Vettolani, G., Garilli, B., et al. 2005, A&A, 439, 845

Leroy, A. K., Bolatto, A., Gordon, K., et al. 2011, ApJ, 737, 12

Li, A., & Draine, B. T. 2001, ApJ, 554, 778

Li, Q., Narayanan, D., & Davé, R. 2019, MNRAS, 490, 1425

Lilly, S. J., Carollo, C. M., Pipino, A., Renzini, A., & Peng, Y. 2013, ApJ, 772, 119

Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415

Magdis, G. E., Daddi, E., Bethermin, M., et al. 2012, ApJ, 760, 6

Magnelli, B., Karim, A., Staguhn, J., et al. 2019, ApJ, 877, 45

Magnelli, B., Lutz, D., Saintonge, A., et al. 2014, A&A, 561, 86

Magnelli, B., Lutz, D., Santini, P., et al. 2012a, A&A, 539, A155

Magnelli, B., Saintonge, A., Lutz, D., et al. 2012b, A&A, 548, A22

Martin, A. M., Papastergis, E., Giovanelli, R., et al. 2010, ApJ, 723, 1359

McKinnon, R., Torrey, P., Vogelsberger, M., Hayward, C. C., & Marinacci, F. 2017, MNRAS, 468, 1505

McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in ASP Conf. Ser. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell (San Francisco, CA: ASP), 127

Ménard, B., & Fukugita, M. 2012, ApJ, 754, 116

Momcheva, I. G., Brammer, G. B., van Dokkum, P. G., et al. 2016, ApJS, 225, 27

Morris, A. M., Kocevski, D. D., Trump, J. R., et al. 2015, AJ, 149, 178

Mortlock, A., Conselice, C. J., Hartley, W. G., et al. 2015, MNRAS, 447, 2

Neistein, E., & Dekel, A. 2008, MNRAS, 388, 1792

Peng, Y.-j., & Maiolino, R. 2014, MNRAS, 443, 3643

Pettini, M., & Pagel, B. E. J. 2004, MNRAS, 348, L59

Planck Collaboration, Abergel, A., Ade, P. A. R., et al. 2011, A&A, 536, A21

Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594 A13

Popping, G., Somerville, R. S., & Galametz, M. 2017, MNRAS, 471, 3152

Pozzi, F., Calura, F., Zamorani, G., et al. 2019, MNRAS, 491, 5073

Rathaus, B., & Sternberg, A. 2016, MNRAS, 458, 3168

Rémy-Ruyer, A., Madden, S. C., Galliano, F., et al. 2014, A&A, 563, A31

Riechers, D. A., Pavesi, R., Sharon, C. E., et al. 2019, ApJS, 872, 7

Rujopakarn, W., Dunlop, J. S., Rieke, G. H., et al. 2016, ApJ, 833, 12

Saintonge, A., Catinella, B., Tacconi, L. J., et al. 2017, ApJS, 233, 22

Saintonge, A., Lutz, D., Genzel, R., et al. 2013, ApJ, 778, 2

Santini, P., Maiolino, R., Magnelli, B., et al. 2014, A&A, 562, 30

Sargent, M. T., Daddi, E., Bournaud, F., et al. 2015, ApJL, 806, L20

Schinnerer, E., Groves, B., Sargent, M. T., et al. 2016, ApJ, 833, 112

Schreiber, C., Pannella, M., Elbaz, D., et al. 2015, A&A, 575, A74

Scoville, N., Aussel, H., Sheth, K., et al. 2014, ApJ, 783, 84

Scoville, N., Lee, N., Vanden Bout, P., et al. 2017, ApJ, 837, 150

Scoville, N., Sheth, K., Aussel, H., et al. 2016, ApJ, 820, 83

Skelton, R. E., Whitaker, K. E., Momcheva, I. G., et al. 2014, ApJS, 214, 24

Tacconi, L. J., Genzel, R., Saintonge, A., et al. 2018, ApJ, 853, 179

Thacker, C., Cooray, A., Smidt, J., et al. 2013, ApJ, 768, 58

Valiante, R., Schneider, R., Bianchi, S., & Andersen, A. C. 2009, MNRAS, 397, 1661

Vijayan, A. P., Clay, S. J., Thomas, P. A., et al. 2019, MNRAS, 489, 4072

Walter, F., Decarli, R., Aravena, M., et al. 2016, ApJ, 833, 67

Whitaker, K. E., Franx, M., Leja, J., et al. 2014, ApJ, 795, 104

Xu, C., Pirzkal, N., Malhotra, S., et al. 2007, AJ, 134, 169