Disappearing Galaxies: The Orientation Dependence of JWST-bright, HST-dark, Star-forming Galaxy Selection

Galaxies that are invisible in deep optical–near-infrared 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 ≳ 4. The depths now achievable with JWST at the longest wavelengths probed by the Hubble Space Telescope (HST), coupled with the transformative resolution at longer wavelengths, are already enabling detailed, spatially resolved characterization 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 project, with multiwavelength images along several lines of sight forward-modeled using radiative transfer. We naturally recover a population of modeled sources that meet commonly used selection criteria (H AB > 27 mag and H AB − 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,4.4 μm ≲ 1 kpc). Orientation is very important: for all but one of the 17 simulated galaxy snapshots with HST-dark sight lines, there exist other sight lines 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 semianalytic 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.

Other populations of OIR-dark galaxies have been selected using optical-infrared colors.Galaxies first identified by red  −  colors and faint −band fluxes were called 'Extremely Red Objects' (EROs; Elston et al. 1988; Hu & Ridgway 1 email: rc3660@columbia.edu1994); 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 −band imaging and IRAC CH2:  AB − [4.5] > 2.25 mag.62% of sources selected in this way were detected down to  870 m > 0.6 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  ≳ 3, though estimates of the contribution vary (e.g.Wang et al. 2016 calculated that −dark sources contribute 15 − 25 % of the SFRD at  = 4 − 5, while Sun et al. 2021  estimated 8 +8  −4 % at  = 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 − or −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  AB > 27 mag sources predicted by an early version of the lgalaxies semi-analytic model (SAM) (Henriques et al. 2015) lies ∼ 2 orders of magnitude below the observed number density.They speculated that our understanding of massivegalaxy 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. 2014Simpson et al. , 2015Simpson et al. , 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. 2018Stach et al. , 2019;;Dudzevičiūtė et al. 2020), Smail et al. (2021) found that 15 ± 2 % of the bright ( 850 m > 3.6 mJy) SMG population is fainter than  AB = 25.3 mag.These −faint galaxies tend to have higher redshifts and dust attenuation than their −detected counterparts (see also the high photometric redshifts derived for the −dark, 1.1 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 (FWHM ∼ 2 ′′ ) 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-Guĳarro 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  = 3 − 6 ( AB > 27 &  AB − 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 (  ∼ 2), massive ( ★ ∼ 10 9 − 10 10.5 M ⊙ ) star-forming galaxies at  ∼ 2 − 8 that lie approximately along the main sequence.Using the same data, Nelson et al. (2022) identified a population of spatially extended ( , 4.4 m > 0.17 ′′ ), 4.4 m-bright, −faint sources, with higher inferred stellar masses, the majority in the range  ★ = 10 10−11 M ⊙ (using EAzY, the median stellar mass is  ★ = 10 10.5 M ⊙ , but they also find systematic differences between SED fitting codes).These 4.4 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 -faint and whole sample of Smail et al. 2021 are  ★ = 10 11.10±0.04M ⊙ and  ★ = 10 11.00±0.06M ⊙ , respectively, and some of the most extreme sub-millimeter bright sources are likely even more massive (e.g. the stellar mass derived for the  850 m = 15.3±0.4mJy,  = 4.26 source presented by Smail et al. 2023 is 10 11.8 M ⊙ ).
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. 2014Hopkins et al. , 2018Hopkins et al. , 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.

FIRE simulations
The FIRE1 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 -body+hydrodynamics code GIZMO2 (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 finestructure processes are modelled from  = 10 K to  = 10 10 K. Star particles form from locally self-gravitating, molecular, Jeans unstable gas above a minimum hydrogen number density   ≥ 1000 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. (2016Feldmann et al. ( , 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  halo ∼ 10 12.5 M ⊙ at  = 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  = 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 × 10 4 M ⊙ , 2.7 × 10 5 M ⊙ , and 2.2 × 10 6 M ⊙ , for A-, B-, and C-series haloes, respectively.For dark matter particles, the respective mass resolutions are: 1.7 × 10 5 M ⊙ , 1.4 × 10 6 M ⊙ , and 1.1 × 10 7 M ⊙ .Convergence tests for these simulations are presented in Cochrane et al. (2022) (see Appendix B).

Modelling observable emission
We model spectral energy distributions (SEDs) and multiwavelength emission maps at every snapshot, between  ∼ 8 and  ∼ 1, for each of the eight simulated halos (in total, 1736 snapshots).Following Cochrane et al. (2019Cochrane et al. ( , 2022Cochrane et al. ( , 2023)); Cochrane & Angles-Alcazar (2023) we use the skirt3 radiative transfer code (version 8; Baes et al. 2011;Camps & Baes 2015) to make predictions for emission between restframe ultraviolet (UV) and far-infrared (FIR) wavelengths along seven lines of sight.For all snapshots, gas and star particles within 0.1 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 reason-3 http://www.skirt.ugent.beable assumption for enriched, massive galaxies like these.We have checked that using the metallicity-dependent dustto-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 K (Draine & Salpeter 1979;Tielens et al. 1994).Following Cochrane et al. (2019Cochrane et al. ( , 2022Cochrane et al. ( , 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 ∼ 100 wavelengths modelled, for seven lines of sight.We transform flux densities in Jy to AB magnitudes (using a zero-point of 3631 Jy).

HST-dark source selection
Following Barrufet et al. (2022), we select sightlines from our snapshots that meet the following −band and color criteria: and Variations in both -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 HSTdark 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 AB all FIRE snapshots subsample, selected subsample, not selected Figure 1: Selection of snapshot orientations detected as HSTdark (dark red), according to the criteria:  AB − F444W AB > 2.3 (solid line) and  AB > 27 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 −bright, or both.Lines of sight for snapshots with no orientations meeting the criteria are plotted in gray. 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 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.

Comments on selection criteria
Here, we comment on the effects of each of the selection criteria on the demographics of the selected population.The −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 −faint sources into the sample.All of our selected snapshots lie within the redshift range 4 <  < 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 modestsized 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 m emission of our sources is fairly faint (F444W AB ≳ 24.8 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 −dark SMGs studied by Smail et al. (2021) is  = 3.44 ± 0.06, the majority lying in the range  = 3 − 4, with a handful below  = 3.Our simulated sources become too −bright below  ∼ 4, but we would expect that dustier, more sub-millimeter-bright sources would be more highly attenuated in the −band and hence meet the −faint criterion (the degeneracy between   and  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 ( 850 m ≲ 1 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. (2019Cochrane et al. ( , 2022Cochrane et al. ( , 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  (the HII region compactness) and  PDR (the covering fraction of the associated photo-dissociation regions) on the predicted SED.They found that the choice of log  had little impact on the UV-optical SED, while a higher  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  PDR could result in more snapshots meeting the HST-dark criteria (due to both lower −band fluxes and redder colours).

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 ∼ 10 9 M ⊙ to ∼ 10 11 M ⊙ , with mean log 10 ( ★ /M ⊙ ) = 9.8.Dust masses range from ∼ 5 × 10 6 M ⊙ to ∼ 5 × 10 7 M ⊙ , with mean 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 ∼ 4 − 7, with a broad range of stellar masses (from ∼ 10 9 M ⊙ to ∼ 10 11 M ⊙ ).log 10 ( dust /M ⊙ ) = 7.0.Star formation rates (defined here as instantaneous) range from a few to a few hundred solar masses per year, with mean SFR = 56 M ⊙ yr −1 .These are in broad agreement with the inferred physical properties of 4.4/4.5 m-selected HST-dark sources (Sun et al. 2021;Barrufet et al. 2022;Nelson et al. 2022).In contrast, the submillimeter-selected −dark galaxies studied by Smail et al. (2021) have higher derived dust masses (> 10 8 M ⊙ ), with stellar masses ∼ 10 11 M ⊙ ; the majority of those sources lie at  < 4.

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   (left) and median -band-normalised dust attenuation   /  (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 <  < 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   , 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   for the different populations.The median   values for selected and not-selected lines of sight of HST-dark snapshots are   = 2.8 and 2.6, respectively, while   = 0.5 for FIRE snapshots in the redshift range  = 4 − 7.

Physical sizes of HST-dark sources
In Figure 4a, 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  1/2,  ★ and  1/2, SFR 10 Myr < 0.5 kpc.Although more extended than the stellar mass distribution, the dust mass is also compact ( 1/2,  dust ∼ 1 kpc).For the whole sample at  = 4 − 7, the median  1/2, SFR 10 Myr / 1/2,  dust = 0.49.For the subsample of sources with one or more lines of sight meeting the HST-dark criteria, the median  1/2, SFR 10 Myr / 1/2,  dust = 0.14, with median  1/2, SFR 10 Myr = 0.12 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.

Multi-wavelength sizes of HST-dark sources
In Figure 4b, we show histograms of half-light radii at 1.6 m, 4.4 m and 870 m, for selected and not-selected subsamples, as well as for all snapshots in the redshift range 4 <  < 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 m, median half-light radii are 0.81 kpc and 0.72 kpc for selected and not-selected lines of sight, respectively, compared to 0.90 kpc for the whole 4 <  < 7 sample.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 <  < 7 (gray).Solid lines and shaded regions show median curves and 1  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   .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 Δ  < 0.22 at all wavelengths.Right: the distribution of   values for each of the samples, with dashed lines showing medians (median   = 2.8, 2.6 for selected and not-selected orientations, respectively, compared to   = 0.5 for the whole FIRE sample at  = 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   values.
for selected and not-selected lines of sight, respectively, compared to 0.77 kpc for the whole 4 <  < 7 sample.The compact sizes of our simulated sources are in line with recent observations: Gómez-Guĳarro 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 ≲ 1.2 kpc.The difference between the spatial extent of the 4.4 m and 1.6 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 m, median half-light radii are 0.17 kpc for both selected and not-selected lines of sight, compared to 0.72 kpc for the whole 4 <  < 7 sample.The particularly compact dust continuum emission of observed −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 4a, right-hand panel) but by particularly compact star formation (see Figure 4a, 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).

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 −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 Myr), and predicted emission maps at observed-frame 1.6 m and 4.4 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  = 6.5.This galaxy is faintest in the observed-frame −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 , which also leads to a bluer  − F444W color.Figure 6 shows galaxy B2 at  = 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  and also too blue.In both examples, there are lines of sight from which the −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).(b) Distributions of half-light sizes at observed-frame 1.6 m (left), 4.4 m (center), and 870 m (right) for subsamples of selected and notselected sightlines (red and blue), compared to those of all FIRE-2 snapshots at 4 <  < 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 m.The half-light sizes (at 1.6 m and 4.4 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 m and 4.4 m emission is biased to larger radii due to central dust obscuration.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.

Implications of line-of-sight variations
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 ∼ 40 times lower than our highest resolution simulations (Pillepich et al. 2018), and that of SIMBA's largest box is ∼ 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 −dropouts and modelled −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.

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 −band magnitudes and red  − 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 ≲  ≲ 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  = 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 <   < 4) compared to others in the same redshift range.Physically, these high   values are associated with substantial dust masses ( dust > 5 × 10 6 M ⊙ ) and compact star formation ( SFR 10 Myr ∼ 0.1 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 semianalytic models and coarse-resolution simulations that do not resolve the detailed geometry of stars and dust.
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 ∼ 4 − 7, with a broad range of stellar masses (from ∼ 10 9 M ⊙ to ∼ 10 11 M ⊙ ).

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

Figure 5 :
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  = 6.5.Bottom panels: SED and attenuation curve for the different lines of sight.Vertical dashed lines indicate −band and 4.4 m.Only one of the simulated lines of sight meets both −band and color selection criteria of HST-dark galaxies.Two others meet the −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 -band and are also too blue.

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

Figure 7 :
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.