The First Billion Years project: constraining the dust attenuation law of star-forming galaxies at z $\simeq$ 5
F. Cullen R. J. McLure S. Khochfar J. S. Dunlop C. Dalla Vecchia




firstpage–lastpage 2016

Secure SWIPT Networks Based on a  Non-linear Energy Harvesting Model 
    Elena Boshkovska,
 Nikola Zlatanov, Linglong Dai,   Derrick Wing Kwan Ng, and Robert Schober E. Boshkovska and R. Schober are with Friedrich-Alexander-University Erlangen-Nürnberg (FAU), Germany.  Linglong Dai is with Tsinghua University, Beijing, China.  L. Dai is supported by the International Science & Technology Cooperation Program of China (Grant No. 2015DFG12760) and the National Natural Science Foundation of China (Grant No. 61571270).
D. W. K. Ng is with The University of New South Wales, Australia. N. Zlatanov is with Monash University, Australia. R. Schober is supported by the AvH Professorship Program of the Alexander von Humboldt Foundation. D. W. K. Ng is supported under Australian Research Council's Scheme Discovery Early Career Researcher Award funding scheme (project
number DE170100137).


    Accepted – . Received 2016 May 30
====================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================







We present the results of a study investigating the dust attenuation law at z≃ 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≃ 5, which have M_1500≤ -18.0 and 7.5 ≤log(M/M_⊙)≤ 10.2, display a mass-dependent α-enhancement, with a median value of [α/Fe]_z=5 ≃ 4 × [α/Fe]_Z_⊙. 
The median Fe/H ratio of the simulated galaxies is 0.14±0.05 which produces steep intrinsic UV continuum slopes; ⟨β_i⟩ = -2.4 ± 0.05. 
Using a set of simple dust attenuation models, in which the wavelength-dependent attenuation is assumed to be of the form A(λ) ∝λ^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 (A_1500=0.5×log(M/M_⊙) - 3.3), and the dust attenuation slope (n) is within the range -0.7 ≤ n ≤-0.3, can successfully reproduce the LF and CMR over a wide range of stellar population synthesis model (SPS) assumptions, including the effects of massive binaries. 
This range of attenuation curves is consistent with a power-law fit to the Calzetti attenuation law in the UV (n=-0.55).
In contrast, 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≃5 IRX-β relation.



galaxies: dust - galaxies: high redshift - galaxies: evolution - 
galaxies: star-forming






§ INTRODUCTION


Tracing the star-formation rate of galaxies across cosmic time, and particularly at z ≳ 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≈8 <cit.>, the lack of far-infrared (FIR) data for galaxies at z ≳ 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 ≳ 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 <cit.> (M99) IRX-β relation (IRX=L_IRL_UV) is used to estimate the intrinsic UV luminosity from the observed UV continuum slope β (f_λ∝λ^β).
This relation was derived from a sample of local starburst galaxies and has been shown to apply for star-forming galaxies out to z ≃ 2 <cit.>.
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 <cit.>, 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 A_UV-β relation measured for local starbursts is what one would expect for a relatively grey reddening law, similar to that derived by <cit.>.
Steeper laws such as the Small Magellanic Cloud (SMC) extinction curve <cit.>, or laws which include a strong 2175 Å absorption feature <cit.> are not compatible with the local starburst galaxies.
Therefore, when applying the <cit.> IRX-β relation at z ≳ 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 <cit.>, 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-β relation, by far the most common assumption is the <cit.> 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 ≃ 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 <cit.>.
Similar results have also previously been reported for young galaxies (< 100 Myr old) <cit.> and sample of gravitationally lensed galaxies <cit.> 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-β relation becomes much shallower).
If dust reddening at high redshift is SMC-like, then previous UV luminosity corrections made assuming the <cit.> IRX-β relation, or <cit.>, will have overestimated the dust attenuation, and hence overestimated the intrinsic luminosities and SFRs at z ≃ 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, <cit.> found a Calzetti-like attenuation curve shape at λ<2600Å for a sample of N=224 star-forming galaxies at 1.4 < z < 2.6.
Using an independent method, <cit.> 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, <cit.> 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 <cit.> 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) <cit.>, 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 <cit.>; 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 <ref> 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 <ref> 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 <ref> 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 <ref>. In Section <ref> we discuss the IR properties of z=5 galaxies predicted by our models and compare these to recent ALMA observations. Finally, in Sections <ref> and <ref> 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: Ω_n =0.265, Ω_b = 0.0448, Ω_Λ =0.735, H_0 =71 km s^-1 Mpc^-1 and σ_8=0.81.



§ 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.



 §.§ SPS models


To construct the SED for each galaxy we used the latest release of `Binary Population and Spectral Synthesis' <cit.>.
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 (≳ 70%) of massive stars undergo a binary interaction (such as mass transfer or a merger) during their lifetimes <cit.>.
Moreover, recent studies of galaxies at z ≃ 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 <cit.>.
We have checked that, for the purposes of this work, BPASSv2 single star models are equivalent to the more commonly adopted STARBURST99 models <cit.>.
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.



 §.§ 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 <cit.>. 
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 <cit.>. 
These simulations reproduce the stellar mass function and star formation rate of galaxies at z≥ 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 <cit.>, and that supernovae-feedback can drive the formation of cores in galaxies <cit.>.

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 <cit.> using tables pre-calculated with CLOUDY v07.02 <cit.>. 
Furthermore, the simulation incorporates full non-equilibrium primordial chemistry networks <cit.> including molecular cooling functions for H_2 and HD. 
Star formation is modelled using the pressure law implementation of <cit.>, which yields results consistent with the Schmidt-Kennicutt law <cit.>. 
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 <cit.>. 
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 (≈ 3.7 × 10^4 Mpc^3 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 <cit.>.
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 <ref>, 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 ≤ 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 10^6 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 M_1500≤ -18.0 where M_1500 was calculated by integrating the rest-frame BPASSv2-100bin SED within a top-hat filter of width 100 Å at a central wavelength of 1500 Å.
This magnitude limit differs slightly according to the assumed SPS model, the limiting value of 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 ≡ SFR / M_*) of our z=5 sample are shown in Fig. <ref>.
The maximum galaxy mass in our sample is 1.4 x 10^10 M_⊙, with a median value of 2.3 x 10^8 M_⊙; the sample is mass complete down to 2.6 x 10^8 M_⊙.
The star-formation rates range from 0.3 - 60 M_⊙yr^-1 with a median of 1.1 M_⊙yr^-1.
The median sSFR is 4.8 Gyr^-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 α/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 <cit.>.
For rising (or constant) star-formation histories the rate of enrichment of α elements from Type II (core-collapse) supernovae, which has a typical time scale of ∼ 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 ∼ 1 Gyr.
Therefore, rapidly star-forming galaxies are expected to have enhanced α/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 α/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. <ref> shows the distributions of UV-weighted Fe metallicity and total metallicity of the star particles (hereafter referred to, for simplicity, as the α element metallicity), relative to solar, for the N=498 galaxies in our sample, where the solar mass fractions have been taken from <cit.>. 
Consistent with the rising star-formation histories, the distributions in Fig. <ref> (left-hand panel) show the difference between the α and Fe-based metallicities.
The median α/Fe ratio in our sample is ≃ 3.9 times larger than the solar value, with the median Fe and α metallicities being Z_Fe/Z_⊙ = 0.14 and Z_α/Z_⊙ = 0.53 respectively.
The right-hand panel of Fig. <ref> shows that, under both definitions, metallicity is an increasing function of stellar mass, however Z_α/Z_⊙ shows a stronger evolution, implying that the α/Fe ratio is also an increasing function of stellar mass.

Non-solar α/Fe ratios have an important effect on the details of the emergent SED.
Recently, <cit.> 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 ≃ 4 - 5 between the stellar and nebular metallicities (Z_*/Z_⊙≃ 0.1 and Z_neb/Z_⊙≃ 0.4 - 0.5).
This is consistent with the offset in our sample (Fig. <ref>); in fact, restricting our sample to galaxies above 10^9 M_⊙ yields a median α/Fe ratio of 4.2.
As argued by <cit.>, this offset is naturally explained if star-formation histories at these redshifts are typically rising <cit.>, 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 α elements.
To a first approximation, any differences between stellar and nebular metallicities is therefore plausibly a result of the same α/Fe offset illustrated in Fig. <ref>; the lower stellar metallicity is driven by the sensitivity of the stellar UV spectrum to Fe/H, whereas the nebular metallicity follows α/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 α/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 <cit.>, 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. <ref>, which shows the average composite BPASSv2-100bin spectra of all N=73 galaxies with log(M/M_⊙) ≥ 9.0, assuming both Z_*=Z_Fe and Z_*=Z_Tot, and will be discussed further in Section <ref>. 




 §.§ 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  regions in the galaxy. 
To model the nebular continuum and emission spectrum we use the photoionization code Cloudy <cit.>.
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  region.
Following the discussion in Section <ref>, we used the UV-weighted alpha-element stellar metallicity to set the nebular abundances, scaling to the <cit.> 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≃2-3 which have illustrated that the typical values of these parameters are evolving from z=0 out to higher redshifts <cit.>.
For example, <cit.> have presented a detailed analysis of N∼ 380 galaxies at z∼2.3 and concluded that the typical value of electron density is ≃ 300 cm^-3 <cit.>, with ionization parameters in the range -3.1 < log(U) < -2.5, consistent with the best-fitting value of log(U) = -2.8 from <cit.> and again representing an increase on typical SDSS values <cit.>.

We made the assumption that the average z∼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 <cit.> (which have medians of 1.0 x 10^10 M_⊙ and 24 M_⊙yr^-1, factors ∼ 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 <cit.> 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 ≃ 5.
Secondly, the metallicities of our sample (both Fe/H and α/H) follow their best-fitting values at z ≃ 2 - 3 in the relevant mass range. 
Given the similarities in sSFRs and metallicities, we made the assumption the other physical conditions in  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 <ref>, how our results are affected under the assumption that physical conditions in  regions continue to evolve towards even more extreme values of these parameters at z ≳ 3.

The effect of adding the nebular continuum to the UV spectrum is illustrated in Fig. <ref> for a typical galaxy in our sample (M_* = 2.8 x 10^8 M_⊙; SFR = 1.3 M_⊙yr^-1; sSFR = 4.6 Gyr^-1; Z_Fe/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 β slopes (1268 < λ < 1950 Å) is ≃ 6 %, and the overall effect of the nebular continuum is to make the spectrum redder (i.e. increase the value of β).
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. <ref>  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 β across the whole sample is given in Section <ref>.

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≃ 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.

    

    



§ INTRINSIC Β AND M_1500 DISTRIBUTIONS


As will be discussed in detail in Sections <ref> and <ref>, 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 (β) and (ii) the volume density of galaxies as a function of their absolute magnitude at λ = 1500 Å (M_1500). 
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.



 §.§ UV continuum slope (β)


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

    f_λ∝λ^β

where β is the power-law index <cit.>.
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.



  §.§.§ Intrinsic stellar UV continuum


To measure β for our synthetic SEDs we performed a linear fit to logf_λ versus logλ over the wavelength range 1268 < λ < 2580, masking out stellar and interstellar absorption features using the windows defined in <cit.>, we refer to this as a `spectroscopic' β measurement.
As we will describe in Section <ref>, we used β slopes of z=5 galaxies reported in <cit.> as our observational comparison.
These authors measured β by fitting a power-law to five photometric bands which sampled the rest-frame UV spectrum in the wavelength range λ≈ 1500 - 2600 Å. 
To ensure consistency, we checked that our spectroscopic β measurements were consistent with the <cit.> photometric method.
We found that both methods are consistent for idealized data (see Fig. <ref>), with the spectroscopic β measurements systematically smaller by a median of Δβ = 0.03 across the whole sample.
For our sample, we chose to use the spectroscopically-measured β values rather than explicitly mimic observational methods.

    

Fig. <ref> shows how the intrinsic stellar UV continuum (β_i) of our sample varies as a function of the UV-weighted age (defined as age weighted by the luminosity at 1500Å of each star particle) and the UV-weighted Fe-based metallicity (Z_Fe).
We note that the data in Fig. <ref> are taken from the SEDs generated with the BPASSv2-100bin models, however, although the normalization of β_i differs slightly, the trends are the same across all SPS models.
As expected, the value of β_i increases with increasing UV-weighted age and UV-weighted metallicity. 
Moreover, we find that, at z=5, β_i is more strongly correlated with the UV-weighted age within the range 10 - 80 Myr probed by our sample, with the scatter in β_i, at fixed UV-weighted age, ≈ 1/3 × the scatter at fixed metallicity.
β_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. <ref>, the median β_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 <ref>). 
Independent of SPS model, the dispersion of β_i values is σ_β_i≈ 0.08. 

Given the correlation between metallicity and β_i shown in Fig. <ref>, 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. <ref>).
We find that the offset between β_i, Fe and β_i, Tot is an increasing function of stellar mass, consistent with the trend shown in the right-hand panel of Fig. <ref>.
For the BPASSv2-100bin model the offset can be as large as δβ_i≈ 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 α/Fe ratios.

    

    

    




  §.§.§ Effect of the nebular continuum


As illustrated in Figs. <ref> and <ref>, 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 β_i across our full sample.
For the remainder of this paper, β refers to the UV continuum slope of the galaxies including the nebular continuum.

Fig. <ref> shows the relationship between β and β_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 β).
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 Δβ = 0.4, but decreases to Δβ≈ 0.05 at β_i≈ -2.35, with a medians of Δβ≈ 0.15 - 0.25 across all models.
In other words, the effect of the nebular continuum on β is of the order ≈ 2 - 15 %, with an average of ≈ 10 %.

It is also evident from Fig. <ref> that, while the difference between binary and single-star models is minor, Δβ is systematically larger for the BPASSv2-300 models. 
At its maximum, the magnitude of this effect represents a difference in Δβ, at fixed β_i, of ≈ 0.05; interestingly, this is large enough to mean that while the median β_i is bluer for the BPASSv2-300 compared to the BPASSv2-100 models, the situation is reversed for the β medians (see Table <ref>).

Finally, across our sample, including the nebular continuum has the effect of reducing the scatter in intrinsic β values by a factor ≈ 2 (see Table <ref>) and significantly flattens the relations between β and metallicity and UV age (Fig. <ref>).
The median β across all SPS models is ≈ -2.40, and we find that the minimum value of β in the case of no dust attenuation at z=5 is β_min = -2.50.
This theoretical dust-free β limit is consistent with the work of <cit.> and <cit.> at z = 6 - 8 and, as we will return to later in the paper, significantly different to the value of β_min = -2.23 assumed for the <cit.> IRX -β relation.



 §.§ UV continuum magnitude


The UV continuum magnitude (M_1500,i) is primarily sensitive to the star-formation history of a galaxy over the last 10 - 100 Myr <cit.>.
As mentioned in Section <ref>, we calculated the intrinsic absolute UV magnitude magnitude (M_1500,i) by integrating the rest-frame SEDs over a top-hat filter of width 100 Å at a central wavelength of 1500 Å.

Since our sample was defined by requiring an intrinsic stellar magnitude (i.e. dust free, no nebular contribution) M_1500,i≤ -18.0, the faintest galaxies across all SPS models have roughly this magnitude (see Section <ref>); the (median, maximum) M_1500,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 ≈ 3.7 × 10^4 Mpc^3 volume, we do not find any galaxies intrinsically brighter than M_1500,i≈ -23.4.

As with the β values, the choice of stellar metallicity also affects the intrinsic stellar M_1500,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 M_1500,i is 0.23 mag.
Finally, the effect of the nebular continuum on M_1500,i is also small; across all models, adding nebular continuum increases the brightness of the SEDs by ≲ 0.1 mag.

    



§ 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 <cit.>.
To do this, we considered two simple models for the dust attenuation: (i) a total birth-cloud extinction model (TEx) and (ii) a <cit.>-like model (CF2000) which we describe in detail below.

    







 §.§ 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 ∼ 10 - 30 Myr <cit.>, values which are in good agreement with current predictions from hydrodynamical simulations <cit.>; while estimates of the timescales for the ejection of OB associations from their birth clouds are of the order 1 - 3 Myr <cit.>, consistent with the youngest ages of galactic OB associations <cit.>.

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 <cit.> 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 <cit.>.
The reddening (in magnitudes) as a function of wavelength in this model can be written as,


    A(λ) =
      ∞     if  t ≤ t_BC,
    ϕλ^n       if  t > t_BC,

where ϕ is a linear function in log(M/M_⊙) given by:


    ϕ = a_1log_10(M/M_⊙) + a_0.


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 <cit.>; 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 Å `bump' as illustrated in Fig. <ref>.
It is worth mentioning that the results of this paper are not dependent of the presence or size of the 2175 Å `bump' feature, since the <cit.> windows are selected to mask out the 2175 Å region, and <cit.> have demonstrated that their observationally derived β values are also insensitive to this feature.



 §.§ CF2000


The CF2000 model follows the prescription of <cit.>, with the addition of a mass-dependent normalization (ϕ) as above.
In this model the wavelength-dependent optical depth experienced by a star particle of age (t) is given by,

    A(t, λ) =
      ϕλ^n       if  t ≤ t_BC,
    μϕλ^n       if  t > t_BC,

where μ defines the fraction of the total dust optical depth contributed by the ambient ISM, n is the exponent of the reddening curve, ϕ 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) ϕ as a linear function of the logarithm of stellar mass as in Eq. <ref>.

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, <cit.> estimate μ≈ 1/3, n ≈ -0.7 and t_BC≈ 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.

    




§ 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 <ref>, we discuss how the predictions of our best-fitting dust models compare to the observed IR properties of high-redshift galaxies.



 §.§ TEx


For the TEx model, the mass normalization ϕ was derived by directly fitting to the z=5 luminosity function of <cit.>.
The method is illustrated in Fig. <ref>: for each birth-cloud age and SPS model, the sample was split into six bins of M_1500, the number density of galaxies in each bin was determined, and the value of A_1500 required to match the luminosity function calculated.
A value of A_1500 = 0 was assigned to bins in which the number density fell below the observed relation. 

We were then able to fit the relationship between A_1500 and the median stellar mass in each bin.
The inset panels of Fig. <ref> show linear fits to A_1500 versus stellar mass (bins with A_1500 = 0 were excluded from the fitting).
Using this relation we can then rewrite Equation <ref> for t > t_BC as,

    A(λ, M) = A_1500(M) (λ/1500)^n,

where M ≡ 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 <cit.> by varying the value of n, applying Equation <ref>, and re-calculating β and M_1500.
For a giving mass normalization, the change in β can be expressed as,

    Δβ = β_obs - β_int = 1.297 (A_1268 - A_2580),

where A_1268 and A_2580 are the attenuations at the wavelengths limits of the <cit.> 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. <ref>).
For each value of n, the sample was binned into six bins of M_1500, and the total χ^2 was then computed,

    χ^2 = ∑_i(model(i) - data(i)/σ(i))^2,

where the summation is over all bins and σ(i)=0.1 is taken as the intrinsic β scatter at the faint end of the CMR estimated by <cit.>, and is representative of the β distributions given in Table <ref>.
The results of this fitting procedure are illustrated in Fig. <ref> and the best-fitting n values and their corresponding χ^2 values are listed in Table <ref>.
Figs. <ref> and <ref> show data in which the nebular continuum correction as described in <ref> has been applied; however, for completeness, in Table <ref> 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σ confidence (χ^2 > χ^2_min + 4).
Combining this with the luminosity functions shown in Fig. <ref>, we can further constrain the birth-cloud age to be ≲ 3 Myr, since at ages ≥ 3 Myr the intrinsic volume density of galaxies at M_1500 > -19.0 falls below the <cit.> 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_neb=-0.46 for the nebular corrected spectra, and n_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≥ 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σ level (Δχ^2 = 1), between the different SPS models.
Taken all together, the results suggest a birth cloud age in the range 0 ≤ t_BC < 3, and a attenuation curve slope in the range -0.6 ≤ n ≤ -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 <ref> 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 ≲ n ≲ -1.0, still below the SMC extinction-law value of n_SMC≈ -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≥ 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 ≲ 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 <cit.>.

    


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 ≲ 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.



 §.§ CF2000


To fit the CF2000 model we applied a more detailed analysis, performing a χ^2 minimization across all parameters as listed in Table <ref>.
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 ≤ t_BC) and ambient ISM (t > t_BC) as given by Equation <ref>, and reddened the spectrum accordingly.
We binned the reddened sample into six bins of M_1500 and performed a χ^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 χ^2 minimization, where the χ^2 value for each set of model parameters was taken as the sum of the CMR and LF values (i.e. χ^2=χ_LF^2+χ_CMR^2).
From the resulting χ^2 we constructed marginalized probability density functions for each parameter.

    

    

Fig. <ref> shows the LF and and CMR for the model parameters at the minimum χ^2 (see Table <ref> 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 <ref> it can be seen that acceptable fits can be achieved for both stellar-only and stellar + nebular models.
The insets in Fig. <ref> 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 χ^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σ uncertainties of all models, the range of acceptable n values is -0.7 ≤ n ≤ -0.3, consistent with both <cit.> and <cit.>, 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 ≤ n ≤ -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≃2 <cit.>.
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 ≲ -0.7 (-0.9) at 1σ (2σ) (see Fig <ref>). 
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 <ref>.



  §.§.§ 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 μ values of either 0.3 or 0.4 depending on the SPS and nebular assumptions, and values within 1σ uncertainties in the range 0.1 < μ < 0.8.
This is consistent with the original <cit.> results of μ≈ 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 (dA_1600/dM) in the range 0.2 < dA_1600/dM < 0.7, with a best-fitting value of dA_1600/dM≃ 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.

    



 §.§ Recommended dust model parameters
















Despite the relatively large range of parameter values over which acceptable can be found, we find that a <cit.> 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 ≃ 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 <ref>.
Given that a Calzetti-like curve is fully consistent within the 1σ 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 λ < 2600 Å (Fig. <ref>), for attenuation curves extending from UV to optical we recommend using the <cit.> relation explicitly or, given its similarity, the <cit.> curve.
In addition to the values listed in Table <ref>, we recommend the following relation between A_1500 and log(M/M_⊙), which provides the normalization of the attenuation curve:

    A_1500 = 0.5 ×log(M/M_⊙) - 3.3.

This relation implies A_1500=0 for log(M/M_⊙)≤ 6.6, and that, for an order of magnitude increase in mass, dust attenuation at 1500Å increases by 0.5 mags.














§ 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 ≃ 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 () 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.:

    L_IR = ∫_912^∞ l_λ, i - l_λ dλ,

where l_λ, i is the intrinsic SED, and l_λ is the reddened SED.
This integral converges at an upper wavelength limit of λ≈ 10,000 Å.
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).



 §.§ ALMA HUDF z=5 source


Recently, <cit.> (D17) have presented a 1.3-mm ALMA mosaic covering the full 4.5 arcmin^2 of the Hubble Ultra-Deep Field (HUDF).
Their final sample consists of 16 sources with point-source flux densities S_1.3 > 120 μJy at a mean redshift ⟨ z ⟩ = 2.15.
Interestingly, one source in their sample (UDF12) has a spectroscopic redshift of z=5.00 <cit.>, 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±0.12, and a sSFR =9.29±4.80, placing it in the high mass and high sSFR tail of our sample (see Fig. <ref>).



  §.§.§ 1.3mm number counts


    

To derive S_1.3, first  was converted into a SFR_IR using the same relation adopted by D17 <cit.>, and S_1.3 was calculated using the conversion given in D17: S_1.3 (in μJy) ≈ 3.3 × SFR_IR (in M_⊙yr^-1).
The number counts as a function of S_1.3 are shown by the histograms in Fig. <ref>.
Depending on model assumptions, we find that between N ≈ 3 - 10 galaxies have flux densities above the threshold S_1.3 > 120 μJy within the ≈ 3.7 × 10^4 Mpc^3 volume of our simulation, or equivalently we find a volume density of N(S_1.3 > 120 μJy) ≈ (1-3) × 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 .

Over the 4.5 arcmin^2 HUDF area, the comoving volume between 4.5<z<5.5 is ≈ 1.3 × 10^4 Mpc^3, therefore, to a rough approximation, the simulation volume is ≈ 3 × the HUDF volume at z=5.
Scaling our number counts by this factor results in a prediction of N ≈ 1 - 3 galaxies with S_1.3 > 120 μ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.

    




  §.§.§ 1.3mm flux density versus UV continuum slope


In Fig. <ref> we directly compare observed quantities tracing the absolute attenuation (S_1.3) and dust law (β) of UDF12.
UDF12 has a 1.3mm flux density of S_1.3=154±40 <cit.> and a UV continuum slope of β=-1.70±0.12, where β 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  (β=-1.71).

Fig. <ref> shows that the S_1.3 and β values extracted from the best-fitting dust models are consistent with the observation.
Although only ≈ 1 % of our sample have comparable flux densities (S_1.3≥ 120 μJy), all fall within < 2σ of UDF12 and the majority fall within < 1σ.
Fig. <ref> also demonstrates that, for the majority of models, the flux densities of the S_1.3≥ 120 μJy sources are biased high with respect to the average S_1.3 of galaxies with similar β slopes.
There is no significant difference between the TEx and CF2000 models, or between the various SPS assumptions; UDF12 is consistent with the 1σ 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 ≤ n ≤ -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≃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.



 §.§ The z=5 IRX-β relation


As discussed in the introduction, the <cit.> (M99) IRX-β relation, and the underlying, implicit, A_UV-β 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 <cit.>.
In this section we explore the IRX-β, and underlying A_UV-β,  relations implied by our best-fitting models.



  §.§.§ IRX -β relation


<cit.> define IRX as the ratio of far infrared flux (F_FIR) to flux at 1600Å (F_1600), where F_1600 is defined as a generalized flux of the form F_λ=λ f_λ, and f_λ is the flux density per wavelength interval.
In their derivation of IRX_1600, 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 IRX_1600 can be written as

    IRX_1600≡L_IR/L_1600 = (10^0.4A_1600-1) ∫_912^∞ l_λ', idλ'/L_1600,i.

Here the first term gives the fraction of L_1600,i absorbed by dust and A_1600 is the attenuation at 1600Å.
For a given dust model, A_1600 was calculated by differencing the magnitudes at 1600Å of the intrinsic SED and the attenuated SED, where the magnitudes were calculated using a top-hat filter centered on 1600Å of width 100Å.
The second term is the ratio of the maximum amount of of luminosity available for dust heating (l_λ', i is the unattenuated luminosity of the SED) divided by the intrinsic 1600Å luminosity.

Assuming the second term is a constant, this equation implies a direct transformation between the attenuation at 1600Å and IRX.
Then, since for a given attenuation curve the attenuation at 1600Å is related to the observed β slope, the position of galaxies in the IRX-β plane can be used to infer the properties of the attenuation curve <cit.>, although one must be wary of the crucial dependence on the intrinsic β slopes as we will discuss below.

In <cit.>, the second term in the equation, denoted BC(1600)_*, was calculated using 100 Myr duration constant star-formation rate SPS models from <cit.>.
They found 1.56 ≤BC(1600)_* ≤ 1.76, adopting a final value BC(1600)_* = 1.66 ± 0.15.
Other authors have adopted the value of 1.75 ± 0.25 derived later by <cit.> <cit.>.
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 ≤BC(1600)_* ≤ 1.78, and for stellar continuum only models we find 1.42 ≤BC(1600)_* ≤ 1.52.
The full range of values is consistent with original derivation of <cit.> but emphasizes how large an effect the nebular continuum can have.
Interestingly, the BC(1600)_* values derived from SEDs including nebular continuum are in much better agreement with previous estimates <cit.>.
For the remainder of this paper we adopt BC(1600)_* = 1.75.
The IRX equation then simply becomes:

    IRX_1600 = 1.75 (10^0.4A_1600-1),

and the relationship between IRX and β becomes dependent only on the A_1600 - β relation.


    




  §.§.§ A_1600-β relation


If one assumes that all sources have the same intrinsic UV slopes (i.e. a constant β_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 A_1600 and β given by combining equations <ref> and <ref>.
For a Calzetti-like reddening law (n=-0.55) this becomes,

    A_1600, Calzetti = 2.09 (β+β_i),

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

    A_1600, SMC = 0.99 (β+β_i).

However, if there is an underlying correlation between β_i and A_1600 (i.e. β_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 β_i and A_1600, such that there is a general β_i - 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 β_i values (see Table <ref>).
The situation is further complicated for `clumpy' geometries like CF2000, since for these models Equation <ref> does not apply, and so simple relations like equations <ref> and <ref> above cannot be derived.
The upshot of all this is that the slope of the A_1600-β 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 β_i in the resulting fits are representative of the minimum values of β_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.



  §.§.§ IRX - β for the best-fitting dust models


    

    

We derived A_1600-β and IRX - β relations for our sample for both the stellar +  nebular and stellar-only cases; the results are shown in Figs. <ref> and <ref>.
For each case we performed a linear least-squares fit to the best-fitting β and A_1600 values from our models (combining the data across both dust models and all SPS models), and substituted the resulting fit into Equation <ref> to give the corresponding IRX -β curve.
The best-fitting A_1600-β relations are as follows (see Fig. <ref>), for the stellar + nebular case:

    A_1600 = 2.10 (β + 2.52),

and for the stellar-only case:

    A_1600 = 1.28 (β + 2.83),

In both cases, the scatter about the relations is σ=0.14.

The relation for the stellar + nebular case is virtually identical to assuming a Calzetti attenuation law for sources with an intrinsic β_i=-2.52 (see Equation <ref>).
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 β_i and best-fitting A_1600.
The effect of the nebular continuum is substantial, up to ΔA_1600≈ 0.5 at the extreme ends of β values.
The scatter of σ=0.14 around the A_1600-β relation derived above only represents the scatter about the best fitting models, to get a more robust estimate from our models we estimated A_1600-β based on the 1σ limits on n (-0.7 < n < -0.3; see Fig. <ref> and Table <ref>).
Taking these limit into account, we estimate the following best-fitting A_1600-β relation for the stellar + nebular case,
 
    A_1600 = 2.10^+1.8_-0.3(β + 2.52).

In reality, a nebular continuum contribution is expected in star-forming galaxies and, though it will likely vary for individual galaxies, we recommend Equation <ref> as the more reliable estimate of the average A_1600-β relationship at z=5.
Futhermore, this equation is reassuringly similar to observational constrainsts based on of samples of Lyman-beak galaxies at z ∼ 3 - 6, especially in terms of the intrinsic UV slopes at these redshits <cit.>.
One immediate consequence of this is that dust corrections to the UV, based on the <cit.> IRX -β relation, will underestimate A_1600 for an observed β, by a median of ≈ 0.6 across the range of observed β values -2.5 < β < -1.6, equivalent to an underestimation of the resulting star-formation rate of a factor ≈ 2.

Finally, in Fig. <ref> we show the predicted position of z ≃ 5 galaxies in the IRX-β plane along with a variety of data from the literature at z>3 <cit.>.
It can be seen that the best-fitting line (using Eq. <ref>) is roughly consistent with the shape of the original <cit.> curve with a shift to bluer intrinsic UV slopes of Δβ_i≃ - 0.2.
Again, this is fully consistent with typical z≃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. <ref> that our predictions are in agreement with some data at high redshift <cit.>, and again are consistent with the ALMA HUDF <cit.>.
However, some tension clearly exists with other recent ALMA observations as we explore in more detail below.



  §.§.§ IRX-β discrepancy with <cit.>


It is clear from Fig. <ref> that there is a significant discrepancy between the predictions of our simulations and the observed positions of the z≈5 galaxies from <cit.> in the IRX -β plane <cit.>.
As discussed previously, under the assumption that the BC(1600)_* factor is constant (1.75), then the position of galaxies in the IRX -β plane is dependent, to first order, only on (i) the slope of the reddening curve and (ii) the intrinsic UV continuum slopes (β_i). 
Therefore, in an attempt to explain this, it is worth asking, given the above discussion on IRX -β, what type of reddening law, and what intrinsic properties of galaxies, are required to find relations compatible with the <cit.> data?

To do this we derived IRX -β curves for two simple scenarios: firstly, the A_1600-β relation was calculated for the <cit.> SMC extinction curve at three values of β_i (-2.5, -2.0, -1.5) and converted into an IRX-β relation assuming BC(1600)_* = 1.75; secondly, the A_1600 - β relation was calculated for a set of empirical attenuation curves of the form A(λ) ∝λ^n at different values of n (-0.55, -1.24, -20.00) and converted into an IRX-β relation assuming a fixed β_i=-2.5 and BC(1600)_* = 1.75.
As highlighted previously, assuming a fixed β_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. <ref>, along with data points showing the three z>5 sources with individual detections, and the stack of undetected sources, from <cit.>.

Fig. <ref> illustrates two important points: firstly, if the reddening law at z=5 resembles the SMC extinction law then the implied β_i of the <cit.> galaxies is β_i≈ -1.5.
This is possible for individual galaxies, although 2/3 of the <cit.> 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 β_i≃ -2.2.
However, it is clearly incompatible with the average β_i of galaxies at this redshift, since the typical observed slopes are bluer than this value ⟨β⟩≈ -2.0 <cit.>.
Secondly, if we assume β_i=-2.5 (i.e. the <cit.> 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 <cit.> and <cit.> 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 <cit.>[This 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 <cit.>], and typically an upper limit of 45K has been assumed <cit.>.
In this case  will be underestimated by a factor ≈ 2 and the resulting IRX underestimated by ≈ 0.3 dex.
The magnitude of this shift is illustrated by the arrow in the bottom right-hand corner of Fig. <ref>.
It would bring the <cit.> and <cit.> 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 , which we will discuss in detail in a future work (McLure et al. in prep).



  §.§.§ SMC-like attenuation at lower redshifts?


In Fig. <ref> we have not included data at lower redshifts (z ∼ 2 - 3) which also support steeper, SMC-like, attenuation laws.
In particular, <cit.> find that galaxies with young mass-weighted ages and typically lower masses (age ≲ 100 Myr; M/M_⊙ < 10^10) have a smaller IRX, at a given β, relative to older galaxies at the same redshift.
Finally, similar conclusions have also been reached using small numbers of individual lensed galaxies <cit.>.

In <cit.>, these young galaxies represented ≈ 13 % of their full star-forming galaxy sample, and in our simulated sample ≈ 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 ≈ 5 galaxy population down to log(M/M_⊙) ≈ 7.5.




§ 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 <cit.>, who performed a similar analysis, but reached different conclusions.

    

 


 §.§ 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Å 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 β for a given A_UV needs to be larger than is currently required by intrinsic simulated galaxy properties.
Therefore, either A_UV is kept the same and the intrinsic β slopes must move further away from the CMR (i.e. become bluer), or the β slopes remain the same and the UV magnitudes of the galaxies must move closer to the LF (i.e. become fainter).



  §.§.§ 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 ≲ t_UV≲ 80 Myr with a median of 33 Myr (see Fig. <ref>), and the median intrinsic stellar UV continuum slope across all models is β_i≈ -2.6 (see Table. <ref>).
The bluest UV slopes possible with the current BPASSv2 models are β_i≈ -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 χ^2 fit with the TEx model was performed as described in Section <ref>.
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 (β_i≲ -2.9) would require that the UV-weighted ages of our galaxies be in the range 0 - 10 Myr <cit.>.
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 × 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. <ref> compares the star-forming main sequence of our simulated sample to the observed z=5 relation taken from <cit.>, 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≃5 <cit.>.
Furthermore, as demonstrated in Fig. <ref>, 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 β_i≈ -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 β slopes at z=5 are being overestimated to the extent that a steep SMC-like reddening law is required to match observations.


    



  §.§.§ Overestimated UV luminosities?


As we discussed previously, the intrinsic M_1500 is relatively insensitive to the metallicity and nebular continuum contribution, but rather tracks the star-formation rate over ≈ 100 Myr timescales.
Therefore, M_1500 will be overestimated if the SFRs of the galaxies are overestimated.
As with the UV continuum slopes, we performed a test whereby the intrinsic M_1500 of the galaxies are increased by an amount δM_1500, and a χ^2 fit with the TEx model carried out; we found that an increase of the order ΔM_1500≈ 1.0 was required to make the data consistent with a steep SMC-like reddening law.
Assuming the <cit.> conversion (in which 1 M_⊙yr^-1 is equivalent to an absolute magnitude of M_1500≃ -18.75) this shift in absolute magnitude corresponds to a factor ≃ 2.5 decrease in SFR, or in log space Δlog(SFR)= -0.4, which would make the data incompatible with the observed main sequence at z=5 (Fig. <ref>).
Again, given the consistency of our data with observed SFRs at the same redshift, it is unlikely that we are overestimating M_1500.
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 ≈ 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 M_1500.




 §.§ Dust properties at 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 ≃ 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 <cit.>. 
This seems plausible, since the low metallicity of the SMC <cit.>, 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. <ref>, where we show theoretical attenuation curves from the 3D radiative transfer models of <cit.>.
In these models the output attenuation curve is dependent on the system geometry (for which we have adopted the parameters  M_s=20 and R_s/R_d=1, which represent a turbulent clumpy median in which the stars are distributed throughout the total volume), the amount of dust (A_V) and the dust properties, for which we have selected, for illustrative purposes, the Milky Way model from <cit.>, and the SMC model of <cit.>.
The value of A_V we have selected is typical of the median of our sample across all best-fitting models (A_V≈ 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σ.
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 <cit.>.



 §.§ Comparison to <cit.>


Recently, <cit.> (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 IRX - β 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 A_1500 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 (M_UV≈ -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. <ref> 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 (β_i≈ −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 M_UV≈ -18.0 (7.6 < log(M/M_⊙) < 7.9) would correspond to N ≈ 3 - 10 individual star particles in their simulation (compared to N ≈ 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. <cit.>), 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 <cit.>, 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.



§ CONCLUSIONS


We have presented an investigation of the dust attenuation law at z=5 using a sample of N=498 galaxies, with M_1500≤ -18.0 and 7.5<log(M/M_⊙)<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 (≃ 0.9 Gyr) are less than the typical timescale for the ISM to be enriched with Fe via SNe Ia (≃ 1 Gyr), we find that the galaxies in our sample have α/Fe ratios enhanced by a factor ≈ 4 relative to the solar value.
This is consistent with the recent observation of galaxies at z=2.4 by <cit.>, 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.


  * Applying a simple dust model (the TEx model) in which stars below a certain age (t_BC) 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(λ)∝λ^n), we find that the observed CMR and LF can only be recreated if t_BC < 3 Myr and -0.6 ≤ n ≤ -0.3.
The constraints on t_BC are consistent with the age of the youngest OB associations observed in local galaxies <cit.>, and the slope of the attenuation curve consistent with a Calzetti-like attenuation curve.


  * Assuming a more complicated dust model, similar to that of <cit.>, in which stars with ages < t_BC experience enhanced optical depths rather than being completely obscured, and again assuming A(λ)∝λ^n, we find that an acceptable match to the observed CMR and LF can only be achieved with -0.7 ≤ n ≤ -0.3.
Although other free parameters within the model (t_BC, μ) are not well constrained using only these observations, the formal best-fitting values are consistent with those originally proposed in <cit.>.


  * 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 <cit.>, and the simulations of <cit.>.
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 ≈ 0.4 dex, this would result in a shift towards more SMC-like curves.


  * 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 <cit.>, 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.


  * The A_1600-β and resulting IRX -β relations of our best-fitting models are inconsistent with the recent results of <cit.> and <cit.>, which require an attenuation curve as steep, or steeper, than the SMC extinction curve.
We find that the best-fitting IRX -β curve is consistent with what would be expected for a simple homogeneous star/dust mixture geometry combined with the <cit.> attenuation law assuming all sources have intrinsic UV continuum slopes of β_i=-2.5.
There is clearly some tension here which can only be resolved with future, deeper, ALMA data at these redshifts.


  * Our best estimate of the A_1600-β relation at z=5 is A_1600 = 2.10^+1.8_-0.3(β + 2.52); this equation implies that previous estimates of dust-corrected SFRs at these redshifts, based on the <cit.> relation, will have underestimated the SFR by a factor ≈ 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 ≥ 5.

 





§ 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 <cit.>, NumPy and SciPy <cit.>, Matplotlib <cit.>, IPython <cit.> and NASA's Astrophysics Data System Bibliographic Services.




mnras                      
       







