Dark-ages Reionization and Galaxy Formation Simulation -- XIX: Predictions of infrared excess and cosmic star formation rate density from UV observations

TL;DR

This study uses a semi-analytic galaxy formation model combined with UV observations to predict infrared excess and cosmic star formation rate density at high redshift, highlighting the impact of dust modeling and star formation activity.

Contribution

It introduces a Bayesian calibration of dust and star formation parameters in a semi-analytic model, providing self-consistent IRX-$eta$ predictions aligned with observations.

Findings

Large intrinsic scatter in IRX-$eta$ driven by specific star formation rate.

Systematic uncertainty of a factor of two in dust-corrected star formation rates.

Best-fit models agree well with UV luminosity functions and color-magnitude relations.

Abstract

We present a new analysis of high-redshift UV observations using a semi-analytic galaxy formation model, and provide self-consistent predictions of the infrared excess (IRX) -- $\beta$ relations and cosmic star formation rate density. We combine the Charlot & Fall dust attenuation model with the Meraxes semi-analytic model, and explore three different parametrisations for the dust optical depths, linked to star formation rate, dust-to-gas ratio and gas column density respectively. A Bayesian approach is employed to statistically calibrate model free parameters including star formation efficiency, mass loading factor, dust optical depths and reddening slope directly against UV luminosity functions and colour-magnitude relations at z ~ 4-7. The best-fit models show excellent agreement with the observations. We calculate IRX using energy balance arguments, and find that the large intrinsic scatter in the IRX -$\beta$ plane is driven by the specific star formation rate. Additionally, the difference among the three dust models suggests a factor of two systematic uncertainty in the dust-corrected star formation rate when using the Meurer IRX - $\beta$ relation at z > 4.

Tables

Table 1. Table 1: The five windows selected from Calzetti et al. ( 1994 ) to fit UV slopes for the on-the-fly calibrations.

	Wavelength range [Å]
1	1342 - 1371
2	1562 - 1583
3	1866 - 1890
4	1930 - 1950
5	2400 - 2580

Table 2. Table 2: Summary of free galaxy and dust parameters.

Parameter	Section		Equation		Description
${\alpha }_{\text{SF}}$	Sec.2.1.1		Eq.3		Star formation efficiency
${\Sigma }_{\text{SF}}$	Sec.2.1.1		Eq.1		Critical mass normalisation
${\eta }_0$	Sec.2.2.1		Eq.9		Mass loading normalisation
${ϵ}_0$	Sec.2.2.1		Eq.10		supernova energy coupling normalisation
${\tau }_{\text{SFR}}^{\text{ISM}}$/${\tau }_{\text{DTG}}^{\text{ISM}}$/${\tau }_{\text{GCD}}^{\text{ISM}}$	Sec.3.1.1/Sec.3.1.2/Sec.3.1.3		Eq.16/Eq.19/Eq.22		Dust optical depth normalisation of ISM
${\tau }_{\text{SFR}}^{\text{BC}}$/${\tau }_{\text{DTG}}^{\text{BC}}$/${\tau }_{\text{GCD}}^{\text{BC}}$	Sec.3.1.1/Sec.3.1.2/Sec.3.1.3		Eq.17/Eq.20/Eq.23		Dust optical depth normalisation of BC
${\gamma }_{\text{SFR}}$/${\gamma }_{\text{DTG}}$/${\gamma }_{\text{GCD}}$	Sec.3.1.1/Sec.3.1.2/Sec.3.1.3		Eq.15/Eq.18/Eq.21		Dust optical depth slope of galaxy property
$n$	Sec.3.1.1/Sec.3.1.2/Sec.3.1.3		Eq.15/Eq.18/Eq.21		Reddening slope
$a$	Sec.3.1.1/Sec.3.1.2/Sec.3.1.3		Eq.15/Eq.18/Eq.21		Dust optical depth redshift dependence
Parameter	Prior scale	Prior range			Best-fit ${}^{\text{a}}$			16/84-th percentiles ${}^{\text{b}}$
		M-SFR	M-DTG	M-GCD	M-SFR	M-DTG	M-GCD	M-SFR	M-DTG	M-GCD
${\alpha }_{\text{SF}}$	log	[0.005, 0.2]	[0.05, 0.18]	[0.04, 0.08]	0.10	0.10	0.05	[0.08, 0.13]	[0.10, 0.11]	[0.05, 0.07]
${\Sigma }_{\text{SF}}$	log	[0.1, 0.8]	[0.001, 0.25]	[0.05, 0.25]	0.19	0.01	0.16	[0.21, 0.42]	[0.007, 0.06]	[0.14, 0.19]
${\eta }_0$	log	[2.0, 12.0]	[2.0, 15.0]	[3.5, 7.5]	4.6	7.0	6.4	[4.0, 7.8]	[6.6, 7.9]	[4.9, 6.1]
${ϵ}_0$	log	[0.35, 0.65]	[0.8, 2.2]	[1.0, 1.7]	0.5	1.5	1.3	[0.4, 0.6]	[1.5, 1.7]	[1.3, 1.5]
${\tau }_{\text{SFR}}^{\text{ISM}}$/${\tau }_{\text{DTG}}^{\text{ISM}}$/${\tau }_{\text{GCD}}^{\text{ISM}}$	linear	[0.5, 2.4]	[0.0, 50.0]	[2.0, 8.0]	1.7	13.5	3.7	[1.4, 1.7]	[9.9, 17.0]	[3.5, 5.3]
${\tau }_{\text{SFR}}^{\text{BC}}$/${\tau }_{\text{DTG}}^{\text{BC}}$/${\tau }_{\text{GCD}}^{\text{BC}}$	linear	[2.0, 10.0]	[0.0, 1000.0]	[25.0, 140.0]	2.5	381.3	69.7	[3.9, 6.6]	[225.1, 476.1]	[60.4, 91.0]
${\gamma }_{\text{SFR}}$/${\gamma }_{\text{DTG}}$/${\gamma }_{\text{GCD}}$	linear	[0.0, 0.6]	[0.4, 2.2]	[1.3, 1.7]	0.19	1.20	1.48	[0.23, 0.32]	[1.05, 1.38]	[1.44, 1.52]
$n$	linear	[-1.00, -0.25]	[-2.5, -0.8]	[-1.6, -1.0]	-0.3	-1.6	-1.3	[-0.5, -0.3]	[-1.7, -1.5]	[-1.4, -1.2]
$a$	linear	[0.00, 0.15]	[0.10, 0.65]	[0.20, 0.55]	0.04	0.34	0.39	[0.02, 0.07]	[0.25, 0.37]	[0.36, 0.42]

Table 3. Table 3: Tabular data of predicted cosmic star formation rate density (SFRD) for the three different dust models. These values are plotted in Figure 7 and are in a unit of log 10 ⁡ ( ρ SFR / ( M ⊙ ​ yr − 1 ​ Mpc − 3 ) ) subscript 10 subscript 𝜌 SFR subscript M direct-product superscript yr 1 superscript Mpc 3 \log_{10}(\rho_{\text{SFR}}/(\,\text{M}_{\odot}\text{yr}^{-1}\text{Mpc}^{-3})) .

z	All galaxies
	M-SFR	M-DTG	M-GCD	M-SFR	M-DTG	M-GCD
4	-1.04	-1.24	-1.28	-1.07	-1.29	-1.32
5	-1.19	-1.44	-1.44	-1.25	-1.51	-1.52
6	-1.38	-1.62	-1.64	-1.46	-1.72	-1.72
7	-1.56	-1.78	-1.83	-1.73	-1.93	-1.97

Equations

$$m_{\text{crit}}={\Sigma }_{\text{SF}}\left(\frac{V_{\text{max}}}{100\text{km/s}}\right)\left(\frac{r_{\text{disk}}}{10\text{kpc}}\right)\times 10^{10}{\text{M}}_{\odot },$$

$$r_{\text{disk}}=3R_{\text{vir}}\frac{\lambda }{\sqrt{2}},$$

$$\Delta m_{\text{star}}={\alpha }_{\text{SF}}\frac{m_{\text{gas}}-m_{\text{crit}}}{t_{\text{dyn,disk}}}\Delta t,$$

$$\Delta m_{\text{reheat}}=\{\begin{array}{ll}\eta \Delta m_{\text{new}},& \Delta E_{\text{SN}}\ge \Delta E_{\text{hot}}\\ \frac{\Delta E_{\text{SN}}}{\text{nicefrac}12V_{\text{vir}}^2},& \Delta E_{\text{SN}}<\Delta E_{\text{hot}}\end{array},$$

$$\Delta E_{\text{hot}}=\frac{1}{2}\eta \Delta m_{\text{new}}{V}_{\text{vir}}^2,$$

$$\Delta m_{\text{eject}}=\frac{\Delta E_{\text{SN}}-\Delta E_{\text{hot}}}{\text{nicefrac}12V_{\text{vir}}^2}.$$

$$\Delta E_{\text{SN}}=ϵ\times {\int }_t^{t+\Delta t}d{t}^{\prime }{\int }_0^{\infty }d\tau \frac{d\epsilon }{d\tau }{\psi }_{}(t^{\prime }-\tau ),$$

$$\Delta m_{\text{new}}=\frac{{\int }_t^{t+\Delta t}d{t}^{\prime }{\int }_0^{\infty }d\tau \frac{d\epsilon }{d\tau }\psi (t^{\prime }-\tau )}{{\int }_0^{\infty }d\tau \frac{d\epsilon }{d\tau }}$$

$$\eta =\{\begin{array}{ll}{\eta }_0{\left(\frac{1+z}{4}\right)}^{{\alpha }_{\text{reheat}}}{\left(\frac{V_{\text{max}}}{60\text{km/s}}\right)}^{-1},& V_{\text{max}}\ge 60\text{km/s}\\ {\eta }_0{\left(\frac{1+z}{4}\right)}^{{\alpha }_{\text{reheat}}}{\left(\frac{V_{\text{max}}}{60\text{km/s}}\right)}^{-3.2},& V_{\text{max}}<60\text{km/s}\end{array},$$

$$ϵ=\{\begin{array}{ll}{ϵ}_0{\left(\frac{1+z}{4}\right)}^{{\alpha }_{\text{eject}}}{\left(\frac{V_{\text{max}}}{60\text{km/s}}\right)}^{-1},& V_{\text{max}}\ge 60\text{km/s}\\ {ϵ}_0{\left(\frac{1+z}{4}\right)}^{{\alpha }_{\text{eject}}}{\left(\frac{V_{\text{max}}}{60\text{km/s}}\right)}^{-3.2},& V_{\text{max}}<60\text{km/s}\end{array},$$

$$\Delta m_{\text{recycle}}={\int }_t^{t+\Delta t}d{t}^{\prime }{\int }_0^{\infty }d\tau \frac{dy}{d\tau }\psi (t^{\prime }-\tau ),$$

$$\Delta m_{\text{reinc}}$$

$$t_{\text{reinc}}$$

$$T_{\lambda }(t)=\{\begin{array}{ll}\exp (-{\tau }_{\lambda }^{\text{ISM}})& t\ge t_{\text{BC}}\\ \exp (-{\tau }_{\lambda }^{\text{ISM}}-{\tau }_{\lambda }^{\text{BC}})& t<t_{\text{BC}}\end{array}.$$

$${\Gamma }_{\lambda }=e^{-az}{\left(\frac{\mathrm{SFR}}{100{\text{M}}_{\odot }/\text{yr}}\right)}^{{\gamma }_{\text{SFR}}}{\left(\frac{\lambda }{1600\text{A˚}}\right)}^n,$$

$${\tau }_{\lambda }^{\text{ISM}}={\tau }_{\mathrm{SFR}}^{\text{ISM}}{\Gamma }_{\lambda },$$

$${\tau }_{\lambda }^{\text{BC}}={\tau }_{\mathrm{SFR}}^{\text{BC}}{\Gamma }_{\lambda },$$

$${\Gamma }_{\lambda }=e^{-az}{\left(\frac{Z_{\text{cold}}}{Z_{\odot }}\right)}^{{\gamma }_{\text{DTG}}}\left(\frac{m_{\text{cold}}}{10^{10}{h}^{-1}{\text{M}}_{\odot }}\right){\left(\frac{r_{\text{disk}}}{h^{-1}\mathrm{kpc}}\right)}^{-2}{\left(\frac{\lambda }{1600\text{A˚}}\right)}^n,$$

$${\tau }_{\lambda }^{\text{ISM}}={\tau }_{\mathrm{DTG}}^{\text{ISM}}{\Gamma }_{\lambda },$$

$${\tau }_{\lambda }^{\text{BC}}={\tau }_{\mathrm{DTG}}^{\text{BC}}{\Gamma }_{\lambda },$$

$${\Gamma }_{\lambda }=e^{-az}{\left(\frac{m_{\text{cold}}}{10^{10}{h}^{-1}{\text{M}}_{\odot }}\right)}^{{\gamma }_{\text{GCD}}}{\left(\frac{r_{\text{disk}}}{h^{-1}\mathrm{kpc}}\right)}^{-2}{\left(\frac{\lambda }{1600\text{A˚}}\right)}^n,$$

$${\tau }_{\lambda }^{\text{ISM}}={\tau }_{\mathrm{GCD}}^{\text{ISM}}{\Gamma }_{\lambda },$$

$${\tau }_{\lambda }^{\text{BC}}={\tau }_{\mathrm{GCD}}^{\text{BC}}{\Gamma }_{\lambda }.$$

$$L_{\lambda }(t)={\int }_0^{t}d\tau {\int }_{{\mathrm{Z}}_{\min }}^{Z_{\max }}dZ\psi (t-\tau ,Z)S_{\lambda }(\tau ,Z)T_{\lambda }(\tau ),$$

$$\begin{array}{rl}\ln \mathcal{L}=& -\frac{1}{2}\sum _i\left[\frac{(n_i^{\mathrm{obs}}-n_i^{\mathrm{model}}{)}^2}{{\sigma }_{\mathrm{LF},\mathrm{i}}^2}+\ln (2\pi {\sigma }_{\mathrm{LF},\mathrm{i}}^2)\right]\\ & -\frac{1}{2}\sum _i\left[\frac{({\beta }_i^{\mathrm{obs}}-{\beta }_i^{\mathrm{model}}{)}^2}{{\sigma }_{\mathrm{CMR},\mathrm{i}}^2}+\ln (2\pi {\sigma }_{\mathrm{CMR},\mathrm{i}}^2)\right].\end{array}$$

$$F_{\text{IR}}={\int }_{912\text{A˚}}^{\infty }\left(L_{\lambda }-L_{\lambda }^{\text{intrinsic}}\right)d\lambda$$

$$\displaystyle F_{\text{UV}}=\lambda L_{\lambda}\big{\rvert}_{\lambda=1600\text{\AA}}$$

$$\text{IRX}=\frac{F_{\text{IR}}}{F_{\text{UV}}}$$

Full text

††pagerange: Dark-ages Reionization and Galaxy Formation Simulation - XIX: Predictions of infrared excess and cosmic star formation rate density from UV observations–12††pubyear: 2017

1 Introduction

One fundamental question in astronomy is to understand the buildup of stars and galaxies from baryonic matter in the early Universe. During this epoch, observations focus mainly on rest-frame UV properties due to cosmic redshift. These include measurements of UV luminosity functions (LFs) (van der Burg et al., 2010; Bouwens et al., 2015; Livermore et al., 2017; Bhatawdekar et al., 2018; Ono et al., 2018), and UV continuum slope to UV magnitude relations (Finkelstein et al., 2012; Bouwens et al., 2014; Rogers et al., 2014), which are also known as the colour-magnitude relations (CMRs). The UV luminosity is a tracer of star formation since most UV photons are emitted by young stars. However, star formation can be heavily obscured by the interstellar dust. One commonly adopted approach to perform dust corrections at high redshifts is to infer the infrared excess (IRX) from the observed UV slopes using a relation calibrated by Meurer et al. (1999) (e.g. Bouwens et al., 2015; Mason et al., 2015; Liu et al., 2016). However, the Meurer et al. (1999) relation is calibrated against local starburst galaxies, and observations of far infrared data is rather challenging at high redshifts. Recent observations at $z\gtrsim 3$ show large scatter in the IRX - $\beta$ relation (Capak et al., 2015; Álvarez-Márquez et al., 2016; Bouwens et al., 2016; Barisic et al., 2017; Fudamoto et al., 2017; Koprowski et al., 2018). For instance, the observed IRX by Bouwens et al. (2016) is much lower than the Meurer et al. (1999) relation, while Koprowski et al. (2018) suggest that the IRX - $\beta$ relation does not evolve with redshift. These observations motivate investigation of the IRX - $\beta$ at high redshifts from theoretical models.

Theoretical studies of dust extinction require intrinsic galaxy properties as input, and one approach is to postprocess the output of a hydrodynamical simulation. This method has been implemented in Safarzadeh et al. (2017) and Narayanan et al. (2018) to investigate the origin of the IRX - $\beta$ relation. At $z\gtrsim 5$, the IRX - $\beta$ relation has been studied by Mancini et al. (2016), Cullen et al. (2017) and Ma et al. (2019). However, their results suggest different extinction curves. Cullen et al. (2017) pointed out that the reason for the disagreement could be due to systematics associated with different simulations.

Semi-analytic models are another popular approach for studying galaxy formation (e.g. Guo et al., 2011; Somerville et al., 2015; Croton et al., 2016; Lacey et al., 2016; Cora et al., 2018; Lagos et al., 2018; Cousin et al., 2019b), and their results can also be used as input to dust models. Semi-analytic models solve a system of differential equations that govern the mass accretion and transition of several key baryonic components of galaxies such as gas and stellar mass. The construction of these models is relatively simple, and hence they are computationally efficient. These models also introduce several free parameters to describe the unknown efficiency or strength of certain physics processes. These parameters bring flexibility, and allow the exploration of different galaxy formation scenarios, which is very useful for identifying which galaxy processes regulate certain observations.

The Dark-ages Reionization And Galaxy Observables from Numerical Simulations (DRAGONS) project111http://dragons.ph.unimelb.edu.au/ introduces the meraxes semi-analytic model (Mutch et al., 2016), which is coupled with the high cadence Tiamat N-body simulation (Poole et al., 2016). The model concentrates on studying galaxy formation at high redshifts. This work utilises meraxes to predict intrinsic galaxy properties, and combines it with a simple and flexible dust attenuation model. The dust optical depths are calculated empirically using relevant galaxy properties. By taking full advantage of the fast computational speed of both the galaxy formation and dust models, we carry out a Bayesian analysis on all the model free parameters, and use UV LFs and CMRs as constraints, which are the most fundamental observables at high redshift. This approach allows these observations to put direct constraints on both galaxy formation and dust parameters, and provides self-consistent predictions of the IRX and star formation rate (SFR).

We organise the paper as follows. Section 2 provides an overview of our meraxes galaxy formation model, and introduces several updates on the model for this work. Section 3 describes the dust models that are integrated into meraxes and the computation of galaxy spectral energy distributions (SEDs). The description of our calibration method can found in Section 4, and the results are discussed in Section 5. We demonstrate the predicted IRX - $\beta$ relations and cosmic star formation rate density (SFRD) in Section 6 and Section 7 respectively. Finally, this work is summarised in Section 8. Throughout the paper, we adopt a flat $\Lambda$CDM cosmology, with $(h,\Omega_{\text{m}},\Omega_{\text{b}},\Omega_{\Lambda},\sigma_{8},n_{\text{s}})=(0.678,0.308,0.0484,0.692,0.815,0.968)$ (Planck Collaboration et al., 2016). Magnitudes are in the AB system (Oke & Gunn, 1983).

2 Galaxy formation model

2.1 Overview

The meraxes semi-analytic model1 is the backbone of the present work. It extends the models of Croton et al. (2006) and Guo et al. (2011) to high redshifts, and is modified to run on high cadence halo merge trees with a delayed supernova feedback scheme. It also implements gas infall, radiative cooling, star formation, supernova feedback, metal enrichment, and reionisation feedback. A detailed description of the model can be found in Mutch et al. (2016), hereafter M16. The active galactic nuclei (AGN) feedback of the model is later introduced by Qin et al. (2017). This work also applies several updates to the model, aiming to improve the predicted gas phase metallicity, which is an input of galaxy SEDs. These will be introduced in Section 2.2.

We utilise the halo merger trees of the Tiamat N-body simulation (Poole et al., 2016, 2017) as input to our semi-analytic model. The simulation contains $2160^{3}$ particles in a $(67.8h^{-1})^{3}\,\text{Mpc}^{3}$ box, with mass resolution $m_{\rm p}=2.64\times 10^{6}h^{-1}\,\,\text{M}_{\odot}$. Halos and friends-of-friends groups are identified using subfind (Springel et al., 2001). The timestep of the simulation is $11.1\,\text{Myr}$ between $z=35$ and $z=5$ and is evenly distributed in dynamical time between $z=5$ and $z=1.8$. The high cadence of the simulation is critical to this study since UV magnitudes are sensitive to starbursts in the recent 100 Myr.

Since this work requires evaluating the model many times, and does not focus on ionising structures, we adopt homogeneous reionisation feedback (Gnedin, 2000) instead of using 21cmfast (Mesinger & Furlanetto, 2007). Both approaches are described in M16 and found to have almost the same predictions on global galaxy properties such as the stellar mass function up to $z\simeq 5$. However, the homogeneous prescription is more computationally efficient.

2.1.1 Star Formation

The star formation model in M16 should be mentioned here, since the free parameters in the model will be calibrated statistically in this work. Our model assumes that gas undergoes shock heating and forms a quasi-static hot halo when it is accreted by the host dark matter halo. The gas can cool and form a cold disk in the central region, which then becomes fuel for star formation. The gas in the hot halo and the cold disk is labelled as hot and cold gas respectively. Following the disk stability argument of Kauffmann (1996), our model assumes that gas can only form stars when its mass is greater than the critical mass

[TABLE]

where $V_{\text{max}}$ is the maximum circular velocity of the host halo. The disk scale radius $r_{\text{disk}}$ is defined by

[TABLE]

where $R_{\text{vir}}$ is the virial radius of the host halo, and $\lambda$ is the spin parameter defined by Bullock et al. (2001). Then, the mass of new formed stars can be calculated from

[TABLE]

where $t_{\text{dyn, disk}}=r_{\text{disk}}/V_{\text{max}}$ is the dynamical time of the disk, $m_{\text{gas}}$ is the mass of cold gas and $\Delta t$ is the timestep. In the model, the normalisation of the critical mass $\Sigma_{\text{SF}}$ and the star formation efficiency $\alpha_{\text{SF}}$ are the two free parameters. Their preferred values will be discussed in Section 5.

2.2 Updates to Meraxes

2.2.1 Supernova Feedback

We update the supernova feedback model with a different treatment of supernova energy, and a different parametrisation of mass loading factor and energy coupling efficiency. Our original supernova model in M16 is a modified version of Guo et al. (2011), taking into account the high cadence of our halo merger trees. As mentioned in Section 2.1.1, our model galaxies have hot and cold gas components, and the effect of supernova feedback is to transfer the gas in the cold disk to the hot halo. The amount of mass that is reheated by supernova can be calculated by

[TABLE]

with

[TABLE]

where $\eta$ is the mass loading factor, $\Delta m_{\text{new}}$ is the mass of new formed stars, $\Delta E_{\text{SN}}$ is the supernova energy that is injected into the interstellar medium (ISM), and $V_{\text{vir}}$ is the virial velocity of the friends-of-friends group. If the amount of reheated mass is $\Delta m$, the energy increase of the hot halo is $\Delta E=\frac{1}{2}\Delta mV_{\text{vir}}^{2}$ after virialisation. This model first estimates the reheated mass by the mass loading factor argument, and reduces the mass if the energy injected by supernova is smaller than the underlying energy increase of the hot halo. Moreover, if $\Delta E_{\text{SN}}\geq\Delta E_{\text{hot}}$, materials can be further ejected from the hot halo. The amount of ejected mass is given by

[TABLE]

The ejected mass is subtracted from the hot gas and put into a separated component.

The injected supernova energy $\Delta E_{\text{SN}}$ plays a important role in the model described above. This quantity is given by

[TABLE]

where $\epsilon$ is the energy coupling efficiency, $t$ is the simulation time, $\Delta t$ is the timestep, $(d\varepsilon/d\tau)d\tau$ is the energy released by type-II supernova from stars with age between $\tau$ to $\tau+d\tau$ per unit mass of stellar population, and $\psi_{\text{}}(t)$ is the star formation rate as a function of the simulation time. The second term on the right hand side of Equation (7) is the total energy released by type-II supernova during a snapshot. M16 uses an analytic fit of star lifetime and an initial mass function to estimate $d\varepsilon/d\tau$, while in this work, we generate $d\varepsilon/d\tau$ using starburst99 (Leitherer et al., 1999; Vázquez & Leitherer, 2005; Leitherer et al., 2010, 2014) with metallicity dependence, assuming a Kroupa (2002) initial mass function (IMF). This treatment provides more reasonable and self-consistent estimates of the supernova energy, and can be generalised to other stellar evolutionary libraries (e.g Saitoh, 2017; Ritter et al., 2018). A similar approach has already been applied in the fire hydrodynamic simulations (Hopkins et al., 2014).

To evaluate the integral in Equation (7), we adopt the same method as M16. meraxes tracks the mass of new formed stars and their metals in four previous snapshots and assumes that they are formed by a single burst in the middle of each corresponding snapshot. Stars formed in earlier snapshots have ages greater than 55 Myr, and typically do not end with a type-II supernova. To tackle the metallicity dependence, we interpolate the table of $d\varepsilon/d\tau$ from starburst99 to a grid in a range of $Z=0.001\sim 0.040$ with resolution $\Delta Z=0.001$, and apply nearest interpolation based on the grid for each starburst.

Since supernova energy is released by stars formed in current and several previous snapshots, $\Delta m_{\text{new}}$ in Equation (4) should have contributions from these stars. This quantity is computed by

[TABLE]

In other words, we use the average star formation history weighted by the supernova energy to calculate $\Delta m_{\text{new}}$. If we assume constant canonical energy for every type-II supernova explosion, the above equation is equivalent to the number-weighted expression given by Equation (16) in M16.

The remaining parameters in the supernova feedback model are the mass loading factor $\eta$ and energy coupling efficiency $\epsilon$. In this work, we adopt different parametrisations from M16. They are given by

[TABLE]

[TABLE]

where $V_{\text{max}}$ is the maximum circular velocity. We force the maximum of $\epsilon$ to be unity due to energy conservation. Muratov et al. (2015) originally obtained a broken power law for the mass loading factor. Their study is based on model galaxies in the fire simulations (Hopkins et al., 2014). This form is subsequently implemented in several semi-analytic models (Hirschmann et al., 2016; Cora et al., 2018; Lagos et al., 2018). The implementation of this form in the present work is primarily motivated by its impact on the metallicity, which is an input of galaxy SEDs. Hirschmann et al. (2016) tested eight different supernova feedback schemes in their semi-analytic model, and found that only explicit redshift-dependent models can lead to evolution of the mass metallicity relation. Collacchioni et al. (2018) demonstrated that a steeper slope of the redshift dependence can result in stronger evolution of the mass metallicity relation using the semi-analytic model of Cora et al. (2018). In this work, we set $\alpha_{\text{reheat}}=2$ according to the optimisation result in Cora et al. (2018), assume no redshift dependence on the energy coupling efficiency (i.e. $\alpha_{\text{eject}}=0$) and leave $\eta_{0}$ and $\epsilon_{0}$ as free parameters.

2.2.2 Mass recycling and metal enrichment

We also apply starburst99 to the mass recycling and metal enrichment. The mass of materials that are produced by type-II supernova and released into the ISM can be obtained by

[TABLE]

where $(dy/d\tau)d\tau$ is the mass produced by type-II supernova from stars with age $\tau$ to $\tau+d\tau$ per unit mass of stellar population, i.e. the yield. This quantity depends on the IMF and varies with different elements. We generate the table of $dy/d\tau$ using starburst99, including metallicity dependence and assuming a Kroupa (2002) IMF. In the present work, we only consider two cases, i.e. the yield of all elements and the yield of all metal elements. The former gives the amount of recycled mass, while the latter introduces metal enrichment. The evaluation of the integral in Equation (11) and the treatment of metallicity dependence follow the same approach as the calculation of the total supernova energy. All recycled materials are added into the cold gas component. This age-dependent mass recycling scheme is introduced due to the short timestep of the halo merger tree, and is more realistic than the commonly adopted constant recycling fraction and yield, particularly at high redshift.

2.2.3 Reincorporation

The ejected gas component mentioned in the previous section can be transferred back to the hot gas halo. We employ the reincorporation model proposed by Henriques et al. (2013)

[TABLE]

where $m_{\text{eject}}$ is the mass of ejected gas, and $M_{\text{vir}}$ is the virial mass of the friends-of-friends group. We also force the reincorporation time scale to be smaller than the halo dynamic time. The statistical analysis of Henriques et al. (2013) indicates that this model provides better fit of the stellar mass functions against observations at $z\leq 3$. We set $\gamma=18\,\text{Gyr}$ as suggested by Henriques et al. (2013). This model is also implemented in Hirschmann et al. (2016), Cora et al. (2018) and Lagos et al. (2018). We note that with this choice of $\gamma$, the reincorporation time scale equals the forced upper limit, i.e. the halo dynamical time, for $M_{\text{vir}}\lesssim 10^{12}M_{\odot}$ at the redshift range of interest in this study. At $z\sim 4-7$, the halo dynamical time is $\sim 100$ Myr. Therefore, reincorporation is more efficient in our model relative to Henriques et al. (2013) at these redshifts, in particular for low mass halos. This behaviour is very different from the original model proposed by Henriques et al. (2013), and may impact on the predicted stellar mass functions at lower redshifts. We defer the exploration of this effect to later works.

3 Dust model and synthetic spectral energy distributions

3.1 Dust Model

We implement the dust model proposed by (Charlot & Fall, 2000). The transmission function due to the ISM is expressed by

[TABLE]

This model takes into account the relative stars-dust geometry of different stellar populations. Photons emitted by young stars are absorbed by an additional component due to the surrounding molecular cloud where the stars form. The birthcloud is assumed to have lifetime $t_{\text{BC}}$, and for stars whose age is older than $t_{\text{BC}}$, their starlight is only absorbed by the diffuse ISM dust. We fix $t_{\text{BC}}=10$ Myr according to previous studies (Charlot & Fall, 2000; da Cunha et al., 2008). The attenuation due to the birth cloud and diffuse ISM dust is described by their optical depths $\tau^{\text{BC}}_{\lambda}$ and $\tau^{\text{ISM}}_{\lambda}$ respectively, which should vary with different galaxies. In this study, we explore three different parametrisations, linked to star formation rate (SFR), dust-to-gas (DTG) ratio and gas column density (GCD). We name them as M-SFR, M-DTG and M-GCD respectively. In general, these properties are indirectly related to the dust. One dust production channel is from the ejecta of supernova (e.g. Dayal & Ferrara, 2018), which is proportional to the SFR. Dust is also mixed with gas. Accordingly, they are expected to have similar properties. We will see that M-DTG and M-GCD have similar results since they primarily depend on the gas density.

3.1.1 Star formation rate model

The dependence of the dust optical depths on SFR is motivated by observations of the CMRs at high redshifts, i.e the relation between UV continuum slope and UV magnitude. These observations suggest that more UV luminous galaxies have redder UV continuum slopes (Finkelstein et al., 2012; Bouwens et al., 2014; Rogers et al., 2014). Since brighter galaxies correspond to higher SFR, one could expect that SFR and dust content are positively correlated. Similar trends have been found in low redshift studies (e.g. da Cunha et al., 2010; Qin et al., 2019). Hence, we assume the following parameterisation

[TABLE]

where $\tau^{\text{ISM}}_{\rm SFR}$, $\tau^{\text{ISM}}_{\rm SFR}$, $\gamma_{\text{SFR}}$, $a$ and $n$ are free parameters. Yung et al. (2019) also use a parametric model to calculate dust attenuation in their semi-analytic model. They adjusted the normalisation of the optical depth to fit the observed UV LFs at individual redshifts. Their results indicate that the normalisation depends on redshift and the trend can be fit by an exponential function. Therefore, for all the three parametrisations proposed in this work, we also include an exponential redshift dependence factor to fit the model against multiple redshifts.

3.1.2 Dust-to-gas ratio model

In the literature, dust optical depths are often linked to the gas column density, which is then converted to the dust column density using the dust-to-gas (DTG) ratio (De Lucia & Blaizot, 2007; Guo et al., 2011; Somerville et al., 2012; Yung et al., 2019). In this model, optical depths are expressed by

[TABLE]

where $Z_{\text{cold}}$ is the metallicity of cold gas, $m_{\text{cold}}$ is the mass of cold gas, and $r_{\text{disk}}$ is the disk scale radius defined in Equation (2). We adopt the solar metallicity as $Z_{\odot}=0.02$. Free parameters are $\tau^{\text{ISM}}_{\rm DTG}$, $\tau^{\text{BC}}_{\rm DTG}$, $\gamma_{\text{DTG}}$, $a$ and $n$.

3.1.3 Gas column density model

We propose an additional gas mass related dust model, which is independent of the metallicity. In M16 and this work, when metals are produced by supernova explosions, we assume that they are first fully mixed with cold gas, and then ejected into the hot gas reservoir. In reality, since the materials produced by supernova have quite different initial velocities from the surrounding gas, the mixing may take some time. Thus, we provide a metallicity independent parametrisation of the dust optical depths

[TABLE]

There are also five free parameters in this model, i.e. $\tau^{\text{ISM}}_{\rm GCD}$, $\tau^{\text{BC}}_{\rm GCD}$, $\gamma_{\text{GCD}}$, $a$ and $n$. This model includes a power law scaling on the cold gas mass, unlike the M-DTG model, where the scaling is on metallicity.

3.2 Synthetic spectral energy distributions

The computation of galaxy spectral energy distributions (SEDs) follows standard stellar population synthesis. The luminosity of a galaxy at time $t$ can be obtained by

[TABLE]

where $\tau$ is the stellar age, $\psi(t-\tau,Z)\,d\tau dZ$ is the mass of stars formed at $t-\tau$ with an age between $\tau$ to $\tau+d\tau$ and metallicity between $Z$ to $Z+dZ$, $S_{\lambda}(\tau,Z)$ is the luminosity of a simple stellar population (SSP) per unit mass, and $T_{\lambda}(\tau)$ is the transmission function of the ISM described in the previous subsection. We generate $S_{\lambda}(\tau,Z)$ using starburst99 (Leitherer et al., 1999; Vázquez & Leitherer, 2005; Leitherer et al., 2010, 2014), assuming a metallicity range from $Z=0.001$ to $Z=0.040$ and a Kroupa (2002) IMF. Nebular continuum emssions are also added using starburst99. To compute UV magnitudes, we apply a tophat filter centred at $\lambda=1600\text{\AA}$ with width $100\text{\AA}$. UV slopes are obtained by a linear fit in the logarithmic flux space using the ten windows proposed by Calzetti et al. (1994). However, for computational speed, we only choose five of them (including the longest wavelength window) for on-the-fly calibrations. The selected windows are given in Table 1. The median errors from this treatment are negligible in the range of the observed CMRs.

We also make a numeric approximation in order to accelerate the speed of evaluating Equation (24). We first compute the intrinsic luminosity in necessary filters. The dust transmission is then applied to the luminosity of the filters using the central wavelength instead of the full SEDs. This approximation is found to have a negligible effect on the results, since all filters used in this work have a simple shape and are relatively narrow.

4 Calibration

An essential part of this work is to determine the free parameters in both the galaxy formation and dust attenuation models introduced in the previous sections. We carry out a Bayesian analysis on these parameters, and use observed UV LFs and CMRs at $z\sim 4-7$ as constraints.

A key goal of a Bayesian analysis is to estimate the posterior distribution of model parameters, which is non-trivial for high dimensional spaces. Kampakoglou et al. (2008) and Henriques et al. (2009) first applied the Markov chain Monte Carlo (MCMC) method to sample the parameter space of semi-analytic models. This approach has been implemented by several subsequent studies (Henriques et al., 2013; Mutch et al., 2013; Henriques et al., 2015). However, the MCMC method has several drawbacks. Firstly, it requires additional evaluations of the model to ensure the final sample reaches a stationary distribution, and it is generally difficult to determine whether a Monte Carlo chain has fully converged (see Cowles & Carlin, 1996, for a review). Moreover, without special treatments, MCMC samplers can encounter difficulties in approaching a stationary distribution when the parameter space contains isolated modes (which is the case in this study), since random walkers can be trapped by a local minimum and fail to jump to other modes (e.g Neal, 1996). A possible improvement to handle multimodal distributions for MCMC methods can be found in Earl & Deem (2005).

In this work, in order to achieve higher sampling efficiency and obtain more stable results on multimodal parameter spaces, we utilise the multimodal nested sampling introduced by Feroz et al. (2009) to estimate the posterior distributions. This algorithm is found to be a competitive alternative to MCMC methods, and addresses the issues mentioned above to some extent. The nested sampling was designed to evaluate the Bayesian evidence (Skilling, 2004). However, the output samples produced by the algorithm can also be used to estimate posterior distributions, which is equivalent to the MCMC method. In difference from MCMC methods, no burn-in phase is required in this algorithm. The stopping criterion of the nested sampling is based on an estimated error of the resulting value of the Bayesian evidence, which is also proposed by Skilling (2004). The sampling efficiency of the original algorithm is improved by Feroz et al. (2009), who use the information of existing sample points to approximate the iso-likelihood surfaces in the parameter space as hyper ellipsoids (see also Mukherjee et al., 2006). Secondly, the algorithm includes a special treatment for multimodal problems. Again using the information of existing sample points, it applies a clustering algorithm to detect multiple modes and splits the parameter space (see also Shaw et al., 2007). This approach has been tested against toy models which contain several equally-high peaks, and is found to have good performance. The reader is referred to Feroz & Hobson (2008), Feroz et al. (2009) and references therein for a detailed description of the algorithm. A comparison between the nested sampling and the MCMC method can found be in Speagle (2019).

The Bayesian posterior distribution is comprised of the likelihood and prior distributions of each free model parameter. We construct the log-likelihood as

[TABLE]

Observational data of LFs ($n^{\rm obs}_{i}$, $\sigma^{2}_{\rm LF,i}$) and CMRs ($\beta^{\rm obs}_{i}$, $\sigma^{2}_{\rm CMR,i}$) are taken from Bouwens et al. (2015) and the biweight mean measurements of Bouwens et al. (2014) respectively. The LFs are defined by the co-moving number density. We convert the dimensionless Hubble constant from $h=0.7$ to $h=0.678$ for these observations in order to be consistent with our model. Due to the limited size of the simulation box, the model is unable to probe the full range of the LFs and CMRs. Therefore, for each LF and CMR bin, we use the observed LF to estimate an expected number of galaxies in the simulation box, and drop the bin if the number is less than five and twenty for the LF and CMR respectively.

The model parameters are from both the semi-analytic model and the dust relations introduced in Section 3.1. We focus on four galaxy formation parameters: the star formation efficiency $\alpha_{\text{SF}}$, normalisation of the critical mass $\Sigma_{\text{SF}}$, mass loading factor $\eta_{0}$, and supernova energy coupling efficiency $\epsilon_{0}$. Their prior distributions are chosen to be uniform in logarithmic space. Three different dust models were introduced in Section 3.1. Each of them has five free parameters. We adopt uniform priors in linear space for them.

The prior ranges of these model parameters are used in the initialisation of nested sampling, and they are chosen based on several experiments. We first run the sampler in a very large parameter space and find the high probability regions. We then shrink the prior ranges accordingly, keeping the posterior distribution at the bounds negligible compared with the high probability regions. There is an exception for the mass loading factor $\eta_{0}$ in the M-SFR model, which is found to have no upper limit. However, this will not affect our main results, since the energy coupling efficiency $\epsilon_{0}$ puts physical upper limit on the strength of the supernova feedback, and this parameter is constrained. This approach of choosing the prior ranges allows the sampler to spend more time on the high probability regions and improves the sampling efficiency. A summary of all model parameters and their prior ranges can be found in Table 2.

In practice, we utilise a modified verion of the open source Python package nestle111https://github.com/kbarbary/nestle. See https://github.com/smutch/nestle for the modified version., which implements the algorithm, and couple it with the meraxes Python interface mhysa (Mutch, in prep.). We set the number of active points to be 300 for the sampler. The stop criterion follows the remaining Bayesian evidence approach suggested by Skilling (2004). The algorithm terminates when the logarithmic change due to the remaining Bayesian evidence is below one, and the convergence requires evaluating the model for 50,000 - 100,000 times.

5 fitting results

For the three different dust models, we obtain $5,000$ - $6,000$ sample points from the nested sampling algorithm, which describe the posterior distributions of both galaxy and dust parameters. The point that has highest value of the posterior distribution is chosen to be the best-fit result. The best-fit parameter values are listed in Table 2, and the corresponding LFs and CMRs are shown in Figure 1 for each dust model. The three models all fit the observational data extremely well. In Figures 10, 11 and 12, we show the posterior distributions of M-SFR, M-DTG and M-GCD respectively. In plotting these figures, we adopt a similar approach with Henriques et al. (2009) and Henriques et al. (2013), i.e. using contours to show the marginalised distributions and colours to indicate the values of the whole posterior distributions. For the M-SFR model, it can be seen from Figure 2 that the marginalised distribution of the mass loading factor $\eta_{0}$ extends to large values, which means that this parameter is less constrained. On the other hand, all parameters are well constrained for the other two models.

An interesting finding is that the derived galaxy formation parameters preferred by these three dust models are quite different. Figure 2 illustrates a comparison of the marginalised distributions for the four galaxy formation parameters. We found that M-DTG and M-GCD suggest similar mass loading factor and supernova energy coupling efficiency. However, M-DTG shows evidence of a more active star formation scenario, with higher star formation efficiency and lower normalisation of the critical mass. For M-SFR, the marginalised distribution of the star formation efficiency overlaps with that of M-DTG. However, M-SFR requires much smaller supernova energy coupling efficiency. Moreover, differences can also be found in the two parameter correlations between the galaxy formation parameters for the three different models. For instance, in the third row and first column of Figure 12, M-GCD shows a strong correlation between the star formation efficiency $\alpha_{\text{SF}}$ and the mass loading factor $\eta_{0}$. However, this correlation cannot be found the in the other two models. The variation among the posterior distributions of the three models implies that these free parameters fit the observational data in a very complex way and the constraints on them depend on the assumptions used to model the dust attenuation.

By comparing the posterior distributions of the three dust models, we find similar correlations among the parameters of the supernova feedback, galaxy property scaling of the dust relation and reddening slope. We demonstrate this in Figure 3. It can be seen that the supernova energy coupling efficiency is positively and inversely correlated with $\gamma_{\text{SFR,\,DTG,\,GCD}}$ and the reddening slope $n$, respectively. While the trends are the weakest for the M-GCD model, the correlation between $\gamma_{\text{SFR,\,DTG,\,GCD}}$ and $n$ is obvious for all the three models. Similar correlations are also found for the mass loading factor $\eta_{0}$. The reader is referred to Figures 10, 11 and 12 for the two parameter marginalised distributions of all model parameters. The dependence between the galaxy formation and dust parameters is important, since it suggests that the observations can put constraints on intrinsic or dust-unattenauted galaxy properties despite the degrees of freedom in the dust models.

In order to understand the correlations mentioned above, we plot the intrinsic LFs and CMRs for the three best-fit models in Figure 1. They are shown as dashed lines. It can be seen that the LFs of the best-fit M-SFR is roughly a factor of two higher than for the other two models. We note that the intrinsic luminosity functions are more sensitive to feedback processes rather than the star formation law due to self-regulation (e.g. Schaye et al., 2010; Lagos et al., 2011). The supernova coupling efficiency of the best-fit M-SFR model is much smaller than the other two models, which is likely to be the main reason for the difference in the intrinsic LFs, irrespective of those star formation parameters. We examine the effect of supernova feedback in Figure 4. For the upper panel, we vary the mass loading factor $\eta_{0}$ at fixed supernova energy coupling efficiency $\epsilon_{0}$ for the best-fit M-DTG model, and compare the resulting LFs with the best-fit results. The y-axis shows the ratio of the logarithmic intrinsic LFs. We find that the number density at fixed UV magnitude decreases with increasing $\eta_{0}$. Since the energy coupling efficiency puts an upper limit on the reheated mass (see Equation 4), the change in the LFs is smaller at larger $\eta_{0}$. The results of varying $\epsilon_{0}$ at fixed $\eta_{0}$ are shown in the lower panel of Figure 4. While higher $\epsilon_{0}$ decreases the LFs, the effect is found to be more significant at the bright end. Since the energy coupling efficiency is assumed to scale as a power law of the maximum circular velocity (see Equation 10), the efficiency can easily reach the maximum value of unity for small galaxies. The median $M_{\text{UV}}$ - $V_{\text{max}}$ relation of the best-fit M-DTG model indicates that the energy coupling efficiency becomes unity at an intrinsic magnitude $M_{\text{UV}}\sim-18$ with $\epsilon_{0}=1.5$. Therefore, the major effect of increasing $\epsilon_{0}$ is to allow more gas to be reheated in galaxies hosted by more massive halos. This explains why this parameter has larger impact at the bright end of the intrinsic LFs. Overall, the above discussion implies that supernova feedback plays an important role in regulating the intrinsic LFs.

We next investigate the effect of the galaxy property scaling of the dust relation. Figure 5 shows the resulting dust-attenuated LFs and CMRs when $\gamma_{\text{SFR,DTG,GCD}}$ is changed by a factor of two. It can be seen that the shape of both the LFs and CMRs are quite sensitive to this parameter. Furthermore, since the dust optical depths are assumed to depend on different galaxy properties, this parameter changes the shape of the LFs and CMRs in different ways.

Combining the discussions of the supernova feedback parameters and $\gamma_{\text{SFR,DTG,GCD}}$ above, we provide an explanation of the correlations between the model parameters in Figure 3. In our dust models, the effective UV optical depth is a function of the galaxy property scaling $\gamma_{\text{SFR,DTG,GCD}}$ and the optical depth normalisations of both the ISM and BC. The galaxy property scaling has a direct impact on the shape of the dust-attenuated LFs and CMRs, and this single parameter is required to satisfy two shapes. Hence, the effective UV optical depth is very sensitive to $\gamma_{\text{SFR,DTG,GCD}}$. The effective optical depth should be degenerate with the intrinsic UV LFs, which are primarily controlled by the supernova feedback parameters. These imply that both $\eta_{0}$ and $\epsilon_{0}$ should be correlated with $\gamma_{\text{SFR,DTG,GCD}}$. On the other hand, the observed UV continuum slope $\beta$ depends on the reddening curve, which is assumed to be a power law of wavelength. Since there is a natural degeneracy between the power law slope and the normalisation, the reddening slope $n$ should be degenerate with the effective UV optical depth, and therefore $\gamma_{\text{SFR,DTG,GCD}}$. The dependence between the supernova feedback parameters and the reddening slope can be derived from the above two correlations. Figure 5 shows that the galaxy property scaling on the SFR changes the dust-attenuated LFs and CMRs differently from the other two models, which may explain why the best-fit M-SFR model requires much weaker supernova feedback. This also explains the shallower reddening slope required by the best-fit M-SFR model.

In addition, we contrast the best-fit models for M-DTG and M-GCD. Whilst their intrinsic LFs and CMRs are almost the same, a difference is found in the metallicity. Figure 6 illustrates the cold gas metallicity at two stellar mass bins as a function of redshift for all best-fit models. It is clear that the cold gas is more metal enriched in the best-fit M-DTG than in M-GCD. We identify that the normalisation of the critical mass $\Sigma_{\text{SF}}$ is the primary driver for the difference since both best-fit models have similar parameters of supernova feedback. To confirm this, we run a model variant, setting $\Sigma_{\text{SF}}=0.2$, with other parameters being the same with the best-fit M-DTG. The resulting metallicity is also shown in Figure 6, which is similar to that of the best-fit M-GCD. This finding is unsurprising since the star formation law affects gas fraction and therefore metallicity (e.g. Schaye et al., 2010; Lagos et al., 2011). In addition, from Figure 6, it is worth noting that the metallicity evolves with redshift in our model, with higher metallicity at lower redshifts. This is expected due to the explicit redshift dependence on the mass loading factor, which is motivated by previous studies (Muratov et al., 2015; Hirschmann et al., 2016; Collacchioni et al., 2018).

6 Infrared excess to UV continuum slope relations

As demonstrated in previous sections, by simultaneously fitting our galaxy formation and dust models to the observed UV LFs and CMRs, we are able to obtain constraints on both the dust attenuation in the UV band and the reddening. This allows estimates of the infrared luminosity $F_{\text{IR}}$ and therefore the infrared excess (IRX) using energy balance arguments, i.e.

[TABLE]

We compute the IRX for galaxies in the best-fit models of M-SFR, M-DTG and M-GCD. The resulting IRX - $\beta$ relations for galaxies with stellar mass greater than $10^{8}\,\text{M}_{\odot}$ are shown in Figure 7 with several observations for comparison. Taking into account intrinsic scatter in the relations, our results cover the observations from Koprowski et al. (2018), who performed a stacking analysis of 4209 Lyman-break galaxies (LBGs) at $3\lesssim z\lesssim 5$, and individual detection from Barisic et al. (2017). We also compare our predictions with the relations calibrated by Meurer et al. (1999) using local starburst galaxies. The Meurer et al. (1999) relation is frequently used to correct dust extinction in both observational and theoretical studies at high redshifts (e.g. Duncan et al., 2014; Bouwens et al., 2015; Mason et al., 2015; Liu et al., 2016; Harikane et al., 2018). It can be seen from Figure 7 that the best-fit M-SFR predicts higher IRX than the Meurer et al. (1999) relation at fixed $\beta$, while the other two best-fit models suggest lower IRX. Thus, our models indicate dust extinction that differs from the Meurer et al. (1999) relation, which implies that a direct application of the relation at high redshifts may lead to systematic errors on the dust corrections. We will discuss the resulting uncertainties on estimations of the cosmic SFRD in Section 7.

6.1 Reddening slope

The best-fit models for M-SFR, M-DTG and M-GCD have quite different reddening slopes $n$, which can be read from Table 2. The best-fit M-SFR model has the shallowest slope of $n=-0.3$, while much steeper slopes are found for the best-fit M-DTG and M-GCD, with $n=-1.6$ and $-1.3$ respectively. This difference is directly reflected on the IRX - $\beta$ plane. In Figure 6, the best-fit M-DTG and M-GCD show a shallower IRX - $\beta$ relation at redder UV slope regime. Similar disagreements can also be found from other studies. For example, Cullen et al. (2017) post-processed the outputs of the FiBY hydrodynamic simulation (Johnson et al., 2013; Paardekooper et al., 2015). They propose a similar dust model to the present work, linking the dust optical depths to the logarithmic stellar mass. The free parameters in their model are adjusted to fit the observed LFs and CMRs from Rogers et al. (2014) at $z\sim 5$. They suggest $n=-0.55^{+0.25}_{-0.15}$. We plot their results as solid red lines in Figure 7, which is more consistent with the best-fit M-SFR than our other models. On the other hand, Mancini et al. (2016) also post-processed a hydrodynamic simulation by Maio et al. (2010), and coupled it with an dust evolution model. Their results reproduce the observed LFs of Bouwens et al. (2015) and CMRs of Bouwens et al. (2014) at $z\sim 5-8$ when using a SMC-like extinction curve. The slope of the SMC-like extinction curve is steeper, and is similar to those of the best-fit M-DTG and M-GCD. Since all our models can well reproduce observed LFs and CMRs, we cannot draw any firm conclusions on the reddening slope. Instead, we treat this as systematic uncertainties arising due to different assumptions in the dust models.

6.2 Intrinsic scatter

At $z\gtrsim 3$, observations show considerable scatter in the IRX - $\beta$ plane (e.g. Capak et al., 2015; Álvarez-Márquez et al., 2016; Bouwens et al., 2016; Barisic et al., 2017; Fudamoto et al., 2017; Koprowski et al., 2018), which might be explained by the large intrinsic scatter in our predicted relations. Hence, it is instructive to examine the main drivers of the scatter. We first notice that from Figure 7, low IRX galaxies vanish in the best-fit M-SFR model. This is due to the nature of our star formation prescription (see Section 2.1.1). The SFR of galaxies whose cold gas mass is below the critical mass is zero. Accordingly, in the M-SFR model, the dust optical depths of these galaxies are also zero, which results in the disappearance of the IRX. This unrealistic feature shows the limitations of this model.

In left and right panels of Figure 8, we illustrate the IRX - $\beta$ relations at $z\sim 5$ as functions of stellar mass and specific star formation rate (sSFR) respectively. The relations at other redshifts show similar trends. For the stellar mass case, we see that massive galaxies form a tight correlation between IRX and $\beta$ in the high IRX and red $\beta$ regions. The trend that more massive galaxies have higher IRX is also observed by Álvarez-Márquez et al. (2016) and Fudamoto et al. (2017). However, we also find several larger stellar mass galaxies which have lower IRX and redder $\beta$. They might explain some of the outliers individually detected by Barisic et al. (2017), as shown in Figure 8. On the other hand, the right panels show that the scatter of the IRX - $\beta$ relation is tightly correlated with sSFR. At fixed IRX, redder galaxies typically have lower sSFR. This trend is consistent with other theoretical studies (Popping et al., 2017b; Safarzadeh et al., 2017; Narayanan et al., 2018; Cousin et al., 2019a). In addition, it is worth noting that although the dust optical depths are related to different galaxy properties for the three dust models, we find similar dependence of the scatter in the IRX - $\beta$ plane on both stellar mass and sSFR.

7 Cosmic star formation rate density

Dust corrections are typically required for the conversion between the UV luminosity and the SFR. As mentioned, in high redshift observational studies of SFR, the Meurer et al. (1999) relation is widely used, though it is calibrated against local galaxies. The previous section has shown that the dust extinction predicted by our models, which reproduce both LFs and CMRs at $z\sim 4-7$, is rather different from the Meurer et al. (1999) relation. In principle, we could derive similar relations based on our results to be used by other studies to perform the dust corrections. However, by using such relations, we should be able to recover the SFR functions of our models given the LFs. Therefore, we directly present the predicted SFRs. Furthermore, the difference among the three dust models allows us to estimate the systematic uncertainties in the observed SFRs.

Figure 9 illustrates the predicted cosmic star formation rate density (SFRD) for the best-fit models of M-SFR, M-DTG and M-GCD. Their values are listed in Table 3. We compare the results with Bouwens et al. (2015), whose estimations are based on the CMRs of Bouwens et al. (2014) and the Meurer et al. (1999) relation. Since the results of Bouwens et al. (2015) and our models use the same observational information, the comparison between them quantifies the systematic errors of using the Meurer et al. (1999) relation with respect to correcting the dust extinction. We note that all our models suggest bluer intrinsic UV continuum slopes than the one used in Meurer et al. (1999) as shown in Figure 1. Figure 7 also illustrates that the best-fit results of M-DTG and M-GCD have shallower IRX - $\beta$ relations. Thus, compared with the Meurer et al. (1999) relation, the dust extinction in these two models is stronger for bluer galaxies but weaker for redder galaxies. On the other hand, the dust attenuation is stronger for all galaxies in the best-fit M-SFR. It can be seen from Figure 9 that the cosmic SFRD of the best-fit M-DTG and M-GCD are consistent with those of Bouwens et al. (2015), while the results of the best-fit M-SFR is roughly a factor of two higher. We also compare our results with Driver et al. (2018). Their dust-corrected SFRs are derived from the energy balance SED-fitting code magphys (da Cunha et al., 2008). Better consistency is found between their measurements and our best-fit models of M-DTG and M-GCD, given the size of the error bars on those points. Overall, Figure 9 suggests that uncertainty in the dust relations introduces at least a factor of two systematic error into the inferred cosmic SFRD at $z\gtrsim 6$.

8 Summary

This work investigates the IRX - $\beta$ relation and cosmic SFRD at $z\sim 4-7$ by combining the the meraxes semi-analytic galaxy formation model (Mutch et al., 2016; Qin et al., 2017) and the Charlot & Fall (2000) dust attenuation model. The supernova feedback model of meraxes is updated using results from previous studies (Muratov et al., 2015; Hirschmann et al., 2016; Cora et al., 2018), which aim to reproduce the redshift evolution of the mass metallicity relation. We introduce three different parametrisations of the dust optical depths, which are related to the star formation rate (M-SFR), dust-to-gas ratio (M-DTG) and gas column density (M-GCD) respectively. These lead to five free parameters in each dust model in additional to those in meraxes.

The determinations on not only the dust parameters but also the meraxes free parameters constitute the primary part of this work. For galaxy formation parameters, we focus on the star formation efficiency, critical mass, mass loading factor and supernova coupling efficiency. We adopt a Bayesian approach, calibrating these parameters against the UV luminosity functions (LFs) of Bouwens et al. (2015) and colour magnitude relations (CMRs) of Bouwens et al. (2014) at $z\sim 4-7$. The posterior distribution of these parameters is estimated using multimodal nested sampling (Feroz et al., 2009). We find that these observations can be fit extremely well by all the three dust models. However, the preferred parameter ranges are quite different among the three dust models. Our analysis indicates that the combination of the LFs and CMRs can put strong constraints on a given dust attenuation model, since the model is required to reproduce the shape of both observations. The differences in our results are due to the different assumptions of the dust models, which results in different relations between UV dust attenuation and intrinsic UV magnitude.

We then demonstrate the predictions of our calibration results. Using energy balance arguments, we estimate the IRX for each model galaxy. We find that the predicted IRX - $\beta$ relations are quite different from the Meurer et al. (1999) relation, and contain large intrinsic scatter, which might explain the current discrepancy among several high redshift observations (e.g. Capak et al., 2015; Álvarez-Márquez et al., 2016; Bouwens et al., 2016; Barisic et al., 2017; Fudamoto et al., 2017; Koprowski et al., 2018). We also confirm the correlation between the intrinsic scatter and sSFR. This finding is consistent with other theoretical studies (Popping et al., 2017b; Safarzadeh et al., 2017; Narayanan et al., 2018; Cousin et al., 2019a). Secondly, we present model predictions for the cosmic SFRD, and compare these with the observations of Bouwens et al. (2015) and Driver et al. (2018). The difference among the three dust models implies at least a factor of two systematic uncertainty in the observed SFRD when corrected using the Meurer IRX - $\beta$ relation.

This work has simultaneously constrained the free parameters of a semi-analytic galaxy formation model and additional dust parameters using observations of UV properties. Within a Bayesian framework, our approach establishes a more direct connection between the model and observations despite the complexity. This approach is particularly useful for studies at high redshifts where UV properties are the most robust observables. This work could be further improved by explicitly modelling the dust evolution (e.g Mancini et al., 2016; Popping et al., 2017a; Dayal & Ferrara, 2018), which might reduce the systematic uncertainties due to different assumptions in the dust models. Additional free parameters (e.g. the time scale of dust growth) in such dust evolution models could also be constrained using our methodology.

acknowledgement

This research was supported by the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013. This work was performed on the OzSTAR national facility at Swinburne University of Technology. OzSTAR is funded by Swinburne University of Technology and the National Collaborative Research Infrastructure Strategy (NCRIS). We thank the referee for providing constructive comments to improve the quailty of this paper.

In addition to those already been mentioned in the paper, we also acknowledge the use of the following software: astropy 111http://www.astropy.org (Astropy Collaboration et al., 2013, 2018), corner (Foreman-Mackey, 2016), cython (Behnel et al., 2011), ipython (Perez & Granger, 2007), matplotlib (Hunter, 2007), numpy (van der Walt et al., 2011), pandas (McKinney, 2010), seaborn222https://github.com/mwaskom/seaborn and scipy (Jones et al., 2001).

Appendix A Posterior distributions

This appendix illustrates the posterior distributions of M-SFR, M-DTG and M-GCD in Figures 10, 11 and 12, respectively. These results are obtained using the methodology introduced in Section 4.