# Dust in High-Redshift Galaxies: Reconciling UV Attenuation and IR Emission

Jianhan Zhao,11footnotetext: Corresponding author. Steven R. Furlanetto

###### Abstract

Dust is a key component of galaxies, but its properties during the earliest
eras of structure formation remain elusive. Here we present a simple semi-
analytic model of the dust distribution in galaxies at $z\gtrsim 5$. We
calibrate the free parameters of this model to estimates of the UV attenuation
(using the IRX-$\beta$ relation between infrared emission and the UV spectral
slope) and to ALMA measurements of dust emission. We find that the observed
dust emission requires that most of the dust expected in these galaxies is
retained (assuming a similar yield to lower-redshift sources), but if the dust
is spherically distributed, the modest attenuation requires that it be
significantly more extended than the stars. Interestingly, the retention
fraction is larger for less massive galaxies in our model. However, the
required radius is a significant fraction of the host’s virial radius and is
larger than the estimated extent of dust emission from stacked high-$z$
galaxies. These can be reconciled if the dust is distributed anisotropically,
with typical covering fractions of $\sim 0.2$–0.7 in bright galaxies and
$\lesssim 0.1$ in fainter ones.

##  1 Introduction

Dust is crucial for understanding galaxy evolution on a variety of levels. Not
only is it a key component in chemical evolution, but it also absorbs a
substantial fraction of UV starlight and re-emits it in the infrared. However,
its role in galaxies during the Cosmic Dawn is mostly unknown, as only
recently have the tools emerged to measure its properties.

JWST has, of course, made crucial contributions by pushing the observational
frontier to higher redshifts and fainter sources (e.g., [1, 2, 3, 4, 5, 6,
7]). Interestingly, these studies have found an excess of UV-luminous galaxy
candidates at $z\sim 10$ compared to HST-calibrated models [8, 9, 10]. The
most popular explanation is burstiness in the star formation histories [8].
However, dust may play a further role in this phenomenon [10, 11] if the
highest-redshift sources are mostly dust-free. [9] also argued that the burst
solution also depends on dust: more efficient dust production is needed to
match the JWST-observed UV colours while maintaining a star formation
efficiency that fits the observed luminosity function. Others have noted that
the dust may have a complex morphological relation to the star-forming regions
[12, 11]. A thorough study of high-redshift dust properties is called for to
answer these questions.

Beyond the luminosity function, an understanding of dust is also crucial to
interpreting galaxy properties during the Cosmic Dawn. Dust grains attenuate
starlight, especially in the UV [13, 14, 15]. Therefore, it is critical for
measuring the intrinsic SFR as well as other properties of the stellar
population. However, most high-$z$ galaxies are only identified
photometrically, and SED modeling is required to estimate the effects of dust.
A common simple method is to use the UV spectral slope $\beta$ as a proxy for
the dust distribution, which appears to be a good assumption for lower-
redshift galaxies, although the details may evolve with redshift [16, 17, 18,
19, 20, 21, 22]. The fundamental motivation behind these empirical relations
is that more dust (and hence a larger infrared luminosity) will redden the
galaxy more (making the spectral slope flatter).

The advent of ALMA has allowed new studies directly targeting dust emission in
the far-infrared from high-$z$ galaxies [21, 23, 24, 25, 26, 12, 27, 28, 29].
In particular, large surveys like ALPINE [30, 31], CRISTAL [32], SERENADE
[33], and REBELS [34] have begun to place statistical constraints on the dust
content of early galaxies. With the shared objective of studying dust and
other ISM properties in high-$z$ galaxies, these surveys scan one or both of
the brightest fine structure lines ([O III] at 88$\mu$m and [C II] at
158$\mu$m) from UV or spectroscopically selected samples, which also allows a
measurement of the nearby dust continuum. These surveys cover a broad range of
redshifts (ALPINE and CRISTAL at $4.4\lesssim z\lesssim 5.9$, SERENADE at
$z\sim 6$, and REBELS at $z\sim 7$).

However, constructing a first-principles model of dust is a difficult task. It
requires a prescription for dust formation in supernova ejecta [35, 36], AGB
winds [37], ISM grain growth [38, 25], and other processes. The dust must then
be distributed through the interstellar and circumgalactic media through winds
and turbulent transport, leaving a complex geometrical distribution [39]. It
then attenuates starlight in a manner that depends on the composition and size
distribution of the dust grains [13, 40]. Recent work has brought a variety of
methods to bear on this problem. Numerical simulations [41, 42, 43, 44] can
track the evolution of dust content in detail but are still dependent on
assumptions about feedback and dust creation and destruction. Analytic models
offer more flexibility but are limited in the number of physical mechanisms
they can include [45, 46, 47]. Many of these models use parameterized semi-
analytic approaches to account for these challenges. [24, 25, 26, 48, 29]

In this paper, we use the new JWST and ALMA observations to understand dust
through a simplified model of galaxy formation. We use the “minimalist” model
of high-$z$ galaxy evolution from [49] as our baseline and incorporate a new
treatment of dust production and morphology. We then use both UV and far-IR
measurements to constrain the properties of the model.

The paper is organized as follows. In section 2, we introduce our galaxy
formation and dust models, building on [49]; in section 3, we compare to UV
dust measurements; in section 4, we appeal to the ALMA far-IR flux observation
to further constrain the model; in section 5, we explore the implications of
our model; and in section 6 we conclude.

The numerical calculations in this work assume a flat $\Lambda$CDM cosmology
with $\Omega_{m}=0.3111$, $\Omega_{\Lambda}=0.6889$, $\Omega_{b}=0.0489$,
$h=0.6766$, consistent with recent results of [50]. The magnitudes in this
paper are expressed in the AB system [51].

##  2 Dust in High-redshift Galaxies

This section introduces our dust prescription, which extends the star
formation model of [49]. Since our objective is to constrain the dust mass and
radius using multi-band observations, we extend their minimalist approach to
include both dust evolution and its observable consequences. Rather than
attempt to model dust from first principles, we will introduce scaling
parameters for the dust mass and radius, to be calibrated by observations in
later sections. The resultant model depends on the host dark matter halo mass
$M_{\rm halo}$ and redshift $z$ of a galaxy, generating its (1) dust mass, (2)
dust radius, (3) UV attenuation, and (4) dust temperature and far-IR continuum
emission flux.

###  2.1 The Minimalist Galaxy Formation Model

Here, we introduce our basis model: a minimalist semi-analytic star formation
model by [49]. Calibrated to the HST-observed UVLF, this model predicts
various properties of a galaxy by balancing its star formation and feedback,
while assuming galaxies grow through matter accretion onto their host halos.
We note that [49] only provides a reasonable fit to the UVLF up to $z\approx
10$ [9], so we must be cautious applying it to earlier eras.

The model of [49] uses a simple approximation to dark matter halo growth. We
first define the halo mass function, by assuming that the comoving number
density $n_{\rm h}$ of halos in the infinitesimal mass range $(M,M+dM)$ at a
given redshift $z$ is

$n_{\rm h}(M,z)=f(\sigma)\frac{\bar{\rho}}{M}\frac{d\ln{(1/\sigma)}}{dM},$
(2.1)

where $\bar{\rho}$ is the comoving matter density of halos, $\sigma$ is the
RMS fluctuation of the linear density field smoothed on a scale $M$, and
$f(\sigma)$ is a dimensionless parametrization usually fit to a simulation.
Here, we take $f(\sigma)$ from [52]’s high-$z$ simulation results, where

$f(\sigma)=0.15\left[1+\left(\frac{\sigma}{2.54}\right)^{-1.36}\right]e^{-1.14/\sigma^{2}}.$
(2.2)

To find the desired halo accretion rate, we then apply a method analogous to
abundance matching [53], where we deduce the halo mass accretion rate by
assuming that halos evolve at constant number density across redshifts. That
is, we assume that the masses $(M_{1},M_{2})$ at $(z_{1},z_{2})$ are related
via

$\int_{M_{1}}^{\infty}dM\,n_{\mathrm{h}}(M|z_{1})=\int_{M_{2}}^{\infty}dM\,n_{\mathrm{h}}(M|z_{2}).$
(2.3)

We obtain the accretion rate $\dot{M}_{\rm h}$ by taking the limit
$z_{1}\rightarrow z_{2}$. We assume that the baryonic accretion rate, which is
the source of star and dust formation, is simply proportional to the total
accretion rate in most halos, except at high masses. Simulations have shown
that high-mass halos are quenched by their virial shocks, so we apply a shock
suppression factor $f_{\mathrm{sh}}$ to the baryonic accretion rate given by
the fit of [54]. Thus the total accretion rate is $f_{\rm sh}\dot{M}_{b}\equiv
f_{\rm sh}(\Omega_{b}/\Omega_{m})\dot{M}_{\rm h}$.

In the minimalist model, we assume the baryonic accretion rate is in
equilibrium with star formation and feedback, such that

$\dot{M}_{b}=\dot{M}_{\star}+\dot{M}_{w},$ (2.4)

where $\dot{M}_{\star}$ is the total star formation rate, and $\dot{M}_{w}$ is
the rate of baryonic expulsion due to feedback. As this feedback is generated
by the stars, we assume that $\dot{M}_{w}=\eta\dot{M}_{\star}$, where the
mass-loading parameter $\eta$ determines the strength of such feedback and is
in general a function of halo mass and redshift. We also define the
instantaneous star formation efficiency as the fraction of the accreted
baryons that turn into stars, such that
$f_{\star}\equiv\dot{M}_{\star}/\dot{M}_{b}$. Additionally, we impose an upper
limit of star formation efficiency denoted $f_{\star,\mathrm{max}}$ to
maintain the continuity of the resultant function. Our finalized star
formation efficiency prescription therefore becomes

$f_{\star}=\frac{f_{\mathrm{shock}}}{f^{-1}_{\star,\mathrm{max}}+\eta(M_{\mathrm{h}},z)}.$
(2.5)

We now require a prescription for the feedback parameter $\eta$. Our
parameterization is based on two limiting scenarios. In the first, we assume
that all of the kinetic energy from supernovae is available to drive the
outflow, so we balance the supernova energy with the binding energy of the
accreting gas,

$\frac{1}{2}\dot{M}_{w}v^{2}_{\mathrm{esc}}=\dot{M}_{\star}\epsilon_{k}\omega_{\mathrm{SN}},$
(2.6)

where $v_{\mathrm{esc}}$ is the halo’s escape velocity and
$\omega_{\mathrm{SN}}=10^{49}\mathrm{erg/M_{\odot}}$ is the fiducial supernova
kinetic energy per unit mass of star formation, where the value here assumes a
typical stellar initial mass function (IMF) and metallicity. This method
assumes that the gas always receives a fixed portion of supernova blast
energy. Alternatively, if most of the energy is lost to radiative cooling, we
could obtain a similar relation by balancing the momentum of the supernova
ejecta with the escape velocity. In both cases, the feedback parameter can be
fit by the general form

$\eta(M_{\mathrm{h}},z)=A\left(\frac{10^{11.5}M_{\odot}}{M_{\mathrm{h}}}\right)^{\xi}\left(\frac{9}{1+z}\right)^{\sigma},$
(2.7)

where $\xi$ and $\sigma$ are the power law indices and $A$ is the
normalization factor. Setting the remaining parameters (including the overall
normalization) with an approximate fit to the $z\sim 6$ UV LF yields
$(A,\,\xi,\,\sigma,\,f_{\rm\star,max})\approx(1,\,2/3,\,1,\,0.1)$ for the
energy-regulated case and $(5,\,1/3,\,1/2,\,0.2)$ for the momentum-regulated
case [55]. We will use the latter as our fiducial model because it is more
consistent with numerical simulations of feedback regulation (e.g., the FIRE
simulations of [56]), but this choice does not significantly alter our
conclusions.

To compare with existing observations, we convert the computed star formation
rate of a galaxy to the corresponding UV luminosity it generates. Throughout
this paper, we adopt the SFR-$L_{\mathrm{UV}}$ conversion from [57], which can
be written

$\dot{M}_{\mathrm{\star}}=\mathcal{K}_{\mathrm{UV}}\times L_{\mathrm{UV}},$
(2.8)

where the proportionality constant is $\mathcal{K}_{\mathrm{UV}}=1.15\times
10^{-28}\mathrm{M_{\odot}yr^{-1}/(erg\,s^{-1}Hz^{-1})}$ and $L_{\mathrm{UV}}$
is the continuum UV luminosity measured at 1500 Å. This conversion assumes a
Salpeter IMF from 0.1 to 100 $M_{\odot}$ and an extended period of continuous
star formation for $\gtrsim 300\ {\rm Myr}$.

We have defined the star formation efficiency as an instantaneous value, but
it will be useful to have the overall (time-averaged) value for a galaxy. We
define $\tilde{f}_{\star}\equiv M_{\star}/M_{b}$, where
$M_{b}=(\Omega_{b}/\Omega_{m})M_{\rm h}$. [49] showed that
$\tilde{f}_{\star}(M_{\mathrm{halo}},z)\approx
f_{\star}(M_{\mathrm{halo}},z)/(1+\xi)$, which we use as an estimate here; the
slight error can be absorbed into our scaling constants.

###  2.2 Dust Production in High-Redshift Galaxies

We next assume that the total dust mass is proportional to the total stellar
mass with a yield $y_{d}$, such that the total dust mass is $y_{d}M_{\star}$.
We will use a fiducial $y_{d}=0.01$ value comparable to a typical disk galaxy
from low redshifts. However, this parameter is uncertain at high redshifts,
because the dust production mechanisms may differ, or a portion of the dust
can be ejected by feedback [11] or destroyed by supernova shock waves [58]. We
therefore define a scaling parameter $f_{M}$ via

$M_{\mathrm{dust}}\equiv f_{M}\times(0.01M_{\star}).$ (2.9)

We naively expect $f_{M}\lesssim 1$, because dust formation may be less
efficient at early times and because dust may be ejected or destroyed. For
example, the disk models of [59] found that only $\sim 0.1$ of the total
baryonic mass associated with a halo was retained in the galaxy’s interstellar
medium, with most of the remainder ejected. We absorb those uncertainties into
the mass scaling parameter $f_{M}$, which we will fix by comparison to
observations in later sections.

We note that we do not explicitly consider the effect of asymptotic giant
branch (AGB) stars’ dust production throughout this paper. Though AGB stars
are the main sources of dust in local systems [60] such as the Milky Way, they
are less likely to constitute a significant portion in higher redshifts
because of their long evolutionary timescales [38]. Their cumulative effects
could, however, also be absorbed into $f_{M}$, although that would not fully
capture the long timescale of dust production.

The dotted curves in Figure 1 show the resulting dust masses, assuming
$f_{M}=1$. We show them as a function of each galaxy’s SFR, because
observations are often presented that way. In our minimalist model, the SFR is
a monotonic function of halo (and stellar) mass, so we provide the translation
along the top axis. The solid curves preview the results of our model
calibration to UV and IR observations – we will find that $f_{M}$ may be a
strong function of halo mass. The overall efficiency of dust formation is thus
$\approx 0.01f_{M}\tilde{f}_{\star}$, which for halos in the mass range
$10^{8}$–$10^{12}\ M_{\odot}$ is typically $\sim 3\times 10^{-4}f_{M}$. The
expected dust mass is smaller in smaller halos, because their star formation
efficiencies are small, and also in very large halos, because our virial shock
treatment suppresses accretion and hence star formation in these systems as
well.

Figure 1: Expected dust masses as a function of SFR and $M_{\star}$ at $z=7$
(bottom and top axes; these are monotonically related in our model). The
dotted curves show the dust mass assuming $f_{M}=1$, while the solid curves
are the calibrated dust masses with $f_{M}$ formulated in section 4.2. The
dash-dotted curve shows [25]’s dust masses at $z=7$.

Before developing our model further, we briefly compare our predicted dust
masses with [25], who also use a semi-analytic dust treatment, albeit with
more details in modelling specific dust evolution processes, such as dust
grain growth, destruction, and astration. Their model uses a dark matter halo
merger tree and traces the gas content associated with such halos to deduce
the dust mass. In Figure 1, we compare our work’s fiducial dust mass (with
$f_{M}=1$; dotted curves) to [25]’s results presented in their equation (11)
using our model derived $M_{\star}$ (shown in the top axis on the left panel).
(In this comparison, we have converted the UV emissivity per unit star
formation assumed by [25] to our value.) Our “initial” dust masses (with
$f_{M}=1$) are typically a few times larger than their reported masses.
However, the expected trends with halo mass are quite comparable. We do note
that this discrepancy also increases with redshift: at $z\sim 10$, [25]’s
model has dust masses more than an order of magnitude smaller than ours (if
$f_{M}=1$). We return to this comparison in more detail in section 5.1.

###  2.3 Dust Absorption

We next use simple models (similar to [46, 47]) to connect to observations
through two observables: the absorption of UV photons and emission in the IR
continuum.

In the UV, the wavelength-dependent optical depth along a line of sight is

$\tau_{\lambda}=\int_{0}^{R_{\mathrm{dust}}}\rho_{\mathrm{dust}}\kappa_{\mathrm{UV}}dr,$
(2.10)

where $R_{\mathrm{dust}}$ is the dust radius, $\rho_{\mathrm{dust}}$ is the
dust density, and $\kappa_{\mathrm{UV}}$ is the UV dust opacity.

The opacity depends on the choice of dust extinction law. As we focus on the
rest-UV regime, we adopt a simplified power law approach from [61], with

$\kappa_{\rm
UV}=\kappa_{1000}\left(\frac{\lambda}{10^{3}\text{\AA}}\right)^{\gamma_{k}},$
(2.11)

where $\gamma_{k}=-1$ is the fiducial power law adopted in their work and
$\kappa_{1000}=10^{5}\ \mathrm{cm^{2}\ g^{-1}}\approx 20\ \mathrm{pc^{2}\
M_{\odot}^{-1}}$ is the dust opacity evaluated at 1000 Å. This prescription is
consistent with the more sophisticated SMC dust extinction law presented by
[14], which is a carbonaceous-silicate grain model (and which appears to match
high-$z$ galaxies reasonably well; [62]). We note that the simple power law
estimate is reasonable for the UV continuum at 1500–1600 Å, but a more
sophisticated dust extinction law may be warranted to consider the effects of
specific atomic and molecular absorption lines (e.g. the 2175 Å “bump”), as
described in detail by [40]. To test the effect of particular extinction laws
on our model, we have estimated the dust optical depth (equation 2.10) with
both [13] and [14]’s results. We found that our model is not sensitive to the
choice of extinction laws, thus motivating the choice of a simple power law.

Because the integral in equation (2.10) extends over the entirety of the dust
distribution, we need to model the dust’s radial extent. We assume that the
dust radius is proportional to the scale radius of a disk galaxy embedded in a
dark matter halo. At high redshifts, the halo virial radius is by [63]

$R_{{\rm halo}}\approx 1.5\left(\frac{M_{\rm
halo}}{10^{8}M_{\odot}}\right)^{1/3}\left(\frac{10}{1+z}\right)\,\mathrm{kpc}.$
(2.12)

Assuming the gas in a galaxy is rotationally supported, the galaxy scale
radius is conventionally written as $R_{\rm gal}=(\lambda/\sqrt{2})R_{\rm
halo}$, where $\lambda$ is the spin parameter. We take $\lambda=0.028$ as a
fiducial value for a typical disk galaxy [64]. While this scale provides a
reasonable estimate for the gas and stellar distributions, we expect that dust
may be more extended– if it is ejected via feedback, for example. We therefore
insert a radius scaling parameter $f_{R}$

$R_{\mathrm{dust}}=f_{R}\frac{\lambda}{\sqrt{2}}R_{\mathrm{halo}},$ (2.13)

which we will later calibrate to observations.

To evaluate equation (2.10), we also need to determine the radial dust
distribution, which is unconstrained. We will consider three simple
morphologies below. But we note that, even after that distribution is
specified, the optical depth for each star’s light will vary depending on its
location within the galaxy. We are ultimately interested in the dust’s effects
averaged over all lines of sight. Therefore, we describe this via the galaxy-
averaged escape probability of UV photons $P^{\rm esc}$,

$L_{\rm obs}=P^{\rm esc}L_{\rm int},$ (2.14)

where $L_{\rm int}$ is the intrinsic luminosity produced by the stars and
$L_{\rm obs}$ is the apparent luminosity after dust processing for a distant
observer.

####  2.3.1 Homogeneous Mixture

We first consider a homogeneous, spherical mixture of stars and dust. This
morphology corresponds to the scenario where dust is retained within the
galaxy and uniformly mixed. In this case, equation (2.10) becomes

$\tau_{\mathrm{mix}}=\frac{3M_{\mathrm{dust}}}{4\pi
R^{2}_{\mathrm{dust}}}\kappa_{\mathrm{UV}}.$ (2.15)

Note that $\tau_{\rm mix}\propto f_{M}/f_{R}^{2}$ ; if only the absorption is
measured, the degeneracy between the dust mass and spatial distribution cannot
be broken. We will fix this combination of parameters by comparison to UV
observations in section 3. We can then use additional observations – dust
emission in the IR – to constrain the individual parameters in section 4.2.

For a homogenous mixture, [65] gives the escape probability of a UV photon as
a function of optical depth as

$P^{\rm
esc}_{\mathrm{mix}}(\tau)=\frac{3}{4\tau}\left[1-\frac{1}{2\tau^{2}}+\left(\frac{1}{\tau}+\frac{1}{2\tau^{2}}\right)e^{-2\tau}\right].$
(2.16)

####  2.3.2 Spherical Shell

Next, we consider a spherical shell surrounding the galaxy. This may occur if
feedback has ejected the dust to moderate distances from the galaxy. In that
case, equation (2.10) becomes

$\tau_{\mathrm{shell}}=\frac{M_{\mathrm{dust}}}{4\pi
R^{2}_{\mathrm{dust}}}\kappa_{\mathrm{UV}}.$ (2.17)

Note that this is three times less than $\tau_{\mathrm{mix}}$, as dust is
farther from the stars (and hence more spread out) than in the mixture
geometry. Note that we have the same $\tau_{\rm shell}\propto f_{M}/f_{R}^{2}$
relation. However, in this case, we have

$P^{\rm esc}_{\mathrm{shell}}(\tau)=e^{-\tau}.$ (2.18)

Although the characteristic optical depth is smaller than the homogeneous
case, the average escape probability of a shell decreases much more rapidly
with optical depth, because the stars are all buried within the shell.

####  2.3.3 Partially Covered Shell

The shell model is motivated by galactic feedback ejecting dust from the
central regions, but in reality, that process is highly anisotropic. We
therefore next generalize the shell model to one with “holes” where UV photons
can escape without any absorption. We define $f_{\mathrm{cov}}$ as the
fraction of the galaxy’s UV luminosity subject to dust absorption. (In a
simple model, this would correspond to the fraction of lines of sight that
intersect the dust, but because the stars may also be inhomogeneously
distributed we leave the definition more open-ended.) This may happen if
feedback is particularly strong in some directions and ejects the dust
entirely out of the galaxy’s host dark matter halo. The escape probability for
UV photons in this geometry is then

$P^{\rm esc}_{\mathrm{hol}}=(1-f_{\mathrm{cov}})+f_{\mathrm{cov}}[P^{\rm
esc}_{\mathrm{shell}}(\tau_{\mathrm{shell}})],$ (2.19)

where $\tau_{\rm shell}$ is defined as in equation (2.17).

Note that this geometry has an additional free parameter $f_{\rm cov}$ and
will therefore require additional assumptions when we return to it in section
4.3.

###  2.4 Dust Emission

Next, we estimate the IR emission from heated dust in order to compare it with
ALMA observations. We will compare the observed continuum fluxes around a
specific wavelength (near the [C II] emission line at 158 $\mu$m, as that is a
common window for observations). The total emission depends on the dust mass
and temperature. We adopt an analytic temperature estimate from [47], which
assumes that the dust is in radiative equilibrium with the absorbed starlight
and CMB heating. Then

$T_{\mathrm{dust}}=\left(\frac{L_{\star}^{\mathrm{abs}}}{C\kappa_{\mathrm{IR,0}}M_{\mathrm{dust}}}+T_{\mathrm{CMB}}^{\beta_{\mathrm{IR}}+4}\right)^{1/(\beta_{\mathrm{IR}}+4)},$
(2.20)

where the infrared absorption coefficient is given by a power law in
frequency,
$\kappa_{158}=\kappa_{\mathrm{IR,0}}(\nu/\nu_{0})^{\beta_{\mathrm{IR}}}$ with
$\kappa_{\mathrm{IR,0}}=16\ \mathrm{cm^{2}\ g^{-1}}$, $\nu_{0}=c/100\ \mu{\rm
m}$ is the normalization frequency, and $\beta_{\mathrm{IR}}=1.6$ as the IR
spectral slope. Both these fiducial values are taken from [66]’s Table 1 and
are appropriate for the Milky Way. Also, $T_{\mathrm{CMB}}$ is the CMB
temperature at a given $z$. The constant $C$ is

$C=\frac{8\pi
k_{B}^{\beta_{\mathrm{IR}}+4}}{c^{2}\nu_{0}h^{\beta_{\mathrm{IR}}+3}}\zeta(\beta_{\mathrm{IR}}+4)\Gamma(\beta_{\mathrm{IR}}+4),$
(2.21)

where $k_{B}$ is the Boltzmann constant, $c$ is the speed of light, $h$ is
Planck’s constant, $\zeta$ is the Riemann zeta function and $\Gamma$ is the
gamma function. The total absorbed UV stellar emission
$L_{\star}^{\mathrm{abs}}$ is computed by integrating the absorbed UV
luminosity across all frequencies, so

$L_{\star}^{\mathrm{abs}}=\int
L_{\mathrm{obs}}\frac{1-P^{\mathrm{esc}}(\tau_{\lambda})}{P^{\mathrm{esc}}(\tau_{\lambda})}d\lambda,$
(2.22)

where $\tau_{\lambda}$ is the wavelength-dependent optical depth. We normalize
$\tau_{\lambda}$ with the effective optical depth at 1600Å (which is uncertain
until we fix the scaling parameters) and assume that
$\tau_{\lambda}\propto\lambda^{-1}$ at other UV wavelengths, as approximated
by the [61] extinction law. To compute this integral, we integrate from the
Lyman break wavelength 912 Å to 6000 Å, though we find that the result is not
sensitive to any choice for the upper limit beyond 2800 Å. We note that
equation (2.20) includes heating from the CMB. Sometimes this effect is
ignored in estimates taken from observations; in that case, one can take the
limit $T_{\mathrm{CMB}}=0\mathrm{K}$ in equation (2.20).

This method assumes that the dust is at a single temperature, so that the
absorbed starlight is distributed evenly across the dust grains. [26] notes
that when the dust is optically thick, grains will be exposed to a range of
intensities so will reach a range of temperatures. They thus compute a
luminosity-weighted temperature, such that the dust population closer to the
source is hotter. We ignore this effect because our opacities are relatively
small; we have found that it has $\lesssim 10\%$ effect. However, such a
correction may be more important in a different morphology or very optically
thick medium.

With the temperature in hand, we can next predict the far-IR continuum flux.
To compare to the largest set of observations, we focus on the far-IR dust
continuum flux around the [C II] emission line – it is likely the strongest
line so is the target of many campaigns, and the dust continuum can be
conveniently estimated elsewhere in the band. Hereafter, we use $F_{158}$ as a
short-hand for the dust far-IR continuum flux in this window. To calculate
$F_{158}$, we adopt the “grey-body” treatment of [45], so

$F_{158}=\frac{1+z}{d_{L}^{2}}\,M_{\mathrm{dust}}\kappa_{158}\left[B_{158}(T_{\mathrm{dust}})-B_{158}(T_{\mathrm{CMB}})\right].$
(2.23)

where the $(1+z)$ factor results from the shift in the bandwidth, $d_{L}$ is
the luminosity distance, and $B_{158}$ is the Planck function at the dust
temperature at the emitted frequency $\nu$. We note again that some authors do
not subtract the CMB contribution, so the last term may not be necessary when
comparing to observations.

##  3 Empirical Dust Attenuation Calibration

Our simple dust model has two free parameters, characterizing the dust’s
spatial extent (through $f_{R}$) and its mass (parameterized by $f_{M}$). In
this section, we will describe how we can use existing UV observations to
place a joint constraint on these parameters.

###  3.1 UV Attenuation Calibration

Ideally, one would estimate the UV dust attenuation via SED modelling, where
the level of extinction is obtained in combination with the stellar
parameters. However, such modelling is not yet available for large numbers of
high-$z$ galaxies, so instead, we follow a common approximate method. This
method assumes that the UV slopes $\beta$ of these galaxies depend on the dust
properties (as dust reddens the spectrum). We use observed estimates of these
slopes together with an empirical relation from [16] that relates the UV slope
to the level of extinction.

First, we require the UV slopes of “typical” high-$z$ galaxies as a function
of magnitude and redshift, which we obtain from two observed samples detailed
later. For simplicity, we assume that the average UV slope
$\langle\beta\rangle$ is linear in the observed magnitude, with the slope
$d\beta/dM_{\mathrm{UV}}$ and intercept $\beta_{M_{0}}$ varying linearly with
redshift as well. Thus we have

$\langle\beta\rangle\left(M_{\mathrm{UV}}^{\mathrm{dust}},z\right)=\frac{d\beta}{dM^{\rm
dust}_{\mathrm{UV}}}(z)\left[M_{\mathrm{UV}}^{\mathrm{dust}}-M_{0}\right]+\beta_{M_{0}}(z),$
(3.1)

where $M_{0}$ is a reference magnitude, and we assume that

$\displaystyle\beta_{\mathrm{M_{0}}}(z)$
$\displaystyle=\frac{d\beta}{dz}z+\beta|_{z=0},$ (3.2)
$\displaystyle\frac{d\beta}{dM^{\mathrm{dust}}_{\mathrm{UV}}}(z)$
$\displaystyle=\frac{d}{dz}\left(\frac{d\beta}{dM^{\mathrm{dust}}_{\mathrm{UV}}}\right)z+\left.\frac{d\beta}{dM^{\mathrm{dust}}_{\mathrm{UV}}}\right|_{z=0}.$
(3.3)

Table 1: The best-fit parameters of $\beta_{\mathrm{M_{0}}}$ and $\frac{d\beta}{dM^{\mathrm{dust}}_{\mathrm{UV}}}$ to [67] & [6] as linear functions of redshift, displayed as purple curves in Figure 2. The first column represents the slope of the linear function, and the second column represents the offset evaluated at $z=0$.  | slope | offset $(z=0)$
---|---|---
$\beta_{\mathrm{M_{0}}}$ | $-0.081\pm 0.016$ | $-1.58\pm 0.115$
$d\beta/dM^{\mathrm{dust}}_{\mathrm{UV}}$ | $0.012\pm 0.006$ | $-0.216\pm 0.040$
Figure 2: The redshift evolution of the UV slope fit parameters,
$\beta_{\mathrm{M_{0}}}(z),\frac{d\beta}{dM^{\rm dust}_{\mathrm{UV}}}(z)$,
measured from the [67] (blue squares) and [6] (orange circles) data sets. We
excluded [67]’s measurement at $z=8$ (open square) from the fit since the
corresponding slope was fixed. We also show measurements from [68] at $\langle
z\rangle\approx 9$ and 11 (green diamonds; the latter was also estimated
without varying the slope). These are not included in our fit. The purple
curve exhibits the joint best fit to [67] and [6] with parameters presented in
Table. 1, whereas the dashed blue and orange curves are the best fit for the
individual datasets. We align all the data to have $M_{0}=-19.5$.

To estimate these parameters, we use galaxies from [67], who measured slopes
from over 4000 HST sources in deep surveys at $z=4$-$8$, and [6], who measured
slopes over $z=6$-$12$ with a large set of sources from the JWST Advanced Deep
Extragalactic Survey (JADES). We summarize the fit parameters in Table. 1 and
plot the resulting relation in Figure 2. We show the measurements at each
redshift with the data points, and our best-fit linear relations for the
parameters via the solid line. The dashed lines show the individual fits to
the different datasets. Note that we have excluded the $z=8$ data point from
[67] for our fit, because its properties were estimated after holding the
slope fixed. The qualitative evolution of the UV slopes is clear, although the
precise form is only weakly constrained.

Figure 2 also includes measurements from [68], who find a steeper evolution in
the zero-point than [6]. This is likely due to the sparseness of the data at
$z\gtrsim 10$. Given this, and the potential for contamination in photometric
surveys [69, 70, 71] and cosmic variance [8], we do not include additional
measurements. However, in Appendix A, we do compare our inferred relation to
other recent JWST measurements, finding reasonable consistency (albeit with a
large scatter amongst individual sources).

We must now transform these UV slopes into estimates of the dust extinction.
To do so, we use the relation from [16], which is an empirical relation
(calibrated to local starburst galaxies) that links the observed slope to its
corresponding infrared excess, which parameterizes the amount of dust emission
(ultimately due to absorbed starlight). This can therefore be converted to the
total UV emission in the absence of dust, $M_{\mathrm{UV}}^{\mathrm{int}}$ by
[42]

$M_{\mathrm{UV}}^{\mathrm{int}}=\left(1-C_{1}\frac{d\beta}{dM_{\mathrm{UV}}^{\mathrm{dust}}}\right)M_{\mathrm{UV}}^{\mathrm{dust}}-C_{0}-C_{1}\beta_{M_{0}}+C_{1}\frac{d\beta}{dM_{\mathrm{UV}}^{\mathrm{dust}}}M_{0},$
(3.4)

where we take constant values $C_{0}=4.43,C_{1}=1.99$ for absorption at 1600 Å
from [16]. We also require that $M_{\mathrm{UV}}^{\mathrm{dust}}\leq
M_{\mathrm{UV}}^{\mathrm{int}}$ to prevent unphysical negative attenuation
(which can happen as a result of our simple UV slope parameterization). This
becomes important for faint galaxies with negligible dust production. We refer
interested readers to the detailed derivation of equation (3.4) in [42]’s
section 3.2.1. Also, we note that we use luminosities evaluated at 1500 Å,
while the [16] relation uses the absorption at 1600 Å. We convert between
these with the power law $L_{\lambda}\propto(\lambda/\lambda_{0})^{\beta}$.

Figure 3: The escape probability of UV photons at 1600 Å in our model. The
upper panel shows the empirical probabilities calculated with the IRX-$\beta$
relations from [16] (solid curves) and [21] (dashed curves), respectively,
whereas the lower panel shows the probabilities calculated from our models
assuming $f_{d}=1$. The violet star symbol represents the brightest galaxy in
our JWST data compilation in Appendix A, 817482 from the COSMOS field, from
[72] at $z=9.89$. It has $P_{\mathrm{esc}}\approx 0.72$ according to the [16]
relation at $z=10$.

This procedure provides the intrinsic luminosity, and hence the probability
that UV photons escape the system to the observer (via equation 2.14). We show
the results in Figure 3. The attenuation increases towards lower redshifts and
towards higher luminosities, as one would naturally expect. Interestingly, for
much of the Cosmic Dawn, the attenuation is negligible except in very luminous
galaxies: even at $z\sim 5$, $P_{\rm esc}\approx 1$ for
$M_{\mathrm{UV}}\gtrsim-18$ 222Hereafter, we use $M_{\rm UV}$ to denote the
dust-processed or observed magnitude $M_{\rm UV}^{\rm dust}$ for conciseness
unless otherwise noted., which corresponds to a halo mass $\sim
10^{11}M_{\odot}$. 333This result motivates our temperature treatment in
section 2.4, where we assume a single mean dust temperature in the optically
thin regime.  As shown in appendix A, most observed systems at very high
redshifts are near this limit, so the overall dust correction to the universal
star formation rate density is expected to be modest during and before
reionization (see section 5.5). To provide some intuition, the star in Figure
3 corresponds to the brightest high-$z$ galaxy in the COSMOS field (see the
Appendix), ID 817482, at $z\approx 9.89$ [72]. In this case, the escape
fraction of star formation is about $72\%$ according to our model.

In the next subsection, we will compare these empirical results to our model.
But before proceeding, we must note some caveats in adopting this empirical
dust law. First, it is not clear that the IRX-$\beta$ relation of [16] can be
applied to higher redshifts. Observations at $z\approx 3$ have found
consistency with [16] for faint galaxies, although more massive galaxies with
more dust may present problems. At $z\approx 5$, [21] finds the dust continuum
in many galaxies to be smaller than that predicted by the [16] relation. One
possible explanation is that the dust and star-forming regions may be
spatially segregated in the high-$z$ regime so that the dust is illuminated by
only a fraction of the UV light [26, 12, 28]. As the IRX-$\beta$ relation
assumes the UV and IR emitting regions to be co-spatial, extrapolation from
low-$z$ may result in an overestimation of dust reddening. We return to a
detailed discussion about the current IRX-$\beta$ relations in section 5.2.

To test the importance of the particular IRX-$\beta$ relation, we have also
estimated the escape probability using the empirical relations from [18, 19,
20, 21]. In Figure 3 we also show the photon escape fractions inferred from
the IRX-$\beta$ relation measured by [21] from a (small) sample of higher
redshift galaxies ($z\sim 5$). The results are qualitatively similar to our
fiducial model, but they show a gentler dependence on halo mass (or
luminosity). Given that these faint galaxies are only now beginning to be
observed by JWST, we stick with the widely-used [16] relation for consistency
with past work. However, for global properties where the faint galaxies play a
significant role (e.g. SFRD; section. 5.5), we will also show an estimate from
[21]’s IRX-$\beta$ relation.

###  3.2 Inferred Dust Properties

With the empirically-derived dust attenuation in hand, we can now begin to
calibrate the free parameters in our model. We will focus on the shell and
mixture morphologies for now, which have parameters for the dust spatial
extent ($f_{R}$, defined relative to the disk scale length) and dust mass
($f_{M}$, defined relative to the total amount of dust expected to be
created). However, because the escape fraction only depends on the dust
optical depth in these simple models, we can only constrain their combination,
$f_{d}\equiv f_{M}/f_{R}^{2}$, with the UV properties.

The lower panel of Figure 3 shows the predicted UV escape fractions in our
shell and mixture models assuming $f_{d}=1$. It is immediately obvious that
the models predict far more absorption than is observed, with the model
optical depths one or two orders of magnitude larger than the empirical
results. This implies that either the dust is more spread out (so that each
line of sight has a smaller optical depth) or that the galaxies have less dust
than naively expected (or both). This appears to be a generic conclusion of
galaxy formation models at these redshifts (e.g., [61, 25, 29, 11]). Note that
our model also predicts an opposite redshift-evolution trend compared to that
of the observations, with the optical depth increasing with redshift in the
models. This is because star formation (at a fixed halo mass) is slightly more
efficient at high redshifts and (more importantly) because galaxies are more
compact.

Figure 4: Contours showing the loci of possible solutions $(f_{M},\,f_{R})$
that gives the empirically-derived $f_{d}$ across the modelled galaxy
populations by our UV calibration at various redshifts. The solid curves are
computed for a halo mass $10^{9}M_{\odot}$ and the dash-dotted curves take
$10^{11}M_{\odot}$. All spaces above each contour are also possible solutions
for galaxies that require a smaller $f_{d}$.

Figure 4 shows the loci of $f_{d}$ values that bring our model in agreement
with the UV observations at a variety of redshifts (assuming that the
IRX-$\beta$ relation from [16] is accurate). As shown, the required $f_{d}$
value decreases towards higher redshifts and lower halo masses as dust
extinction decreases. Because the mass and radius scaling parameters are
degenerate, a range of solutions is possible at any given redshift. To isolate
their values, we require additional observations of the dust emission, which
we will explore in the next section.

##  4 Constraining the Model with Dust Emission

With its ability to detect the rest-frame far-IR dust continuum of high-$z$
galaxies, the Atacama Large Millimetre Array (ALMA) has become our heaviest
artillery for understanding dust in early cosmic times. Specifically, ALMA
surveys targeting the [C II] line can measure (or constrain) the nearby dust
continuum. In this section, we will use the IR emission to break the
degeneracy between the dust parameters.

###  4.1 The ALMA Comparison Sample

We focus on two ALMA surveys: the ALMA Large Program to INvestigate [C II] at
Early times [ALPINE, 30] carried out in Cycle-5 and Cycle-6, and the
Reionization-Era Bright Emission Line Survey [REBELS, 34] carried out in
Cycle-7. With the shared objective of providing large statistical samples of
far-IR continuum observations of UV-selected bright galaxies, ALPINE observed
the [C II] continuum emission of 23 galaxies in the redshift range
$4.4\lesssim z\lesssim 5.9$, while REBELS observed 14 galaxies with a mean
redshift around $\langle z\rangle\approx 7$. Complementary studies from both
surveys supplement UV and optical properties of the observed galaxy
populations, presented by [31] for ALPINE and [34] and [26] for REBELS,
respectively.

Although the analyses of these surveys present many properties of the observed
galaxies, we must be careful comparing them to our model. For example, the
star formation rate of these galaxies inferred from their UV luminosity
depends on assumptions about the IMF, metallicity, etc., while the reported
stellar masses depend on the details of the SED modelling. We will generally
report our results as a function of the SFR, because it is more easily related
to observables. In order to ensure a meaningful comparison, we adopt the
conversion in equation (2.8), applied to the reported luminosities. For
ALPINE, we convert the UV magnitudes at 1500 Å from [31] to the SFR. For
REBELS, we convert [26]’s total SFR to the observed-UV SFR using their derived
transmissivity factor $T_{\mathrm{1500}}$, while also replacing the
SFR-$L_{\mathrm{UV}}$ conversion factor in their equation (1) with our choice
for $\mathcal{K}_{\rm UV}$.

One additional complexity is that stellar mass inferences depend on
assumptions about the stellar populations and on the photometric library used
for the measurement. The assumptions of these codes typically do not match the
results of high-$z$ galaxy models (e.g., [9]). The SED-fitting process is
complex enough that we cannot easily translate between the inferences and
models, so we will not attempt to do so here. Instead, we will mostly focus on
comparisons as a function of SFR, which are at least transparent to convert.
When we report stellar masses, we will simply use the results of the ALPINE
and REBELS groups without modification.

However, we note that our minimalist model assumes a one-to-one correspondence
between halo mass, stellar mass, and SFR. Real galaxies are of course
variable, both because of the differing halo environments and because of
temporal variations in individual systems. Because the stellar mass is a
cumulative measure, it is relatively stable when averaging over this
variability – while the SFR against which we will compare is much more
sensitive to it. There are indications that galaxies at $z\gtrsim 8$ are
extremely bursty [8, 9, 73], so we will return to this issue in section 5.3.
For now, we note that we should not take the comparison too seriously for any
individual system.

Figure 5 provides a rough comparison between the model and observed samples in
terms of their “star-forming main sequence” relating SFR and stellar mass. The
model results are shown by the solid curves, while the individual galaxies in
the surveys are shown by the circles with error bars. Visually, the ALPINE
data are more or less consistent with the model, but the REBELS points lay
significantly above the models. This is confirmed by the dashed curves: these
are fit to the data points, assuming that they follow the same power-law trend
as the model but varying the normalization. The REBELS fit is driven higher by
the three leftmost points, which have large UV luminosities (and hence SFRs)
but small inferred stellar masses. These could be vigorous starbursts in small
galaxies, or they could have “hidden” mass components from old stellar
populations. (We note that the SED-based stellar masses did not have JWST data
available, so are mostly insensitive to such old stars.)

Figure 5: The star-forming main sequence at $z=5$ and 7 as seen in our models
(solid curves) and the ALMA surveys. The orange points are from the REBELS
survey [26] at $z\approx 7$, corrected to our assumptions about the stellar
populations. The green points are from the ALPINE survey [31] at $z\approx
4.5$–5.5. The dashed curves show fits to these data sets if we constrain the
slope to follow that expected in the model. The darker star symbols show
results from stacking analyses that include sources undetected by these
surveys [21, 28]. We also show the conversion between stellar mass and UV
magnitude in our model along the top axis.

However, the individual survey detections are biased, as they only include
galaxies with bright enough dust emission to be detected. Both surveys also
have significant numbers of undetected sources. The population properties
must, therefore, be constrained through stacking analyses that include these
fainter sources. The dark stars in Figure 5 show the results of such analyses
for both surveys [21, 28]. Again, ALPINE is reasonably consistent with our
model, albeit over a limited range in stellar masses, but we see that now
REBELS is also in reasonable agreement. We also note that both surveys find a
much flatter relationship between $M_{\star}$ and SFR than our model provides;
whether this is a selection effect is unclear, given the limited dynamic range
in stellar mass.

###  4.2 Implications for the Dust Properties

We are now able to leverage the ALMA dust emission measurements to further
constrain our models. In section 3, we used the UV attenuation inferred from
the IRX-$\beta$ relation to estimate the overall optical depth, and from that,
we placed a joint constraint on the mass and radius scaling factors. The dust
emission also depends on the dust mass, so adding it can break the degeneracy
between mass and radius.

Figure 6: Dust scaling factors in our model, inferred through the combination
of UV attenuation and dust emission. The green points correspond to ALPINE
galaxies (using the dust continua reported by [30], while the orange points
are from the REBELS survey [26]. The dashed curves show power-law fits to each
survey’s results, while the solid curve shows a power-law fit to their
combination. _Left:_ The mass scaling parameter $f_{M}$. _Right:_ The radius
scaling parameters in the shell (top) and mixture (bottom) morphologies. In
each of these panels, we also show the resolution of the individual REBELS
detections as the solid red curve ($\approx 6.3\mathrm{kpc}$) and the
estimated sizes from stacking sources from both surveys with stars [27].

Figure 6 shows our main results. The left panel shows the inferred value of
$f_{M}$ for each galaxy (which is independent of the assumed dust morphology
in the context of our simple model), while the right panels show the inferred
values of $f_{R}$ for the shell and mixed dust morphologies. Interestingly, we
find that relatively high dust masses – with $f_{M}\sim 1$ – are required to
produce the overall dust luminosity. This implies that the systems are able to
retain most of the dust they produce. However, we have already seen that the
optical depth is much lower than expected in this case, so we also need
$f_{R}\sim 10$ for each model. The physical dust radius is then $\sim
f_{R}\lambda R_{\rm vir}\sim 0.3R_{\rm vir}$. Thus, the dust would mostly be
mixed through the circumgalactic media of these halos.

We also see some potential trends in the data. The REBELS points tend to have
slightly smaller $f_{M}$ and larger $f_{R}$ than the ALPINE points, but the
differences are small and within the scatter of the points at each redshift.
For a simple estimate, we therefore take a redshift-independent average fit.

The trends with SFR are more obvious. In contrast to intuition from feedback-
based models, $f_{M}$ appears to _decrease_ as the star formation rate
increases (as previously found by [25, 29]). If we follow the minimalist
galaxy formation model and assume that the SFR increases monotonically with
halo mass, this means that small halos would be more likely to retain their
metals than massive halos – contrary to the expectations of galaxy evolution
models based on feedback. One possibility is enhanced dust destruction in
massive systems, possibly due to the increasing amount of SNe shocks [58]. On
the other hand, $f_{R}$ also decreases with SFR – so that the dust that
remains is closer to the galaxy’s central regions in more massive systems.

We must be cautious in interpreting these trends however, as we would expect
similar qualitative results from burstiness: galaxies that have undergone
recent starburst episodes that have now faded (so lay at small SFRs) may have
residual dust from their previous event (so appear to have a large $f_{M}$).
Going forward, we therefore consider two cases: one in which we take a single
average value for each parameter (mass- and redshift-independent) and one in
which we assume a power-law dependence on SFR. Fits for the latter are shown
in Figure 6. The redshift-independent fit (where we vary both the
normalization and power-law slope) is given by

$f_{M}(\dot{m}_{\star})=42.4\,\left(\frac{\dot{m}_{\star}}{M_{\odot}/{\rm
yr}}\right)^{-1.17},$ (4.1)

If we instead assume a constant value, we find the median $f_{M}$ value across
the galaxy sample to be $\approx 0.75$.

Our model therefore requires that the dust be quite spread out across the
circumgalactic medium, a result that is in tension with recent observations.
[12] found that the REBELS galaxies typically have dust radii below the ALMA
resolution, which was approximately $R_{\mathrm{max}}\approx 6.3\
\mathrm{kpc}$ for their survey (shown by the red curve in the right panels of
Figure 7). The radii required by our models are dangerously close to this
limit, with some even above it. Moreover, [27] have measured the extent of the
dust emission in these sources with a stacking analysis, shown by the stars in
Figure 6. Stacking UV and dust emission is quite challenging, as it requires
the sources to be spatially aligned, whereas the dust is often offset (by an
unpredictable amount) relative to the UV emission in any given source [26,
12]. Nevertheless, the observed dust radii are several times smaller than
expected in our model. This motivates a more complicated geometry that we
consider in the next section.

Figure 7: Comparison of model and observed dust continuum fluxes for the
ALPINE and REBELS galaxies. The solid curves show the model fluxes using the
power law fit of equation 4.1, while the dashed curves take $f_{M}=1$. The
green points are observed by ALPINE with $4.4\lesssim z\lesssim 5.9$, while
the blue points are observed by REBELS with $\langle z\rangle\approx 7$.

Before turning to that, however, we first examine a couple of other aspects of
the model framework. Figure 7 compares the resulting model fluxes (solid
lines) to the observations (points with error bars). Because we have taken a
redshift-independent fit, the model underestimates the ALPINE sources and
overestimates those of REBELS sources by a modest amount. This also
illustrates another potential limitation of the surveys: because they are
limited by the survey depths, they cannot determine how wide the distribution
of dust luminosities is for a given SFR. The dotted lines show our model with
$f_{M}=1$, which decreases at small SFR – unfortunately, the surveys would not
be sensitive to such a decline.

Figure 8: Dust temperatures (neglecting the CMB) in our model. The solid
curves show the ALMA-calibrated temperature with the $f_{M}$ fit provided by
equation 4.1, while the dotted curves show the temperatures assuming
$f_{M}=1$.

As described above, our model also estimates the dust temperature. This has
not been directly measured for these sources, because they have only been
observed over a very narrow range of wavelengths, but we show the model
results in Figure 8; here we take $T_{\rm CMB}=0\ {\rm K}$ in equation (2.20),
as often done when interpreting observations. The true dust temperatures are
somewhat smaller, because the CMB contributes to the emission as well. The
solid curves show our model with the SFR-dependent fits, while the dotted
curves are with constant parameters. Since $f_{M}$ decreases towards large
${\rm SFR}$, the calibrated dust temperature is higher than the original on
the luminous end and lower on the faint end. However, for the luminosity range
that concerns us most, the temperatures do not change much.

###  4.3 Gaps in the Dust Distribution?

In the previous subsection, we found that our model required quite extended
dust distributions in order to account simultaneously for the large dust
luminosity and modest extinction. This occurs in both the limiting cases of
dust mixed uniformly with the stars and dust expelled to a shell outside of
the central galaxy, but it is in tension with observations [27, 12] of the
radius of the dust emitting region. One potential resolution is to drop the
assumption of spherical symmetry in both our shell and mixture morphologies:
we now turn to the partially covered shell model, in which dust appears only
along a fraction $f_{\rm cov}$ of the lines of sight.

Such a model does seem reasonable in light of other observations. [12] also
found a spatial offset between the UV and IR emission peaks in the REBELS
survey, suggesting along with [28]’s stacking analysis that the UV and IR
emitting regions may not be co-spatial. [74] concludes the same as a part of
their specific SFR analysis joining UV and dust emission. In our “partial
shell” model, this could be interpreted as a shell with a small covering
factor and large radius, so that the dust region is confined to a small area
displaced from the stars. More generally, this could occur if dust coincides
with only a fraction of the star-forming regions.

Figure 9: The covering fraction $f_{\mathrm{cov}}$ as a function of UV
magnitude at various redshifts in our partially covered shell model. The
curves are calculated from a scenario with the median $f_{M}$ value of 0.75
from the previous fit and a fiducial scaling $f_{R}=3$, motivated by [27]’s
stacking analysis.

However, a morphology with “holes” introduces a new free parameter into the
system, which is hard to constrain given the existing observations. We
therefore take a simplified model in which: _(i)_ we take the median $f_{M}$
value of $0.75$ from Figure 6 and _(ii)_ we set $f_{R}=3$ to match the stacked
measurements of [27]. We then fix the covering fraction by equating the left
side of equation (2.19) to the empirically-determined UV escape probability.
We assume that all these parameters are mass-independent. Note that the
covering fraction does not affect the dust emission, because we hold its mass
fixed (see equation 2.17), implicitly assuming that the dust remains in the
system but is confined to a smaller region. Our results are not terribly
sensitive to these choices, so long as we use values close to those suggested
by the previous section.

We show the resulting covering fractions in Figure 9. The small dust
attenuation in faint galaxies implies that they are nearly clear of dust, but
the covering fraction becomes large in luminous systems. More luminous
galaxies have $0.2\lesssim f_{\rm cov}\lesssim 0.7$, implying that the dust is
more uniformly distributed in such systems. It seems reasonable that their
longer star formation histories and larger SFRs mean less stochastic feedback
and more time for metal mixing. Conversely, dust in smaller galaxies could be
heavily impacted by bursty feedback, resulting in “holes” in the feedback-
affected region and dust clumps in the unaffected region.

##  5 Discussion

Here we examine some other aspects of the dust model and its consequences.

###  5.1 Comparison to Past Work

Other groups have also considered the dust content of high-$z$ galaxies. Here
we compare to a few particularly relevant studies. [25] used their DELPHI
semi-analytic model to estimate dust properties. They calibrated their star
formation parameters to the observed UVLF and stellar mass function (using HST
surveys). They then have physically motivated prescriptions for dust
production, destruction, and escape via feedback, following its abundance in
each system over time; while their dust production parameters are comparable
to ours, we do not explicitly model destruction or escape. Figure 1 compares
our dust mass estimates to theirs (dash-dotted line). By comparing to ours, we
find that their model corresponds to an effective $f_{M}\sim 0.3$. They
compare to REBELS sources and find that their model underpredicts the dust
mass inferred from infrared emission measurements, which we also find for the
dominant low-mass galaxy population. Similarly, [29] find that their own semi-
analytic dust model underpredicts the emission, arguing that the tension could
be alleviated by a top-heavy IMF (amongst other possibilities).

[26] and [75] have also estimated dust masses and temperatures from the REBELS
sample. [26] uses more sophisticated dust distributions and radiative transfer
models, including a disk geometry and allowing for the dust to be at a range
of temperatures. They use the $\beta$ slope and stellar population synthesis
to break the degeneracy between dust mass and radius. Interestingly, their
dust radii derived from $F_{158}$ are consistently lower than ours but their
dust masses are consistent with ours. Thus, like our model and those of [25,
29], they also find a high IR-to-UV ratio, concluding that this could possibly
suggest a non-uniform dust morphology. To compare their results to our
temperatures, we convert their luminosity-weighted temperature
($T_{d}^{\prime}$) in Table 2 back to the mean temperature without the CMB
effect ($\bar{T}_{d}$) described by their equation (15). (We remind the reader
that our model finds that most high-$z$ galaxies are in the optically thin
regime, so this has only a modest effect.) After this conversion, [26] finds
temperatures of $30$–$60\ \mathrm{K}$ in the observed sources. Meanwhile, [75]
find $40$–$58\ \mathrm{K}$. Even with rather different models, the predicted
temperatures all roughly agree; this is not too surprising given that the dust
luminosity is so sensitive to temperature.

[76] used an even simpler approach than ours to estimating dust continuum
fluxes, using a fixed star formation efficiency, assuming that all halo gas is
retained, and imposing a fixed metallicity and dust-to-metal ratio. They take
a much more sophisticated model for the dust properties, however, and then
predict the far IR luminosities of very high-$z$ sources. Their fiducial model
is comparable to ours, though they show that the dust emission can vary
significantly depending on their assumptions about its microscopic properties.

###  5.2 The IRX-$\beta$ Relation

Figure 10: The effect of the assumed IRX-$\beta$ relation on the UV photon
escape probability at $z=5$ and 10. The solid curves are calculated from [16]
derived from local starburst galaxies. The dotted curves use an SMC-like
extinction curve calibrated to galaxies at $z\sim 2$–3 [20]. The dashed curves
use a “grey attenuation” curve calibrated to galaxies at $z\sim 2$–3 [19]. The
dash-dotted curves are estimated from [21] using the ALPINE survey at $z\sim
5$.

A key part of our framework is the empirical calibration of the UV
attenuation. However, this quantity is not measured directly in the existing
samples but is inferred from the UV spectral indices, based on a calibration
to a sample at much lower redshifts. As noted in section 3, recent work shows
that this calibration may break down at high redshifts [21]. To further
examine the effects of particular IRX-$\beta$ relations on our model, we show
in Figure 10 how different IRX-$\beta$ calibrations affect the UV photon
escape probabilities. The solid lines show our fiducial model (from local
starburst galaxies; [16]), the dashed and dotted lines show estimates from
galaxies at $z\sim 2-3$ [19, 20] that adopt two types of extinction curves
(the “grey” curve and the “SMC-like” curve, respectively), and the dash-dotted
curves show an estimate from ALPINE galaxies at $z\sim 5$ [21].

Interestingly, our fiducial model finds the steepest dependence of attenuation
on galaxy luminosity: all the models calibrated to higher redshift galaxies
have gentler slopes, with the flattest resulting from calibration to ALPINE
galaxies [21]. The estimates at $z\sim 2$–3 are in between these limits. The
gentler slopes mean that small galaxies still have a modest amount of dust,
though it is still a very small effect on the overall UV luminosity. We re-run
the $f_{R}$ analysis using [21] and show the results in Figure 11. Note that
the $f_{M}$ results are independent of the IRX-$\beta$ choice (because they
are derived from the IR emission), so we only show the change in $f_{R}$. It
is obvious that a gentler IRX-$\beta$ slope increases $f_{R}$ compared to our
fiducial model, increasing the tension with the dust radii inferred from
stacking [27].

Figure 11: Estimated radius scaling parameters if we use [21] to infer the
level of UV attenuation (c.f. the right panels of Fig. 6, which used the
relation from [16], shown here by the dash-dotted curves).

We note as well that the recent study of [68] found a steep evolution of the
UV slope at $z\gtrsim 9$ (Figure 2). Along the same lines, [20] has suggested
the intrinsic UV slope could become bluer towards higher redshift. If so, this
would compromise the extrapolation of any of these results to very high
redshifts during the reionization era. A full understanding of the dust will
therefore hinge on a better understanding of the intrinsic UV slopes of these
stellar populations.

###  5.3 Dust and Bursty Galaxies

Another limitation of our model is the assumption of a one-to-one
correspondence between halo mass and star formation rate. In reality,
starbursts appear to be a crucial part of high-$z$ galaxy evolution [8, 9].
The ultimate cause of this burstiness is not clear, but one possibility is
that the short dynamical times of high-$z$ galaxies allow star formation to
“overshoot” the expectations from feedback regulation [64, 59, 77]. The excess
feedback energy from the runaway star formation can then shut down that
process for a period by driving gas (and presumably dust) outside of the host
halo, until the cycle begins anew. This process can occur in small galaxies
[59] and large ones [77].

As a result, the mapping between star formation and stellar mass can have
large fluctuations, which we have thus far ignored. In this section, we will
estimate the effects of such scatter in a simple phenomenological way. We
assume that the total dust mass is unaffected by burstiness, because it
depends on the stellar mass in our model (which is integrated over all past
episodes of star formation so is not as strongly affected by the current SFR).
However, the instantaneous luminosity – and hence the energy input into the
dust phase – does depend on the current SFR.

For a toy model, we assume that the intrinsic luminosity of galaxies varies by
a fixed amount of $\pm 1$ or $2$ magnitudes. The former is comparable to the
scatter expected in supernova-delay induced burstiness [59], while the latter
is close to the required level to explain the bright end of the luminosity
function at high redshifts [8, 9]. We then use this luminosity in equation
(2.20) to compute the dust temperature, from which the dust emission follows.
(Note that here we assume that the fraction of absorbed starlight does not
change from the fiducial calculation; see below.)

Figure 12: The burstiness-induced scatter in [C II] continuum flux at various
redshifts as a function of $\mathrm{SFR_{UV}}$. The solid curves are the
average continuum flux. The shaded region highlights the scattering induced in
these quantities at a given redshift with (1, 2) magnitudes. The darker shades
are the 1 magnitude scatters, and shallower shades are the 2 magnitude
scatters. The green circles are the ALPINE results from [30]. The blue circles
are the REBELS results from [26].

Figure 12 shows the resulting spread in the observed flux for several
redshifts. We find that a moderate amount of burstiness can easily explain the
observed scatter in the REBELS sample. The ALPINE galaxies show a wider range,
however. We note that scatter of $\pm 1$ mag causes a scatter in the dust
temperatures of $\pm 5$ K as well.

Here we have assumed for simplicity that the dust spatial distribution is
unaffected by burstiness. If the dust distribution is primarily a function of
the total stellar mass, then our UV-calibration procedure would imply that the
escape fraction of UV photons also has scatter at each UV luminosity. Also, we
note that it may actually help explain gaps in the dust morphology, as we
appear to require (see also [12, 27]). The increased radiation pressure during
a starburst can drive dust away from the stars, segregating it from the UV
emission, especially in small systems [78, 11].

###  5.4 Dust and the UVLF

One of the key surprises of JWST’s early surveys is the unexpected abundance
of luminous systems at $z\gtrsim 10$. While burstiness provides one potential
explanation, dust may contribute as well [10, 11]. The key idea is that the
abundance of UV-luminous systems at $z\sim 7$ is suppressed by their
relatively large amounts of dust, but that at higher redshifts the amount of
dust decreases, partly because it may be ejected by strong feedback [11]. Thus
if we normalize our model parameters to match the low-$z$ data, the galaxy
abundance will decrease less than expected at higher redshifts because the
decreasing dust obscuration partially counteracts the decrease in the halo
mass function.

Figure 13: The effect of dust on the high-$z$ UVLF. The dotted curves show the
UVLF computed from [49], which ignores dust. The solid curves show the
_intrinsic_ UVLF without dust-processing, assuming the IRX-$\beta$ relation
from [16]. The dashed curves are LFs recovered from the _intrinsic_ by
applying the [16] IRX-$\beta$ relation. Note that the effects of dust
obscuration are quite modest by $z\sim 10$.

We now consider this explanation in the context of our model. Figure 13 shows
two sets of luminosity functions: one in which we ignore dust entirely (and is
thus equivalent to the minimalist model of [49] as the “observed” LFs; dotted
curves) and one in which we include our dust model with the IRX-$\beta$
relation from [16]. Here we have obtained the intrinsic LF at $z\sim 7$ by
correcting to the total SFR using the [16] relation. The solid curves show
this intrinsic LF from $z=7$–10; note that this is _not_ the observable UVLF,
because it includes star formation hidden by dust. The dashed curves then show
the expected LF after dust attenuation is included.

We see that, as expected, dust attenuation is significant at the bright end at
$z\sim 7$, but by $z\sim 10$ it has a much smaller impact. Importantly, the
dashed and dotted curves are near each other at all redshifts – that is, the
overall redshift evolution is comparable when we ignore and include dust. We
do not find that it can account entirely for the overabundance of bright
galaxies at high redshifts, although it does help to some extent. This
emphasizes how dust evolution can complement other explanations for the excess
of sources, and it may even be physically associated with burstiness (as
described in the previous subsection).

###  5.5 The Fraction of Total Obscured SFR

So far, we have focused on dust in individual galaxies. We now shift to the
cumulative effect of dust across the entire galaxy population by considering
the star formation rate density (SFRD), which is given by

$\rho_{\rm tot}(z)=\int dm\,\dot{M}_{\rm\star,tot}(m,z)n(m,z).$ (5.1)

where the integration variable $m$ is the halo mass and $n$ is the halo
comoving number density in the mass range $(m,m+dm)$ at redshift $z$. We
calculate both the total SFRD and the obscured SFRD $\rho_{\rm obs}$, so that
the obscured fraction is $f_{\rm obsc}=\rho_{\rm obsc}/\rho_{\rm tot}$. We
calculate this quantity for the entire halo population that is able to form
stars. We suppose that the lowest mass limit for a halo to form stars is given
by its virial temperature $T_{\rm vir}\gtrsim 10^{4}{\rm K}$. (We take the
virial temperature formulation from [79]’s equation (26), which roughly gives
a halo mass limit of $10^{8}M_{\odot}$.) We also calculate the SFRD for a
subset of massive halos (with a minimum at $10^{11}\ M_{\odot}$). The latter
corresponds roughly to $M_{\rm UV}\lesssim-19$, so it better reflects the
fraction that is directly observable.

Figure 14: Effect of dust obscuration on the global SFRD as a function of
redshift. In the upper panel, the orange curves show the SFRD including all
halos able to form stars (roughly $10^{8}-10^{13}\ M_{\odot}$ in our model),
while blue curves only include massive galaxies ($10^{11}-10^{13}M_{\odot}$,
which corresponds to $M_{\rm UV}\lesssim-19$ at $z\approx 5$ so are more
easily observable). The solid curves are the intrinsic SFRD, the dash-dotted
curves are the dust-obscured SFRD from [16], and the dotted curves are the
obscured SFRD computed from [21]’s IRX-$\beta$ relation for all halo masses.
In the lower panel, we show $f_{\rm obsc}$ directly. The purple symbols are
[28]’s stacking analysis of the obscured fraction in observed galaxies, with
the downward triangle representing their higher stellar mass stack and the
upward symbol representing their lower stellar mass stack. The green squares
are the ALPINE stacking analysis from [23].

Figure 14 shows the results, with the SFRD in the top panel and the obscured
fraction in the bottom panel. We see that at the low end of the redshift range
a substantial fraction of the star formation is hidden by dust, especially in
the massive galaxies. But our model predicts that by $z\sim 10$ obscuration is
a modest effect overall. In general, obscuration is minimal in faint galaxies
across this redshift interval, so if they dominate the photon budget for,
e.g., reionization [80], dust plays only a relatively minor role. We note that
this remains true even using the IRX-$\beta$ relation from [21], which has a
gentler dependence of the dust mass on luminosity. Even though the faint
galaxies have more dust according to their prescription, it is not enough to
increase the overall SFRD significantly.

We also compare our results to estimates based on observed sources from [23]
and [28]. They are quite consistent with our model, at least for massive
galaxies, but we emphasize that the (dominant) low-mass galaxies can differ
substantially, especially at high redshifts. We also note that all of these
results (including our own) are based on UV-selected sources. [81] has
suggested that the presence of UV-faint but massive galaxies, undetected by
HST, may contribute significantly to the obscured cosmic SFRD at $z\gtrsim 4$,
resulting in an underestimation of that quantity. Spectroscopic surveys across
diverse galaxy populations will ultimately be needed to better constrain the
obscured contribution.

##  6 Conclusions

We have presented a semi-analytic model of dust in high-$z$ galaxies, aiming
to constrain its parameters by simultaneously matching the observed UV and IR
properties of galaxies during reionization. Our dust model pairs a simple
galaxy evolution model [49] with a streamlined treatment of the dust
distribution [46, 61] to generate self-consistent predictions for the UV
attenuation and far-IR continuum emission. The key input parameters are the
dust mass and radial distribution. We constrain these quantities with the UV
attenuation inferred from the IRX-$\beta$ relation [16] and from ALMA
measurements of the far-IR continuum flux. We find:

_(i)_ The dust attenuation increases with galaxy mass and decreases with
redshift, as expected [26, 22]. In our model, very few galaxies have
significant attenuation beyond $z\sim 10$. The observed galaxies at $z\lesssim
10$, which are of course relatively bright, suffer from non-negligible
attenuation. Nevertheless, our model predicts that these modest extinctions
are far less than one expects from normal dust production models, if the dust
and stars are confined to rotation-supported disks.

_(ii)_ By simultaneously calibrating to the UV attenuation and IR emission,
our model requires that the total dust mass be close to the amount produced by
normal stellar populations but that the dust be significantly more extended
than a rotationally-supported disk. This conclusion holds regardless of
whether we assume the dust to be distributed in a shell surrounding the galaxy
or mixed throughout the stellar population. Interestingly, we find that dust
may be better retained (albeit with a larger relative extent) in the smallest
galaxies, as previously suggested by [25, 74, 29] but in contrast to
expectations from standard feedback models.

_(iii)_ However, the inferred radii of the dust distribution are larger than
those measured by [27] by a factor of a few. One way to reconcile these
results – while maintaining the relatively large dust luminosities – is if the
dust has a more complicated distribution, with the dust and young stars
segregated from each other. In that case, if the dust is distributed on the
spatial scale suggested by [27], we find typical dust covering fractions
$f_{\rm cov}\sim 0.2$–$0.7$ for massive galaxies (similar to many ALMA-derived
models; e.g. [23, 25, 26, 28]), while faint galaxies have $f_{\rm cov}\approx
0$. This could result from feedback ejecting the dust locally from star-
forming regions but not from the halo.

Our results have implications not only for dust evolution but potentially for
other aspects of early galaxy formation. For example, reionization requires a
substantial fraction of ionizing photons to escape from early galaxies. If
indeed feedback has cleared out channels in the dust distribution, it may also
have cleared out neutral gas – allowing these ionizing photons to escape more
easily. Additionally, the large dust luminosities suggest that most of the
dust is retained in the galaxy’s halo, which means feedback may not be
effective at expelling metals to large distances (c.f. [82]).

However, there are several caveats to our model. Most importantly, our galaxy
formation model associates halos to star formation histories (and dust
properties) on a one-to-one basis, such that the SFR is a monotonically
increasing function of halo mass. Real halos have a substantial scatter in
their SFRs, which makes a comparison with observations challenging. Using a
simple treatment of burstiness, we found that it can explain the scatter in
the observations, but it is not clear if it is responsible for the apparent
decrease in the retention fraction of dust as a function of SFR. Additionally,
we have assumed that the IRX–$\beta$ relation from low-$z$ galaxies [16]
applies to high-$z$ sources. We have checked that our qualitative conclusions
hold when other relations, based on higher-redshift sources [19, 20, 21], are
used, but the quantitative results are sensitive to this relation. Finally,
while our results are suggestive of a complex dust morphology, we do not
attempt to model it or explore the physical mechanism behind it. Future
observations will no doubt shed light on all of these issues.

##  Appendix A The UV Slope–Luminosity Relation in Early JWST Surveys

Figure 15: The $\beta$-$M_{\mathrm{UV}}$ relation at redshifts $6.5\leq z\leq
12.5$ separated into six redshift bins and compared to measurements from the
literature (blue points; see text for a full list of sources). The purple
curve represents the best fit relation to the [67] and [6] measurements (as
used in the main text), while the light blue curves show the best fit to data
displayed here in each redshift interval.

In the main text, we based our UV slope fit on two large homogeneous samples
[67, 6], but other measurements are available in the literature. To check for
any obvious biases, here we compare our fits to other samples, which we have
assembled from ERO and ERS NIRCam imaging and ground-based near-infrared
imaging of the COSMOS field [72, 83], with additional pointings in the Hubble
ultra deep field (HUDF) field and GOODS-South field [84], the UNCOVER survey
[1], the JEMS survey [3], the NGDEEP field [2], and the CEERS survey [4, 7].
Figure 15 shows that our fit is largely consistent with these other
measurements, although the scatter is large, even at redshifts beyond the
$z\sim 10$ limit used in this work.

## Acknowledgments

This work was supported by by NASA through award 80NSSC22K0818 and by the
National Science Foundation through awards AST-1812458 and AST-2205900. We
thank Hiddo Algera, Yana Khusanova, and Michael Topping for providing their
data, as well as Sahil Hegde and Adam Trapp for useful discussions. JZ
acknowledges the fellowship support from the James and Meredith Henry
Endowment administered through the UCLA Undergraduate Research Scholars
Program (URSP) under the UCLA Undergraduate Research Center - Sciences. JZ
thanks Meadowridge School for hospitality while much of this work was
completed, along with Dale C. Stevenson for helpful conversations. This work
has made extensive use of NASA’s Astrophysics Data System
(http://ui.adsabs.harvard.edu/) and the arXiv e-Print service
(http://arxiv.org), as well as the following softwares: matplotlib [85], numpy
[86], Pandas [87], and scipy [88].

## References

  * [1] H. Atek, I. Chemerynska, B. Wang, L. Furtak, A. Weibel, P. Oesch et al., _JWST UNCOVER: Discovery of $z>9$ Galaxy Candidates Behind the Lensing Cluster Abell 2744_, _arXiv e-prints_ (2023) arXiv:2305.01793 [2305.01793].
  * [2] D. Austin, N.J. Adams, C.J. Conselice, T. Harvey, K. Ormerod, J. Trussler et al., _A Large Population of Faint 8 <z<16 Galaxies Found in the First JWST NIRCam Observations of the NGDEEP Survey_, _arXiv e-prints_ (2023) arXiv:2302.04270 [2302.04270].
  * [3] R.J. Bouwens, M. Stefanon, G. Brammer, P.A. Oesch, T. Herard-Demanche, G.D. Illingworth et al., _Evolution of the UV LF from z 15 to z 8 using new JWST NIRCam medium-band observations over the HUDF/XDF_ , _MNRAS,_ 523 (2023) 1036 [2211.02607].
  * [4] R. Endsley, D.P. Stark, L. Whitler, M.W. Topping, Z. Chen, A. Plat et al., _A JWST/NIRCam Study of key contributors to reionization: The star-forming and ionizing properties of UV-faint z 7 - 8 galaxies_ , _MNRAS,_ (2023)  [2208.14999].
  * [5] S. Tacchella, B.D. Johnson, B.E. Robertson, S. Carniani, F. D’Eugenio, N. Kumari et al., _JWST NIRCam + NIRSpec: interstellar medium and stellar populations of young galaxies with rising star formation and evolving gas reservoirs_ , _MNRAS,_ 522 (2023) 6236 [2208.03281].
  * [6] M.W. Topping, D.P. Stark, R. Endsley, L. Whitler, K. Hainline, B.D. Johnson et al., _The UV Continuum Slopes of Early Star-Forming Galaxies in JADES_ , _arXiv e-prints_ (2023) arXiv:2307.08835 [2307.08835].
  * [7] L. Whitler, R. Endsley, D.P. Stark, M. Topping, Z. Chen and S. Charlot, _On the ages of bright galaxies 500 Myr after the big bang: insights into star formation activity at z $\gtrsim$ 15 with JWST_, _MNRAS,_ 519 (2023) 157 [2208.01599].
  * [8] C.A. Mason, M. Trenti and T. Treu, _The brightest galaxies at cosmic dawn_ , _MNRAS,_ 521 (2023) 497 [2207.14808].
  * [9] J. Mirocha and S.R. Furlanetto, _Balancing the efficiency and stochasticity of star formation with dust extinction in z $\gtrsim$ 10 galaxies observed by JWST_, _MNRAS,_ 519 (2023) 843 [2208.12826].
  * [10] A. Ferrara, A. Pallottini and P. Dayal, _On the stunning abundance of super-early, luminous galaxies revealed by JWST_ , _MNRAS,_ 522 (2023) 3986 [2208.00720].
  * [11] F. Ziparo, A. Ferrara, L. Sommovigo and M. Kohandel, _Blue monsters. Why are JWST super-early, massive galaxies so blue?_ , _MNRAS,_ 520 (2023) 2445 [2209.06840].
  * [12] H. Inami, H.S.B. Algera, S. Schouws, L. Sommovigo, R. Bouwens, R. Smit et al., _The ALMA REBELS Survey: dust continuum detections at z > 6.5_, _MNRAS,_ 515 (2022) 3126 [2203.15136].
  * [13] D. Calzetti, _The Dust Opacity of Star-forming Galaxies_ , _PASP_ 113 (2001) 1449 [astro-ph/0109035].
  * [14] J.C. Weingartner and B.T. Draine, _Dust Grain-Size Distributions and Extinction in the Milky Way, Large Magellanic Cloud, and Small Magellanic Cloud_ , _ApJ,_ 548 (2001) 296 [astro-ph/0008146].
  * [15] S. Schouws, M. Stefanon, R. Bouwens, R. Smit, J. Hodge, I. Labbé et al., _Significant Dust-obscured Star Formation in Luminous Lyman-break Galaxies at z 7-8_ , _ApJ,_ 928 (2022) 31 [2105.12133].
  * [16] G.R. Meurer, T.M. Heckman and D. Calzetti, _Dust Absorption and the Ultraviolet Luminosity Density at z ~3 as Calibrated by Local Starburst Galaxies_, _ApJ,_ 521 (1999) 64 [astro-ph/9903054].
  * [17] S.L. Finkelstein, C. Papovich, B. Salmon, K. Finlator, M. Dickinson, H.C. Ferguson et al., _Candels: The Evolution of Galaxy Rest-frame Ultraviolet Colors from z = 8 to 4_ , _ApJ,_ 756 (2012) 164 [1110.3785].
  * [18] C.M. Casey, N.Z. Scoville, D.B. Sanders, N. Lee, A. Cooray, S.L. Finkelstein et al., _Are Dusty Galaxies Blue? Insights on UV Attenuation from Dust-selected Galaxies_ , _ApJ,_ 796 (2014) 95 [1410.0702].
  * [19] R.J. McLure, J.S. Dunlop, F. Cullen, N. Bourne, P.N. Best, S. Khochfar et al., _Dust attenuation in 2 < z < 3 star-forming galaxies from deep ALMA observations of the Hubble Ultra Deep Field_, _MNRAS,_ 476 (2018) 3991 [1709.06102].
  * [20] N.A. Reddy, P.A. Oesch, R.J. Bouwens, M. Montes, G.D. Illingworth, C.C. Steidel et al., _The HDUV Survey: A Revised Assessment of the Relationship between UV Slope and Dust Attenuation for High-redshift Galaxies_ , _ApJ,_ 853 (2018) 56 [1705.09302].
  * [21] Y. Fudamoto, P.A. Oesch, A. Faisst, M. Béthermin, M. Ginolfi, Y. Khusanova et al., _The ALPINE-ALMA [CII] survey. Dust attenuation properties and obscured star formation at z $\sim$ 4.4-5.8_, _A &A_ 643 (2020) A4 [2004.10760].
  * [22] R.A.A. Bowler, H. Inami, L. Sommovigo, R. Smit, H.S.B. Algera, M. Aravena et al., _The ALMA REBELS survey: obscured star formation in massive Lyman-break galaxies at z = 4-8 revealed by the IRX- $\beta$ and $M_{\star}$ relations_, _arXiv e-prints_ (2023) arXiv:2309.17386 [2309.17386].
  * [23] Y. Khusanova, M. Bethermin, O. Le Fèvre, P. Capak, A.L. Faisst, D. Schaerer et al., _The ALPINE-ALMA [CII] survey. Obscured star formation rate density and main sequence of star-forming galaxies at z > 4_, _A &A_ 649 (2021) A152 [2007.08384].
  * [24] L. Sommovigo, A. Ferrara, S. Carniani, A. Zanella, A. Pallottini, S. Gallerani et al., _Dust temperature in ALMA [C II]-detected high-z galaxies_ , _MNRAS,_ 503 (2021) 4878 [2102.08950].
  * [25] P. Dayal, A. Ferrara, L. Sommovigo, R. Bouwens, P.A. Oesch, R. Smit et al., _The ALMA REBELS survey: the dust content of z 7 Lyman break galaxies_ , _MNRAS,_ 512 (2022) 989 [2202.11118].
  * [26] A. Ferrara, L. Sommovigo, P. Dayal, A. Pallottini, R.J. Bouwens, V. Gonzalez et al., _The ALMA REBELS Survey. Epoch of Reionization giants: Properties of dusty galaxies at z $\approx$ 7_, _MNRAS,_ 512 (2022) 58 [2202.07666].
  * [27] Y. Fudamoto, R. Smit, R.A.A. Bowler, P.A. Oesch, R. Bouwens, M. Stefanon et al., _The ALMA REBELS Survey: Average [C II] 158 $\mu$m Sizes of Star-forming Galaxies from z 7 to z 4_, _ApJ,_ 934 (2022) 144 [2206.01886].
  * [28] H.S.B. Algera, H. Inami, P.A. Oesch, L. Sommovigo, R.J. Bouwens, M.W. Topping et al., _The ALMA REBELS survey: the dust-obscured cosmic star formation rate density at redshift 7_ , _MNRAS,_ 518 (2023) 6142 [2208.08243].
  * [29] M. Palla, I. De Looze, M. Relaño, S. van der Giessen, P. Dayal, A. Ferrara et al., _Metal and dust evolution in ALMA REBELS galaxies: insights for future JWST observations_ , _arXiv e-prints_ (2023) arXiv:2311.16071 [2311.16071].
  * [30] M. Béthermin, Y. Fudamoto, M. Ginolfi, F. Loiacono, Y. Khusanova, P.L. Capak et al., _The ALPINE-ALMA [CII] survey: Data processing, catalogs, and statistical source properties_ , _A &A_ 643 (2020) A2 [2002.00962].
  * [31] A.L. Faisst, D. Schaerer, B.C. Lemaux, P.A. Oesch, Y. Fudamoto, P. Cassata et al., _The ALPINE-ALMA [C II] Survey: Multiwavelength Ancillary Data and Basic Physical Measurements_ , _ApJS,_ 247 (2020) 61 [1912.01621].
  * [32] I. Mitsuhashi, K.-i. Tadaki, R. Ikeda, R. Herrera-Camus, M. Aravena, I. De Looze et al., _The ALMA-CRISTAL survey: Widespread dust-obscured star formation in typical star-forming galaxies at z=4-6_ , _arXiv e-prints_ (2023) arXiv:2311.17671 [2311.17671].
  * [33] I. Mitsuhashi, Y. Harikane, F.E. Bauer, T. Bakx, A. Ferrara, S. Fujimoto et al., _SERENADE II: An ALMA Multi-Band Dust-Continuum Analysis of 28 Galaxies at $5<z<8$ and the Physical Origin of the Dust Temperature Evolution_, _arXiv e-prints_ (2023) arXiv:2311.16857 [2311.16857].
  * [34] R.J. Bouwens, R. Smit, S. Schouws, M. Stefanon, R. Bowler, R. Endsley et al., _Reionization Era Bright Emission Line Survey: Selection and Characterization of Luminous Interstellar Medium Reservoirs in the z > 6.5 Universe_, _ApJ,_ 931 (2022) 160 [2106.13719].
  * [35] T. Nozawa, T. Kozasa, H. Umeda, K. Maeda and K. Nomoto, _Dust in the Early Universe: Dust Formation in the Ejecta of Population III Supernovae_ , _ApJ,_ 598 (2003) 785 [astro-ph/0307108].
  * [36] E.S. Brooker, S.M. Stangl, C.M. Mauney and C.L. Fryer, _Dependence of Dust Formation on the Supernova Explosion_ , _ApJ,_ 931 (2022) 85 [2103.12781].
  * [37] R. Valiante, R. Schneider, S. Bianchi and A.C. Andersen, _Stellar sources of dust in the high-redshift Universe_ , _MNRAS,_ 397 (2009) 1661 [0905.1691].
  * [38] M. Mancini, R. Schneider, L. Graziani, R. Valiante, P. Dayal, U. Maio et al., _The dust mass in z > 6 normal star-forming galaxies._, _MNRAS,_ 451 (2015) L70 [1505.01841].
  * [39] D. Anglés-Alcázar, C.-A. Faucher-Giguère, D. Kereš, P.F. Hopkins, E. Quataert and N. Murray, _The cosmic baryon cycle and galaxy mass assembly in the FIRE simulations_ , _MNRAS,_ 470 (2017) 4698 [1610.08523].
  * [40] B.T. Draine, _Interstellar Dust Grains_ , _ARA &A_ 41 (2003) 241 [astro-ph/0304489].
  * [41] D. Narayanan, R. Davé, B.D. Johnson, R. Thompson, C. Conroy and J. Geach, _The IRX- $\beta$ dust attenuation relation in cosmological galaxy formation simulations_, _MNRAS,_ 474 (2018) 1718 [1705.05858].
  * [42] M. Vogelsberger, D. Nelson, A. Pillepich, X. Shen, F. Marinacci, V. Springel et al., _High-redshift JWST predictions from IllustrisTNG: dust modelling and galaxy luminosity functions_ , _MNRAS,_ 492 (2020) 5167 [1904.07238].
  * [43] M. Mushtaq, D. Ceverino, R.S. Klessen, S. Reissl and P.H. Puttasiddappa, _Dust attenuation in galaxies at cosmic dawn from the FirstLight simulations_ , _arXiv e-prints_ (2023) arXiv:2304.10150 [2304.10150].
  * [44] D. Narayanan, J.D.T. Smith, B.S. Hensley, Q. Li, C.-Y. Hu, K. Sandstrom et al., _A Framework for Modeling Polycyclic Aromatic Hydrocarbon Emission in Galaxy Evolution Simulations_ , _ApJ,_ 951 (2023) 100 [2301.07136].
  * [45] E. da Cunha, B. Groves, F. Walter, R. Decarli, A. Weiss, F. Bertoldi et al., _On the Effect of the Cosmic Microwave Background in High-redshift (Sub-)millimeter Observations_ , _ApJ,_ 766 (2013) 13 [1302.0844].
  * [46] N. Imara, A. Loeb, B.D. Johnson, C. Conroy and P. Behroozi, _A Model Connecting Galaxy Masses, Star Formation Rates, and Dust Temperatures across Cosmic Time_ , _ApJ,_ 854 (2018) 36 [1801.01499].
  * [47] A.K. Inoue, T. Hashimoto, H. Chihara and C. Koike, _Radiative equilibrium estimates of dust temperature and mass in high-redshift galaxies_ , _MNRAS,_ 495 (2020) 1577 [2004.12612].
  * [48] V. Mauerhofer and P. Dayal, _The dust enrichment and detectability of early galaxies in the JWST and ALMA era_ , _arXiv e-prints_ (2023) arXiv:2305.01681 [2305.01681].
  * [49] S.R. Furlanetto, J. Mirocha, R.H. Mebane and G. Sun, _A minimalist feedback-regulated model for galaxy formation during the epoch of reionization_ , _MNRAS,_ 472 (2017) 1576 [1611.01169].
  * [50] Planck Collaboration, N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi et al., _Planck 2018 results. VI. Cosmological parameters_ , _A &A_ 641 (2020) A6 [1807.06209].
  * [51] J.B. Oke and J.E. Gunn, _Secondary standard stars for absolute spectrophotometry._ , _ApJ,_ 266 (1983) 713.
  * [52] H. Trac, R. Cen and P. Mansfield, _SCORCH I: The Galaxy-Halo Connection in the First Billion Years_ , _ApJ,_ 813 (2015) 54 [1507.02685].
  * [53] A. Vale and J.P. Ostriker, _Linking halo mass to galaxy luminosity_ , _MNRAS,_ 353 (2004) 189 [astro-ph/0402500].
  * [54] C.-A. Faucher-Giguère, D. Kereš and C.-P. Ma, _The baryonic assembly of dark matter haloes_ , _MNRAS,_ 417 (2011) 2982 [1103.0001].
  * [55] A.C. Trapp and S.R. Furlanetto, _A flexible analytic model of cosmic variance in the first billion years_ , _MNRAS,_ 499 (2020) 2401 [2009.05059].
  * [56] P.F. Hopkins, D. Kereš, J. Oñorbe, C.-A. Faucher-Giguère, E. Quataert, N. Murray et al., _Galaxies on FIRE (Feedback In Realistic Environments): stellar feedback explains cosmologically inefficient star formation_ , _MNRAS,_ 445 (2014) 581 [1311.2073].
  * [57] P. Madau and M. Dickinson, _Cosmic Star-Formation History_ , _ARA &A_ 52 (2014) 415 [1403.0007].
  * [58] I. Cherchneff and E. Dwek, _The Chemistry of Population III Supernova Ejecta. II. The Nucleation of Molecular Clusters as a Diagnostic for Dust in the Early Universe_ , _ApJ,_ 713 (2010) 1 [1002.3060].
  * [59] S.R. Furlanetto and J. Mirocha, _Bursty star formation during the Cosmic Dawn driven by delayed stellar feedback_ , _MNRAS,_ 511 (2022) 3895 [2109.04488].
  * [60] M. Matsuura, P.M. Woods and P.J. Owen, _The global gas and dust budget of the Small Magellanic Cloud_ , _MNRAS,_ 429 (2013) 2527 [1212.1468].
  * [61] J. Mirocha, C. Mason and D.P. Stark, _Effects of self-consistent rest-ultraviolet colours in semi-empirical galaxy formation models_ , _MNRAS,_ 498 (2020) 2645 [2005.07208].
  * [62] R.J. Bouwens, M. Aravena, R. Decarli, F. Walter, E. da Cunha, I. Labbé et al., _ALMA Spectroscopic Survey in the Hubble Ultra Deep Field: The Infrared Excess of UV-Selected z = 2-10 Galaxies as a Function of UV-Continuum Slope and Stellar Mass_ , _ApJ,_ 833 (2016) 72 [1606.05280].
  * [63] A. Loeb and S.R. Furlanetto, _The First Galaxies in the Universe_ (2013).
  * [64] C.-A. Faucher-Giguère, _A model for the origin of bursty star formation in galaxies_ , _MNRAS,_ 473 (2018) 3717 [1701.04824].
  * [65] F. Városi and E. Dwek, _Analytical Approximations for Calculating the Escape and Absorption of Radiation in Clumpy Dusty Environments_ , _ApJ,_ 523 (1999) 265 [astro-ph/9905029].
  * [66] S. Bianchi, V. Casasola, M. Baes, C.J.R. Clark, E. Corbelli, J.I. Davies et al., _Dust emissivity and absorption cross section in DustPedia late-type galaxies_ , _A &A_ 631 (2019) A102 [1909.12692].
  * [67] R.J. Bouwens, G.D. Illingworth, P.A. Oesch, I. Labbé, P.G. van Dokkum, M. Trenti et al., _UV-continuum Slopes of >4000 z ~4-8 Galaxies from the HUDF/XDF, HUDF09, ERS, CANDELS-South, and CANDELS-North Fields_, _ApJ,_ 793 (2014) 115 [1306.2950].
  * [68] F. Cullen, D.J. McLeod, R.J. McLure, J.S. Dunlop, C.T. Donnan, A.C. Carnall et al., _Evidence for the emergence of dust-free stellar populations at z > 10_, _arXiv e-prints_ (2023) arXiv:2311.06209 [2311.06209].
  * [69] R.P. Naidu, P.A. Oesch, D.J. Setton, J. Matthee, C. Conroy, B.D. Johnson et al., _Schrodinger’s Galaxy Candidate: Puzzlingly Luminous at $z\approx 17$, or Dusty/Quenched at $z\approx 5$?_, _arXiv e-prints_ (2022) arXiv:2208.02794 [2208.02794].
  * [70] S.R. Furlanetto and J. Mirocha, _On the expected purity of photometric galaxy surveys targeting the Cosmic Dawn_ , _MNRAS,_ 523 (2023) 5274 [2208.12828].
  * [71] J.A. Zavala, V. Buat, C.M. Casey, S.L. Finkelstein, D. Burgarella, M.B. Bagley et al., _Dusty Starbursts Masquerading as Ultra-high Redshift Galaxies in JWST CEERS Observations_ , _ApJ,l_ 943 (2023) L9 [2208.01816].
  * [72] F. Cullen, R.J. McLure, D.J. McLeod, J.S. Dunlop, C.T. Donnan, A.C. Carnall et al., _The ultraviolet continuum slopes ( $\beta$) of galaxies at z $\simeq$ 8-16 from JWST and ground-based near-infrared imaging_, _MNRAS,_ 520 (2023) 14 [2208.04914].
  * [73] J.B. Muñoz, J. Mirocha, S. Furlanetto and N. Sabti, _Breaking degeneracies in the first galaxies with clustering_ , _MNRAS,_ 526 (2023) L47 [2306.09403].
  * [74] M.W. Topping, D.P. Stark, R. Endsley, R.J. Bouwens, S. Schouws, R. Smit et al., _The ALMA REBELS Survey: specific star formation rates in the reionization era_ , _MNRAS_ 516 (2022) 975 [2203.07392].
  * [75] L. Sommovigo, A. Ferrara, A. Pallottini, P. Dayal, R.J. Bouwens, R. Smit et al., _The ALMA REBELS Survey: cosmic dust temperature evolution out to z 7_ , _MNRAS,_ 513 (2022) 3122 [2202.01227].
  * [76] M.E. De Rossi and V. Bromm, _Probing the Nature of the First Galaxies with JWST and ALMA_ , _ApJL_ 946 (2023) L20 [2303.06328].
  * [77] A. Dekel, K.C. Sarkar, Y. Birnboim, N. Mandelker and Z. Li, _Efficient formation of massive galaxies at cosmic dawn by feedback-free starbursts_ , _MNRAS,_ 523 (2023) 3201 [2303.04827].
  * [78] D. Tsuna, Y. Nakazato and T. Hartwig, _A Photon Burst Clears the Earliest Dusty Galaxies: Modelling Dust in High-redshift Galaxies from ALMA to JWST_ , _arXiv e-prints_ (2023) arXiv:2309.02415 [2309.02415].
  * [79] R. Barkana and A. Loeb, _In the beginning: the first sources of light and the reionization of the universe_ , _Physics Reports_ 349 (2001) 125 [astro-ph/0010468].
  * [80] C. Simmonds, S. Tacchella, K. Hainline, B.D. Johnson, W. McClymont, B. Robertson et al., _Low-mass bursty galaxies in JADES efficiently produce ionising photons and could represent the main drivers of reionisation_ , _arXiv e-prints_ (2023) arXiv:2310.01112 [2310.01112].
  * [81] C.C. Williams, S. Alberts, Z. Ji, K.N. Hainline, J. Lyu, G. Rieke et al., _The galaxies missed by Hubble and ALMA: the contribution of extremely red galaxies to the cosmic census at 3 <z<8_, _arXiv e-prints_ (2023) arXiv:2311.07483 [2311.07483].
  * [82] N. Yamaguchi, S.R. Furlanetto and A.C. Trapp, _The extent of intergalactic metal enrichment from galactic winds during the Cosmic Dawn_ , _MNRAS,_ 520 (2023) 2922 [2209.09345].
  * [83] C.T. Donnan, D.J. McLeod, J.S. Dunlop, R.J. McLure, A.C. Carnall, R. Begley et al., _The evolution of the galaxy UV luminosity function at redshifts z $\simeq$ 8 - 15 from deep JWST and ground-based near-infrared imaging_, _MNRAS,_ 518 (2023) 6011 [2207.12356].
  * [84] C.T. Donnan, D.J. McLeod, R.J. McLure, J.S. Dunlop, A.C. Carnall, F. Cullen et al., _The abundance of z $\gtrsim$ 10 galaxy candidates in the HUDF using deep JWST NIRCam medium-band imaging_, _MNRAS,_ 520 (2023) 4554 [2212.10126].
  * [85] J.D. Hunter, _Matplotlib: A 2d graphics environment_ , _Computing in Science & Engineering_ 9 (2007) 90.
  * [86] C.R. Harris, K.J. Millman, S.J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau et al., _Array programming with NumPy_ , _Nature_ 585 (2020) 357.
  * [87] T. pandas development team, _pandas-dev/pandas: Pandas_ , Feb., 2020.  10.5281/zenodo.3509134.
  * [88] P. Virtanen, R. Gommers, T.E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau et al., _SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python_ , _Nature Methods_ 17 (2020) 261.