Predictions for deep galaxy surveys with JWST from $\Lambda$CDM

TL;DR

This paper uses a physical galaxy formation model within the $\\Lambda$CDM framework to predict galaxy observations with JWST, including luminosity functions, counts, and sizes across redshifts, validated against existing data.

Contribution

It provides the first comprehensive predictions for JWST galaxy surveys based on a calibrated, physically motivated galaxy formation model covering high redshifts up to 16.

Findings

Predicts ~1 galaxy per JWST field at z~11 with 10^4 s exposure.

Higher reionization redshift model predicts 5x more galaxies at z~11.

MIRI is unlikely to detect galaxies beyond z~6.

Abstract

We present predictions for the outcome of deep galaxy surveys with the $James$ $Webb$ $Space$ $Telescope$ ($JWST$) obtained from a physical model of galaxy formation in $\Lambda$CDM. We use the latest version of the GALFORM model, embedded within a new ($800$ Mpc)$^{3}$ dark matter only simulation with a halo mass resolution of $M_{\rm halo}>2\times10^{9}$ $h^{-1}$ M$_{\odot}$. For computing full UV-to-mm galaxy spectral energy distributions, including the absorption and emission of radiation by dust, we use the spectrophotometric radiative transfer code GRASIL. The model is calibrated to reproduce a broad range of observational data at $z\lesssim6$, and we show here that it can also predict evolution of the rest-frame far-UV luminosity function for $7\lesssim z\lesssim10$ which is in good agreement with observations. We make predictions for the evolution of the luminosity function from $z=16$ to $z=0$ in all broadband filters on the Near InfraRed Camera (NIRCam) and Mid InfraRed Instrument (MIRI) on $JWST$ and present the resulting galaxy number counts and redshift distributions. Our fiducial model predicts that $\sim1$ galaxy per field of view will be observable at $z\sim11$ for a $10^4$ s exposure with NIRCam. A variant model, which produces a higher redshift of reionization in better agreement with $Planck$ data, predicts number densities of observable galaxies $\sim5\times$ greater at this redshift. Similar observations with MIRI are predicted not to detect any galaxies at $z\gtrsim6$. We also make predictions for the effect of different exposure times on the redshift distributions of galaxies observable with $JWST$, and for the angular sizes of galaxies in $JWST$ bands.

Figures

Tables

Table 1. Table 1: Changes between parameter values presented in Lacey et al. ( 2016 ) and those used in this work (and discussed further in Baugh et al. in preparation). The galaxy formation parameters are listed in the bottom part of the table.

Parameter	Description	Lacey et al. (2016)	This work
Cosmological parameters		Komatsu et al. (2011)	Planck Collaboration (2016)
${\Omega }_m$	Matter density	$0.272$	$0.307$
${\Omega }_{\Lambda }$	Vacuum energy density	$0.728$	$0.693$
${\Omega }_b$	Baryon density	$0.0455$	$0.0483$
$h$	Hubble Parameter	$0.704$	$0.678$
${\sigma }_8$	Fluctuation amplitude	$0.810$	$0.829$
$N$-body simulation parameters
$M_{\mathrm{halo},\min }$	Minimum halo mass	$1.87\times {10}^{10}$ $h^{-1}$ M☉	$2.12\times {10}^9$ $h^{-1}$ M☉
Galaxy merger timescale		Jiang et al. (2008)	Simha & Cole (2017)
Galaxy formation parameters
${\alpha }_{\mathrm{ret}}$	Gas reincorporation timescale factor	$0.64$	$1.00$
${\gamma }_{\mathrm{SN}}$	Slope of SN feedback mass loading	$3.2$	$3.4$

Table 2. Table 2: Adopted values for adjustable parameters in grasil . See the text in Section  2.3 for their definitions.

Parameter	Value
$h_z/h_R$ (disc)	$0.1$
$h_z/h_R$ (burst)	$0.5$
$h_z\left(\mathrm{dust}\right)/h_z\left(\mathrm{stars}\right)$	$1$
$\eta$	$1.0$
$f_{\mathrm{cloud}}$	$0.5$
$m_{\mathrm{cloud}}/r_{\mathrm{cloud}}^2$	${10}^6$ M${}_{☉}{}^{}/(16$ pc$){}^2$
$t_{\mathrm{esc}}$	$1$ Myr
${\beta }_b$	$1.5$

Table 3. Table 3: Adopted sensitivities for JWST filters based on 10 ​ σ 10 𝜎 10\sigma point source and 10 4 superscript 10 4 10^{4}  s exposure.

Instrument	Filter	${\lambda }_{\mathrm{eff}}$ ($\mu$m)	Sensitivity ($\mu$Jy)
NIRCam	F070W	$0.70$	$20.9\times {10}^{-3}$
	F090W	$0.90$	$13.1\times {10}^{-3}$
	F115W	$1.15$	$11.8\times {10}^{-3}$
	F150W	$1.50$	$9.6\times {10}^{-3}$
	F200W	$2.00$	$7.9\times {10}^{-3}$
	F277W	$2.77$	$11.5\times {10}^{-3}$
	F356W	$3.56$	$11.1\times {10}^{-3}$
	F444W	$4.44$	$17.6\times {10}^{-3}$
MIRI	F560W	$5.6$	$0.2$
	F770W	$7.7$	$0.28$
	F1000W	$10.0$	$0.7$
	F1130W	$11.3$	$1.7$
	F1280W	$12.8$	$1.4$
	F1500W	$15.0$	$1.8$
	F1800W	$18.0$	$4.3$
	F2100W	$21.0$	$8.6$
	F2550W	$25.5$	$28.0$
Note: Adapted from https://jwst.stsci.edu/files/live/sites/jwst/files/home/science%20planning/Technical%20documents/JWST-PocketBooklet_January17.pdf

Table 4. Table 4: Adopted JWST instrument fields of view (FoV).

Instrument	FoV (arcmin2)
NIRCam	$2\times 2.2\times 2.2$
MIRI	$1.23\times 1.88$
Note: From https://jwst.stsci.edu/files/live/sites/jwst/files/home/science%20planning/Technical%20documents/JWST-PocketBooklet_January17.pdf

Equations

$$R_{\mathrm{mol}}=\frac{{\Sigma }_{\mathrm{mol}}}{{\Sigma }_{\mathrm{atom}}}={\left[\frac{P}{P_0}\right]}^{{\alpha }_{\mathrm{P}}},$$

$${\Sigma }_{\mathrm{SFR}}={\nu }_{\mathrm{SF}}{\Sigma }_{\mathrm{mol}}={\nu }_{\mathrm{SF}}{f}_{\mathrm{mol}}{\Sigma }_{\mathrm{cold}},$$

$${\psi }_{\mathrm{burst}}={\nu }_{\mathrm{SF},\mathrm{burst}}{M}_{\mathrm{cold},\mathrm{burst}},$$

$${\tau }_{\star ,\mathrm{burst}}=\max [f_{\mathrm{dyn}}{\tau }_{\mathrm{dyn},\mathrm{bulge}},{\tau }_{\mathrm{burst},\min }].$$

$${\dot{M}}_{\mathrm{eject}}=\beta (V_{\mathrm{c}})\psi ={\left(V_{\mathrm{c}}/V_{\mathrm{SN}}\right)}^{-{\gamma }_{\mathrm{SN}}}\psi .$$

$${\dot{M}}_{\mathrm{return}}={\alpha }_{\mathrm{ret}}\frac{M_{\mathrm{res}}}{{\tau }_{\mathrm{dyn},\mathrm{halo}}},$$

$$\beta (V_{\mathrm{c}},z)=\{\begin{array}{ll}[V_{\mathrm{c}}/V_{\mathrm{SN}}^{\prime }(z){]}^{-{\gamma }_{\mathrm{SN}}^{\prime }}& V_{\mathrm{c}}\le V_{\mathrm{thresh}}\\ [V_{\mathrm{c}}/V_{\mathrm{SN}}(z){]}^{-{\gamma }_{\mathrm{SN}}}& V_{\mathrm{c}}>V_{\mathrm{thresh}},\end{array}$$

$$V_{\mathrm{SN}}(z)=\{\begin{array}{ll}{V}_{SN2}& z>z_{SN2}\\ c_{0}z+c_1& z_{SN2}\le z\le z_{SN1}\\ V_{SN1}& z<z_{SN1},\end{array}$$

$$\mathcal{R}(z)=\frac{{\int }_z^{\infty }ϵ(z^{\prime })\mathrm{d}{z}^{\prime }}{n_{\mathrm{H}}},$$

$$\mathcal{R}(z_{\mathrm{re},\mathrm{full}})=\frac{1+N_{\mathrm{rec}}}{f_{\mathrm{esc}}},$$

$$\frac{{\mathrm{d}}^3\eta }{\mathrm{d}\ln S_{\nu }\mathrm{d}z\mathrm{d}\Omega }=\frac{\mathrm{d}n}{\mathrm{d}\ln L_{\nu }}\frac{{\mathrm{d}}^{2}V}{\mathrm{d}z\mathrm{d}\Omega },$$

Full text

††pubyear: 2017††pagerange: Predictions for deep galaxy surveys with JWST from $\Lambda$CDM–A.4

1 Introduction

The James Webb Space Telescope (JWST) is scheduled for launch in October 2018 and is expected to significantly advance our understanding of the high-redshift ($z\gtrsim 7$) Universe (e.g. Gardner et al., 2006). Two of its on-board instruments, the Near InfraRed Camera (NIRCam) and the Mid InfraRed Instrument (MIRI), are dedicated to obtaining broadband photometry over the wavelength range $0.7-25.5$ $\mu$m with unprecedented sensitivity and angular resolution. This wavelength coverage will enable JWST to probe the rest-frame UV/optical/near-IR spectral energy distributions (SEDs) of high-redshift ($z\gtrsim 7$) galaxies, opening up a hitherto unexplored regime of galaxy formation and evolution.

An early breakthrough in the study of galaxies in the high-redshift Universe came from the identification of galaxies at $z\sim 3$ using the Lyman-break technique (e.g. Steidel & Hamilton, 1993; Steidel et al., 1996). This study took advantage of the break in galaxy SEDs produced at the Lyman limit ($912$ Å) to identify galaxies at $z\sim 3$ by searching for ‘dropouts’ in a set of broadband photometric filters. The significance of this development in the context of galaxy formation and evolution, in particular, the implications for the cosmic star formation rate density and the formation of massive galaxies in the $\Lambda$CDM cosmological model, was discussed in Baugh et al. (1998, see also Mo & Fukugita 1996 and Mo et al. 1999). A further advance came with the installation of the Advanced Camera for Surveys (ACS) on the Hubble Space Telescope which, using the $z$-band, pushed the Lyman-break technique selection to $z\sim 6$ (e.g. Bouwens et al., 2003; Stanway et al., 2003). At these redshifts the Lyman-break technique makes use of the fact that neutral hydrogen in the intergalactic medium (IGM) effectively absorbs radiation with wavelengths shorter than the Lyman $\alpha$ transition ($1216$ Å), resulting in a strong break in the galaxy SED at the observer-frame wavelength of this transition. Installation of the Wide-Field Camera 3 (WFC3) with near-IR filters increased the number of galaxies that could be identified at $z\sim 7$ (e.g. Bouwens et al. 2010; Wilkins et al. 2010), pushing the samples of galaxies at these redshifts into the thousands, with a few examples at $z\sim 10$. These advances have been complemented by ground-based telescopes, such as the Visible and Infrared Survey Telescope for Astronomy (VISTA), that typically provide a larger field of view than their space-based counterparts; this has allowed the bright end of the rest-frame far-UV luminosity function to be probed robustly at $z\sim 7$ (e.g. Bowler et al., 2014).

As observations in the near-IR with Hubble have identified the highest-redshift galaxies to date, a wealth of further information regarding galaxy properties at intermediate redshifts ($z\sim 3$) has come from surveys with the Spitzer Space Telescope in the same wavelength range that will be probed by JWST (e.g. Labbé et al., 2005; Caputi et al., 2011; Caputi et al., 2015), though JWST will have greater angular resolution and sensitivity than Spitzer. As a result, JWST is expected to greatly increase the number of observed galaxies at $z\gtrsim 7$, providing important information about their SEDs which can help characterise their physical properties, whilst also extending observations of the high-redshift Universe towards the first luminous objects at the end of the so-called cosmic dark ages.

Understandably, in recent years a number of studies have made predictions for galaxy formation at the high redshifts expected to be probed by JWST. Numerical hydrodynamical simulations such as the First Billion Years simulation suites (e.g. Paardekooper et al., 2013), the BlueTides simulation (e.g. Wilkins et al., 2016; Wilkins et al., 2017), the Renaissance simulations suite (e.g. Xu et al., 2016) and others (e.g. Dayal et al., 2013; Shimizu et al., 2014) have typically focused on the earliest stars and galaxies as potential sources of reionization and have made predictions for the galaxy rest-frame UV luminosity function. These calculations are generally only run to very high redshift ($z\gtrsim 6$), as the computational expense of adequately resolving the physical processes involved becomes prohibitive towards later times. As such, there is considerable uncertainty as to whether such simulations would be able to reproduce the galaxy population at $z=0$. It should be noted, however, that some cosmological hydrodynamical simulations are able to reproduce the galaxy population at $z=0$ (e.g. Vogelsberger et al., 2014; Schaye et al., 2015).

Simple empirical models (e.g. Behroozi & Silk, 2015; Mason et al., 2015; Mashian et al., 2016) have also been used to make predictions for the high-redshift galaxy rest-frame UV luminosity function. These models are much less computationally expensive than the hydrodynamical schemes mentioned above and as such can be run to $z=0$. However, they ignore most of the physical processes of galaxy formation and instead rely on arbitrary scalings to compute a small number of galaxy properties from those of the host halo. As such they have a limited predictive power and a physical interpretation of their predictions is foregone. Nevertheless, these models can reproduce evolution of the rest-frame far-UV luminosity function in reasonable agreement with observations for $z\lesssim 8$ (though they are often calibrated on these data at some redshifts), and suggest small numbers of galaxies will be observable with future JWST galaxy surveys at $z\gtrsim 10$.

A powerful technique for studying the formation and evolution of galaxies is semi-analytical modelling (see the reviews by Baugh 2006 and Benson 2010). In such models, the complex physical processes of galaxy formation are fully accounted for and are described by simplified prescriptions that are based on either theoretical arguments or observational or simulation data. This makes the galaxy formation and evolution calculation more computationally tractable, whilst still encapsulating its intrinsic complexity. The free parameters introduced as a result of these simplified prescriptions are then calibrated against a predetermined set of observational data, often requiring that any viable model reproduce the galaxy population observed at $z=0$. Once this has been done the model is fully specified and can be used to make genuine predictions for a wide range of other observable properties at any redshift. An advantage of semi-analytical models is that their predictions can then be readily interpreted in terms of the modelling and interplay of the physical processes involved, and comparing their predictions to observational data provides a test of our understanding of these processes.

Clay et al. (2015) made predictions for the evolution of the rest-frame far-UV luminosity function for $z\sim 4-7$ using the semi-analytical model of Henriques et al. (2015). However, this model underpredicted the bright end of the observed luminosity function and relied on an ad hoc scaling with redshift of the dust optical depth in galactic discs. Liu et al. (2016) achieved a better fit to the observed rest-frame far-UV luminosity function using the semi-analytical model meraxes (Mutch et al., 2016). However, this model only provides predictions for $z\gtrsim 5$ and does not account for feedback from an active galactic nucleus (AGN) [though see Qin et al. (2017) for an updated version of this model that addresses these shortcomings]. Additionally, neither of these works attempt to model dust emission and thus their predictions are restricted to the rest-frame UV/optical/near-IR.

Here we present theoretical predictions for deep galaxy surveys with JWST NIRCam and MIRI, in the form of luminosity functions, number counts and redshift distributions from a semi-analytical model of hierarchical galaxy formation within $\Lambda$CDM (Lacey et al., 2016). The model provides a physically-motivated computation of galaxy formation and evolution from $z\gtrsim 20$ to $z=0$. For computing galaxy SEDs the model is coupled with the spectrophotometric code grasil (Silva et al., 1998), which takes into account the absorption and re-emission of stellar radiation by interstellar dust by solving the equations of radiative transfer in an assumed geometry. This broadens the predictive capability of the model to the full wavelength range that will be probed by JWST. The Lacey et al. model is calibrated to reproduce a broad range of observational data at $z\lesssim 6$, these include the optical and near-IR luminosity functions at $z=0$, the evolution of the rest-frame near-IR luminosity functions for $z=0-3$, far-IR/sub-mm galaxy number counts and redshift distributions, and the evolution of the rest-frame far-UV luminosity function for $z=3-6$. The predictions of this model presented in this work thus represent an exciting opportunity to test the modelling and interplay of the physical processes of galaxy formation against JWST observations at higher redshifts than those at which the model was calibrated. At the same time, they can potentially inform future JWST galaxy survey strategies.

A shortcoming of the fiducial Lacey et al. model, however, is that it does not reproduce the reionization redshift of $z=8.8_{-1.4}^{+1.7}$ inferred from cosmic microwave background (CMB) data by Planck Collaboration et al. (2016). This is an important constraint for high-redshift predictions of the galaxy population. The model produces too few ionizing photons at early times, reionizing the Universe at $z=6.3$ (Hou et al., 2016).

A simple and effective solution to this shortcoming was proposed by Hou et al. (2016) who, motivated by the dynamical supernova (SN) feedback model of Lagos et al. (2013), allowed the strength of SN feedback in the Lacey et al. (2016) model to vary as a function of redshift. Reducing the strength of SN feedback at high redshift meant that the model could produce more ionizing photons at this epoch. The evolving feedback also enabled this model to reproduce the $z=0$ luminosity function of the Milky Way satellites, as well as their metallicity–stellar mass relation. These further successes in matching observational data do not come at the expense of the agreement of the model with the data against which it was originally calibrated at $z\lesssim 6$, but it does introduce new parameters to describe the effects of SN feedback.

SN feedback is an extremely important physical process in galaxy evolution (e.g. Larson, 1974; White & Rees, 1978; White & Frenk, 1991; Cole, 1991). However, its precise details, for example, exactly how energy input from supernovae (SNe) should couple to the interstellar medium (ISM), are still poorly understood. This is mainly due to the difficulty of fully resolving individual star-forming regions in hydrodynamical simulations spanning a cosmologically significant time period and volume (e.g. Vogelsberger et al., 2014; Schaye et al., 2015). It is hoped that comparing the predictions of phenomenological models of SN feedback, such as those presented here, with future observations from JWST, will lead to a greater understanding of the efficiency of this crucial process.

This paper is structured as follows: in Section 2 we present some of the pertinent details of our galaxy formation model and the evolving feedback variant, the radiative transfer code used for the computation of UV-to-mm galaxy SEDs and some information regarding the coupling of these two models. In Section 3 we present our main results111Some of the model data presented here will be made available at http://icc.dur.ac.uk/data/. For other requests please contact the first author.; these include galaxy luminosity functions, number counts and redshift distributions for varying exposures, and angular sizes in each of the NIRCam and MIRI broadband filters. We also present predictions for the evolution of some of the physical properties of the model galaxies (e.g. stellar masses, star formation rates) and compare some model predictions to available high-redshift ($z\gtrsim 7$) observational data. We conclude in Section 4. A brief discussion of the dependence of our high-redshift predictions on some of our model assumptions is given in Appendix A.

Throughout we assume a flat $\Lambda$CDM cosmology with cosmological parameters consistent with recent Planck satellite results (Planck Collaboration et al., 2016)222$\Omega_{\rm m}=0.307$, $\Omega_{\Lambda}=0.693$, $h=0.678$, $\Omega_{\rm b}=0.0483$, $\sigma_{8}=0.829$. All magnitudes are presented in the absolute bolometric (AB) system (Oke, 1974).

2 The Theoretical Model

In this Section we introduce our galaxy formation model, which combines a dark matter only $N$-body simulation, a semi-analytical model of galaxy formation (galform) and the spectrophotometric radiative transfer code grasil (Silva et al., 1998) for computing UV-to-mm galaxy SEDs.

2.1 GALFORM

The Durham semi-analytic model of hierarchical galaxy formation, galform, was introduced in Cole et al. (2000), building on ideas outlined earlier by White & Rees (1978), White & Frenk (1991) and Cole et al. (1994). Galaxy formation is modelled ab initio, beginning with a specified cosmology and a linear power spectrum of density fluctuations, and ending with predicted galaxy properties at different redshifts.

Galaxies are assumed to form from baryonic condensation within the potential wells of dark matter halos, with their subsequent evolution being controlled in part by the merging history of the halo. Here, these halo merger trees are extracted directly from a dark matter only $N$-body simulation (e.g. Helly et al., 2003; Jiang et al., 2014) as this approach allows us to predict directly the spatial distribution of the galaxies. We use a new ($800$ Mpc)3 Millennium-style simulation (Springel et al., 2005) with cosmological parameters consistent with recent Planck satellite results (Planck Collaboration et al., 2016), henceforth referred to as P–Millennium (Baugh et al. in preparation; McCullagh et al. 2017). This large volume (800 Mpc)3 gives the bright end of our predicted luminosity functions a greater statistical precision than could be achieved using dark matter only simulations with a better halo mass resolution but smaller volume.

The halo mass resolution of this simulation is $2.12\times 10^{9}$ $h^{-1}$ M*☉, where a halo is required to have at least $20$ dark matter particles and is defined according to the ‘DHalo’ algorithm (Jiang et al., 2014). This mass resolution is approximately an order of magnitude better than previous dark matter simulations that were used with this galaxy formation model. For example, the MR7 simulation (Springel et al., 2005; Guo et al., 2013) in which the Lacey et al. (2016) model was originally implemented had a halo mass resolution of $1.87\times 10^{10}$ $h^{-1}$ M☉*. This improved resolution is particularly important for predictions of the high-redshift Universe where, due to the hierarchical nature of structure formation in $\Lambda$CDM, galaxy formation takes place in lower mass halos. This halo mass resolution is in the regime where ignoring baryonic effects on the dark matter, an implicit assumption of the semi-analytical technique, is still a reasonable one. The ‘back-reaction’ due to baryonic effects, such as feedback processes, on the dark matter is expected to reduce the mass of dark matter halos by only $\sim 30$ per cent at the mass limit of the P–Millennium simulation (e.g. Sawala et al., 2013).

We have tested that the results presented in this paper have converged with respect to the halo mass resolution used in the P-Millennium simulation and that any artificial features this introduces into our predicted luminosity functions are at luminosities fainter than those studied here. For example, at $z=10$ we find a halo mass resolution ‘turn-over’ in our predicted rest-frame far-UV ($1500$ Å) luminosity function at $M_{\rm AB}-5\log_{10}h\sim-14$, which is approximately one magnitude fainter than the sensitivity of a $10^{6}$ s exposure with the NIRCam–F150W filter at this redshift.

Baryonic physics in galform are included as a set of coupled differential equations which track the exchange of mass and metals between between the stellar, cold disc gas and hot halo gas components in a given halo. These equations comprise simplified prescriptions for the physical processes (e.g. gas cooling, star formation and feedback) understood to be important for galaxy formation. We discuss some of the main features of the model below and refer the interested reader to Lacey et al. (2016) for more details.

2.1.1 The star formation law and stellar initial mass function

Star formation in the galactic disc is based on the surface density of molecular gas. Cold disc gas is partitioned into molecular and atomic components based on an empirical relation involving the mid-plane gas pressure, $P$, proposed by Blitz & Rosolowsky (2006) based on observations of nearby galaxies:

[TABLE]

where $R_{\rm mol}$ is the ratio of molecular to atomic gas; $\alpha_{\rm P}=0.8$ and $P_{0}=1.7\times 10^{4}$ cm*-3* K based on the local observations of Leroy et al. (2008). It is assumed that gas and stars are distributed in a exponential disc, the radial scalelength of which is predicted by galform (see Section 2.1.4). The star formation rate surface density is then given by

[TABLE]

where $f_{\rm mol}=R_{\rm mol}/(1+R_{\rm mol})$ and the parameter $\nu_{\rm SF}=0.74$ Gyr*-1*, based on the observations of Bigiel et al. (2011). This expression is then integrated over the whole disc to yield the global star formation rate, $\psi$. For further details of this star formation law we refer the reader to Lagos et al. (2011). For star formation in the galactic disc a Kennicutt (1983) stellar initial mass function (IMF) is assumed. This IMF is described by $x=0.4$ in $\mathrm{d}N/\mathrm{d}\ln m\propto m^{-x}$ for $m<1$ M*⊙* and $x=1.5$ for $m>1$ M*⊙* [for reference, a Salpeter (1955) IMF has an unbroken slope of $x=1.35$].

Star formation in bursts, triggered by a dynamical process (see Section 2.1.3), takes place in a forming galactic bulge. It is assumed that $f_{\rm mol}\approx 1$ and the star formation rate depends on the dynamical timescale of the bulge

[TABLE]

where $\nu_{\rm SF,burst}=1/\tau_{\star\rm,burst}$ and

[TABLE]

Here $\tau_{\rm dyn,bulge}$ is the dynamical time of the bulge and $f_{\rm dyn}$ and $\tau_{\rm burst,min}$ are model parameters. This means that for large dynamical times the star formation rate scales with the dynamical time, but has a floor value when the dynamical time of the bulge is short. Here $f_{\rm dyn}=20$ and $\tau_{\rm burst,min}=100$ Myr (Lacey et al., 2016).

For star formation in bursts, it is assumed that stars form with a top-heavy stellar initial mass function (IMF), described by a slope of $x=1$ in $\mathrm{d}N/\mathrm{d}\ln m\propto m^{-x}$. This assumption is primarily motivated by the requirement that the model reproduce the observed far-IR/sub-mm galaxy number counts and redshift distributions (e.g. Baugh et al., 2005; Cowley et al., 2015; Lacey et al., 2016, see also Fontanot 2014 for a study of the effects of IMF variation in semi-analytical models). It should be noted that the slope in this new model is much less top-heavy than the one suggested by Baugh et al. (2005), where $x=0$ was assumed.

The assumption of a top-heavy IMF for starburst galaxies is often seen as controversial. For example, in their review of observational studies Bastian et al. (2010) argue against significant IMF variation in the local Universe. However, Gunawardhana et al. (2011) infer an IMF for nearby star-forming galaxies that becomes more top-heavy with increasing star formation rate, reaching a slope of $x\approx 0.9$, and a similar IMF slope was inferred for a star-forming galaxy at $z\sim 2.5$ by Finkelstein et al. (2011). Both of these studies utilise modelling of a combination of nebular emission and broadband photometry to infer an IMF slope. More recently, Romano et al. (2017) inferred an IMF slope of $x=0.95$ in nearby starburst galaxies through modelling the observed CNO isotopic ratios. Thus whilst the issue of a varying IMF is far from resolved, there are a number of observational studies that support both this assumption and the adopted value of $x=1$.

2.1.2 Feedback processes

The model includes three modes of feedback from stars and AGN on the galaxy formation process.

Photoionization feedback: the IGM is reionized and photo-heated by ionizing photons produced by stars. This inhibits star formation through (i) preventing gas accretion onto low-mass haloes through an increased IGM pressure and (ii) continued photo-heating reducing the cooling rate of gas already within haloes. Here a simple scheme is implemented that assumes that after the IGM is reionized at a fixed redshift, $z_{\rm reion}$, no cooling of gas occurs in haloes with circular velocities $V_{\rm vir}<V_{\rm crit}$. Here we assume $z_{\rm reion}=10$ (Dunkley et al., 2009)333This value of $z_{\rm reion}$ is slightly different to the one predicted by the models, however varying this parameter within the range suggested by the model predictions has a negligible effect on our results [see Fig. 13(c)]. and $V_{\rm crit}=30$ km s*-1*, based on hydrodynamical simulations (e.g. Hoeft et al., 2006; Okamoto et al., 2008). Whilst this model is very simple it is based on a self-consistent calculation of reionization in galform described by Benson et al. (2002), and it was shown by Font et al. (2011) to reproduce results from more detailed treatments (e.g. Muñoz et al., 2009; Busha et al., 2010) of this process.

Supernova (SN) feedback: the injection of energy into the ISM from SNe ejects gas from the disc to beyond the virial radius of the halo at a rate, $\dot{M}_{\rm eject}$. As SNe are short-lived this rate is proportional to the star formation rate, $\psi$, according to a ‘mass loading’ factor, $\beta$, such that

[TABLE]

Here $V_{\rm c}$ is the circular velocity of the disc; $\psi$ is the star formation rate; and $V_{\rm SN}$ and $\gamma_{\rm SN}$ are adjustable parameters. We assume $V_{\rm SN}=320$ km s*-1* (Lacey et al., 2016) and $\gamma_{\rm SN}=3.4$ (Baugh et al. in preparation, see Section 2.1.5). The ejected gas accumulates in a reservoir of mass $M_{\rm res}$, and then falls back within the virial radius at a rate

[TABLE]

where $\tau_{\rm dyn,halo}$ is the halo dynamical time and $\alpha_{\rm ret}=1.0$ (Baugh et al., in preparation, see also Section 2.1.5).

AGN feedback: the model implements a hot-halo mode AGN feedback, first implemented into galform by Bower et al. (2006). Energy released by the direct accretion of hot gas from the halo onto the supermassive black hole (SMBH) powers relativistic jets that deposit thermal energy into the hot halo gas and thus inhibit further cooling. In the model, gas cooling is turned off if: (i) the gas is cooling quasi-statically (i.e. the cooling time is long compared to the free-fall time), and (ii) the SMBH is massive enough such that the power required to balance the radiative cooling luminosity of the gas is below some fraction of its Eddington luminosity.

2.1.3 Dynamical processes

Morphological transitions occur, and starbursts are triggered, through dynamical processes. These are either galaxy mergers, where the orbit of a satellite galaxy in a dark matter halo has decayed through dynamical friction such that it merges with the central galaxy, or disc instabilities, in which the galactic disc becomes sufficiently self-gravitating that it is unstable to bar formation [using the criterion of Efstathiou et al. (1982), which is based on simulations of isolated disc galaxies].

Major galaxy mergers (and minor mergers above a baryonic mass ratio) and all disc instabilities trigger bursts of star formation. In these, all of the cold gas in the disc is transferred to a forming bulge/spheroid and forms stars according to the star formation law for bursts and assuming a top-heavy IMF as is described earlier.

When we refer to starburst galaxies throughout this paper, we are referring to this dynamically triggered star formation rather than, for example, a galaxy’s position on the specific star formation rate–stellar mass plane. This distinction is discussed in more detail in Cowley et al. (2017).

2.1.4 Galaxy sizes

In galform it is assumed that a disc with an exponential radial profile is formed from cold gas once it has had sufficient time to cool and fall to the centre of the dark matter halo potential well. The size of the disc is calculated by assuming conservation of angular momentum and centrifugal equilibrium (Cole et al., 2000).

Galaxy bulges/spheroids are assumed to have a projected $r^{1/4}$ density profile and are formed through a dynamical process, either a disc instability or a galaxy merger. The size of the bulge is determined by the conservation of energy for the components involved i.e. baryons and dark matter in the disc and bulge of the galaxies (Cole et al., 2000).

2.1.5 Changes to the Lacey et al. (2016) model

This work assumes different cosmological parameters from those assumed by Lacey et al. (2016), and utilises an $N$-body simulation with a better halo mass resolution. The model used here also incorporates an improved prescription for the merger timescale of satellite galaxies (Simha & Cole, 2017), which was first introduced into galform by Campbell et al. (2015), but was not considered by Lacey et al. This new treatment accounts for the effects of both dynamical friction and tidal stripping on the sub-halo and thus more closely follows the underlying $N$-body simulation than the analytical prescription used in galform previously (Lacey & Cole, 1993; Jiang et al., 2008). Additionally, the earlier prescription for the merger timescale resulted in a radial distribution of satellite galaxies that was too centrally concentrated (Contreras et al., 2013).

As a result of these changes, it is necessary to adjust some of the galaxy formation parameters in the fiducial model so that it can still reproduce certain pre-specified observational datasets to the desired accuracy. The adjustments will be discussed in more detail in Baugh et al. (in preparation); however, we briefly summarise the main ideas here. A minor reduction in the number of bright galaxies at $z=0$ required the gas reincorporation timescale factor, $\alpha_{\rm ret}$ (equation 6), to be increased from $0.64$ to $1.00$, thus returning gas ejected by SN feedback to the hot halo faster. Additionally, the change in the halo mass resolution resulted in the number of faint galaxies being slightly overpredicted, so it was necessary to increase the strength of the SN feedback through increasing the value of the parameter $\gamma_{\rm SN}$ (equation 5) from $3.2$ to $3.4$ to mitigate this.

We summarise these minor adjustments to the model presented in Lacey et al. (2016) in Table 1.

2.2 Evolving supernova feedback and the redshift of reionization

As mentioned earlier, a shortcoming of the fiducial Lacey et al. (2016) model is that it does not reionize the Universe at a redshift as high as implied by recent Planck data, as it does not produce enough ionizing photons at early enough times. Here we discuss the variant feedback model of Hou et al. (2016) which provides a simple and effective solution to this shortcoming.

In the fiducial galform model, gas outflows due to SN feedback are implemented according to equation 5. A dynamical model of SN feedback, which followed the evolution of pressurized SNe bubbles in a multiphase ISM was implemented into the galform framework by Lagos et al. (2013). Whilst this SN feedback model is not complete as it only considers gas escaping from the galactic disc, and not from the halo, it suggested that the dependence of the mass loading, $\beta$, solely on galaxy circular velocity may be an oversimplification of this physical process. This standard parametrization of $\beta$ is motivated by the fact that gas outflows should depend on the depth of the gravitational potential well, for which $V_{\rm c}$ is a commonly used proxy. However, it is reasonable to propose that it may also depend on properties such as the gas density, the gas metallicity and the molecular gas fraction. For example, the local gas density and metallicity determine the local gas cooling rate in the ISM and in turn the fraction of the injected SN energy that can be used to launch outflows; and dense molecular gas may not be affected by SNe explosions and thus not ejected in such outflows. These additional physical parameters will evolve with redshift, and may not be well described by a power law with $V_{\rm c}$.

In order to produce more ionizing photons, and thus reionize the Universe earlier than the fiducial model, Hou et al., motivated by the dynamical SN feedback model of Lagos et al. (2013), introduced a break into the power-law parametrisation of the mass loading factor and also a redshift dependence into its normalisation, such that

[TABLE]

where $V_{\rm thresh}$ and $\gamma_{\rm SN}^{\prime}$ are additional adjustable parameters [$V_{\rm SN}^{\prime}(z)$ is set by the condition that $\beta$ be a continuous function at $V_{\rm c}=V_{\rm thresh}$]. The redshift evolution of the normalisation is parametrised as

[TABLE]

where $V_{\rm SN2}$, $z_{\rm SN2}$ and $z_{\rm SN1}$ are additional adjustable parameters [the constants $c_{0}$ and $c_{1}$ are set by the condition that $V_{\rm SN}(z)$ be a continuous function]. This form parametrizes our ignorance of the precise physical mechanisms at play, whilst allowing for the dependencies of the mass outflow rates on physical properties other than $V_{\rm c}$, as discussed above, to be described. Though we acknowledge that a detailed physical interpretation of this variant feedback model is somewhat lacking, it provides a tractable approximation that is calibrated not only on the reionization redshift suggested by Planck data, but also on the luminosity function and metallicity–stellar mass relation of $z=0$ Milky Way satellites. These independent observational data provide strong constraints on the form that the mass-loading factor for SN feedback can take, as is discussed in Hou et al. (2016).

Here we use the same values for the additional adjustable parameters in this variant feedback model as Hou et al.: $V_{\rm thresh}=50$ km s*-1*, $\gamma_{\rm SN}^{\prime}=1.0$, $V_{\rm SN2}=180$ km s*-1*, $z_{\rm SN1}=4$ and $z_{\rm SN2}=8$, without any further calibration, although we remind the reader that the value for $\gamma_{\rm SN}$ is different to the one used by Hou et al. Additionally, we adopt $V_{\rm SN1}=V_{\rm SN}$, as was done by Hou et al.

We show the predicted redshift of reionization for both the fiducial model (lc16) and the evolving feedback variant (lc16.EvolFB) in Fig. 1. Following Hou et al. we calculate the ratio, $\mathcal{R}(z)$, of ionizing photons produced before redshift $z$, to the number density of hydrogen nuclei as

[TABLE]

where $\epsilon(z^{\prime})$ is the number of hydrogen-ionizing photons produced per unit comoving volume per unit redshift at redshift $z^{\prime}$, and $n_{\rm H}$ is the comoving number density of hydrogen nuclei. The Universe is assumed to be fully ionized at redshift $z_{\rm re,\,full}$, for which,

[TABLE]

where $N_{\rm rec}$ is the mean number of recombinations per hydrogen atom up to reionization, and $f_{\rm esc}$ is the fraction of ionizing photons that can escape into the IGM from the galaxy producing them. Here we adopt $N_{\rm rec}=0.25$ and $f_{\rm esc}=0.2$ as was done by Hou et al. This gives a threshold for reionization of $\mathcal{R}(z_{\rm re,full})=6.25$.

Observations of the CMB (e.g. Planck Collaboration et al., 2016) directly constrain the electron scattering optical depth to recombination, which is then converted to a reionization redshift by assuming a simple model for the redshift dependence of reionization (e.g. Appendix B of Lewis et al., 2008). The redshift of reionization is commonly expressed in terms of the redshift, $z_{\rm re,\,half}$, at which half of the IGM is reionized. Here we assume $\mathcal{R}_{\rm re,\,half}=0.5\,\mathcal{R}_{\rm re,\,full}$ as was done by Hou et al. The value of $\mathcal{R}_{\rm re,\,half}$ is shown as the horizontal dot-dashed line in Fig. 1. We can see that the evolving feedback model predicts $z_{\rm re,\,half}=8.9$, in good agreement with the $68$ per cent confidence interval inferred from Planck satellite data (Planck Collaboration et al., 2016), $z_{\rm re,\,half}=8.8_{-1.4}^{+1.7}$. For the fiducial model the reionization redshift turns out to be lower, $z_{\rm re,\,half}=6.9$, which is discrepant by $\sim 1.5\sigma$ with the Planck data.

2.3 The Dust Model

We use the spectrophotometric radiative transfer code grasil (Silva et al., 1998) to compute model galaxy SEDs. Using the star formation and metal enrichment histories, gas masses and geometrical parameters predicted by galform, and assuming a composition and geometry for interstellar dust, grasil computes the SEDs of the model galaxies, accounting for dust extinction (absorption and scattering) of radiation and its subsequent re-emission. In this Section, we briefly describe the grasil model. For further details we refer the reader to Silva et al. (1998) and Granato et al. (2000).

Here grasil assumes that stars exist in a disc + bulge system, as is the case in galform. The disc has a radial and vertical exponential profile with scale lengths, $h_{R}$ and $h_{z}$, and the bulge is described by an analytic King model profile, $\rho\propto(r^{2}+r_{\rm c}^{2})^{-3/2}$ out to a truncation radius, $r_{\rm t}$. The half-mass radii, $r_{\rm disc}$ and $r_{\rm bulge}$, are predicted by galform. By definition, given the assumed profiles, the bulge core radius is related to the half-mass radius by $r_{\rm c}=r_{\rm bulge}/14.6$ whilst the radial disc scale-length, $h_{\rm R}$, is related to the half-mass disc radius by $h_{R}=r_{\rm disc}/1.68$. Star formation histories are calculated separately for the disc and bulge by galform. For galaxies undergoing a starburst, the burst star formation, as well as the associated gas and dust, are assumed to also be in an exponential disc but with a half-mass radius, $r_{\rm burst}=\eta r_{\rm bulge}$, rather than $r_{\rm disc}$, where $\eta$ is an adjustable parameter. The disc axial ratio, $h_{z}/h_{R}$, is a parameter of the grasil model; for starburst galaxies, the axial ratio of the burst is allowed to be different from that of discs in quiescent galaxies.

The gas and dust exist in an exponential disc, with the same radial scale-length as the disc stars but in general with a different scale-height, so $h_{z}\mathrm{(dust)}/h_{z}\mathrm{(stars)}$ is an adjustable parameter. The gas and dust are assumed to exist in two components: (i) giant molecular clouds in which stars form, escaping on some time scale, $t_{\rm esc}$, and (ii) a diffuse cirrus ISM. The total gas mass, $M_{\rm cold}$, and metallicity, $Z_{\rm cold}$, are calculated by galform. The fraction of gas in molecular clouds is determined by the parameter $f_{\rm cloud}$. The cloud mass, $m_{\rm cloud}$, and radius, $r_{\rm cloud}$, are also parameters, though the results of the model depend only on the ratio, $m_{\rm cloud}/r_{\rm cloud}^{2}$, which determines (together with the gas metallicity) the optical depth of the clouds.

The dust is assumed to consist of a mixture of graphite and silicate grains and polycyclic aromatic hydrocarbons (PAHs), each with a distribution of grain sizes. The grain mix and size distribution were determined by Silva et al. so that the extinction and emissivity properties of the local ISM are reproduced using the optical properties of the dust grains tabulated by Draine & Lee (1984). At long wavelengths ($\lambda>30$ $\mu$m) this results in a dust opacity that approximates $\kappa_{\rm d}\propto\lambda^{-2}$. However, in galaxies undergoing a starburst this is modified (for $\lambda>100$ $\mu$m) such that $\kappa_{\rm d}\propto\lambda^{-\beta_{\rm b}}$, where $\beta_{\rm b}$ is treated as an adjustable parameter. Laboratory measurements suggest that values in the range $\beta_{\rm b}=1.5-2$ are acceptable (Agladze et al. 1996). Here a value of $\beta_{\rm b}=1.5$ is adopted (Lacey et al., 2016). The total dust mass in a galaxy is proportional to the cold gas mass and metallicity, both of which are predicted by galform.

The adopted values of adjustable grasil parameters are summarised in Table 2.

For the parameters which are analogous to those in the dust model used by Lacey et al. (2016): $f_{\rm cloud}$, $m_{\rm cloud}/r_{\rm cloud}^{2}$, $t_{\rm esc}$ and $\beta_{\rm b}$, we use the values chosen by Lacey et al. For other parameters specific to the grasil model, we use the values chosen by by Baugh et al. (2005, see also Lacey et al. 2008, Swinbank et al. 2008 and Lacey et al. 2011), which was the last time a published version of galform was coupled with grasil in the manner presented here.

The luminosities of the stellar components are calculated assuming the Maraston (2005) evolutionary population synthesis model, as is done in Lacey et al. (2016). grasil then calculates the radiative transfer of the stellar radiation through the interstellar dust. For molecular clouds, a full radiative transfer calculation is performed. For the diffuse cirrus the effects of scattering are included approximately by using an effective optical depth for the absorption $\tau_{\rm abs,eff}=[\tau_{\rm abs}(\tau_{\rm abs}+\tau_{\rm scat})]^{1/2}$. The dust-attenuated stellar radiation field can be calculated at any point inside or outside the galaxy. grasil then computes the final galaxy SED by calculating the absorption of stellar radiation, thermal balance and the re-emission of radiation for each grain species and size at every point in the galaxy.

Examples of predicted star formation histories and the resulting galaxy UV-to-mm SEDs computed by grasil are shown in Fig. 2. One can see that the star formation histories are extremely ‘bursty’ at early times when the Universe is a few Gyr old. Significant dust extinction and re-emission is evident for each of the galaxy SEDs shown. There are also a number of interesting features in the galaxy SEDs. These include: (i) Lyman-continuum breaks in the galaxy SEDs at $912$ Å; (ii) a prominent $4000$ Å break for the $z=0$ galaxy, indicative of an old stellar population (which would be expected from the smoothly declining star formation history of this galaxy); (iii) dust emission approximating a modified blackbody that peaks at $\lambda_{\rm rest}\approx 100$ $\mu$m, indicative of cold ($\sim 30$ K) dust, though the peak of the emission shifts to shorter wavelengths with increasing redshift suggesting hotter dust, and (iv) PAH emission lines in the cirrus dust at $\lambda_{\rm rest}=3.3$, $6.2$, $7.7$, $8.6$, and $11.3$ $\mu$m.

Once an SED has been computed, luminosities in specified bands are calculated by convolving the SED (redshifted into the observer frame) with the filter transmission of interest. We use the Meiksin (2005) prescription for attenuation of radiation in the IGM due to neutral hydrogen, also shown in Fig. 2.

2.4 Coupling galform and grasil

Here we briefly describe how the galform and grasil models are used in conjunction. For further details, we refer the reader to Granato et al. (2000).

Due to the computational expense of running grasil ($\sim 3-5$ CPU mins per galaxy) it is not feasible to compute an SED for each galaxy in the simulation volume, as has been discussed in previous studies (e.g. Granato et al., 2000; Almeida et al., 2010; Lacey et al., 2011). However, for the purposes of constructing luminosity functions, it is possible to circumvent this by running grasil on a sample of galaxies, from which the luminosity function can be constructed if the galaxies in question are weighted appropriately. We choose to sample galaxies according to their stellar mass such that $\sim 10^{3}$ galaxies per dex of stellar mass are sampled. We use a lower mass limit of $10^{6}$ $h^{-1}$ M*☉*, which we choose so that any artificial features it introduces into our predicted luminosity functions (see Section 3.2) are at fainter luminosities than those investigated here. This represents a factor of $\sim 10$ increase over the number of galaxies sampled by Granato et al. (2000).

The procedure that we use to construct luminosity functions in a given band at each output redshift is as follows: (i) run galform to the redshift of interest; (ii) create a subsample of galaxies; (iii) re-run galform to output the star formation and metal enrichment history for each of the sampled galaxies; (iv) run grasil on each of the sampled galaxies to produce a predicted SED; (v) convolve the output SED with the relevant broadband filter response and IGM attenuation curve (Meiksin, 2005) and (vi) construct the galaxy luminosity function using the weights from the initial sampling and luminosities from the previous step.

We have made a number of improvements to steps (iii) to (v) above, which allow us to run grasil for samples of $\sim 10^{5}$ galaxies for each model, spread over $25$ output redshifts from $z=16$ to $z=0$. For each model, this takes $\sim 7\times 10^{3}$ CPU hours, approximately $95$ per cent of which is spent by grasil, with the remaining time being taken by galform to calculate the necessary star formation histories.

3 Results

In this Section, we present our main results. In Section 3.1 we present predictions for the evolution of physical properties of the galaxy population as well as a comparison of our predictions with available high-redshift ($z\gtrsim 7$) observational data. In Section 3.2 we present the predicted evolution of the galaxy luminosity function for the NIRCam–F200W and MIRI–F560W filters. We make such predictions for each NIRCam and MIRI broadband filter but only show these two in this paper for brevity; results for other filters will be made available online. In Section 3.3 we present predictions for galaxy number counts and redshift distributions (for a $10^{4}$ s exposure) observable by JWST in each NIRCam and MIRI band; we also show predictions for the redshift distributions of galaxies observable with longer ($10^{5}$ and $10^{6}$ s) exposures. Finally, in Section 3.4 we present predictions for the angular sizes of galaxies for the NIRCam–F200W and MIRI–F560W filters, again we make such predictions for all NIRCam filters but show only these two here for brevity. Throughout we show predictions for our fiducial model ‘lc16’ and the variant ‘lc16.EvolFB’ that adopts the evolving feedback model presented in Hou et al. (2016) and is discussed in Section 2.2. The dependence of our high-redshift predictions on some assumptions made in the model is discussed briefly in Appendix A.

3.1 The Lacey et al. (2016) model at high redshift

In this Section, we present model predictions for the evolution of some physical properties of the galaxy population and compare our predictions at $z\gtrsim 7$ to available observational data. In Fig. 3 we show predictions of the fiducial and evolving feedback variant models for the evolution of: (a) the galaxy stellar mass function; (b) the galaxy star formation rate function (for $M_{\star}>10^{6}$ $h^{-1}$ M*☉* galaxies); and (c) the fraction of bulge-dominated (i.e. with bulge-to-total stellar mass ratios of $B/T>0.5$) galaxies as a function of stellar mass, from $z=15.1$ to $z=0$.

The stellar mass function [Fig. 3 (a)] evolves rapidly at $z\gtrsim 2$ in both models. At lower redshifts, further evolution is predominantly at the high-mass end. It is easily seen that (for $z\gtrsim 2$) the evolving feedback model results in both more massive galaxies and a greater abundance of galaxies at a given stellar mass (for $M_{\star}\gtrsim 10^{6}$ $h^{-1}$ M*☉, as galaxies with a lower stellar mass are not included in our grasil sampling) by factors of up to $\sim 10$. For $z<4$, the normalisation of the supernova feedback strength is the same in both models and the differences between their stellar mass functions begin to disappear. At the low mass end ($M_{\star}\lesssim 10^{8}$ $h^{-1}$ M☉), however, the break in the power law for the mass loading factor (at $V_{\rm thresh}=50$ km s-1*) in the evolving feedback model results in a greater abundance of galaxies at these stellar masses than in the fiducial model. At the high mass end ($M_{\star}\gtrsim 10^{11}$ $h^{-1}$ M*☉*), an increase in stellar mass at low redshift due to the reduced feedback strength at higher redshift is apparent.

The distributions of star formation rates [Fig. 3 (b)] tell a similar story. For $z<4$ the distributions predicted by both models are essentially identical, except at low star formation rates (SFRs$\,\lesssim 10^{-2}$ $h^{-1}$ M*☉* yr*-1*) where the break in the evolving feedback model results in this model having a greater abundance of galaxies. At higher redshifts $z>4$ the differences in the star formation rate distributions are greater due to the different normalisations of feedback, with the evolving feedback variant having significantly more galaxies with SFRs$\,\gtrsim 3\times 10^{-2}$ $h^{-1}$ M*☉* yr*-1*. The apparent peak seen in each SFR distribution is mostly due to the imposed stellar mass limit of $10^{6}$ $h^{-1}$ M*☉*, if lower stellar mass galaxies were included it would shift to lower star formation rates according to the (approximately) constant relation between specific star formation rate and stellar mass predicted by the model (e.g. Mitchell et al., 2014; Cowley et al., 2017).

Fig. 3 (c) shows the evolution in the fraction of galaxies with a bulge-to-total stellar mass ratio of $B/T>0.5$, as a function of total stellar mass. In galform, bulges are created by a dynamical process, either a galaxy merger or a disc instability. The transition from a disc-dominated to a bulge-dominated galaxy population is relatively sharp, occurring over roughly one dex in stellar mass in most cases. In the evolving feedback model, this transition generally occurs at lower stellar masses. At higher redshifts (and thus lower stellar masses), the shape of this relationship is different for the evolving feedback variant, which predicts a much smoother transition. We caution against over-interpreting the predicted $B/T$ as a proxy for the morphological type. The instabilities that create bulges in galform do not necessarily create slowly rotating bulges, and so defining bulges as slow rotators would give different results to those presented here.

Having established some predicted physical properties of galaxies in the two models, we now compare predictions of the models to observational data at $z\gtrsim 7$. We note that none of the observational data considered here were used to calibrate model parameters [Lacey et al. (2016) only considered rest-frame far-UV luminosity functions at $z\lesssim 6$ in their model calibration].

We compare the predictions of the models for the evolution of the rest-frame far-UV luminosity function to observational data over the redshift interval $7\lesssim z\lesssim 10$ in Fig. 4. We can see that both models provide reasonable agreement with the observed data, and appear to ‘bracket’ the data for $M_{\rm AB}(1500\leavevmode\nobreak\ \AA)-5\log_{10}h\gtrsim-18$. However, at brighter magnitudes, the predictions of the two models converge. This is due to dust extinction becoming the limiting factor in a galaxy’s intrinsic brightness at far-UV wavelengths. To illustrate this, we show the predictions of the two models, without dust attenuation, in the $z=10$ panel. These predictions resemble the star formation rate distributions in Fig. 3 (b), as the star formation rate of a galaxy is essentially traced by the rest-frame far-UV.

Finally, we compare predictions for the angular sizes of galaxies to observational data in the redshift range $7\lesssim z\lesssim 9$ in Fig 5. The stellar component of the model galaxies is assumed to be a composite system, consisting of an exponential disc and a bulge with a projected $r^{1/4}$ density profile (Cole et al., 2000). We compute the half-light radii for our model galaxies by weighting the density profile of each component by their predicted rest-frame far-UV ($1500$ Å) luminosity, dividing the half-light radii of the disc by a factor of $1.34$ to account for inclination effects (Lacey et al., 2016), and interpolating to find the half-light radius of the composite system. We then bin the galaxies according to their flux, $S_{\nu}$. The symbols in Fig. 5 show the median size in each flux bin, with the error bars representing the $16-84$ percentile scatter in each bin. We show this for the whole galaxy population and also for starburst and quiescent galaxies. The differences between the predictions of the two models are small and they both show reasonable agreement with data from Ono et al. (2013) and Shibuya et al. (2015), who use GALFIT (Peng et al. 2002) to derive sizes from Hubble Space Telescope imaging. For the Ono et al. data we present their stacked image results. For the Shibuya et al. data we bin their sizes for individual galaxies into bins of $1$ mag width and present the median size in each bin. The error bars presented represent the $16-84$ percentile scatter of sizes within these bins. For reference, we also show the diffraction limit of JWST. The models predict that JWST should be able to resolve most galaxies in the rest-frame far-UV at these redshifts.

In summary, the predictions of both models show good agreement with the evolution of the rest-frame far-UV ($1500$ Å) luminosity function and observed galaxy sizes at high redshift ($z\gtrsim 7$). We re-iterate that these high-redshift data were not considered when calibrating the model.

3.2 Luminosity functions observable with JWST

In this Section we present predictions for the evolution of the galaxy luminosity function in the JWST NIRCam and MIRI bands. These are listed in Table 3, with their sensitivities (for a $10^{4}$ s exposure), and the field of view (FoV) for each instrument is shown in Table 4. In Fig. 6 we show the predicted luminosity functions for the NIRCam–F200W and MIRI–F560W bands. We make such predictions for all broadband NIRCam and MIRI filters, but show only these two here for brevity. The predictions for other filters will be made available online.

In the top panels of Fig. 6 we can see that at high redshifts the difference between the two models is similar to that seen in Fig. 4, and that the models predict similar luminosity functions for $z<4$, when the normalisation of the feedback strength is the same in both models.

In the bottom panels, we show the predicted luminosity function at $z=11$ for NIRCam–F200W (bottom left panel), and at $z=6$ for MIRI–F560W (bottom right panel). We choose these values as they are the redshifts at which we predict JWST will see $\sim 1$ object per field of view (FoV) for a $10^{4}$ s exposure, as is discussed below. Here we show the contribution to the luminosity function predicted by the fiducial model from quiescent and starburst galaxies. We can see that the bright end of the luminosity function is dominated by galaxies undergoing a burst of star formation. As mentioned earlier, the definition of starburst here refers to a dynamical process, either a galaxy merger or disc instability, triggering a period of enhanced star formation. In this case, the majority of the bursts are triggered by disc instabilities, as mergers appear to be inefficient at boosting the specific star formation rates of galaxies in this model, as is also discussed in Cowley et al. (2017). We also show predictions of the fiducial model without dust and can see that the bright end of the luminosity functions at these redshifts is composed of heavily dust-attenuated objects. We, therefore, expect such observations to provide a further constraint on the way dust absorption is accounted for in galaxy formation models.

For reference, we have also shown the sensitivity limits of the filters based on $10^{4}$ and $10^{5}$ s exposures as the vertical dashed and dotted lines respectively. Our adopted sensitivities for a $10^{4}$ s exposure are summarised in Table 3. We derive sensitivities for other exposures assuming they scale as $t^{-1/2}$.

In conjunction, we also show the abundance at which the instrument will see one object per FoV per unit redshift at this redshift. Our adopted fields of view are summarised in Table 4. Objects that are in the upper right quadrant of each plot would be observable with a $10^{4}$ s exposure in a single FoV. Therefore, the fiducial model predicts that $\sim 1$ object will be observable at $z=11$ by NIRCam–F200W, and $\sim 2$ will be observable at $z=6$ by MIRI–F560W. We recognise that single FoV observations will be sensitive to field-to-field variance. We hope to make direct predictions for the field-to-field variance by creating lightcone catalogues from our simulation in a future work.

3.3 Galaxy number counts and redshift distributions observable with JWST

The simplest statistic of a galaxy population that can be derived from an imaging survey is their number counts. Here we present the predictions for the cumulative number counts observable with NIRCam (Fig. 7) and MIRI (Fig. 8). We also show the corresponding redshift distributions (for a $10^{4}$ s exposure) in Fig. 9 (NIRCam) and Fig. 10 (MIRI). We obtain the number counts and redshift distributions by integrating the predicted luminosity functions according to

[TABLE]

where $\eta$ is the surface density of galaxies projected on the sky, $n$ is the number density of galaxies and $\mathrm{d}^{2}V/\mathrm{d}z\,\mathrm{d}\Omega$ is the comoving volume element per unit solid angle. We show the contribution to the predicted number counts and redshift distributions from quiescent and starburst galaxies. For the NIRCam filters, the counts are dominated by quiescent galaxies. This is because they are dominated by galaxies at low redshift, for which starbursts are not a significant population at these wavelengths. This is also why the predicted number counts from the fiducial and evolving feedback variant models are so similar, as at low redshifts the feedback normalisations are equal, though the lc16.EvolFB model does predict slightly more galaxies at faint fluxes. For the MIRI number counts, we see the burst population becoming important at brighter fluxes in bands $\lambda_{\rm obs}\gtrsim 10$ $\mu$m. These wavelengths also correspond to a shift from the number counts being dominated by dust-attenuated stellar light to dust emission. Again, these number counts are dominated by relatively low-redshift galaxies, for which the MIRI filters probe the dust emission from the rest-frame mid-IR.

The redshift distributions in Figs 9 and 10 exhibit a more discernible difference between the two models, particularly in the NIRCam bands at high redshift. For instance, in the NIRCam–F200W filter, the redshift at which one object per FoV per unit redshift is predicted to be observable with a $10^{4}$ s exposure is $z\sim 11$. For the evolving feedback variant $\sim 5$ times more galaxies are predicted to be observable at this redshift. From our predictions, it appears that very few galaxies will be observable at $z\gtrsim 10$ with NIRCam and at $z\gtrsim 6$ with MIRI, although we stress that this is the case for a single FoV and a $10^{4}$ s exposure. Additionally, we note that we have not considered effects such as gravitational lensing, which would allow surveys to probe fainter galaxies at higher redshifts (e.g. Infante et al., 2015).

Various features in the predicted MIRI redshift distributions can be related to PAH emission. For example, the peaks at $z\sim 2.5$ in the MIRI–F1130W distribution and at $z\sim 3.6$ in the MIRI-F1500W distribution correspond to the $3.3$ $\muup$m PAH feature.

We briefly consider the possibility that nebular emission lines may affect our predicted broadband photometry (e.g. Smit et al., 2015), as they are not included in our galaxy SEDs. For this we focus on the MIRI–F560W filter at $z\sim 7$ as the H $\alpha$ emission line is redshifted across the filter. The luminosity of the H $\alpha$ line is calculated assuming that all photons emitted with wavelengths shorter than $912$ Å will ionize a hydrogen atom in the gas surrounding the star. We then assume ‘Case B’ recombination i.e. we ignore recombinations directly to the ground state ($n=1$), as these just produce another ionizing photon. Thus only recombinations to $n>1$ are counted. The fraction of such recombinations that produce an H $\alpha$ photon ($n=2\rightarrow 1$) is taken from Osterbrock (1974). We apply the dust extinction factor predicted by grasil at the wavelength of the line to the line luminosity. We find that the predicted equivalent widths (EWs) of the line are $\sim 400$ Å, significantly narrower than the width of the MIRI–F560W filter $\sim 1.2$ $\muup$m. As a result, the line luminosity has a minor effect on the broadband photometry. For example, at $z=7.5$ in both models $95$ per cent of the sampled galaxies have their MIRI–F560W luminosity increased by less than $\sim 10$ per cent and $90$ per cent by less than $\sim 7$ per cent. This results in a negligible difference in the luminosity functions if H $\alpha$ emission is included. Thus we conclude that a more detailed inclusion of nebular emission lines (e.g. Panuzzo et al., 2003) is unlikely to affect the results presented here (see also Bisigello et al. 2016 for an investigation of the effect of nebular emission lines on MIRI photometry).

We now consider the predicted redshift distributions of galaxies that would be observable with longer exposures than considered in Figs 9 and 10. In Fig. 11 we show predictions for $10^{4}$, $10^{5}$ and $10^{6}$ s exposures, for the NIRCam–F200W and MIRI–F560W filters. For the fiducial model a $10^{6}$ s exposure will increase the number of observable objects in the NIRCam–F200W filter at $z\sim 11$ from $1$ per FoV to $\sim 10$ per Fov, and will increase the highest redshift at which an object is observable in a single FoV from $z\sim 11$ to $z\sim 13$. For the evolving feedback model, the highest redshift will be $z\sim 14.5$. Thus, we expect that long ($>10^{4}$ s) exposures with JWST will provide better constraints on the effectiveness of supernova feedback in galaxies at high redshift.

3.4 Sizes of galaxies in JWST bands

Finally, we present predictions for the angular sizes of galaxies for the NIRCam–F200W and MIRI-F560W filters in Fig. 12. We make such predictions for all NIRCam and MIRI filters but show only these two here for brevity, the predictions for other filters will be made available online. The sizes in each band are calculated as described in Section 3.1.

We can see that the predicted sizes are $\sim 0.1$ arcsec, with the evolving feedback variant generally predicting slightly smaller sizes. By comparison to the diffraction limits for JWST, shown here as dashed horizontal lines, it is evident that NIRCam will be able to resolve the majority of detected galaxies whereas this will not be the case for MIRI (for $z\gtrsim 2$).

4 Summary

The James Webb Space Telescope (JWST) is scheduled for launch in October 2018 and is expected to significantly advance our understanding of the high-redshift ($z\gtrsim 7$) Universe.

Here we present predictions for deep galaxy surveys with JWST. To do so we couple the hierarchical galaxy formation model galform (Lacey et al., 2016), with the spectrophotometric code grasil (Silva et al., 1998) for computing galaxy SEDs. grasil calculates the absorption and re-emission of stellar radiation by interstellar dust by solving the equations of radiative transfer in an assumed geometry. This allows us to produce UV-to-mm galaxy SEDs, broadening the predictive power of the model to cover the full wavelength range that will be probed by JWST. The galaxy formation model is implemented within a dark matter only $N$-body simulation using Planck cosmological parameters (Planck Collaboration et al., 2016). Adjustable parameters in the model are calibrated against a broad range of observational data such as optical and near-IR luminosity functions at $z=0$, the evolution of the rest-frame near-IR luminosity functions for $z=0-3$, far-IR galaxy number counts and redshift distributions, and the evolution of the rest-frame far-UV luminosity function for $z=3-6$ (Lacey et al. 2016; Baugh et al. in preparation). Here we have shown that the model predicts evolution of the rest-frame far-UV luminosity function for $7\lesssim z\lesssim 10$, and galaxy sizes for $7\lesssim z\lesssim 9$, in good agreement with observations.

We also present predictions for an evolving feedback variant model, in which the strength of supernova feedback is allowed to vary as a function of redshift (Hou et al., 2016). This adjustment allows the model to reproduce the reionization redshift inferred from Planck data (Planck Collaboration et al., 2016), as well as the $z=0$ luminosity function of the Milky Way satellites and their metallicity–stellar mass relation.

We present predictions for JWST in the form of luminosity functions, number counts, redshift distributions and angular sizes for each of the broadband filters on NIRCam and MIRI on JWST, for both the fiducial model ‘lc16’ and the evolving feedback variant ‘lc16.EvolFB’.

We find that for a $10^{4}$ s exposure the fiducial model predicts that JWST will be able to observe a single galaxy per field of view at $z\sim 11$ in the NIRCam–F200W filter; though the evolving feedback model predicts number surface densities factors of $\sim 5$ greater. The model predicts that similar exposures with MIRI will not detect any galaxies at $z\gtrsim 6$ (in a single FoV). Longer integration times will increase the number of galaxies that are observable, for example, a $10^{6}$ s integration will increase the number of galaxies predicted by the fiducial model to be observable in a single FoV by a factor of $\sim 10$. A similar effect may be achieved by utilising strong gravitational lenses; however, we do not consider such an effect here. We consider a simple model for calculating H $\alpha$ emission and conclude that nebular emission lines will have a negligible effect on these results.

The predicted sizes of high-redshift galaxies observable with JWST are $\sim 0.1$ arcsec, and as such we expect NIRCam to be capable of resolving the majority of detected galaxies.

We hope that the predictions presented here will help inform galaxy survey strategies for JWST. In the future, we plan to make our results public for such a purpose and to further develop our methodology to produce realistic mock galaxy catalogues for NIRCam and MIRI. This will allow us to make direct predictions for the field-to-field variance. We envisage that observations with JWST will provide a wealth of information on physical processes important for galaxy formation, such as the effectiveness of supernova feedback in galaxies at high redshift.

Acknowledgements

We thank the anonymous referee for their comprehensive and constructive comments. The authors would like to thank Alessandro Bressan, Gian-Luigi Granato and Laura Silva for use of, and discussions relating to, the grasil code. Additionally, we would like to acknowledge Karina Caputi, Lydia Heck, John Helly, Hou Jun, Sarah Kendrew, John Pye, Luiz Felippe Rodrigues, Martin Ward and Rogier Windhorst for many helpful discussions during the development of this work. This work was supported by the Science and Technology Facilities Council [ST/K501979/1, ST/L00075X/1]. CMB acknowledges the receipt of a Leverhulme Trust Research Fellowship. CSF acknowledges an ERC Advanced Investigator Grant, COSMIWAY [GA 267291] and the Science and Technology Facilities Council [ST/F001166/1, ST/I00162X/1]. This work used the DiRAC Data Centric system at Durham University, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by BIS National E-infrastructure capital grant ST/K00042X/1, STFC capital grant ST/H008519/1, and STFC DiRAC Operations grant ST/K003267/1 and Durham University. DiRAC is part of the National E-Infrastructure.

Appendix A Dependence of high-redshift predictions on model assumptions

The galaxy formation model used in this work assumes prescriptions for the physical processes involved in galaxy formation that are in some cases calibrated on observational, or simulation, data at $z\lesssim 6$. In this Appendix, we briefly discuss the dependence of our results on some of the assumptions made and their applicability at the high redshifts explored in this paper. We note that this is not intended to be an exhaustive discussion (which we consider to be beyond the scope of this work), and emphasize that the primary focus of this study is to present the predictions for the outcomes of future JWST galaxy surveys for the model described in Lacey et al. (2016), rather than to investigate to what extent these predictions are sensitive to the various prescriptions for physical processes within the model. We envisage that this will form the basis of future study once the predictions presented earlier are confronted with the corresponding observations.

Some results from a brief parameter space exploration at high redshift are shown in Fig. 13, for a more detailed exploration of the parameter space of this model we refer to reader to Appendix C of Lacey et al. (2016). For ease of computation, the results presented in Fig. 13 have used the model for the absorption and re-emission of radiation by interstellar dust described in Lacey et al. (2016, see also Cowley et al., 2017), rather than grasil. However, for the rest-frame wavelength investigated in Fig. 13 ($1500$ Å) the difference between the predictions of these two models is negligible.

A.1 Formation of molecular hydrogen and star formation efficiency

The model assumes a star formation law for galactic discs in which the surface density of star formation is related to the surface density of molecular gas (see equation 2). The molecular gas fraction, $f_{\rm mol}$, is computed according to the mid-plane gas pressure in the disc using the empirical relation of Blitz & Rosolowsky (2006). Any possible dependence of the molecular gas fraction on physical properties such as the gas metallicity or the interstellar far-UV radiation field as proposed in some models (e.g. Gnedin et al., 2009; Krumholz et al., 2009) are not explicitly considered. The model for the formation of molecular hydrogen of Krumholz et al. (2009) was implemented into galform by Lagos et al. (2011), however, Lagos et al. concluded that the Blitz & Rosolowsky (2006) relation produced better agreement with observations at low redshift.

The value for the inverse timescale of star formation from molecular gas, sometimes referred to as the star formation efficiency, $\nu_{\rm SF}$ (see equation 2), is chosen based on the value measured for this parameter by Bigiel et al. (2011). Those authors presented the $1\sigma$ range for this quantity as being $0.25-0.74$ Gyr*-1*, so in practice, this parameter is only varied within this range during the model calibration. Doing so has a negligible effect on our predictions for the rest-frame far-UV ($1500\leavevmode\nobreak\ {\AA}$) luminosity function at $z=10$ [see Fig. 13(a)].

A.2 Feedback processes

Throughout this work, we have explored a variant model in which the strength of SN feedback is allowed to vary as a function of redshift. This variant allows the model to reproduce a higher redshift of reionization in better agreement with that inferred from Planck data, as well as the luminosity function and metallicity–stellar mass relation of $z=0$ Milky Way satellite galaxies, though it does introduce some new parameters in the SN feedback model (Hou et al., 2016). We find that this modification increases the numbers of galaxies observable with JWST at high-redshift by factors of $\sim 5$.

The prescription for AGN feedback in our model was first introduced by Bower et al. (2006) and has a minimal impact on model predictions at high redshifts. In this prescription, the conditions for gas cooling in the halo to be turned off as a result of thermal energy deposited by relativistic radio jets are: (i) that the gas is cooling quasi-statically (i.e. the cooling time is long compared to the free-fall time); and (ii) that the SMBH is massive enough such that the power required to balance the radiative cooling luminosity is below some fraction of its Eddington luminosity. These conditions are rarely met at high-redshift as: (i) gas cooling times are generally shorter as the physical density of the Universe is higher ($\tau_{\rm cool}\propto\rho^{-1}$); and (ii) typically the SMBHs have not had time to grow massive enough to satisfy the second criterion. As such, turning off AGN feedback has a negligible impact on our prediction for the rest-frame UV luminosity function at $z=10$ [see Fig. 13(b)]. However, we stress that turning off AGN feedback does not result in a viable model of galaxy formation, as it fails to reproduce many properties of the galaxy population at $z=0$ (see e.g. Fig. C2 of Lacey et al., 2016).

Photoionization feedback is implemented such that for $z<z_{\rm reion}$, no cooling of gas occurs in halos with circular velocities $V_{\rm c}<V_{\rm crit}$ (Benson et al., 2002; Font et al., 2011). Adjusting $z_{\rm reion}$ has a negligible impact on our high-redshift predictions at luminosities observable by JWST, and adjusting $V_{\rm crit}$ mainly shifts the ‘break’ at the faint end of the luminosity function [see Fig. 13(c)]. It should be noted, however, that these parameters are considered ‘fixed’ (see Table 1 of Lacey et al., 2016) and are not varied in the model calibration. The value for $z_{\rm reion}$ is chosen based on WMAP7 data (Dunkley et al., 2009) and the value for $V_{\rm crit}$ is based on cosmological gas simulations (Hoeft et al., 2006; Okamoto et al., 2008).

A.3 Disc instabilities and top-heavy IMF

The bright ends of the predicted NIRCam and MIRI luminosity functions presented in this work are predicted to be dominated by galaxies undergoing a disc instability-triggered starburst for $z\gtrsim 2$, so they will be somewhat sensitive to the treatment of this process. For example, turning off disc instabilities shifts the bright end of the rest-frame far-UV luminosity function towards fainter magnitudes by $\sim 2$ magnitudes [see Fig. 13(b)]. We stress very strongly, however, that turning this process off does not result in a viable model of galaxy formation, as it no longer reproduces the observed far-IR/sub-mm galaxy number counts and redshift distributions, which are predominantly composed of disc-instability triggered starburst galaxies in the redshift range $z\sim 1-4$ (Lacey et al., 2016).

A related issue is the assumption of a mildly top-heavy IMF for star formation during a starburst episode. This assumption is made so that the model can reproduce the observed far-IR/sub-mm galaxy number counts and redshift distributions (Baugh et al., 2005; Cowley et al., 2015; Lacey et al., 2016), though the slope for the IMF adopted here is much less top-heavy than was advocated by Baugh et al. Baugh et al. suggested a slope of $x=0$ in $\mathrm{d}N/\mathrm{d}\ln m\propto m^{-x}$, however in the updated galform model used here, $x=1$ was found to produce better agreement with the observed far-IR/sub-mm galaxy number counts and redshift distributions (Lacey et al., 2016).

This assumption has only a fairly minor impact on our model predictions in the rest-frame far-UV, as is shown in Fig. 13(d), where we compare the predictions of the fiducial model to a model that adopts a Kennicutt (1983) IMF for both disc and burst star formation (the fiducial model adopts a Kennicutt IMF for disc star formation only). The increase in intrinsic UV flux as a result of forming more massive stars is counterbalanced by the increase in dust mass (which absorbs the UV radiation) due to the increased rate of SNe returning metals to the ISM. We emphasize this point by also showing the model predictions without dust attenuation in this panel, the difference between the two predicted luminosity functions is greater without dust extinction. As a result, the rest-frame far-UV luminosity function is fairly insensitive to the assumption of a top-heavy IMF for burst star formation, which has a much greater effect in far-IR/sub-mm bands that trace the dust emission, against which the model has been calibrated.

A.4 Stellar population synthesis model

The models presented here make use of the stellar population synthesis (SPS) model of Maraston (2005). This includes an empirical calibration for the light produced by thermally-pulsating asymptotic giant branch stars (TP-AGB). This is difficult to model accurately from purely theoretical stellar evolution models, so Maraston (2005) calibrate this using observations of star clusters. This phase is important for the rest-frame near-IR luminosities of stellar populations with ages $\sim 0.1-1$ Gyr. We refer the interested reader to Conroy (2013) for a review of the Maraston, and other, SPS models. A comprehensive investigation of the effects of implementing different SPS models within galform was performed by Gonzalez-Perez et al. (2014), who found that SPS models with a larger contribution from TP-AGB stars produced stronger evolution in the rest-frame near-IR luminosity function. Lacey et al. (2016) found that implementing different SPS models made only a modest change to the predicted rest-frame $K$-band luminosity function, but that the Maraston (2005) models produced better agreement with observations of this quantity up to $z\sim 3$ (see Fig. C16 of Lacey et al., 2016).

Recently, Stanway et al. (2016) proposed an SPS model that includes the effects of binary stars (bpass), which produces higher hydrogen-ionizing luminosities than single-star SPS models by up to $\sim 60$ per cent at sub-solar ($0.1-0.2$ Z*⊙*) metallicities. Wilkins et al. (2016) found that the increased ionizing radiation that results from assuming this SPS model also increases the nebular emission by up to $\sim 25$ per cent. It is conceivable that assuming the Stanway et al. models would have a similar effect on predictions from our simple model for nebular emission. However, given that we estimate the impact of nebular emission on our predictions for the broadband luminosity functions presented in this work to be negligible, it is unlikely that an increase of only $25$ per cent in nebular emission would produce a significant effect on our predicted broadband luminosity functions.