# Disappearing galaxies: the orientation dependence of JWST-bright, HST-dark,
star-forming galaxy selection

R. K. Cochrane Department of Astronomy, Columbia University, New York, NY
10027, USA  Center for Computational Astrophysics, Flatiron Institute, 162
Fifth Avenue, New York, NY 10010, USA  D. Anglés-Alcázar Department of
Physics, University of Connecticut, 196 Auditorium Road, U-3046, Storrs, CT
06269-3046, USA
Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue,
New York, NY 10010, USA  F. Cullen Institute for Astronomy, University of
Edinburgh, Royal Observatory, Blackford Hill, Edinburgh, EH9 3HJ, UK  C. C.
Hayward Center for Computational Astrophysics, Flatiron Institute, 162 Fifth
Avenue, New York, NY 10010, USA

###### Abstract

Galaxies that are invisible in deep optical-NIR imaging but detected at longer
wavelengths have been the focus of several recent observational studies, with
speculation that they could constitute a substantial missing population and
even dominate the cosmic star formation rate density at $z\gtrsim 4$. The
depths now achievable with JWST at the longest wavelengths probed by HST,
coupled with the transformative resolution at longer wavelengths, are already
enabling detailed, spatially-resolved characterisation of sources that were
invisible to HST, often known as ‘HST-dark’ galaxies. However, until now,
there has been little theoretical work to compare against. We present the
first simulation-based study of this population, using highly-resolved
galaxies from the Feedback in Realistic Environments (FIRE) project, with
multi-wavelength images along several lines of sight forward-modelled using
radiative transfer. We naturally recover a population of modelled sources that
meet commonly-used selection criteria ($H_{\rm{AB}}>27\,\rm{mag}$ and
$H_{\rm{AB}}-\rm{F444W}>2.3$). These simulated HST-dark galaxies lie at high
redshifts ($z=4-7$), have high levels of dust attenuation ($A_{V}=2-4$), and
display compact recent star formation ($R_{1/2,\,\rm{4.4\,\mu\rm{m}}}\lesssim
1\,\rm{kpc}$). Orientation is very important: for all but one of the 17
simulated galaxy snapshots with HST-dark sightlines, there exist other
sightlines that do not meet the criteria. This result has important
implications for comparisons between observations and models that do not
resolve the detailed star-dust geometry, such as semi-analytic models or
coarsely-resolved hydrodynamical simulations. Critically, we demonstrate that
HST-dark sources are not an unexpected or exotic population, but a subset of
high-redshift, highly-dust-attenuated sources viewed along certain lines of
sight.

galaxies: evolution – galaxies: high-redshift – galaxies: ISM – radiative
transfer

figuret

111email<EMAIL_ADDRESS>
##  1 Introduction

It has been known for decades that not all bright sub-millimeter galaxies
(SMGs) have detectable optical/near-infrared (OIR) counterparts (Lilly et al.,
1999; Smail et al., 1999, 2000, 2002a, 2004; Bertoldi et al., 2000; Frayer et
al., 2000, 2004; Dannerbauer et al., 2002; Wang et al., 2007). This can hamper
redshift determination; indeed, the brightest source in the Hubble Deep Field
at $850\,\mu\rm{m}$, HDF850.1, showed no clear OIR counterpart (Hughes et al.,
1998; Cowie et al., 2009), and it took over a decade for its redshift to be
identified ($z=5.183$; Walter et al., 2012).
Other populations of OIR-dark galaxies have been selected using optical-
infrared colors. Galaxies first identified by red $R-K$ colors and faint
$K-$band fluxes were called ‘Extremely Red Objects’ (EROs; Elston et al.,
1988; Hu & Ridgway, 1994); some, but not all, of these sources were also
identified in the sub-millimeter (Cimatti et al., 1998; Dey et al., 1999).
Later, Spitzer-IRAC imaging (Fazio et al., 2004) enabled efficient surveys of
new red populations (Wilson et al., 2004; Huang et al., 2011) and better
separation of passive and star-forming galaxies. Wang et al. (2016) introduced
a color selection to identify massive, high-redshift galaxies using a
combination of deep HST $H-$band imaging and IRAC CH2:
$H_{\rm{AB}}-[4.5]>2.25\,\rm{mag}$. $62\%$ of sources selected in this way
were detected down to $S_{870\,\mu\rm{m}}>0.6\,\rm{mJy}$ (Wang et al., 2019).
OIR-faint sources have also been detected at radio wavelengths (Kondapally et
al., 2021; Talia et al., 2021; Enia et al., 2022; van der Vlugt et al., 2022;
Behiri et al., 2023).
OIR-faint galaxies are important for several reasons. Firstly, several works
have suggested that they may contribute significantly to the cosmic star
formation rate density (SFRD) at $z\gtrsim 3$, though estimates of the
contribution vary (e.g. Wang et al. 2016 calculated that $H-$dark sources
contribute $15-25\,\%$ of the SFRD at $z=4-5$, while Sun et al. 2021 estimated
$8^{+8}_{-4}\,\%$ at $z=3-5$; see also Williams et al. 2019; Yamaguchi et al.
2019; Xiao et al. 2023). Hence, excluding OIR-faint sources from censuses of
star formation (e.g. by using $H-$ or $K-$selected samples, or by requiring
OIR counterparts for sources selected in other ways) will lead to a redshift-
dependent underestimate of the SFRD. Secondly, degeneracies between the
signatures of redshift, age and dust attenuation in photometry result in the
possibility of dusty sources contaminating surveys of high-redshift quiescent
galaxies (e.g. Smail et al., 2002b; Dunlop et al., 2007; Simpson et al.,
2017a). Thirdly, OIR-faint sources may help constrain theoretical models of
galaxy formation and evolution. Wang et al. (2019) showed that the number
density of $H_{\rm{AB}}>27\,\rm{mag}$ sources predicted by an early version of
the l-galaxies semi-analytic model (SAM) (Henriques et al., 2015) lies $\sim
2$ orders of magnitude below the observed number density. They speculated that
our understanding of massive-galaxy formation may require substantial
revision.
The advent of ALMA has facilitated statistical studies of large samples of
sub-millimeter bright galaxies, including locating and resolving sources
identified with single dish telescopes (e.g. Hodge et al., 2013; Karim et al.,
2013; Simpson et al., 2014, 2015, 2020; Stach et al., 2018). This has enabled
more detailed exploration of the number density and physical properties of the
subset of sub-millimeter-bright galaxies that are also OIR-dark. Using ALMA
follow-up of SCUBA-2-identified sources within the UKIDSS Ultra Deep Survey
(Simpson et al., 2017b; Stach et al., 2018, 2019; Dudzevičiūtė et al., 2020),
Smail et al. (2021) found that $15\pm 2$ % of the bright
($S_{850\,\mu\rm{m}}>3.6\,\rm{mJy}$) SMG population is fainter than
$K_{\rm{AB}}=25.3\,\rm{mag}$. These $K-$faint galaxies tend to have higher
redshifts and dust attenuation than their $K-$detected counterparts (see also
the high photometric redshifts derived for the $HST-$dark,
$1.1\,\rm{mm}$-selected sources from Franco et al. 2018, and for OIR-dark
galaxies selected at longer wavelengths: Williams et al. 2019; Casey et al.
2021; Manning et al. 2022). However, given the relatively low angular
resolution of IRAC ($\rm{FWHM}\sim 2^{\prime\prime}$) and the typically small
angular sizes of OIR-dark sources, characterising the spatial distribution of
rest-frame optical emission was not possible until the launch of JWST.
JWST has provided the sensitivity and wavelength coverage required to gather
high angular resolution (compared to IRAC), multi-wavelength imaging and hence
constrain the properties of sources that were invisible to HST (e.g. Fujimoto
et al., 2023; Gómez-Guijarro et al., 2023; Pérez-González et al., 2023;
Rodighiero et al., 2023; Smail et al., 2023). Employing a similar selection
used in previous studies to identify dusty galaxies at $z=3-6$
($H_{\rm{AB}}>27$ & $H_{\rm{AB}}-\rm{F444W_{AB}}>2.3$; Caputi et al. 2012;
Wang et al. 2016), Barrufet et al. (2022) selected a sample of 33 HST-dark
sources within CEERS NIRCam imaging (Finkelstein et al., 2022). They showed
that galaxies selected using these criteria tend to be fairly dust-obscured
($A_{V}\sim 2$), massive ($M_{\star}\sim 10^{9}-10^{10.5}\,\rm{M_{\odot}}$)
star-forming galaxies at $z\sim 2-8$ that lie approximately along the main
sequence. Using the same data, Nelson et al. (2022) identified a population of
spatially extended ($R_{e,\,4.4\,\mu\rm{m}}>0.17^{\prime\prime}$),
$4.4\,\mu\rm{m}$-bright, $H-$faint sources, with higher inferred stellar
masses, the majority in the range $M_{\star}=10^{10-11}\,\rm{M_{\odot}}$
(using EAzY, the median stellar mass is $M_{\star}=10^{10.5}\,\rm{M_{\odot}}$,
but they also find systematic differences between SED fitting codes). These
$4.4\,\mu\rm{m}$-identified sources appear to have lower stellar masses, on
average, than sub-millimeter-identified OIR-dark sources. For example, the
median stellar masses of the $K$-faint and whole sample of Smail et al. 2021
are $M_{\star}=10^{11.10\pm 0.04}\,\rm{M_{\odot}}$ and $M_{\star}=10^{11.00\pm
0.06}\,\rm{M_{\odot}}$, respectively, and some of the most extreme sub-
millimeter bright sources are likely even more massive (e.g. the stellar mass
derived for the $S_{850\,\mu\rm{m}}=15.3\pm 0.4\,\rm{mJy}$, $z=4.26$ source
presented by Smail et al. 2023 is $10^{11.8}\,\rm{M_{\odot}}$).
Despite substantial observational interest in HST-dark galaxies, until now
there have been no focused theoretical studies on the nature of HST-dark
galaxies. In this paper, we study the physical properties of galaxies modelled
with highly-resolved zoom-in simulations using the Feedback in Realistic
Environments (FIRE; Hopkins et al. 2014, 2018, 2023) model that meet
observational HST-dark selection criteria. In Section 2, we describe the
simulations and the radiative transfer calculations performed to forward-model
observable emission. We describe the selection of HST-dark galaxies. In
Section 3, we study the demographics, dust attenuation properties and sizes of
our simulated HST-dark galaxies, and highlight a substantial line-of-sight
dependence of the source section. We draw our conclusions in Section 4.

##  2 Simulations and HST-dark source selection

###  2.1 FIRE simulations

The FIRE222https://fire.northwestern.edu project is a suite of state-of-the-
art hydrodynamical cosmological zoom-in simulations described fully in Hopkins
et al. (2014, FIRE-1), Hopkins et al. (2018, FIRE-2) and Hopkins et al. (2023,
FIRE-3). In this paper, we study galaxies modelled with FIRE-2, which uses the
“meshless finite mass” mode of the $N$-body+hydrodynamics code
GIZMO333http://www.tapir.caltech.edu/~phopkins/Site/GIZMO.html (Hopkins,
2015); gravitational forces are computed following the methods presented in
Hopkins et al. (2013), using an improved version of the parallel TreeSPH code
GADGET-3 (Springel et al., 2005). Cooling and heating processes including
free-free, photoionization/recombination, Compton, photo-electric, metal-line,
molecular and fine-structure processes are modelled from $T=10\,\rm{K}$ to
$T=10^{10}\,\rm{K}$. Star particles form from locally self-gravitating,
molecular, Jeans unstable gas above a minimum hydrogen number density
$n_{H}\geq 1000\,\rm{cm}^{-3}$. Each star particle represents a single stellar
population with known mass, age, and metallicity, injecting feedback locally
in the form of mass, momentum, energy, and metals from Type Ia and Type II
Supernovae (SNe), stellar winds, photoionization and photoelectric heating,
and radiation pressure, with all feedback quantities and their time dependence
taken directly from the starburst99 population synthesis model (Leitherer et
al., 1999).
We study the central galaxies of eight massive haloes originally selected and
simulated by Feldmann et al. (2016, 2017) as part of the MassiveFIRE-1 suite.
The same haloes were studied in Cochrane et al. (2022) and Cochrane et al.
(2023). The first four haloes are drawn from the ‘A-series’ (A1, A2, A4 and
A8); these haloes were selected to have dark matter halo masses of
$M_{\rm{halo}}\sim 10^{12.5}\,\rm{M_{\odot}}$ at $z=2$. The ‘A-series’ halos
studied in this paper are drawn from Anglés-Alcázar et al. (2017), who re-
simulated them down to $z=1$ with the upgraded FIRE-2 physics model (Hopkins
et al., 2018). We supplement these haloes by re-running four more haloes from
Feldmann et al. (2017), with the updated FIRE-2 physics. Two of the haloes are
drawn from their ‘B-series’ (B1 and B2) and two from the ‘C-series’ (Cm1:0,
hereafter C1, and C2:0, hereafter C2). Haloes B2, C1 and C2 tend to be more
massive than the A-series haloes. The mass resolution for both gas and star
particles is $3.3\times 10^{4}\,\rm{M_{\odot}}$, $2.7\times
10^{5}\,\rm{M_{\odot}}$, and $2.2\times 10^{6}\,\rm{M_{\odot}}$, for A-, B-,
and C-series haloes, respectively. For dark matter particles, the respective
mass resolutions are: $1.7\times 10^{5}\,\rm{M_{\odot}}$, $1.4\times
10^{6}\,\rm{M_{\odot}}$, and $1.1\times 10^{7}\,\rm{M_{\odot}}$. Convergence
tests for these simulations are presented in Cochrane et al. (2022) (see
Appendix B).

###  2.2 Modelling observable emission

We model spectral energy distributions (SEDs) and multi-wavelength emission
maps at every snapshot, between $z\sim 8$ and $z\sim 1$, for each of the eight
simulated halos (in total, $1736$ snapshots). Following Cochrane et al. (2019,
2022, 2023); Cochrane & Angles-Alcazar (2023) we use the
skirt444http://www.skirt.ugent.be radiative transfer code (version 8; Baes et
al., 2011; Camps & Baes, 2015) to make predictions for emission between rest-
frame ultraviolet (UV) and far-infrared (FIR) wavelengths along seven lines of
sight. For all snapshots, gas and star particles within $0.1R_{\rm{vir}}$ are
drawn directly from FIRE-2 simulation data. Dust particles are assumed to
follow the distribution of the gas particles, with a dust-to-metals mass ratio
of $0.4$ (Dwek, 1998; James et al., 2002). This is a reasonable assumption for
enriched, massive galaxies like these. We have checked that using the
metallicity-dependent dust-to-metal ratios draw from the relations derived by
Rémy-Ruyer et al. (2014) would yield similar values. We assume dust
destruction at $>10^{6}\,{\rm K}$ (Draine & Salpeter, 1979; Tielens et al.,
1994). Following Cochrane et al. (2019, 2022, 2023); Cochrane & Angles-Alcazar
(2023), we model a mixture of graphite, silicate and PAH grains using the
Weingartner & Draine (2001) Milky Way dust prescription. Star particles are
assigned Bruzual & Charlot (2003) SEDs based on their ages and metallicities.
We perform the radiative transfer on an octree dust grid, in which cell sizes
are adjusted according to the dust density distribution, with the condition
that no dust cell may contain more than $0.0001\%$ of the total dust mass of
the galaxy. skirt parameter convergence tests are presented in Cochrane et al.
(2022).
The output from skirt comprises predictions for global galaxy SEDs as well as
maps of the resolved emission at each of the $\sim 100$ wavelengths modelled,
for seven lines of sight. We transform flux densities in Jy to AB magnitudes
(using a zero-point of $3631\,\rm{Jy}$).

###  2.3 HST-dark source selection

Following Barrufet et al. (2022), we select sightlines from our snapshots that
meet the following $H-$band and color criteria:

$H_{\rm{AB}}>27\,\rm{mag}$ (1)

and

$H_{\rm{AB}}-\rm{F444W}_{\rm{AB}}>2.3.$ (2)

Variations in both $H$-band magnitude and color with orientation result in
snapshots for which some lines of sight meet the criteria and some do not. We
define two subsets of lines of sight of the HST-dark subsample: ‘selected’ are
orientations that meet the criteria; ‘not-selected’ are orientations that do
not meet the criteria, where there exist other orientations of the same
snapshot (i.e. the same galaxy and redshift) that do meet the criteria. In
total, we obtain $17$ independent galaxy snapshots that meet the selection
criteria at one or more viewing angles. We simulate a total of $119$ viewing
angles of these sometimes-HST-dark snapshots, where $46$ meet the selection
and $73$ do not. We study the viewing angle dependence in more detail in
Section 3.5.
In Figure 1, we show the positions in the color-magnitude plane of selected
and not-selected lines of sight of snapshots that are classified as HST-dark
from one or more viewing angles. In gray, we show all seven lines of sight of
FIRE snapshots with no HST-dark sightlines. The selected HST-dark subsample
occupies a similar region of this parameter space to the observed sample of
Barrufet et al. (2022), giving confidence in the realism of our simulated
galaxies. Lines of sight in the not-selected sample tend to lie just outside
the boundaries of the selection, either because they are slightly too bright
at $1.6\,\mu\rm{m}$, slightly too blue, or both. Only one simulated galaxy
snapshot is HST-dark from all seven modelled viewing angles. This implies that
HST-dark sources selected in this way are not an entirely unique population;
instead, they are sources viewed from preferential orientations that would be
visible by HST if viewed from other directions. In Section 3, we explore the
physical properties of the selected HST-dark sources and the drivers of this
orientation dependence.

Figure 1: Selection of snapshot orientations detected as HST-dark (dark red),
according to the criteria: $H_{\rm{AB}}-\rm{F444W}_{\rm{AB}}>2.3$ (solid line)
and $H_{\rm{AB}}>27\,\rm{mag}$ (dashed line). Lines of sight that do not meet
this criteria, when one or more lines of sight of the same galaxy snapshot do,
are plotted in blue. These tend to be slightly too blue, slightly too
$H-$bright, or both. Lines of sight for snapshots with no orientations meeting
the criteria are plotted in gray.

###  2.4 Comments on selection criteria

Here, we comment on the effects of each of the selection criteria on the
demographics of the selected population. The $H-$dark criterion favours the
selection of higher redshift sources, whereas the color criterion favours
intermediate redshifts. A looser color criterion (i.e. extending to bluer
colors) would lead to the inclusion of higher redshift $H-$faint sources into
the sample. All of our selected snapshots lie within the redshift range
$4<z<7$, but the redshift distribution and range will depend on the details of
the mass assembly and metal enrichment of the simulated galaxies. While a
modest-sized sample constructed in this way does not allow predictions for the
redshift distribution of observed sources, we can note some general trends.
Firstly, the selection criteria identify simulated sources within the general
redshift range of the HST-dark, IRAC/F444W-bright sources reported by
observational studies (e.g. Wang et al., 2019; Barrufet et al., 2022; Nelson
et al., 2022). The modelled $4.4\,\mu\rm{m}$ emission of our sources is fairly
faint ($\rm{F444W_{AB}}\gtrsim 24.8\,\rm{mag}$) compared to observed HST-dark
sources, but still within the range measured by Barrufet et al. (2022).
Importantly, the OIR-dark, sub-millimeter-bright population tends to extend to
lower redshifts: the median photometric redshift of the $K-$dark SMGs studied
by Smail et al. (2021) is $z=3.44\pm 0.06$, the majority lying in the range
$z=3-4$, with a handful below $z=3$. Our simulated sources become too
$H-$bright below $z\sim 4$, but we would expect that dustier, more sub-
millimeter-bright sources would be more highly attenuated in the $H-$band and
hence meet the $H-$faint criterion (the degeneracy between $A_{V}$ and $z$ is
also discussed in Caputi et al. 2012; Wang et al. 2016; Sun et al. 2021). Due
to the limited dynamic range of halos simulated here, we do not model these
sub-millimeter-bright, highly attenuated sources. Instead, we focus on the
sub-millimeter-faint ($S_{850\,\mu\rm{m}}\lesssim 1\,\rm{mJy}$) HST-dark
population.
Finally, we comment on some modelling choices that affect the predicted fluxes
of our simulated galaxies, and therefore the selection. In this work,
following the methods of Cochrane et al. (2019, 2022, 2023); Cochrane &
Angles-Alcazar (2023), we do not implement any models for unresolved, small-
scale dust. Ma et al. (2019) argues that the FIRE galaxies are sufficiently
well-resolved that this is not necessary. Tests have, however, been performed
using the MAPPINGS III models (Groves et al., 2008) to model emission from the
warm dust associated with the unresolved birth clouds of young star clusters.
Liang et al. (2020) studied the impact of varying two MAPPINGS parameters,
$\log C$ (the HII region compactness) and $f_{\rm{PDR}}$ (the covering
fraction of the associated photo-dissociation regions) on the predicted SED.
They found that the choice of $\log C$ had little impact on the UV-optical
SED, while a higher $f_{\rm{PDR}}$ yields a higher dust optical depth, and
therefore reduced UV-optical fluxes and a steeper dust attenuation curve.
Although the effects on the SED are modest, we note that implementing the
MAPPINGS models with a non-zero $f_{\rm{PDR}}$ could result in more snapshots
meeting the HST-dark criteria (due to both lower $H-$band fluxes and redder
colours).

Figure 2: The redshift evolution of stellar mass (top left), dust mass
(center) and star formation rate (right), for all eight FIRE halos in this
study (grey). Snapshots with HST-dark sightlines are plotted in red. Simulated
galaxies with HST-dark sightlines span redshifts $\sim 4-7$, with a broad
range of stellar masses (from $\sim 10^{9}\,\rm{M_{\odot}}$ to $\sim
10^{11}\,\rm{M_{\odot}}$).

##  3 The nature of HST-dark sources

###  3.1 HST-dark source demographics

In Figure 2, we show the redshift evolution of stellar mass, dust mass and
star formation rate for the eight FIRE halos studied. We briefly summarise the
demographics of snapshots with HST-dark lines of sight here. Stellar masses of
the galaxies in selected snapshots range from $\sim 10^{9}\,\rm{M_{\odot}}$ to
$\sim 10^{11}\,\rm{M_{\odot}}$, with mean
$\log_{10}(M_{\star}/\rm{M_{\odot}})=9.8$. Dust masses range from $\sim
5\times 10^{6}\,\rm{M_{\odot}}$ to $\sim 5\times 10^{7}\,\rm{M_{\odot}}$, with
mean $\log_{10}(M_{\rm{dust}}/\rm{M_{\odot}})=7.0$. Star formation rates
(defined here as instantaneous) range from a few to a few hundred solar masses
per year, with mean $\rm{SFR}=56\,\rm{M_{\odot}\rm{yr}^{-1}}$. These are in
broad agreement with the inferred physical properties of
$4.4/4.5\,\mu\rm{m}$-selected HST-dark sources (Sun et al., 2021; Barrufet et
al., 2022; Nelson et al., 2022). In contrast, the sub-millimeter-selected
$K-$dark galaxies studied by Smail et al. (2021) have higher derived dust
masses ($>10^{8}\,\rm{M_{\odot}}$), with stellar masses $\sim
10^{11}\,\rm{M_{\odot}}$; the majority of those sources lie at $z<4$.

Figure 3: Left: dust attenuation curves for our subsamples of selected and
not-selected sightlines (red and blue) of HST-dark sources, compared to those
of all FIRE-2 snapshots at $4<z<7$ (gray). Solid lines and shaded regions show
median curves and $1\,\sigma$ range for all lines of sight that fall into a
given subsample, respectively. Center: dust attenuation curves as in the left-
hand panel, now normalised by $A_{V}$. We overplot the canonical Calzetti et
al. (2000) law. In both panels, the average attenuation curves of selected and
not-selected lines of sight are similar, with the median not-selected curve
lying below the selected curve by $\Delta A_{V}<0.22$ at all wavelengths.
Right: the distribution of $A_{V}$ values for each of the samples, with dashed
lines showing medians (median $A_{V}=2.8,\,2.6$ for selected and not-selected
orientations, respectively, compared to $A_{V}=0.5$ for the whole FIRE sample
at $z=4-7$). Snapshots that are selected as HST-dark along one or more lines
of sight are distinguished from the overall population by particularly high
$A_{V}$ values.

###  3.2 Dust attenuation of HST-dark sources

By comparing intrinsic and attenuated emission as a function of wavelength, we
can study the effective attenuation curves of the simulated galaxies at each
snapshot. In Figure 3, we plot median dust attenuation $A_{\lambda}$ (left)
and median $V$-band-normalised dust attenuation $A_{\lambda}/A_{V}$ (center)
versus wavelength, for the subsample of lines of sight selected as HST-dark
and red (red line) and non-selected lines of sight of the same snapshots (blue
line). For comparison, we also show curves generated using all snapshots in
the redshift range within which our HST-dark-selected sources fall, $4<z<7$.
It is clear from the left-hand panel that both selected and not-selected lines
of sight of the HST-dark snapshots are significantly more attenuated, on
average, than other sources in the same redshift range. Once normalised by
$A_{V}$, though, the shapes of the attenuation curves are very similar (center
panel) and grayer than the curve derived by Calzetti et al. (2000) (note,
though, that we have used the same dust grain model for all simulated sources
and snapshots: different choices of dust grain mix and size distributions
would change the details of the attenuation and possibly also change the
numbers of simulated sources that meet the HST-dark selection criteria). In
the right-hand panel, we show a histogram of $A_{V}$ for the different
populations. The median $A_{V}$ values for selected and not-selected lines of
sight of HST-dark snapshots are $A_{V}=2.8$ and $2.6$, respectively, while
$A_{V}=0.5$ for FIRE snapshots in the redshift range $z=4-7$.

(a) Distributions of radii containing half the stellar mass (left), SFR over
the last $10\,\rm{Myr}$ (center) and dust mass (right), for sources with HST-
dark sightlines (red) compared to all FIRE-2 snapshots in the redshift range
$4<z<7$ (grey). Snapshots with HST-dark sightlines display extremely compact
stellar mass, star formation ($R_{1/2,\,M_{\star}}$ &
$R_{1/2,\,\rm{SFR}}<0.5\,\rm{kpc}$), and dust mass
($R_{1/2,\,M_{\rm{dust}}}\sim 1\,\rm{kpc}$), compared to other snapshots at
similar redshifts.

(b) Distributions of half-light sizes at observed-frame $1.6\,\mu\rm{m}$
(left), $4.4\,\mu\rm{m}$ (center), and $870\,\mu\rm{m}$ (right) for subsamples
of selected and not-selected sightlines (red and blue), compared to those of
all FIRE-2 snapshots at $4<z<7$ (grey). Differences in the distributions of
half-light radii of selected and not-selected sightlines of snapshots with
HST-dark sightlines are minor. HST-dark snapshots typically display more
compact half-light radii than others at the same redshift at all three
wavelengths, with the largest difference at $870\,\mu\rm{m}$. The half-light
sizes (at $1.6\,\mu\rm{m}$ and $4.4\,\mu\rm{m}$) of HST-dark snapshots and the
general population are more similar than the half-mass and half-SFR radii.
This is because, for the more highly attenuated sources, more of the
$1.6\,\mu\rm{m}$ and $4.4\,\mu\rm{m}$ emission is biased to larger radii due
to central dust obscuration.

Figure 4: Physical (top panel) and observable (bottom panel) radii of HST-dark
galaxies compared to all FIRE snapshots within a similar redshift range.

###  3.3 Physical and light-based sizes of HST-dark sources

####  3.3.1 Physical sizes of HST-dark sources

In Figure 4(a), we present distributions of the half-mass, half-SFR and half-
dust mass sizes of galaxies with and without HST-dark sightlines. These sizes
are calculated from particle data, and are hence mass-weighted rather than
light-weighted. The stellar mass (both of all stars and recently-formed stars)
is extremely compact for the HST-dark sources, with $R_{1/2,\,M_{\star}}$ and
$R_{1/2,\,\rm{SFR_{10\,Myr}}}<0.5\,\rm{kpc}$. Although more extended than the
stellar mass distribution, the dust mass is also compact
($R_{1/2,\,M_{\rm{dust}}}\sim 1\,\rm{kpc}$). For the whole sample at $z=4-7$,
the median $R_{1/2,\,\rm{SFR_{10\,Myr}}}/R_{1/2,\,M_{\rm{dust}}}=0.49$. For
the subsample of sources with one or more lines of sight meeting the HST-dark
criteria, the median
$R_{1/2,\,\rm{SFR_{10\,Myr}}}/R_{1/2,\,M_{\rm{dust}}}=0.14$, with median
$R_{1/2,\,\rm{SFR_{10\,Myr}}}=0.12\,\rm{kpc}$. Thus, our simulated HST-dark
sources feature particularly compact star formation embedded within a more
extended dust mass distribution. This drives large optical depths.

####  3.3.2 Multi-wavelength sizes of HST-dark sources

In Figure 4(b), we show histograms of half-light radii at $1.6\,\mu\rm{m}$,
$4.4\,\mu\rm{m}$ and $870\,\mu\rm{m}$, for selected and not-selected
subsamples, as well as for all snapshots in the redshift range $4<z<7$. Note
that these half-light sizes are calculated using a curve-of-growth technique,
rather than from Sérsic profile fits to PSF-convolved images (as was performed
in Parsotan et al. 2021; Cochrane & Angles-Alcazar 2023). It is immediately
clear that both selected and not-selected lines of sight of HST-dark snapshots
are more compact, on average, than other simulated snapshots within the same
redshift range. At $1.6\,\mu\rm{m}$, median half-light radii are
$0.81\,\rm{kpc}$ and $0.72\,\rm{kpc}$ for selected and not-selected lines of
sight, respectively, compared to $0.90\,\rm{kpc}$ for the whole $4<z<7$
sample. At $4.4\,\mu\rm{m}$, median half-light radii are $0.52\,\rm{kpc}$ and
$0.55\,\rm{kpc}$ for selected and not-selected lines of sight, respectively,
compared to $0.77\,\rm{kpc}$ for the whole $4<z<7$ sample. The compact sizes
of our simulated sources are in line with recent observations: Gómez-Guijarro
et al. (2023) found that the F444W-derived effective radii of optically-
dark/faint galaxies are 30% smaller than the average star-forming galaxy, at
fixed stellar mass and redshift, with most $\lesssim 1.2\,\rm{kpc}$. The
difference between the spatial extent of the $4.4\,\mu\rm{m}$ and
$1.6\,\mu\rm{m}$ emission is greatest for the HST-dark snapshots; this is
driven by preferential dust attenuation in the inner regions of the galaxy
affecting shorter wavelength light most strongly (see Cochrane & Angles-
Alcazar, 2023) and is in qualitative agreement with the recent observational
study of Suess et al. (2022).
The most significant differences in size are seen at longer wavelengths. At
$870\,\mu\rm{m}$, median half-light radii are $0.17\,\rm{kpc}$ for both
selected and not-selected lines of sight, compared to $0.72\,\rm{kpc}$ for the
whole $4<z<7$ sample. The particularly compact dust continuum emission of
observed $K-$dark SMGs was also noted by Smail et al. (2021). In our
simulations, the compact emission is not only driven by a compact dust mass
distribution (see Figure 4(a), right-hand panel) but by particularly compact
star formation (see Figure 4(a), central panel; this effect was also discussed
in Cochrane et al. 2019, who showed that compact dust emission is driven by
steep dust temperature gradients associated with compact star formation).

###  3.4 Understanding the orientation-dependence of HST-dark source selection

In Section 2.3, we noted that for a given simulated galaxy snapshot, there
exist lines of sight that meet the $H-$band magnitude and color criteria and
others that do not. To illustrate this, in Figures 5 and 6 we show two
examples of snapshots selected as HST-dark along just one of the modelled
lines of sight. We plot projected dust mass, recently-formed stars (with age
$<100\,\rm{Myr}$), and predicted emission maps at observed-frame
$1.6\,\mu\rm{m}$ and $4.4\,\mu\rm{m}$. Contours of the projected dust mass are
shown on both emission maps. We also show predicted SEDs and dust attenuation
curves for all lines of sight.
Figure 5 shows galaxy C1 at $z=6.5$. This galaxy is faintest in the observed-
frame $H-$band in the approximately face-on orientation (top panel). This is
the only orientation for which both criteria are met. Viewing the galaxy from
other angles (subsequent panels), the dust is less well-orientated to cover
the emission from the young stars; hence, the galaxy becomes brighter in $H$,
which also leads to a bluer $H-\rm{F444W}$ color. Figure 6 shows galaxy B2 at
$z=5.4$. This galaxy meets the criteria along just one of the modelled lines
of sight; from all others, the galaxy appears too bright in $H$ and also too
blue. In both examples, there are lines of sight from which the $H-$band
emission appears to escape from the edges of the dust distribution. This is in
line with observations of dusty high-redshift galaxies (e.g. Hodge et al.,
2016; Cochrane et al., 2021).

Figure 5: Top panels: Projected distributions of dust mass, young stars, and
predicted emission (intrinsic and observed) along different lines of sight for
galaxy C1 at $z=6.5$. Bottom panels: SED and attenuation curve for the
different lines of sight. Vertical dashed lines indicate $H-$band and
$4.4\,\mu\rm{m}$. Only one of the simulated lines of sight meets both $H-$band
and color selection criteria of HST-dark galaxies. Two others meet the
$H-$band criterion but not the color criterion. Along other lines of sight,
the dust does not cover the young stellar emission as completely; hence, they
appear brighter in the $H$-band and are also too blue.

Figure 6: The same as Fig 5, for galaxy B2 at $z=5.4$. Only one line of sight
meets the HST-dark criteria; the others are brighter in $H$ and hence too
blue.

###  3.5 Implications of line-of-sight variations for comparison with models

One key result of this paper is the clear dependence of HST-dark source
selection on line of sight. Only one of the 17 selected HST-dark galaxies is
selected from all seven modelled orientations; the rest are selected along
some lines of sight but not along others. To investigate this further, we
repeat the radiative transfer calculation for each of the $17$ selected
snapshots, this time modelling $50$ lines of sight, with solid angle evenly
sampled. The number of lines of sight meeting the HST-dark selection criteria
ranges from $1$ to $45$ (see Figure 7). In total, $268$ of $850$ modelled
sightlines meet the criteria ($32\%$). This has implications for the
measurement of number densities of HST-dark sources from observations: for
each galaxy selected as HST-dark, there are, on average, two others with the
same physical properties that would not meet the selection criteria, simply
because of sky orientation. More massive, sub-millimeter-bright sources, not
simulated here, may be obscured along more lines of sight.
The high resolution of the FIRE simulations is key in enabling this result. In
coarser resolution simulations, resolving the detailed geometry of gas and
stars might not be possible. For example, the baryonic particle resolution of
IllustrisTNG-100 is $\sim 40$ times lower than our highest resolution
simulations (Pillepich et al., 2018), and that of SIMBA’s largest box is $\sim
14$ times lower (Davé et al., 2019). Furthermore, simulations that do not
explicitly resolve stellar feedback (e.g. using effective equation of state
models) like IllustrisTNG will artificially suppress small-scale ISM clumping.
Semi-analytic models tend to model galaxy geometry in a simplified manner. For
example, Wang et al. (2019) compared the relative contributions of observed
$H-$dropouts and modelled $H-$dropouts in the l-galaxies SAM (2015 version;
Henriques et al. 2015) to the cosmic SFRD. However, the SAM models stars only
as bulge, disc and intracluster light components. Optical depths are
calculated separately for the diffuse ISM and molecular birth clouds,
following De Lucia & Blaizot (2007). The overall extinction curve is derived
by assuming an ‘inclined slab’ geometry implementation for the diffuse ISM.
Although some geometric effects will be included as a result of the
inclination angle, this is a simplistic implementation that does not allow for
a clumpy ISM (or scattering into the line-of-sight). We caution that coarse-
resolution hydrodynamical simulations and semi-analytic models may not resolve
the crucial dust-star geometry sufficiently for robust comparisons with
observed HST-dark populations.

Figure 7: We re-ran the radiative transfer procedure for each snapshot with
one or more previously-identified HST-dark lines of sight. We used 50 lines of
sight, evenly sampling solid angle. Here, we show the distribution of the
percentage of snapshots that meet the HST-dark selection criteria. This ranges
from $2$ to $90$ per cent.

##  4 Conclusions

Using highly-resolved zoom-in simulations from the FIRE suite, we perform the
first detailed modelling of HST-dark galaxies, a population of galaxies
identified observationally by their faint $H-$band magnitudes and red
$H-\rm{F444W}$ colors. The high resolution of these simulations enables us to
infer realistic distributions of dust and stars for galaxies spanning a range
of stellar masses and redshifts. We run radiative transfer calculations on
eight massive galaxies across $1\lesssim z\lesssim 8$, and make predictions
for the observed-frame UV-FIR emission. Applying selection criteria identical
to those used in observational searches for IRAC/JWST-bright HST-dark
galaxies, we naturally recover $17$ galaxy snapshots that meet HST-dark
selection criteria at one or more viewing angles, all in the redshift range
$z=4-7$. The four most massive halos all pass through a HST-dark stage at some
point in their evolution. This is because those halos build enough stellar and
dust mass early enough to meet the colour selection. We then study the
physical properties of these galaxies in light of the selections applied.
Galaxy snapshots selected as HST-dark show high levels of dust attenuation
($2<A_{V}<4$) compared to others in the same redshift range. Physically, these
high $A_{V}$ values are associated with substantial dust masses
($M_{\rm{dust}}>5\times 10^{6}\,\rm{M_{\odot}}$) and compact star formation
($R_{\rm{SFR\,10\,Myr}}\sim 0.1\,\rm{kpc}$). Modelling the observable emission
along additional lines of sight for each snapshot enables us to test the role
of geometry in HST-dark galaxy selection. Importantly, galaxies that are HST-
dark along some sightlines do not meet the criteria along others. We infer
that the observational selection of HST-dark galaxies is subject to a strong
viewing angle dependence: rather than a particularly special population, HST-
dark galaxies are a subset of high-redshift galaxies viewed along lines of
sight with particularly high dust attenuation. This result has implications
for comparisons of observations with semi-analytic models and coarse-
resolution simulations that do not resolve the detailed geometry of stars and
dust.

## Acknowledgements

We thank the anonymous reviewer for helpful suggestions on an early version of
this paper. The Flatiron Institute is supported by the Simons Foundation. The
simulations presented in this work were run on the Flatiron Institute’s
research computing facilities (Gordon-Simons, Popeye, and Iron compute
clusters), supported by the Simons Foundation. DAA acknowledges support by NSF
grants AST-2009687 and AST-2108944, CXO grant TM2-23006X, Simons Foundation
Award CCA-1018464, and Cottrell Scholar Award CS-CSA-2023-028 by the Research
Corporation for Science Advancement. FC acknowledges support from a UKRI
Frontier Research Guarantee Grant (PI Cullen; grant reference EP/X021025/1).

## Data availability

The FIRE-2 simulations are publicly available (Wetzel et al., 2023) at
http://flathub.flatironinstitute.org/fire. Additional FIRE simulation data is
available at https://fire.northwestern.edu/data. A public version of the Gizmo
code is available at http://www.tapir.caltech.edu/~phopkins/Site/GIZMO.html.

## References

  * Anglés-Alcázar et al. (2017) Anglés-Alcázar, D., Faucher-Giguère, C. A., Quataert, E., et al. 2017, MNRAS, 472, L109, doi: 10.1093/mnrasl/slx161
  * Baes et al. (2011) Baes, M., Verstappen, J., De Looze, I., et al. 2011, ApJS, 196, doi: 10.1088/0067-0049/196/2/22
  * Barrufet et al. (2022) Barrufet, L., Oesch, P. A., Weibel, A., et al. 2022, arxiv:2207.14733.  https://arxiv.org/abs/2207.14733
  * Behiri et al. (2023) Behiri, M., Talia, M., Cimatti, A., et al. 2023.  https://arxiv.org/abs/2309.00050
  * Bertoldi et al. (2000) Bertoldi, F., Carilli, C. L., Menten, K. M., et al. 2000, A&A, 360, 92
  * Bruzual & Charlot (2003) Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000.  https://arxiv.org/abs/0309134v1
  * Calzetti et al. (2000) Calzetti, D., Armus, L., Bohlin, R., et al. 2000, ApJ, 533, 682
  * Camps & Baes (2015) Camps, P., & Baes, M. 2015, Astronomy and Computing, 9, 20, doi: 10.1016/j.ascom.2014.10.004
  * Caputi et al. (2012) Caputi, K. I., Dunlop, J. S., McLure, R. J., et al. 2012, ApJL, 750, 1, doi: 10.1088/2041-8205/750/1/L20
  * Casey et al. (2021) Casey, C. M., Zavala, J. A., Manning, S. M., et al. 2021, ApJ, 923, 215, doi: 10.3847/1538-4357/ac2eb4
  * Cimatti et al. (1998) Cimatti, A., Andreani, P., Röttgering, H., & Tilanus, R. 1998, Nature, 392, 895, doi: 10.1038/31872
  * Cochrane & Angles-Alcazar (2023) Cochrane, R., & Angles-Alcazar, D. 2023, in prep.
  * Cochrane et al. (2021) Cochrane, R. K., Best, P. N., Smail, I., et al. 2021, MNRAS, 503, 2622
  * Cochrane et al. (2022) Cochrane, R. K., Hayward, C. C., & Angles-Alcazar, D. 2022, ApJL, 939, L27
  * Cochrane et al. (2023) Cochrane, R. K., Hayward, C. C., Angles-Alcazar, D., & Somerville, R. S. 2023, MNRAS, 518, 5522
  * Cochrane et al. (2019) Cochrane, R. K., Hayward, C. C., Anglés-Alcázar, D., et al. 2019, MNRAS, 488, 1779, doi: 10.1093/mnras/stz1736
  * Cowie et al. (2009) Cowie, L. L., Barger, A. J., Wang, W. H., & Williams, J. P. 2009, ApJ, 697, 1, doi: 10.1088/0004-637X/697/2/L122
  * Dannerbauer et al. (2002) Dannerbauer, H., Lehnert, M. D., Lutz, D., et al. 2002, ApJ, 573, 473
  * Davé et al. (2019) Davé, R., Anglés-Alcázar, D., Narayanan, D., et al. 2019, MNRAS, 486, 2827, doi: 10.1093/mnras/stz937
  * De Lucia & Blaizot (2007) De Lucia, G., & Blaizot, J. 2007, MNRAS, 375, 2, doi: 10.1111/j.1365-2966.2006.11287.x
  * Dey et al. (1999) Dey, A., Graham, J. R., Ivison, R. J., et al. 1999, ApJ, 519, 610, doi: 10.1086/307395
  * Draine & Salpeter (1979) Draine, B. T., & Salpeter, E. E. 1979, ApJ, 231, 77, doi: 10.1086/157165
  * Dudzevičiūtė et al. (2020) Dudzevičiūtė, U., Smail, I., Swinbank, A. M., et al. 2020, MNRAS, 494, 3828.  https://arxiv.org/abs/arXiv:1910.07524v1
  * Dunlop et al. (2007) Dunlop, J. S., Cirasuolo, M., & McLure, R. J. 2007, MNRAS, 376, 1054, doi: 10.1111/j.1365-2966.2007.11453.x
  * Dwek (1998) Dwek, E. 1998, ApJ, 1, 643
  * Elston et al. (1988) Elston, R., Rieke, G. H., & Rieke, M. J. 1988, ApJ, 331, L77
  * Enia et al. (2022) Enia, A., Talia, M., Pozzi, F., et al. 2022, ApJ, 927, 204, doi: 10.3847/1538-4357/ac51ca
  * Fazio et al. (2004) Fazio, G. G., Hora, J. L., Allen, L. E., et al. 2004, ApJS, 154, 10, doi: 10.1086/422843
  * Feldmann et al. (2016) Feldmann, R., Hopkins, P. F., Quataert, E., Faucher-Giguère, C. A., & Kerěs, D. 2016, MNRAS, 458, L14, doi: 10.1093/mnrasl/slw014
  * Feldmann et al. (2017) Feldmann, R., Quataert, E., Hopkins, P. F., Faucher-Giguére, C. A., & Kereš, D. 2017, MNRAS, 470, 1050, doi: 10.1093/mnras/stx1120
  * Finkelstein et al. (2022) Finkelstein, S. L., Bagley, M., Song, M., et al. 2022, ApJ, 928, 52, doi: 10.3847/1538-4357/ac3aed
  * Franco et al. (2018) Franco, M., Elbaz, D., Béthermin, M., et al. 2018, A&A, 620, A152, doi: 10.1051/0004-6361/201832928
  * Frayer et al. (2004) Frayer, D. T., Reddy, N. A., Armus, L., et al. 2004, ApJ, 127, 728
  * Frayer et al. (2000) Frayer, D. T., Smail, I., Ivison, R. J., & Scoville, N. Z. 2000, ApJ, 120, 1668
  * Fujimoto et al. (2023) Fujimoto, S., Bezanson, R., Labbe, I., et al. 2023, eprint arXiv:2309.07834.  https://arxiv.org/abs/2309.07834
  * Gómez-Guijarro et al. (2023) Gómez-Guijarro, C., Magnelli, B., Elbaz, D., et al. 2023, A&A, 677, 1, doi: 10.1051/0004-6361/202346673
  * Groves et al. (2008) Groves, B., Dopita, M. A., Sutherland, R. S., et al. 2008, ApJS, 176, 438, doi: 10.1086/528711
  * Henriques et al. (2015) Henriques, B. M. B., White, S. D. M., Thomas, P. A., et al. 2015, MNRAS, 451, 2663, doi: 10.1093/mnras/stv705
  * Hodge et al. (2013) Hodge, J. A., Karim, A., Smail, I., et al. 2013, ApJ, 768, 91, doi: 10.1088/0004-637X/768/1/91
  * Hodge et al. (2016) Hodge, J. A., Swinbank, A. M., Simpson, J. M., et al. 2016, doi: 10.3847/1538-4357/833/1/103
  * Hopkins (2015) Hopkins, P. F. 2015, MNRAS, 450, 53, doi: 10.1093/mnras/stv195
  * Hopkins et al. (2014) Hopkins, P. F., Kereš, D., Oñorbe, J., et al. 2014, MNRAS, 445, 581, doi: 10.1093/mnras/stu1738
  * Hopkins et al. (2013) Hopkins, P. F., Narayanan, D., & Murray, N. 2013, MNRAS, 432, 2647, doi: 10.1093/mnras/stt723
  * Hopkins et al. (2018) Hopkins, P. F., Wetzel, A., Keres, D., et al. 2018, MNRAS, 480, 800.  https://arxiv.org/abs/1702.06148
  * Hopkins et al. (2023) Hopkins, P. F., Wetzel, A., Wheeler, C., et al. 2023, MNRAS, 519, 3154.  https://arxiv.org/abs/2203.00040
  * Hu & Ridgway (1994) Hu, E. M., & Ridgway, S. E. 1994, ApJ, 107
  * Huang et al. (2011) Huang, J. S., Zheng, X. Z., Rigopoulou, D., et al. 2011, ApJL, 742, 6, doi: 10.1088/2041-8205/742/1/L13
  * Hughes et al. (1998) Hughes, D. H., Serjeant, S., Dunlop, J., et al. 1998, Nature, 394, 241, doi: 10.1038/28328
  * James et al. (2002) James, A., Dunne, L., Eales, S., & Edmunds, M. G. 2002, MNRAS, 335, 753, doi: 10.1046/j.1365-8711.2002.05660.x
  * Karim et al. (2013) Karim, A., Swinbank, A. M., Hodge, J. A., et al. 2013, MNRAS, 432, 2, doi: 10.1093/mnras/stt196
  * Kondapally et al. (2021) Kondapally, R., Best, P. N., Hardcastle, M. J., et al. 2021, A&A, 648, A3, doi: 10.1051/0004-6361/202038813
  * Leitherer et al. (1999) Leitherer, C., Schaerer, D., Goldader, J. D., et al. 1999, ApJS, 123, 3
  * Liang et al. (2020) Liang, L., Feldmann, R., Hayward, C. C., et al. 2020, 32, 1.  https://arxiv.org/abs/2009.13522
  * Lilly et al. (1999) Lilly, S., Eales, S. A., Gear, W. K., et al. 1999, ApJ, 518, 614
  * Ma et al. (2019) Ma, X., Hayward, C. C., Casey, C. M., et al. 2019, MNRAS, 487, 1844.  https://arxiv.org/abs/1902.10152
  * Manning et al. (2022) Manning, S. M., Casey, C. M., Zavala, J. A., et al. 2022, ApJ, 925, 23, doi: 10.3847/1538-4357/ac366a
  * Nelson et al. (2022) Nelson, E. J., Suess, K. A., Bezanson, R. S., et al. 2022, arxiv:2208.01630.  https://arxiv.org/abs/2207.14733
  * Parsotan et al. (2021) Parsotan, T., Cochrane, R. K., Hayward, C. C., et al. 2021, MNRAS, 501, 1591, doi: 10.1093/mnras/staa3765
  * Pérez-González et al. (2023) Pérez-González, P. G., Barro, G., Annunziatella, M., et al. 2023, ApJL, 946, L16.  https://arxiv.org/abs/2211.00045
  * Pillepich et al. (2018) Pillepich, A., Nelson, D., Hernquist, L., et al. 2018, MNRAS, 475, 648, doi: 10.1093/mnras/stx3112
  * Rémy-Ruyer et al. (2014) Rémy-Ruyer, A., Madden, S. C., Galliano, F., et al. 2014, A&A, 563, 31, doi: 10.1051/0004-6361/201322803
  * Rodighiero et al. (2023) Rodighiero, G., Bisigello, L., Iani, E., et al. 2023, MNRAS, 518, L19, doi: 10.1093/mnrasl/slac115
  * Simpson et al. (2017a) Simpson, J. M., Smail, I., Wang, W.-H., Riechers, D., & Dunlop, J. S. 2017a, ApJL, 844, L10.  https://arxiv.org/abs/arXiv:1704.03868v1
  * Simpson et al. (2014) Simpson, J. M., Swinbank, A. M., Smail, I., et al. 2014, ApJ, 788, doi: 10.1088/0004-637X/788/2/125
  * Simpson et al. (2015) Simpson, J. M., Smail, I., Swinbank, A. M., et al. 2015, ApJ, 799, 81, doi: 10.1088/0004-637X/799/1/81
  * Simpson et al. (2017b) —. 2017b, ApJ, 839, 58, doi: 10.3847/1538-4357/aa65d0
  * Simpson et al. (2020) Simpson, J. M., Smail, I., Dudzevičiūtė, U., et al. 2020, MNRAS, 495, 3409, doi: 10.1093/mnras/staa1345
  * Smail et al. (2004) Smail, I., Chapman, S. C., Blain, A. W., & Ivison, R. J. 2004, ApJ, 617, 64, doi: 10.1086/425171
  * Smail et al. (2002a) Smail, I., Ivison, R., Blain, A., & Kneib, J.-P. 2002a, MNRAS, 331, 495
  * Smail et al. (2000) Smail, I., Ivison, R., Owen, F. N., Blain, A. W., & Kneib, J.-P. 2000, ApJ, 528, 612
  * Smail et al. (1999) Smail, I., Ivison, R. J., Kneib, J. P., et al. 1999, MNRAS, 308, 1061, doi: 10.1046/j.1365-8711.1999.02819.x
  * Smail et al. (2002b) Smail, I., Owen, F. N., Morrison, G. E., et al. 2002b, ApJ, 581, 844, doi: 10.1086/344440
  * Smail et al. (2021) Smail, I., Dudzevičiūtė, U., Stach, S. M., et al. 2021, MNRAS, 502, 3426\.  https://arxiv.org/abs/2010.02250
  * Smail et al. (2023) Smail, I., Dudzeviciute, U., Gurwell, M., et al. 2023, eprint arXiv:2306.16039.  https://arxiv.org/abs/2306.16039
  * Springel et al. (2005) Springel, V., White, S. D. M., Jenkins, A., et al. 2005, Nature, 435, 629, doi: 10.1038/nature03597
  * Stach et al. (2018) Stach, S. M., Smail, I., Swinbank, A. M., et al. 2018, ApJ, 860, 161, doi: 10.3847/1538-4357/aac5e5
  * Stach et al. (2019) Stach, S. M., Dudzevičiūtė, U., Smail, I., et al. 2019, MNRAS, 487, 4648, doi: 10.1093/mnras/stz1536
  * Suess et al. (2022) Suess, K. A., Bezanson, R., Nelson, E. J., et al. 2022, ApJ, 937, L33, doi: 10.3847/2041-8213/ac8e06
  * Sun et al. (2021) Sun, F., Egami, E., Pérez-González, P. G., et al. 2021, ApJ, 922, 114, doi: 10.3847/1538-4357/ac2578
  * Talia et al. (2021) Talia, M., Cimatti, A., Giulietti, M., et al. 2021, ApJ, 909, 23, doi: 10.3847/1538-4357/abd6e3
  * Tielens et al. (1994) Tielens, A. G. G. M., McKee, C. F., Seab, C. G., & Hollenbach, D. J. 1994, ApJ, 431, 321, doi: 10.1086/174488
  * van der Vlugt et al. (2022) van der Vlugt, D., Hodge, J. A., Algera, H. S. B., et al. 2022, ApJ, 941, 10, doi: 10.3847/1538-4357/ac99db
  * Walter et al. (2012) Walter, F., Decarli, R., Carilli, C., et al. 2012, Nature, 486, 233.  https://arxiv.org/abs/arXiv:1011.1669v3
  * Wang et al. (2016) Wang, T., Elbaz, D., Schreiber, C., et al. 2016, ApJ, 816, 84, doi: 10.3847/0004-637x/816/2/84
  * Wang et al. (2019) Wang, T., Schreiber, C., Elbaz, D., et al. 2019, Nature, 572, 211, doi: 10.1038/s41586-019-1452-4
  * Wang et al. (2007) Wang, W.-H., Cowie, L. L., van Saders, J., Barger, A. J., & Williams, J. P. 2007, ApJ, 670, L89
  * Weingartner & Draine (2001) Weingartner, J. C., & Draine, B. T. 2001, ApJ, 548, 296, doi: 10.1086/318651
  * Wetzel et al. (2023) Wetzel, A., Hayward, C. C., Sanderson, R. E., et al. 2023, ApJS, 265, 44.  https://arxiv.org/abs/2202.06969
  * Williams et al. (2019) Williams, C. C., Labbe, I., Spilker, J., et al. 2019, ApJ, 884, 154, doi: 10.3847/1538-4357/ab44aa
  * Wilson et al. (2004) Wilson, G., Huang, J., Perez-Gonzalez, P. G., et al. 2004, ApJS, 154, 107, doi: 10.1086/422716
  * Xiao et al. (2023) Xiao, M., Oesch, P., Elbaz, D., et al. 2023, arXiv:2309.02492, 17.  https://arxiv.org/abs/2309.02492
  * Yamaguchi et al. (2019) Yamaguchi, Y., Kohno, K., Hatsukade, B., et al. 2019, ApJ, 878, 73, doi: 10.3847/1538-4357/ab0d22