Predictions for the spatial distribution of the dust continuum emission in 1<z<5 star-forming galaxies

TL;DR

This study uses radiative transfer simulations to predict the spatial distribution of dust emission in high-redshift star-forming galaxies, revealing wavelength-dependent morphologies and sizes consistent with observations.

Contribution

First detailed simulation-based predictions of dust continuum spatial distribution in galaxies at 1<z<5, linking morphology to dust and star formation properties.

Findings

Dust emission is more compact than cold gas and dust mass but more extended than stars.

Simulated 870-$b5$m sizes match ALMA observations.

Dust morphology varies significantly across wavelengths.

Abstract

We present the first detailed study of the spatially-resolved dust continuum emission of simulated galaxies at 1<z<5. We run the radiative transfer code SKIRT on a sample of submillimeter-bright galaxies drawn from the Feedback in Realistic Environments (FIRE) project. These simulated galaxies reach Milky Way masses by z=2. Our modelling provides predictions for the full rest-frame far-ultraviolet-to-far-infrared spectral energy distributions of these simulated galaxies, as well as 25-pc-resolution maps of their emission across the wavelength spectrum. The derived morphologies are notably different in different wavebands, with the same galaxy often appearing clumpy and extended in the far-ultraviolet yet an ordered spiral at far-infrared wavelengths. The observed-frame 870-$\mu$m half-light radii of our FIRE-2 galaxies are ~0.5-4kpc, consistent with existing ALMA observations of galaxies with similarly high redshifts and stellar masses. In both simulated and observed galaxies, the dust continuum emission is generally more compact than the cold gas and the dust mass, but more extended than the stellar component. The most extreme cases of compact dust emission seem to be driven by particularly compact recent star-formation, which generates steep dust temperature gradients. Our results confirm that the spatial extent of the dust continuum emission is sensitive to both the dust mass and SFR distributions.

Figures

Tables

Table 1. Table 1 : Properties of the four simulated FIRE-2 haloes at z = 2 𝑧 2 z=2 . M halo subscript 𝑀 halo M_{\rm{halo}} denotes the total mass of the central halo at z = 2 𝑧 2 z=2 . M ∗ subscript 𝑀 M_{\rm{*}} , M gas subscript 𝑀 gas M_{\rm{gas}} , and SFR SFR \rm{SFR} denote the stellar mass, gas mass, and star-formation rate of the halo’s central galaxy at z = 2 𝑧 2 z=2 , all calculated within 0.1 ​ R vir 0.1 subscript 𝑅 vir 0.1R_{\rm{vir}} . r 1 / 2 subscript 𝑟 1 2 r_{1/2} is the half-mass radius, calculated using the stellar mass within 0.1 ​ R vir 0.1 subscript 𝑅 vir 0.1R_{\rm{vir}} .

Name	${\log }_{10}{M}_{\mathrm{halo}}/M_{\odot }$	${\log }_{10}{M}_{\ast }/M_{\odot }$	${\log }_{10}{M}_{\mathrm{gas}}/M_{\odot }$	$\mathrm{SFR}/M_{\odot }{\mathrm{yr}}^{-1}$	$r_{1/2}/\mathrm{kpc}$
A1	$12.45$	$11.24$	$10.22$	65	0.73
A2	$12.56$	$11.46$	$10.51$	168	0.98
A4	$12.49$	$11.10$	$10.33$	63	0.81
A8	$12.41$	$10.85$	$10.45$	88	0.91

Equations

$$f_{850}=0.81\,\rm{mJy}\times\Bigg{(}\frac{SFR}{100M_{\odot}yr^{-1}}\Bigg{)}^{0.43}\Bigg{(}\frac{M_{\rm{dust}}}{10^{8}M_{\odot}}\Bigg{)}^{0.54},$$

Full text

††pubyear: 2019††volume: in press††pagerange: Predictions for the spatial distribution of the dust continuum emission in $\mathbf{1<z<5}$ star-forming galaxies–References

1 Introduction

Observations suggest that the physical properties of star-forming (SF) galaxies at the peak of cosmic star-formation ($z\sim 2$) differ greatly from those of the ordered disks and ellipticals in the local Universe. Galaxies at these redshifts display high star-formation rates, believed to be driven by large molecular gas reservoirs (Tacconi et al., 2010; Tacconi et al., 2013, 2017; Papovich et al., 2016; Falgarone et al., 2017; Jiménez-Andrade et al., 2018) that arise due to steady accretion of cold gas along filaments of the cosmic web (Kereš et al., 2005; Dekel et al., 2009; Faucher-Giguère et al., 2011; Martin et al., 2016; Kleiner et al., 2017). Structurally, high-redshift galaxies are less ordered than their low-redshift counterparts, with star-formation taking place within turbulent disks (Genzel et al., 2008; Kassin et al., 2012; Guo et al., 2015; Tadaki et al., 2018) that often harbour massive ultraviolet (UV)-bright clumps (Elmegreen et al., 2013; Guo et al., 2017; Soto et al., 2017). However, our measurements of the high-redshift Universe are largely reliant on data at rest-frame optical and UV wavelengths, which can be biased towards the least dust-obscured galaxies.

Only in the last few years have new facilities such as ALMA had the resolving power to resolve and probe the morphology of longer-wavelength emission from highly star-forming galaxies. The angular resolution of previous instruments such as SCUBA (Holland et al., 1999) was low, so it was not possible to determine the structural properties of high-redshift galaxies. Source confusion has also been a hindrance in the identification of fainter sources (e.g. with Herschel; Oliver et al., 2012; Scudder et al., 2016). ALMA has the potential to be particularly fruitful in identifying high-redshift galaxies, due to the so-called ‘negative k-correction’ at rest-frame far-infrared (FIR) fluxes (observed $\lambda\sim 850$-µm flux from galaxies with similar intrinsic spectra remains approximately constant across the redshift range $z\sim 1-6$, as we trace further up the Rayleigh-Jeans tail at higher redshift). The most interesting physical insights will likely come from the combination of these new millimeter/sub-millimeter (mm/sub-mm) data with shorter-wavelength imaging. To this end, Dunlop et al. (2017) present $1.3\rm{mm}$ ALMA imaging of the Hubble UltraDeep Field (HUDF), previously mapped with the Wide Field Camera 3/IR on the Hubble Space Telescope (HST) to an unprecedented $5\sigma$ depth of $30\rm{AB\,mag}$ (Bouwens et al., 2010; Oesch et al., 2010; Illingworth et al., 2013; Dunlop et al., 2013; Ellis et al., 2013), and also with the HST Advanced Camera for Surveys (ACS; Beckwith et al., 2006), over an area of $4.5\rm{arcmin}^{2}$. Combining these new ALMA data with Herschel and Spitzer 24-µm  photometry and fitting to a template spectral energy distribution (SED), Dunlop et al. (2017) find that $\sim 85\%$ of the total star-formation at $z\sim 2$ is enshrouded in dust (see also Bourne et al. 2017; Whitaker et al. 2017). They show that for high-mass galaxies ($M_{*}>2\times 10^{10}M_{\odot}$), which host $\sim 65\%$ of the total star-formation at this epoch, the star-formation rate derived from long-wavelength emission is an extraordinary $200$ times that derived from unobscured light. Bowler et al. (2018) demonstrate that dust-obscured star-formation could be substantial even as early as $z\sim 7$.

In addition to finding substantially different measurements of star-formation rates compared to dust-uncorrected short-wavelength light, studies of galaxies at longer wavelengths also present a different view of the morphologies of high-redshift galaxies. Unobscured emission (probed at short wavelengths) tends to be significantly more extended and clumpier than the rest-frame FIR emission. Barro et al. (2016), for example, find that the 345-GHz (ALMA Band 7) dust continuum emission of a $M_{*}=10^{10.9}M_{\odot},\,\,\rm{SFR}=500M_{\odot}yr^{-1},\,\,z=2.45$ galaxy has a half-light radius which is half that of the rest-frame optical emission probed by HST. Hodge et al. (2016) imaged 16 $z\sim 2.5$ similarly massive, highly star-forming, luminous sub-mm galaxies at $0.16^{\prime\prime}$ in the same ALMA band. Many of these galaxies display clumpy structures in HST’s $H_{160}$ and $I_{814}$ bands, but their dust emission appears substantially smoother and more compact. The 870-µm  radii obtained are small (median $1.8\pm 0.2\,\rm{kpc}$), with no convincing evidence for clumpy dust emission at the ALMA resolution probed.

Molecular gas reservoirs have also been mapped for a handful of high-redshift, FIR-luminous galaxies. High-spatial-resolution studies show that molecular gas is compact, though it tends to be slightly more extended than the dust continuum emission. Tadaki et al. (2017) imaged two $z=2.5$ galaxies at 345 GHz, obtaining 870-µm  dust-emission radii of $1.2\pm 0.1\,\rm{kpc}$ and $1.3\pm 0.1\,\rm{kpc}$, around half the size of the CO(J=3-2) emission. Consistent results have been found by other studies (e.g. Strandet, 2017; Calistro Rivera et al., 2018; Tadaki et al., 2018), though sample sizes remain small.

A physical understanding of these differences in spatial extent of emission in different wavebands is critical in order to make the best use of the unmatched sensitivity and high spatial resolution of ALMA. This is currently difficult observationally, as only small samples of galaxies have been resolved at high resolution with multiple instruments, and such samples are often biased towards either the least dust-obscured systems (if selected in the UV) or the most compact, FIR-bright systems (if selected by FIR surface brightness). Interpretations are further complicated by uncertainty in what FIR/sub-mm fluxes actually probe. Frequently used relations between observed FIR/sub-mm fluxes and SFR (e.g. Kennicutt & Evans, 2012) do not fold in the shape of the full dust SED, which should reflect not only SFR but also dust mass, dust temperature, and the geometry of the source (e.g. Misselt et al., 2001; Hayward et al., 2011, 2012; Lanz et al., 2014; Safarzadeh et al., 2016; Kirkpatrick et al., 2017).

The primary aim of this paper is to understand which physical properties of high-redshift galaxies are probed by their dust continuum emission, in a spatially-resolved manner. We address this using simulations to study individual galaxies at high spatial resolution across the rest-frame UV-FIR spectrum. We select a sample of the most FIR-bright galaxies from the Feedback In Realistic Environments 2 (FIRE-2) simulations (Hopkins et al., 2018b)111http://fire.northwestern.edu presented in Anglés-Alcázar et al. (2017b). We perform three-dimensional continuum radiative transfer on selected galaxy snapshots from these simulations, modelling the effects of dust attenuation and re-emission to predict the spatially-resolved multi-wavelength emission of these high-redshift FIR/sub-mm-bright simulated galaxies. We then compare the spatial extent of the dust continuum emission to various intrinsic physical properties of our simulated galaxies.

The structure of this paper is as follows. In Section 2, we discuss our method for selecting sub-mm-bright snapshots from the FIRE-2 simulations. We describe the radiative transfer modelling used to post-process these galaxy snapshots and present their simulated rest-frame UV-FIR spectral energy distributions in Section 3. We present predictions for the spatial extent of dust continuum emission in Section 4 and also compare our predictions to observational results. We also study which physical properties of high-redshift, dusty star-forming galaxies are best probed by the spatial extent of the dust continuum emission. We summarise our conclusions in Section 5.

2 A sample of simulated high-redshift galaxies

2.1 Galaxies in the FIRE-2 simulations

The Feedback In Realistic Environments (FIRE) project (Hopkins et al., 2014, 2018b) is a set of state-of-the-art hydrodynamical cosmological zoom-in simulations that explore the role of stellar feedback in galaxy formation and evolution. Stellar feedback must play an important role in galaxy formation. Without it, the galactic ISM would collapse on dynamical time-scales, leading to gravitational collapse, fragmentation and accelerated star-formation. While galaxies simulated without stellar feedback thus rapidly convert all available gas into stars (e.g. Hopkins et al. 2011; see the review by Somerville & Davé 2015), the tight locus of observed galaxies on the Kennicutt-Schmidt relation (Kennicutt, 1998) implies that gas consumption timescales in real galaxies are long. Furthermore, both galaxy stellar mass functions (e.g. Ilbert et al., 2013; Muzzin et al., 2013; Davidzon et al., 2017) and the stellar mass-halo mass relation (e.g. Moster et al. 2010; Behroozi et al. 2013; Cochrane et al. 2018) imply that galaxies convert only a small fraction of the universal baryon fraction into stars. Galactic outflows are therefore needed to regulate the mass of galaxies over time (Kereš et al., 2009; Anglés-Alcázar et al., 2014; Muratov et al., 2015; Anglés-Alcázar et al., 2017a). These outflows are likely responsible for the observed enrichment of the circumgalactic medium (CGM) and intergalactic medium (IGM), and are also observed directly (e.g. Weiner et al. 2009; Steidel et al. 2010; see the review by Rupke 2018). These pieces of observational evidence demand that stellar feedback must be at play (Schaye et al., 2003; Oppenheimer & Davé, 2006; Faucher-Giguère et al., 2015; Hafen et al., 2017). Stellar feedback is believed to be particularly important in low-stellar-mass galaxies, below the peak of the stellar mass-halo mass relation.

Various stellar feedback processes are thought to contribute. These processes include supernovae, protostellar jets, photo-heating, stellar mass loss from O- and AGB-stars and radiation pressure (see Dale, 2015, for a review). Importantly, these processes are believed to act non-linearly, and so modelling the stellar processes of even a single galaxy becomes a complex computational challenge. Only recently have cosmological zoom simulations achieved sufficient resolution to model these feedback processes directly.

The FIRE project employs two main techniques to explicitly model multi-channel stellar feedback. Firstly, the FIRE simulations resolve the formation of giant molecular clouds (GMCs), and star-formation takes place only in self-gravitating (according to the Hopkins et al. 2013 criterion), self-shielding molecular gas (following Krumholz & Gnedin 2011) at high densities ($n_{H}>1000\,\rm{cm}^{-3}$ in the simulations we use here). Secondly, FIRE includes models for both energy and momentum return from the main stellar feedback processes, using directly the predictions of stellar population synthesis models without the extensive parameter tuning employed in other simulations. Specifically, once a star particle forms, the simulations explicitly follow several different stellar feedback mechanisms as detailed in Hopkins et al. 2018a, including (1) local and long-range momentum flux from radiation pressure (in both the initial UV/optical single-scattering regime and re-radiated light in the IR); (2) energy, momentum, mass and metal injection from SNe (Types Ia and II) and stellar mass loss (both OB and AGB); and (3) photo-ionization and photo-electric heating. Every star particle is treated as a single stellar population with known mass, age, and metallicity, and then all feedback event rates, luminosities and energies, mass-loss rates, and all other quantities are tabulated directly from stellar evolution models (starburst99; Leitherer et al. 1999), assuming a Kroupa (2001) initial mass function. The FIRE simulations succeed in generating galactic winds self-consistently (Muratov et al., 2015; Anglés-Alcázar et al., 2017a), without relying on sub-grid hydrodynamic decoupling or delayed cooling techniques, and broadly reproducing many observed galaxy properties, including stellar masses, star-formation histories and the ‘main sequence’ of star-forming galaxies (see Hopkins et al. 2014; Sparre et al. 2017), metallicities and abundance ratios (Ma et al., 2016; van de Voort et al., 2015), as well as morphologies and kinematics of both thin and thick disks (Ma et al., 2017).

For this paper, we study the central galaxies of four massive haloes originally selected and simulated by Feldmann et al. (2016); Feldmann et al. (2017) with the original FIRE model (Hopkins et al., 2014) as part of the MassiveFIRE suite. The haloes were selected to have dark matter halo masses of $M_{\rm{halo}}\sim 10^{12.5}M_{\odot}$ at $z=2$. The central galaxies of these haloes have stellar masses of $7\times 10^{10}-3\times 10^{11}M_{\odot}$ at $z=2$ (detailed in Table 1), with a variety of formation histories; see Feldmann et al. (2017) for details. The simulations in the present study are drawn from Anglés-Alcázar et al. (2017b), who re-simulated these four massive haloes down to $z=1$ with the upgraded FIRE-2 physics model (Hopkins et al., 2018b) and with a novel on-the-fly treatment for the seeding and growth of supermassive black holes (SMBHs). Note that these simulations do not include feedback from the accreting SMBHs. Compared to FIRE, FIRE-2 simulations are run with a new, more accurate hydrodynamics solver (a mesh-free Godunov solver implemented in the gizmo222http://www.tapir.caltech.edu/~phopkins/Site/GIZMO.html code; Gaburov & Nitadori 2011; Hopkins 2015). They also feature improved treatments of cooling and recombination rates, gravitational softening and numerical feedback coupling, and they adopt a higher density threshold for star-formation (Hopkins et al., 2018a). The mass resolution of our simulations is $3.3\times{}10^{4}$ $M_{\odot}$ for gas and star particles and $1.7\times{}10^{5}$ $M_{\odot}$ for dark matter particles.

2.2 Selection of sub-mm-bright galaxy snapshots at $z>2$

We wish to simulate galaxies that are representative of those typically observed with ALMA at high spatial resolution at high redshifts, which implies that we should select those that are likely to have high 850-µm  flux densities ($f_{850}\ga 1\,{\rm mJy}$). Performing radiative transfer on each of the $\sim 300$ redshift snapshots available to predict sub-mm fluxes and then selecting the brightest would be unnecessarily computationally intensive. The first step in our analysis is therefore to select redshift snapshots for each of the four galaxies for which we expect the sub-mm flux to be particularly bright, using simply the SFR and dust mass at each redshift. We adopt the following equation, derived from fits to the observed-frame 850-µm flux densities of simulated galaxies computed via dust radiative transfer and presented in Hayward et al. (2013):

[TABLE]

where we estimate the dust mass for the present purposes using $M_{\rm{dust}}=0.01M_{\rm{gas}}$, where $M_{\rm{gas}}$ is the total gas mass within $0.1R_{\rm{vir}}$. The agreement between these predictions and the 850-µm flux densities that result from the radiative transfer modelling is correct to within a factor of $\sim 2$ for fluxes above $\sim 1\rm{mJy}$. Future work will involve running radiative transfer over a larger number of snapshots and re-fitting this formula. We select the $\sim 20$ snapshots with the highest predicted $f_{850}$ for each simulated halo (named A1, A2, A4 & A8)333Our analysis focuses on the central galaxies of each halo. Throughout, we will refer to the central galaxy of e.g. halo A1 simply as ‘galaxy A1’., within the redshift range $1<z<5$. A couple of snapshots were excluded from the analysis after performing radiative transfer due to poorly defined centres and extremely diffuse dust emission. We do not attempt to produce a complete sample of sub-mm bright galaxies from the FIRE-2 simulations; instead, our selection is sufficient to yield a sample of sub-mm bright snapshots for which we can perform radiative transfer and study multi-wavelength properties.

3 Radiative transfer methods and results

3.1 Performing radiative transfer with skirt

Modelling dust and its emission in galaxies is a difficult computational problem (see Steinacker et al., 2013, for a comprehensive review). The process of radiative transfer is non-local in space (photons can propagate long distances before interacting with dust), and it is coupled in terms of both direction and wavelength. The distribution of dust in (both real and simulated) galaxies is far from a simple screen; instead, it is necessary to model the three-dimensional distribution of sources of radiation (stars and AGN) and dust.

In this work, we make use of the Stellar Kinematics Including Radiative Transfer (skirt)444http://www.skirt.ugent.be Monte Carlo radiative transfer code (Baes et al., 2011; Camps & Baes, 2015), which has also been used to model dust attenuation and emission in the EAGLE simulations (see Trayford et al., 2017; Camps et al., 2018; McAlpine et al., 2019). Monte Carlo radiative transfer codes like skirt treat the radiation field from stars (and sometimes AGN) as a flow of photons through the dusty medium of a galaxy to compute the effects of dust absorption, scattering, and re-emission of the absorbed light, including dust self-absorption. We are able to perform the radiative transfer in post-processing because the light-crossing time is short compared to dynamical times on resolved scales (such that the dust and stellar geometry does not change appreciably over a light-crossing time).

We extract gas and star particles from the FIRE-2 simulations at each of our chosen snapshots. The coordinate system is rotated to align with the angular momentum vector of the gas particles within $0.1R_{\rm{vir}}$ prior to input to skirt, so that for disk galaxies, a viewing angle of 0 degrees corresponds to face-on. For gas particles with temperature $<10^{6}\,\rm{K}$, we compute dust masses using the metallicity of the gas particles and a dust-to-metals mass ratio of $0.4$ (Dwek, 1998; James et al., 2002). We assume that dust is destroyed in gas particles with temperature $>10^{6}\,{\rm K}$ (Draine & Salpeter, 1979; Tielens et al., 1994). We use a Weingartner & Draine (2001) Milky Way dust prescription to model a mixture of graphite, silicate and PAH grains. Star particles are assigned Bruzual & Charlot (2003) SEDs according to their ages and metallicities (note that our results are unchanged if we instead use the starburst99 SED templates presented by Leitherer et al. 1999). We use an octree dust grid, in which cell sizes are adjusted according to the dust density distribution. We impose the condition that no dust cell may contain more than $0.0001\%$ of the total dust mass of the galaxy, which yields excellent convergence of the integrated SED.

We also specify a number of properties relating to the SED output. We define a wavelength grid with $\sim 100$ discrete wavelengths, spaced uniformly in log(wavelength) between rest-frame UV and FIR wavelengths and including each of the ALMA bands. To model the flux that would be received by an observer on Earth, we place seven detectors at $z=0$, uniformly spaced at different inclinations with respect to the disk plane of the galaxy. These detectors have pixel sizes corresponding to a proper length of $25\,\rm{pc}$, and the field of view is set to 10% of the virial radius for each galaxy snapshot studied. Note that the pixel scale of the images should not be confused with the resolution of the simulation; this varies across the galaxy, and is higher than $25\,\rm{pc}$ in denser regions but lower in general ISM gas.

The focus on resolved dust emission in the present work is complementary to other radiative transfer studies of FIRE galaxies with skirt. Liang et al. (2018) analyzed how well galaxy-integrated, single-band dust continuum emission traces molecular gas. The implications of defining galaxy-integrated dust temperatures in different ways, for instance with regard to various scaling relations, are studied in Liang et al. (2019). Finally, Ma et al. (2019) make predictions for the dust extinction and emission properties of $z\ga 5$ galaxies using a sample of 34 FIRE-2 haloes, including those first presented in Ma et al. (2018).

3.2 Morphology across the wavelength spectrum

Our radiative transfer modelling enables us to track the emission from each of the galaxies at multiple epochs. We compute images of each of the galaxy snapshots at every wavelength on our grid, spanning the rest-frame far-UV to far-IR. We show an example of the high quality of our maps of the sub-mm flux in Figure 1, to illustrate the wide range of morphologies displayed by a single galaxy evolving in the redshift range $z=2.41-1.45$. We also find that the same galaxy can look vastly different in the different wavebands. We illustrate this qualitatively with five representative wavelengths for the central galaxy of halo A1 (hereafter ‘galaxy A1’) at $z=4.38$ in Figure 2 and galaxy A2 at $z=2.95$ in Figure 3. We also show for comparison the spatial distributions of four key galaxy properties: total stellar mass, cool, dense gas mass ($T<300\,{\rm K},n_{H}>10\,\rm{particles}/\rm{cm}^{-3}$), dust mass, and recently formed stars ($\rm{age}<100\,\rm{Myr}$).

Galaxy A1 at $z=4.38$ has a skirt-predicted observed-frame 850-µm  flux density of $0.79\,\rm{mJy}$. At this snapshot, the galaxy exhibits very clumpy FUV emission (rest-frame $1500\mbox{\AA}$), spanning a few kpc. This emission is aligned with the cool gas and exhibits peaks where this gas is densest. The optical emission (rest-frame $4000\mbox{\AA}$) is also clumpy. Longwards of 24-µm  the emission becomes more regular and centrally concentrated, resembling the total and recently formed stellar mass distributions more closely.

Galaxy A2 at $z=2.95$ is brighter, with a skirt-predicted observed-frame 850-µm  flux density of $1.45\,\rm{mJy}$. It is also substantially more extended in all wavebands. The rest-frame $1500\mbox{\AA}$ and $7000\mbox{\AA}$ emission is again clumpy, bearing little resemblance to the ordered bulge plus spiral structure that is clear from the SFR and stellar, gas and dust mass maps. Interestingly, the peak of the short-wavelength emission occurs in a region to the right of centre that is largely free of dust (see Figure 4). It appears that the clumpiness of the FUV emission is driven by the structure of the dust, with FUV emission tracing holes in the dust. This is just one example, and such offsets are common in the sample, as in observations (e.g. Chen et al., 2017; Calistro Rivera et al., 2018, Cochrane et al. in prep). Future work will explore these offsets, as well as the resolved multi-wavelength emission on a pixel-by-pixel basis.

4 Quantifying the spatial extent of dust continuum emission

4.1 The dependence of size on FIR wavelength

In this section, we quantify the sizes of sub-mm-bright FIRE-2 galaxies as a function of wavelength. An example of an observed-frame SED predicted by skirt, overplotted with the wavelengths corresponding to eight ALMA bands, is shown in Figure 5 (left-hand panel). We extract the skirt image at each of these wavelengths and derive an ‘effective radius’ for the predicted emission. This is the radius that contains half of the total flux (calculated within a $\sim 12\rm{kpc}$ radius), derived using a circular aperture centred on the flux-weighted centre of the emission in each band. In each case, an error bar is derived from the standard deviation of the effective radius measurements at seven different galaxy/detector inclinations.

In the right-hand panel of Figure 5, we show an example of our results. The effective radius of the emission varies with wavelength, with longer wavelength FIR emission spanning a greater spatial extent. We also overplot the effective radii of several key physical quantities of the galaxy: stellar mass, dust mass and cool, dense gas mass. We define cool, dense gas using the criteria $T<300\,{\rm K},n_{H}>10\,\rm{cm}^{-3}$. This has been shown to be a reasonable proxy for molecular gas by Orr et al. (2018). In this example, and in general, the cool gas mass and the dust mass are more spatially extended than the stellar mass. The sizes of the dust continuum emission tend to be between those of the dust mass and the stellar mass. Shorter-wavelength FIR emission, corresponding to hotter dust, tends to be more compact, whereas longer-wavelength FIR emission better traces the extended cold gas and dust.

The 870-µm  (345-GHz) observed-frame emission probed by ALMA Band 7 is frequently used to study the dust continuum emission of high-redshift galaxies (e.g. Barro et al., 2016; Chen et al., 2017; Simpson et al., 2017). For simplicity, from here on, we consider only the 870-µm  flux in our discussion of the spatial extent of the dust continuum emission. However, we note that we consistently find that emission at longer rest-frame FIR wavelengths is more extended, as shown in Figure 5.

4.2 Agreement with compact dust continuum observations

In Figure 6, we compare the sizes of both our predicted 870-µm  dust emission and the cool, dense gas to the following ALMA observations of galaxies at similar redshifts: 860-µm emission and $12\rm{CO}(J=4-3)$ line flux for a $\rm{SFR}>1000M_{\odot}/yr$, $M_{*}\sim 10^{11}M_{\odot}$ galaxy at $z=4.3$, from Tadaki et al. (2018) (grey star); 870-µm  emission and $12\rm{CO}(J=3-2)$ line flux for two massive ($M_{*}\sim 10^{11}M_{\odot}$) galaxies at $z=2.5$, from Tadaki et al. (2017) (black stars); the effective radii for stacked ALMA maps of 870-µm  and $12\rm{CO}(J=3-2)$ emission for 16 ALESS galaxies at $z=2.5\pm 0.2$, from Calistro Rivera et al. (2018) (black circle); 870-µm  emission and $12\rm{CO}(J=3-2)$ line flux for a $\rm{SFR}\sim 500M_{\odot}/yr$, $M_{*}\sim 2\times 10^{11}M_{\odot}$ galaxy at $z=2.12$, from Chen et al. (2017) (black square); FIR emission (from APEX/LABOCA and Herschel) and $12\rm{CO}(J=3-2)$ at $z=6.9$, from Strandet (2017) (black triangle); stacked 870-µm  and [CII] emission for $z\sim 4.5$ galaxies, from Cooke et al. (2018) (purple star); $1.1\rm{mm}$ (from NOEMA) and $\rm{H}\alpha$ emission, mapped for a $M_{*}\sim 7\times 10^{10}M_{\odot}$ at $z=1.25$, from Nelson et al. (2019) (purple circle). Our derived 870-µm  effective radii are $\sim 0.5-4\rm{kpc}$. This is in excellent agreement with observations at a range of redshifts (e.g. Ikarashi et al., 2015; Iono et al., 2016; Simpson et al., 2015; Hodge et al., 2016). In line with the observations, the snapshots of each of the four FIRE-2 galaxies tend to lie below the 1-1 line (dashed black), i.e. the dust continuum emission is more compact than the ‘molecular’ gas. Our predictions agree fairly well with the observational results for the small number of high-redshift galaxies that have been observed in both dust continuum and $\rm{CO}$ at high spatial resolution. However, the few observations that do exist tend to have slightly more compact dust emission at fixed molecular gas effective radius than our predictions. This could be due to the definition of ‘molecular’ gas within the simulation (different temperature and density cuts yield slightly different ‘molecular’ gas sizes, and hence slopes; ‘very dense gas’ is much more compact, as shown in Figure 7), or due to the selection of observational samples.

4.3 How does the dust continuum emission trace the physical properties of galaxies?

The dust continuum emission is frequently used as an indicator of both SFR and dust mass. However, until now, the spatial extent of these physical and observable quantities has not been studied consistently with simulated galaxies. Motivated by the clear differences in effective radii of the dust continuum emission, stellar mass, dust mass, and gas mass found for individual FIRE-2 galaxies (see the right-hand panel of Figure 5), we set out to identify which physical properties are best reflected by the spatial extent of the dust continuum emission.

In Figure 7, we plot the effective radius of the simulated observed-frame 870-µm  dust continuum emission versus those of eight different physical quantities derived directly from the FIRE-2 simulations. The 870-µm  emission is more compact than the total gas component, but only slightly more compact than the cool, dense gas. The dust emission is also more compact than the dust mass distribution for all FIRE-2 snapshots. We also consider how the spatial extent of dust emission correlates with that of the total stellar mass, recently formed stellar mass (within 10 Myr and 100 Myr) and instantaneous SFR. We find that the emission from dust is more extended than all of these stellar quantities and than the densest gas ($n_{H}>100\,\rm{particles}/\rm{cm}^{3}$). The extent of the dust emission appears to correlate more tightly with the extent of the most recently formed stars than the total stellar mass. This reflects the role of young stars as the primary heating source for the dust grains that reprocess their short-wavelength light, as will be discussed in detail in the next two subsections. Note also that stellar mass appears to assemble inside-out, with the more recently formed stars spanning a larger half-mass radius than the total stellar component.

4.4 The role of star-formation in determining the extent of the dust continuum emission

The observed-frame submm emission depends on the effective dust temperature, which is sensitive to the luminosity absorbed by dust, the dust mass, and possibly the geometry of the system (e.g. Misselt et al., 2001; Hayward et al., 2011; Safarzadeh et al., 2016). Consequently, we expect that the spatial extent of the dust continuum will depend on those of the young stars and dust. To quantify this dependence, we use a Principal Component Analysis (PCA). PCA is a statistical technique used to describe the variance within a dataset. Variables are converted into a set of uncorrelated, orthogonal principal components. The first component reveals the direction of maximum variance, which describes the key correlation between variables within the dataset. Successive components account for less of the variance of the population. Some latter components may be dominated by noise, leaving the data decomposed into fewer dimensions. PCA has been used in a number of recent galaxy evolution studies, including in studies of the mass-metallicity relation (Lagos et al., 2016; Bothwell et al., 2016; Hashimoto et al., 2018) and to study the quenching of galaxies within the EAGLE simulations (Cochrane & Best, 2018).

Here, we use PCA to study the relationship between the effective radii of the dust component, SFR and dust continuum emission. We use the PCA python tool scikit.learn to perform the PCA analysis. Each variable is normalized to its mean and scaled to unit variance before performing the PCA. For each snapshot studied, we construct a vector of the form [$r_{\rm{eff,\,\,dust\,\,emission}}$, $r_{\rm{eff,\,\,dust}}$, $r_{\rm{eff,\,\,stars\,\,<100Myr}}$] and supply all of these vectors to the PCA.

We find that the first principal component, [0.62, 0.56, 0.55], almost entirely describes the dataset, accounting for $83\%$ of the variance of the sample. This primary correlation indicates that all three variables – the effective radii of dust emission, dust mass and recently formed stars – are positively correlated. The effective radius of the dust emission exhibits a strong correlation with the effective radius of the dust mass itself and with the effective radius of the recently formed stars. This result confirms our physical intuition that the sub-mm flux density (and thus the effective radius of the dust continuum emission) is sensitive to both the SFR and dust mass (Hayward et al., 2011; Safarzadeh et al., 2016). The second principal component, [0.03, 0.69,-0.73], accounts for $15\%$ of the scatter, essentially all of the rest. This component represents the scatter in the relation between the effective radii of the dust mass and the recently formed stars. Note that there is, therefore, little scatter in the strong correlations between the effective radius of the dust emission and the dust mass and between the effective radius of the dust emission and the recently formed stars.

4.5 The physical drivers of compact dust emission: dust temperature gradients

In Section 4.4, we found that the effective radius of the dust continuum emission is correlated with both the effective radius of the dust mass itself, and the effective radius of the recently formed stars. In this final section, we investigate the physical drivers of these correlations, in particular the role of dust heating. In Figure 8, we plot the effective radius of the 870-µm  dust continuum emission against that of the recently formed stars, normalising both quantities by the effective radius of the dust mass. If the dust continuum emission traced the dust mass directly, all points would lie at $r_{e,\rm{dust\,emission}}/r_{e,\rm{dust}}=1$. While $r_{e,\rm{dust\,emission}}/r_{e,\rm{dust}}$ is close to 1 for most snapshots, we find a clear slope in this relation: galaxies with low $r_{e,\rm{dust\,emission}}/r_{e,\rm{dust}}$ (which have particularly compact dust emission with respect to the total dust mass distribution) also have recent star-formation which is much more compact than the dust. This reflects the dual role of the dust mass and star-formation in determining the spatial extent of the dust continuum emission.

We have identified that the spatial scales of dust and recently formed stars play roles in determining the extent of the dust continuum emission, with very compact dust emission appearing to correlate with compact recent star-formation. Now, we use the dust mass-weighted temperatures of the grid cells output by skirt to explore the physical drivers of particularly compact dust emission. The temperature of the dust in our simulated galaxies varies significantly as a function of distance from the galaxy’s centre. Typical dust temperatures are $\sim 20-40\rm{K}$ outside the central kpc, and $\sim 40-70\rm{K}$ inside the central kpc. Here, we study the relationship between dust temperature gradients and the spatial extent of dust emission. For each snapshot, we construct thin concentric shells about the halo centre, out to the effective radius calculated for the dust mass. For each shell, we calculate the mean dust temperature, $T_{\rm{dust}}$, and mean dust particle position, $r_{\rm{dust}}$. We then fit the following relation: $\log_{10}T_{\rm{dust}}=\nabla T_{\rm{dust}}\log_{10}r_{\rm{dust}}+C$, where the units of $T_{\rm{dust}}$ are K and the units of $r_{\rm{dust}}$ are kpc. This yields a dust temperature gradient (strictly a power-law index), $\nabla T_{\rm{dust}}$, for each galaxy snapshot. The snapshots in Figure 8(a) are colour-coded by this gradient. It is clear that those galaxies with compact dust emission and recent star-formation also have the steepest dust temperature gradients. The other panels of Figure 8 show ratios of the effective radii of the dust mass, dust continuum emission and recently formed stars. For gradients flatter than $\nabla T_{\rm{dust}}\sim-0.2$, the dust continuum emission broadly traces the extent of both recently formed stars and dust mass. For steeper gradients, below $\nabla T_{\rm{dust}}\sim-0.2$, the recently formed stars are substantially more compact than the dust mass. This causes the ratio of absorbed luminosity to dust mass to be higher in the galaxy centre than at larger radii, leading to the steep dust temperature gradients. At these low values of $\nabla T_{\rm{dust}}$, the spatial extent of the dust continuum emission is not a good approximation of the distribution of dust mass or recently formed stars.

These results confirm that the spatial extent of the dust continuum emission is sensitive to the scales of recent star-formation because of the effects of dust heating. Thus, inferring spatial scales of star-formation or dust mass from observations of dust continuum emission is non-trivial. This may be done more robustly using measurements of dust temperature gradients within galaxies. Studies of local galaxies show that dust temperature distributions can be derived from spatially-resolved imaging in multiple bands (e.g. Galametz et al., 2012).

5 Conclusions

We have performed radiative transfer modelling on a subset of rest-frame FIR-bright redshift snapshots of four massive galaxies drawn from the FIRE-2 simulations presented in Anglés-Alcázar et al. (2017b). The simulated galaxies have stellar masses $7\times 10^{10}-3\times 10^{11}M_{\odot}$, and reside in dark matter haloes of mass $M_{\rm{halo}}\sim 10^{12.5}M_{\odot}$ at $z=2$. Our modelling yields full rest-frame FUV-FIR SEDs as well as maps of the emission in $\sim 100$ wavebands, resolved to $25\rm{pc}$ scales. We find clear differences between the morphologies of the same galaxies in the different wavebands, with shorter-wavelength emission (UV/optical) frequently appearing clumpy and extended. We find that this is due to the emitted short-wavelength light tracing ‘holes’ in the dust distribution. At rest-frame FIR wavelengths, the galaxies tend to have more regular morphologies.

The primary focus of our study is the spatial extent of the dust continuum emission. For the snapshots we study, the FIRE-2 galaxies have observed-frame 870-µm  fluxes of up to $\sim 2\,\rm{mJy}$ at $z\sim 1-5$. These simulated FIR-bright galaxies exhibit very compact dust continuum emission, with effective radii of $\sim 0.5-4\rm{kpc}$, in line with existing observations of star-forming galaxies at these redshifts. At rest-frame FIR wavelengths (which we can probe with ALMA for high-redshift galaxies), longer-wavelength emission tends to be more extended because the dust tends to be hotter in the central regions of galaxies, and the shorter-wavelength emission is dominated by hotter grains.

We also compare the spatial extent of the 870-µm  dust continuum emission to that of several key intrinsic physical quantities, including the dust, gas and stellar components. In both simulated and observed galaxies, the dust continuum emission is more compact than the cool, dense gas and the dust, but more extended than the stellar component. Extracting only recently formed ($\rm{age}<100\rm{Myr}$) stars from the simulations enables us to study the role of recent star-formation in determining the spatial extent of the dust emission. We find that in some snapshots, the simulated galaxies exhibit extremely compact dust emission ($\sim 0.5\rm{kpc}$). This emission appears to be driven by particularly compact recent star-formation. Overall, the spatial extent of the dust continuum emission correlates with those of both the dust mass and the recently formed stars. Physically, this makes sense: the dust that emits strongly in the FIR is predominantly heated by young massive stars (as these simulations do not contain AGN), and when the recent star-formation is particularly compact, the central regions of the galaxy have steeper dust temperature gradients and consequently more compact emission. In such systems, constraints on the dust temperature gradient are necessary to infer the spatial extents of the young stars and dust reliably. This work thus motivates multiband ALMA observations to constrain the dust temperature gradients of observed galaxies.

Acknowledgements

This work was initiated as a project for the Kavli Summer Program in Astrophysics held at the Center for Computational Astrophysics of the Flatiron Institute in 2018. The program was co-funded by the Kavli Foundation and the Simons Foundation. We thank them for their generous support. RKC acknowledges funding from STFC via a studentship and thanks Philip Best for helpful discussions and comments. DAA acknowledges support from a Flatiron Fellowship. The Flatiron Institute is supported by the Simons Foundation. DK was supported by NSF grant AST-1715101 and the Cottrell Scholar Award from the Research Corporation for Science Advancement. RF acknowledges financial support from the Swiss National Science Foundation (grant number 157591). CAFG was supported by NSF through grants AST-1517491, AST-1715216, and CAREER award AST-1652522; by NASA through grant 17-ATP17-0067; by CXO through grant TM7-18007; by STScI through grant HST-AR-14562.001; and by a Cottrell Scholar Award from the Research Corporation for Science Advancement. The simulations were run using XSEDE, supported by NSF grant ACI-1053575, and Northwestern University’s compute cluster ‘Quest’.