# Cosmic evolution of the H2 mass density and the epoch of molecular gas

T. K. Garratt<EMAIL_ADDRESS>Centre for Astrophysics Research,
University of Hertfordshire, Hatfield, AL10 9AB, UK.  K. E. K. Coppin  Centre
for Astrophysics Research, University of Hertfordshire, Hatfield, AL10 9AB,
UK.  J. E. Geach  Centre for Astrophysics Research, University of
Hertfordshire, Hatfield, AL10 9AB, UK.  O. Almaini  School of Physics and
Astronomy, University of Nottingham, University Park, Nottingham, NG7 2RD, UK
W. G. Hartley  Department of Astronomy, University of Geneva, ch. d’Écogia 16,
CH-1290 Versoix, Switzerland  D. T. Maltby  School of Physics and Astronomy,
University of Nottingham, University Park, Nottingham, NG7 2RD, UK  C. J.
Simpson  Gemini Observatory, Hilo, HI 96720, USA  A. Wilkinson  Sterrenkundig
Observatorium, Universiteit Gent, Krijgsl aan 281 S9, 9000 Gent, Belgium  C.
J. Conselice  School of Physics and Astronomy, University of Nottingham,
University Park, Nottingham, NG7 2RD, UK  M. Franco  Centre for Astrophysics
Research, University of Hertfordshire, Hatfield, AL10 9AB, UK.  R. J. Ivison
European Southern Observatory, Karl Schwarzschild Strasse 2, Garching, Germany
M. P. Koprowski  Institute of Astronomy, Faculty of Physics, Astronomy and
Informatics, Nicolaus Copernicus University, Grudziadzka 5, 87-100 Torun,
Poland  C. C. Lovell  Centre for Astrophysics Research, University of
Hertfordshire, Hatfield, AL10 9AB, UK.  A. Pope  Department of Astronomy,
University of Massachusetts, 710 North Pleasant Street Amherst, MA 01003, USA
D. Scott  Dept. of Physics & Astronomy, University of British Columbia,
Vancouver, Canada  P. van der Werf  Leiden Observatory, Leiden University,
P.O. Box 9513, NL-2300 RA Leiden, The Netherlands

###### Abstract

We present new empirical constraints on the evolution of $\rho_{\rm H_{2}}$,
the cosmological mass density of molecular hydrogen, back to $z\approx 2.5$.
We employ a statistical approach measuring the average observed $850\,\mu$m
flux density of near-infrared selected galaxies as a function of redshift. The
redshift range considered corresponds to a span where the $850\,\mu$m band
probes the Rayleigh-Jeans tail of thermal dust emission in the rest-frame, and
can therefore be used as an estimate of the mass of the interstellar medium
(ISM). Our sample comprises of ${\approx}150,000$ galaxies in the UKIDSS-UDS
field with near-infrared magnitudes $K_{\rm AB}\leq 25$ mag and photometric
redshifts with corresponding probability distribution functions derived from
deep $12$-band photometry. With a sample approximately $2$ orders of magnitude
larger than in previous works we significantly reduce statistical
uncertainties on $\rho_{\rm H_{2}}$ to $z\approx 2.5$. Our measurements are in
broad agreement with recent direct estimates from blank field molecular gas
surveys, finding that the epoch of molecular gas coincides with the peak epoch
of star formation with $\rho_{\rm H_{2}}\approx 2\times 10^{7}\,{\rm
M_{\odot}}\,{\rm Mpc^{-3}}$ at $z\approx 2$. We demonstrate that $\rho_{\rm
H_{2}}$ can be broadly modelled by inverting the star-formation rate density
with a fixed or weakly evolving star-formation efficiency. This “constant
efficiency” model shows a similar evolution to our statistically derived
$\rho_{\rm H_{2}}$, indicating that the dominant factor driving the peak star
formation history at $z\approx 2$ is a larger supply of molecular gas in
galaxies rather than a significant evolution of the star-formation rate
efficiency within individual galaxies.

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

##  1 Introduction

Three intimately linked observational tracers broadly characterise the cosmic
evolution of galaxies: the volume averaged star-formation rate density
$\rho_{\dot{M}_{\star}}(z)$, the stellar mass density $\rho_{M_{\star}}(z)$,
and the molecular gas density $\rho_{M_{\rm H_{2}}}(z)$. Our current
understanding of galaxy evolution is largely driven by comprehensive
measurements of the former two (see Madau & Dickinson, 2014, for a review),
with a clear empirical picture emerging of an evolution of star formation,
which rises rapidly to a peak around $z\approx 2$ and then decays to the
present day. Completing the triptych is important since the evolution of the
molecular gas content of galaxies encodes several important pieces of
astrophysics: gas consumption in star formation; gas recycling via feedback;
and fresh gas accretion. Ultimately, it is the evolution of molecular gas that
drives galaxy evolution as it is the fuel from which stars are assembled.
Measurements of molecular gas in galaxies are therefore needed to complete the
picture, and to resolve a key outstanding question: Was the peak of star
formation history driven by a larger supply of molecular gas or because
galaxies formed stars more efficiently (e.g., driven by galaxy
mergers/instabilities etc.), or both?

The bulk of the cold gas reservoir in the Universe is comprised of hydrogen
gas in the form of atomic hydrogen (HI) and molecular hydrogen (H2). In the
current model of galaxy formation gas is delivered into galaxies via hot- or
cold-mode accretion (e.g., Birnboim & Dekel, 2003). The cooling gas must form
H2 for star formation to occur. The two main routes to H2 formation in
galaxies are via the gas phase reaction ${\rm H+e^{-}\rightarrow
H^{-}+\gamma}$, ${\rm H^{-}+H\rightarrow H_{2}+e^{-}}$, and via a dust phase,
where H2 forms on the surface of dust grains via efficient three-body
reactions (Gould & Salpeter, 1963).

H2 radiates poorly in typical ISM conditions due to the lack of a permanent
dipole moment and a minimum rotational excitation temperature that is
significantly higher (${\approx}500$ K) than typical temperatures of the cold
star-forming ISM (Wakelam et al., 2017). However, H2 can be indirectly traced
through its interactions with CO, which traces the same cold, dense ISM and
has a low dipole moment enabling its excitation in regions of low density
($n_{\rm crit}\sim 10^{2}\,{\rm cm}^{-3}$). Consequently, 12CO, the second
most abundant molecule in the ISM, is commonly used as a tracer of the
available reservoir of molecular gas in galaxies (e.g., Solomon & Vanden Bout,
2005; Carilli & Walter, 2013). The ground state transition CO ($J=1\rightarrow
0$) is a reliable tracer of total molecular gas, with the conversion factor
from CO luminosity to H2 gas mass ($\alpha_{\rm
CO}=M_{\rm{H}_{2}}/L^{\prime}_{\rm CO}$) calibrated locally (see Bolatto et
al., 2013, for a review). Observing the ground-state transition line avoids
additional uncertainties inherent in observations of higher-J CO transitions,
which require a correction for gas excitation to derive the equivalent CO
($J=1\rightarrow 0$) luminosity.

Until recently measurements of the cosmological molecular gas mass density
were hampered by a paucity of observational data. Over the past few years
direct measurements of the cold molecular gas reservoirs of individual
galaxies have increased rapidly, with surveys primarily targeting star-forming
and lensed galaxies (e.g., Frayer et al., 1998, 1999; Coppin et al., 2007;
Tacconi et al., 2010; Ivison et al., 2011; Thomson et al., 2012; Bothwell et
al., 2013; Riechers et al., 2013; Stach et al., 2017; Oteo et al., 2018;
Gómez-Guijarro et al., 2019; Lenkić et al., 2019). However, as these surveys
rely on observationally-expensive detections of faint spectral lines,
measurements of molecular gas mass are still dwarfed in number in comparison
to the samples for which star-formation rate (SFR) and stellar mass estimates
are derived. Moreover, to properly assess the cosmological evolution of the
cold gas content of galaxies a blank field survey approach is required to
measure the gas mass function, rather than targeted (and therefore biased)
observations of high-z galaxies as has generally been the case for cold gas
observations outside the local volume.

Recently surveys using a blank field molecular line scan strategy have emerged
as an alternative to targeted observations. These surveys evade many of the
biases towards massive star-forming galaxies inherent in targeted approaches.
The inaugural blank field CO survey employed the Plateau de Bure
Interferometer in observations of the Hubble Deep Field North (Walter et al.,
2014; Decarli et al., 2014). This was followed more recently by the ALMA
Spectroscopic Survey in the Hubble Ultra Deep Field (ASPECS; Walter et al.,
2016; Decarli et al., 2016a, b; Aravena et al., 2016a, b; Bouwens et al.,
2016; Carilli et al., 2016), the ASPECS Large Program (ASPECS LP; Decarli et
al., 2019; González-López et al., 2019; Boogaard et al., 2019; Popping et al.,
2019; Decarli et al., 2020) and the CO Luminosity Density at High Redshift
survey (COLDZ; Pavesi et al., 2018; Riechers et al., 2019). These surveys have
presented results setting out valuable new constraints on the evolution of
$\rho_{\rm H_{2}}$ over a redshift range $0\lesssim z\lesssim 7$, obtained
through blank field observations of CO line emission. However, due to low
number statistics and the small survey areas (which are prone to strong
clustering-enhanced sample variance) used to derive the CO luminosity
functions, these measurements are hampered by large statistical uncertainties.

To combat the shortfall in direct measurements of molecular gas Scoville
(2013); Scoville et al. (2014, 2016, 2017) employed a complementary approach
that utilises submillimeter observations of the long wavelength dust continuum
as a measure of the molecular gas mass in galaxies. The Rayleigh-Jeans (RJ)
tail is nearly always optically thin, and consequently measurements of dust
emission can be used as a direct probe of molecular gas mass (e.g., Eales et
al., 2012; Magdis et al., 2012). Whilst ordinarily a conversion from dust to
gas mass would require dust emissivity and dust-to-gas abundance to be
constrained, Scoville et al. (2014, 2016, 2017) circumvent this by deriving an
empirically calibrated RJ luminosity-to-gas mass ratio using CO
($J=1\rightarrow 0$) and submillimeter continuum observations of a sample of
normal star-forming and starburst galaxies at low-$z$ and submillimeter
galaxies (SMGs) at high-$z$. This approach requires assumptions about dust
temperature and the evolution of the gas-to-dust mass ratio, but provides
molecular gas mass ($M_{\rm mol}$) estimates within factor of ${\approx}2$
accuracy (e.g., Scoville et al., 2016; Kaasinen et al., 2019). Since dust
continuum measurements can be made in minutes (in contrast to CO line
observations which can take multiple hours, e.g., Bothwell et al., 2013;
Tacconi et al., 2013) this method can be used to derive molecular gas
measurements for much larger samples of galaxies.

The Scoville et al. (2016) RJ luminosity-to-gas mass calibration has been used
to estimate the molecular gas mass for $\sim 700$ ALMA-detected galaxies from
the COSMOS field, with Scoville et al. (2017) deriving molecular gas masses
for individual galaxies at redshifts $0.3<z<4.5$ and Liu et al. (2019)
extending this approach to redshifts of $z\approx 6$. This method has also
been used in combination with stacking methodologies to estimate average
molecular gas masses for large samples of galaxies. Millard et al. (2020) use
the Scoville et al. (2016) calibration and apply this to a sample of $63,658$
galaxies to derive the gas mass fraction out to $z\approx 5$. Magnelli et al.
(2020) use a method similar to Scoville et al. (2016, and which is equivalent
at solar metallicity) and apply this to a sample of $555$ galaxies to derive
the molecular gas mass density to $z\approx 3$.

In this paper we contribute to the picture of cosmic galaxy evolution by
building on the approach of Scoville (2013); Scoville et al. (2014, 2016,
2017) estimating the evolution of the cosmological mass density of molecular
hydrogen to $z\approx 2.5$ via the average submillimeter continuum emission of
a sample of $150{,}000$ galaxies selected from a deep near-infrared survey in
the well studied UKIDSS-UDS field. The Ultra-Deep Survey (UDS) is the deepest
component of the UK InfraRed Telescope (UKIRT) Infrared Deep Sky Survey
(UKIDSS; Lawrence et al., 2007). We limit our estimate of the molecular gas
mass density to $z\approx 2.5$ as the Scoville et al. (2016) calibration has
only been shown to be robust out to this redshift. Adopting a statistical
approach allows us to take advantage of a near-infrared selected sample which
is an order of magnitude larger than in surveys that measure dust emission
(e.g., Scoville et al., 2017; Liu et al., 2019; Magnelli et al., 2020) or CO
spectral line emission (e.g., Walter et al., 2016; Riechers et al., 2019;
Decarli et al., 2019; Kaasinen et al., 2019) for individual sources. Our
method differs from previous stacking approaches (e.g., Millard et al., 2020;
Magnelli et al., 2020) as we do not use a combination of spectroscopic and
photometric redshifts for our binning. Instead, in the absence of a sample
complete with spectroscopic redshifts, we utilise the full photometric
redshift probability distribution functions for all our sources. Our method is
complementary to previous works in this field (e.g., Decarli et al., 2020; Liu
et al., 2019; Magnelli et al., 2020; Riechers et al., 2019) and allows us to
reduce the statistical uncertainties on the cosmological molecular gas mass
density out to $z\approx 2.5$.

This paper is organized as follows: in Section $2$, we define the maps and
catalogs used; in Section $3$, we present a $3$-dimensional stacking method
which we employ to measure the average (stacked) observed $850\,\mu$m flux
densities for near-infrared selected galaxies as a function of redshift; in
Section $4$, we show that the approach of Scoville et al. (2016) can be
applied to our stacked $850\,\mu$m flux densities to derive the cosmological
molecular gas density to $z\approx 2.5$. We also demonstrate that the cosmic
molecular gas density can be broadly modelled by $2$ complementary approaches
(i) from the halo mass function assuming a constant halo mass range, and
employing stellar-halo mass and ISM-stellar mass ratios, and (ii) inverting
the star-formation rate density assuming a “constant efficiency” model, and in
Section $5$ we interpret the overall evolution of the cosmic molecular gas
mass density in the context of our results and in comparison to previous
works. We present our conclusions in Section $6$. We assume a Planck $2015$
cosmology, where $\Omega_{\rm m}=0.31,\,\Omega_{\Lambda}=0.69,\,H_{0}=68\,{\rm
km\,s^{-1}\,Mpc^{-1}}$ (Planck Collaboration et al., 2016), and a Chabrier
(2003) Initial Mass Function. The AB magnitude system is used throughout.

##  2 Data

###  2.1 SCUBA-2 Cosmology Legacy Survey

The UKIDSS-UDS field was mapped at $850\,\mu$m as part of the Submillimetre
Common-User Bolometer Array 2 (SCUBA-2) Cosmology Legacy Survey (Geach et al.,
2017). The full details of the data collection, reduction and map properties
are given in Geach et al. (2017). Briefly, the beam-convolved map spans
approximately $1\,{\rm{deg}}^{2}$ covering the bulk of the multi-wavelength
coverage of this field, with a uniform (instrumental) noise of
$\sigma_{850}=0.9$ mJy beam-1. Geach et al. (2017) estimate the SCUBA-2
confusion limit to be $\sigma_{\rm conf}=0.8$ mJy beam-1. The beam full width
half maximum (FWHM) is approximately 15′′, with a full analytic description of
the point spread function (PSF) given by Geach et al. (2017).

###  2.2 UKIDSS-UDS ultraviolet–optical–mid-infrared imaging and catalog

The UDS Data Release 11 (DR$11$) $12$-band matched catalog is K-band selected
with the $95\%$ completeness limit estimated to be $K_{\rm AB}=25$ mag. The
full details of this catalog will be comprehensively provided in Almaini et
al., (in prep.) and Hartley et al., (in prep.), and only a summary is given
here. The catalog provides photometry in 12 bands (U, B, V, R, $i^{\prime}$,
$z^{\prime}$, Y, J, H, K, $3.6\,\mu$m and $4.5\,\mu$m), where available.

The J, H, and K photometry is taken from the DR11 release of UKIDSS-UDS. The
UKIDSS project, described in Lawrence et al. (2007) utilises the UKIRT Wide
Field CAMera (WFCAM; Casali et al., 2007). The photometric system and
calibration are outlined in Hewett et al. (2006) and Hodgkin et al. (2009),
respectively, and the pipeline processing and science archive are described in
Irwin et al., (in prep.) and Hambly et al. (2008). UKIDSS-UDS covers an area
of $0.8\,{\rm deg}^{2}$, reaching median depths of J$=25.6$, H$=25.1$, and
K$=25.3$ ($5\sigma$, AB, estimated from $2^{\prime\prime}$ apertures in source
free areas; Almaini et al. in prep.).

The B, V, R, $i^{\prime}$, and $z^{\prime}$ optical imaging is from the
Subaru/XMM-Newton Deep Survey, which utilises Suprime-Cam on the Subaru
Telescope (Furusawa et al., 2008). U-band data are from the Canada-France-
Hawaii Telescope Megacam instrument (Almaini et al. in prep.) and Y-band
imaging is obtained from the VISTA Deep Extragalactic Observations survey
(VISTA-VIDEO: Jarvis et al., 2013). The InfraRed Array Camera (IRAC) imaging
at $3.6\,\mu$m and $4.5\,\mu$m is from the Spitzer UKIDSS Ultra Deep Survey
(SpUDS: PI Dunlop), combined with deeper data from the Spitzer Extended Deep
Survey (SEDS: Ashby et al., 2013). To expand the coverage to outer regions of
the field, shallower data are also used from the SIRTF Wide-area InfraRed
Extragalactic survey (SWIRE: Lonsdale et al., 2003).

UDS DR11 provides image masks, with masked regions corresponding to image
boundaries, artefacts, and bright stars. We employ the UDS binary mask for
“good” regions which has an unmasked area of $0.64\,{\rm deg}^{2}$. This
binary mask combines the masked regions of the photometry images detailed
above (not including the deeper SEDS or SpUDS IRAC images).

We also utilise the subsets feature of the UKIDSS-UDS catalog and for our
galaxy sample chose the catalog-defined “good galaxy” subset, which comprises
$217{,}429$ sources. These sources have full $12$-band photometry and lie
within the corresponding “good” mask regions, are not cross-talk sources (for
which JHK photometry is likely compromised), and are not classified as stars.

UDS DR11 also includes photometric redshifts derived using the code eazy (Easy
and Accurate $Z_{\rm phot}$ from Yale; Brammer et al., 2008). To estimate the
photometric redshift for each source the 12-band broadband photometry was fit
with a spectral energy distribution template producing a redshift probability
distribution (Hartley et al. in prep.). We utilise both the maximum-likelihood
photometric redshifts and redshift probability distributions provided with UDS
DR11. eazy performs well compared to other commonly used photometric redshift
codes (e.g., zphot, HyperZ, Rainbow; Dahlen et al., 2013), with the resulting
normalized mean absolute deviation between eazy derived photometric redhifts
and spectroscopic redshifts found to be only $\sigma_{\rm nmad}\approx 0.02$
(Hartley et al. in prep.).

##  3 Methods

We employ a $3$-dimensional stacking approach based on the simultaneous
stacking algorithm simstack (presented in detail in Viero et al., 2013). This
method allows for the simultaneous fitting of the average observed flux
density for multiple populations that contribute to the flux density in the
observed map (such as a population of galaxies split into bins of redshift).
Importantly, this method takes into account the (usually) large beam in
single-dish submillimeter maps with simulations demonstrating that this method
returns an unbiased estimate of the average observed flux density for beam
sizes ranging from ${\rm FWHM}=15$–$35^{\prime\prime}$ (Viero et al., 2013).
This approach also mitigates against boosting of stacking signals from
clustered galaxies (e.g., Chary & Pope, 2010; Alberts et al., 2014).

Our goal is to find the average observed submillimeter flux densities at given
redshift intervals for a population of near-infrared selected galaxies that
best fit the observed flux density in the SCUBA-2 map, taking into account the
convolution of point sources with the large beam. In this work, rather than
binning galaxies by discrete photometric redshift values, we split our sample
across redshift intervals according to the redshift probability distribution
of each source.

First, we define our sample, performing a selection in observed $K$-band total
magnitude, $K_{\rm AB}\leq 25$ mag, with the faint-end corresponding to the
$95\%$ completeness limit of the UKIDSS-UDS catalog, giving us a sample of
153,399 galaxies. At this limiting magnitude the $95\%$ stellar mass
completeness is $\approx 10^{9.5}\,{\rm M}_{\odot}$ at $z=2.5$ (Wilkinson et
al. in prep.). The redshift probability distribution, $\mathcal{P}(z)$, for
each source is discretized in bins of $\Delta z$ (Hartley et al. in prep.). We
make a completeness correction to the redshift probability distribution of
each source, such that $\mathcal{P}(z)$ of a source of magnitude $K$
integrates to $C(K)^{-1}$, where $C(K)$ is the catalog completeness at $K$
(Hartley et al. in prep.). We assume there is no systematic redshift bias in
$C(K)$ for this correction.

With the sample defined we consider a sky model in which each galaxy
contributes a flux density that can be described as

$\displaystyle S_{\nu}=\int_{0}^{\infty}\mathcal{S}_{\nu}(z)\mathcal{P}(z){\rm
d}z.$ (1)

where $\mathcal{P}(z)$ is the normalized redshift probability distribution
function and $\mathcal{S}_{\nu}(z)$ is the flux density “weighting” at
redshift $z$. In practice we have discrete redshift probability distributions
defined over $R$ bins such that, for a population of $N$ galaxies, the flux
density in the $ij$th pixel of a map can be written

$\displaystyle
M_{ij}=\sum^{N_{ij}}_{p}\sum^{R}_{q}\mathcal{S}_{\nu}(z_{q})_{p,ij}\mathcal{P}(z_{q})_{p,ij}\Delta
z$ (2)

Because of the PSF, the flux contribution of each galaxy is distributed over
many pixels according to the convolution

$\displaystyle\mathcal{M}=M\otimes{\rm PSF}.$ (3)

In effect equation 1 uses $\mathcal{P}(z)$ to split each of the $N$ galaxies
in our $K$-selected sample into $R$ redshift bins and assumes that the
galaxies in each redshift bin can be represented by an average observed flux
density, $\langle S_{\nu}(z)\rangle$. This is effectively the “stacked” flux
density.

With the model sky defined we consider an optimization problem where the set
of average observed flux densities, $\langle S_{\nu}(z)\rangle$, per redshift
interval in equation 1 are unknown coefficients. A key decision in defining
our sky model is in the binning of $\mathcal{P}(z)$. The UKIDSS-UDS
$\mathcal{P}(z)$ are binned in non-linear steps of
$1+z_{(n+1)}=1.001(1+z_{n})$. This would result in hundreds of free parameters
across the redshift range of interest, which is computationally impractical as
well as unnecessary given the photometric redshift uncertainties. Instead we
bin each $\mathcal{P}(z)$ to $\Delta z=0.5$, giving $20$ equally-sized bins
across the redshift range $0<z\lesssim 10$.

We aim to find the optimal set of average flux densities that minimizes the
square of the residual flux between the model in equation 3 and the observed
beam-convolved map, weighted by the noise. We use the Markov chain Monte Carlo
(MCMC) sampler emcee (Foreman-Mackey et al., 2013) to estimate the best fit
flux densities and their uncertainties. We minimize a negative log likelihood
$\ln(\mathcal{L})=-0.5\chi^{2}$, with

$\displaystyle\chi^{2}=\sum_{ij}\left(\frac{\mathcal{O}_{ij}-\mathcal{M}_{ij}}{\sigma_{{\rm
rms},{ij}}}\right)^{2}$ (4)

where $\mathcal{O}$ is the observed map and $\sigma_{\rm rms}$ is the
instrumental noise map. We initialise $1000$ “walkers” with an uninformative
prior, such that each walker is set with a vector of flux densities
(representing $S_{\nu}(z)$) with each flux density drawn from a Gaussian
distribution of mean $0.5$ mJy and width $0.05$ mJy. The sampler runs for
$1000$ iterations with the first $500$ iterations discarded as burn-in. The
best fitting flux densities and the $1\sigma$ bounds are estimated from the
$16{\rm th}$, $50{\rm th}$ and $84{\rm th}$ percentiles of accepted samples
for $500$ iterations. In Figures 1a, 1b and 1c we show the emcee corner plot
of the posterior distributions for all our free parameters.

Figure 1a: Standard emcee corner plot showing the one- and two-dimensional
posterior distributions for our parameters (the average observed flux density
in $\mu$Jy of galaxies in each redshift bin) for the redshift intervals
$\Delta z\,(0.25)$–$\Delta{z}\,(4.75)$. The density of the points and the
contours correlate with the posterior probability distributions from a
$1000$-step run (with $500$ steps for burn-in discarded) based on our sky
model and SCUBA-$2$ maps, and employing the “delete one” jackknife technique
with the map segment area corresponding to $A=1$ deleted to take into account
sample variance (Tukey, 1958). The vertical red lines show the average
$850\,\mu$m flux density of galaxies for each redshift interval, with the
dashed red lines showing the associated $1\sigma$ uncertainties (these values
are not corrected for the influence of the CMB, see Table 1 for CMB corrected-
estimates). As evidenced by this plot there are only very weak correlations
between our parameters for the redshifts intervals $\Delta
z\,(0.25)$–$\Delta{z}\,(4.75)$. Figure 1b: Standard emcee corner plot showing
the one- and two-dimensional posterior distributions for our parameters (the
average observed flux density in $\mu$Jy of galaxies in each redshift bin) in
the redshift intervals $\Delta z\,(5.25)$–$\Delta{z}\,(9.75)$. Detailed
description as in Figure 1a. As evidenced by this plot there are only very
weak correlations between our parameters for the redshifts intervals $\Delta
z\,(5.25)$–$\Delta{z}\,(9.75)$. Figure 1c: Standard emcee corner plot showing
the two-dimensional posterior distributions for our parameters (the average
observed flux density in $\mu$Jy of galaxies in each redshift bin) for
redshift intervals $\Delta z\,(0.25)$–$\Delta{z}\,(4.75)$ and $\Delta
z\,(5.25)$–$\Delta{z}\,(9.75)$. Detailed description as in Figure 1a. As
evidenced by this plot there are only very weak correlations between our
parameters for the redshifts intervals $\Delta z\,(0.25)$–$\Delta{z}\,(4.75)$
and $\Delta z\,(5.25)$–$\Delta{z}\,(9.75)$.

To estimate the additional uncertainty on the stacked flux densities due to
sampling variance we employ the “delete one” jackknife technique (Tukey,
1958), splitting the map into $A=21$ approximately equal area sectors and
running the MCMC fit for each jackknife. We find that the sampler chains
converge quickly (within $200$ steps), and tests indicate that the best fit
parameters are insensitive to the initialisation parameters. The covariance
matrix is given by

$\displaystyle\mathcal{C}_{ij}=\frac{A-1}{A}\sum_{i=1}^{A}\left(S^{k}_{i}-\bar{S}_{i}\right)\left(S^{k}_{j}-\bar{S}_{j}\right)$
(5)

where $S^{k}_{i}$ is the average flux density in the $i$th redshift bin,
eliminating the $k$th sample and $\bar{S}_{i}$ is the average over all
samples. The $1\sigma$ uncertainties on the stacked fluxes are estimated by
the square root of the diagonal elements of $\mathcal{C}$.

At high-$z$ the increase in the cosmic microwave background (CMB) temperature
affects the measurement of submillimeter dust continuum in two ways (see da
Cunha et al., 2013, for a detailed discussion). Firstly, the CMB provides an
additional source of dust heating increasing the intrinsic dust temperature as
shown in equation 6 (da Cunha et al., 2013):

$\displaystyle
T_{\mathrm{dust}}(z)=\Big{(}(T^{z=0}_{\mathrm{dust}})^{\beta+4}+(T^{z=0}_{\mathrm{CMB}})^{\beta+4}[(1+z)^{\beta+4}-1]\Big{)}^{\frac{1}{\beta+4}}$
(6)

Secondly, submillimeter observations of dust emission are always measured
against the background of the CMB. At low-z $T_{\rm dust}(z)>>T_{\rm cmb}(z)$,
so essentially all the intrinsic flux is detected against the CMB. However, at
high-z, as $T_{\rm cmb}(z)$ approaches $T_{\rm dust}(z)$, the fraction of
submillimeter flux detected against the CMB decreases. Equation 7 (da Cunha et
al., 2013) shows the fraction of the intrinsic dust emission from a galaxy
measured at a given frequency, $v_{\mathrm{obs}}=v_{\mathrm{rest}}/(1+z)$
against the CMB background:

$\displaystyle\frac{F^{\mathrm{\,obs\,against\,CMB}}_{v_{\mathrm{obs}}}}{F^{\mathrm{intrinsic}}_{v_{\mathrm{obs}}}}=1-\frac{B_{v}[T_{\mathrm{CMB}}(z)]}{B_{v}[T_{\mathrm{dust}}(z)]}$
(7)

Assuming $T^{z=0}_{\mathrm{CMB}}=2.73\,$K, $T^{z=0}_{\mathrm{dust}}=25\,$K (in
line with the $T_{\mathrm{dust}}$ adopted in Scoville et al., 2016) and
$\beta=2$ we derive the fraction of submillimeter flux observed against the
CMB for the redshift range $0<z<10$ at $v_{\mathrm{obs}}{=}353$ GHz
($\lambda_{\mathrm{obs}}{=}850\,\mu$m), including the extra heating
contributed by the CMB. We apply this correction to our average observed
$850\,\mu$m flux densities in all redshift bins to account for the impact of
the CMB on our estimates.

##  4 Results

###  4.1 Estimating molecular gas mass: RJ luminosity-to-gas mass relation

In Table 1 we present the average observed $850\,\mu$m flux densities for our
galaxy sample as a function of redshift. We quote the uncertainties to
$1\sigma$ and include the additional uncertainty due to sample variance. We
note that at $z\gtrsim 6$ the UDS redshifts are untested. However, as our
sample is binned according to the $\mathcal{P}(z)$ for each source, every
galaxy effectively contributes to the flux in each redshift interval. Hence,
we show the average $850\,\mu$m flux density estimates for our galaxy sample
to $z=10$.

We sum $\mathcal{P}(z)$ (which is completeness corrected) in each redshift
bin, giving us the galaxy “weighting” for each redshift interval. The integral
of the summed $\mathcal{P}(z)$ across all redshift intervals should be
approximately equal to the total number of galaxies in our sample. We
calculate this to be $154,839$, which is consistent with our galaxy sample of
$153,399$ sources (taking into account the completeness corrections). With the
summed $\mathcal{P}(z)$ and taking the area of our sample as the unmasked
region of the SCUBA-$2$ $850\,\mu$m map (which corresponds to the UDS binary
mask for “good galaxy” regions) we calculate the number density of galaxies as
a function of redshift. By combining this with our average flux density (see
Table 1 column $4$ for CMB corrected values) we calculate the summed flux
density for our galaxy sample in each redshift interval. In Figure 2 we
present the number density and summed flux density for galaxies in our sample
as a function of redshift. The distribution of the summed flux densities with
redshift is broadly comparable to the redshift distribution found for SMGs,
which peaks at $z\approx 2$ (e.g., Blain et al., 2002; Chapman et al., 2005;
Simpson et al., 2014; Miettinen et al., 2017; Zavala et al., 2018), whilst the
number density distribution generally declines with increasing redshift as
expected. The difference in the evolution of these distributions demonstrates
that our derived flux densities are not biased by the number density of
galaxies in each redshift bin.

Figure 2: Number density of galaxies in each redshift bin (blue bars)
calculated from the galaxy “weighting” (the sum of the completeness corrected
$\mathcal{P}(z)$ in each redshift bin) and the area of the unmasked region of
the SCUBA-$2$ image. The black points show the summed flux density of galaxies
as a function of redshift derived by combining the number density and the
average flux density of galaxies (see Table 1 column $4$ for average flux
density values) in each redshift bin. This figure shows the contrasting
distributions, with the summed flux density mirroring the redshift
distribution found for SMGs with a peak at $z\approx 2$, whereas the number
density generally shows a steady decline. At $z\approx 7$ there is an
unexpected upturn in the number density and the summed flux density of
galaxies, this is likely due to a combination of uncertainties in the
photometric redshift fitting at high-$z$ (UDS photometric redshifts are
untested at $z\gtrsim 6$) and the small number statistics in these bins with
the galaxy weighting being $\lesssim 100$ per redshift interval. Table 1:
Redshift intervals, galaxy “weighting” in each redshift interval, average 850
$\mu$m flux density, and average 850 $\mu$m CMB corrected flux density.

$\Delta z$ |  Galaxy weighting |  $\langle S_{850}\rangle$ |  CMB corr. $\langle S_{850}\rangle$
---|---|---|---
|  | [$\mu$Jy dz-1] | [$\mu$Jy dz-1]
$0.25$ | $90003$ | $2.1\pm 2.0$ | $2.1\pm 2.0$
$0.75$ | $79917$ | $32.9\pm 8.2$ | $33.0\pm 8.2$
$1.25$ | $64475$ | $67.5\pm 10.3$ | $68.0\pm 10.3$
$1.75$ | $47922$ | $146.9\pm 12.8$ | $148.5\pm 12.9$
$2.25$ | $31448$ | $135.0\pm 17.9$ | $137.2\pm 18.2$
$2.75$ | $18368$ | $224.9\pm 21.5$ | $230.3\pm 22.0$
$3.25$ | $13071$ | $104.6\pm 23.3$ | $108.3\pm 24.1$
$3.75$ | $5271$ | $236.5\pm 45.1$ | $248.4\pm 47.4$
$4.25$ | $2413$ | $195.0\pm 52.9$ | $209.1\pm 56.7$
$4.75$ | $848$ | $130.8\pm 109.6$ | $144.3\pm 120.9$
$5.25$ | $336$ | $251.2\pm 109.8$ | $288.4\pm 126.1$
$5.75$ | $79$ | $422.9\pm 362.7$ | $511.9\pm 439.0$
$6.25$ | $86$ | $74.2\pm 95.7$ | $96.2\pm 124.1$
$6.75$ | $37$ | $1606.2\pm 868.8$ | $2276.3\pm 1231.3$
$7.25$ | $49$ | $60.2\pm 77.3$ | $95.2\pm 122.1$
$7.75$ | $38$ | $106.0\pm 127.0$ | $190.8\pm 228.6$
$8.25$ | $91$ | $63.2\pm 73.7$ | $132.2\pm 154.1$
$8.75$ | $121$ | $29.0\pm 36.5$ | $71.8\pm 90.3$
$9.25$ | $73$ | $155.2\pm 152.5$ | $460.6\pm 452.5$
$9.75$ | $70$ | $289.5\pm 161.3$ | $1039.0\pm 578.9$

Note. — Col. 1. Mid-point of each redshift interval ($\Delta z$); col. 2. the
galaxy “weighting” which is the sum of the completeness corrected
$\mathcal{P}(z)$ in each redshift bin. The summed $\mathcal{P}(z)$ across all
redshift intervals integrates to $154,839$, which is consistent with the
number of galaxies in our sample ($153,399$) taking into account the
completeness corrections; col. 3. average (stacked) $850\,\mu$m flux density
as a function of redshift with the uncertainty quoted to $1\sigma$; col. 4.
average $850\,\mu$m flux density corrected for the impact of the CMB as a
function of redshift with the uncertainty quoted to $1\sigma$.

We adopt the approach of Scoville et al. (2016, hereafter S16) and, utilising
our flux density measurements, estimate the average molecular gas mass for our
galaxy sample in each redshift interval. The full details of this approach are
given in S16; however, we provide a brief description here.

The long wavelength RJ tail of dust emission is nearly always optically thin
($\tau\ll 1$) and consequently this provides a direct probe of the total dust
mass and hence the molecular gas mass. S16 utilise this to obtain an
empirically calibrated RJ luminosity-to-gas mass ratio

$\displaystyle\left(\frac{M_{\rm mol}}{\rm
M_{\odot}}\right)=\frac{1}{\alpha_{850}}\left(\frac{L_{\rm 850,\ rest}}{\rm
erg\,s^{-1}\,Hz^{-1}}\right){\rm for}\ \lambda_{\rm rest}\gtrsim 250\,\mu{\rm
m}$ (8)

with $\alpha_{850}=6.7\pm 1.7\times 10^{19}\,{\rm erg}\,{\rm s}^{-1}\,{\rm
Hz}^{-1}\,{\rm M_{\odot}}^{-1}$. The restriction $\lambda_{\rm rest}\gtrsim
250\,\mu{\rm m}$ is required to ensure that at an observed wavelength of
$850\,\mu$m the rest-frame emission stays on the RJ tail. S16 demonstrate that
this luminosity-to-mass ratio is relatively constant for high-stellar mass
($M_{\rm stellar}=(2-40)\times 10^{10}\,{\rm M}_{\odot}$) normal star-forming
and star-bursting galaxies, both locally and at high-z.

We estimate the average rest-frame $850\,\mu$m luminosity density of galaxies
as a function of redshift using the average flux densities detailed in Table
1, assuming a mass-weighted dust temperature of $25$ K 111The S16 calibration
uses a mass-weighted temperature of $25$ K, rather than a luminosity-weighted
dust temperature (see Appendix A.$2$ of S16) and employing the relation from
S16:

$L_{\nu 850}=S_{\nu}[\rm Jy]\ \times\ 1.19\ \times 10^{27}\ \times\
\left(\frac{\nu(850\,\mu{\rm m})}{\nu_{obs}(1+z)}\right)^{3.8}\\\ \times\
\frac{(d_{\rm L}[{\rm Mpc}])^{2}}{1+z}\ \times\ \frac{\Gamma_{\rm
RJ}(25,\nu_{850\,\mu{\rm m}},0)}{\Gamma_{\rm RJ}(25,\nu_{\rm obs},z)}\ [{\rm
erg\,s^{-1}\,{\rm Hz}^{-1}}]$ (9)

The $\Gamma_{\rm RJ}$ term in equation 10 corrects for departures from the RJ
$\nu^{2}$ dependence as the observed emission approaches the spectral energy
distribution (SED) peak in the rest frame and where $T_{\rm dust}$ is the
mass-weighted temperature characterizing the RJ dust emission

$\displaystyle\Gamma_{\rm RJ}(T_{\rm dust},\nu_{\rm obs},z)=\frac{h\nu_{\rm
obs}(1+z)/kT_{\rm dust}}{e^{h\nu_{obs}(1+z)/kT_{\rm dust}}-1}$ (10)

Table 2: Redshift intervals, average rest-frame $850\mu$m luminosity, average
molecular gas mass, co-moving molecular gas mass density as a function of $z$.

$\Delta z$ |  $\langle L_{\rm 850,\,rest}\rangle$ |  $\langle M_{\rm H_{2}}\rangle$ |  $\rho_{\rm H_{2}}$ |  $\Omega_{\rm H_{2}}$
---|---|---|---|---
|  [$10^{27}$ erg s-1 Hz-1 dz-1] |  [$10^{7}\,\rm{\rm M_{\odot}}$ dz-1] |  [$10^{6}\,{\rm M_{\odot}}\,{\rm Mpc}^{-3}$] |  [$10^{-7}$]
$0.25$ | $1.64\pm 0.02$ | $1.80\pm 0.03$ | $1.92\pm 1.80$ | $115.66\pm 108.42$
$0.75$ | $81.97\pm 3.89$ | $89.95\pm 4.26$ | $17.12\pm 4.39$ | $569.51\pm 146.13$
$1.25$ | $217.79\pm 18.32$ | $239.01\pm 20.10$ | $23.22\pm 4.09$ | $430.70\pm 75.93$
$1.75$ | $514.61\pm 64.72$ | $564.77\pm 71.02$ | $34.67\pm 5.32$ | $380.54\pm 58.41$
$2.25$ | $492.12\pm 85.19$ | $540.08\pm 93.49$ | $20.74\pm 4.60$ | $143.48\pm 31.84$

Note. — Col. 1. Mid-point of each redshift interval ($\Delta z$); cols. 2 and
3. the average rest-frame $850\,\mu$m luminosity and average molecular gas
mass for galaxies in each redshift interval; col. 4. molecular gas mass
density as a function of redshift and col. 5. molecular gas mass density in
terms of the critical mass density. We restrict our results to $z\lesssim 2.5$
to ensure that the observed $850\,\mu$m dust emission is tracing the RJ tail.

We restrict our estimates of the average rest-frame $850\,\mu$m luminosity to
$z\lesssim 2.5$ to ensure that rest-frame emission stays on the RJ tail. With
the average rest-frame $850\,\mu$m luminosity density derived (detailed in
Table 2) we use the RJ luminosity-to-gas mass ratio from equation 8 to
estimate the average molecular gas mass as a function of redshift to $z\approx
2.5$. This calibration includes a factor of 1.36 to account for the associated
mass of heavy elements (mostly Helium at $8\%$ by number), so we correct our
results by a factor $1/1.36$ ($M_{\rm mol}$) to obtain $M_{\rm H_{2}}$.

Since the summed photometric redshift probability distributions inform us
about the galaxy “weighting” in each redshift interval and the UDS binary mask
for “good” regions gives us the unmasked area of the SCUBA-$2$ $850\,\mu$m
map, we combine this information with the average molecular gas mass and
differential co-moving volume element to estimate the co-moving volume density
of molecular gas

$\displaystyle\rho_{\rm H_{2}}=\Omega\int_{z-\Delta z/2}^{z+\Delta
z/2}N(z)\langle M_{\rm H_{2}}(z)\rangle\frac{dV}{dzd\Omega}dz.$ (11)

We present our values for $\rho_{\rm H_{2}}(z)$ in Table 2 as a function of
redshift, also giving this in terms of the critical mass density $\Omega_{\rm
H_{2}}=\rho_{\rm H_{2}}(z)/\rho_{\rm crit}(z)$. We use a Monte Carlo analysis
to calculate the uncertainties for our values of ${L_{\rm 850,\,rest}}$,
$M_{\rm H_{2}}$ and $\rho_{\rm H_{2}}(z)$, first drawing random values for
$S_{\nu}$ from a Gaussian distribution where the mean is the average flux
density and width the uncertainty on the average flux density, and then
drawing values for a mass-weighted $T_{\rm dust}$ from a Gaussian distribution
with a mean of $25$ K (corresponding to the constant $T_{\rm dust}$ assumed by
S16), and width $3$ K. Observations (e.g., ; , ) and simulations (e.g., Liang
et al., 2018, 2019) find that a mass-weighted $T_{\rm dust}$ shows little
variation with galaxy $L_{\rm 850}$ or redshift (cf. Behrens et al., 2018),
and by utilising a temperature distribution with $\sigma=\pm 3$ K we recognise
this minimal variance in our uncertainty calculations. We use these values to
estimate ${L_{\rm 850,\,rest}}$, $M_{\rm H_{2}}$ and $\rho_{\rm H_{2}}(z)$
from equations 9, 8, and 11 respectively for $1000$ runs, with the uncertainty
being taken as the standard deviation across these trials.

We note that the galaxy sample we use to derive the results in Tables 1 and 2
includes all galaxies in the “good galaxy” subset of the UDS DR11 catalog,
regardless of the reliability of photometric redshifts for individual sources.
If we apply a $\chi^{2}$ cut to exclude galaxies with the least reliable
redshifts (omitting galaxies with a reduced $\chi^{2}$ value for the
photometric redshift of $>10$) and repeat the process outlined in Sections 3
and 4 above we find a less than $2\%$ variation in our results with estimates
consistent with those in Table 2 within the uncertainties.

We plot $\rho_{\rm H_{2}}(z)$ in Figure 3, compared to direct CO line
estimates from ASPECS (Decarli et al., 2016b, 2020), COLDZ (Riechers et al.,
2019) and VLASPECS (Riechers et al., 2020), as well as values derived using
far-infrared and UV photometry (Berta et al., 2013). We fit a function of the
same form as the star-formation rate density function presented in Madau &
Dickinson (2014) to the log of our results, and derive the best fit parameters
for our data using emcee. This yields:

${\rm log_{10}}\left(\frac{\rho_{\rm H_{2}}}{\rm M_{\odot}\,{\rm
Mpc^{-3}}}\right)=\\\ (6.59\pm 0.30)\times\frac{(1+z)^{0.38\pm
0.16}}{1+[(1+z)/5.57\pm 1.69]^{1.78\pm 0.61}}$ (12)

which we plot in Figure 3. Our results show a peak $\rho_{\rm H_{2}}(z)$ at
$z\approx 2$ mirroring existing constraints.

###  4.2 Deriving molecular gas mass density from the halo mass function.

Using an alternative approach we derive $\rho_{\rm H_{2}}(z)$ from first
principles using the halo mass function from Murray et al. (2013) and assuming
a constant halo mass range of $10^{11.5}$–$10^{15}\,{\rm M_{\odot}}$. We
estimate the molecular gas mass density as a function of halo mass (for
redshifts $0\leq z\leq 7$) using the stellar-halo mass ratio from Moster et
al. (2013) and the ISM-stellar mass relation from Scoville et al. (2017). The
ISM-stellar mass relation is calibrated using a sample of high mass galaxies
($M_{\rm{stellar}}\gtrsim 10^{10}\,{\rm M}_{\odot}$), therefore we adopt a
halo mass range for which the corresponding stellar masses are comparable with
the Scoville et al. (2017) calibration sample. Integrating these estimates
with respect to halo mass gives the total molecular gas density as a function
of redshift, which we present in Figure 3.

###  4.3 Estimating molecular gas mass density using a “constant efficiency”
model

We also estimate $\rho_{\rm H_{2}}(z)$ from the star-formation rate density,
$\rho_{\dot{M_{\star}}}(z)$, assuming a corresponding volume averaged star-
formation “efficiency”, $\eta(z)=\rho_{\rm H_{2}}(z)/\psi_{\star}(z)$. We use
the functional fit of Wilkins et al. (2019), a recalibration of the well-known
Madau & Dickinson (2014) cosmic star formation history. We make the assumption
that $\eta(z)$ is constant and that the total molecular gas mass per galaxy
can be related to on-going star formation as $\xi M_{\rm H_{2}}={\rm
SFR}/\epsilon$ (Geach & Papadopoulos, 2012). Here $\xi$ is the ratio of dense,
actively star-forming molecular gas to the total molecular reservoir with
$\xi\approx 0.04$ for quiescent disks and $\xi>0.5$ for starbursts (e.g.,
Papadopoulos & Geach, 2012), while the factor $\epsilon$ describes the rate at
which the dense molecular gas forms stars. Figure 3 shows the predicted
$\rho_{\rm H_{2}}(z)$ inferred from the Wilkins et al. (2019) fit, assuming a
constant “average” $\eta(z)=0.3$ Gyr corresponding to $\xi=0.1$ and
$\epsilon=37$ Gyr-1 (e.g., Geach & Papadopoulos, 2012).

This value for $\eta(z)$ is similar to the typical values of $t_{\rm dep}$
($M_{\mathrm{H}_{2}}/\mathrm{SFR}$) $\approx 1$ Gyr (e.g., Tacconi et al.,
2018) for main-sequence galaxies. Tacconi et al. (2018) find a relatively weak
dependence of $t_{\rm dep}$ with redshift, $t_{\rm dep}\propto(1+z)^{-0.57}$
to $z<2.5$, consistent with our picture of a common mode of star formation in
normal galaxies, at least out to the peak epoch. We use this relation from
Tacconi et al. (2018) and the Wilkins et al. (2019) fit to derive an estimate
of $\rho_{\rm H_{2}}(z)$ which incorporates a weakly evolving star-formation
efficiency. We present our predicted $\rho_{\rm H_{2}}(z)$ in Figure 3.

Figure 3: Values for $\rho_{\rm H_{2}}(z)$ (CMB corrected) derived using a
$3$D stacking method and the RJ luminosity-to-gas mass ratio of S16. The upper
x-axis shows the rest-frame wavelength of observed $850\,\mu$m emission for
the redshift range shown, illustrating the range at which rest-frame emission
traces the long-wavelength RJ tail (i.e. $\lambda_{\rm rest}\gtrsim 250\mu$m)
and the S16 calibration can be reliably applied. Our estimates are represented
by the red points (uncertainties shown to $2\sigma$) with the solid red line
showing the best-fit function derived using emcee and a function of the same
form as the star-formation rate density function from Madau & Dickinson
(2014). We note that at $z\gtrsim 2.5$ our results should be considered
speculative as the observed $850\,\mu$m emission no longer traces the rest-
frame RJ tail. Hence we only show an average of our estimates after this
point. We do not include any estimates for $z\gtrsim 6$ as the UDS photometric
redshifts beyond this are untested and as such highly uncertain. Alongside our
values we show results from ASPECS (off white rectangles; Decarli et al.,
2016b), ASPECS LP (blue rectangles; Decarli et al., 2020), COLDz (green
rectangles; Riechers et al., 2019) and VLASPECS (green diamond; Riechers et
al., 2020) which are derived from direct measurements of CO line emission. The
grey points show values from Berta et al. (2013), estimated using deep far-
infrared and UV data, and assuming either typical gas depletion times (Tacconi
et al., 2013) or from IR luminosity and obscuration properties (Nordon et al.,
2013). We plot our “constant efficiency” models derived using the re-
calibrated star formation history from Wilkins et al. (2019) and assuming
either a corresponding constant (e.g., Geach & Papadopoulos, 2012) or weakly
evolving (e.g., Tacconi et al., 2018) volume averaged star-formation
“efficiency” to infer $\rho_{\rm H_{2}}(z)$. The former is shown as the dark
blue dot dash line and the latter illustrated by the dark blue circles. We
also plot $\rho_{\rm H_{2}}(z)$ to $z\approx 7$ derived from the halo mass
function (Murray et al., 2013), assuming the stellar-halo mass ratio from
Moster et al. (2013) and ISM-stellar-mass relation from Scoville et al.
(2017). The dotted dark blue line corresponds to a halo mass range
$10^{11.5}$–$10^{15}\,{\rm M_{\odot}}$, with the minimum stellar masses
derived for this range (${\approx}10^{9.5}\,{\rm M_{\odot}}$) being consistent
with the lowest stellar masses probed in ASPECS LP (Boogaard et al., 2019).

##  5 Discussion

Our results appear to be in reasonable agreement with existing empirical
constraints, indicating that the epoch of molecular gas coincided with the
peak epoch of star formation at $z\approx 2$ . So what does this mean in terms
of the evolving molecular gas budget? We might ask what is the complete
picture of $\rho_{\rm H_{2}}(z)$, or rather, what galaxies host the majority
of the cosmic molecular gas budget across cosmic time? In the following
discussion we interpret the overall evolution of the cosmic molecular gas
density, in the context of our results, within the established framework of
star formation in galaxies from the cosmic dawn to the present day.

###  5.1 Evolution of cosmic molecular gas mass density at $0\leq z\leq 2.5$

The S16 RJ luminosity-to-gas mass ratio has been shown to provide molecular
gas mass estimates accurate to within a factor of around $2$ when compared
with measurements made via direct CO ($J=1\rightarrow 0$) line observations
(e.g., Scoville et al., 2017; Kaasinen et al., 2019), with variations in the
dust emissivity index, temperature, and gas-to-dust ratios being accountable
for the deviations. This factor of $2$ accuracy is based on samples of
galaxies with high stellar masses $M_{\rm stellar}=(2-40)\times 10^{10}\,{\rm
M_{\odot}}$ as these galaxies are likely to have near-solar metallicity
(Tremonti et al., 2004). This avoids probing low metallicity sources for which
the dust-to-gas abundance ratio is likely to drop or the CO gas fraction is
low (Bolatto et al., 2013).

In Figure 4 we show the marginalised stellar mass estimates for UDS galaxies
(Almaini et al. in prep.) and the corresponding $95\%$ stellar mass
completeness (derived using the method of Pozzetti et al., 2010, Wilkinson et
al. in prep.) for the UDS catalog. As can be seen in Figure 4 our galaxy
sample includes a proportion of galaxies with stellar masses lower than those
used to derive the S16 RJ luminosity-to-gas mass ratio, with these sources
being more abundant in lower redshift bins.

The dust-to-gas relation has been found to be relatively consistent for nearby
galaxies with $M_{\rm stellar}>10^{9}\,{\rm M_{\odot}}$ (e.g., Groves et al.,
2015), but drops for galaxies with lower stellar masses (hence lower
metallicities). Cosmological galaxy formation simulations have shown that
deviations from this relation become significant ($\gtrsim 0.5\,{\rm dex}$) at
$L_{\rm 850}\lesssim 10^{28}$ erg s-1 Hz-1 in the redshift range $0<z<9.5$
(e.g., Privon et al., 2018). As shown in Figure 4 at $z\gtrsim 1$ the majority
of our sample are likely to have $M_{\rm stellar}>10^{9}\,{\rm M_{\odot}}$ and
Table 2 shows that the mean rest-frame $850\,\mu$m luminosity for our sample
in all but the lowest redshift interval ($\Delta z=0.25$) is $L_{\rm 850\mu
m}\gtrsim 10^{28}$ erg s-1 Hz-1. Therefore, whilst the RJ luminosity-to-gas
mass ratio has been calibrated on high stellar mass galaxies ($M_{\rm
stellar}=(2-40)\times 10^{10}\,{\rm M_{\odot}}$), we make the assumption that
applying this calibration to our sample at redshifts $z\gtrsim 0.5$ is likely
to result in comparable uncertainties (i.e. a factor of $2$). However, in the
redshift bin $\Delta z=0.25$ the mean rest-frame luminosity is $\langle L_{\rm
850\mu m}\rangle=1.6\times 10^{27}$ erg s-1 Hz-1, and as such our results are
likely to be under-predicted by $\gtrsim 0.5$ dex (e.g., Privon et al., 2018)
due to the abundance of lower mass (low metallicity) galaxies in this redshift
interval.

Figure 4: Marginalised stellar mass estimates (Almaini et al. in prep.) for
UDS galaxies in our sample as a density plot, with darker colours
corresponding to higher number densities of galaxies. The red line corresponds
to the UDS catalog $95\%$ stellar mass completeness (derived using the method
of Pozzetti et al., 2010, see Wilkinson et al. in prep.). This figure
demonstrates that at high redshifts we are only sensitive to the most massive
galaxies.

###  5.2 Comparison of the evolution of molecular gas mass density to other
studies in the literature

In Figure 3 we compare our results, which are revised to account for the
influence of the CMB, to those from direct CO line surveys (e.g., Decarli et
al., 2016b, 2020; Riechers et al., 2019, 2020). We limit our discussion to the
results from these surveys at $z\lesssim 2.5$ as we are restricted to this
redshift range due to the wavelength of our observations
($\lambda_{\mathrm{obs}}=850\,\mu\mathrm{m}$). Whilst the results from
ASPECS/ASPECS LP (Decarli et al., 2016b, 2020) and COLDz (Riechers et al.,
2019, 2020) were not corrected for the influence of the CMB, we note that this
is not necessary at $z\lesssim 4.5$ as the effect of the CMB on measurements
of molecular gas mass density from direct CO line observations is minimal
($\lesssim 15\%$, e.g., Decarli et al., 2019) and as such this does not impact
on our analysis here.

We find that our results are broadly consistent with the estimates from direct
CO line surveys within uncertainties and show notably good agreement with
results obtained through observations of the ground state CO line (ASPECS LP
at $\Delta z\simeq 0.25$ and COLDZ at $\Delta z\simeq 2.25$). Albeit, we
caution that our results at $\Delta z\simeq 0.25$ are likely under estimated
due to the abundance of low stellar mass galaxies in this redshift bin.

When compared with the ASPECS LP survey (Decarli et al., 2020) our results
generally trace the lower boundaries of their estimates between $0.75\lesssim
z\lesssim 1.75$. The S16 RJ luminosity-to-gas mass ratio is calibrated using
the ground state CO (${J=1\rightarrow 0}$) line, whereas at $z>0.75$ the
ASPECS LP results are derived from observations of higher state excitation CO
lines. Therefore, this offset could be explained by the uncertainties
associated with translating higher excitation CO lines observations to ground
state CO (${J=1\rightarrow 0}$) luminosities. However, in the redshift
interval $\Delta z=1.25$, even if the extreme case of thermalized gas is
assumed, our upper estimate (taking into account $1\sigma$ uncertainties)
falls a factor of $1.42$ below the lower boundary of the ASPECS LP survey
(e.g., $38.90\times 10^{6}\,{\rm M_{\odot}}\,{\rm Mpc}^{-3}$, Table A3.
Decarli et al., 2019). As such the uncertainties in CO line ratios do not
fully account for the offset we see.

Building on previous studies Liu et al. (2019) derive the molecular gas mass
density using a dataset comprised of ${\approx}700$ ALMA continuum detected
galaxies and ${\approx}1000$ galaxies with CO observations (taken from the
literature). To derive molecular gas masses for the continuum detected
galaxies Liu et al. (2019) employ the Hughes et al. (2017) luminosity-to-gas
mass calibration. Liu et al. (2019) estimate a SMF (stellar mass function)
integrated molecular gas mass density based on the SMF integrated to $M_{\rm
stellar}=10^{9}\,{\rm M_{\odot}}$ and using a gas fraction function derived
from their composite sample of ${\approx}1700$ galaxies. Their results trace
the upper boundaries of the molecular gas mass density derived from the most
recent blank field CO line surveys (e.g., Decarli et al., 2019; Riechers et
al., 2019) and are ${\approx}1$ dex higher than our estimates. This offset
with our estimates of molecular gas mass density could be in part due to
assumptions made in Liu et al. (2019) to derive an SMF integrated molecular
gas mass density (i.e. that all star-forming galaxies are on the main-
sequence) or potentially differences in sample selection. The majority
(${\approx}800$) of the CO detected sources in the Liu et al. (2019) composite
sample are in the Local Universe (i.e $z<0.3$), so at $z>0.3$ their dataset is
dominated by ALMA continuum detected galaxies which are preferentially massive
and dust-rich (and hence, using a dust-to-gas mass conversion, gas-rich). In
contrast our galaxy sample is near-infrared selected and as such our selection
is less likely to sample these luminous dust-rich SMGs, with previous studies
finding that $\approx{20}\%$ of SMGs are missed in optical/new-infrared
surveys (e.g. Dudzevičiūtė et al., 2019).

Magnelli et al. (2020) use a stacking approach to measure the comoving gas
mass density of a sample of $555$ near-infrared selected galaxies, with
galaxies split into bins of $z$ and ${M}_{\rm stellar}$. Their stacking method
accounts for the metallicity of galaxies in these bins (inferred using the
stellar mass-metallicity from Tacconi et al., 2018) and is equivalent to the
S16 calibration at solar metallicity. Our results trace the lower boundaries
of their estimates, but are inconsistent (within $1\sigma$ uncertainties) at
$z>1$. This discrepancy may be in part due to our method not accounting for
the metallicity of low mass galaxies, resulting in an under-estimation of the
molecular gas mass for galaxies with $M_{\rm stellar}<10^{9}\,{\rm
M_{\odot}}$. However, if this was the sole reason for the difference in our
results we would expect this to have more of an impact at $z<1$ where this
effect will be more prominent.

We caution that our results also rely solely on photometric redshifts, which
despite the high quality 12 band photometry of the UDS catalog, cannot compete
with the accuracy of redshifts derived via spectroscopic surveys. By utilising
the redshift probability distributions in our 3D stacking approach we aim to
provide mitigation against these uncertainties. However, whilst our estimates
rely exclusively on the use of photometric redshifts, the results obtained in
both ASPECS LP (Decarli et al., 2020; Magnelli et al., 2020) and Liu et al.
(2019) benefit from the inclusion of sources with more reliable spectroscopic
redshifts. This may also play a part in deviations seen when we compare our
estimates with these previous surveys.

###  5.3 Contribution of the brightest submillimeter sources to the cosmic
evolution of the molecular gas mass density

In order to present the most complete view of the evolution of the molecular
gas mass density we stack all sources in our near-infrared selected sample,
including counterparts to the bright submillimeter sources in the SCUBA-$2$
UDS map. To test the contribution of these galaxies to our results we repeat
our stacking analysis with the SCUBA$2$ UDS source subtracted map. As expected
excluding the $\sim 1000$ UDS submillimeter sources reduces the average
observed $850\,\mu$m flux in our redshift intervals, which propagates to our
estimate of the comoving molecular gas mass density. At $z<1.5$ the exclusion
of the UDS submillimeter sources has a minimal impact on our estimates of
$\rho_{\rm H_{2}}(z)$ and these remain consistent within the $1\sigma$
uncertainties. However, at $z>1.5$ our estimates of $\rho_{\rm H_{2}}(z)$ drop
by a factor of $2.05$ and $2.33$ in the redshift intervals $\Delta z=1.75$ and
$\Delta z=2.25$ respectively. This coincides with the peak number density of
SMGs at $z\approx 2$. This indicates that approximately $50\%$ of the
molecular gas mass density at the peak of the star formation rate density is
locked in dust-rich SMGs. We note that our inferred contribution of SMGs is
also likely under-estimated as we expect that approximately $20\%$ of SMGs are
undetected in our near-infrared selected sample (e.g., Dudzevičiūtė et al.,
2019). Our finding is in keeping with Zavala et al. (2021) who find that
bright SMGS ($L_{\rm IR}>10^{12}L_{\odot}$) dominate the obscured star
formation rate density at $z\approx 2$ and also Magnelli et al. (2020) who
find that the bulk of dust and gas in galaxies is locked in massive star-
forming galaxies.

###  5.4 Additional constraints on the evolution of molecular gas mass density

We have added further valuable constraints to this picture of cosmic molecular
gas evolution using two alternative approaches.

We estimate $\rho_{\rm H_{2}}(z)$ from the halo mass function (Murray et al.,
2013) assuming a constant halo mass range of $10^{11.5}$–$10^{15}\,{\rm
M_{\odot}}$, and using the stellar-halo mass ratio from Moster et al. (2013)
and the ISM-stellar mass relation from Scoville et al. (2017). For the latter
relation we make the assumption that all galaxies are on the star-forming main
sequence (e.g., $\rm{sSFR/sSFR_{MS}}=1$). The evolution of our halo mass
derived $\rho_{\rm H_{2}}(z)$(shown in Figure 3) follows a similar shape to
the star-formation rate density, rising to a peak at $1\lesssim z\lesssim 3$
and decreasing to the present day. The minimum halo mass we assume corresponds
to stellar masses of ${\approx}10^{9.5}\,{\rm M_{\odot}}$ (e.g., Moster et
al., 2013), equivalent to the lowest stellar masses probed in the ASPECS LP
survey (Decarli et al., 2020). Our $\rho_{\rm H_{2}}(z)$ estimates show good
agreement with the ASPECs/ASPECS LP surveys (Decarli et al., 2016a, 2020) at
$z\gtrsim 0.6$. However, as shown in Figure 3 at $z\lesssim 0.6$ our estimate
of $\rho_{\rm H_{2}}(z)$ lies above the lowest redshift bins from the ASPECS
LP survey (Decarli et al., 2020) and is ${\approx}1\,\rm{dex}$ higher than our
estimate of $\rho_{\rm H_{2}}(z)$ derived from measurements of observed
$850\,\mu$m flux. To obtain an estimate of $\rho_{\rm H_{2}}(z)$ from the halo
mass function we make the assumption that all galaxies are star-forming. As
such our estimate of $\rho_{\rm H_{2}}(z)$ derived from the halo mass function
can be seen as an upper limit $\rho_{\rm H_{2}}(z)$ for the stellar mass range
sampled. It follows that we see a more significant deviation between our
estimate and observationally derived results at lower redshifts as the
fraction of passive galaxies is higher at later epochs.

We also estimate $\rho_{\rm H_{2}}(z)$ from the star-formation rate density
(Wilkins et al., 2019), assuming a constant (Geach & Papadopoulos, 2012) and
weakly evolving (Tacconi et al., 2018) star-formation efficiency. These
“constant efficiency” models predict a co-moving molecular gas mass density in
good agreement with both measurements of molecular gas mass via observations
of direct CO line emission (Decarli et al., 2016b, 2020; Riechers et al.,
2019, 2020) and our results derived from measurements of the long-wavelength
dust emission, out to a peak at $z\approx 2$. A simple conclusion is that the
peak epoch of star formation at $z\approx 2$ is not driven by significantly
more efficient (or starburst-like) star formation in galaxies, but by a higher
abundance of molecular fuel in galaxies. We note that the estimate derived
from weakly evolving star-formation is ${\approx}1\,\rm{dex}$ higher than our
results at $\Delta z\simeq 0.25$. This is likely a consequence of the latter
being under estimated due to the abundance of low stellar mass galaxies in
this redshift bin.

We recognise that our assumption of a “constant efficiency” model is at odds
with Scoville et al. (2017), who argue that whilst cold molecular gas
reservoirs increase with $z$ (as $(1+z)^{1.84}$), the star-formation rate
increases more rapidly (as $(1+z)^{2.9}$), indicating that the peak of star
formation is a consequence of both increased molecular gas content in galaxies
and higher star-formation efficiency. We also note that at $z\gtrsim 1$ early-
type galaxies have been shown to be more compact for a given stellar mass than
their local counterparts (e.g., Daddi et al., 2005; Cappellari et al., 2009),
which taken in combination with the “Kennicutt-Schmidt” relation (a power-law
relation between star-formation rate and gas surface densities, Kennicutt,
1998; Schmidt, 1959) implies that star formation may be more efficient at
$z\gtrsim 1$.

The $3$D stacking approach we use derives the average properties for galaxies
in our sample as a function of redshift, and thus we do not measure the
molecular gas mass and star-formation rates for individual sources. Whilst the
UDS DR11 catalog does include Mstellar estimates (which are evaluated at the
peak maximum likelihood redshift) for individual galaxies, our $3$D stacking
method bins galaxies according to the discretized redshift probability
distribution ($\mathcal{P}(z)$), and as such each galaxy in our sample
effectively contributes to the flux in all redshift intervals. Hence, using
this $3$D stacking method precludes a Mstellar selection relative to our
redshift bins. Therefore, we are not able to repeat the analysis from Scoville
et al. (2017) to test their assertion of an evolving star-formation
efficiency.

###  5.5 The epoch of molecular gas

Although we cannot quantify the contribution of higher star-formation
efficiencies to the peak of star-formation rate density at $z\approx 2$, the
symmetry between our “constant efficiency” models with our statistically
derived $\rho_{\rm H_{2}}(z)$ indicates a star formation history which is
predominantly driven by an increased supply of molecular gas in galaxies,
rather than a significant evolution in star-formation efficiency (consistent
with the findings of Decarli et al., 2020; Magnelli et al., 2020). With this
in mind we now turn to the formation of H2 itself.

Cazaux & Spaans (2004) combine a microscopic model for the relative rates of
gas-phase and dust H2 production with a cosmological model to show the more
efficient dust-phase production becomes the dominant route to H2 formation at
$z\approx 3$–$6$ for reasonable assumptions about the conditions of the
interstellar medium of early galaxies. Therefore, there is a perfect storm for
massive galaxy growth at $z\approx 2$: not only is the cosmic accretion rate
at its peak, massive halos have had time to grow, galaxies have increased gas
densities, and previous generations of stars in the progenitors of these
systems have provided the metal enrichment that accelerates the formation of
H2, which, as the fuel for star formation, drives galaxy growth; this could be
described as the epoch of molecular gas.

###  5.6 Estimating the evolution of cosmic molecular gas mass density at
$z\gtrsim 2.5$

S16 intentionally restrict their calibration sample to galaxies at $z\leq 3$
to ensure observed $850\,\mu$m emission is from the rest-frame long wavelength
RJ tail, where dust is optically thin and emission is dominated by the
contribution of cold dust (which is well represented by a mass-weighted
$T_{\rm dust}=25$ K). In Figure $3$ we have shown an average of our $\rho_{\rm
H_{2}}(z)$ estimates at $2.5\lesssim z\lesssim 6$ (the UDS redshifts are
untested at earlier epochs as there are no UDS galaxies with spectroscopic
redshifts at $z\gtrsim 6.5$), but note that at these redshifts estimates are
less reliable due to large uncertainties in the RJ correction (see equation
10) as rest-frame emission approaches the peak of the SED.

In the optically thick regime (as rest-frame dust emission moves off the long-
wavelength RJ tail) the rest-frame emission no longer correlates with the
total dust mass of a galaxy and probes only the surface dust, which using the
approach of S16 would result in under-estimation of $L_{\rm 850}$ and hence
the molecular gas mass. However, as the rest-frame emission approaches the
peak of the SED we are increasingly sensitive to the dense, warm dust
component, which significantly boosts the luminosity (with only a small mass
fraction) and dominates the emission close to the SED peak. Consequently rest-
frame dust emission at high-$z$ is not well represented by a mass-weighted
$T_{\rm dust}=25$ K, which would result in a over-estimate of the dust and gas
mass.

In addition to these competing effects, we are also likely to be missing a
significant population of lower mass galaxies at $z>2.5$. As shown in Figure 4
the $95\%$ stellar mass completeness at $z\approx 2.5$ is predicted to be
$\simeq 10^{9.5}\,{\rm M_{\odot}}$, so we are are simply not sensitive to the
majority of low mass galaxies at the highest redshifts. In addition, although
relatively rare (with number counts $\mathrm{N}(>3.5\,\mathrm{mJy)\simeq 3000\
deg}^{-2}$; Geach et al., 2017) SMGs are dust-rich ($M_{\rm dust}\sim
10^{9}\,{\rm M_{\odot}}$; e.g., da Cunha et al., 2015; Magnelli et al., 2019)
and about $20\%$ are undetected in optical/near-infrared surveys (e.g.,
Dudzevičiūtė et al., 2019). This non-detection of SMGs is unlikely to have a
significant impact on our estimates at low $z$. However, at $z\gtrsim 2.5$
since we are significantly under-sampling the galaxy population and as the
number of galaxies in our redshift bins fall the non-detection of dust-rich
SMGs becomes more statistically significant, further contributing to an under-
estimation of the molecular gas mass density at $z\gtrsim 2.5$.

The overall impact of the above is difficult to quantify. However, as shown by
Figure 3 our results at $z\gtrsim 2.5$ are systematically lower than the
estimates obtained via direct CO line emission, which suggests that the use of
this method past $z\approx 2.5$ (when $\lambda_{obs}=850\,\mu$m no longer
probes the rest-frame RJ tail) results in an under-estimation of the molecular
gas mass density. In consequence, whilst our results are highly uncertain at
$z>2.5$ we suggest that to $z\lesssim 6$ these can be seen as providing a
lower-limit to the molecular gas mass density.

##  6 CONCLUSIONS

We employ a $3$-dimensional stacking method (Viero et al., 2013) and an
empirically calibrated RJ luminosity-to-gas mass ratio (S16) to derive the
average molecular gas mass as a function of redshift utilising a sample of
${\approx}150,000$ galaxies in the UKIDSS-UDS field. By combining these
techniques we are able to reduce the statistical uncertainties on the
evolution of the molecular gas mass density, $\rho_{\rm H_{2}}(z)$, within the
redshift range $0.5\lesssim z\lesssim 2.5$. We find that:

• $\rho_{\rm H_{2}}(z)$ shows a clear evolution over cosmic time which traces
that of the star-formation rate density with a peak $\approx 2\times
10^{7}\,{\rm M_{\odot}}\,{\rm Mpc^{-3}}$ at $z\approx 2$.

• Our results are consistent with those of blank field CO line surveys, albeit
our estimates are systematically lower than those derived using observations
of higher excitation CO lines. This may in part be a consequence of the line
ratios used to translate higher excitation CO line luminosity to ground state
CO line luminosity.

• Our results are an order of magnitude lower than those derived by Liu et al.
(2019) who use the Hughes et al. (2017) luminosity-to-gas mass calibration to
estimate molecular gas masses for the ALMA continuum detected galaxies in
their sample. This difference in results may be in part due to selection
effects, as their ALMA-selected sample preferentially selects dust-rich (and
consequently gas-rich), sources, whereas by using a NIR selection we are
likely to miss $\approx{20}\%$ of these dust-rich SMGs.

• $\rho_{\rm H_{2}}(z)$ can be broadly modelled by inverting the star-
formation rate density (Wilkins et al., 2019) with a constant (Geach &
Papadopoulos, 2012) or weakly evolving (Tacconi et al., 2018) volume averaged
star-formation efficiency. Our “constant efficiency” models closely align to
our statistically derived $\rho_{\rm H_{2}}(z)$.

• $\rho_{\rm H_{2}}(z)$ can be derived from first principles from the halo
mass function (Murray et al., 2013) in conjunction with stellar-halo mass
(Moster et al., 2013) and ISM-stellar mass ratios (Scoville et al., 2017). To
obtain this estimate we make the assumption that all galaxies are star-forming
and hence this can be seen as an upper limit for $\rho_{\rm H_{2}}(z)$ with
respect to the stellar mass range sampled.

We have demonstrated that by applying a statistical method and the approach of
Scoville et al. (2016) we can provide robust, statistically significant
constraints to the cosmological gas mass density to $z\lesssim 2.5$. Our
results show an evolution that mirrors that of the star-formation rate density
indicating that the peak of the star formation history is primarily driven by
an increased supply of molecular gas rather than a significantly increased
star-formation efficiency. We have shown that at $z\gtrsim 2.5$ we detect dust
emission from high mass galaxies, even with our near-infrared selected sample.
Hence, in the future there is potential for this approach to be extended to
provide improved constraints at higher-$z$ through $1$ mm/$3$ mm wide-field
surveys with facilities such as the Large Millimeter Telescope.

We thank the anonymous referee for their helpful comments. TKG acknowledges
support from a UK Science and Technology Facilities Council studentship. KEKC
is supported by a Royal Society Leverhulme Senior Research Fellowship
(SRF/R1/191013). KEKC and MF acknowledge the support from STFC (grant number
ST/R000905/1). JEG is supported by a Royal Society University Research
Fellowship. CCL acknowledges support from the Royal Society under grant
RGF/EA/181016. MPK acknowledges support from the First TEAM grant of the
Foundation for Polish Science No. POIR.04.04.00-00-5D21/18-00. MPK also
acknowledges support from the Polish National Agency for Academic Exchange
under Grant No. PPI/APM/2018/1/00036/U/001. The James Clerk Maxwell Telescope
has historically been operated by the Joint Astronomy Centre on behalf of the
Science and Technology Facilities Council of the United Kingdom, the National
Research Council of Canada and the Netherlands Organisation for Scientific
Research. The S2CLS map data were taken as part of Program ID MJLSC02.
Additional funds for the construction of SCUBA-2 were provided by the Canada
Foundation for Innovation. The UKIDSS project is defined in Lawrence et al.
(2007). Further details on the UDS can be found in Almaini et al. (in prep.).
UKIDSS uses the UKIRT Wide Field Camera (WFCAM; Casali et al. 2007). The
photometric system is described in Hewett et al. (2006), and the calibration
is described in Hodgkin et al. (2009). The pipeline processing and science
archive are described in Irwin et al (in prep.) and Hambly et al. (2008).

## References

  * Alberts et al. (2014) Alberts, S., Pope, A., Brodwin, M., et al. 2014, MNRAS, 437, 437, doi: 10.1093/mnras/stt1897
  * Aravena et al. (2016a) Aravena, M., Decarli, R., Walter, F., et al. 2016a, ApJ, 833, 71, doi: 10.3847/1538-4357/833/1/71
  * Aravena et al. (2016b) —. 2016b, ApJ, 833, 68, doi: 10.3847/1538-4357/833/1/68
  * Ashby et al. (2013) Ashby, M. L. N., Willner, S. P., Fazio, G. G., et al. 2013, ApJ, 769, 80, doi: 10.1088/0004-637X/769/1/80
  * Behrens et al. (2018) Behrens, C., Pallottini, A., Ferrara, A., Gallerani, S., & Vallini, L. 2018, MNRAS, 477, 552, doi: 10.1093/mnras/sty552
  * Berta et al. (2013) Berta, S., Lutz, D., Nordon, R., et al. 2013, A&A, 555, L8, doi: 10.1051/0004-6361/201321776
  * Birnboim & Dekel (2003) Birnboim, Y., & Dekel, A. 2003, MNRAS, 345, 349, doi: 10.1046/j.1365-8711.2003.06955.x
  * Blain et al. (2002) Blain, A. W., Smail, I., Ivison, R. J., Kneib, J.-P., & Frayer, D. T. 2002, Phys. Rep., 369, 111, doi: 10.1016/S0370-1573(02)00134-5
  * Bolatto et al. (2013) Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013, ARA&A, 51, 207, doi: 10.1146/annurev-astro-082812-140944
  * Boogaard et al. (2019) Boogaard, L. A., Decarli, R., González-López, J., et al. 2019, ApJ, 882, 140, doi: 10.3847/1538-4357/ab3102
  * Bothwell et al. (2013) Bothwell, M. S., Smail, I., Chapman, S. C., et al. 2013, MNRAS, 429, 3047, doi: 10.1093/mnras/sts562
  * Bouwens et al. (2016) Bouwens, R. J., Aravena, M., Decarli, R., et al. 2016, ApJ, 833, 72, doi: 10.3847/1538-4357/833/1/72
  * Brammer et al. (2008) Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, The Astrophysical Journal, 686, 1503, doi: 10.1086/591786
  * Cappellari et al. (2009) Cappellari, M., di Serego Alighieri, S., Cimatti, A., et al. 2009, ApJ, 704, L34, doi: 10.1088/0004-637X/704/1/L34
  * Carilli & Walter (2013) Carilli, C. L., & Walter, F. 2013, ARA&A, 51, 105, doi: 10.1146/annurev-astro-082812-140953
  * Carilli et al. (2016) Carilli, C. L., Chluba, J., Decarli, R., et al. 2016, ApJ, 833, 73, doi: 10.3847/1538-4357/833/1/73
  * Casali et al. (2007) Casali, M., Adamson, A., Alves de Oliveira, C., et al. 2007, A&A, 467, 777, doi: 10.1051/0004-6361:20066514
  * Cazaux & Spaans (2004) Cazaux, S., & Spaans, M. 2004, ApJ, 611, 40, doi: 10.1086/422087
  * Chabrier (2003) Chabrier, G. 2003, Publications of the Astronomical Society of the Pacific, 115, 763, doi: 10.1086/376392
  * Chapman et al. (2005) Chapman, S. C., Blain, A. W., Smail, I., & Ivison, R. J. 2005, ApJ, 622, 772, doi: 10.1086/428082
  * Chary & Pope (2010) Chary, R.-R., & Pope, A. 2010, arXiv e-prints, arXiv:1003.1731.  https://arxiv.org/abs/1003.1731
  * Coppin et al. (2007) Coppin, K. E. K., Swinbank, A. M., Neri, R., et al. 2007, ApJ, 665, 936, doi: 10.1086/519789
  * da Cunha et al. (2013) da Cunha, E., Groves, B., Walter, F., et al. 2013, ApJ, 766, 13, doi: 10.1088/0004-637X/766/1/13
  * da Cunha et al. (2015) da Cunha, E., Walter, F., Smail, I. R., et al. 2015, ApJ, 806, 110, doi: 10.1088/0004-637X/806/1/110
  * Daddi et al. (2005) Daddi, E., Renzini, A., Pirzkal, N., et al. 2005, ApJ, 626, 680, doi: 10.1086/430104
  * Dahlen et al. (2013) Dahlen, T., Mobasher, B., Faber, S. M., et al. 2013, ApJ, 775, 93, doi: 10.1088/0004-637X/775/2/93
  * Decarli et al. (2014) Decarli, R., Walter, F., Carilli, C., et al. 2014, ApJ, 782, 78, doi: 10.1088/0004-637X/782/2/78
  * Decarli et al. (2016a) Decarli, R., Walter, F., Aravena, M., et al. 2016a, ApJ, 833, 70, doi: 10.3847/1538-4357/833/1/70
  * Decarli et al. (2016b) —. 2016b, ApJ, 833, 69, doi: 10.3847/1538-4357/833/1/69
  * Decarli et al. (2019) Decarli, R., Walter, F., Gónzalez-López, J., et al. 2019, arXiv e-prints, arXiv:1903.09164.  https://arxiv.org/abs/1903.09164
  * Decarli et al. (2020) Decarli, R., Aravena, M., Boogaard, L., et al. 2020, The Astrophysical Journal, 902, 110, doi: 10.3847/1538-4357/abaa3b
  * Dudzevičiūtė et al. (2019) Dudzevičiūtė, U., Smail, I., Swinbank, A. M., et al. 2019, arXiv e-prints, arXiv:1910.07524.  https://arxiv.org/abs/1910.07524
  * Eales et al. (2012) Eales, S., Smith, M. W. L., Auld, R., et al. 2012, ApJ, 761, 168, doi: 10.1088/0004-637X/761/2/168
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * Frayer et al. (1998) Frayer, D. T., Ivison, R. J., Scoville, N. Z., et al. 1998, ApJ, 506, L7, doi: 10.1086/311639
  * Frayer et al. (1999) —. 1999, ApJ, 514, L13, doi: 10.1086/311940
  * Furusawa et al. (2008) Furusawa, H., Kosugi, G., Akiyama, M., et al. 2008, ApJS, 176, 1, doi: 10.1086/527321
  * Geach & Papadopoulos (2012) Geach, J. E., & Papadopoulos, P. P. 2012, ApJ, 757, 156, doi: 10.1088/0004-637X/757/2/156
  * Geach et al. (2017) Geach, J. E., Dunlop, J. S., Halpern, M., et al. 2017, MNRAS, 465, 1789, doi: 10.1093/mnras/stw2721
  * Gómez-Guijarro et al. (2019) Gómez-Guijarro, C., Riechers, D. A., Pavesi, R., et al. 2019, ApJ, 872, 117, doi: 10.3847/1538-4357/ab002a
  * González-López et al. (2019) González-López, J., Decarli, R., Pavesi, R., et al. 2019, arXiv e-prints, arXiv:1903.09161.  https://arxiv.org/abs/1903.09161
  * Gould & Salpeter (1963) Gould, R. J., & Salpeter, E. E. 1963, ApJ, 138, 393, doi: 10.1086/147654
  * Groves et al. (2015) Groves, B. A., Schinnerer, E., Leroy, A., et al. 2015, ApJ, 799, 96, doi: 10.1088/0004-637X/799/1/96
  * Hambly et al. (2008) Hambly, N. C., Collins, R. S., Cross, N. J. G., et al. 2008, MNRAS, 384, 637, doi: 10.1111/j.1365-2966.2007.12700.x
  * Hewett et al. (2006) Hewett, P. C., Warren, S. J., Leggett, S. K., & Hodgkin, S. T. 2006, Monthly Notices of the Royal Astronomical Society, 367, 454, doi: 10.1111/j.1365-2966.2005.09969.x
  * Hodgkin et al. (2009) Hodgkin, S. T., Irwin, M. J., Hewett, P. C., & Warren, S. J. 2009, Monthly Notices of the Royal Astronomical Society, 394, 675, doi: 10.1111/j.1365-2966.2008.14387.x
  * Hughes et al. (2017) Hughes, T. M., Ibar, E., Villanueva, V., et al. 2017, MNRAS, 468, L103, doi: 10.1093/mnrasl/slx033
  * Ivison et al. (2011) Ivison, R. J., Papadopoulos, P. P., Smail, I., et al. 2011, MNRAS, 412, 1913, doi: 10.1111/j.1365-2966.2010.18028.x
  * Jarvis et al. (2013) Jarvis, M. J., Bonfield, D. G., Bruce, V. A., et al. 2013, Monthly Notices of the Royal Astronomical Society, 428, 1281, doi: 10.1093/mnras/sts118
  * Kaasinen et al. (2019) Kaasinen, M., Scoville, N., Walter, F., et al. 2019, ApJ, 880, 15, doi: 10.3847/1538-4357/ab253b
  * Kennicutt (1998) Kennicutt, Robert C., J. 1998, ApJ, 498, 541, doi: 10.1086/305588
  * Lawrence et al. (2007) Lawrence, A., Warren, S. J., Almaini, O., et al. 2007, MNRAS, 379, 1599, doi: 10.1111/j.1365-2966.2007.12040.x
  * Lenkić et al. (2019) Lenkić, L., Bolatto, A. D., Förster Schreiber, N. M., et al. 2019, arXiv e-prints, arXiv:1908.01791.  https://arxiv.org/abs/1908.01791
  * Liang et al. (2018) Liang, L., Feldmann, R., Faucher-Giguère, C.-A., et al. 2018, MNRAS, 478, L83, doi: 10.1093/mnrasl/sly071
  * Liang et al. (2019) Liang, L., Feldmann, R., Kereš, D., et al. 2019, MNRAS, 489, 1397, doi: 10.1093/mnras/stz2134
  * Liu et al. (2019) Liu, D., Schinnerer, E., Groves, B., et al. 2019, arXiv e-prints, arXiv:1910.12883.  https://arxiv.org/abs/1910.12883
  * Lonsdale et al. (2003) Lonsdale, C. J., Smith, H. E., Rowan-Robinson, M., et al. 2003, PASP, 115, 897, doi: 10.1086/376850
  * Madau & Dickinson (2014) Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415, doi: 10.1146/annurev-astro-081811-125615
  * Magdis et al. (2012) Magdis, G. E., Daddi, E., Béthermin, M., et al. 2012, ApJ, 760, 6, doi: 10.1088/0004-637X/760/1/6
  * Magnelli et al. (2019) Magnelli, B., Karim, A., Staguhn, J., et al. 2019, ApJ, 877, 45, doi: 10.3847/1538-4357/ab1912
  * Magnelli et al. (2020) Magnelli, B., Boogaard, L., Decarli, R., et al. 2020, ApJ, 892, 66, doi: 10.3847/1538-4357/ab7897
  * Miettinen et al. (2017) Miettinen, O., Delvecchio, I., Smolčić, V., et al. 2017, A&A, 606, A17, doi: 10.1051/0004-6361/201730762
  * Millard et al. (2020) Millard, J. S., Eales, S. A., Smith, M. W. L., et al. 2020, Monthly Notices of the Royal Astronomical Society, 494, 293–315, doi: 10.1093/mnras/staa609
  * Moster et al. (2013) Moster, B. P., Naab, T., & White, S. D. M. 2013, MNRAS, 428, 3121, doi: 10.1093/mnras/sts261
  * Murray et al. (2013) Murray, S. G., Power, C., & Robotham, A. S. G. 2013, Astronomy and Computing, 3, 23, doi: 10.1016/j.ascom.2013.11.001
  * Nordon et al. (2013) Nordon, R., Lutz, D., Saintonge, A., et al. 2013, ApJ, 762, 125, doi: 10.1088/0004-637X/762/2/125
  * Oteo et al. (2018) Oteo, I., Ivison, R. J., Dunne, L., et al. 2018, ApJ, 856, 72, doi: 10.3847/1538-4357/aaa1f1
  * Papadopoulos & Geach (2012) Papadopoulos, P. P., & Geach, J. E. 2012, ApJ, 757, 157, doi: 10.1088/0004-637X/757/2/157
  * Pavesi et al. (2018) Pavesi, R., Sharon, C. E., Riechers, D. A., et al. 2018, ApJ, 864, 49, doi: 10.3847/1538-4357/aacb79
  * Planck Collaboration et al. (2011) Planck Collaboration, Abergel, A., Ade, P. A. R., et al. 2011, A&A, 536, A25, doi: 10.1051/0004-6361/201116483
  * Planck Collaboration et al. (2016) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A13, doi: 10.1051/0004-6361/201525830
  * Popping et al. (2019) Popping, G., Pillepich, A., Somerville, R. S., et al. 2019, ApJ, 882, 137, doi: 10.3847/1538-4357/ab30f2
  * Pozzetti et al. (2010) Pozzetti, L., Bolzonella, M., Zucca, E., et al. 2010, A&A, 523, A13, doi: 10.1051/0004-6361/200913020
  * Privon et al. (2018) Privon, G. C., Narayanan, D., & Davé, R. 2018, ApJ, 867, 102, doi: 10.3847/1538-4357/aae485
  * Riechers et al. (2013) Riechers, D. A., Bradford, C. M., Clements, D. L., et al. 2013, Nature, 496, 329, doi: 10.1038/nature12050
  * Riechers et al. (2019) Riechers, D. A., Pavesi, R., Sharon, C. E., et al. 2019, ApJ, 872, 7, doi: 10.3847/1538-4357/aafc27
  * Riechers et al. (2020) Riechers, D. A., Boogaard, L. A., Decarli, R., et al. 2020, The Astrophysical Journal, 896, L21, doi: 10.3847/2041-8213/ab9595
  * Schmidt (1959) Schmidt, M. 1959, ApJ, 129, 243, doi: 10.1086/146614
  * Scoville et al. (2014) Scoville, N., Aussel, H., Sheth, K., et al. 2014, ApJ, 783, 84, doi: 10.1088/0004-637X/783/2/84
  * Scoville et al. (2016) Scoville, N., Sheth, K., Aussel, H., et al. 2016, The Astrophysical Journal, 820, 83, doi: 10.3847/0004-637x/820/2/83
  * Scoville et al. (2017) Scoville, N., Lee, N., Bout, P. V., et al. 2017, The Astrophysical Journal, 837, 150, doi: 10.3847/1538-4357/aa61a0
  * Scoville (2013) Scoville, N. Z. 2013, Evolution of star formation and gas, ed. J. Falcón-Barroso & J. H. Knapen, 491
  * Simpson et al. (2014) Simpson, J. M., Swinbank, A. M., Smail, I., et al. 2014, ApJ, 788, 125, doi: 10.1088/0004-637X/788/2/125
  * Solomon & Vanden Bout (2005) Solomon, P. M., & Vanden Bout, P. A. 2005, ARA&A, 43, 677, doi: 10.1146/annurev.astro.43.051804.102221
  * Stach et al. (2017) Stach, S. M., Swinbank, A. M., Smail, I., et al. 2017, ApJ, 849, 154, doi: 10.3847/1538-4357/aa93f6
  * Tacconi et al. (2010) Tacconi, L. J., Genzel, R., Neri, R., et al. 2010, Nature, 463, 781, doi: 10.1038/nature08773
  * Tacconi et al. (2013) Tacconi, L. J., Neri, R., Genzel, R., et al. 2013, ApJ, 768, 74, doi: 10.1088/0004-637X/768/1/74
  * Tacconi et al. (2018) Tacconi, L. J., Genzel, R., Saintonge, A., et al. 2018, ApJ, 853, 179, doi: 10.3847/1538-4357/aaa4b4
  * Thomson et al. (2012) Thomson, A. P., Ivison, R. J., Smail, I., et al. 2012, MNRAS, 425, 2203, doi: 10.1111/j.1365-2966.2012.21584.x
  * Tremonti et al. (2004) Tremonti, C. A., Heckman, T. M., Kauffmann, G., et al. 2004, ApJ, 613, 898, doi: 10.1086/423264
  * Tukey (1958) Tukey, J. W. 1958, Ann. Math. Statist, 29, 614
  * Viero et al. (2013) Viero, M. P., Moncelsi, L., Quadri, R. F., et al. 2013, ApJ, 779, 32, doi: 10.1088/0004-637X/779/1/32
  * Wakelam et al. (2017) Wakelam, V., Bron, E., Cazaux, S., et al. 2017, Molecular Astrophysics, 9, 1, doi: 10.1016/j.molap.2017.11.001
  * Walter et al. (2014) Walter, F., Decarli, R., Sargent, M., et al. 2014, ApJ, 782, 79, doi: 10.1088/0004-637X/782/2/79
  * Walter et al. (2016) Walter, F., Decarli, R., Aravena, M., et al. 2016, ApJ, 833, 67, doi: 10.3847/1538-4357/833/1/67
  * Wilkins et al. (2019) Wilkins, S. M., Lovell, C. C., & Stanway, E. R. 2019, MNRAS, 490, 5359, doi: 10.1093/mnras/stz2894
  * Zavala et al. (2018) Zavala, J. A., Aretxaga, I., Dunlop, J. S., et al. 2018, MNRAS, 475, 5585, doi: 10.1093/mnras/sty217
  * Zavala et al. (2021) Zavala, J. A., Casey, C. M., Manning, S. M., et al. 2021, The Evolution of the IR Luminosity Function and Dust-obscured Star Formation in the Last 13 Billion Years.  https://arxiv.org/abs/2101.04734


* [100] M. Honda, M. Sajjad Athar, T. Kajita, K. Kasahara, and S. Midorikawa Phys. Rev. D 92 no. 2, (2015) 023004, arXiv:1502.03916 [astro-ph.HE].
  * [101] J. A. Formaggio and G. P. Zeller Rev. Mod. Phys. 84 (2012) 1307–1341, arXiv:1305.7513 [hep-ex].
  * [102] IceCube Collaboration, M. G. Aartsen et al. Eur. Phys. J. C 80 no. 1, (2020) 9, arXiv:1902.07771 [hep-ex].
  * [103] KM3Net Collaboration, S. Adrian-Martinez et al. J. Phys. G 43 no. 8, (2016) 084001, arXiv:1601.07459 [astro-ph.IM].
  * [104] KM3NeT Collaboration, L. Nauta et al. PoS ICRC2021 (2021) 1123.
  * [105] A. Aurisano, “Recent Results from MINOS and MINOS+,” June, 2018.  https://zenodo.org/record/1286760.
  * [106] NOvA Collaboration, M. A. Acero et al. Phys. Rev. Lett. 123 no. 15, (2019) 151803, arXiv:1906.04907 [hep-ex].
  * [107] KM3NeT Collaboration, S. Aiello et al. Eur. Phys. J. C 82 no. 1, (2022) 26, arXiv:2103.09885 [hep-ex].
  * [108] G. Cowan, K. Cranmer, E. Gross, and O. Vitells Eur. Phys. J. C 71 (2011) 1554, arXiv:1007.1727 [physics.data-an]. [Erratum: Eur.Phys.J.C 73, 2501 (2013)].
  * [109] IceCube-Gen2 Collaboration, M. G. Aartsen et al. Phys. Rev. D 101 no. 3, (2020) 032006, arXiv:1911.06745 [hep-ex].
  * [110] Hyper-Kamiokande Collaboration, K. Abe et al. arXiv:1805.04163 [physics.ins-det].
  * [111] DUNE Collaboration, A. Abud Abed et al. Phys. Rev. D 105 no. 7, (2022) 072006, arXiv:2109.01304 [hep-ex].
  * [112] H. Nunokawa, S. J. Parke, and R. Zukanovich Funchal Phys. Rev. D 72 (2005) 013009, arXiv:hep-ph/0503283.
  * [113] L. Fields, A. D. Marino, J. P. Ochoa-Ricoux, and J. Spitz arXiv:2209.07480 [hep-ex].
  * [114] S. P. Mikheyev and A. Yu. Smirnov Sov. J. Nucl. Phys. 42 (1985) 913–917. [,305(1986)].
  * [115] H. A. Bethe Phys. Rev. Lett. 56 (1986) 1305.
  * [116] S. P. Rosen and J. M. Gelb Phys. Rev. D 34 (1986) 969.
  * [117] G. J. Feldman and R. D. Cousins Phys. Rev. D 57 (1998) 3873–3889, arXiv:physics/9711021.
  * [118] A. Cabrera et al. arXiv:2008.11280 [hep-ph].
  * [119] KM3NeT, JUNO Collaboration, S. Aiello et al. arXiv:2108.06293 [hep-ex].
  * [120] M. Blennow and T. Schwetz JHEP 09 (2013) 089, arXiv:1306.3988 [hep-ph].
  * [121] S. T. Petcov and M. Piai Phys. Lett. B 533 (2002) 94–106, arXiv:hep-ph/0112074.
  * [122] DUNE Collaboration, A. A. Abud et al. in 2022 Snowmass Summer Study.  3, 2022.  arXiv:2203.06100 [hep-ex].
  * [123] A. Lowe, K. Majumdar, K. Mavrokoridis, B. Philippou, A. Roberts, C. Touramanis, and J. Vann Instruments 4 no. 4, (2020) 35, arXiv:2011.02292 [physics.ins-det].
  * [124] S. Para et al. in 2022 Snowmass Summer Study.  3, 2022.  arXiv:2203.07501 [hep-ex].
  * [125] D. A. Dwyer et al. JINST 13 no. 10, (2018) P10007, arXiv:1808.02969 [physics.ins-det].
  * [126] S. Kubota et al. arXiv:2203.12109 [hep-ex].
  * [127] A. Mastbaum, F. Psihas, and J. Zennamo arXiv:2203.14700 [physics.ins-det].
  * [128] Theia Collaboration, M. Askins et al. in 2022 Snowmass Summer Study.  2, 2022.  arXiv:2202.12839 [hep-ex].
  * [129] M. Ball et al.
  * [130] J. Eldred, S. Nagaitsev, V. Shiltsev, A. Valishev, R. Zwaska, and M. Syphers in 2022 Snowmass Summer Study.  3, 2022.  arXiv:2203.08276 [physics.acc-ph].
  * [131] S. Nagaitsev and V. Lebedev arXiv:2111.06932 [physics.acc-ph].
  * [132] R. Ainsworth et al. arXiv:2106.02133 [physics.acc-ph].
  * [133] A. M. Suliga arXiv:2207.09632 [astro-ph.HE].
  * [134] Hyper-Kamiokande Collaboration, K. Abe et al. PTEP 2018 no. 6, (2018) 063C01, arXiv:1611.06118 [hep-ex].
  * [135] P. F. De Salas, S. Gariazzo, O. Mena, C. A. Ternes, and M. Tórtola Front. Astron. Space Sci. 5 (2018) 36, arXiv:1806.11051 [hep-ph].
  * [136] K. Møller, A. M. Suliga, I. Tamborra, and P. B. Denton JCAP 05 (2018) 066, arXiv:1804.03157 [astro-ph.HE].
  * [137] IceCube Collaboration, M. G. Aartsen et al. Science 342 (2013) 1242856, arXiv:1311.5238 [astro-ph.HE].
  * [138] P. B. Denton, D. Marfatia, and T. J. Weiler JCAP 08 (2017) 033, arXiv:1703.09721 [astro-ph.HE].
  * [139] IceCube Collaboration, M. G. Aartsen et al. Astrophys. J. 849 no. 1, (2017) 67, arXiv:1707.03416 [astro-ph.HE].
  * [140] IceCube Collaboration, M. G. Aartsen et al. Astrophys. J. 809 no. 1, (2015) 98, arXiv:1507.03991 [astro-ph.HE].
  * [141] IceCube Collaboration, J. Stachurska EPJ Web Conf. 207 (2019) 02005, arXiv:1905.04237 [hep-ex].
  * [142] IceCube-Gen2 Collaboration, M. G. Aartsen et al. J. Phys. G 48 no. 6, (2021) 060501, arXiv:2008.04323 [astro-ph.HE].
  * [143] KATRIN Collaboration, M. Aker et al. Phys. Rev. Lett. 123 no. 22, (2019) 221802, arXiv:1909.06048 [hep-ex].
  * [144] L. Gastaldo et al. Eur. Phys. J. ST 226 no. 8, (2017) 1623–1694.
  * [145] B. Alpert et al. Eur. Phys. J. C 75 no. 3, (2015) 112, arXiv:1412.5060 [physics.ins-det].
  * [146] Project 8 Collaboration, A. Ashtari Esfahani et al. J. Phys. G 44 no. 5, (2017) 054004, arXiv:1703.02037 [physics.ins-det].
  * [147] E. Di Valentino, S. Gariazzo, and O. Mena Phys. Rev. D 104 no. 8, (2021) 083504, arXiv:2106.15267 [astro-ph.CO].
  * [148] S. Gariazzo et al. arXiv:2205.02195 [hep-ph].
  * [149] R. X. Adhikari et al. arXiv:2209.11726 [hep-ph].
  * [150] KamLAND-Zen Collaboration, S. Abe et al. arXiv:2203.02139 [hep-ex].
  * [151] A. J. Long, C. Lunardini, and E. Sabancilar JCAP 08 (2014) 038, arXiv:1405.7654 [hep-ph].
  * [152] E. Roulet and F. Vissani JCAP 10 (2018) 049, arXiv:1810.00505 [hep-ph].
  * [153] ESSnuSB Collaboration, A. Alekou et al. Eur. Phys. J. C 81 no. 12, (2021) 1130, arXiv:2107.07585 [hep-ex].
  * [154] A. Alekou et al. in 2022 Snowmass Summer Study.  3, 2022.  arXiv:2203.08803 [physics.acc-ph].
  * [155] H. Abele et al. arXiv:2211.10396 [physics.ins-det].
  * [156] MEMPHYS Collaboration, L. Agostino, M. Buizza-Avanzini, M. Dracos, D. Duchesneau, M. Marafini, M. Mezzetto, L. Mosca, T. Patzak, A. Tonazzo, and N. Vassilopoulos JCAP 01 (2013) 024, arXiv:1206.6665 [hep-ex].
  * [157] ICAL Collaboration, S. Ahmed et al. Pramana 88 no. 5, (2017) 79, arXiv:1505.07380 [physics.ins-det].
  * [158] S. Sahoo, A. Kumar, and S. K. Agarwalla arXiv:2110.13207 [hep-ph].
  * [159] A. Bogacz et al. in 2022 Snowmass Summer Study.  3, 2022.  arXiv:2203.08094 [hep-ph].
  * [160] J. De Blas et al. in 2022 Snowmass Summer Study.  3, 2022.  arXiv:2203.07261 [hep-ph].
  * [161] C. Aimè et al. in 2022 Snowmass Summer Study.  3, 2022.  arXiv:2203.07256 [hep-ph].
  * [162] nuSTORM Collaboration, L. A. Ruso et al. in 2022 Snowmass Summer Study.  3, 2022.  arXiv:2203.07545 [hep-ex].
  * [163] S. Vihonen in 22nd International Workshop on Neutrinos from Accelerators.  2, 2022.  arXiv:2202.05509 [hep-ex].
  * [164] J. Cao et al. Phys. Rev. ST Accel. Beams 17 (2014) 090101, arXiv:1401.8125 [physics.acc-ph].
  * [165] J. Tang, S. Vihonen, and T.-C. Wang JHEP 12 (2019) 130, arXiv:1909.01548 [hep-ph].
  * [166] A. V. Akindinov et al. arXiv:1902.06083 [physics.ins-det].
  * [167] M. Perrin-Terrin in Very Large Volume Neutrino Telescope Workshop.  2021\.  https://indico.ific.uv.es/event/3965/contributions/14787/.
  * [168] M. Perrin-Terrin Eur. Phys. J. C 82 no. 5, (2022) 465, arXiv:2112.12848 [hep-ex].

##  Appendix A LOIs tagged in NF01, “Neutrino Oscillations”

Neutrino experiments and facilities
---
1 |  Ongoing Science Program of Super-Kamiokande
2 |  Legacy of the Daya Bay Reactor Neutrino Experiment
3 |  Reactor and Geo Neutrinos at SNO+
4 |  T2K Experiment: future plans and capabilities
5 |  The NOvA Physics Program through 2025
6 |  Expected Final Sensitivity of the NOvA Experiment to 3-Flavor Neutrino Oscillations
7 |  Development of a joint oscillation analysis by the NOvA and T2K collaborations
8 |  The JUNO Experiment
9 |  Long-Baseline Physics in DUNE
10 |  Atmospheric $\nu_{\tau}$ Appearance in the Deep Underground Neutrino Experiment
11 |  Neutrino oscillations with IceCube-DeepCore and the IceCube Upgrade
12 |  Neutrino Detection and Ranging
13 |  Neutrino physics with muon-decay medium-baseline neutrino beam facility (MOMENT)
Simulation and reconstruction
14 |  Fast Simulations for Noble Liquid Experiments
15 |  Differentiable Simulators for HEP
16 |  The Future of Machine Learning in Rare Event Searches
17 |  $\nu_{\tau}$ Reconstruction in the Deep Underground Neutrino Experiment
18 |  Scalable, End-to-End Optimizable Data Reconstruction and Physics Inference Techniques for Large-scale Particle Imaging Neutrino Detectors
Theory and phenomenology
19 |  Direct Probes of the Matter Effect in Neutrino Oscillations
20 |  Role of higher order maxima of oscillation probabilities at long baseline neutrino experiments
21 |  Role of high energy beam tunes in optimizing the sensitivity to current unknowns at DUNE
22 |  Leptonic Sum Rules
23 |  Tau Neutrino Physics
24 |  Physics with Sub-GeV Atmospheric Neutrinos
25 |  Joint Experimental Oscillation Analyses in Search of Sterile Neutrinos
26 |  Three Sterile Neutrinos in E6
27 |  Ultralight dark matter and neutrinos
28 |  Neutrino Non-Standard Interactions
29 |  Neutrino Minimal Standard Model — a unified theory of microscopic and cosmic scales
Supporting program (theory and experiment)
30 |  Neutrino-induced Shallow- and Deep-Inelastic Scattering
31 |  Electron scattering and neutrino programs
32 |  Nucleon Form Factors for Neutrino Physics
33 |  Physics Opportunities in ANNIE

##  Appendix B Snowmass whitepapers related to NF01

1 |  Snowmass Neutrino Frontier: DUNE Physics Summary
---|---
2 |  A Gaseous Argon-Based Near Detector to Enhance the Physics Capabilities of DUNE
3 |  Electron Scattering and Neutrino Physics
4 |  High Energy Physics Opportunities Using Reactor Antineutrinos
5 |  Design Considerations for Fermilab Multi-MW Proton Facility in the DUNE/LBNF era
6 |  A Cost-Effective Upgrade Path for the Fermilab Accelerator Complex
7 |  An Upgrade Path for the Fermilab Accelerator Complex
8 |  The Physics Case for a Neutrino Factory
9 |  Theia: Summary of physics program
10 |  Physics Opportunities for the Fermilab Booster Replacement
11 |  Japan’s Strategy for Future Projects in High Energy Physics
12 |  Tau Neutrinos in the Next Decade: from GeV to EeV