The IRX-Beta Dust Attenuation Relation in Cosmological Galaxy Formation Simulations

TL;DR

This study uses high-resolution cosmological simulations with dust radiative transfer to analyze the IRX-beta relation in galaxies, revealing key factors affecting scatter and deviations, and providing a fitting relation to unify diverse galaxy populations.

Contribution

It offers a detailed physical understanding of the IRX-beta relation's scatter and deviations, and introduces a fitting relation to unify galaxy populations across redshifts.

Findings

Galaxies with young stars and Milky Way-like dust follow the standard IRX-beta relation.

Scatter is mainly caused by stellar age, geometry, and dust extinction curve variations.

A new fitting relation reduces scatter and explains high-redshift galaxy behaviors.

Abstract

We utilise a series of high-resolution cosmological zoom simulations of galaxy formation to investigate the relationship between the ultraviolet (UV) slope, beta, and the ratio of the infrared luminosity to UV luminosity (IRX) in the spectral energy distributions (SEDs) of galaxies. We employ dust radiative transfer calculations in which the SEDs of the stars in galaxies propagate through the dusty interstellar medium. Our main goals are to understand the origin of, and scatter in the IRX-beta relation; to assess the efficacy of simplified stellar population synthesis screen models in capturing the essential physics in the IRX-beta relation; and to understand systematic deviations from the canonical local IRX-beta relations in particular populations of high-redshift galaxies. Our main results follow. Galaxies that have young stellar populations with relatively cospatial UV and IR emitting regions and a Milky Way-like extinction curve fall on or near the standard Meurer relation. This behaviour is well captured by simplified screen models. Scatter in the IRX-beta relation is dominated by three major effects: (i) older stellar populations drive galaxies below the relations defined for local starbursts due to a reddening of their intrinsic UV SEDs; (ii) complex geometries in high-z heavily star forming galaxies drive galaxies toward blue UV slopes owing to optically thin UV sightlines; (iii) shallow extinction curves drive galaxies downward in the IRX-beta plane due to lowered NUV/FUV extinction ratios. We use these features of the UV slopes of galaxies to derive a fitting relation that reasonably collapses the scatter back toward the canonical local relation. Finally, we use these results to develop an understanding for the location of two particularly enigmatic populations of galaxies in the IRX-beta plane: z~2-4 dusty star forming galaxies, and z>5 star forming galaxies.

Figures

Tables

Table 1. Table 1: Summary of Model Galaxies. M ∗ subscript 𝑀 M_{*} and M halo subscript 𝑀 halo M_{\rm halo} refer to the stellar and halo mass at z = 2 𝑧 2 z=2 .

Model Name	$M_{\ast }$ (z=zfinal)	$M_{\mathrm{halo}}$ (z=zfinal)	Final Redshift	Notes
	$M_{\odot }$	$M_{\odot }$
mz0	$1.7\times {10}^{11}$	$2.95\times {10}^{13}$	2
mz0q	$6\times {10}^{10}$	$8.4\times {10}^{12}$	2.2	Quenching model included
mz5	$1.5\times {10}^{11}$	$7.65\times {10}^{12}$	2
mz10	$3.3\times {10}^{10}$	$4.4\times {10}^{12}$	2
mz45	$1.2\times {10}^{10}$	$9.95\times {10}^{11}$	2
z0mz401	$2.6\times {10}^{10}$	$2.05\times {10}^{12}$	0

Equations

$$\frac{d{M}_{\ast }}{dt}={ϵ}_{\ast }\frac{\rho f_{H2}}{t_{\mathrm{dyn}}}$$

$$\eta=3.55\left(\frac{M_{*}}{10^{10}\mbox{M${}_{\sun}$ }}\right)^{-0.351}$$

$$\mathrm{IRX}\equiv \frac{L_{\mathrm{IR}}}{L_{\mathrm{UV}}}$$

$$\mathrm{IR}{\mathrm{X}}_{\mathrm{ref}}=1.68\times \left[10^{0.4[(3.36)+2.04{\beta }_{\mathrm{ref}}]}-1\right]$$

$$\Delta {\beta }_{\mathrm{ref}}=-0.72-0.27\mathrm{lo}{\mathrm{g}}_{10}\left(\frac{\mathrm{SFR}}{M_{\odot }{\mathrm{yr}}^{-1}}\right)+\exp \left[1.47\times \mathrm{lo}{\mathrm{g}}_{10}\left(\frac{t_{\mathrm{age}}}{\mathrm{Gyr}}\right)\right]$$

$${\tau }_{1600}=1.47+1.98\times ({\mathrm{IRX}}^{0.23})$$

Full text

††pagerange: The IRX-$\beta$ Dust Attenuation Relation in Cosmological Galaxy Formation Simulations–LABEL:lastpage††pubyear: 2015

1 Introduction

Measuring the evolution of the cosmic star formation rate density remains both a challenge and fundamental goal of galaxy evolution astrophysics (see Madau & Dickinson, 2014, for a recent review). Essentially all tracers of a galaxy’s star formation rate (SFR) measure the luminosity (either directly or reprocessed) from massive, short-lived stars. Because the ultraviolet (UV) continuum from actively star-forming galaxies is dominated by massive stars, direct measurements of UV starlight from galaxies (alongside an assumption of a stellar initial mass function) serves as one of the more straight forward methods for deriving a galaxy’s SFR (e.g. Kennicutt & Evans, 2012). This said, obscuring dust in the interstellar medium of galaxies preferentially absorbs UV radiation, and therefore complicates the interpretation of SFRs calibrated from UV measurements.

The observed relationship between the rest-frame UV continuum slope from galaxies ($\beta$, where $f_{\lambda}\propto\lambda^{\beta}$), and the ratio of the infrared luminosity to the UV luminosity (IRX) has been commonly employed as a tool for accounting for the dust obscuration of UV and optical radiation. In seminal work, Calzetti (1997) and Meurer et al. (1999) found that local starburst galaxies lie on a well-defined sequence in the IRX-$\beta$ plane. In principle, this notional relationship is quite powerful: even when infrared data is unavailable to constrain the true transfer of energy from UV photons to thermal dust continuum, an IRX-$\beta$ relationship can clarify the degree of dust obscuration from observations where only UV measurements are available.

Indeed, a large number of studies have employed this technique to infer the dust content of \mbox{z\sim }2-6 galaxies for a range of applications. Amongst many others, these include: constraining the cosmic star formation rate density (Bouwens et al., 2014, 2015), cosmic star formation rate functions (e.g. Bothwell et al., 2011; Smit et al., 2012, 2015), the evolution of specific star formation rates (e.g. González et al., 2010, 2014), the growth of galaxy dust content (Finkelstein et al., 2012), physical properties of individual galaxy populations (e.g. Stark et al., 2009; Ho et al., 2010; Williams et al., 2014), and SFR tracers themselves (e.g. Treyer et al., 2007; Kennicutt et al., 2009; Hao et al., 2011). The rationale for using the IRX-$\beta$ relation across a diverse range of observations has been, in part, bolstered by a number of observations that show the relation holding for samples of galaxies between $z\sim$ 2-4 (e.g. Reddy et al., 2008; Reddy et al., 2012; Pannella et al., 2009; Seibert et al., 2002; Heinis et al., 2013; To et al., 2014; Bourne et al., 2016).

At the same time, a number of observations have called into question the relationship between $\beta$ and the infrared excess. For example, at low-redshift, Goldader et al. (2002) and Howell et al. (2010) find that luminous infrared galaxies (LIRGs; LIR$>10^{11}$ L*☉) and ultraluminous infrared galaxies (ULIRGs; LIR$>10^{12}$ L☉*) are offset from the canonical IRX-$\beta$ relations, with either an excess of IRX at a fixed $\beta$, or bluer UV spectral slopes at a fixed IRX (i.e. they lie above the relationship).

This offset extends to high-redshift dusty galaxies as well. For example, while Bourne et al. (2016) found that lower luminosity galaxies at $z\sim$ 2 present similarly as the Meurer et al. (1999) sample in the IRX-$\beta$ plane, higher luminosity dusty star forming galaxies tend to exhibit bluer UV continuum slopes than the canonical relationship at a given IRX value. A similar trend is observed by Reddy et al. (2010), Penner et al. (2012), Oteo et al. (2013); Watson et al. (2015), and Casey et al. (2014b) in high redshift ($z\ga 2$) dusty star forming galaxies. Casey et al. (2014b) explicitly demonstrate that galaxies of increasing infrared luminosity at $z\ga 0.5$ have increasingly blue UV SEDs. We will return to this point later in this paper.

Offsets from the Meurer et al. (1999) local relation are not exclusive to dusty star forming galaxies. Metal-poor systems such as the Small and Large Magallenic Clouds (SMC and LMC respectively) have redder colours than more metal rich galaxies on the IRX-$\beta$ plane (e.g. Bell et al., 2002; Buat et al., 2005). Similarly, some high-$z$ Lyman Break Galaxies are more consistent with the SMC/LMC IRX-$\beta$ curve than the Meurer et al. (1999)-defined starburst locus (e.g. Reddy et al., 2010; Koprowski et al., 2016). Even some so-called “normal” galaxies at high-redshift (those not currently undergoing a starburst event, with $\sim$solar metallicity) have redder UV continuum slopes than the canonical Meurer relation, and lie in between the Meurer relation and SMC/LMC curves (e.g. Buat et al., 2005; Seibert et al., 2005; Boquien et al., 2009; Casey et al., 2014b; Pope et al., 2017).

With the advent of the Atacama Large Millimetre Array (ALMA), dust continuum detections (or sensitive upper limits) have now become routine at the highest redshifts (e.g. Dunlop et al., 2017; Rujopakarn et al., 2016). For example, Capak et al. (2015) and Bouwens et al. (2016) employed $0.8-1.1$ mm observations of $z\ga 4$ galaxies in order to infer their far infrared luminosities (this requires an assumed dust temperature; we will return to this point in § 5.2). From this data, Capak et al. and Bouwens et al. find that these high-redshift galaxies are systematically low in their dust content, lying below even the SMC/LMC IRX-$\beta$ curves.

In this paper, we aim to develop a physical model for the origin of and deviations from the IRX-$\beta$ relationship in galaxies. To do this, we will employ a series of high-resolution cosmological zoom galaxy formation simulations that are coupled with both stellar population synthesis and dust radiative transfer models. This will allow us to in effect “observe” the simulations, and therefore relate the physical properties of the model galaxies to their observed UV continuum and infrared luminosity properties.

We have three principle goals in this paper:

To understand the origin of and scatter in the IRX-$\beta$ relation 2. 2.

To asses the efficacy of simplified (i.e. screen) models in capturing the essential physics driving the IRX-$\beta$ relation 3. 3.

To utilise these models to understand systematic deviations from the canonical IRX-$\beta$ relations in particular populations of galaxies.

To this end, we organize the paper as follows. After describing our numerical methodology and simulations (§ 2), we deconstruct the IRX-$\beta$ relation, assessing the physical drivers of dispersion in § 3 and § 4. We then discuss individual galaxy populations that appear to have systematic deviations from the IRX-$\beta$ relation in § 5. Throughout these sections, we will also include comparisons to simplified models in which we place a dust screen in front of a simple stellar population (hereafter “screen models”). In § 6, we provide discussion, where we discuss a methodology to minimise the uncertainty resultant from systematic deviations from the fiducial IRX-$\beta$ relation, compare to other theoretical models, and assess the $\tau_{\rm UV}-\beta$ and $\tau_{\rm UV}$-IRX relationships. Finally, we summarise in § 7

2 Numerical Methods

2.1 Overview

Our overall aim is to simulate the rest-frame UV-millimetre wave SED of a realistic sample of model galaxies in order to investigate the IRX-$\beta$ relation in galaxies. To do this, we model the formation and evolution of galaxies utilising the cosmological “zoom” technique, wherein individual galaxies are tracked in a cosmological simulation at a higher resolution than the bulk of the cosmic volume. This technique, while computationally expensive, allows us to maintain a relatively high mass resolution while still following the realistic evolution of galaxies in their cosmic environment.

We couple these simulations with stellar population synthesis models to model the spectral emission properties of the stars, and dust radiative transfer simulations in order to calculate the dust attenuation the stellar radiation encounters as it leaves the galaxy. In the remainder of this section, we describe the details of these galaxy formation simulations and dust radiative transfer calculations, though note that the reader interested primarily in the main results can skip this section without loss of continuity.

2.2 Cosmological Zoom Galaxy Formation Simulations

We conduct all of our galaxy formation simulations with a modified version of the hydrodynamic code gizmo (Hopkins, 2014), which builds off of much of the code base of gadget-3 (Springel, 2005). We simulate a $50h^{-1}$ Mpc box from an initial redshift $z_{\rm init}=249$ down to $z_{\rm final}=0$, with initial conditions generated with music(Hahn & Abel, 2011). We employ a cosmology $\Omega_{\Lambda}=0.7$, $\Omega_{\rm b}=0.048$, $H_{0}=68$ km s*-1* Mpc*-1*, and $\sigma_{8}=0.82$. Our initial simulation includes dark matter only, and included $512^{3}$ particles, resulting in a dark matter mass resolution of $7.8\times 10^{8}h^{-1}M_{\odot}$.

We identify five arbitrary dark matter halos at redshift $z=z_{\rm final}$, where $z_{\rm final}$ is the final redshift of the zoom. We have two categories of simulations: massive galaxies intended to model a diverse range of physical properties at high redshift, as well as a lower mass galaxy intended to represent ’normal’ Milky Way like galaxies in the present epoch. The former category of simulations are selected at (and eventually run to) $z_{\rm final}\approx 2$, while the latter $z_{\rm final}=0$. The massive simulations span more than a decade in mass, with $z=2$ halo masses of M_{\rm DM}\approx 1\times 10^{12}-3\times 10^{13}\mbox{M{}_{\sun} }, while the low-mass (low redshift) simulation has a final halo mass comparable to the Milky Way’s (e.g. Boylan-Kolchin et al., 2013).

Utilising caesar (Thompson, 2015), we follow the standard procedure from Hahn & Abel (2011) in which we build an ellipsoidal mask around all particles within $2.5\times$ the radius of the maximum distance (i.e. the farthest away) dark matter particle in the $z=z_{\rm final}$ halo, and define this as the Lagrangian high-resolution region to be re-simulated at a higher-resolution. When re-run to $z=z_{\rm final}$, we have zero low-resolution particles contaminating our main halos within three virial radii.

We operate gizmo in the meshless finite mass (MFM) mode, which evolves the fluid in a mass-conserving manner within individual mesh nodes. We utilise a cubic spline kernel with 64 neighbours in the MFM hydrodynamics, which is used to define the volume partition between gas elements, and therefore the faces over which the hydrodynamics is solved via the Riemann solver. Davé et al. (2016) have shown that relatively minimal differences exist when using the MFM solvers in galaxy formation problems as a pressure-entropy SPH (Hopkins, 2013).

Our simulations use the suite of physics described in the mufasa cosmological hydrodynamic simulations (Davé et al., 2016, 2017a, 2017b), and indeed our parent dark matter simulation is generated with identical initial conditions as the $50$ Mpc $h^{-1}$ mufasa simulation. In this, stars form only in dense molecular gas, where the H2 fraction is dependent on the gas metallicity and surface density following the Krumholz et al. (2009) prescription (e.g. Thompson et al., 2014). We impose a minimum metallicity for star formation of $Z=10^{-3}Z_{\odot}$.

Star formation proceeds following a Schmidt (1959) relation:

[TABLE]

where $f_{\rm H2}$ is the mass fraction of gas in a given element that is molecular, $t_{\rm dyn}$ is the local dynamical time scale, and $\epsilon_{*}$ is the dimensionless star formation efficiency per free fall time. We set this latter quantity to be $\epsilon_{*}=0.02$ as motivated by observations and theoretical modeling (Kennicutt, 1998; Narayanan et al., 2008, 2012; Hopkins et al., 2013).

Feedback from young stars are modeled using a decoupled, two-phase wind. The wind physics are described in significant detail in Davé et al. (2016), and we refer the reader to that work, summarising only the most relevant details here. The modeled stellar winds have a probability for ejection which is modeled as a fraction ($\eta$) times the star formation rate probability. $\eta$ is modeled from the best-fit relation by Muratov et al. (2015), from the Feedback in Realistic Environments (fire) simulation suite:

[TABLE]

which is additionally motivated by analytic arguments (Hayward & Hopkins, 2017). The stellar mass of the galaxy is determined using a fast on-the-fly friends of friends finder (Davé et al., 2016). The ejection velocity scales in a modified way with the circular velocity of the galaxy as well, following the scaling relations derived by Muratov et al. (2015). 30% of the winds are ejected hot, with a thermal energy given by the difference between the supernova energy and the kinetic energy of the launch. These winds are decoupled from hydrodynamic forces or cooling until one of three conditions are met: ($i$) its relative velocity compared to the surrounding gas is less than half the local sound speed; ($ii$) the wind reaches a density limit of $1\%$ of the critical density for star formation; ($iii$) $2\%$ of the Hubble time at launch has elapsed.

As we will discuss in § 2.3, the dust masses in our radiative transfer simulations are tied to the metal content in the ISM. Feedback from longer-lived stars (Asymptotic Giant Branch [AGB], and Type Ia supernovae) are included as well, following Bruzual & Charlot (2003) stellar evolution tracks with a Chabrier (2003) initial mass function. These delayed feedback sources deposit metals into the interstellar medium. We track the evolution of 11 elements: H, He, C, N, O, Ne, Mg, Si, S, Ca and Fe. The type II Supernovae yields follow the Nomoto et al. (2006) parameterisations as a function of metallicity, though reduced by a factor $50\%$, owing to studies that show these scalings result in metallicities a factor $\sim 2$ too high as compared to the observed mass-metallicity relation (Davé et al., 2011). For Type Ia supernovae, we employ the yields from Iwamoto et al. (1999), assuming 1.4 M*☉*  of returned metals per SN1a event. Chemical enrichment from AGB stars follows the Oppenheimer & Davé (2008) lookup tables where the yields are a function of age and metallicity.

The default mufasa model includes a heuristic quenching scheme in which all gas that is not self-shielded is heated in massive halos (where the threshold mass for quenching is a function of redshift Davé et al., 2016). For the bulk of our models we do not include this scheme as we endeavour to simulate a broad range of physical conditions in star-forming galaxies, and including the quenching model would limit the parameter space that we can study. This said, in order to build a substantial population of evolved stars (that will enable us to understand the relationship between old stellar populations and the IRX-$\beta$ relation), we include one model in which this quenching model is included.

In Table 1, we summarise the model zooms employed in this study. The dark matter particle masses are $M_{\rm DM}=1\times 10^{6}h^{-1}M_{\odot}$, and baryon particle masses $M_{\rm b}=1.9\times 10^{5}h^{-1}M_{\odot}$. We employ adaptive gravitational softening of all particles throughout the simulation (Hopkins, 2015), with minimum force softening lengths of $12,3,$ and $3$ pc for dark matter, gas and stars respectively. Each model galaxy is evolved to $z=2$, outputting $85$ snapshots from $z\sim 30-2$. In Figure 1, we show the gas surface density projections for one of our model galaxies (mz5) as it evolves from $z\sim 5-2$.

2.3 Dust Radiative Transfer

With a galaxy sample in hand, we proceed to calculating the UV-far infrared SED. Our methodology contains two major steps: generating the unattenuated stellar SEDs for all star particles in the simulation, and then propagating these stellar SEDs through the dusty interstellar medium to calculate the resultant attenuation. This process is performed on each snapshot at $z<10$ in post-processing.

These combined calculations are performed with powderday, an open-source dust radiative transfer package first described in Narayanan et al. (2015) that wraps a number of publicly available codes as described below. The stellar SEDs for each snapshot of each modeled galaxy are calculated with the Flexible Stellar Population Synthesis (fsps) code (Conroy et al., 2009, 2010; Conroy & Gunn, 2010) via python-fsps111https://github.com/dfm/python-fsps, a set of Python hooks for fsps. The stellar SEDs are calculated as simple stellar populations with their ages and metallicities taken directly from the hydrodynamic simulation. We assume a Kroupa (2002) stellar initial mass function (IMF), and the Padova isochrones (Marigo & Girardi, 2007; Marigo et al., 2008).

We calculate the attenuation suffered by these stellar SEDs by performing 3D Monte Carlo dust radiative transfer simulations. The metal mass from the galaxy formation simulations are projected onto an adaptive grid with an octree memory structure. The grid encompasses $25$ kpc around the center of the galaxy. Functionally, the octree is constructed by placing the entire set of particles associated with the halo within the radiative transfer simulation region into a single cell, and then recursively subdividing this cell into octs until a threshold maximum number of gas particles are in the cell (we employ $n_{\rm subdivide,thresh}=64$). The physical properties of the gas are projected using a spline smoothing kernel, and functionally done with software within yt (Turk et al., 2011, M.J. Turk et al. in prep.).

The radiative transfer happens in a Monte Carlo fashion utilising the dust radiative transfer code hyperion (Robitaille, 2011). hyperion uses the Lucy (1999) iterative algorithm for determining the radiative equilibrium temperature in dust grains. In practice, photons are emitted from stellar sources, and absorbed, scattered, and re-emitted by dust. This process is iterated upon until the dust temperature and radiation field are converged. We determine convergence when the energy absorbed by $99\%$ of the cells has changed by less than $1\%$ between iterations. Unless otherwise specified (i.e. the tests done in § 4.2.3), we assume a Weingartner & Draine (2001) grain size distribution, with $R\equiv A_{v}/E(B-V)=3.15$ (where $A_{v}$ is the visual extinction, and $E(B-V)$ is the difference between the $B$ and $V$ band extinctions. The emissivities are updated to include an approximation for polycyclic aromatic hydrocarbons (PAHs) following Robitaille et al. (2012). The dust mass is assumed to be a constant fraction (here, $40\%$) of the metal mass in each oct cell, where this fraction is taken from constraints across a range of epochs (Dwek, 1998; Vladilo, 1998; Watson, 2011).

The result of the powderday calculations is a simulated SED for each galaxy snapshot, spanning the wavelength range $912\mbox{\AA}$ to $1$ mm. It is these SEDs that we use as both a comparison to observations, as well as an interpretative tool for understanding the IRX-$\beta$ relation in galaxies. Example SEDs can be seen in Figure 4.

3 IRX-$\beta$ Fundamentals

3.1 IRX-$\beta$ Calculations

The infrared excess (IRX) is defined as:

[TABLE]

We define the infrared luminosity as the bolometric luminosity between $8$ and $1000$ µm, though note that our results are not significantly changed if we use a minimum wavelength of $\lambda=1$µm in our definition of $L_{\rm IR}$. Here, we use the flux at the closest wavelength output from our models to $1600\mbox{\AA}$, which is $1601.03\mbox{\AA}$.

In order to calculate the UV continuum slope, $\beta$, we fit a simple power law $f(\lambda)\propto\lambda^{\beta}$ between $\lambda=1000-3000\mbox{\AA}$ (where the flux is in units of erg s*-1* cm*-2* $\mbox{\AA}^{-1}$). We choose this relatively broad range to marginalize the effects of the $2175\mbox{\AA}$ dust feature, when present.

3.2 The Reference Relations

Throughout the paper we will compare to a set of ’reference’ relations to help put our results into context. These include: (i) the original Meurer et al. (1999) relation, which is a fit to observations of local starburst galaxies; (ii) the Casey et al. (2014b) relation, which is additionally a fit to low redshift galaxies, but includes a larger dynamic range of SFRs as well as aperture-corrected data from heterogeneous samples; (iii) recent calibrations by Takeuchi et al. (2012) and Overzier et al. (2011); and (iv) SMC curves calculated by Siana et al. (2009). For the SMC curves, Siana et al. utilise constant star formation history Bruzual & Charlot (2003) population synthesis models ($t_{\rm age}=300$ Myr) with an SMC extinction curve.

3.3 The Parent Sample in the IRX-$\beta$ Plane

In Figure 2, we show the location in the IRX-$\beta$ plane of all snapshots of our simulated galaxies (orange points). Alongside our simulated galaxies, we show a compilation of observed data (blue), comprised of the data presented in Casey et al. (2014b), which in itself is a compilation of new observations, and the samples of Gil de Paz et al. (2007); Howell et al. (2010); Reddy et al. (2012); Penner et al. (2012) and Heinis et al. (2013). We supplement these with more recent detections by Capak et al. (2015) and Bouwens et al. (2016). We have omitted upper limits for the time being, but will return to these later in this paper. To guide the eye, we additionally plot the reference relations (c.f. § 3.2). The comparison of our model galaxies to the observed IRX-$\beta$ plane presents a number of salient points.

First, the simulated galaxies naturally lie on a locus similar to the local reference relations, though there are a number of deviations both at high IRX, and low IRX. The dispersion is significant, and owes to a range of physical effects that we will explore in turn in this paper.

Second, while our galaxies follow the same generic trend as observed galaxies, they do not fully cover the same extent in IRX and $\beta$ as observed galaxies. We quantify this in Figure 3, where we show histograms of the $L_{\rm IR}$, $L_{\rm UV}$, IRX and $\beta$ of both our model galaxies, and the observed comparison samples. While our galaxies are comparably UV-luminous, we do not achieve infrared luminosities much beyond L_{\rm IR}\geq 10^{13}\mbox{L{}_{\sun}}, and consequently do not have model galaxies with as extreme IRX values as the most extreme outliers in the observed sample. Of course a lack of the dustiest systems in our model sample also results in model galaxies restricted to a a $\beta$ range $-2<\beta<1$, therefore falling short of the very reddest UV-continuum slopes observed.

The lack of exact range matching is not a failure of the models, but rather to be taken as a comparison of sample selections: our models contain, on average, galaxies of roughly $5-10\times$ lower infrared luminosity than the average observed galaxy that we compare to. This said, despite the fact that the simulations do not cover the entire range of values seen in the observations, significant insight can still be gained.

In the remainder of this paper, we will dissect Figure 2, aiming to understand both the origin of the IRX-$\beta$ relation in galaxies, as well as trends that result in deviations from the mean relation. To do this, we will no longer present the entire sample, but rather individual sub-samples restricted by either physical or observable properties as it pertains to the relevant discussion.

4 Deconstructing the IRX-$\beta$ Relation

4.1 The Origin of the IRX-$\beta$ Relation

We begin by understanding the origin of the reference IRX-$\beta$ relations; we then proceed to deconstruct large deviations from these relations.

The simplest place to start is with stellar population synthesis experiments. In the top left panel of Figure 4, we have created a simple stellar population with fsps that is covered with a foreground dust screen and a stellar age $t_{\rm age}=10^{-3}$ Gyr. These are denoted with the star symbols (we will discuss the circles, as well as the other panels in the plot shortly). For maximal simplicity, we employ this dust screen to all stars, and manually increase the normalisation (i.e. optical depth) of the attenuation curve to study the location of the emergent SED in the IRX-$\beta$ plane. In our default population synthesis model shown here, old stars ($t_{\rm age}>10^{7}$ yr) do not see any dust, although in this particular example we fix all stellar ages to much less than this threshold value. The dust attenuation curve is a powerlaw with index $-0.7$ (i.e. $\tau(\lambda)\propto\lambda^{-0.7}$), stellar initial mass function of a Kroupa (2002) form, and solar metallicity. The colours of the points denote the monochromatic ($1600\mbox{\AA}$) UV optical depth. As the optical depth for UV photons increases, the UV SED reddens, while at the same time power is transferred from the UV to the IR. This idealized example serves as a control experiment in that it is unmuddied by the complicating effects of a diverse stellar population, and the complex star-dust geometry characteristic of real galaxies.

We further expand on this experiment by conducting a slightly more complex numerical exploration utilising our model galaxies. We now turn to the filled circles in the top left panel of Figure 4. Here, we have selected three arbitrary galaxies from our sample (all selected at the same redshift [$z\sim 2.25$] from galaxies mz0, mz5 and mz10), and plot their location using our fiducial parameters on the IRX-$\beta$ relation. We then manually decrease and increase the dust mass of these galaxies by adjusting our dust to metals ratio and show the location of these galaxies in the IRX-$\beta$ plane. (We remind the reader that our fiducial dust to metals ratio is $40\%$). Alongside this, we show the panchromatic SEDs for these points (top right), the UV optical depth (bottom left; normalised by their $1600\mbox{\AA}$ optical depth, and the un-normalised optical depth (bottom right).

As the dust column increases in Figure 4 the optical depth seen by UV photons increases, and the UV SED flattens. At the same time, power is transferred from the UV into the infrared, resulting in increased IRX values. Of course, Figure 4 represents contrived scenarios to show how galaxies may move along the reference IRX-$\beta$ relations in the face of increasing dust column. However, as is evident both from the observed data, as well as our full suite of numerical simulations (e.g. Figure 2), there is significant dispersion about these reference relations. In the remainder of this paper, we explore the origin of these deviations from the reference IRX-$\beta$ relations.

4.2 Dispersion in the IRX-$\beta$ Relation

4.2.1 Can Simplified Screen Models Capture IRX-$\beta$ Dispersion? The Role of Geometry

Prior to understanding the origin of the dispersion in the IRX-$\beta$ relation, it is first worth considering whether a simple screen model at a single stellar age can capture a similar range of IRX and $\beta$ values as both the cosmological zoom galaxy formation simulations, and the observations.

In Figure 5 we plot the results of fsps models akin to those presented in the top left panel of Figure 4. These models are a simple stellar population with a uniform dust screen with a powerlaw attenuation curve. In the left panel of Figure 5, we vary just the dust attenuation powerlaw index (from [$-0.2,-1.2$]) alongside the attenuation law normalisation. While this model is able to span arbitrarily large $\beta$ values, the very blue UV colours of high IRX galaxies are not encapsulated in this model. In other words, simple screen geometries in population synthesis models fail to capture the complex physics underlying the locations of galaxies on the IRX-$\beta$ plane.

In the right panel of Figure 5, we allow a fraction of young O and B stars ($t_{\rm age}<10^{7}$ yr) in the SPS models to run away from their birth clouds, as well as a fraction of old ($t_{\rm age}>10^{7}$ yr) stars to be decoupled from the diffuse dust emission. The goal here is to decouple the sources of UV luminosity from the obscuring dust. We allow both the fraction of runaway OB stars and the fraction of decoupled old stars from dust to vary between [$0,1$]. In this somewhat extreme experiment, it is clear that the complex star-dust geometry is able to drive a large range of IRX-$\beta$ values, and that the observed dispersion is due to a decoupling of stars from the dust attenuation screen.

4.2.2 The Impact of Older Stellar Populations

We now return to Figure 2, and highlight the locus of galaxies that lie below the reference IRX-$\beta$ relations. These galaxies predominantly arise from our halo0q (quenched) model, and reflect quenched galaxies with older stellar populations. In general, older stellar populations can drive red UV SEDs, even while maintaining relatively low IRX values.

To demonstrate this point more explicitly, in Figure 6, we show the results of a controlled numerical experiment in which we compare the mz0 and mz0_q models. The primary difference in these models is that the latter includes our heuristic quenching model, and therefore has a significant population of evolved stars. The points are coloured by their median stellar age. As is clear, the galaxies with the oldest stellar populations preferentially lie well below the reference IRX-$\beta$ relations.

This effect owes to the reddening of UV SEDs in older stellar populations. To show this, in the left panel of Figure 7, we conduct a numerical experiment with fsps, where we evolve a simple stellar population behind a dust screen. The example is the same as in Figure 4, though we fix the dust optical depth and vary the stellar ages. In the right panel of the same figure, we show the UV and optical SEDs for this aging population. As stellar populations age, their UV SEDs shift toward redder $\beta$ slopes, as demonstrated by the right panel of Figure 7. At the same time, they are less likely to be obscured by dust (as a reminder, in these fsps calculations, stellar populations with $t_{\rm age}>10^{7}$ yr are unobscured), and their IRX ratios drop. (We note there is a slight increase in IRX as the populations continue to age owing to dramatic decreases in $L_{\rm UV}$).

Of course, using the median stellar age as a sole parameter for describing deviations below the IRX-$\beta$ relation is a crude parameter. What actually matters is the stellar population dominating the UV SED, which is a function of the stellar age distribution. Because of that, there is significant dispersion in the ages of galaxies that lie below the reference IRX-$\beta$ relations. This said, the general point stands that, on average, older stellar populations deviate toward redder UV colours. This point has been discussed by Kong et al. (2004), who parameterised deviations from the IRX-$\beta$ relation in terms of a stellar birthrate parameter, comparing the present SFR to the time averaged past.

4.2.3 Dust Composition and the $2175\mbox{\AA}$ Bump

A number of studies have underscored the importance of the underlying extinction curve in driving how galaxies present in the IRX-$\beta$ plane. For example, as previously discussed, Siana et al. (2009) utilise Bruzual & Charlot (2003) calculations to demonstrate that an LMC and SMC-like dust reddening curve result in IRX-$\beta$ relations that lie, for the most part, below the standard Meurer et al. (1999) relation. Similarly, Bell et al. (2002) demonstrate that sightlines in the LMC all lie below the Meurer et al. (1999) galaxies in the IRX-$\beta$ plane.

To understand how the intrinsic extinction curve impacts the observed IRX-$\beta$ relation, we have conducted two numerical experiments on top of our fiducial model (which, as a reminder, utilises a Milky Way Weingartner & Draine (2001) $R_{\rm v}=3.1$ extinction curve). In our first experiment, we utilise the same curve, though have eliminated the $2175\mbox{\AA}$  extinction ’bump’. We remove the bump by linearly interpolating the extinction curve between $1600-4000\mbox{\AA}$ in log-space. This results in reduced attenuation in this wavelength range. For the second experiment, we utilise the Weingartner & Draine (2001) SMC dust curve. In Figure 8, we show these extinction curves (red lines).

In Figure 9, we show the results from these models in the IRX-$\beta$ plane. Galaxies without a $2175\mbox{\AA}$ bump (but otherwise similar dust extinction properties to the Milky Way) has both lower IRX and larger $\beta$ values than our fiducial model that includes the UV extinction bump. The SMC dust curve which has even further reduced extinction in the UV presents with both even lower IRX and redder $\beta$.

To understand the origin of the lower IRX and redder $\beta$ values for the model without a UV bump and the SMC dust curve, we now turn to the the UV SEDs (blue lines) presented in Figure 8. At a fixed far ultraviolet (FUV) opacity, the Milky Way dust curve with no bump has a lower near ultraviolet (NUV) opacity, and the SMC curve a yet lower NUV opacity. As this NUV extinction is reduced, the observed UV SED shows larger powerlaw SED slopes owing to increased transparencies in the NUV/optical bands. Of course, the location of galaxies with the SMC dust curve is degenerate with the impact of older stellar populations (c.f. Figure 7). One must therefore exercise caution in the interpretation of galaxies that lie below the IRX-$\beta$ relation.

5 Application to Observations

5.1 Blue Dusty Star Forming Galaxies

Dusty star forming galaxies are the most luminous, heavily star-forming galaxies in the Universe (e.g. Smail et al., 1997; Barger et al., 1998; Hughes et al., 1998; Hayward et al., 2011; Hayward et al., 2013; Narayanan et al., 2009; Narayanan et al., 2010; Casey et al., 2014a). A number of works recently have observed that dusty star forming galaxies at a range of redshifts exhibit particularly blue UV SEDs given their IRX values as compared to the traditional IRX-$\beta$ relations. In other words, the UV slopes, $\beta$ are too low for DSFGs given their infrared excesses. This was neatly summarised for a large compiled dataset of \mbox{z\sim }0-4 DSFGs by Casey et al. (2014b) who demonstrated that galaxies become bluer (i.e. larger departures from the reference IRX-$\beta$ relations toward lower $\beta$ slopes) with increasing galaxy infrared luminosity. Similar effects have been observed by Penner et al. (2012) and Oteo et al. (2013).

In Figure 10, we investigate the origin of blue infrared bright galaxies. In the top left panel, we plot our parent sample, with galaxies colour-coded by their infrared luminosity. Here, we only plot galaxies at $z>2$ to remain consistent with the bulk of the observations in these samples. Immediately it is clear that more infrared-bright systems present both above and to the left of the reference IRX-$\beta$ relations. This is further quantified in the top right panel of Figure 10, where we show the departure from the reference IRX-$\beta$ relations (here, we use the Casey et al. (2014b) relation as the reference from which we measure $\beta$ departures) as a function of galaxy infrared luminosity. The light blue points show observations, while the coloured points are our cosmological zoom simulations. We ignore the actual colours of the simulated points for the time being. Regardless, it is evident that in both observed galaxies, as well as our models, more infrared-luminous systems indeed have bluer UV SEDs than the reference IRX-$\beta$ relation. Moreover, the magnitude of this departure in our simulations shows good agreement with what is observed.

We first explore the possibility that frosting of UV bright regions toward the outskirts of galaxies causes the blue offset for high infrared luminosity galaxies at high-redshift. Massive galaxies (i.e. the kind that typically present as DSFGs at high-$z$) live in complex environments, with significant substructure surrounding the central galaxy (e.g. Davé et al., 2010; Narayanan et al., 2015). As demonstrated by Geach et al. (2016), many of these star-forming subhalos are metal poor, and UV-bright. In the top right panel of Figure 10, we now highlight the colour-coding, which maps to the UV-luminosity weighted radius of our model galaxies. Indeed, this increases by a factor $\sim 5$ over the dynamic range modeled, with more infrared-luminous galaxies broadly having larger UV radii.

However, while this likely plays some role in driving the blue offset of DSFGs from the reference IRX-$\beta$ relations, in detail the rise in UV disk sizes for more infrared luminous galaxies is tempered by the fact that the galaxies with the highest infrared luminosity also tend to be the most massive, and hence the largest. More colloquially, bigger things are bigger. We demonstrate this in the bottom right panel of Figure 10, where we have now normalised the UV luminosity weighted radii of our model galaxies by their virial radii. Indeed, the galaxies with the largest $R_{\rm UV}/R_{\rm vir}$ ratios reside in the bluest, most infrared-luminous locus in $\Delta\beta$-LIR plane. However, there are a number of systems without notably large UV radii. Hence, frosting plays a role in driving some DSFGs toward bluer UV colours, though does not sufficiently paint the entire picture.

Instead, also important are low optical depth sightlines reflective of a complex star-dust geometry in these systems. In Figure 1, we showed the gas phase morphology for an arbitrary sample of model galaxies. The geometries are complex, and result in significantly lower UV optical depths from decoupled UV and IR emission sites than what would be expected for a maximally obscured geometry. We demonstrate this quantitatively in the bottom-left panel of Figure 10, where we show the monochromatic UV ($1600\ \mbox{\AA}$) optical depth as a function of the infrared luminosity (normalised by galaxy stellar mass) for our model galaxies. The rise in $\tau_{1600}$ is shallow with increasing $L_{\rm IR}/M_{*}$. As a comparison, we revisit the fsps numerical experiment conducted in the top left panel of Figure 4, where we examined a simple stellar population behind a screen of dust with increasing optical depth, and plot the optical depth of this population as a function of LIR/$M_{*}$ alongside our simulations. Here, we have slightly modified the numerical experiment from Figure 4, in that we enforce that dust dust obscures all stars, both young and old; this represents maximal obscuration. There is a dramatic divergence in the $1600\mbox{\AA}$ optical depths at high specific infrared luminosity between the simplified maximally obscured stellar population synthesis model, and our hydrodynamic simulations that exhibit more complex geometries. At the highest specific infrared luminosities, the maximal optical depth is $\tau_{1600}\sim 4-5$, reflective of the large number of sightlines toward UV sources that are relatively unobscured.

In short, we find that dusty star forming galaxies at high-$z$ exhibit blue UV slopes due to low optical depth sightlines toward UV-bright regions. A similar origin is hypothesized by Casey et al. (2014b), though we note that Safarzadeh et al. (2017) suggest that frosting (i.e. recent star formation toward the outskirts of galaxies) may instead dominate the origin of blue DSFGs.

5.2 Galaxies in the First Billion Years

We now turn our attention to the IRX-$\beta$ relation in galaxies at $z\ga 5$ (for recent reviews on UV slopes at these redshifts, see Stark, 2016; Finkelstein, 2016). The advent of ALMA has allowed for rest-frame far infrared detections of galaxies in the early Universe that, coupled with a number of assumptions, allows for FIR luminosity measurements. As a result, recent years have seen a number of constraints on the IRX-$\beta$ relation in low metallicity systems in this epoch. By and large, these galaxies appear to lie below the locus of more metal rich galaxies at $z\sim 0-2$ (and even below the canonical SMC IRX-$\beta$ curve).

As an example, Capak et al. (2015) surveyed nine galaxies (between $1-4\ L_{*}$ and UV slopes $-1.4<\beta<-0.7$) at $z\approx 5-6$ in ALMA’s band 7 ($0.8-1.1$ mm observed frame) and found that seven of them lie below the SMC IRX-$\beta$ line (the other two being consistent with the SMC curve). Similarly, Bouwens et al. (2016) utilised the ALMA Spectroscopic Survey ($2\sigma$ detections) in the Hubble Ultra Deep Field (ASPECS; Walter et al., 2016; Aravena et al., 2016a, b; Decarli et al., 2016) to constrain the IRX-$\beta$ relation in galaxies from $z\sim 2-10$. Like the Capak et al. (2015) study, Bouwens et al. (2016) found many galaxies in this redshift range to be roughly consistent with the SMC curve (the more massive systems), or to lie below it.

In Figure 11 we summarise both the detections (and non-detections) in the aforementioned Capak et al. (2015) and Bouwens et al. (2016) studies, as well as show the location of our own high-redshift model galaxies on the IRX-$\beta$ plane. Here, we show our model systems between $5\leq z\leq 7$, and for reference, include the Casey et al. (2014b) compilation of $z\sim 0-2$ detections. We restrict our analysis to models mz0, mz5 and mz10 as these models best overlap the stellar mass range observed by Capak et al. (2015) and Bouwens et al. (2016) (M_{*}>10^{9}\mbox{M{}_{\sun} } at $z>5$). As is clear, our model galaxies all lie on a relatively tight locus near (or above) the $z\sim 0$ relations, while many of the $z\ga 5$ observations have much lower IRX values.

At face value, the observed lack of dust content in $z\ga 5$ galaxies presents a challenge to our models, which predict that these galaxies ought to lie on a similar IRX-$\beta$ relation as local systems (albeit toward the low IRX and $\beta$ end). However, one possibility is that the observed galaxies have underestimated infrared luminosities (and hence, IRX values) due to relatively cold assumed dust temperatures.

To elaborate, in deriving the infrared luminosity, Capak et al. (2015) and Bouwens et al. (2016) both required an assumed $T_{\rm dust}$ due to the fact that only one photometric FIR data point was available. Capak et al. (2015) assumed 25 K $\leq T_{\rm dust}\leq$ 45 K, and Bouwens et al. (2016) assumed as a fiducial value, $T_{\rm dust}\sim 35$ K. Of course if the true dust temperature is warmer than this, then the derived $L_{\rm IR}$ will increase, and the observed galaxies will move closer to the local IRX-$\beta$ relation.

In the right panel of  11, we plot the distribution of dust temperatures from our model galaxies, showing both the entire $z>2$ sample, as well as just the $z>5$ sample222To ensure that we don’t miss any cold dust outside our fiducial $25$ kpc box, we utilise a larger $200$ kpc box for the radiative transfer calculations for our $T_{\rm dust}$ derivation.. The dust temperatures are derived simply from the peak of the modeled infrared SED. While lower redshift ($z\sim 2$) galaxies have dust temperatures comparable to the assumed observed values of $35-60$ K, at $z\ga 5$, the dust temperatures rise dramatically, with typical dust temperatures $50<T_{\rm dust}<75$ K. This owes both to lower dust masses (so UV photons have lower optical depths), as well as harder radiation fields in the lower metallicity stars.

On average, dust temperatures drop individual galaxies with decreasing redshift. We show this general case in Figure 12, where we show the dust temperature for all of the galaxies used in Figure 11, binned by their stellar mass at the redshift shown. A generic trend is that with increasing redshifts, the average dust temperature of galaxies rises. Indeed, Bouwens et al. (2016) find that their upper limits would be $\sim 0.4$ dex higher if they had assumed a $T_{\rm dust}=45-50$ K, which is still well below the maximum dust temperatures we find at $z\gtrsim 5$ of $\sim 75$ K. Even an increase of $\sim 0.4$ dex would be enough to bring the bulk of the upper limits to values consistent with the SMC curve.

In summary, we therefore find that the typical dust temperature of galaxies at $z\gtrsim 5$ is substantially larger than what is observed at lower ($z\sim 2$) redshifts. When accounting for this increased $T_{\rm dust}$, inferred infrared luminosities from these systems rise and bring IRX values of observed galaxies consistent with local IRX-$\beta$ relations.

6 Discussion

6.1 Secondary (and tertiary) parameters in the Relation

While the IRX-$\beta$ relation appears to hold over a range of physical conditions, the numerical experiments conducted thus far have demonstrated a number of scenarios in which galaxies may deviate strongly from the reference IRX-$\beta$ relations. These deviations have been observed in the literature for some time.

As a result, a number of authors have investigated the possibility of “second parameters” in the IRX-$\beta$ relation. For example, Kong et al. (2004) and Mao et al. (2014) have studied the role of star formation history in driving dispersion in the UV colours in the IRX-$\beta$ relation. Kong et al. (2004) defined the ’birthrate parameter’, proportional to ratio of the present to past average SFR. Galaxies with lower birthrate parameters move farther from the reference IRX-$\beta$ relations in a perpendicular manner. At the same time, however, Johnson et al. (2007a, 2013) find little correlation with the dispersion in $\beta$ and the star formation history for as parameterised by the $4000\mbox{\AA}$ break for blue cloud galaxies (though substantial deviations from the reference IRX-$\beta$ relations for red-sequence galaxies). Similar results were obtained by Seibert et al. (2005); Dale et al. (2007) and Cortese et al. (2008). Grasha et al. (2013) instead find a better relation between the mean stellar age and deviations from the reference IRX-$\beta$ relations. Indeed, our own models suggest that the median stellar age is quite relevant in driving red UV slopes, and thereby significant dispersion about the reference IRX-$\beta$ relations (e.g. Figure 6).

Similarly, a number of authors have investigated the dust attenuation law as a possible driver for dispersion in the IRX-$\beta$ relation (e.g. Johnson et al., 2007b; Burgarella et al., 2005; Boquien et al., 2009). The dust attenuation law in galaxies encapsulates both the loss of photons along the line of sight, as well as scattering back into the line of sight, and therefore is a proxy for both the extinction law, as well as radiative transfer and geometric effects.

In our simulations we can identify three major physical drivers of deviations from the IRX-$\beta$ relation. These are:

Old stellar populations. 2. 2.

Complex geometries of highly star-forming galaxies (and the resultant decoupled UV and IR emission sites) 3. 3.

Dust extinction curves that deviate substantially from a Milky Way-like curve.

In Figure 13, we summarise the direction of these deviations schematically.

Because the variation of dust extinction curves with galaxy physical property is relatively unconstrained (e.g. Kriek & Conroy, 2013), we will concentrate on the impact of older stellar populations and star-ISM geometry in the remainder of our discussion of additional parameters in the IRX-$\beta$ relation. While observationally resolving the complexity of the star-ISM geometry in galaxies (especially at high-$z$) is infeasible for most systems, Figure 10 shows that the $\beta$ deviation from the reference IRX-$\beta$ relations ($\beta_{\rm ref}$) is well parameterised by the total infrared luminosity of the galaxy. It is apparent that the \Delta\beta_{\rm ref}\propto\mbox{L{}_{\rm IR}}. Because the star formation rate is reasonably well represented by the infrared luminosity in galaxies (within limits; see, e.g. Younger et al., 2009; Hayward et al., 2014; Narayanan et al., 2015), then we can fit $\Delta\beta_{\rm ref}$ as a function of the galaxy SFR. Numerically, we utilise the best fit relationship from the compendium by Casey et al. (2014b) as the reference IRX-$\beta$ relation from which we measure deviations in $\Delta\beta_{\rm ref}$:

[TABLE]

and find a relatively good fit with:

[TABLE]

In Figure 14, we show the IRX-$\beta$ relation after applying the correction from Equation 5 (green hexagons). We find that our model parameterisation for $\Delta\beta$ (SFR, $t_{\rm age}$) does a reasonable job at collapsing the deviations from our reference Casey et al. (2014b) IRX-$\beta$ relation back toward the reference, and substantially reducing the dispersion. We additionally apply the same fit from Equation 5 to the observational data in Figure 14, noting that the dispersion about the reference relations in observations as well is significantly reduced. For comparison, in the same figure, we show the fiducial IRX-$\beta$ relation from our simulations (orange circles), without accounting for the blue $\beta$-offset in heavily star forming galaxies333For the observed galaxies, we employ the Murphy et al. (2009) LIR-SFR calibrations as summarised in Kennicutt & Evans (2012). We assume half of the radiation from star formation is reprocessed in the infrared. Because we have no information about the stellar ages, we assume $250$ Myr stellar ages for all galaxies.. By adding the $\Delta\beta$ derived from Equation 5 to the observed value, one can significantly reduce the dispersion and uncertainty inherent in the IRX-$\beta$ relation.

6.2 IRX or $\beta$?

At this point, a natural question is: does IRX or $\beta$ serve as a better proxy for the UV optical depth? Indeed, this point has been investigated by a number of observational studies in recent years (e.g. Cortese et al., 2008; Kennicutt et al., 2009; Hao et al., 2011)

In Figure 15, we plot $\tau_{1600}$ against both $\beta$ and IRX, colour-coding the points by $\Delta\beta_{\rm ref}$, where, again we use the Casey et al. (2014b) derived relation as our reference relation. As is clear, $\beta$ generally does a poor job of representing the monochromatic UV optical depth. The reason for this is discussed in § 5.1. In short, strongly varying star-dust geometries complicate the interpretation of $\beta$.

At the same time, IRX serves as a reasonable proxy for $\tau_{1600}$. We find a relation:

[TABLE]

well describes the $\tau_{1600}$-IRX relation, and we include this fit in the right panel (solid blue line) of Figure 15.

6.3 Comparison to Other Theoretical Models

A number of works in recent years have attempted to understand the IRX-$\beta$ relation in galaxies from a theoretical standpoint. By and large, these methods have employed complementary techniques to our tool of choice: cosmological hydrodynamic zoom simulations.

In seminal work, Granato et al. (2000) coupled the galform semi-analytic galaxy formation model with dust radiative transfer calculations to develop a model for infrared-luminous galaxies (with a particular eye toward local ULIRGs). Semi-analytic models necessarily require a simplified star-dust geometry as the structure of galaxies are evolved in an analytic method. To wit, Granato et al. (2000) assume that the gas and dust are distributed in exponential disks, with the dust distributed both in molecular clouds and diffuse (cirrus) ISM. Stars are allowed to live within molecular clouds within some time scale $t_{\rm esc}$. These galaxies are evolved through cosmic time (with star formation rates/histories determined by the physics inherent in the analytic prescriptions combined with the modeled halo growth rates), and processed through grasil dust radiative transfer. grasil, similar to powderday, calculates the stellar SEDs, and computes the radiative transfer through the dusty ISM.

Granato et al. (2000) found that starbursting galaxies within the galform semi-analytic framework nearly always lie close to the fiducial Meurer et al. (1999) IRX-$\beta$ relationship when the stellar population is young and obscured. These authors find that the time scale for which young stars are obscured by their birth clouds, $t_{\rm esc}$, is one of the more important parameters in dictating deviations from the reference IRX-$\beta$ relations. These results appear to be in good agreement with our own model interpretations. Young stars that are well hidden by a screen of dust (c.f. Figure 4) tend to move along the reference IRX-$\beta$ relations. When the UV and IR emitting regions become decoupled (c.f. Figure 5), then the population exhibits large deviations in the IRX-$\beta$ plane.

More recent work by Safarzadeh et al. (2017) utilized a set of idealized hydrodynamic galaxy evolution simulations coupled with sunrise dust radiative transfer to predict the UV-infrared SED properties of both idealized disk galaxies, and 1:1 major galaxy mergers. Here, the star-dust geometry is determined by the hydrodynamic evolution of the galaxy/galaxies and has no interaction with cosmologically infalling gas (this said, it is important to note that substantial radial variations in the ISM properties of idealized dics and mergers can still occur Torrey et al., 2011; Narayanan et al., 2009, 2011). Amongst other issues, these authors investigated both the origin of blue DSFGs in their simulations. Safarzadeh et al. (2017) show that their models are able to produce galaxies with blue UV colours, akin to those in the Casey et al. (2014b) study, though attribute the origin to ’frosting’ of young stars (i.e. young stars that lie outside the central dust-obscured nucleus). This is in contrast to our model, which demonstrates (Figure 10) that, when normalised by galaxy virial radius, frosting does not dominate the blue UV colours of highly star-forming galaxies in a cosmological context.

Safarzadeh et al. (2017) additionally investigated the usage of an SMC dust curve on the location of their idealized model galaxies in the IRX-$\beta$ plane. Similar to the results found here, these authors showed that the lack of NUV attenuation drives galaxies toward lower IRX and redder $\beta$ slopes for the bulk of their star-forming lives.

Ferrara et al. (2016) developed an analytic model for dust growth in galaxies with the purpose of understanding Capak et al. (2015) and Bouwens et al. (2016) $z\ga 5$ galaxies that lie below the reference IRX-$\beta$ relations. This model derives an equilibrium dust temperature by balancing the rate of energy absorption with that of emission by dust grains. Ferrara et al. (2016) find that galaxies with significant diffuse dust exposed to the UV radiation fields typical of $z\sim 5-6$ Lyman-break galaxies will have warm dust temperatures, and in effect have true infrared luminosities much larger than what is inferred when assuming colder dust temperatures typical of lower redshift galaxies. This is consistent with our interpretation of galaxies at $z\ga 5$. At the same time, these authors point out that if the clouds in $z\ga 5$ galaxies are denser, then the majority of the dust will be shielded from the UV radiation field, and at lower temperatures. These galaxies will have low IRX values, consistent with the Capak et al. (2015) and Bouwens et al. (2016) interpretations.

6.4 Possible Tensions with Observations

We now turn to areas where our models may conflict with observational results. Reddy et al. (2010) examined a sample of $90$ Lyman break galaxies at $z\sim 2$ that had both infrared (Spitzer) and UV measurements. With the aid of SED modeling, these authors determined that young galaxies ($t_{\rm age}<100$) Myr have systematically lower IRX and redder $\beta$ values than the older systems ($t_{\rm age}>100$ Myr in their sample. At face value, this is in conflict with our model results presented in Figures 6 and  7.

One possibility is that the young galaxies in the Reddy et al. (2010) sample have an intrinsic SMC-type extinction curve. When these authors calculate IRX and $\beta$ for their young galaxies assuming an SMC-like dust curve, the observed young galaxies all lie relatively close to the SMC curve, which could then be accommodated by our model (c.f. Figure 9). While it is unknown what drives variations in extinction and attenuation laws in galaxies, certainly a broad range of attenuation laws are observed in galaxies at low and high-redshift. Indeed, Kriek & Conroy (2013) find that galaxies with larger specific star formation rates tend to have shallower dust attenuation curves. If the younger systems in the Reddy et al. (2010) sample have shallower attenuation curves (i.e. comparable to an SMC curve), then the tension between our model and the Reddy et al. (2010) observations may be reduced. This said, Reddy et al. (2006) and Reddy et al. (2010) suggest that attenuation laws are steeper for younger UV-selected galaxies.

7 Summary

We have developed a theoretical model for the origin of and variations in the IRX-$\beta$ dust attenuation relation in galaxies. To do this, we have combined cosmological zoom hydrodynamic galaxy formation simulations with stellar population synthesis models and 3D dust radiative transfer calculations. Additionally, we compare these to a variety of stellar population synthesis models. Our main results follow:

Galaxies with relatively young stellar populations, a Milky Way-like extinction curve, and relatively cospatial IR and UV emitting regions tend to lie on or near the standard Meurer et al. (1999) or Casey et al. (2014b) local relations. As the dust content in galaxies increase, galaxies move along these relations. 2. 2.

Substantial variations can be present from these reference relations. These generally owe their origin to the following effects:

• •

Older stellar populations tend to lie below the canonical IRX-$\beta$ relations, due to changes in the UV-optical SEDs of evolved stars.

• •

Complex star-dust geometries tend to drive galaxies above the canonical relations, due to low optical-depth sightlines that cause galaxies to have bluer $\beta$ values than those that have more cospatial IR and UV emitting regions. These complex geometries are common in high-redshift heavily star-forming galaxies.

• •

Galaxies with SMC-like extinction curves, or those with a Milky Way-like extinction curve (but without a $2175\mbox{\AA}$ UV bump) lie below the canonical IRX-$\beta$ relations. 3. 3.

We have used these results to understand the origin of deviations from the IRX-$\beta$ relation in both high-redshift dusty star forming galaxies, as well as those detected at $z>5$. The former class of galaxies have relatively blue UV SEDs due to complex star-dust geometries, and low optical depth sightlines toward UV bright regions. The latter class of galaxies tend to lie on or even above the canonical IRX-$\beta$ relation. This said, their dust temperatures are quite warm ($50-70$ K), so IRX inferences based on a single long-wavelength photometric point are subject to systematically underestimating the LIR if assuming a lower $T_{\rm dust}$. 4. 4.

We use the results from these simulations to derive two fitting relations to maximize the utility of the IRX-$\beta$ relation in galaxies.

• •

We derive a fitting relation (Equation 5) between $\Delta\beta_{\rm ref}$, $t_{\rm age}$ and the SFR (i.e. the deviation in $\beta$ from the reference relations due to complex geometries [characterized by the galaxy SFR] and the stellar ages). The usage of this correction factor reduces the dispersion inherent in the IRX-$\beta$ relation.

• •

We derive a fit between the UV optical depth ($\tau_{1600}$), and IRX (which serves as a good proxy for the UV optical depth); Equation 6.

In general, our models suggest that no single IRX-$\beta$ relation exists for galaxies, and that some caution must be exercised when correcting for dust obscuration using this method. We find that it is possible to correct for some of these deviations, though a reasonable amount of scatter in the relationship is still present.

Acknowledgements

The authors thank Caitlin Casey, Alex Pope, Naveen Reddy, Nick Scoville, Brian Siana and Dan Stark for helpful conversations. We are additionally appreciative to Rychard Bouwens, Caitlin Casey, Brian Siana, and Fabian Walter for providing the results from their observational studies in easily digestable digital format. D.N. thanks the organizers of the “Physical Characteristics of Normal Galaxies at z$$>2” conference in Leiden, where many of the ideas for this paper were generated. The simulations published here were run on the University of Florida HiPerGator supercomputing facility, and the authors are grateful to the University of Florida Research Computing for providing computational resources and support that have contributed to the research results reported in this publication. Partial support for DN was provided by NSF AST-1009452, AST-1445357, AST-1724864, NASA HST AR-13906.001 from the Space Telescope Science Institute, which is operated by the Association of University for Research in Astronomy, Incorporated, under NASA Contract NAS5-26555, and a Cottrell College Science Award, awarded by the Research Corporation for Science Advancement. JEG acknowledges the Royal Society for a University Research Fellowship.