The First Billion Years project: constraining the dust attenuation law of star-forming galaxies at z $\simeq$ 5

TL;DR

This study investigates the dust attenuation law at redshift 5 using synthetic galaxy spectra from simulations, finding that a grey, mass-dependent attenuation curve best explains observed galaxy properties and is consistent with recent ALMA data.

Contribution

It introduces a simple, mass-dependent dust attenuation model at z~5 that aligns with observations and rules out steeper extinction curves like the SMC.

Findings

A mass-dependent UV attenuation law reproduces the luminosity function and color-magnitude relation.

The attenuation curve slope n is between -0.7 and -0.3, consistent with Calzetti law.

Steeper curves like the SMC are ruled out by the data.

Abstract

We present the results of a study investigating the dust attenuation law at $z\simeq 5$, based on synthetic spectral energy distributions (SEDs) calculated for a sample of N=498 galaxies drawn from the First Billion Years (FiBY) simulation project. The simulated galaxies at $z\simeq 5$, which have M$_{1500} \leq -18.0$ and $7.5 \leq \rm{log(M/M}_{\odot}\rm{)} \leq 10.2$, display a mass-dependent $\alpha$-enhancement, with a median value of $[\alpha/\rm{Fe}]_{z=5}~\simeq~4~\times~[\alpha/\rm{Fe}]_{Z_{\odot}}$. The median Fe/H ratio of the simulated galaxies is $0.14\pm0.05$ which, even including the effects of nebular continuum, produces steep intrinsic UV continuum slopes; $\langle \beta_{i} \rangle = -2.4 \pm 0.05$. Using a set of simple dust attenuation models, in which the wavelength-dependent attenuation is assumed to be of the form $A(\lambda) \propto \lambda^{n}$, we explore the parameter values which best reproduce the observed $z=5$ luminosity function (LF) and colour-magnitude relation (CMR). We find that a simple model in which the absolute UV attenuation is a linearly increasing function of log stellar mass, and the dust attenuation slope ($n$) is within the range $-0.7 \leq n \leq-0.3$, can successfully reproduce the LF and CMR over a wide range of stellar population synthesis model (SPS) assumptions. This range of attenuation curves is consistent with a power-law fit to the Calzetti attenuation law in the UV ($n=-0.55$), and other similarly `grey' star-forming galaxy attenuation curves recently derived at $z\simeq2$. In contrast, attenuation curves as steep as the Small Magellanic Cloud (SMC) extinction curve ($n=-1.24$) are formally ruled out. Finally, we show that our models are consistent with recent 1.3mm ALMA observations of the Hubble Ultra Deep Field (HUDF), and predict the form of the $z\simeq5$ IRX$-\beta$ relation.

Figures

Tables

Table 1. Table 1: The medians and standard deviations of the intrinsic UV stellar continuum slope ( β i subscript 𝛽 𝑖 \beta_{i} ), the stellar + nebular UV continuum slope ( β 𝛽 \beta ) and Δ ​ β Δ 𝛽 \Delta\beta ( = β − β i absent 𝛽 subscript 𝛽 𝑖 =\beta-\beta_{i} ) for all four SPS models (all in the case of zero dust). Δ ​ β Δ 𝛽 \Delta\beta is not the same as the difference in medians because the β 𝛽 \beta distributions are skewed.

SPS Model	${\beta }_i$	$\beta$	$\Delta \beta$
BPASSv2-100bin	$-2.60\pm 0.08$	$-2.42\pm 0.04$	0.17 $\pm$ 0.05
BPASSv2-100	$-2.55\pm 0.09$	$-2.42\pm 0.04$	0.14 $\pm$ 0.06
BPASSv2-300bin	$-2.63\pm 0.08$	$-2.39\pm 0.03$	0.23 $\pm$ 0.07
BPASSv2-300	$-2.58\pm 0.08$	$-2.39\pm 0.03$	0.19 $\pm$ 0.07

Table 2. Table 2: The best-fitting values of the reddening law slope ( n 𝑛 n ) for the TEx model, including both the stellar + nebular continuum and stellar continuum only cases.

SPS Model	Birth Cloud Age / Myr	$n_{\mathrm{neb}}$	${\chi }_{\mathrm{neb}}^2$	$n_{\mathrm{stellar}}$	${\chi }_{\mathrm{stellar}}^2$
BPASSv2-100bin	0	$-0.45$	0.33	$-0.60$	1.94
BPASSv2-100bin	1	$-0.47$	0.19	$-0.60$	3.95
BPASSv2-100bin	3	$-0.51$	0.49	$-0.62$	3.24
BPASSv2-100bin	5	$-0.60$	1.83	$-0.69$	4.58
BPASSv2-100bin	10	$-0.73$	12.84	$-0.86$	10.47
BPASSv2-100bin	15	$-0.53$	20.99	$-0.87$	17.57
BPASSv2-100	0	$-0.48$	0.12	$-0.58$	1.95
BPASSv2-100	1	$-0.52$	0.74	$-0.59$	3.80
BPASSv2-100	3	$-0.55$	1.54	$-0.60$	3.17
BPASSv2-100	5	$-0.67$	5.75	$-0.67$	4.09
BPASSv2-100	10	$-0.63$	13.66	$-0.75$	8.69
BPASSv2-100	15	$-0.33$	17.34	$-0.70$	15.00
BPASSv2-300bin	0	$-0.39$	1.04	$-0.58$	2.03
BPASSv2-300bin	1	$-0.43$	0.23	$-0.58$	3.37
BPASSv2-300bin	3	$-0.50$	0.47	$-0.64$	3.65
BPASSv2-300bin	5	$-0.64$	2.15	$-0.72$	5.10
BPASSv2-300bin	10	$-0.74$	17.17	$-0.93$	12.45
BPASSv2-300bin	15	$-0.45$	25.39	$-0.90$	21.55
BPASSv2-300	0	$-0.42$	0.42	$-0.56$	2.25
BPASSv2-300	1	$-0.47$	0.29	$-0.57$	3.38
BPASSv2-300	3	$-0.56$	1.87	$-0.63$	3.64
BPASSv2-300	5	$-0.70$	7.49	$-0.71$	4.95
BPASSv2-300	10	$-0.61$	17.59	$-0.80$	10.90
BPASSv2-300	15	$-0.25$	19.72	$-0.70$	18.94

Table 3. Table 3: The range of parameter values used in the fitting on the CF2000 model.

Parameter	Min	Max	Step Size
$n$	$-1.5$	$0.0$	0.05
$\mu$	$0.0$	$1.0$	0.10
$a_0$	$-4.5$	$-1.0$	0.25
$a_1$ (d$A_{1600}$/dM)	$0.1$	$0.8$	0.05
$t_{BC}$	$0.0$	$15$	1.00

Table 4. Table 4: The best-fitting n 𝑛 n values for the CF2000 model, including both the stellar + nebular continuum and stellar continuum only cases. The 1 σ 𝜎 \sigma confidence interval is obtained from solutions with Δ ​ χ 2 = 1 Δ superscript 𝜒 2 1 \Delta\chi^{2}=1 .

Parameter	$n_{neb}$	$\Delta n_{neb}$	${\chi }_{neb}^2$	$n_{stellar}$	$\Delta n_{stellar}$	${\chi }_{stellar}^2$
BPASSv2-100bin	$-0.35$		1.36	$-0.50$		0.95
BPASSv2-100	$-0.45$		1.53	$-0.45$		1.03
BPASSv2-300bin	$-0.40$		2.87	$-0.50$		0.81
BPASSv2-300	$-0.35$		2.84	$-0.50$		1.31

Table 5. Table 5: Recommend values for the CF2000 model, provided to allow a simple comparison with our model.

Parameter	Value
$n$	$-0.55$
$\mu$	$0.30$
$t_{BC}$	$12$ Myr

Equations

$$f_{\lambda }\propto {\lambda }^{\beta }$$

$$A(\lambda )=\{\begin{array}{ll}\infty & \text{if }t\le t_{BC},\\ \varphi {\lambda }^n& \text{if }t>t_{BC},\end{array}$$

$$\varphi =a_1{\log }_{10}(M/{\mathrm{M}}_{\odot })+{\mathrm{a}}_0.$$

$$A(t,\lambda )=\{\begin{array}{ll}\varphi {\lambda }^n& \text{if }t\le t_{BC},\\ \mu \varphi {\lambda }^n& \text{if }t>t_{BC},\end{array}$$

$$\mathrm{A}(\lambda ,\mathrm{M})={\mathrm{A}}_{1500}(\mathrm{M}){\left(\frac{\lambda }{1500}\right)}^{\mathrm{n}},$$

$$\Delta \beta ={\beta }_{\mathrm{obs}}-{\beta }_{\mathrm{int}}=1.297({\mathrm{A}}_{1268}-{\mathrm{A}}_{2580}),$$

$${\chi }^2=\sum _i{\left(\frac{\mathrm{model}(\mathrm{i})-\mathrm{data}(\mathrm{i})}{\sigma (i)}\right)}^2,$$

$${\mathrm{A}}_{1500}=0.5\times \log (M/{\mathrm{M}}_{\odot })-\mathrm{3.3.}$$

$$L_{\mathrm{IR}}={\int }_{912}^{\infty }{l}_{\lambda ,i}-l_{\lambda }d\lambda ,$$

$${\mathrm{IRX}}_{1600}\equiv \frac{{\mathrm{L}}_{\mathrm{IR}}}{{\mathrm{L}}_{1600}}=(10^{0.4{\mathrm{A}}_{1600}}-1)\frac{{\int }_{912}^{\infty }{\mathrm{l}}_{{\lambda }^{\prime },\mathrm{i}}\mathrm{d}{\lambda }^{\prime }}{{\mathrm{L}}_{1600,\mathrm{i}}}.$$

$${\mathrm{IRX}}_{1600}=1.75(10^{0.4{\mathrm{A}}_{1600}}-1),$$

$${\mathrm{A}}_{1600,\mathrm{Calzetti}}=2.09(\beta +{\beta }_{\mathrm{i}}),$$

$${\mathrm{A}}_{1600,\mathrm{SMC}}=0.99(\beta +{\beta }_{\mathrm{i}}).$$

$$A_{1600}=2.10(\beta +2.52),$$

$$A_{1600}=1.28(\beta +2.83),$$

$$A_{1600}=2.10_{-0.3}^{+1.8}(\beta +2.52).$$

Full text

††pagerange: The First Billion Years project: constraining the dust attenuation law of star-forming galaxies at $\mn@boldsymbol{z\simeq 5}$–The First Billion Years project: constraining the dust attenuation law of star-forming galaxies at $\mn@boldsymbol{z\simeq 5}$††pubyear: 2016

1 Introduction

Tracing the star-formation rate of galaxies across cosmic time, and particularly at $z\gtrsim 3$ (prior to the peak epoch of cosmic star formation), remains one of the most important focuses of observational cosmology, key to shedding light on the build up of stellar mass and metal content at early times, and clarifying the process of cosmic reionization. However, despite the substantial progress that has been made in constraining the star-formation rate density (SFRD) out to $z\approx 8$ (e.g. Madau & Dickinson, 2014), the lack of far-infrared (FIR) data for galaxies at $z\gtrsim 3$ means that, at these redshifts, star-formation rate estimates are mainly derived solely from rest-frame ultraviolet (UV) to optical observations. Therefore, our knowledge of the SFRD at $z\gtrsim 3$ is affected by inherent uncertainties regarding the correct dust corrections to apply at high redshifts.

When only rest-frame UV data are available, the observed UV luminosities must be dust-corrected based on correlations derived at lower redshifts. Most commonly, the Meurer et al. (1999) (M99) IRX-$\beta$ relation (IRX=$\nicefrac{{L_{IR}}}{{L_{UV}}}$) is used to estimate the intrinsic UV luminosity from the observed UV continuum slope $\beta$ ($f_{\lambda}\propto\lambda^{\beta}$). This relation was derived from a sample of local starburst galaxies and has been shown to apply for star-forming galaxies out to $z\simeq 2$ (e.g. Reddy et al., 2010; Reddy et al., 2012a). Under the assumption that the dust heating is dominated by a young stellar population, IRX can be related directly to the attenuation in the UV (AUV, see Meurer et al., 1999), allowing one to relate the absolute UV attenuation to the observed shape of the UV continuum slope, assuming the intrinsic shape is known. The form of the subsequent AUV-$\beta$ relation measured for local starbursts is what one would expect for a relatively grey reddening law, similar to that derived by Calzetti et al. (2000). Steeper laws such as the Small Magellanic Cloud (SMC) extinction curve (e.g. Gordon et al., 2003), or laws which include a strong $2175\rm{\AA}$ absorption feature (e.g., Milky Way; Cardelli et al., 1989) are not compatible with the local starburst galaxies. Therefore, when applying the Meurer et al. (1999) IRX-$\beta$ relation at $z\gtrsim 3$, one is implicitly assuming Calzetti-like reddening. Furthermore, in cases where rest-frame UV plus optical data is available, the intrinsic UV luminosities are typically derived via detailed SED-fitting (e.g. McLure et al., 2011; Bowler et al., 2015), and the subsequent UV dust-correction depends on which reddening laws are assumed for the fitting. Again, and motivated in part by its compatibility with the local IRX-$\beta$ relation, by far the most common assumption is the Calzetti et al. (2000) reddening law. Thus, this assumption is implicit in our current understanding of the UV luminosity function and SFRD evolution at high redshifts.

However, with the advent of the Atacama Large Millimeter Array (ALMA), the form of the attenuation law at high redshifts has begun to be questioned. Recent ALMA observations have indicated that the IR luminosities of ‘typical’ star-forming galaxies at $z\simeq 3$ are lower than expected assuming a Calzetti law, and have been interpreted as evidence that the attenuation curve is as steep, or steeper, than the SMC extinction curve at these redshifts (Capak et al., 2015; Bouwens et al., 2016). Similar results have also previously been reported for young galaxies ($<$ 100 Myr old) (e.g. Reddy et al., 2010; Reddy et al., 2012a) and sample of gravitationally lensed galaxies (e.g. Baker et al., 2001; Siana et al., 2008; Siana et al., 2009) at these redshifts. For a steeper, SMC-like, law a given amount of reddening in the UV is achieved with less absolute attenuation (i.e. the $A_{UV}-\beta$ relation becomes much shallower). If dust reddening at high redshift is SMC-like, then previous UV luminosity corrections made assuming the Meurer et al. (1999) IRX-$\beta$ relation, or Calzetti et al. (2000), will have overestimated the dust attenuation, and hence overestimated the intrinsic luminosities and SFRs at $z\simeq 3$. On the other hand, direct measurements of the attenuation curve at high redshifts have yet to find strong evidence for the steep SMC-like curves inferred from IR measurements. For example, Reddy et al. (2015) found a Calzetti-like attenuation curve shape at $\lambda<2600\rm{\AA}$ for a sample of N=224 star-forming galaxies at $1.4<z<2.6$. Using an independent method, Zeimann et al. (2015) found an average attenuation curve at $1.90<z<2.35$ similar to the Calzetti curve, and perhaps even greyer at lower masses, with a tentative detection of a mass dependence of the attenuation curve slope. Furthermore, Scoville et al. (2015) performed a similar analysis with a sample of N=266 galaxies at $2.0<z<6.5$ and again found no evidence for a deviation from a Calzetti-like attenuation curve. Interestingly, the Scoville et al. (2015) analysis is the only example of a direct attenuation curve measurement at the same redshifts probed by recent ALMA observations, and does not support the existence of a steep SMC-like curve.

In light of these recent results, the primary motivation of this paper is to use state-of-the-art hydrodynamical simulations of galaxies at $z=5$, from the First Billion Years project (FiBY) (Paardekooper et al., 2013, Khochfar et al. in prep), to constrain the form of the attenuation curve, by comparing the SEDs of simulated galaxies to the luminosity function (LF) and colour magnitude relation (MR) at $z=5$. This specific redshift was chosen primarily because it is the highest redshift at which there is good consistency between measurements of the LF and CMR across a variety of independent studies (e.g Bouwens et al., 2012; Rogers et al., 2014; Bowler et al., 2015; Bouwens et al., 2015); furthermore, it is at this approximate redshift that a number of recent ALMA observations have been made claiming evidence for an evolution of dust properties in the early Universe.

The structure of this paper is as follows: in Section 2 we describe how we extracted the relevant stellar data for galaxies in our sample from the FiBY simulation, and outline the method we employed to generate synthetic SEDs, including a description of SPS models and photoionization model assumptions. In Section 3 we review the intrinsic properties of our synthetic SEDs and describe how these properties are affected by our stellar population synthesis (SPS) and photoionization model choices. In Section 4 we describe the two dust models we have employed to constrain various properties of the attenuation law, before presenting the results of fitting our sample of simulated galaxies to the observed LF and CMR in Section 5. In Section 6 we discuss the IR properties of $z=5$ galaxies predicted by our models and compare these to recent ALMA observations. Finally, in Sections 7 and 8 we provide a discussion and summary of our main conclusions. The FiBY simulation adopts the following cosmological parameters, that we will also adopt throughout this paper: $\Omega_{n}=0.265$, $\Omega_{b}=0.0448$, $\Omega_{\Lambda}=0.735$, $H_{0}=71$ km s*-1* Mpc*-1* and $\sigma_{8}=0.81$.

2 Generating Synthetic SEDs

In this section we discuss how the synthetic SEDs were generated from the FiBY simulation data, including a detailed discussion of our SPS and photoionization modelling assumptions, and an overview of the physical properties of our final simulated galaxy sample.

2.1 SPS models

To construct the SED for each galaxy we used the latest release of ‘Binary Population and Spectral Synthesis’ (BPASSv2; Eldridge & Stanway, 2016; Stanway et al., 2016). The BPASSv2 models were chosen because, along with conventional single-star models, they also incorporate the evolution of massive stars in binary systems. Observations in the local Universe have shown that a potentially substantial fraction ($\gtrsim 70\%$) of massive stars undergo a binary interaction (such as mass transfer or a merger) during their lifetimes (e.g. Sana et al., 2012). Moreover, recent studies of galaxies at $z\simeq 2-3$ indicate that the BPASSv2 models are best able to replicate the observed nebular and stellar properties of star-forming galaxies at these redshifts (Steidel et al., 2016; Strom et al., 2017). We have checked that, for the purposes of this work, BPASSv2 single star models are equivalent to the more commonly adopted STARBURST99 models (Leitherer et al., 1999). In the context of this paper, the main effect of including massive binary evolution in SPS models is to generate harder ionizing spectra, particularly at low metallicities, and to boost UV output at older stellar ages relative to single-star populations.

We considered four sets of models which we refer to according to the upper mass cutoff of the initial mass function (IMF) and whether or not binary evolution is included. Each galaxy is assigned four separate SEDs corresponding to each of the four models. BPASSv2-100bin models include binary star evolution with an IMF cutoff of 100 M*⊙, and BPASSv2-100 are the equivalent single-star evolution models; similarly, BPASSv2-300bin models include binary evolution with an IMF cutoff of 300 M⊙, and BPASSv2-300 are the equivalent single-star evolution models. All models have an IMF index of $-1.3$ between $0.1-0.5$ M⊙* and $-2.35$ above 0.5 M*⊙* and cover the following metallicities: Z*∗* = 0.001, 0.002, 0.003, 0.004, 0.006, 0.008, 0.010, 0.014, 0.020, 0.030, 0.040.

2.2 FiBY simulation data

Here, we give a brief summary of the FiBY simulations but refer the reader to other papers for a more detailed description (e.g. Johnson et al., 2013; Paardekooper et al., 2015). The FiBY simulation suite is a set of high-resolution cosmological hydrodynamics simulations using a modified version of the GADGET code used in the Overwhelmingly Large Simulations (OWLS) project (Schaye et al., 2010). These simulations reproduce the stellar mass function and star formation rate of galaxies at $z\geq 6$ (Khochfar et al. in prep.), and at the same time also recover the right trends in the metallicity evolution of galaxies (Dalla Vecchia et al. in prep.). Furthermore, the simulations recover that low mass galaxies can reionize the Universe (Paardekooper et al., 2013, 2015), and that supernovae-feedback can drive the formation of cores in galaxies (Davis et al., 2014).

The code tracks metal pollution for 11 elements: H, He, C, N, O, Ne, Mg, Si, S, Ca and Fe and calculates the cooling of gas based on line-cooling in photoionization equilibrium for these elements (Wiersma et al., 2009) using tables pre-calculated with CLOUDY v07.02 (Ferland et al., 1998). Furthermore, the simulation incorporates full non-equilibrium primordial chemistry networks (Abel et al., 1997; Galli & Palla, 1998; Yoshida et al., 2006) including molecular cooling functions for H2 and HD. Star formation is modelled using the pressure law implementation of Schaye & Dalla Vecchia (2008), which yields results consistent with the Schmidt-Kennicutt law (Schmidt, 1959; Kennicutt, 1998). The threshold density above which stars form is set to $n=10$ cm*-3*. The simulations include feedback from stars is by injecting thermal energy in the neighboring particles (Dalla Vecchia & Schaye, 2012). For each Pop II supernova $10^{51}$ erg is injected once a star particle has reached an age of 30 Myr, corresponding to the maximum lifetime of stars that end their lives as core-collapse supernovae.

In this work we will focus on the FiBY_L and FiBY_XL simulations which cover co-moving volumes of (16 Mpc)3 and (32 Mpc)3 respectively ($\approx 3.7\times 10^{4}$ Mpc3 combined comoving volume). The individual gas and star particle masses in the simulations are log(M/M*⊙*) = 5.68 for the FiBY_XL and 4.81 for the FiBY_L. To build the spectral energy distributions of galaxies at $z=5$, halos at this redshift were identified in the simulation using the SUBFIND algorithm (Springel et al., 2001; Dolag et al., 2009). Then, for each halo, we extracted the masses, ages and metallicities of each star particle. The mass of each star particle is defined as the total initial mass as opposed to the current mass (i.e. at $z=5$). since the effects of mass-loss with stellar age are accounted for in the BPASSv2 stellar population synthesis models we use to generate the SEDs. The age of each star particle was calculated using the formation redshift and assumed cosmological parameters. Finally, we extracted the mass in metals for each star particle; for reasons which will be described in Section 2.2.1, we tracked both the total mass in metals and the mass in Fe. To ease computation time, we binned star particles with ages $>20$ Myr into bins of width 10 Myr, taking averages of mass, age and metallicity for each bin; star particles with ages $\leq 20$ Myr were stored individually.

To build a galaxy SED from a given model we looped over each star particle associated with the galaxy and assigned the stellar spectral energy distribution which best matched the age and metallicity of the particle. Since the BPASSv2 stellar SEDs are modelled as an instantaneous starburst with a total mass of 1 x 106 M*⊙*, we additionally scaled the flux to the mass of the star particle. We constructed the final galaxy SED by summing the SEDs of all the individual star particles.

For our final sample, we selected halos at $z=5$ with an absolute UV magnitude magnitude $\rm{M}_{1500}\leq-18.0$ where $\rm{M}_{1500}$ was calculated by integrating the rest-frame BPASSv2-100bin SED within a top-hat filter of width 100 $\rm{\AA}$ at a central wavelength of $1500\rm{\AA}$. This magnitude limit differs slightly according to the assumed SPS model, the limiting value of $\rm{M}_{1500}=-18.0$ for the BPASSv2-100bin model corresponds to $-17.89$, $-18.08$ and $-17.95$ for the BPASSv2-100, BPASSv2-300bin and BPASSv2-300 models respectively. However, changing the SPS model for which the magnitude cut is defined does not have a significant effect on the physical properties of our sample, or any of the results presented in this paper.

The final sample contains N=498 galaxies; the distributions of current mass (i.e. at $z=5$), star-formation rate (SFR) (averaged over the past 100 Myr), and specific star-formation rate (sSFR $\equiv$ SFR / M*∗) of our $z=5$ sample are shown in Fig. 1. The maximum galaxy mass in our sample is 1.4 x 1010* M*⊙, with a median value of 2.3 x 108 M⊙; the sample is mass complete down to 2.6 x 108 M⊙. The star-formation rates range from 0.3 - 60 M${}_{\odot}\rm{yr}^{-1}$ with a median of 1.1 M${}_{\odot}\rm{yr}^{-1}$. The median sSFR is 4.8 Gyr-1*.

2.2.1 Stellar metallicity

The BPASSv2 models assume solar element abundance ratios, however it is not obvious that the abundance ratios at $z=5$ should be solar; in particular, given the typical star-formation histories at $z=5$, there could be differences in the Fe/H and $\alpha$/H ratios relative to solar values.

The star-formation histories of the galaxies in the FiBY simulation are typically rising (Khochfar et. al. in prep), consistent with other observations and simulations at $z>2$ (e.g. Papovich et al., 2011; Finlator et al., 2011; Reddy et al., 2012b). For rising (or constant) star-formation histories the rate of enrichment of $\alpha$ elements from Type II (core-collapse) supernovae, which has a typical time scale of $\sim$ 10 Myr, will exceed that rate of Fe enrichment, since Fe is mainly returned to the ISM via SNe Ia, which have a typical timescale of $\sim$ 1 Gyr. Therefore, rapidly star-forming galaxies are expected to have enhanced $\alpha/\rm{Fe}$ ratios relative to the solar value. Indeed, it is not strictly necessary to appeal to star-formation histories to understand this offset at $z=5$, since the median age of the oldest star particle across our galaxy sample is 0.74 Gyr, with a maximum of 0.94 Gyr, implying that the majority of galaxies have not had time to undergo significant SNIa enrichment, regardless of their star-formation histories. As we discuss below, this $\alpha/\rm{Fe}$ enhancement can have a significant effect on the shape of the UV spectrum.

We were able to calculate the detailed abundance ratios of all galaxies in our sample because the mass fractions for 11 individual elements (H, He, C, N, O, Ne, Mg, Si, S, Ca and Fe) are tracked in the FiBY simulations, along with the total mass fraction in all metals. Fig. 2 shows the distributions of UV-weighted Fe metallicity and total metallicity of the star particles (hereafter referred to, for simplicity, as the $\alpha$ element metallicity), relative to solar, for the N=498 galaxies in our sample, where the solar mass fractions have been taken from Asplund et al. (2009). Consistent with the rising star-formation histories, the distributions in Fig. 2 (left-hand panel) show the difference between the $\alpha$ and Fe-based metallicities. The median $\alpha/\rm{Fe}$ ratio in our sample is $\simeq$ 3.9 times larger than the solar value, with the median Fe and $\alpha$ metallicities being ZFe/Z*⊙* = 0.14 and Zα/Z*⊙* = 0.53 respectively. The right-hand panel of Fig. 2 shows that, under both definitions, metallicity is an increasing function of stellar mass, however Zα/Z*⊙* shows a stronger evolution, implying that the $\alpha/$Fe ratio is also an increasing function of stellar mass.

Non-solar $\alpha/$Fe ratios have an important effect on the details of the emergent SED. Recently, Steidel et al. (2016) performed detailed stellar population synthesis and photoionization modelling of the composite rest-frame UV and optical spectra of 30 star-forming galaxies at $z=2.4$ with log(M/M*⊙)$>9.0$. They found that simultaneously fitting (i) the shape of the stellar spectrum in the UV, and (ii) the observed UV + optical emission line ratios, implied a difference of a factor $\simeq 4-5$ between the stellar and nebular metallicities (Z∗/Z${}_{\odot}\simeq 0.1$ and Zneb/Z${}_{\odot}\simeq 0.4-0.5$). This is consistent with the offset in our sample (Fig. 2); in fact, restricting our sample to galaxies above 109 M⊙* yields a median $\alpha/\rm{Fe}$ ratio of 4.2. As argued by Steidel et al. (2016), this offset is naturally explained if star-formation histories at these redshifts are typically rising (arguments for rising SFHs at these redshifts have been made in a number of studies e.g. Papovich et al., 2011; Reddy et al., 2012a), and by the fact that stellar opacity in the UV is dominated by Fe, whereas the nebular metallicity is indicative of the abundances of the important coolants in ionized gas (e.g O, N, C), and hence traces the $\alpha$ elements. To a first approximation, any differences between stellar and nebular metallicities is therefore plausibly a result of the same $\alpha$/Fe offset illustrated in Fig. 2; the lower stellar metallicity is driven by the sensitivity of the stellar UV spectrum to Fe/H, whereas the nebular metallicity follows $\alpha$/H. Taken together, all these considerations imply that, when applying SPS models with solar abundance sets to galaxies in which one suspects there will be an enhanced $\alpha/\rm{Fe}$ ratio, the stellar metallicity should be defined as Fe/H, rather than the mass fraction of all metals as is commonly assumed, at least when studying the properties of the UV spectrum. Ideally, one would use SSPs with non-solar abundance sets however, to our knowledge, there is no publicly available SPS code which accounts for both binary star evolution and also allows non-solar element abundances. As emphasized by Steidel et al. (2016), this should be a key goal for any future SPS models.

To summarize, as a consequence of the above considerations, we used the mass fraction in Fe, relative to solar, rather than the total mass fraction in metals, to select the appropriate metallicity BPASSv2 model when constructing the galaxy SEDs. The effect of this choice on the UV SED is illustrated in Fig. 3, which shows the average composite BPASSv2-100bin spectra of all N=73 galaxies with log(M/M*⊙*) $\geq 9.0$, assuming both $\rm{Z}_{*}=\rm{Z}_{\rm{Fe}}$ and $\rm{Z}_{*}=\rm{Z}_{\rm{Tot}}$, and will be discussed further in Section 3.

2.3 Photoionization models

Since the BPASSv2 models only account for the stellar component of each spectrum, we separately modelled the nebular contribution for each galaxy. The nebular continuum is a result of free-free, free-bound and two-photon emission from H ii regions in the galaxy. To model the nebular continuum and emission spectrum we use the photoionization code Cloudy (v13.03; Ferland et al., 2013). We ran Cloudy assuming a plane-parallel geometry and considered the following four key input parameters: (i) the shape of the incident ionizing radiation field, (ii) the metallicity of the nebular gas, (iii) the electron density and (iv) the ionization parameter. Of these, we directly set (i) and (ii) from the FiBY simulation data. For the shape of the incident radiation field we input the BPASSv2 stellar continuum of each galaxy, thus the nebular spectrum is modelled by assuming the entire galaxy acts as a single H ii region. Following the discussion in Section 2.2.1, we used the UV-weighted alpha-element stellar metallicity to set the nebular abundances, scaling to the Asplund et al. (2009) solar abundance set. This explicitly assumes that the gas phase metallicity traces the metallicity of the most recently formed stars. To set appropriate values of electron density and ionization parameter at $z=5$ we used recent near-IR spectroscopic observations of star-forming galaxies at $z\simeq 2-3$ which have illustrated that the typical values of these parameters are evolving from $z=0$ out to higher redshifts (e.g Steidel et al., 2014; Sanders et al., 2016; Cullen et al., 2016; Strom et al., 2017). For example, Strom et al. (2017) have presented a detailed analysis of $\rm{N}\sim 380$ galaxies at $z\sim 2.3$ and concluded that the typical value of electron density is $\simeq 300$ cm*-3* (representing, roughly, an order of magnitude increase over typical SDSS star-forming galaxies; see also Sanders et al., 2015), with ionization parameters in the range $-3.1<\rm{log}(\rm{U})<-2.5$, consistent with the best-fitting value of $\rm{log}(\rm{U})=-2.8$ from Steidel et al. (2016) and again representing an increase on typical SDSS values (log(U) $\simeq-3.2$; Liang et al., 2006; Kewley & Ellison, 2008).

We made the assumption that the average $z\sim 2.3$ values of these parameters are appropriate to use for our $z=5$ galaxies based on the following reasoning. Firstly, though the range of masses and SFRs of our sample are not representative of the galaxies from Strom et al. (2017) (which have medians of 1.0 x 1010 M*⊙* and 24 M${}_{\odot}\rm{yr}^{-1}$, factors $\sim$ 40 and 20 larger respectively) the typical values of sSFR are within a factor 2 (4.8 Gyr*-1* for our sample compared to 2.4 Gyr*-1*), and one of the key results presented in Strom et al. (2017) is that the sSFR is correlated with the degree of excitation and the evolution in optical line ratios they observe. Unfortunately, no explicit relations between sSFR and either log(U) or $n_{e}$ exist in the literature to our knowledge, and we do not attempt to second guess the potential evolution of these parameters at $z\simeq 5$. Secondly, the metallicities of our sample (both Fe/H and $\alpha$/H) follow their best-fitting values at $z\simeq 2-3$ in the relevant mass range. Given the similarities in sSFRs and metallicities, we made the assumption the other physical conditions in H ii regions at $z=5$ are also plausibly equivalent. Therefore, we adopted an electron density of $n_{e}=300$ cm*-3* and assumed log(U)$=-2.8$. We briefly discuss, in Section 7, how our results are affected under the assumption that physical conditions in H ii regions continue to evolve towards even more extreme values of these parameters at $z\gtrsim 3$.

The effect of adding the nebular continuum to the UV spectrum is illustrated in Fig. 4 for a typical galaxy in our sample (M*∗* = 2.8 x 108 M*⊙; SFR = 1.3 M${}_{\odot}\rm{yr}^{-1}$; sSFR = 4.6 Gyr-1*; ZFe/Z*⊙* = 0.06). In this particular example we show a BPASSv2-300bin model SED; the median nebular contribution over the wavelength range within which we measured $\beta$ slopes ($1268<\lambda<1950\rm{\AA}$) is $\simeq 6\%$, and the overall effect of the nebular continuum is to make the spectrum redder (i.e. increase the value of $\beta$). We note that although nebular emission is included in our photoionization modelling, nebular emission lines have no effect on our method of constraining the form of the attenuation law, therefore in Fig. 4 we do not show emission lines, since only the continuum is relevant in this context. A detailed discussion of the effect of the nebular continuum on $\beta$ across the whole sample is given in Section 3.

To summarize, we have modelled the nebular continuum of all galaxies in our sample by using the BPASSv2 stellar SEDs to define the shape of the incident spectrum, assumed the nebular metallicity to be equal to the UV-weighted total stellar metallicity of the galaxy, and used values for electron density and ionization parameter typical of star-forming galaxies at $z\simeq 2-3$. Therefore, for each galaxy, we have both the intrinsic stellar spectrum, and the combined stellar $+$ nebular spectrum. As we will discuss later in the paper, the pure stellar spectrum provides a useful limit for our study of plausible attenuation curves, since it represents the ‘maximally blue’ SEDs of galaxies at $z=5$ from the FiBY simulations, based purely on their star-formation histories and metallicities.

3 Intrinsic $\mn@boldsymbol{\beta}$ and M1500 Distributions

As will be discussed in detail in Sections 4 and 5, in order to constrain the form of the dust attenuation curve, we have compared our simulations to two observed properties of $z=5$ galaxies: (i) the UV continuum slope ($\beta$) and (ii) the volume density of galaxies as a function of their absolute magnitude at $\lambda=1500\rm{\AA}$ (M1500). In this section we describe the intrinsic values we measured for these parameters from the synthetic SEDs, as well as discussing how the various different assumptions we have made (e.g. stellar metallicity, nebular modelling) systematically affect these observables.

3.1 UV continuum slope ($\beta$)

The UV continuum slope of a galaxy is defined as a power-law expressed as

[TABLE]

where $\beta$ is the power-law index (e.g. Bouwens et al., 2010; Dunlop et al., 2012; Dunlop et al., 2013). In the absence of dust, the intrinsic UV continuum is sensitive to the age and metallicity of the stellar population, with older and more metal rich populations having redder continuum slopes.

3.1.1 Intrinsic stellar UV continuum

To measure $\beta$ for our synthetic SEDs we performed a linear fit to log$f_{\lambda}$ versus log$\lambda$ over the wavelength range $1268<\lambda<2580$, masking out stellar and interstellar absorption features using the windows defined in Calzetti et al. (1994), we refer to this as a ‘spectroscopic’ $\beta$ measurement. As we will describe in Section 5, we used $\beta$ slopes of $z=5$ galaxies reported in Rogers et al. (2014) as our observational comparison. These authors measured $\beta$ by fitting a power-law to five photometric bands which sampled the rest-frame UV spectrum in the wavelength range $\lambda\approx 1500-2600\rm{\AA}$. To ensure consistency, we checked that our spectroscopic $\beta$ measurements were consistent with the Rogers et al. (2014) photometric method. We found that both methods are consistent for idealized data (see Fig. 5), with the spectroscopic $\beta$ measurements systematically smaller by a median of $\Delta\beta=0.03$ across the whole sample. For our sample, we chose to use the spectroscopically-measured $\beta$ values rather than explicitly mimic observational methods.

Fig. 6 shows how the intrinsic stellar UV continuum ($\beta_{i}$) of our sample varies as a function of the UV-weighted age (defined as age weighted by the luminosity at $1500\rm{\AA}$ of each star particle) and the UV-weighted Fe-based metallicity (ZFe). We note that the data in Fig. 6 are taken from the SEDs generated with the BPASSv2-100bin models, however, although the normalization of $\beta_{i}$ differs slightly, the trends are the same across all SPS models. As expected, the value of $\beta_{i}$ increases with increasing UV-weighted age and UV-weighted metallicity. Moreover, we find that, at $z=5$, $\beta_{i}$ is more strongly correlated with the UV-weighted age within the range $10-80$ Myr probed by our sample, with the scatter in $\beta_{i}$, at fixed UV-weighted age, $\approx 1/3\times$ the scatter at fixed metallicity. $\beta_{i}$ also correlates positively with mass-weighted age, although with a larger dispersion due to the mass-weighted age probing older stellar populations which contribute less to the UV SED. For the BPASSv2-100bin model SEDs shown in Fig. 6, the median $\beta_{i}$ across all galaxies is $-2.60$, for the other SPS models the medians are $-2.55$, $-2.63$ and $-2.58$ for BPASSv2-100, BPASSv2-300bin and BPASSv2-300 respectively (see Table 1). Independent of SPS model, the dispersion of $\beta_{i}$ values is $\sigma_{\beta_{i}}\approx 0.08$.

Given the correlation between metallicity and $\beta_{i}$ shown in Fig. 6, it is clear that adopting Fe-based rather than total metallicities for generating stellar SEDs can have a strong impact on the shape of the UV continuum (see also Fig. 3). We find that the offset between $\beta_{i,\rm{Fe}}$ and $\beta_{i,\rm{Tot}}$ is an increasing function of stellar mass, consistent with the trend shown in the right-hand panel of Fig. 2. For the BPASSv2-100bin model the offset can be as large as $\delta\beta_{i}\approx 0.3$ at the highest masses, with a median of 0.13; the median offsets across all other models are 0.02, 0.23 and 0.09 for BPASSv2-100, BPASSv2-300bin and BPASSv2-300 respectively. This emphasizes the importance, in future work (including SED fitting to observations) of adopting the correct stellar metallicity when fitting the UV continuum, which, as we have previously discussed, is not equal to the total mass-weighted metallicity in the case of variable $\alpha$/Fe ratios.

3.1.2 Effect of the nebular continuum

As illustrated in Figs. 4 and 6, the nebular continuum has a non-negligible effect on the UV continuum shape of individual galaxies. Here we discuss the effect of the nebular continuum on $\beta_{i}$ across our full sample. For the remainder of this paper, $\beta$ refers to the UV continuum slope of the galaxies including the nebular continuum.

Fig. 7 shows the relationship between $\beta$ and $\beta_{i}$ across all four BPASSv2 SPS models. The addition of the nebular continuum acts to redden the intrinsic stellar slopes (i.e. leads to larger values of measured $\beta$). The effect is strongest for the galaxies with the steepest (or bluest) stellar UV continuum slopes, or, equivalently, the galaxies with the lowest metallicities and youngest UV-weighted ages. This is unsurprising since young, low metallicity stellar populations generate harder ionizing spectra and hence a stronger nebular continuum. For the bluest stellar continua the effect can be as large as $\Delta\beta=0.4$, but decreases to $\Delta\beta\approx 0.05$ at $\beta_{i}\approx-2.35$, with a medians of $\Delta\beta\approx 0.15-0.25$ across all models. In other words, the effect of the nebular continuum on $\beta$ is of the order $\approx 2-15\%$, with an average of $\approx 10\%$.

It is also evident from Fig. 7 that, while the difference between binary and single-star models is minor, $\Delta\beta$ is systematically larger for the BPASSv2-300 models. At its maximum, the magnitude of this effect represents a difference in $\Delta\beta$, at fixed $\beta_{i}$, of $\approx 0.05$; interestingly, this is large enough to mean that while the median $\beta_{i}$ is bluer for the BPASSv2-300 compared to the BPASSv2-100 models, the situation is reversed for the $\beta$ medians (see Table 1).

Finally, across our sample, including the nebular continuum has the effect of reducing the scatter in intrinsic $\beta$ values by a factor $\approx 2$ (see Table 1) and significantly flattens the relations between $\beta$ and metallicity and UV age (Fig. 6). The median $\beta$ across all SPS models is $\approx-2.40$, and we find that the minimum value of $\beta$ in the case of no dust attenuation at $z=5$ is $\beta_{\rm{min}}=-2.50$. This theoretical dust-free $\beta$ limit is consistent with the work of Dayal & Ferrara (2012) and Wilkins et al. (2013) at $z=6-8$ and, as we will return to later in the paper, significantly different to the value of $\beta_{\rm{min}}=-2.23$ assumed for the Meurer et al. (1999) IRX $-\beta$ relation.

3.2 UV continuum magnitude

The UV continuum magnitude (M1500,i) is primarily sensitive to the star-formation history of a galaxy over the last 10 - 100 Myr (Kennicutt & Evans, 2012). As mentioned in Section 2.2, we calculated the intrinsic absolute UV magnitude magnitude (M1500,i) by integrating the rest-frame SEDs over a top-hat filter of width $100\rm{\AA}$ at a central wavelength of $1500\rm{\AA}$.

Since our sample was defined by requiring an intrinsic stellar magnitude (i.e. dust free, no nebular contribution) M${}_{1500,i}\leq-18.0$, the faintest galaxies across all SPS models have roughly this magnitude (see Section 2.2); the (median, maximum) M1500,i across all models are $(-18.88,-23.25)$, $(-18.77,-23.16)$, $(-18.99,-23.36)$ and $(-18.87,-23.26)$ for the BPASSv2-100bin, BPASSv2-100, BPASSv2-300bin and BPASSv2-300 models respectively. In the $\approx 3.7\times 10^{4}$ Mpc3 volume, we do not find any galaxies intrinsically brighter than M${}_{1500,i}\approx-23.4$.

As with the $\beta$ values, the choice of stellar metallicity also affects the intrinsic stellar M1500,i distributions, in the sense that a galaxy with a low SPS model metallicity will be brighter in the UV, due to reduced photospheric line blanketing. Again the effect is strongest for the binary models and models with a higher mass IMF cutoff, but not strong enough to affect the results presented in this paper; for the BPASSv2-300 models, the median offset in M1500,i is 0.23 mag. Finally, the effect of the nebular continuum on M1500,i is also small; across all models, adding nebular continuum increases the brightness of the SEDs by $\lesssim 0.1$ mag.

4 Modelling the dust attenuation

Having defined the intrinsic properties of simulated galaxies at $z=5$, we now turn to exploring models for dust attenuation, which should map these intrinsic properties to their observed values. In particular, we are interested in whether deviations from a Calzetti-like attenuation curve are required to achieve this match, since evidence (both observed and simulated) for a steeper attenuation curve, more like the SMC extinction curve, at high redshifts has been reported in the recent literature (e.g. Capak et al., 2015; Bouwens et al., 2016; Mancini et al., 2016, see further discussion below). To do this, we considered two simple models for the dust attenuation: (i) a total birth-cloud extinction model (TEx) and (ii) a Charlot & Fall (2000)-like model (CF2000) which we describe in detail below.

4.1 TEx

The total birth-cloud extinction model (TEx) assumes that, below a critical age, star particles are completely embedded in their birth clouds and suffer total extinction, while at larger ages either their birth clouds have been photo-evaporated by OB star associations, or the OB associations have drifted away or been blown-out of the birth-cloud, and hence the star particles suffer only the global ambient ISM attenuation of the galaxy.

Local estimates of the timescale for molecular cloud destruction are of the order $\sim$ 10 - 30 Myr (Blitz & Shu, 1980; Blitz et al., 2007; McKee & Ostriker, 2007), values which are in good agreement with current predictions from hydrodynamical simulations (e.g. Colín et al., 2013; Körtgen et al., 2016); while estimates of the timescales for the ejection of OB associations from their birth clouds are of the order 1 - 3 Myr (e.g. Israel, 1978; Leisawitz & Hauser, 1988), consistent with the youngest ages of galactic OB associations (e.g. Massey et al., 1995).

Given these observations, we chose to model the following six birth cloud age thresholds $t_{BC}=$ 0, 1, 3, 5, 10 and 15 Myr. The upper limit of 15 Myr is chosen to be representative of the timescale for molecular cloud dispersion, and to mirror the work of Jimenez et al. (2000) who find, using a similar model, that applying this age threshold successfully reproduces the shape of the UV spectrum of the $z=2.73$ Lyman-break galaxy 1512-cB58 (Pettini et al., 2000). The reddening (in magnitudes) as a function of wavelength in this model can be written as,

[TABLE]

where $\phi$ is a linear function in log(M/M*⊙*) given by:

[TABLE]

In this prescription, the overall normalization (or ‘amount’ of dust) is a linearly increasing function of the logarithm of stellar mass. It is well established in the local Universe that effective dust optical depth increases with stellar mass (e.g. Garn & Best, 2010); we adopt a linear model, both for simplicity, and because there is no strong evidence at high redshift to favour a more complex functional form. We modelled the wavelength dependence of the attenuation as a power-law with exponent $n$; this simple parameterization is a good approximation to many observationally derived attenuation curves which don’t include a $2175\rm{\AA}$ ‘bump’ as illustrated in Fig. 8. It is worth mentioning that the results of this paper are not dependent of the presence or size of the $2175\rm{\AA}$ ‘bump’ feature, since the Calzetti et al. (1994) windows are selected to mask out the $2175\rm{\AA}$ region, and Rogers et al. (2014) have demonstrated that their observationally derived $\beta$ values are also insensitive to this feature.

4.2 CF2000

The CF2000 model follows the prescription of Charlot & Fall (2000), with the addition of a mass-dependent normalization ($\phi$) as above. In this model the wavelength-dependent optical depth experienced by a star particle of age (t) is given by,

[TABLE]

where $\mu$ defines the fraction of the total dust optical depth contributed by the ambient ISM, $n$ is the exponent of the reddening curve, $\phi$ is the overall normalization and $t_{BC}$ is the age of transition between birth cloud and ambient ISM optical depth. Again we parameterized the normalization (or ‘amount’ of dust) $\phi$ as a linear function of the logarithm of stellar mass as in Eq. 3.

The difference here is that, instead of suffering total extinction, stars within their birth cloud suffer additional attenuation with respect to stars in the ambient ISM of a galaxy. Using a sample of local starburst galaxies, Charlot & Fall (2000) estimate $\mu\approx 1/3$, $n\approx-0.7$ and $t_{BC}\approx 10$ Myr, however, as will be described in more detail below, we kept all three parameters free when attempting to fit our simulation data to the observed LF and CMR.

5 Fitting to the LF and CMR

In this section we present the results of using these two dust models to fit to the observed UV galaxy luminosity function and CMR at $z=5$. The focus of this section is to constrain the dust model parameters (specifically the reddening slope $n$) before, in Section 6.2, we discuss how the predictions of our best-fitting dust models compare to the observed IR properties of high-redshift galaxies.

5.1 TEx

For the TEx model, the mass normalization $\phi$ was derived by directly fitting to the $z=5$ luminosity function of Bowler et al. (2015). The method is illustrated in Fig. 9: for each birth-cloud age and SPS model, the sample was split into six bins of M1500, the number density of galaxies in each bin was determined, and the value of A1500 required to match the luminosity function calculated. A value of A1500 = 0 was assigned to bins in which the number density fell below the observed relation.

We were then able to fit the relationship between A1500 and the median stellar mass in each bin. The inset panels of Fig. 9 show linear fits to A1500 versus stellar mass (bins with A1500 = 0 were excluded from the fitting). Using this relation we can then rewrite Equation 2 for $t>t_{\rm{BC}}$ as,

[TABLE]

where M $\equiv$ log(M/M*⊙*), leaving the slope of the attenuation law ($n$) as the only model parameter to be fitted.

To constrain the slope of the attenuation law, the CMR of our sample was compared to the observed CMR of Rogers et al. (2014) by varying the value of $n$, applying Equation 5, and re-calculating $\beta$ and M1500. For a giving mass normalization, the change in $\beta$ can be expressed as,

[TABLE]

where $A_{1268}$ and $A_{2580}$ are the attenuations at the wavelengths limits of the Calzetti et al. (1994) windows.

The value of $n$ was varied over the range $-1.50<n<0.00$, in steps of 0.01, to encompass the range of $n$ values consistent with observed attenuation/extinction curves (see Fig. 8). For each value of $n$, the sample was binned into six bins of M1500, and the total $\chi^{2}$ was then computed,

[TABLE]

where the summation is over all bins and $\sigma(i)=0.1$ is taken as the intrinsic $\beta$ scatter at the faint end of the CMR estimated by Rogers et al. (2014), and is representative of the $\beta$ distributions given in Table 1. The results of this fitting procedure are illustrated in Fig. 10 and the best-fitting $n$ values and their corresponding $\chi^{2}$ values are listed in Table 2. Figs. 9 and 10 show data in which the nebular continuum correction as described in 3.1.2 has been applied; however, for completeness, in Table 2 we also report the best-fitting values without applying the nebular correction, this provides a useful comparison for the limiting case of no nebular contribution (i.e. a 100$\%$ escape fraction, and the bluest possible UV slopes).

The results show that, across all SPS models, TEx dust models with birth-cloud ages $>5$ Myr, are ruled out at $>2\sigma$ confidence ($\chi^{2}>\chi^{2}_{min}+4$). Combining this with the luminosity functions shown in Fig. 9, we can further constrain the birth-cloud age to be $\lesssim 3$ Myr, since at ages $\geq 3$ Myr the intrinsic volume density of galaxies at $M_{1500}>-19.0$ falls below the Bowler et al. (2015) $z=5$ LF. Formally, across all SPS models, the best-fitting solutions corresponds to either the 0 or 1 Myr models, with an average best-fitting values of $n_{\rm{neb}}=-0.46$ for the nebular corrected spectra, and $n_{\rm{stellar}}=-0.59$ for the stellar continuum only models. In general, the nebular-corrected spectra require greyer reddening slopes because the UV continuum slopes are already reddened with respect to the intrinsic stellar spectrum; however, this distinction gradually disappears at $t_{BC}\geq 10$, as the strength of the nebular continuum decreases, and the best-fitting $n$ values converge. Finally, we find no distinction, even at the 1$\sigma$ level ($\Delta\chi^{2}=1$), between the different SPS models. Taken all together, the results suggest a birth cloud age in the range $0\leq t_{BC}<3$, and a attenuation curve slope in the range $-0.6\leq n\leq-0.4$.

It is clear that, for the TEx model, ‘greyer’ reddening laws with Calzetti-like slopes are favoured over steeper SMC-like slopes. From Table 2 it can be seen that, in general, the steepness of the best-fitting reddening slope increases with increasing birth cloud age, however even at $t_{BC}=10-15$ Myr the values are in the range $-0.8\lesssim n\lesssim-1.0$, still below the SMC extinction-law value of $n_{\rm{SMC}}\approx-1.24$. The increase in $n$ with $t_{BC}$ occurs because, as the birth cloud age increases, the intrinsic SEDs become fainter in $M_{1500}$, thus a smaller $A_{1500}$ is required to fit to the luminosity function, and a steeper $n$ is required to simultaneously redden the SEDs to match the CMR. Although this effect is partially compensated for by the intrinsic SEDs becoming redder at larger $t_{BC}$, the compensation is not sufficient to keep the reddening law as grey as a Calzetti law. It is possible that increasing the birth cloud age further would lead to SMC-like values of $n$, however this is really a moot point since, at $t_{BC}\geq 3$ Myr, the data become incompatible with the observed LF. Indeed, an acceptable match to the observations can only be achieved if the UV output from stars $\lesssim$ 3 Myr old is visible in the integrated spectra of galaxies; therefore, OB star associations should be able to drift away from, or photo-evaporate, their parent molecular clouds within this timescale. Such a scenario is consistent with the youngest ages of OB star associations observed in the Milky Way (e.g. Massey et al., 1995).

In summary, comparing the results of the TEx model to the observed LF and CMR at $z=5$ leads us to conclude that: (1) no satisfactory match is achieved unless the UV output from stars $\lesssim$ 3 Myr old is visible in the integrated UV spectra of galaxies; (2) the models show no evidence for a departure from Calzetti-like attenuation laws; and (3) there is no strong preference for binary or single-star BPASSv2 models.

5.2 CF2000

To fit the CF2000 model we applied a more detailed analysis, performing a $\chi^{2}$ minimization across all parameters as listed in Table 3. This allowed us to get a better estimate of the formal error on the attenuation curve slope by marginalizing of all free parameters. The ranges of parameter values were chosen after some experimentation using coarser step sizes. For each combination of parameter value, we calculated the optical depth as a function of wavelength for the birth cloud ($t\leq t_{BC}$) and ambient ISM ($t>t_{BC}$) as given by Equation 4, and reddened the spectrum accordingly. We binned the reddened sample into six bins of M1500 and performed a $\chi^{2}$ fit to both the CMR and LF. For the fit to the luminosity function, we defined the error in each bin as the Poisson error on the number counts in that bin. The best-fitting model parameters were determined via $\chi^{2}$ minimization, where the $\chi^{2}$ value for each set of model parameters was taken as the sum of the CMR and LF values (i.e. $\chi^{2}=\chi_{\rm{LF}}^{2}+\chi_{\rm{CMR}}^{2}$). From the resulting $\chi^{2}$ we constructed marginalized probability density functions for each parameter.

Fig. 11 shows the LF and and CMR for the model parameters at the minimum $\chi^{2}$ (see Table 4 for the corresponding parameter values). The figure shows the case including nebular continuum emission but the values for the case of stellar-only SEDs (i.e. maximally blue) are also given in the table. From Table 4 it can be seen that acceptable fits can be achieved for both stellar-only and stellar $+$ nebular models. The insets in Fig. 11 show the normalized probability density functions for the slope of the reddening curve, $p(n)$. Reassuringly, the $p(n)$ curves are consistent with the formally best-fitting $n$ values from our $\chi^{2}$ grid. The formally best-fitting reddening slopes are again slightly greyer than the Calzetti law for models including a nebular continuum; however, within the estimated $1\sigma$ uncertainties of all models, the range of acceptable $n$ values is $-0.7\leq n\leq-0.3$, consistent with both Calzetti et al. (2000) and Reddy et al. (2015), but incompatible with a steep reddening law resembling the SMC extinction curve. Furthermore, the results are fully consistent with the acceptable range of $n$ values ($-0.6\leq n\leq-0.4$) obtained from the simpler TEx model.

Overall, across the TEx and CF2000 models, the majority of formal best-fitting values of $n$ are marginally greyer than a Calzetti law. This is not unprecedented, since evidence for reddening curves shallower than Calzetti has been found in some direct observational studies of normal star forming galaxies at $z\simeq 2$ (e.g. Zeimann et al., 2015). However, given the SPS and nebular continuum modelling uncertainties, it is not possible to say definitively that a shallower law is favoured by our data. The combined analysis does, however, definitively rule out attenuation curves steeper than $n\lesssim-0.7(-0.9)$ at $1\sigma(2\sigma)$ (see Fig 12). Crucially, we find no evidence for an attenuation curve in the UV as steep as the SMC extinction curve ($n=-1.24$), we will return to this issue in Section 7.

5.2.1 Other CF2000 parameters

Although the focus of this analysis is on the slope of the reddening curve, the constraints placed on the other CF2000 model parameters were also investigated. We found best-fitting $\mu$ values of either 0.3 or 0.4 depending on the SPS and nebular assumptions, and values within $1\sigma$ uncertainties in the range $0.1<\mu<0.8$. This is consistent with the original Charlot & Fall (2000) results of $\mu\approx 1/3$ with a large scatter. For the birth cloud age ($t_{BC}$) we find formally best-fitting values in the range $11-13$ Myr, however statistically-acceptable solutions at any value of $t_{BC}$ in the full parameter range ($0-15$ Myr) i.e. the birth-cloud age is not constrained by fitting to the LF and CMR within the CF2000 parameterization. Finally, we find values of optical depth versus mass slopes (d$\rm{A}_{1600}$/dM) in the range $0.2<\rm{d}\rm{A}_{1600}/\rm{dM}<0.7$, with a best-fitting value of $\rm{d}\rm{A}_{1600}/\rm{dM}\simeq 0.5$. Of course, there is a degeneracy between these parameters since they all directly affect the change in $M_{1500}$ which is mainly constrained by fitting to the LF.

5.3 Recommended dust model parameters

Despite the relatively large range of parameter values over which acceptable can be found, we find that a Charlot & Fall (2000) type model, with parameter values not too dissimilar from the original values quoted, provides a simple and adequate framework for describing the dust attenuation law at $z\simeq 5$ from the FiBY simulation. In the interest of future work, we provide a list of our recommend values based on the above analysis in Table 5. Given that a Calzetti-like curve is fully consistent within the $1\sigma$ uncertainties of the model, we recommend the power-law index $n=-0.55$ for consistency with this commonly adopted attenuation curve. However, as this slope is only fit to the Calzetti curve at $\lambda<2600\rm{\AA}$ (Fig. 8), for attenuation curves extending from UV to optical we recommend using the Calzetti et al. (2000) relation explicitly or, given its similarity, the Reddy et al. (2015) curve. In addition to the values listed in Table 5, we recommend the following relation between $\rm{A}_{1500}$ and log(M/M*⊙*), which provides the normalization of the attenuation curve:

[TABLE]

This relation implies $\rm{A}_{1500}=0$ for $\rm{log(M/M_{\odot})}\leq 6.6$, and that, for an order of magnitude increase in mass, dust attenuation at $1500\rm{\AA}$ increases by 0.5 mags.

6 Dust model predictions

In this section we discuss the infrared (IR) predictions of our best-fitting dust models, and how these compare to the latest ALMA observations at $z\simeq 5$. Since we have derived the dust attenuation models by comparing the simulated SEDs to observed UV properties, it is useful to ask whether these same models are consistent with the, admittedly sparse, IR data at the same redshift, which trace the dust emission. In addition, we can use our models to make predictions for future IR surveys. The total IR luminosity (LIR) of each SED was calculated by applying the best-fitting dust model, subtracting the reddened SED from the intrinsic SED, and integrating the difference across all wavelengths i.e.:

[TABLE]

where $l_{\lambda,i}$ is the intrinsic SED, and $l_{\lambda}$ is the reddened SED. This integral converges at an upper wavelength limit of $\lambda\approx 10,000\rm{\AA}$. Throughout this section we include the IR predictions of each SPS model for both the best-fitting TEx and CF2000 dust model, and for both the nebular continuum and stellar continuum only cases (i.e. a total of 16 models encompassing the various SED and dust prescriptions).

6.1 ALMA HUDF $z=5$ source

Recently, Dunlop et al. (2017) (D17) have presented a 1.3-mm ALMA mosaic covering the full 4.5 arcmin2 of the Hubble Ultra-Deep Field (HUDF). Their final sample consists of 16 sources with point-source flux densities $S_{1.3}>120\mu$Jy at a mean redshift $\langle z\rangle=2.15$. Interestingly, one source in their sample (UDF12) has a spectroscopic redshift of $z=5.00$ (Hathi et al., 2008), and therefore provides a useful comparison to the number counts expected over a similar area from our simulations. UDF12 has a stellar mass log(M/M*⊙*) $=9.6\pm 0.12$, and a sSFR $=9.29\pm 4.80$, placing it in the high mass and high sSFR tail of our sample (see Fig. 1).

6.1.1 1.3mm number counts

To derive $S_{1.3}$, first LIR was converted into a SFRIR using the same relation adopted by D17 (from Murphy et al., 2011), and $S_{1.3}$ was calculated using the conversion given in D17: $S_{1.3}$ (in $\mu$Jy) $\approx 3.3\times$ SFRIR (in $M_{\odot}$yr*-1*). The number counts as a function of $S_{1.3}$ are shown by the histograms in Fig. 13. Depending on model assumptions, we find that between N $\approx 3-10$ galaxies have flux densities above the threshold $S_{1.3}>120\mu$Jy within the $\approx 3.7\times 10^{4}$ Mpc3 volume of our simulation, or equivalently we find a volume density of N($S_{1.3}>120\mu$Jy) $\approx(1-3)\times 10^{-4}$ Mpc*-3*. The number counts are larger for the models which include a nebular contribution due to the additional luminosity provided by the nebular continuum and emission lines in these SEDs, leading to an increase in the output LIR.

Over the 4.5 arcmin2 HUDF area, the comoving volume between $4.5<z<5.5$ is $\approx 1.3\times 10^{4}$ Mpc3, therefore, to a rough approximation, the simulation volume is $\approx 3\times$ the HUDF volume at $z=5$. Scaling our number counts by this factor results in a prediction of N $\approx 1-3$ galaxies with $S_{1.3}>120\mu$Jy within the HUDF ALMA mosaic, consistent with the observations. Clearly, a comparison to number counts over such a small area of sky is hindered by cosmic variance, nevertheless these numbers provide some reassurance that we are not seriously under or over-predicting infrared flux densities.

6.1.2 1.3mm flux density versus UV continuum slope

In Fig. 14 we directly compare observed quantities tracing the absolute attenuation ($S_{1.3}$) and dust law ($\beta$) of UDF12. UDF12 has a 1.3mm flux density of $S_{1.3}=154\pm 40$ (Dunlop et al., 2017) and a UV continuum slope of $\beta=-1.70\pm 0.12$, where $\beta$ was measured by fitting a power-law to the $I_{814}$, $z_{850}$, $Y_{105}$, $J_{125}$ and $H_{160}$ photometric data points, consistent with the value obtained using a standard SED fitting procedure, in which the intrinsic SEDs are modelled as a power-laws with varying slopes ($\beta=-1.71$).

Fig. 14 shows that the $S_{1.3}$ and $\beta$ values extracted from the best-fitting dust models are consistent with the observation. Although only $\approx 1\%$ of our sample have comparable flux densities ($S_{1.3}\geq 120$ $\mu$Jy), all fall within $<2\sigma$ of UDF12 and the majority fall within $<1\sigma$. Fig. 14 also demonstrates that, for the majority of models, the flux densities of the $S_{1.3}\geq 120$ $\mu$Jy sources are biased high with respect to the average $S_{1.3}$ of galaxies with similar $\beta$ slopes. There is no significant difference between the TEx and CF2000 models, or between the various SPS assumptions; UDF12 is consistent with the $1\sigma$ scatter on the running averages of all models.

In summary, the simple best-fitting dust models with a Calzetti-like reddening law slope ($-0.7\leq n\leq-0.3$) derived from UV observations are simultaneously compatible with the infrared number counts and properties of sources in the ALMA HUDF image. In other words, for at least this $z\simeq 5$ ALMA source, the models account for the both the absolute attenuation and reddening across the UV-IR wavelength baseline. However, although this is clearly encouraging, larger sample sizes at $z>5$ are needed in order to perform a more robust comparison.

6.2 The $\mn@boldsymbol{z=5}$ IRX-$\mn@boldsymbol{\beta}$ relation

As discussed in the introduction, the Meurer et al. (1999) (M99) IRX-$\beta$ relation, and the underlying, implicit, AUV-$\beta$ relation, has been commonly used to dust correct UV data at high redshifts, and has recently been used to argue that reddening at $z>5$ follows a SMC extinction-like curve, rather than a Calzetti curve (Capak et al., 2015; Bouwens et al., 2016). In this section we explore the IRX-$\beta$, and underlying AUV-$\beta$, relations implied by our best-fitting models.

6.2.1 IRX $-\beta$ relation

Meurer et al. (1999) define IRX as the ratio of far infrared flux ($F_{\rm{FIR}}$) to flux at 1600$\rm{\AA}$ ($F_{1600}$), where $F_{1600}$ is defined as a generalized flux of the form $F_{\lambda}=\lambda f_{\lambda}$, and $f_{\lambda}$ is the flux density per wavelength interval. In their derivation of IRX1600, they include a term (BC(FIR)Dust) to convert the total bolometric IR flux to a FIR flux. In this analysis we work with luminosities, rather than fluxes, and ignore this BC(FIR)Dust correction term since we can calculated the total bolometric IR luminosities of our SEDs. In this case IRX1600 can be written as

[TABLE]

Here the first term gives the fraction of L1600,i absorbed by dust and $A_{1600}$ is the attenuation at $1600\rm{\AA}$. For a given dust model, $A_{1600}$ was calculated by differencing the magnitudes at $1600\rm{\AA}$ of the intrinsic SED and the attenuated SED, where the magnitudes were calculated using a top-hat filter centered on $1600\rm{\AA}$ of width 100$\rm{\AA}$. The second term is the ratio of the maximum amount of of luminosity available for dust heating ($l_{\lambda^{\prime},i}$ is the unattenuated luminosity of the SED) divided by the intrinsic $1600\rm{\AA}$ luminosity.

Assuming the second term is a constant, this equation implies a direct transformation between the attenuation at $1600\rm{\AA}$ and IRX. Then, since for a given attenuation curve the attenuation at $1600\rm{\AA}$ is related to the observed $\beta$ slope, the position of galaxies in the IRX-$\beta$ plane can be used to infer the properties of the attenuation curve (e.g. Meurer et al., 1999), although one must be wary of the crucial dependence on the intrinsic $\beta$ slopes as we will discuss below.

In Meurer et al. (1999), the second term in the equation, denoted BC(1600)∗, was calculated using 100 Myr duration constant star-formation rate SPS models from Leitherer & Heckman (1995). They found $1.56\leq\rm{BC}(1600)_{*}\leq 1.76$, adopting a final value $\rm{BC}(1600)_{*}=1.66\pm 0.15$. Other authors have adopted the value of $1.75\pm 0.25$ derived later by Calzetti et al. (2000) (e.g. Reddy et al., 2006; Bouwens et al., 2016). It is worth comparing this to the BC(1600)∗ values of our SEDs, across all SPS models. For the range of BPASSv2 models adopted in this paper, the median BC(1600)∗ across all SEDs was calculated for the stellar + nebular and stellar-only cases. For the models including nebular continuum we find $1.70\leq\rm{BC}(1600)_{*}\leq 1.78$, and for stellar continuum only models we find $1.42\leq\rm{BC}(1600)_{*}\leq 1.52$. The full range of values is consistent with original derivation of Meurer et al. (1999) but emphasizes how large an effect the nebular continuum can have. Interestingly, the $\rm{BC}(1600)_{*}$ values derived from SEDs including nebular continuum are in much better agreement with previous estimates (e.g. Calzetti et al., 2000). For the remainder of this paper we adopt $\rm{BC}(1600)_{*}=1.75$. The IRX equation then simply becomes:

[TABLE]

and the relationship between IRX and $\beta$ becomes dependent only on the A${}_{1600}-\beta$ relation.

6.2.2 A${}_{1600}-\beta$ relation

If one assumes that all sources have the same intrinsic UV slopes (i.e. a constant $\beta_{i}$), then for a given reddening law slope ($n$) and a simple homogeneous star/dust mixture geometry (e.g. similar to the TEx model), there is a simple mapping between A1600 and $\beta$ given by combining equations 5 and 6. For a Calzetti-like reddening law ($n=-0.55$) this becomes,

[TABLE]

whereas for an SMC-like law ($n=-1.24$) the slope of the relation is much shallower,

[TABLE]

However, if there is an underlying correlation between $\beta_{i}$ and A1600 (i.e. $\beta_{i}$ is not assumed to be constant across the population), then these relations do not hold. Indeed, we find that more massive galaxies, with older UV-weighted ages, have larger values of both $\beta_{i}$ and A1600, such that there is a general $\beta_{i}-\rm{A}_{1600}$ correlation in our sample. This correlation is much stronger for the stellar-only models, since adding nebular continuum significantly reduces the range of $\beta_{i}$ values (see Table 1). The situation is further complicated for ‘clumpy’ geometries like CF2000, since for these models Equation 5 does not apply, and so simple relations like equations 12 and 13 above cannot be derived. The upshot of all this is that the slope of the A${}_{1600}-\beta$ relation, derived from our models, does not simply map to the power-law slope of the reddening law as described above. Furthermore, the values of $\beta_{i}$ in the resulting fits are representative of the minimum values of $\beta_{i}$ of the SEDs, rather than the median value across the population. Nevertheless, these equations are good approximations in the case of a small scatter in intrinsic UV slopes.

6.2.3 $\rm{IRX}-\beta$ for the best-fitting dust models

We derived A${}_{1600}-\beta$ and $\rm{IRX}-\beta$ relations for our sample for both the stellar + nebular and stellar-only cases; the results are shown in Figs. 15 and 16. For each case we performed a linear least-squares fit to the best-fitting $\beta$ and A1600 values from our models (combining the data across both dust models and all SPS models), and substituted the resulting fit into Equation 11 to give the corresponding IRX $-\beta$ curve. The best-fitting A${}_{1600}-\beta$ relations are as follows (see Fig. 15), for the stellar + nebular case:

[TABLE]

and for the stellar-only case:

[TABLE]

In both cases, the scatter about the relations is $\sigma=0.14$.

The relation for the stellar + nebular case is virtually identical to assuming a Calzetti attenuation law for sources with an intrinsic $\beta_{i}=-2.52$ (see Equation 12). In the stellar-only case, the slope of the relation is significantly shallower, this is a result of both the steeper best-fitting attenuation curves, and the stronger correlation between $\beta_{i}$ and best-fitting $A_{1600}$. The effect of the nebular continuum is substantial, up to $\Delta$A${}_{1600}\approx 0.5$ at the extreme ends of $\beta$ values. The scatter of $\sigma=0.14$ around the $\rm{A}_{1600}-\beta$ relation derived above only represents the scatter about the best fitting models, to get a more robust estimate from our models we estimated $\rm{A}_{1600}-\beta$ based on the $1\sigma$ limits on $n$ ($-0.7<n<-0.3$; see Fig. 12 and Table 3). Taking these limit into account, we estimate the following best-fitting $\rm{A}_{1600}-\beta$ relation for the stellar + nebular case,

[TABLE]

In reality, a nebular continuum contribution is expected in star-forming galaxies and, though it will likely vary for individual galaxies, we recommend Equation 16 as the more reliable estimate of the average A${}_{1600}-\beta$ relationship at $z=5$. Futhermore, this equation is reassuringly similar to observational constrainsts based on of samples of Lyman-beak galaxies at $z\sim 3-6$, especially in terms of the intrinsic UV slopes at these redshits (e.g. de Barros et al., 2014; Castellano et al., 2014). One immediate consequence of this is that dust corrections to the UV, based on the Meurer et al. (1999) IRX $-\beta$ relation, will underestimate A1600 for an observed $\beta$, by a median of $\approx 0.6$ across the range of observed $\beta$ values $-2.5<\beta<-1.6$, equivalent to an underestimation of the resulting star-formation rate of a factor $\approx 2$.

Finally, in Fig. 16 we show the predicted position of $z\simeq 5$ galaxies in the $\rm{IRX}-\beta$ plane along with a variety of data from the literature at $z>3$ (taken from: Coppin et al., 2015; Dunlop et al., 2017; Schaerer et al., 2015; Capak et al., 2015; Watson et al., 2015; Bouwens et al., 2016). It can be seen that the best-fitting line (using Eq. 16) is roughly consistent with the shape of the original Meurer et al. (1999) curve with a shift to bluer intrinsic UV slopes of $\Delta\beta_{i}\simeq-0.2$. Again, this is fully consistent with typical $z\simeq 5$ star-forming galaxies experiencing, on average, Calzetti-like attenuation, similar to local star-forming galaxies, yet having bluer intrinsic UV slopes based on the typically rising star-formation histories at these redshifts. It can be seen from Fig. 16 that our predictions are in agreement with some data at high redshift (e.g Coppin et al., 2015; Watson et al., 2015), and again are consistent with the ALMA HUDF (Dunlop et al., 2017). However, some tension clearly exists with other recent ALMA observations as we explore in more detail below.

6.2.4 IRX-$\beta$ discrepancy with Capak

et al. (2015)

It is clear from Fig. 16 that there is a significant discrepancy between the predictions of our simulations and the observed positions of the $z\approx 5$ galaxies from Capak et al. (2015) in the IRX $-\beta$ plane (see also Bouwens et al., 2016). As discussed previously, under the assumption that the BC(1600)∗ factor is constant (1.75), then the position of galaxies in the IRX $-\beta$ plane is dependent, to first order, only on (i) the slope of the reddening curve and (ii) the intrinsic UV continuum slopes ($\beta_{i}$). Therefore, in an attempt to explain this, it is worth asking, given the above discussion on IRX $-\beta$, what type of reddening law, and what intrinsic properties of galaxies, are required to find relations compatible with the Capak et al. (2015) data?

To do this we derived IRX $-\beta$ curves for two simple scenarios: firstly, the A${}_{1600}-\beta$ relation was calculated for the Gordon et al. (2003) SMC extinction curve at three values of $\beta_{i}$ ($-2.5$, $-2.0$, $-1.5$) and converted into an IRX-$\beta$ relation assuming BC(1600)∗ = 1.75; secondly, the A${}_{1600}-\beta$ relation was calculated for a set of empirical attenuation curves of the form A$(\lambda)\propto\lambda^{n}$ at different values of $n$ ($-0.55$, $-1.24$, $-20.00$) and converted into an IRX-$\beta$ relation assuming a fixed $\beta_{i}=-2.5$ and BC(1600)∗ = 1.75. As highlighted previously, assuming a fixed $\beta_{i}$ is a simplification, but nevertheless it is a reasonable assumption for models including nebular continuum, and adequate for the purposes of this discussion. The derived curves are shown in Fig. 17, along with data points showing the three $z>5$ sources with individual detections, and the stack of undetected sources, from Capak et al. (2015).

Fig. 17 illustrates two important points: firstly, if the reddening law at $z=5$ resembles the SMC extinction law then the implied $\beta_{i}$ of the Capak et al. (2015) galaxies is $\beta_{i}\approx-1.5$. This is possible for individual galaxies, although 2/3 of the Capak et al. (2015) galaxies with ALMA detections have masses in the range probe by our simulations (log(M/M*⊙*) = 9.67, 10.17, 10.39) and we do not find intrinsic slopes redder than $\beta_{i}\simeq-2.2$. However, it is clearly incompatible with the average $\beta_{i}$ of galaxies at this redshift, since the typical observed slopes are bluer than this value $\langle\beta\rangle\approx-2.0$ (Dunlop et al., 2012; Dunlop et al., 2013; Rogers et al., 2014). Secondly, if we assume $\beta_{i}=-2.5$ (i.e. the Capak et al. (2015) galaxies are compatible with our simulations) then the slope of the reddening law required to explain the data is much steeper than the SMC-extinction law. For example, the position of the stack of undetected galaxies would require $n<-20$, much steeper than any attenuation, or even extinction law, ever observed.

Given this discrepancy, it is worth considering the potential systematics that could be affecting the Capak et al. (2015) and Bouwens et al. (2016) results. On possibility is an underestimation of the dust temperature, which is taken to be 35K in both cases. It has been suggested that dust temperature may increase with redshift due to an increase in the intensity of the radiation field heating the dust (e.g. Béthermin et al., 2015)111This increase in radiation field intensity is also consistent with the increase in ionizing photons in star-forming galaxies inferred from the evolution of optical emission line ratios (e.g. Steidel et al., 2014; Kewley et al., 2015; Cullen et al., 2016), and typically an upper limit of 45K has been assumed (Schaerer et al., 2015). In this case LIR will be underestimated by a factor $\approx 2$ and the resulting IRX underestimated by $\approx 0.3$ dex. The magnitude of this shift is illustrated by the arrow in the bottom right-hand corner of Fig. 16. It would bring the Capak et al. (2015) and Bouwens et al. (2016) data into better agreement with our model, however it is not enough to fully account for the discrepancy. There are other potential systematic issues affecting the determination of the UV continuum slopes, and estimated upper limits on LIR, which we will discuss in detail in a future work (McLure et al. in prep).

6.2.5 SMC-like attenuation at lower redshifts?

In Fig. 16 we have not included data at lower redshifts ($z\sim 2-3$) which also support steeper, SMC-like, attenuation laws. In particular, Reddy et al. (2010) find that galaxies with young mass-weighted ages and typically lower masses (age $\lesssim 100$ Myr; M/M${}_{\odot}<10^{10}$) have a smaller IRX, at a given $\beta$, relative to older galaxies at the same redshift. Finally, similar conclusions have also been reached using small numbers of individual lensed galaxies (e.g. Baker et al., 2001; Siana et al., 2009).

In Reddy et al. (2010), these young galaxies represented $\approx 13\%$ of their full star-forming galaxy sample, and in our simulated sample $\approx 6\%$ of galaxies have ages $<100$ Myr. However, since our method involves comparing to the observed luminosity function, it can only be applied to the full galaxy population complete down to a given UV magnitude, and cannot be applied to populations split by age. Therefore, we cannot rule out the possibility that certain populations (e.g. the youngest galaxies), or any given individual galaxy, will follow a steeper attenuation curve. We also note that, since the absolute attenuations decrease towards fainter, lower mass galaxies, the uncertainties on the form of the attenuation law will correspondingly increase, as the difference in absolute attenuation between curves becomes more difficult to distinguish. Therefore, our conclusions are most robust at high masses. Nevertheless, our current analysis suggests that, on average, a Calzetti-like attenuation curve is sufficient to describe the observed properties of the $z\approx 5$ galaxy population down to log(M/M${}_{\odot})\approx 7.5$.

7 Discussion

In this section, we discuss in more detail why a reddening law similar to the SMC extinction law is ruled out by our models, and how this relates to the type of dust present at high redshifts; finally we compare our results to the recent study of Mancini et al. (2016), who performed a similar analysis, but reached different conclusions.

7.1 Why are SMC-like laws ruled out?

It is worth asking, what differences in the intrinsic properties of the simulated SEDs are required such that a steeper attenuation law would be favored over a Calzetti-like attenuation law? The main ways in which our simulations could be biasing our results in favour a shallower attenuation law are: (i) the intrinsic UV continuum slopes are too red; ii) the intrinsic luminosity at $1500\rm{\AA}$ is overestimated; (iii) the number density of galaxies at a given mass is overestimated; or some combination of all three. The guiding principle here is that, for a steep law to be necessary, the change in $\beta$ for a given AUV needs to be larger than is currently required by intrinsic simulated galaxy properties. Therefore, either AUV is kept the same and the intrinsic $\beta$ slopes must move further away from the CMR (i.e. become bluer), or the $\beta$ slopes remain the same and the UV magnitudes of the galaxies must move closer to the LF (i.e. become fainter).

7.1.1 Steeper intrinsic UV continuum slopes?

Firstly, if the the UV continuum slopes are too red then either the UV-weighted ages, or Fe/H values, of the galaxies are being overestimated by the simulations. The range of UV-weighted ages of the FiBY galaxies is $20\lesssim t_{UV}\lesssim 80$ Myr with a median of 33 Myr (see Fig. 6), and the median intrinsic stellar UV continuum slope across all models is $\beta_{i}\approx-2.6$ (see Table. 1). The bluest UV slopes possible with the current BPASSv2 models are $\beta_{i}\approx-3.1$ (a starburst of age 1 Myr and metallicity 1/10 solar). Therefore, as a simple test we fixed the intrinsic slopes of all SEDs at this value, and a $\chi^{2}$ fit with the TEx model was performed as described in Section 5. Indeed, in this case the best-fitting values of $n$ are in the range $-1.0<n<-1.2$ (though the fits are formally significantly poorer), nevertheless it serves to demonstrate that if the intrinsic slopes were close to the steepest allowed by current SPS models, SMC-like reddening laws would become more acceptable.

However, to achieve intrinsic slopes approaching this value ($\beta_{i}\lesssim-2.9$) would require that the UV-weighted ages of our galaxies be in the range $0-10$ Myr (e.g. Stanway et al., 2016). Naturally this would lead to a boost in the SFRs and sSFRs; for example, in the extreme scenario of all stars being formed within the past 10 Myr the sSFR is 100 Gyr*-1* ($20\times$ the median value of our sample). However, the SFRs and sSFR of galaxies in the FiBY simulation are in good agreement with observations. For example, Fig. 18 compares the star-forming main sequence of our simulated sample to the observed $z=5$ relation taken from Schreiber et al. (2015), and clearly there is excellent consistency here. Similarly, the median sSFR of our sample (4.8 Gyr*-1*) is consistent with a number observations at $z\simeq 5$ (e.g. Stark et al., 2013; González et al., 2014; Mármol-Queraltó et al., 2016). Furthermore, as demonstrated in Fig. 7, the effect of the nebular continuum increases as the intrinsic stellar continuum becomes bluer, such that the UV continuum slopes are reddened to roughly the same values $\beta_{i}\approx-2.4$. Finally, for completeness, we briefly note that, if we have underestimated either the ionization parameter or electron density in our photoionization modelling, then we will have underestimated the contribution of the nebular continuum to the galaxy SEDs. Assuming more extreme values of these parameters would make the intrinsic slopes redder, not bluer. In summary, it seems unlikely that the intrinsic $\beta$ slopes at $z=5$ are being overestimated to the extent that a steep SMC-like reddening law is required to match observations.

7.1.2 Overestimated UV luminosities?

As we discussed previously, the intrinsic M1500 is relatively insensitive to the metallicity and nebular continuum contribution, but rather tracks the star-formation rate over $\approx$ 100 Myr timescales. Therefore, M1500 will be overestimated if the SFRs of the galaxies are overestimated. As with the UV continuum slopes, we performed a test whereby the intrinsic M1500 of the galaxies are increased by an amount $\delta\rm{M}_{1500}$, and a $\chi^{2}$ fit with the TEx model carried out; we found that an increase of the order $\Delta\rm{M}_{1500}\approx 1.0$ was required to make the data consistent with a steep SMC-like reddening law. Assuming the Kennicutt & Evans (2012) conversion (in which 1 M*⊙yr-1* is equivalent to an absolute magnitude of M${}_{1500}\simeq-18.75$) this shift in absolute magnitude corresponds to a factor $\simeq 2.5$ decrease in SFR, or in log space $\Delta$log(SFR)$=-0.4$, which would make the data incompatible with the observed main sequence at $z=5$ (Fig. 18). Again, given the consistency of our data with observed SFRs at the same redshift, it is unlikely that we are overestimating M1500. Nevertheless, it is worth bearing in mind that a shift in SFR of this magnitude would be required if the observed stellar masses at these redshifts were overestimated by $\approx 0.4$ dex. However, given the consistency of our data with observed mass $-$ SFR relations at the same redshift, and assuming the observed masses at these redshifts are not systematically overestimated, it is unlikely that we are overestimating M1500.

7.2 Dust properties at $\mn@boldsymbol{z=5}$

What do our results mean for the ‘type’ of dust in high redshift galaxies? Results in the literature, which have pointed towards an attenuation curve at $z\simeq 5$ similar in shape to the extinction curve of the SMC, have been used to claim that the properties of the dust grains themselves at high redshift are more SMC-like (Capak et al., 2015). This seems plausible, since the low metallicity of the SMC ($\simeq 0.25$ Z*⊙* Pagel et al., 1978; Yozin & Bekki, 2014), is more representative of high-redshift galaxies, and hence the physical properties of the dust grains may also be comparable. However, it is not necessarily true that SMC-like dust will result in an attenuation curve similar to the SMC extinction curve, due to the significant difference between extinction and attenuation.

This point is illustrated in Fig. 19, where we show theoretical attenuation curves from the 3D radiative transfer models of Seon & Draine (2016). In these models the output attenuation curve is dependent on the system geometry (for which we have adopted the parameters $M_{s}=20$ and Rs/R${}_{d}=1$, which represent a turbulent clumpy median in which the stars are distributed throughout the total volume), the amount of dust (AV) and the dust properties, for which we have selected, for illustrative purposes, the Milky Way model from Weingartner & Draine (2001), and the SMC model of Witt & Gordon (2000). The value of AV we have selected is typical of the median of our sample across all best-fitting models (A${}_{V}\approx 0.5$). It can be seen that, even assuming SMC dust properties, the output attenuation curve is much shallower than the SMC extinction curve, though still steeper than the Calzetti curve. However, it is fully consistent with our best fitting model within $1\sigma$. Assuming Milky Way dust gives better agreement with the Calzetti curve, though it is slightly shallower. Therefore, assuming this scenario, a Calzetti-like attenuation curve at $z=5$ is consistent with dust properties somewhere intermediate between SMC and Milky Way. Moreover, given the inherent uncertainties in radiative transfer modelling, it is not implausible to think that SMC-like dust can result in a Calzetti-like attenuation curve; indeed, as has been claimed in the past (e.g. Witt & Gordon, 2000).

7.3 Comparison to Mancini

et al. (2016)

Recently, Mancini et al. (2016) (M16) performed an analysis similar to the one presented here, comparing the intrinsic CMR and LF of a sample of simulated galaxies at $5<z<8$ to observations, in order to constrain the attenuation law and subsequent $\rm{IRX}-\beta$ relation. Contrary to our empirical prescription, M16 estimate the absolute attenuation (or the normalization of the reddening curve) using a semi-analytic dust evolution model which attempts to explicitly track dust growth / destruction in stars and supernovae. Combining this with a variety of attenuation laws and a CF2000-like model, they conclude that matching their simulations to observations re quires a steep SMC-like curve.

In the method adopted by M16, the A1500 versus stellar mass relation is fixed by the dust model, and their relation is steeper than the one we derive by fitting to the LF (private communication). This can be seen by inspecting Fig. 2 of M16, where the number counts at faint intrinsic magnitudes ($\rm{M_{UV}}\approx-18.0$) from their adopted simulation are much closer to the observed number counts than the number counts at faint magnitudes from FiBY (see Fig. 9 of this paper). In other words, M16 require much smaller dust corrections at lower masses. Given this, and the fact that the intrinsic UV slopes in M16 appear very similar to the median values we find in our sample ($\beta_{i}\approx−2.4$, Table 1), it is unsurprising that they require a steeper attenuation law.

One probable explanation of this discrepancy is the lower particle mass resolution of the M16 simulation (log(M/M*⊙) = 7.11 compared to 5.68 and 4.81 for the two components of the FiBY simulation) since the stellar mass range of the galaxies in our sample with $\rm{M_{UV}}\approx-18.0$ (7.6 $<$ log(M/M⊙*) $<$ 7.9) would correspond to N $\approx 3-10$ individual star particles in their simulation (compared to N $\approx 600-1000$ in FiBY_L). This low resolution at faint magnitudes may result in an underestimation of the number counts (a well known effect at the resolution limit, e.g. Finlator et al. (2011)), and hence lead to the perception of very little dust in low mass galaxies. Indeed, comparing Fig 2. of M16 to the intrinsic luminosity functions of Dayal et al. (2013), who analyse a higher resolution box from the same simulation, would appear to support this explanation. Without being able to perform a direct comparison to the M16 data it is difficult to identify the source of the discrepancy conclusively. Nevertheless, we believe that the main driver is the difference in number counts at the faint end of the intrinsic LF and that this is most likely a result of the difference in particle mass resolution of the two simulations.

8 Conclusions

We have presented an investigation of the dust attenuation law at $z=5$ using a sample of N=498 galaxies, with M${}_{1500}\leq-18.0$ and $7.5<\rm{log(M/M}_{\odot}\rm{)}<10.2$, taken from the First Billion Years (FiBY) simulation. Synthetic SEDs were generated for these galaxies using the BPASSv2 SPS models, both including and excluding the effects of massive binaries and nebular emission. The intrinsic UV properties of these SEDs were then compared to the observed $z=5$ LF and CMR to constrain the properties of dust attenuation. Our main findings can be summarized as follows:

Due to the fact that the galaxies in our sample exhibit mainly rising star-formation histories, and since the age of the oldest star particles in our $z=5$ simulation ($\simeq$ 0.9 Gyr) are less than the typical timescale for the ISM to be enriched with Fe via SNe Ia ($\simeq$ 1 Gyr), we find that the galaxies in our sample have $\alpha$/Fe ratios enhanced by a factor $\approx 4$ relative to the solar value. This is consistent with the recent observation of galaxies at $z=2.4$ by Steidel et al. (2016), and suggests care must be taken when modelling the SEDs of high-redshift galaxies, or any galaxies with a rising star-formation history, since the shape of the UV spectrum is most sensitive to the Fe/H ratio. 2. 2.

Applying a simple dust model (the TEx model) in which stars below a certain age (tBC) are assumed to be completely obscured by their natal birth clouds, and the remaining stars experience dust attenuation assuming a homogeneous mixture of stars and dust/gas (A$(\lambda)\propto\lambda^{n}$), we find that the observed CMR and LF can only be recreated if t${}_{BC}<3$ Myr and $-0.6\leq n\leq-0.3$. The constraints on tBC are consistent with the age of the youngest OB associations observed in local galaxies (Massey et al., 1995), and the slope of the attenuation curve consistent with a Calzetti-like attenuation curve. 3. 3.

Assuming a more complicated dust model, similar to that of Charlot & Fall (2000), in which stars with ages $<$ tBC experience enhanced optical depths rather than being completely obscured, and again assuming $\rm{A}(\lambda)\propto\lambda^{n}$, we find that an acceptable match to the observed CMR and LF can only be achieved with $-0.7\leq n\leq-0.3$. Although other free parameters within the model (tBC, $\mu$) are not well constrained using only these observations, the formal best-fitting values are consistent with those originally proposed in Charlot & Fall (2000). 4. 4.

Across both dust models, a steep attenuation law resembling the SMC extinction law is ruled out by our simulations. This is in contrast to claims based on some recent ALMA observations at high redshift (e.g. Capak et al., 2015; Bouwens et al., 2016), and the simulations of Mancini et al. (2016). In order for the simulations to become more compatible with a steeper reddening curve, the intrinsic UV magnitudes must be fainter and/or the intrinsic UV continuum slopes must be steeper. However, we argue that, based on the good agreement between the simulation data and observations of the star-forming main sequence of galaxies at $z=5$, it is unlikely that we are overestimating either of these quantities. Nevertheless, we note that if the observed stellar masses of galaxies at $z=5$ are being overestimated by $\approx 0.4$ dex, this would result in a shift towards more SMC-like curves. 5. 5.

Comparing the IR predictions of our best-fitting models to the properties of the one $z=5$ source detected in the recent ALMA HUDF imaging of Dunlop et al. (2017), we find consistency in both the measured physical properties and implied number counts. Nevertheless, better statistics are clearly needed to make a more robust comparison. 6. 6.

The A${}_{1600}-\beta$ and resulting IRX $-\beta$ relations of our best-fitting models are inconsistent with the recent results of Capak et al. (2015) and Bouwens et al. (2016), which require an attenuation curve as steep, or steeper, than the SMC extinction curve. We find that the best-fitting IRX $-\beta$ curve is consistent with what would be expected for a simple homogeneous star/dust mixture geometry combined with the Calzetti et al. (2000) attenuation law assuming all sources have intrinsic UV continuum slopes of $\beta_{i}=-2.5$. There is clearly some tension here which can only be resolved with future, deeper, ALMA data at these redshifts. 7. 7.

Our best estimate of the $A_{1600}-\beta$ relation at $z=5$ is $A_{1600}=2.10^{+1.8}_{-0.3}(\beta+2.52)$; this equation implies that previous estimates of dust-corrected SFRs at these redshifts, based on the Meurer et al. (1999) relation, will have underestimated the SFR by a factor $\approx 2$. We recommend using this relation for any young galaxy population expected to have intrinsically blue UV continuum slopes.

Overall, we find no evidence for any significant variation from the Calzetti-like attenuation curve at $z=5$. However, more observations will be needed, both in the rest-frame UV and IR, to settle the question of the shape of the dust attenuation law at $z\geq 5$.

9 Acknowledgments

FC acknowledges the support of the Science and Technology Facilities Council (STFC). RJM acknowledges the support of the European Research Council via the award of a Consolidator Grant (PI McLure). JSD acknowledges the support of the European Research Council via the award of an Advanced Grant, and the contribution of the EC FP7 SPACE project ASTRODEEP (Ref. No. 312725). This research made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration et al., 2013), NumPy and SciPy (Oliphant, 2007), Matplotlib (Hunter, 2007), IPython (Pérez & Granger, 2007) and NASA’s Astrophysics Data System Bibliographic Services.