# The ALMA REBELS Survey: The First Infrared Luminosity Function Measurement
at $\mathbf{z\sim 7}$

L. Barrufet1, P. A. Oesch1,2, R. Bouwens3, H. Inami4, L. Sommovigo5, H.
Algera4,6, E. da Cunha7<EMAIL_ADDRESS>M. Aravena8, P.
Dayal9, A. Ferrara5, Y. Fudamoto10,11, V. Gonzalez12,13 L. Graziani14,15 A.
Hygate3 I. de Looze16,17, T. Nanayakkara18, A. Pallottini5, R. Schneider14,15
M. Stefanon19,20 M. Topping21 and P. van Der Werf3
1Department of Astronomy, University of Geneva, Chemin Pegasi 51, 1290
Versoix, Switzerland
2Cosmic Dawn Center (DAWN), Niels Bohr Institute, University of Copenhagen,
Jagtvej 128, København N, DK-2200, Denmark
3Leiden Observatory, Leiden University, PO Box 9500, 2300 RA Leiden, The
Netherlands
4 Hiroshima Astrophysical Science Center, Hiroshima University, 1-3-1,
Kagamiyama-Higashi, Hiroshima 739-8526, Japan
5 Scuola Normale Superiore, Piazza dei Cavalieri, 7, 56126 Pisa, Italy
6 National Astronomical Observatory of Japan, 2-21-1, Osawa, Mitaka, Tokyo,
Japan
7 International Centre for Radio Astronomy Research, University of Western
Australia, Stirling Hwy, Crawley, 26WA, 6009, Australia
8 Nucleo de Astronomia Facultad de Ingenieria y Ciencias, Universidad Diego
Portales, Av Ejercito 441, Santiago, Chile
9 Kapteyn Astronomical Institute, University of Groningen, P.O. Box 800, 9700,
AV Groningen, The Netherlands
10 Waseda Research Institute for Science and Engineering, Faculty of Science
and Engineering, Waseda University, 3-4-1 Okubo, Shinjuku, Tokyo 169-8555,
Japan
11National Astronomical Observatory of Japan, 2-21-1, Osawa, Mitaka, Tokyo,
Japan
12 Departmento de Astronomia, Universidad de Chile, Casilla 36-D, Santiago
7591245, Chile
13Centro de Astrofisica y Tecnologias Afines (CATA), Camino del Observatorio
1515, Las Condes, Santiago 7591245, Chile
14 Dipartimento di Fisica, Sapienza, Universita di Roma, Piazzale Aldo Moro 5,
I-00185 Roma, Italy
15 INAF/Osservatorio Astronomico di Roma, via Frascati 33, I-00078 Monte
Porzio Catone, Roma, Italy
16 Sterrenkundig Observatorium, Ghent University, Krijgslaan 281-59, 9000,
Gent, Belgium
17 Dept. of Physics & Astronomy, University College London, Gower Street,
London WC1E 6BT, United Kingdom
18 Centre for Astrophysics & Supercomputing, Swinburne University of
Technology, PO Box 218, Hawthorn, VIC 3112, Australia
19 Departament d’Astronomia i Astrofisica, Universitat de Valencia, C. Dr.
Moliner 50, E-46100 Burjassot, Valencia, Spain
cia, C. Dr. Moliner 50, E-46100 Burjassot, Valencia, Spain
20 Unidad Asociada CSIC "Grupo de Astrofisica Extragalactica y Cosmologia"
(Instituto de Fisica de Cantabria - Universitat de Valencia)
21 Steward Observatory University of Arizona, 933 N Cherry Ave, Tucson, AZ
85721, United States

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We present the first observational infrared luminosity function (IRLF)
measurement in the Epoch of Reionization (EoR) based on a UV-selected galaxy
sample with ALMA spectroscopic observations. Our analysis is based on the ALMA
large program Reionization Era Bright Emission Line Survey (REBELS), which
targets 42 galaxies at $\mathrm{z=6.4-7.7}$ with [CII] 158$\micron$ line
scans. 16 sources exhibit a dust detection, 15 of which are also
spectroscopically confirmed through the [CII] line. The IR luminosities of the
sample range from $\log L_{IR}/L_{\odot}=11.4$ to 12.2. Using the UVLF as a
proxy to derive the effective volume for each of our target sources, we derive
IRLF estimates, both for detections and for the full sample including IR
luminosity upper limits. The resulting IRLFs are well reproduced by a
Schechter function with the characteristic luminosity of $\log
L_{*}/L_{\odot}=11.6^{+0.2}_{-0.1}$. Our observational results are in broad
agreement with the average of predicted IRLFs from simulations at $z\sim 7$.
Conversely, our IRLFs lie significantly below lower redshift estimates,
suggesting a rapid evolution from $z\sim 4$ to $z\sim 7$, into the
reionization epoch. The inferred obscured contribution to the cosmic star-
formation rate density at $z\sim 7$ amounts to
$\mathrm{log(SFRD/M_{\odot}/yr/Mpc^{3})=-2.66^{+0.17}_{-0.14}}$ which is at
least $\sim$10% of UV-based estimates. We conclude that the presence of dust
is already abundant in the EoR and discuss the possibility of unveiling larger
samples of dusty galaxies with future ALMA and JWST observations.

###### keywords:

Galaxies: high-redshift, luminosity function. Infrared: galaxies

††pubyear: 2022††pagerange: The ALMA REBELS Survey: The First Infrared
Luminosity Function Measurement at $\mathbf{z\sim 7}$–The ALMA REBELS Survey:
The First Infrared Luminosity Function Measurement at $\mathbf{z\sim 7}$

##  1 Introduction

It is still a crucial open question in astrophysics when the first galaxies
formed and how they built up their mass. The continuous discovery of higher
redshift galaxies is pushing the boundaries of our knowledge of galaxy
evolution (e.g., Dunlop, 2013; Stark, 2016; Dayal & Ferrara, 2018; Schaerer et
al., 2022; Naidu et al., 2022b; Atek et al., 2022; Adams et al., 2022). In
particular, the discovery of a significant population of luminous and massive
galaxies at $\mathrm{z>9}$ has posed questions about the speed of early
stellar mass production (e.g. Oesch et al., 2016; Laporte et al., 2021; Naidu
et al., 2022a; Labbe et al., 2022).

Until recently, the knowledge of galaxies at $\mathrm{z>7}$ was mainly based
on rest-frame ultraviolet (UV) observations (Oesch et al., 2018a; Bouwens et
al., 2021a). These samples might not be complete, however, as they might miss
extremely dust obscured, but highly star-forming galaxies (e.g., Casey et al.,
2019).

From an observational point of view, the Atacama Large Millimeter Array (ALMA)
is the most powerful tool to study dust at high redshift (e.g., Capak et al.,
2015; Bouwens, 2016; Bowler et al., 2018; Béthermin et al., 2020). However,
the cost to obtain statistical samples of galaxies in the EoR results in the
fact that only a modest number of galaxies have been characterized in detail
so far (e.g., Watson et al., 2015; Smit et al., 2018; Laporte et al., 2019;
Faisst et al., 2020; Harikane et al., 2021; Schouws et al., 2022).
Furthermore, the study of dust at $\mathrm{2<z<6}$ was for a long time limited
to bright dusty galaxies such as submillimetre galaxies (SMGs; e.g., Gruppioni
et al., 2013; Wang et al., 2019b; Barrufet et al., 2020). However, ALMA is
bridging the gap between these extreme dusty massive galaxies and more
moderate star-forming galaxies (see Hodge & da Cunha 2020 for a review).

The recent observational improvements have allowed the discovery of the
emergence of high-z dusty galaxies at $\mathrm{z>6}$. In particular, Fudamoto
et al. (2021) has serendipitously detected two dusty galaxies at
$\mathrm{z_{spec}\sim 7}$ near massive neighbors at the same redshifts. This
shows that dusty galaxies in the EoR could be more common than previously
thought, which leads to the question of whether the number of dusty galaxies
at $\mathrm{z>6}$ is higher than expected (see also Barrufet et al., 2022;
Nelson et al., 2022; Rodighiero et al., 2022).

The possible underestimation of the number of dusty galaxies would have a
direct impact on the obscured Star Formation Rate Density (SFRD), which
remains uncertain at $\mathrm{z>3}$ (Casey et al., 2019). Several studies have
calculated the obscured SFRD at $\mathrm{z>5}$ based on serendipitous sources
resulting in largely differing conclusions (e.g Gruppioni et al., 2020;
Fudamoto et al., 2021; Talia et al., 2021; Casey et al., 2021; Viero et al.,
2022). While some studies find that 2mm selected, dusty galaxies contribute
$\mathrm{\sim 30\%}$ to the integrated star-formation rate density between
$\mathrm{3<z<6}$ (Casey et al., 2021), others report a significantly larger
obscured SFRD that remains constant over redshift (e.g., Gruppioni et al.,
2020; Talia et al., 2021). An approach to clarify the contribution of dust-
obscured star formation to the cosmic star formation history is to measure the
infrared luminosity function (IRLF) all the way into the EoR. The shape and
scale of the IRLF are crucial to understanding the abundance of dusty galaxies
and how rapidly dust is formed in the early universe. This directly affects
the fraction of star formation that is obscured in forming galaxies, and
thereby the formation (or rise) of metals.

Due to the wealth of rest-frame UV observations, the UV luminosity function
(UVLF) is well constrained up to $\mathrm{z\sim 9}$ (e.g., Bouwens et al.,
2007; Bouwens et al., 2015; Oesch et al., 2018b; Bowler et al., 2020; Bouwens
et al., 2021b), and we even have some information at $\mathrm{z\sim 9-10}$
(Oesch et al., 2018a; Harikane et al., 2022) and beyond now with JWST (e.g.
Naidu et al., 2022a; Donnan et al., 2022; Atek et al., 2022; Adams et al.,
2022; Finkelstein et al., 2022). In contrast, the IRLF is still quite
uncertain at high redshifts. Current measurements of the IRLF rely on small
numbers of dusty sources at $\mathrm{z>3.5}$ (e.g., Wang et al., 2019a;
Gruppioni et al., 2020). This leads to large uncertainties in the IRLF
parameters, including the faint-end slopes, and disagreements between
different survey results (e.g., Gruppioni et al., 2013; Koprowski et al.,
2017; Lim et al., 2020; Popping et al., 2020; Gruppioni et al., 2020).

The recent study of Zavala et al. (2021) compiled the results of several
surveys and combined those with semi-empirical modelling to constrain the
evolution of the IRLF out to $\mathrm{z>5}$, albeit with significant
uncertainties. However, an IRLF at $\mathrm{z\sim 7}$ has not been measured
directly using dust continuum observations yet. In this context, we use the
data from the Reionization Era Bright Emission Line Survey (REBELS), an ALMA
large program aimed at obtaining a statistical sample of normal star-forming
galaxies at $\mathrm{z>6.4}$ (see Bouwens et al. 2022 for details). REBELS has
increased the number of spectroscopically observed massive galaxies in the EoR
by a factor $\mathrm{\times\sim 4-5}$ compared to the previous literature
(Bouwens et al., 2021a). The same strategy of the REBELS selection was tested
in a pilot program presented in Schouws et al. (2022). This study showed the
potential of ALMA as a high redshift ‘machine’ and the six pilot galaxies are
also included in the main REBELS sample (Smit et al., 2018; Schouws et al.,
2021, 2022). While observations from the REBELS program were just recently
completed and analysis of the full data set now underway, its data have
already been used for a number of scientific analyses, including the discovery
of serendipitous dust-obscured sources at $z\sim 7$ (Fudamoto et al., 2021),
modelling the dust and ISM properties of $z>6$ galaxies (e.g., Sommovigo et
al., 2022; Dayal et al., 2022; Ferrara et al., 2022), measuring their detailed
specific SFRs (Topping et al., 2022), calculating their SFRD Algera et al.
(2022), estimating Ly$\alpha$ transmission around luminous sources in
overdense $z\sim 7$ environments (Endsley et al., 2022), and constraining the
neutral gas fraction out to the EoR (Heintz et al., 2022).

In this paper, we use this survey to calculate – for the first time – an IRLF
at $\mathrm{z\sim 7}$. In Section 2, we describe the ALMA observations and the
infrared luminosity calculations used in this work. The methodology for
calculating the IRLF and their values is described in Section 3. We present
the results on the obscured SFRD of REBELS galaxies in Section 4. We discuss
our results in Section 5 and present a summary and our conclusions in Section
6.

##  2 REBELS observations

###  2.1 ALMA observations and catalogue

In this work, we use data from REBELS (Bouwens et al., 2021a) which is a Cycle
7 ALMA large program of $\mathrm{\sim 40}$ UV bright galaxies at
$\mathrm{z>6.4}$. The selection was based on UV brightness
($\mathrm{-23<M_{UV}<-21.3}$) and photometric redshifts for galaxies
identified over a combined area of $\mathrm{\sim 7deg^{2}}$ in several fields
(see Bouwens et al. 2021a for details). This survey of spectral scan
observations identifies bright ISM cooling lines ([CII], [OIII]) while
simultaneously probing the dust-continuum in bands $\mathrm{158\ \mu m}$ and
$\mathrm{88\ \mu m}$, respectively, which is essential to derive the infrared
luminosity ($\mathrm{L_{IR}}$). Given its selection, the REBELS sample only
spans a limited range in redshift and UV luminosities. Even though it is UV
selected, the sample is representative of massive star-forming galaxies at
$\mathrm{z\sim 7}$, providing an extensive probe of ISM reservoirs in the EoR
(Bouwens et al., 2022; Ferrara et al., 2022).

In this work, we only focus on galaxies that were scanned for [CII], i.e.,
sources with $\mathrm{z_{phot}=6.4-7.7}$. The total sample used in this study
contains 42 galaxies with [CII] scanned, 16 of which with a dust continuum
detection at more than $3\sigma$. Notably, 15 of these 16 sources also do have
a significant [CII] emission line detection and thus a robust spectroscopic
redshift measurement (Inami et al., 2022).

###  2.2 Infrared luminosity from REBELS survey

In this section, we describe the infrared luminosity measurements from Inami
et al. (2022) and the average properties of the REBELS galaxies.

When deriving the infrared luminosities of our sample, we have to make an
assumption about the dust temperature. Estimating this based on a few
photometric detections in the far-infrared is very challenging. Sommovigo et
al. (2021) solve this difficulty using $\mathrm{L_{[CII]}}$ as a proxy for the
dust mass and the underlying continuum to constrain the dust temperature. This
is particularly useful for the REBELS survey, given that [CII] estimates (or
upper limits) are available for the full sample. Using these measurements,
Sommovigo et al. (2022) find an average dust temperature of
$\mathrm{T_{d}=46K}$ for the REBELS sample. Hence, Inami et al. (2022) assumed
a Spectral Energy Distribution (SED) with dust temperature and emissivity from
Sommovigo et al. (2022) ($\mathrm{T_{d}=46K}$ and $\mathrm{\beta=2}$
respectively) to calculate the infrared luminosity based on the ALMA dust
continuum flux. For the galaxies without dust continuum detection a 3$\sigma$
upper limit was derived both for the continuum flux and the corresponding
infrared luminosity. A cosmic microwave background correction was applied for
all galaxies, with and without dust detection. The correction depends on the
exact redshift, but lies in the range of $\mathrm{8-14\%}$ (see Inami et al.
2022 for details).

Figure 1:  Infrared luminosity against UV absolute magnitude with the redshift
colour-coded for the REBELS (filled symbols) and ALPINE (empty symbols)
samples for both $\mathrm{3\sigma}$ detections (dots) and upper limits
(triangles). The REBELS sample does not show significant differences between
detections and upper limits. $\mathrm{L_{IR}}$ does not depend on
$\mathrm{M_{UV}}$ or redshift. The small $\mathrm{L_{IR}}$ dynamic range and
the flatness are comparable with the ALPINE sample at $\mathrm{4.5<z<6}$
although ALPINE extends to fainter UV galaxies (empty triangles and dots for
upper limits and detections respectively). The ALPINE relation presented in
Khusanova et al. (2021) is shown in the black dashed line.

Using the derived IR luminosity measurements, we plot in Figure 1 the relation
between UV and IR-luminosities. Given the selection of UV luminous sources,
the dynamic range both in UV and IR luminosities is limited. The REBELS sample
only probes the most massive, UV-luminous galaxies at these redshifts. It is
composed of luminous infrared galaxies (LIRGs;
$\mathrm{10^{11}<L_{IR}/L_{\odot}<10^{12}}$) except for REBELS-25, the
brightest galaxy in our sample with $\mathrm{log(L_{IR})\sim 12.2L_{\odot}}$
(see Hygate et al. 2022 for details). The fact that we found only one ultra
luminous infrared galaxy (ULIRG; $\mathrm{L_{IR}>10^{12}L_{\odot}}$) in the
REBELS sample could be due to the UV bright selection of REBELS galaxies with
$\mathrm{-23<M_{UV}<-21.3}$. We discuss this further in a later section.

We compare the IR luminosities from REBELS with the sample from the ALMA Large
Program to INvestigate [CII] at Early times (ALPINE, Le Fèvre et al., 2020)
which targets UV-selected sources at lower redshifts at $\mathrm{4.5<z<6}$.
The ALPINE sample spans a wider $\mathrm{M_{UV}}$ range
($\mathrm{-23.3<M_{UV}<-20}$) but is also mostly composed of LIRGs (see Figure
1) finding also in general dusty galaxies (Pozzi et al., 2021, Sommovigo et
al. 2022b in prep). Our REBELS sample shows that UV-selected galaxies at
$\mathrm{z\sim 7}$ have comparable infrared luminosities to UV-selected
galaxies at lower redshift ($\mathrm{4.5<z<6}$) (see Section 5 for
Discussion).

##  3 Infrared luminosity function at z $\sim$ 7

In this section, we explain the procedure to calculate the luminosity function
(LF). The main complication in computing a luminosity function using a
targeted survey such as REBELS is that it is not straightforward to derive a
selection volume for each source. This can be overcome by basing our volume
estimates on the UV luminosity function as a proxy, as was successfully
demonstrated in Yan et al. (2020) who used the ALPINE UV targeted sample to
derive the [CII] luminosity function. Here, we closely follow their approach.

###  3.1 Calculation of the luminosity function

Our derivation is based on the $\mathrm{z\sim 7}$ UVLF from Bouwens et al.
(2021a). This is used to derive a representative volume for the UV-selected
sources. In practice, we use the UVLF to compute the number of expected
galaxies in bins of UV luminosity assuming a volume-limited survey over the
full selection area of the REBELS sample of 7 deg2 and $\mathrm{z=6.4-7.7}$
(see Fig. 2). This is given by:

$\mathrm{N_{exp}=\phi_{UV}(M_{UV})\,\Delta M_{UV}\,V_{tot}}$ (1)

where $\mathrm{\phi_{UV}(M)}$ is the UVLF from Bouwens et al. (2021a) per
magnitude bin $\mathrm{\Delta M_{UV}}$, and $\mathrm{V_{tot}}$ is the total
survey volume over which REBELS sources were selected. REBELS only targets a
very small sub-sample of all galaxies expected in such a large survey. We can
compute a correction factor to account for this sampling incompleteness in
each UV luminosity bin as $\mathrm{f_{UV}=N_{exp}/N_{obs}}$, where
$\mathrm{N_{obs}}$ is the number of targeted REBELS galaxies in each
$\mathrm{M_{UV}}$ bin.

Figure 2:  Number of $\mathrm{L_{IR}}$ detections against the UV absolute
magnitude. The histogram shows the detected sources in red and the non-
detections in grey with the fraction of detections/total indicated in the
lower numbers. Also shown is the UVLF from (Bouwens et al., 2021a) as a dashed
line. This is used to compute the representative volume for each of our
targets. The small numbers above the LF indicate how many galaxies are
expected per MUV bin in a volume-limited survey spanning the REBELS target
selection area of 7 deg2. Clearly, REBELS only targets a very small fraction
of the full galaxy population at faint UV luminosities, which we account for
in our analysis (see main text).

While the correction factor above is derived for a volume-limited survey, the
requirement of a dust continuum detection can further introduce a reduction in
the survey volume for each source. Namely, it can limit the maximum redshift
up to which a given source would remain detected. This is accounted for by
computing the so-called maximum comoving volume $\mathrm{V_{max,i}}$ for each
galaxy i (see Schmidt, 1968). Specifically,
$\mathrm{V_{max,i}}=\int_{z_{min}}^{z_{max,i}}\nicefrac{{d^{2}V}}{{dz\,d\Omega}}\,\Omega\,dz$,
where $\mathrm{z_{max,i}}$ is either the upper edge of the redshift bin of the
LF, or, if smaller, the maximum redshift up to which source i would remain
continuum detected at $>3\sigma$. $\Omega$ is the survey volume. In practice,
$\mathrm{z_{max,i}}=7.7$ for most galaxies, except for the faintest few
sources in the sample.

We now have all quantities to calculate the IR luminosity function
$\mathrm{\phi_{IR}}$ in bins of $\mathrm{L{IR}}$. This is given by:

$\mathrm{\phi_{IR}(\log L_{IR})=\frac{1}{\Delta logL_{IR}}\sum_{i\in
bin}\frac{f_{UV,i}}{V_{max,i}}}$ (2)

where i runs over all sources in a given IR luminosity bin $\mathrm{\log
L_{IR}\pm\Delta\log L_{IR}/2}$ (see Eq. 3 in Yan et al. 2020). The
uncertainties on the IRLF bins are computed as the Poisson errors in each
$\mathrm{L_{IR}}$ bin.

Note that this calculation is independent of the assumed survey area $\Omega$,
since both $\mathrm{{V_{max}}}$ and $\mathrm{{f_{UV}}}$ are directly
proportional to it.

We repeat the above calculation twice. In the first case, we only consider
continuum detected galaxies (16 sources); in the second case, we include the
full REBELS sample (42 sources), treating non-detections as upper-limits. The
completeness factors $\mathrm{f_{UV}}$ are computed separately for both cases.
The resulting IRLFs are in very good agreement, as discussed in the next
section.

###  3.2 The infrared luminosity function at $z\sim 7$

Figure 3:  Infrared luminosity function at $\mathrm{z\sim 7}$ for the REBELS
sample (red dots and lines) compared with simulations (dashed lines) and
observations (solid lines). The IRLF was calculated both only using the
galaxies with dust continuum detections (16 galaxies, empty dots) as well as
using the full sample including upper limits (42 galaxies, filled red dots).
The red line shows the Schechter 1976 fit for the total sample. The shaded
area shows the uncertainty of the luminosity function Schechter function fit
with the total sample which is larger at the low luminosity end due to the
lack of data. The rest of the lines show both theoretical and observational
IRLF studies in several fields. Our study is in agreement with Li et al. in
prep (dark purple line) which predicts a similar number of dusty galaxies in a
broad range of luminosities. The dark grey line is the IRLF at $\mathrm{z\sim
7}$ from Zavala et al. (2021) and predicts a larger number of galaxies than
our study for the bright end with luminosities
($\mathrm{12.5<log(L_{IR}/L_{\odot}<13}$) whereas our luminosity function does
not predict a significant number of galaxies at $\mathrm{z\sim 7}$ with
$\mathrm{log(L_{IR}/L_{\odot}>12.5}$. TNG simulations at $\mathrm{z\sim 6}$
from Shen et al. (2021) show a systematic shift with respect to our fitting,
but consistent in shape (blue dashed line). Dayal et al. (2022) and Lagos et
al. (2020) simulations at $\mathrm{z\sim 7}$ (light blue and grey line
respectively) present a 1 dex difference in the lower luminosity with our
result in between them. The yellow line and dots indicate the IRLF at
$\mathrm{z\sim 5.25}$ predicted by the serendipitous galaxies found in the
ALPINE survey presented in Gruppioni et al. (2020), whereas the orange symbols
show Wang et al. (2019a) results at similar redshift.

####  3.2.1 The Step-Wise IRLF

In this section, we first present the step-wise LF by using the methodology
described in the previous subsection, before we derive parametric Schechter
function fits. Figure 3 shows the resulting LFs in three equidistant
luminosity bins $\mathrm{log(L_{IR}/L_{\odot}}$: [11.3-11.6], [11.6-11.9] and
[11.9-12.2], both for our detections-only and our full sample. The derived
stepwise LFs are in excellent agreement, showing that the detection-only
sample is not biased significantly. In the rest of the paper, we use the total
sample as a baseline.

For the detection-only sample, we further test the possible impact of
uncertainties in the IR luminosity estimates. Specifically, we use a Monte
Carlo technique in which we perturb the initial $L_{IR}$ measurements by their
statistical (Gaussian) uncertainties 10,000 times and rederive the IRLF in
each case. We then use the median and 16th and 84th percentiles, respectively,
as the uncertainties. We do not find significant differences in the resulting
LF values, but the uncertainties are increased as can also be seen in Figure
3.

####  3.2.2 Schechter Function Fits

We now derive a parametric estimate of the IRLF based on the classic Schechter
function from Schechter 1976, commonly used both in the local and the high-z
Universe (Johnston, 2011). The three parameters that define the Schechter
function are $\mathrm{\phi^{*}}$, $\mathrm{L_{*}}$ and $\mathrm{\alpha}$; the
normalization factor of the overall density of galaxies, the characteristic
luminosity, and the faint-end luminosity slope, respectively. Due to the lack
of data at low $\mathrm{L_{IR}}$, we have restricted $\mathrm{\alpha}$ taking
into account the faint-end slope values found in the literature (see Section 5
for details). We fix the slope to $\mathrm{\alpha=-1.3}$ in our fitting, which
is the value derived for the ALPINE high-z IRLF in Gruppioni et al. (2020).

We use a Bayesian Monte Carlo Markov Chain (MCMC) approach to derive the
posterior distribution of the Schechter function parameters. Hence, we compute
the $\mathrm{\phi_{IR}}$, $\mathrm{L_{*}}$, while keeping the slope fixed at
$\mathrm{\alpha=-1.3}$. We have set these initial parameters centered at the
values obtained by minimizing the error function first
($\mathrm{log(\phi_{IR})=-3.5}$, $\mathrm{log(L^{*})=11.7}$), and then use
non-informative Gaussian priors. We then perform 20,000 MCMC iterations and
ensure that these are converged. We find that posterior distribution of the
parameters is similar in both cases, either including the total sample
(considering upper limits) or only detections. Therefore, we only present the
Schechter function with uncertainties for total sample in Figure 3. The
$\mathrm{1\sigma}$ uncertainty of the fit function was also calculated from
the MCMC chains computing the 16th and 84th percentiles of the posterior
distributions. The $\mathrm{\phi_{IR}}$ uncertainties in the fainter end are
$\mathrm{\sim 0.5\ dex}$, while at the brighter end they are $\mathrm{<0.2\
dex}$. The IRLF is best constrained between
$\mathrm{11.5<log(L_{IR}/L_{\odot})<12}$, and shows that the density of
sources drops quickly ($\mathrm{log(\phi_{IR})<-6.5dex^{-1}Mpc^{-3}}$) at
luminosities above $\mathrm{log(L_{IR}/L_{\odot})>12.3}$.

The resulting Schechter function parameters are
$\mathrm{log(\phi_{IR})=-4.38^{+0.38}_{-0.35}dex^{-1}Mpc^{-3}}$ and
$\mathrm{log(L_{*}/L_{\odot})=11.60^{+0.23}_{-0.13}}$ with a fixed
$\mathrm{\alpha=-1.3}$ (see Table 1 for the summary of the main parameters).
Our analysis shows a $\mathrm{z\sim 7}$ IRLF with a considerable number of
LIRGs that drops in the ULIRG range suggesting a limit in luminosity at
$\mathrm{log(L_{IR}/L_{\odot})\sim 12.3}$. This is in general agreement with
some theoretical studies. The IRLF at $\mathrm{L_{IR}<11.5L_{\odot}}$ is
uncertain and a larger study with fainter galaxies should be carried out to
accurately measure the IRLF at the fainter luminosity end.

We compared our results to both theoretical and observational IRLF studies at
similar redshifts (see dashed and continuous lines respectively in Figure 3).
Generally, our results are in broad agreement with some simulated IRLFs at
similar redshift. When comparing to lower redshift observations at $z\sim
5-6$, however, we find that our IRLF is more than an order of magnitude lower.
Finally, our IRLF shows an interesting evolution with redshift, compared with
the literature, not only in number density (as was previously shown in
Koprowski et al. (2020); Fujimoto et al. (2023)), but also in
$\mathrm{L_{*}}$. This could be due to our UV-selected sample being biased to
bright sources and further study with a similar selection at different
redshift should be carried out to confirm the possibility of evolution with
$\mathrm{L_{*}}$.

We discuss the points above in more detail in Section 5. We also discuss in
subsection 5.3 the importance that our data is UV-bright selected which cannot
take into account extremely dust-obscured sources that are faint in the UV.

##  4 Obscured star formation rate density

In this section, we calculate the obscured SFRD directly through the IRLF
derived in the previous section. We calculate the SFRD in two different ways:
1) by simply summing up the step-wise infrared densities for the data in the
REBELS sample and 2) by integrating the Schechter IRLF over the luminosity
range $\mathrm{10.5<log(L_{IR}/L_{\odot})<13}$. These limits were selected in
the range over which we can define the Schechter function. Note that the
integration limits are narrow but, due to the luminosity bins, there is no
data to constrain a lower-limit integration. Further analysis is produced in
section 5. In both cases we use a conversation factor
$\mathrm{\kappa=10^{-10}M_{\odot}/yr/L_{\odot}}$.

For the step-wise estimates, we considered both the total sample and
detections. We find $\mathrm{log(SFRD/(M_{\odot}/yr/Mpc^{3}))=-3.21\pm 0.18}$
taking only into account the dust continuum detections, which is slightly
lower than for the total sample with
$\mathrm{log(SFRD/M_{\odot}/yr/Mpc^{3})=-2.93\pm 0.20}$. This SFRD estimate
needs to be considered as a lower limit, since it only takes into account the
three luminosity bins.

To extrapolate to fainter luminosities, we have calculated the SFRD for the
Schechter LFs. In particular, we use the MCMC chains to derive the median
posterior SFRD and the associated uncertainties. We find
$\mathrm{log(SFRD/M_{\odot}/yr/Mpc^{3})=-2.66^{+0.17}_{-0.14}}$ where the
uncertainties correspond to the 16-84th percentile (see Figure 4 ). As
expected, this SFRD is larger than the SFRD calculated from the observations,
since it is integrated over the full luminosity range (
$\mathrm{10.5<log(L_{IR}/L_{\odot})<13}$). Notice that REBELS is a UV-selected
sample and the obscured SFRD needs to be taken into account as a robust lower
limit (see caveats in Section 5.3). Finally, the SFRD was computed adding the
serendipitous sources from the REBELS sample presented in Fudamoto et al.
(2021). The sum of the two points, UV-selected galaxies and serendipitous
’dark’ systems, is
$\mathrm{log(SFRD/(M_{\odot}/yr/Mpc^{3}))=-2.53^{+0.17}_{-0.14}}$.

$\mathrm{\alpha}$ | $\mathrm{log(L^{*})}$ | $\mathrm{log(\phi_{IR})}$ | $\mathrm{log(SFRD)}$
---|---|---|---
| $\mathrm{[L_{\odot}]}$ | $\mathrm{[dex^{-1}Mpc^{-3}]}$ | $\mathrm{[M_{\odot}/yr/Mpc^{3}]}$
Schechter Function Fit
-1.3 (fix) | $\mathrm{11.60}^{+0.23}_{-0.13}$ | $\mathrm{-4.38^{+0.38}_{-0.35}}$ | $\mathrm{-2.66^{+0.17}_{-0.14}}$
Total sample | $\mathrm{11.15}$ | $\mathrm{-4.3^{+0.1}_{-0.1}}$ | $\mathrm{-2.93\pm 0.20}$
| $\mathrm{11.75}$ | $\mathrm{-4.6^{+0.2}_{-0.2}}$ |
| $\mathrm{12.05}$ | $\mathrm{-5.5^{+0.4}_{-0.5}}$ |
Detections | $\mathrm{11.15}$ | $\mathrm{-4.4^{+0.2}_{-0.2}}$ | $\mathrm{-3.21\pm 0.18}$
| $\mathrm{11.75}$ | $\mathrm{-4.6^{+0.3}_{-0.3}}$ |
| $\mathrm{12.05}$ | $\mathrm{-5.1^{+0.2}_{-0.5}}$ |
Table 1: Summary of the main parameters of this study. The first column shows
the faint luminosity slope ($\mathrm{\alpha}$), and the second column shows
the luminosity function at the determined luminosity bin (third column).
Finally, the fourth column shows the obscured star formation rate density
taking into account the three luminosity bins. The first row shows the best
fit Schechter function parameters for a fixed slope of $\mathrm{\alpha=-1.3}$,
while the subsequent rows show the total sample and only with detections.

We compare our results with previous studies in the literature for both
similar samples to REBELS and other dusty galaxies at high redshift. Our
derived obscured SFRD of the REBELS sample is $\mathrm{\sim 13\pm 1\%}$ of the
total CSFRD at $\mathrm{z\sim 7}$ from Madau & Dickinson (2014) and
$\mathrm{9\%}$ of the unobscured SFRD estimate from Bouwens et al. (2022).
This is in agreement with the range of obscured SFRD predictions of Zavala et
al. (2021), who use a compilation of several surveys to derive a model of the
IRLF evolution. Our resulting obscured SFRD lies in the upper part of their
inferred SFRD range being the first result at $\mathrm{z\sim 7}$ calculated
through [CII] spectroscopic scans. In an accompanying paper, Algera et al.
(2022) also derived the SFRD for the REBELS sample using the stellar mass as a
proxy to calculate the SFRD through a stacking analysis. While our best
estimates are a factor $\sim 2.5\times$ lower, the measurements are consistent
within the $\mathrm{1\sigma}$ uncertainties.

Figure 4:  Star formation rate density against redshift for the REBELS sample
at $\mathrm{z\sim 7}$ and several works in the literature. The black line
shows the total SFRD from Madau & Dickinson (2014) whereas the orange shaded
region shows the obscured SFRD (Zavala et al., 2021). Our results show a
moderate SFRD calculated from the fitted IRLF (red triangle) which increases
if the two serendipitous normal dusty REBELS galaxies from Fudamoto et al.
(2021) are taken into account (orange dot). Similarly, Algera et al. (2022)
obtains a larger contribution to the obscured star formation but in agreement
within $\mathrm{1\sigma}$ error (dark orange dot). DSFGs from the ALMA 2 mm
photometric blind survey show a decrease in SFRD over redshift (purple
squares; Casey et al., 2021). The 1.3 mm ALMA blind survey presented in Dunlop
et al. (2017) shows a obscured SFRD at $\mathrm{1<z<4.5}$ that decreases at
$\mathrm{z>2}$ (purple diamonds). Khusanova et al. (2021) shows the SFRD from
the ALPINE survey at $\mathrm{z\sim 5}$ (brown area). Also from ALPINE,
Gruppioni et al., 2020 present a larger obscured SFRD which is decreasing at
$\mathrm{z>3}$ (pink area) with the last redshift bin at $\mathrm{z>4}$
containing only one source (dashed pink area). Similarly, Wang et al. (2019a)
shows a decreasing SFRD (light purple area) with large uncertainty in the last
bin at $\mathrm{z\sim 4}$ (dashed light purple area). Koprowski et al. (2020)
presented a constrained SFRD up to $\mathrm{z\sim 4}$ (purple area). REBELS
results shows the presence of dust at $\mathrm{z\sim 7}$ even in UV-selected
galaxies.

In Figure 4 we also present the obscured SFRD for several studies showing the
lack of consensus at $\mathrm{z>3}$ on the obscured SFRD. Our SFRD result is
comparable to DSFGs selected at 2mm from Casey et al. (2021), who reports a
decrease in the obscured SFRD over $\mathrm{4<z<6}$. In contrast to these
findings, the SFRD from serendipitous sources found in the ALPINE survey
present a non-evolving SFRD across the whole redshift range of the sample
($\mathrm{1<z<5.5}$). Their calculated SFRD is over two orders of magnitude
more than our results at $\mathrm{z\sim 7}$. Similarly, longer wavelength
studies support a flatter evolution of the SFRD at $\mathrm{3<z<6}$, albeit
with more moderate SFRD (Talia et al., 2021). In contrast, our results show
lower SFRD at $\mathrm{z\sim 7}$, which, when compared to literature at lower
redshifts, supports a non-flat SFRD across redshift (see section 5 for
discussion).

##  5 Discussion

In this section, we compare our IRLF results with observational and
theoretical studies. However, due to the underlying assumptions, IRLFs from
simulations are not directly comparable. As a result, our findings broadly
concur with theoretical research. On the observational side, the literature
shows a large range of IRLF suggesting SFRD discrepancies of $\mathrm{\sim 2}$
orders of magnitude. We also explore the causes for the different results in
the literature and compare to our IRLF and SFRD.

###  5.1 Comparison to Literature

Some theoretical IRLFs at $\mathrm{z\sim 6-7}$ agree quite well with our
findings. For example, Li et al. in prep. show a similar IRLF over the
luminosity range $\mathrm{10.5<log(L_{IR}/L_{\odot})<12.5}$, as do the TNG+300
simulations shown in Shen et al. (2021). But throughout the whole infrared
luminosity range, the latter exhibits larger number densities by $\mathrm{\sim
0.5dex}$. A plausible explanation for this shift is the difference in redshift
($\mathrm{\Delta z\sim 1}$) between our results and those of Shen et al.
(2021), as the IRLF is expected to decrease in number density at increasing
redshift (see e.g. (Koprowski et al., 2017; Fujimoto et al., 2023)).

Our results contrast with those from Lagos et al. (2020) which themselves
differ by $\mathrm{\sim 0.5dex}$ despite the fact that both utilise semi-
analytical models based on merger trees. Over the full range of our directly
observed luminosities ($\mathrm{log(L_{IR}/L_{\odot})>11.5}$), our results are
higher than both of these estimates.

Although the simulations described above are based on different assumptions,
the theoretical work does not contain a UV selected sample bias. This suggests
that, according to simulations, our IRLF estimate is not missing a significant
number of extremely luminous, UV-undetected galaxies at $\mathrm{z\sim 7}$
(for potential caveats, see Section 5.3).

We continue by contrasting with semi-empirical models from Zavala et al.
(2021) at $\mathrm{z\sim 7}$. Their IRLF changes very little at
$\mathrm{12<log(L_{IR}/L_{\odot})<12.5}$, whereas our IRLF sharply declines.
Our study shows an IRLF an order of magnitude higher for LIRGs and a
negligible number of galaxies with $\mathrm{log(L_{IR}/L_{\odot})>12.3}$.
Thus, we find a different distribution also for the bright luminosity end.
These differences in IRLF could be explained by the different methodology, due
to the lack of observational data at $\mathrm{z\sim 7}$, that leads to an
extrapolation of their IRLF at higher redshifts. To do that, it is necessary
to assume two different slopes for the LIRGs and the ULIRGs that might lead to
different outcomes between our study and Zavala et al. (2021).

Finally, we compare our results with IRLFs derived from observations. In
particular, we contrast with the ALPINE IRLF, since it is an analogous survey
to REBELS, but at lower redshift (see section 2 for details). Using the ALPINE
data, Gruppioni et al. (2020) provide the IRLF at $\mathrm{z\sim 5}$ for
serendipitous galaxies. Their IRLF agrees with ours for the lower luminosity
bin, but the overall normalisation is significantly higher. The reason of the
difference is the IRLF rely on several factors. Firstly, the redshift
difference ($\mathrm{\Delta z\sim 2}$) is an obvious reason for the density to
be lower. Furthermore, the REBELS sample was UV-selected, implying a selection
effect that is nonexistent in a blind survey (see section 5.3 for caveats).
Another cause for the disparity with Gruppioni et al. (2020) might the
difference it redshift calculation. Their redshifts were calculated with
multi-band photometry and with only three galaxies at $\mathrm{z\sim 5}$.
Finally, the differing dust temperature assumptions and the SED fitting may
lead to different infrared luminosities, but further analysis is required to
ensure that the differences are significant.

In order to continue the observational comparison, we contrast the IRLF
calculated with the maximum redshift observed to yet in Wang et al. (2019a).
This analysis presents an IRLF with bright infrared galaxies selected with
Herschel Space Observatory Pilbratt et al. (2010) at $\mathrm{z=5.5}$. At same
redshift, their results have a 2 dex greater luminosity function than ours at
the bright end, but a smaller overall luminosity function than the one stated
in Gruppioni et al. (2020). Again, the expected difference is caused by the
disparity in redshift, as does the bias to select massive galaxies with
Herschel.

###  5.2 What IRLF is needed to reproduce extreme SFRD?

This section discusses how changes in the IRLF impact the SFRD. Since there is
lack of consensus about obscured SFRDs at $\mathrm{z>5}$, we evaluate the key
variables that influence the SFRD computation: the IRLF faint end slope, the
$\mathrm{L_{IR}}$ integration limits, and the conversion factor between
$\mathrm{L_{IR}}$ and SFRD. To do that, we compute the SFRD derived for
extreme $\mathrm{\alpha}$ and integration limits to determine whether the most
extreme SFRD described in the literature could be reproduced. We also discuss
the likely causes of these variances.

First, we investigate changes in the IRLF slope. Lower redshift studies
frequently find a slope of $\mathrm{\alpha=-1.3}$, including more galaxies
with lower infrared luminosities (Hammer et al., 2012), but some high redshift
studies report shallower faint-end slopes of $\mathrm{\alpha=-0.4}$ (Koprowski
et al., 2017; Zavala et al., 2021). In Figure 5, we compute the IRLF for these
two extreme cases by using $\mathrm{\alpha=-2}$ and $\mathrm{\alpha=-0.4}$,
respectively. Additionally, we used a wider luminosity range for the
integration than in previous sections of this work, allowing for
$\mathrm{8<log(L_{IR}/L_{\odot})<13}$ as in (Gruppioni et al., 2020).
Nevertheless, we cannot recreate values close to their SFRD, even in the most
extreme scenario ($\mathrm{\alpha=-2}$), yielding a $\mathrm{SFRD\sim 6\cdot
10^{-3}M_{\odot}/Mpc^{3}/yr}$.

This SFRD is, however, consistent with the findings of Talia et al. (2021)
($\mathrm{SFRD\sim 5\cdot 10^{-3}M_{\odot}/yr/Mpc^{3}}$ at $\mathrm{z\sim
5}$). It should be noted that the analysis of Talia et al. (2021) was
conducted using radio galaxies with median $\mathrm{L_{IR}=2.3\pm 0.5\times
10^{12}L_{\odot}}$, and is thus based on a different set of assumptions than
our IR-based estimates.

Despite the fact that it is common to compute the obscured SFRD using the
IRLF, some studies directly calculate it by using the individual SFRs. For
instance, the MORA survey performed blind 2mm ALMA observations (Casey et al.,
2021), and identified a number of $\mathrm{z\sim 4-6}$ DSFGs. They find
$\mathrm{SFRD\sim 10^{-3}\ M_{\odot}/yr/Mpc^{3}}$ at $\mathrm{z\sim 6}$, which
is far lower than the previously mentioned studies such as Talia et al. (2021)
or Gruppioni et al. (2020). The key distinction is that their photometric
redshift estimates are based on submillimetre data, which can be degenerate
with dust temperature. Generally, however, the findings of Casey et al. (2021)
are in good agreement with ours, and their obscured SFRD is compatible with a
$\mathrm{z\sim 6}$ extension of our SFRD at $\mathrm{z\sim 7}$. This agreement
also extends to the 1.3 mm ALMA serendipitous sources at $\mathrm{z<4.5}$ from
Dunlop et al. (2017). Both Dunlop et al. (2017) and Casey et al. (2021)
present a decrease of obscured SFRD at $\mathrm{z>3}$ which likely continues
beyond $\mathrm{z>6}$, as suggested by our data.

Even if several obscured SFRD present large values at $\mathrm{z\sim 5}$ (i.e.
Wang et al. (2019a); Gruppioni et al. (2020); Khusanova et al. (2021)), we
also notice that the highest redshift bin in both Wang et al. (2019a) and
Gruppioni et al. (2020) have larger uncertainty than the rest to the low
number of sources (as shown the hatched areas in Figure 4). Given these larger
uncertainties, a declining SFRD cannot be excluded from these analyses. Hence,
although not in agreement, our results are not in contradiction with the
studies that show large SFRDs and the highest redshift surveys. Studies
including larger samples at $\mathrm{4<z<7}$ would be needed to corroborate
this hypothesis.

Figure 5:  The SFRD depends on the IRLF shape and the luminosity range used in
the integration. The faint end slope $\alpha$ assumed in the low luminosity
end is key for the resulting SFRD. This plot shows the best fit IRLF for two
extreme slopes: $\mathrm{\alpha=-2}$ (red line) and $\mathrm{\alpha=-0.4}$
(orange line). The difference between slopes increases in IRLF being
$\mathrm{\sim 4}$ orders of magnitude higher at
$\mathrm{L_{IR}=10^{9}L_{\odot}}$. The inner plot shows the SFRD for these two
extreme cases which shows an order of magnitude difference depending on the
slope assumed with the same integration luminosity
($\mathrm{10^{8}<L_{IR}/L_{\odot}<10^{13}}$). The dark red dots show the total
REBELS sample for the three luminosity bins. The dark red line shows the
Schechter fit with $\mathrm{\alpha=-1.3}$ (dark red line) as presented
previously in Section 3.

###  5.3 Possible Caveats

In this section, we assess the importance of our data being based on a UV-
bright target selection. This directly implies that our study cannot account
for extremely dust-obscured sources, such as SMGs, that are faint in the UV.
However, given that there are several verified SMGs at $\mathrm{z>4}$, we know
that such galaxies are 100$\times$ less common than UV-based Lyman Break
Galaxies, given the SMG sky surface density of 0.01 arcmin-2 (e.g., Riechers
et al., 2013, 2017; Marrone et al., 2018). Furthermore, extremely dusty high
redshift galaxies have only been discovered up to a maximum $\mathrm{z=6.34}$
(Riechers et al., 2013). All of these findings are based on large surveys
conducted with the South Pole Telescope (SPT), SCUBA-2, or Herschel Space
Observatory.

The serendipitous detection of two dust-obscured galaxies in the REBELS
dataset with similar masses and SFRs as the main sample clearly shows that the
primary target sample of REBELS is not complete (Fudamoto et al., 2021).
While, the contribution of this class of galaxies to the SFRD is still very
uncertain, Fudamoto et al. (2021) estimate a value of $\mathrm{1.2\times
10^{-3}\,M_{\odot}/yr/Mpc^{3}}$, i.e. comparable to our estimate from the
IRLF. This would suggest that UV-undetected galaxies could contribute a
similar, but additional amount of obscured SFR as UV-bright galaxies.

Similar conclusions have been reached from recent JWST observations. The first
deep NIRCam observations revealed the existence of UV-undetected, dusty
galaxies at $\mathrm{z>6}$. Barrufet et al. (2022), in particular, present the
SFRD for high-z dusty galaxies, finding a
$\mathrm{log(SFRD/M_{\odot}/yr/Mpc^{3})\sim-3}$ at $\mathrm{z\sim 7}$ for
highly attenuated galaxies. We thus conclude that the galaxies we are missing
in UV selections might contribute the same order of magnitude as the REBELS
sample itself.

To compute a more complete IRLF it would be necessary to perform a deep, but
blind survey to probe galaxies at $z\sim 7$ at several wavelengths. For the
present, a good first step is to obtain results based on the UV-selected
REBELS galaxies. These results represent a firm lower limit on the total
obscured SFRD at $\mathrm{z\sim 7}$.

##  6 Summary and Conclusions

In this work, we have exploited the data from the REBELS survey, which
consists of ALMA spectroscopic data of UV-bright galaxies in the EoR. Our
sample consists of 42 galaxies at $\mathrm{6.4<z<7.7}$. 16 have revealed
significant dust continuum emission at rest-frame $\sim 158\micron$, and all
but one of these are spectroscopically confirmed through their [CII] emission
lines. This sample was used to:

  * •

We have calculated the Infrared Luminosity Function (IRLF) at $\mathrm{z\sim
7}$ for the first time using a spectroscopically confirmed sample. We find a
$\mathrm{log(\phi_{IR})\sim-4.2\pm 0.2\ dex^{-1}Mpc^{-3}}$ in our faintest
luminosity bin of $\mathrm{log(L_{IR}/L_{\odot})\sim 11.5}$. At higher
luminosities, the IRLF decreases considerably.

  * •

We have fit a Schechter (1976) function with a fix slope of
$\mathrm{\alpha=-1.3}$ for the low luminosity end finding the best fitting
values $\mathrm{log(\phi_{IR})\sim-4.38\ dex^{-1}Mpc^{-3}}$ and
$\mathrm{log(L_{IR}/L_{\odot})=11.6}$. Our results indicate that extremely
luminous galaxies with $\mathrm{log(L_{IR}/L_{\odot})>12.3}$ are extremely
rare at $z\sim 7$, with number densities
$\mathrm{log(\phi_{IR})<-6.5dex^{-1}Mpc^{-3}}$.

  * •

We have derived the obscured Star Formation Rate Density through the IRLF.
From the observations we calculate a lower limit of
$\mathrm{log(SFRD/M_{\odot}/yr/Mpc^{3})=-2.93\pm 0.20}$ at $\mathrm{z\sim 7}$
which represents $\mathrm{\sim 13\%}$ of the total SFRD. When integrating over
the luminosity range $\mathrm{10.5<log(L_{IR}/L_{\odot})<13}$ we infer a
larger value of
$\mathrm{log(SFRD/M_{\odot}/yr/Mpc^{3})=-2.66^{+0.17}_{-0.14}}$.

  * •

Our IRLF is broadly consistent with some simulations at $\mathrm{z\sim 7}$.
The inferred SFRD is a robust lower limit that shows a significant
contribution of obscured star formation at $\mathrm{z\sim 7}$.

We conclude that our results imply a significant amount of obscured SFR at
$\mathrm{z\sim 7}$ of at least
$\mathrm{log(SFRD/M_{\odot}/yr/Mpc^{3})\sim-3}$. Comparing with ALMA blind
surveys, our results suggest a steep evolution of the obscured SFRD over
redshift that continues to $\mathrm{z\sim 7}$, at least.

## Acknowledgements

We acknowledge the constructive feedback of the referee (MB) for his
constructive feedback that helped in the improvement of this paper. We
acknowledge support from: the Swiss National Science Foundation through the
SNSF Professorship grant 190079 (LB and PAO). The Cosmic Dawn Center (DAWN) is
funded by the Danish National Research Foundation under grant No. 140. PD
acknowledges support from the European Research Council’s starting grant ERC
StG-717001 (“DELPHI"), from the NWO grant 016.VIDI.189.162 (“ODIN") and the
European Commission’s and University of Groningen’s CO-FUND Rosalind Franklin
program. AF and AP acknowledge support from the ERC Advanced Grant
INTERSTELLAR H2020/740120. Generous support from the Carl Friedrich von
Siemens-Forschungspreis der Alexander von Humboldt-Stiftung Research Award is
kindly acknowledged. YF acknowledges support from NAOJ ALMA Scientific
Research Grant number 2020-16B. VG gratefully acknowledges support by the ANID
BASAL projects ACE210002 and FB210003.

## References

  * Adams et al. (2022) Adams N. J., et al., 2022, arXiv e-prints, p. arXiv:2207.11217
  * Algera et al. (2022) Algera H., et al., 2022, MNRAS, pp 6142–6157
  * Atek et al. (2022) Atek H., Shuntov M., Furtak L. J., Richard J., Kneib J.-P., Mahler Adi Zitrin G., McCracken Clotilde Laigle Stéphane Charlot H. J., 2022, arXiv e-prints, p. arXiv:2207.12338
  * Barrufet et al. (2020) Barrufet L., et al., 2020, A&A, 641, A129
  * Barrufet et al. (2022) Barrufet L., et al., 2022, arXiv e-prints, p. arXiv:2207.14733
  * Béthermin et al. (2020) Béthermin M., et al., 2020, A&A, 643, A2
  * Bouwens (2016) Bouwens R., 2016, in Mesinger A., ed., Astrophysics and Space Science Library Vol. 423, Astrophysics and Space Science Library. p. 111 (arXiv:1511.01133), doi:10.1007/978-3-319-21957-8_4
  * Bouwens et al. (2007) Bouwens R. J., Illingworth G. D., Franx M., Ford H., 2007, ApJ, 670, 928
  * Bouwens et al. (2015) Bouwens R. J., et al., 2015, ApJ, 803, 34
  * Bouwens et al. (2021a) Bouwens R. J., et al., 2021a, arXiv e-prints, p. arXiv:2106.13719
  * Bouwens et al. (2021b) Bouwens R. J., et al., 2021b, AJ, 162, 47
  * Bouwens et al. (2022) Bouwens R. J., Illingworth G. D., Ellis R. S., Oesch P. A., Stefanon M., 2022, arXiv e-prints, p. arXiv:2205.11526
  * Bowler et al. (2018) Bowler R. A. A., Bourne N., Dunlop J. S., McLure R. J., McLeod D. J., 2018, MNRAS, 481, 1631
  * Bowler et al. (2020) Bowler R. A. A., Jarvis M. J., Dunlop J. S., McLure R. J., McLeod D. J., Adams N. J., Milvang-Jensen B., McCracken H. J., 2020, MNRAS, 493, 2059
  * Capak et al. (2015) Capak P. L., et al., 2015, Nature, 522, 455
  * Casey et al. (2019) Casey C., et al., 2019, BAAS, 51, 212
  * Casey et al. (2021) Casey C. M., et al., 2021, ApJ, 923, 215
  * Dayal & Ferrara (2018) Dayal P., Ferrara A., 2018, Phys. Rep., 780, 1
  * Dayal et al. (2022) Dayal P., et al., 2022, MNRAS, 512, 989
  * Donnan et al. (2022) Donnan C. T., et al., 2022, arXiv e-prints, p. arXiv:2207.12356
  * Dunlop (2013) Dunlop J. S., 2013, in Wiklind T., Mobasher B., Bromm V., eds, Astrophysics and Space Science Library Vol. 396, The First Galaxies. p. 223 (arXiv:1205.1543), doi:10.1007/978-3-642-32362-1_5
  * Dunlop et al. (2017) Dunlop J. S., et al., 2017, MNRAS, 466, 861
  * Endsley et al. (2022) Endsley R., et al., 2022, arXiv e-prints, p. arXiv:2202.01219
  * Faisst et al. (2020) Faisst A. L., et al., 2020, ApJS, 247, 61
  * Ferrara et al. (2022) Ferrara A., et al., 2022, arXiv e-prints, p. arXiv:2202.07666
  * Finkelstein et al. (2022) Finkelstein S. L., et al., 2022, arXiv e-prints, p. arXiv:2207.12474
  * Fudamoto et al. (2021) Fudamoto Y., et al., 2021, Nature, 597, 489
  * Fujimoto et al. (2023) Fujimoto S., Kohno K., Ouchi M., 2023, ApJ
  * Gruppioni et al. (2013) Gruppioni C., et al., 2013, MNRAS, 432, 23
  * Gruppioni et al. (2020) Gruppioni C., et al., 2020, A&A, 643, A8
  * Hammer et al. (2012) Hammer D. M., Hornschemeier A. E., Salim S., Smith R., Jenkins L., Mobasher B., Miller N., Ferguson H., 2012, ApJ, 745, 177
  * Harikane et al. (2021) Harikane Y., et al., 2021, arXiv e-prints, p. arXiv:2112.09141
  * Harikane et al. (2022) Harikane Y., et al., 2022, arXiv e-prints, p. arXiv:2208.01612
  * Heintz et al. (2022) Heintz K. E., et al., 2022, ApJ, 934, L27
  * Hodge & da Cunha (2020) Hodge J. A., da Cunha E., 2020, Royal Society Open Science, 7, 200556
  * Inami et al. (2022) Inami H., et al., 2022, arXiv e-prints, p. arXiv:2203.15136
  * Johnston (2011) Johnston R., 2011, A&ARv, 19, 41
  * Khusanova et al. (2021) Khusanova Y., et al., 2021, A&A, 649, A152
  * Koprowski et al. (2017) Koprowski M. P., Dunlop J. S., Michałowski M. J., Coppin K. E. K., Geach J. E., McLure R. J., Scott D., van der Werf P. P., 2017, MNRAS, 471, 4155
  * Koprowski et al. (2020) Koprowski M. P., et al., 2020, MNRAS, 492, 4927
  * Labbe et al. (2022) Labbe I., et al., 2022, arXiv e-prints, p. arXiv:2207.12446
  * Lagos et al. (2020) Lagos C. d. P., da Cunha E., Robotham A. S. G., Obreschkow D., Valentino F., Fujimoto S., Magdis G. E., Tobar R., 2020, MNRAS, 499, 1948
  * Laporte et al. (2019) Laporte N., et al., 2019, MNRAS, 487, L81
  * Laporte et al. (2021) Laporte N., Meyer R. A., Ellis R. S., Robertson B. E., Chisholm J., Roberts-Borsani G. W., 2021, MNRAS, 505, 3336
  * Le Fèvre et al. (2020) Le Fèvre O., et al., 2020, A&A, 643, A1
  * Lim et al. (2020) Lim C.-F., et al., 2020, ApJ, 889, 80
  * Madau & Dickinson (2014) Madau P., Dickinson M., 2014, ARA&A, 52, 415
  * Marrone et al. (2018) Marrone D. P., et al., 2018, Nature, 553, 51
  * Naidu et al. (2022a) Naidu R. P., et al., 2022a, arXiv e-prints, p. arXiv:2207.09434
  * Naidu et al. (2022b) Naidu R. P., et al., 2022b, arXiv e-prints, p. arXiv:2208.02794
  * Nelson et al. (2022) Nelson E. J., et al., 2022, arXiv e-prints, p. arXiv:2208.01630
  * Oesch et al. (2016) Oesch P. A., et al., 2016, ApJ, 819, 129
  * Oesch et al. (2018a) Oesch P. A., et al., 2018a, ApJS, 237, 12
  * Oesch et al. (2018b) Oesch P. A., Bouwens R. J., Illingworth G. D., Labbé I., Stefanon M., 2018b, ApJ, 855, 105
  * Pilbratt et al. (2010) Pilbratt G. L., et al., 2010, A&A, 518, L1
  * Popping et al. (2020) Popping G., et al., 2020, ApJ, 891, 135
  * Pozzi et al. (2021) Pozzi F., et al., 2021, A&A, 653, A84
  * Riechers et al. (2013) Riechers D. A., et al., 2013, Nature, 496, 329
  * Riechers et al. (2017) Riechers D. A., et al., 2017, ApJ, 850, 1
  * Rodighiero et al. (2022) Rodighiero G., Bisigello L., Iani E., Marasco A., Grazian A., Sinigaglia F., Cassata P., Gruppioni C., 2022, arXiv e-prints, p. arXiv:2208.02825
  * Schaerer et al. (2022) Schaerer D., Marques-Chaves R., Barrufet L., Oesch P., Izotov Y. I., Naidu R., Guseva N. G., Brammer G., 2022, arXiv e-prints, p. arXiv:2207.10034
  * Schechter (1976) Schechter P., 1976, ApJ, 203, 297
  * Schmidt (1968) Schmidt M., 1968, ApJ, 151, 393
  * Schouws et al. (2021) Schouws S., et al., 2021, arXiv e-prints, p. arXiv:2105.12133
  * Schouws et al. (2022) Schouws S., et al., 2022, arXiv e-prints, p. arXiv:2202.04080
  * Shen et al. (2021) Shen X., Vogelsberger M., Nelson D., Tacchella S., Hernquist L., Springel V., Marinacci F., Torrey P., 2021, arXiv e-prints, p. arXiv:2104.12788
  * Smit et al. (2018) Smit R., et al., 2018, Nature, 553, 178
  * Sommovigo et al. (2021) Sommovigo L., Ferrara A., Carniani S., Zanella A., Pallottini A., Gallerani S., Vallini L., 2021, MNRAS, 503, 4878
  * Sommovigo et al. (2022) Sommovigo L., et al., 2022, MNRAS,
  * Stark (2016) Stark D. P., 2016, ARA&A, 54, 761
  * Talia et al. (2021) Talia M., Cimatti A., Giulietti M., Zamorani G., Bethermin M., Faisst A., Le Fèvre O., Smolçić V., 2021, ApJ, 909, 23
  * Topping et al. (2022) Topping M. W., et al., 2022, MNRAS, 516, 975
  * Viero et al. (2022) Viero M. P., Sun G., Chung D. T., Moncelsi L., Condon S. S., 2022, arXiv e-prints, p. arXiv:2203.14312
  * Wang et al. (2019a) Wang T., et al., 2019a, Nature, 572, 211
  * Wang et al. (2019b) Wang L., Pearson W. J., Cowley W., Trayford J. W., Béthermin M., Gruppioni C., Hurley P., Michałowski M. J., 2019b, A&A, 624, A98
  * Watson et al. (2015) Watson D., Christensen L., Knudsen K. K., Richard J., Gallazzi A., Michałowski M. J., 2015, Nature, 519, 327
  * Yan et al. (2020) Yan L., et al., 2020, ApJ, 905, 147
  * Zavala et al. (2021) Zavala J. A., et al., 2021, ApJ, 909, 165