Simulating galaxies in the reionization era with FIRE-2: galaxy scaling relations, stellar mass functions, and luminosity functions

TL;DR

This paper uses cosmological simulations to study galaxy properties during reionization, revealing how stellar mass, luminosity functions, and dust attenuation evolve from redshift 5 to 12, and providing fitting functions for key relations.

Contribution

It presents new high-resolution simulations covering a wide halo mass range, offering detailed predictions and fitting functions for galaxy scaling relations during reionization.

Findings

Stellar mass-halo mass relation remains relatively stable from z=5 to 12.

Faint-end luminosity function steepens with increasing redshift.

Dust obscuration significantly affects UV luminosity at z~6.

Abstract

We present a suite of cosmological zoom-in simulations at z>5 from the Feedback In Realistic Environments project, spanning a halo mass range M_halo~10^8-10^12 M_sun at z=5. We predict the stellar mass-halo mass relation, stellar mass function, and luminosity function in several bands from z=5-12. The median stellar mass-halo mass relation does not evolve strongly at z=5-12. The faint-end slope of the luminosity function steepens with increasing redshift, as inherited from the halo mass function at these redshifts. Below z~6, the stellar mass function and ultraviolet (UV) luminosity function slightly flatten below M_star~10^4.5 M_sun (fainter than M_1500~-12), owing to the fact that star formation in low-mass halos is suppressed by the ionizing background by the end of reionization. Such flattening does not appear at higher redshifts. We provide redshift-dependent fitting functions for the SFR-M_halo, SFR-M_star, and broad-band magnitude-stellar mass relations. We derive the star formation rate density and stellar mass density at z=5-12 and show that the contribution from very faint galaxies becomes more important at z>8. Furthermore, we find that the decline in the z~6 UV luminosity function brighter than M_1500~-20 is largely due to dust attenuation. Approximately 37% (54%) of the UV luminosity from galaxies brighter than M_1500=-13 (-17) is obscured by dust at z~6. Our results broadly agree with current data and can be tested by future observations.

Tables

Table 1. Table 1: Simulation details. Each simulation contains several galaxies in the zoom-in region. Properties below refer to the most massive (or the ‘target’) galaxy.

Name	$M_{\mathrm{halo}}$ ($z=5$)	$M_{\ast }$ ($z=5$)	$M_{\mathrm{halo}}$ ($z=10$)	$M_{\ast }$ ($z=10$)	$m_b$	${ϵ}_b$	${ϵ}_{\mathrm{star}}$	$m_{\mathrm{DM}}$	${ϵ}_{\mathrm{DM}}$
	[$M_{☉}$]	[$M_{☉}$]	[$M_{☉}$]	[$M_{☉}$]	[$M_{☉}$]	[pc]	[pc]	[$M_{☉}$]	[pc]
z5m09a	2.4e+09	8.0e+05	2.3e+08	2.3e+04	119.3	0.14	0.7	6.5e+02	10
z5m09b	2.8e+09	5.9e+05	1.5e+08	5.2e+03	119.3	0.14	0.7	6.5e+02	10
z5m10a	6.7e+09	1.0e+07	1.4e+09	3.5e+05	954.4	0.28	1.4	5.2e+03	21
z5m10b	1.2e+10	1.6e+07	1.6e+09	1.1e+06	954.4	0.28	1.4	5.2e+03	21
z5m10c	1.3e+10	1.1e+07	1.0e+09	3.5e+05	954.4	0.28	1.4	5.2e+03	21
z5m10d	1.9e+10	2.0e+07	2.9e+08	2.4e+04	954.4	0.28	1.4	5.2e+03	21
z5m10e	2.4e+10	1.7e+07	6.8e+08	3.8e+05	954.4	0.28	1.4	5.2e+03	21
z5m10f	3.2e+10	1.1e+08	1.6e+09	3.4e+05	954.4	0.28	1.4	5.2e+03	21
z5m11a	4.1e+10	2.8e+07	1.4e+09	3.3e+05	954.4	0.28	1.4	5.2e+03	21
z5m11b	4.0e+10	9.2e+07	3.3e+09	3.4e+06	890.8	0.28	1.4	4.9e+03	21
z5m11c	7.5e+10	4.5e+08	8.5e+09	1.2e+07	7126.5	0.42	2.1	3.9e+04	42
z5m11d	1.4e+11	9.9e+08	2.5e+10	1.0e+08	7126.5	0.42	2.1	3.9e+04	42
z5m11e	2.4e+11	1.1e+09	1.3e+10	5.2e+07	7126.5	0.42	2.1	3.9e+04	42
z5m12a	4.4e+11	3.0e+09	2.3e+10	4.7e+07	7126.5	0.42	2.1	3.9e+04	42
z5m12b	8.5e+11	1.5e+10	3.2e+10	1.0e+08	7126.5	0.42	2.1	3.9e+04	42

Table 2. Table 2: Best-fit parameters for the magnitude–stellar mass relation and the stellar mass–magnitude relation (see Section 3.4 for details).

Magnitude–stellar mass relation (Equation 6)
Band		$a$	$c$	$d$
1500 Å	median	-2.81	-5.61	-22.38
1500 Å	1$\sigma$ lower	-2.61	-6.83	-23.06
1500 Å	1$\sigma$ upper	-2.74	-3.87	-21.42
B band	median	-2.63	-3.36	-22.46
B band	1$\sigma$ lower	-2.59	-5.17	-22.89
B band	1$\sigma$ upper	-2.64	-2.52	-22.05
J band	median	-2.61	-2.63	-22.69
J band	1$\sigma$ lower	-2.61	-3.86	-23.20
J band	1$\sigma$ upper	-2.63	-2.15	-22.31
Stellar mass–magnitude relation (Equation 7)
Band		$a^{\prime }$	$c^{\prime }$	$d^{\prime }$
1500 Å	median	-0.39	-2.59	8.77
1500 Å	1$\sigma$ lower	-0.42	-1.65	8.38
1500 Å	1$\sigma$ upper	-0.36	-2.99	9.16
B band	median	-0.38	-2.17	8.95
B band	1$\sigma$ lower	-0.41	-1.59	8.66
B band	1$\sigma$ upper	-0.37	-1.40	9.09
J band	median	-0.38	-1.85	8.90
J band	1$\sigma$ lower	-0.41	-1.61	8.67
J band	1$\sigma$ upper	-0.38	-0.90	9.01

Table 3. Table 3: Stellar mass functions and luminosity functions at rest-frame 1500 Å, B, and J band from z = 5 𝑧 5 z=5 –12.

Stellar mass function		Luminosity function
$\log M_{\ast }$	$\log {\varphi }_{\ast }$	M1500	$\log {\varphi }_{1500}$	MB	$\log {\varphi }_B$	MJ	$\log {\varphi }_J$
($M_{☉}$)	(Mpc-3 dex-1)	(mag)	(Mpc-3 mag-1)	(mag)	(Mpc-3 mag-1)	(mag)	(Mpc-3 mag-1)
$z=5$
3.83	0.76	-22.37	-3.96	-22.18	-3.95	-22.60	-4.05
4.50	0.62	-20.68	-3.01	-20.52	-2.94	-20.89	-3.04
5.17	0.33	-18.99	-2.42	-18.85	-2.37	-19.17	-2.36
5.84	-0.04	-17.30	-1.84	-17.18	-1.73	-17.45	-1.76
6.51	-0.47	-15.61	-1.14	-15.51	-1.05	-15.73	-1.08
7.18	-0.92	-13.92	-0.69	-13.84	-0.67	-14.01	-0.66
7.85	-1.62	-12.23	-0.38	-12.17	-0.22	-12.30	-0.20
8.52	-2.35	-10.54	-0.09	-10.50	0.11	-10.58	0.15
9.19	-2.55	-8.85	0.20	-8.83	0.42	-8.86	0.47
9.85	-4.02
$z=6$
3.83	0.87	-22.34	-4.10	-22.01	-4.18	-22.34	-4.19
4.48	0.62	-20.66	-3.35	-20.37	-3.35	-20.66	-3.35
5.13	0.35	-18.97	-2.64	-18.72	-2.51	-18.97	-2.64
5.78	-0.09	-17.28	-1.84	-17.07	-1.81	-17.28	-1.86
6.43	-0.56	-15.59	-1.25	-15.42	-1.27	-15.59	-1.29
7.08	-1.12	-13.91	-0.80	-13.77	-0.70	-13.91	-0.70
7.74	-1.66	-12.22	-0.41	-12.12	-0.30	-12.22	-0.27
8.39	-2.32	-10.53	-0.06	-10.47	0.11	-10.53	0.13
9.04	-2.88	-8.84	0.24	-8.83	0.42	-8.84	0.48
9.69	-3.98
$z=7$
3.81	0.97	-21.39	-3.89	-21.07	-4.14	-21.40	-4.34
4.42	0.60	-19.82	-3.17	-19.54	-3.05	-19.82	-3.23
5.04	0.32	-18.24	-2.54	-18.00	-2.45	-18.24	-2.47
5.65	-0.16	-16.67	-1.92	-16.46	-1.88	-16.67	-1.92
6.26	-0.63	-15.09	-1.26	-14.92	-1.25	-15.09	-1.32
6.88	-1.20	-13.52	-0.75	-13.38	-0.74	-13.52	-0.75
7.49	-1.68	-11.94	-0.37	-11.85	-0.26	-11.94	-0.25
8.11	-2.32	-10.36	-0.04	-10.31	0.05	-10.37	0.13
8.72	-2.99	-8.79	0.26	-8.77	0.46	-8.79	0.50
9.33	-4.00
$z=8$
3.86	0.87	-21.46	-4.31	-20.90	-4.22	-20.91	-4.22
4.57	0.48	-19.39	-3.21	-18.92	-2.95	-18.93	-2.96
5.29	-0.01	-17.32	-2.25	-16.93	-2.27	-16.94	-2.23
6.00	-0.61	-15.25	-1.43	-14.95	-1.45	-14.95	-1.45
6.72	-1.31	-13.18	-0.72	-12.96	-0.70	-12.97	-0.68
7.43	-1.89	-11.11	-0.17	-10.98	-0.14	-10.98	-0.11
8.15	-2.64	-9.04	0.32	-8.99	0.45	-8.99	0.47
8.86	-3.89
$z=9$
3.83	0.77	-20.85	-4.95	-20.27	-4.93	-20.33	-4.93
4.50	0.39	-18.87	-3.14	-18.39	-3.16	-18.43	-3.17
5.17	-0.11	-16.90	-2.38	-16.50	-2.25	-16.53	-2.27
5.84	-0.69	-14.92	-1.52	-14.61	-1.51	-14.64	-1.56
6.50	-1.38	-12.94	-0.80	-12.72	-0.83	-12.74	-0.84
7.17	-2.00	-10.97	-0.20	-10.83	-0.22	-10.85	-0.22
7.84	-2.69	-8.99	0.20	-8.94	0.36	-8.95	0.36
8.51	-3.66
$z=10$
3.80	0.68	-19.25	-3.73	-18.50	-3.60	-18.44	-3.58
4.41	0.32	-17.52	-2.90	-16.88	-2.79	-16.84	-2.70
5.01	-0.12	-15.79	-2.18	-15.27	-2.15	-15.23	-2.23
5.62	-0.74	-14.06	-1.41	-13.65	-1.40	-13.62	-1.40
6.22	-1.39	-12.33	-0.75	-12.04	-0.78	-12.02	-0.78
6.83	-1.98	-10.60	-0.21	-10.42	-0.23	-10.41	-0.18
7.43	-2.54	-8.87	0.16	-8.81	0.29	-8.81	0.27
8.04	-3.44

Table 4. Table 4: Table 3 — continued.

Stellar mass function		Luminosity function
$\log M_{\ast }$	$\log {\varphi }_{\ast }$	M1500	$\log {\varphi }_{1500}$	MB	$\log {\varphi }_B$	MJ	$\log {\varphi }_J$
($M_{☉}$)	(Mpc-3 dex-1)	(mag)	(Mpc-3 mag-1)	(mag)	(Mpc-3 mag-1)	(mag)	(Mpc-3 mag-1)
$z=11$
3.82	0.52	-19.10	-3.85	-18.36	-3.82	-18.07	-3.78
4.46	0.21	-17.08	-3.09	-16.48	-2.95	-16.24	-2.91
5.10	-0.44	-15.06	-2.00	-14.60	-2.00	-14.41	-1.83
5.74	-1.10	-13.05	-1.21	-12.71	-1.32	-12.58	-1.27
6.38	-1.87	-11.03	-0.38	-10.83	-0.42	-10.75	-0.39
7.03	-2.79	-9.01	0.06	-8.94	0.16	-8.92	0.15
7.67	-2.90
$z=12$
3.81	0.37	-18.64	-4.03	-17.85	-3.84	-17.52	-3.79
4.42	0.07	-16.71	-3.22	-16.06	-3.09	-15.79	-3.02
5.04	-0.52	-14.77	-2.10	-14.27	-2.15	-14.06	-2.04
5.65	-1.33	-12.84	-1.38	-12.48	-1.41	-12.33	-1.39
6.26	-1.97	-10.90	-0.49	-10.69	-0.54	-10.60	-0.47
6.88	-2.81	-8.97	-0.08	-8.90	0.02	-8.87
7.49	-3.25

Equations

$$\log {\mathrm{M}}_{\ast }=\alpha (\log {\mathrm{M}}_{\mathrm{halo}}-10)+\beta$$

$${\sigma }_{\log {\mathrm{M}}_{\ast }}=\exp [\gamma (\log {\mathrm{M}}_{\mathrm{halo}}-10)+\delta ],$$

$$(\alpha ,\beta ,\gamma ,\delta )=(1.58,7.10,-0.14,-1.10).$$

$$\log \mathrm{SFR}=\alpha (\log {\mathrm{M}}_{\mathrm{halo}}-10)+\gamma \log \left(\frac{1+z}{6}\right)+\delta$$

$$\log \mathrm{SFR}={\alpha }^{\prime }(\log {\mathrm{M}}_{\ast }-10)+{\gamma }^{\prime }\log \left(\frac{1+z}{6}\right)+{\delta }^{\prime },$$

$${\mathrm{M}}_{\mathrm{AB},\mathrm{band}}=a(\log {\mathrm{M}}_{\ast }-10)+c\log \left(\frac{1+z}{6}\right)+d,$$

$$\log {\mathrm{M}}_{\ast }=a^{\prime }({\mathrm{M}}_{\mathrm{AB},\mathrm{band}}+20)+c^{\prime }\log \left(\frac{1+z}{6}\right)+d^{\prime },$$

$${\mathrm{A}}_{1500}=(0.0306\pm 0.0002)({\mathrm{M}}_{1500}+15{)}^2,$$

Full text

††pagerange: Simulating galaxies in the reionization era with FIRE-2: galaxy scaling relations, stellar mass functions, and luminosity functions–C

1 Introduction

High-redshift galaxies are believed to be the dominant sources contributing to cosmic reionization (e.g. Faucher-Giguère et al. 2008; Haardt & Madau 2012; Kuhlen & Faucher-Giguère 2012; Robertson et al. 2013; Robertson et al. 2015; however, see Madau & Haardt 2015). Current deep surveys using the Hubble Space Telescope have already put reliable constraints on the $z\geq 5$ ultraviolet (UV) luminosity functions for galaxies brighter than $\rm M_{UV}=-17$ (e.g. McLure et al., 2013; Schenker et al., 2013; Bouwens et al., 2015; Finkelstein et al., 2015a), but the faint-end behavior of the UV luminosity function remains highly uncertain. These faint galaxies contribute a non-trivial fraction of the ionizing photons needed for reionization (e.g. Finkelstein et al., 2012; Kuhlen & Faucher-Giguère, 2012; Robertson et al., 2013), although their abundances are poorly understood.

Recently, Livermore et al. (2017) reported the detection of very faint galaxies of $\rm M_{UV}=-12.5$ at $z\sim 6$ that are highly magnified by galaxy clusters in the Hubble Frontier Fields, after performing a novel analysis to remove the cluster light. They found a steep UV luminosity function down to $\rm M_{UV}=-13$ at $z\geq 6$, implying sufficient numbers of faint galaxies to account for cosmic reionization. However, Bouwens et al. (2017b); Bouwens et al. (2017a) later pointed out that the uncertain size distribution of high-redshift galaxies and the uncertain magnification model of the lensing clusters can have a large impact on the inferred faint-end luminosity functions in the Hubble Frontier Fields. The faint-end slope of the UV luminosity function fainter than $\rm M_{UV}=-15$ thus remains poorly constrained.

Great efforts have also been made to measure the galaxy stellar mass functions at these redshifts (e.g. González et al., 2011; Duncan et al., 2014; Grazian et al., 2015; Song et al., 2016; Stefanon et al., 2017). The stellar masses of high-redshift galaxies are usually derived from single-band photometry using empirical relations. Such relations are calibrated against spectral energy distribution (SED) fitting using limited rest-frame optical data for a small sample of galaxies at these redshifts. These relations tend to have large intrinsic scatter and suffer from systematic uncertainties of the underlying stellar population synthesis model. Therefore, the stellar mass functions reported by different authors have considerable discrepancies (e.g. figure 9 in Song et al. 2016).

Consequently, the stellar mass–halo mass relation and the star formation efficiencies inferred from the stellar mass measurements at these redshifts are also very uncertain. For example, Finkelstein et al. (2015b) reported an increasing stellar mass to halo mass ratio with increasing redshift, whereas Stefanon et al. (2017) found no evolution of this ratio at these redshifts. Another related question is to understand the stellar mass growth histories of galaxies at these redshifts. This is not only useful for constraining the total ionizing photon emissivity at the epoch of reionization, but also essential for understanding galaxy populations at lower redshift – both dwarf galaxy abundances in the Local Group (e.g. Boylan-Kolchin et al., 2015) and stellar mass functions in local galaxy clusters (e.g. Lu et al., 2014a).

The James Webb Space Telescope (JWST, scheduled for launch in 2020) and the next generation of ground-based telescopes will make it possible to study $z\geq 5$ galaxies in more detail. Future observations of galaxies in the reionization era will provide substantial data for high-spatial-resolution deep imaging at the rest-frame optical bands, as well as spectroscopic measurements probing the physical conditions of the interstellar medium (ISM) in these galaxies. This may help resolve many current open questions in the field, such as the faint-end slope of the luminosity function, more robust determination of stellar mass, understanding the stellar populations in high-redshift galaxies and their contribution to cosmic reionization (Leitherer et al., 2014; Topping & Shull, 2015; Choi et al., 2017; Stanway, 2017), etc. Therefore, it is necessary from a theoretical point of view to make more realistic predictions of galaxy properties at these redshifts.

Currently there are two broad categories of cosmological simulations of galaxy formation at the epoch of reionization. High-resolution cosmological radiation-hydrodynamic simulations, with a detailed set of baryonic physics, including primordial chemistry and molecular networks, can simultaneously model the formation of first stars and galaxies and the local reionization history (e.g. Wise et al., 2014; Chen et al., 2014; O’Shea et al., 2015; Paardekooper et al., 2015). Such calculations are usually computationally expensive and thus carried out in a small cosmological volume. They generally focus on the formation of Population III (Pop III) stars and low-mass galaxies (in halos below ${\rm M}_{\rm halo}\sim 10^{9}\,{\rm M}_{\sun}$) at relatively high redshifts ($z\gtrsim 10$). These types of simulations have been used to predict the scaling relations of high-redshift, low-mass galaxies (e.g. the stellar mass–halo mass relation, gas fraction, mass–metallicity relation, etc.; Chen et al., 2014), ionizing photon escape fractions from these small galaxies and their importance for cosmic reionization (e.g. Paardekooper et al., 2015; Xu et al., 2016), their spectral properties and detectability with JWST (e.g. Barrow et al., 2017), and the faint-end ($\rm M_{UV}>-14$) UV luminosity functions at these redshifts (e.g. O’Shea et al., 2015).

On the other hand, there are also large-volume cosmological simulations at relatively low resolution using empirically-calibrated models of star formation and stellar feedback (e.g. Feng et al., 2016; Gnedin, 2016; Ocvirk et al., 2016; Finlator et al., 2017; Pawlik et al., 2017). Simulations of this nature broadly reproduce the observed galaxy populations, stellar mass functions, UV luminosity functions (e.g. Gnedin, 2016; Wilkins et al., 2017), and the global reionization histories (e.g. Ocvirk et al., 2016; Pawlik et al., 2017). Forward modeling of galaxies in these simulations provide large samples of mock images and spectra that can be directly confronted with JWST (e.g. Wilkins et al., 2016; Zackrisson et al., 2017). However, these simulations tend to have mass resolution $\gtrsim 10^{5}\,{\rm M}_{\sun}$. Therefore, they are not able to capture the small-scale physics and the detailed structures in galaxies, which can be important for questions such as understanding the escape of ionizing photons (e.g. Ma et al., 2015). Also, some galaxy formation models calibrated to observations in the local universe struggle to reproduce observed galaxy properties at intermediate redshifts ($z\sim 2$–3), such as star formation histories, metallicities, etc. (e.g. Ma et al., 2016a; Davé et al., 2016). This is also a known problem in semi-analytic models of galaxy formation (e.g. Lu et al., 2014b).

In this work, we introduce a new suite of cosmological ‘zoom-in’ simulations at $z\geq 5$ in the $z=5$ halo mass range ${\rm M}_{\rm halo}\sim 10^{8}$–$10^{12}\,{\rm M}_{\sun}$. We mainly focus on relatively massive (above the atomic cooling limit), Population II (Pop II) star-dominated galaxies in the redshift range $z=5$–12. Our simulations cover a range of galaxies that can be well probed by future observations using JWST and next-generation ground-based telescopes. The cosmological zoom-in technique allows us to simulate galaxies in a broad mass range without being limited to a fixed simulation volume. The resolution is adaptively chosen based on the mass of the system, but always much better than that of large-volume simulations. These are not the first cosmological zoom-in simulations at $z\geq 5$: previous works using a similar technique have studied the escape fraction of ionizing photons (e.g. Kimm & Cen, 2014), galaxy properties and scaling relations (e.g. Ceverino et al., 2017), and the importance of stellar feedback for shaping these galaxies (e.g. Yajima et al., 2017). Our work builds on these recent studies by increasing the resolution, expanding sample size, and most importantly including more detailed treatments for stellar feedback.

Our high-resolution cosmological zoom-in simulations use a full set of physically motivated models of the multi-phase ISM, star formation, and stellar feedback from the Feedback in Realistic Environments (FIRE) project111http://fire.northwestern.edu. In a series of previous papers, these models have shown to successfully reproduce a variety of observed galaxy properties at lower redshifts (e.g. Hopkins et al., 2017, and references therein). Therefore, the new simulations presented in this paper are complementary to other state-of-the-art simulations in the field of galaxies in the reionization era.

This paper is the first in a series based on these new simulations, focusing on galaxy properties, scaling relations, stellar mass functions, and luminosity functions at $z>5$. Our results complement previous predictions on the same topics using semi-analytic models (e.g. Clay et al., 2015; Liu et al., 2016; Cowley et al., 2018) and cosmological simulations (e.g. Jaacks et al., 2012; O’Shea et al., 2015; Yajima et al., 2015; Gnedin, 2016; Ocvirk et al., 2016; Xu et al., 2016; Wilkins et al., 2017). In Sections 2.1 and 2.2, we describe the initial conditions and the physical ingredients used in the code. In Sections 2.3 and 2.4, we construct the simulated catalog. In Section 3, we present the general properties of our simulated galaxies. In Sections 4 and 5, we predict the stellar mass functions and luminosity functions from $z=5$–12. We discuss our results in Section 6 and conclude in Section 7.

We adopt a standard flat $\Lambda$CDM cosmology with Planck 2015 cosmological parameters $H_{0}=68{\rm\,km\,s^{-1}\,Mpc^{-1}}$, $\Omega_{\Lambda}=0.69$, $\Omega_{m}=1-\Omega_{\Lambda}=0.31$, $\Omega_{b}=0.048$, $\sigma_{8}=0.82$, and $n=0.97$ (Planck Collaboration et al., 2016). In this paper, we adopt a Kroupa (2002) initial mass function (IMF) from 0.1–$100\,{\rm M}_{\sun}$, with IMF slopes of $-1.30$ from 0.1–$0.5\,{\rm M}_{\sun}$ and $-2.35$ from 0.5–$100\,{\rm M}_{\sun}$. All magnitudes are in the AB system (Oke & Gunn, 1983).

2 The Simulations

The simulations presented in this paper form a subsample of the FIRE project (version 2.0, which we refer as FIRE-2; Hopkins et al., 2017). FIRE-2 is an updated version of the feedback implementations studied in a number of previous papers, which we refer as FIRE-1 (Hopkins et al., 2014).

All FIRE-2 simulations are run using an identical version of the gizmo code (Hopkins, 2015)222http://www.tapir.caltech.edu/~phopkins/Site/GIZMO.html. We use the meshless finite-mass (MFM) method in gizmo to solve the hydrodynamic equations. We refer to Hopkins et al. (2017) for details of the numerical methods and convergence tests of the FIRE-2 simulations, as well as their differences from FIRE-1 simulations. Other FIRE-2 simulations have already been presented and studied in recent papers (e.g. Wetzel et al., 2016; Fitts et al., 2017; Garrison-Kimmel et al., 2017; El-Badry et al., 2018). We describe the initial conditions of our simulated sample in Section 2.1 and review the baryonic physics adopted in FIRE-2 briefly in Section 2.2.

2.1 Initial conditions

We run a set of dark matter-only cosmological boxes at low resolution to $z=5$, select target halos from the $z=5$ snapshots, and re-simulate these halos and the regions around them at much higher resolution with baryons using the well developed multi-scale cosmological ‘zoom-in’ techniques (Katz & White, 1993; Oñorbe et al., 2014). The initial conditions of the parent boxes and the zoom-in simulations are generated at $z=99$ using the music code (Hahn & Abel, 2011), with Planck 2015 cosmological parameters.

We use three dark matter-only cosmological boxes of side-length 11, 22, and 43 comoving Mpc, respectively. We use the spherical overdensity-based Amiga Halo Finder (ahf; Knollmann & Knebe, 2009) to identify halos in the $z=5$ snapshots, applying the redshift-dependent virial parameter from Bryan & Norman (1998), which leads to a virial overdensity $\Delta_{\rm vir}\approx 177$ (relative to background) for the redshift range we consider in this paper. We also checked the results against the six-dimensional phase-space halo finder rockstar (Behroozi et al., 2013a) and found good agreements in halo mass functions. We randomly select target halos in the $z=5$ halo mass range ${\rm M}_{\rm halo}=2\times 10^{9}$–$10^{12}\,{\rm M}_{\sun}$, requiring that there is no more massive halo within $5R_{\rm vir}$ from the target halo. This selection excludes 1/3 of the halos in the box333We also include other well-resolved halos in the zoom-in regions in our analysis (see Section 2.3). These halos live in the vicinity of a more massive halo (the target halo in the zoom-in region) by design. This will partially compensate the selection bias due to the isolation criteria above..

We identify zoom-in regions based on particles within $\sim 3$–$5R_{\rm vir}$ of the targeted halo, and iterate to ensure zero mass contamination from low-resolution particles within $2R_{\rm vir}$ and less than 1% contamination within $3R_{\rm vir}$ at $z=5$. There may be more than one halo in the zoom-in region, but the target halo is the most massive one by design. In Table 1, we list all of our target halos studied in this paper, along with the halo mass and stellar mass of the central galaxy (see Section 2.3 for details) at $z=5$ and $z=10$, and initial particle masses and minimum Plummer-equivalent force softening lengths of baryonic and high-resolution dark matter particles.

2.2 Baryonic physics

We briefly review the baryonic physics here, but refer to Hopkins et al. (2017, sections 2.3–2.5 and appendix B–E) for details. In the simulations, gas follows an ionized-atomic-molecular cooling curve from $10$–$10^{10}$ K, including metallicity-dependent fine-structure and molecular cooling at low temperatures and high-temperature metal-line cooling for 11 separately tracked species (H, He, C, N, O, Ne, Mg, Si, S, Ca, and Fe; see Wiersma et al., 2009). We do not include a primordial chemistry network nor Pop III star formation, but apply a metallicity floor of $Z=10^{-4}\,Z_{\sun}$, which corresponds crudely to the metallicity expected after enrichment by the first supernovae (SNe) from Pop III stars (e.g. Bromm et al., 2003; Wise et al., 2012). This is a reasonable treatment since we mainly focus on relatively massive galaxies at $z\lesssim 15$, which are dominated by Pop II stars. At each timestep, the ionization states and cooling rates H and He are calculated following Katz et al. (1996), and cooling rates from heavier elements are computed from a compilation of cloudy runs (Ferland et al., 2013), applying a uniform but redshift-dependent photo-ionizing background from Faucher-Giguère et al. (2009)444The ionizing background starts at $z=10.6$, with the ionization rate and heating rate increasing with time until the simulations end at $z=5$. We note that both rates show a sharp increase just below $z\sim 7$. A tabulated version of the background is available at http://galaxies.northwestern.edu/uvb/., and an approximate model for H ii regions generated by local sources. Gas self-shielding is accounted for with a local Jeans-length approximation, which is consistent with the radiative transfer calculations in Faucher-Giguère et al. (2010). The on-the-fly calculation of ionization states is broadly consistent with more accurate post-processing radiative transfer calculations (Ma et al., 2015).

We follow the star formation criteria in Hopkins et al. (2013) and allow star formation to take place only in dense, molecular, and locally self-gravitating regions with hydrogen number density above a threshold $n_{\rm th}=1000\,{\rm cm}^{-3}$. Stars form at 100% efficiency per free-fall time when the gas meets these criteria, and there is no star formation elsewhere. Note that star-forming particles can reach densities much higher than $n_{\rm th}$ following the self-gravitating criterion. The simulations include several different stellar feedback mechanisms, including (1) local and long-range momentum flux from radiative pressure, (2) energy, momentum, mass and metal injection from SNe and stellar winds, and (3) photo-ionization and photo-electric heating. Every star particle is treated as a single stellar population with known mass, age, and metallicity, assuming a Kroupa (2002) IMF from 0.1–$100\,{\rm M}_{\sun}$. All feedback quantities are directly calculated from starburst99 (Leitherer et al., 1999).

Note that starburst99 is a single-star stellar population model555Note that the stellar population models used in the simulations do not including stellar rotation, which is another key ingredient in stellar population synthesis (e.g. Leitherer et al., 2014; Choi et al., 2017) and could have important implications for reionization (e.g. Topping & Shull, 2015)., which assumes each star evolves independently, but most massive stars are expected to interact with a companion during their lifetimes. This will have significant effects on the SED of young populations, especially at low metallicities (e.g. Stanway, 2017). It has been suggested that massive binaries can lead to high escape fractions of ionizing photons from high-redshift metal-poor galaxies, and thus have important implications for understanding the sources dominating cosmic reionization (Ma et al., 2016b; Götberg et al., 2017). Nonetheless, binarity only has weak effects on most stellar feedback quantities, such as bolometric luminosities (within 0.05 dex in the first 200 Myr since a stellar population is born) and Type-II SNe rates (e.g. Xiao & Eldridge, 2015), so we do not expect binary interaction to have significant dynamical effects666Binary models do produce more ionizing photons (see Section 6.4), which are likely to enhance photo-ionization feedback, but we checked that this only has sub-dominant effects on gas dynamics (Ma et al., 2015).. For these reasons, we only consider binary stellar population models in post-processing. In this paper, we use the Binary Population and Spectral Synthesis (BPASS) models (version 2.0; Eldridge et al., 2008; Stanway et al., 2016)777http://bpass.auckland.ac.nz to compute the SED of each star particle from its age and metallicity. The BPASS models include both single-stellar and binary stellar population synthesis models. Their single-star models agree well with starburst99. Their binary models take into account mass transfer, common envelope phase, binary mergers, and quasi-homogeneous evolution at low metallicities. Also, the BPASS binary models appear to explain the nebular emission line properties observed in $z\sim 2$–3 galaxies (e.g. Steidel et al., 2016; Strom et al., 2017). In this paper, we mainly consider stellar continuum emission, while detailed modeling of dust extinction and nebular line emission will be the subject of future studies.

2.3 Halo selection and definitions

We run ahf on every snapshot to identify halos and subhalos in the zoom-in region. In general, most stars of the central (satellite) galaxy in a halo (subhalo) reside in $\frac{1}{3}R_{\rm max}$ from the halo center, where $R_{\rm max}$ is the radius at which the maximum circular velocity $V_{\rm max}$ is reached ($R_{\rm max}$ is already computed by ahf). We thus define a galaxy by including all star particles within $\frac{1}{3}R_{\rm max}$ after removing contributions from subhalos outside $\frac{1}{5}R_{\rm max}$. This excludes star particles at large distances from the halo center (corresponding to diffuse stellar distributions) and allow us to mask satellites. For each galaxy, we obtain a list of star particles with known position, age, and metallicity, from which we can compute a number of galaxy properties, such as stellar mass, star formation history, broad-band luminosities and magnitudes, surface brightness, galaxy size, etc. In this paper, we primarily focus on central galaxies, which dominate the stellar light: at a given stellar mass, only a few per cent of the galaxy population are satellites. We have also confirmed that they do not differ significantly from centrals at similar stellar masses in most properties we study in this paper. We restrict our analysis below to central halos that have zero contamination from low-resolution particles within $R_{\rm vir}$ and contain more than $10^{4}$ particles in total888This excludes most halos below ${\rm M}_{\rm halo}\sim 10^{8.6}$ ($10^{7.7}$) ${\rm M}_{\sun}$ in simulations at resolution $m_{b}\sim 7000$ (900) ${\rm M}_{\sun}$. The minimum number of dark matter particles for halos in our simulated catalog is $\sim 5600$.. Our target halos are guaranteed to meet this criteria by construction, but we also consider other halos in the zoom-in regions in our analysis.

In Figure 1, we show the number of halos that meet the above criteria in all zoom-in regions at several redshifts. Our simulations are able to capture (at least some) halo-to-halo variance below ${\rm M}_{\rm halo}=10^{11}\,{\rm M}_{\sun}$ at $z\sim 6$ and below ${\rm M}_{\rm halo}=10^{10.5}\,{\rm M}_{\sun}$ at $z\sim 10$, where the simulations include a few to more than 200 halos in a given halo mass bin. Moreover, these galaxies always have ‘bursty’ star formation histories (see Section 3.3), which leads to significant time variability in their properties (e.g. Muratov et al., 2015; Sparre et al., 2017; Feldmann et al., 2017; Ma et al., 2017b; Faucher-Giguère, 2018). Hence a galaxy tends to move above and below the median of certain scaling relations (see also the discussion in Section 3.2). To account for the scatter due to bursty star formation (as well as galaxy mergers and other time-variable phenomena), we make use of 58 snapshots saved for each simulation at redshifts $z=5$–12 (about 20 Myr apart between two successive snapshots) to build a catalog of over 34,000 simulated halo ‘snapshots’. By doing so, we sample the same halos multiple times in the catalog and treat them as statistically equal in our analysis. Figure 1 essentially shows the number of independent halos in the simulated sample. At lower masses (e.g. ${\rm M}_{\rm halo}\leq 10^{11}\,{\rm M}_{\sun}$ at $z\sim 6$), we are able to account for the scatter both from halo-to-halo variance and time variability within single halos. A priori, it is not clear which effect dominates the scatter for a given scaling relation. At higher masses (${\rm M}_{\rm halo}\geq 10^{11}\,{\rm M}_{\sun}$), our sample only contains 1–2 independent halos at a given redshift, so we are only able to account for the variance due to time variability of individual galaxies. We caution that we may therefore underestimate the scatter of certain scaling relations at the high-mass end. We have also checked that excluding a randomly selected 1/2–2/3 of the snapshots from our analysis (sampling each galaxy at sparser time steps) does not change the results of this paper. In other words, our time-sampling is sufficient for statistically converged results.

2.4 Halo abundances

Since our simulated catalog is constructed from 15 cosmological zoom-in regions, it does not contain information about the halo abundance at a given halo mass and redshift. Therefore, we assign every simulated halo ‘snapshot’ a weight to recover the appropriate number density of halos at its mass and redshift in the Universe. We briefly summarize the method here and refer the readers to Appendix A for details. We use HMFcalc (Murray et al., 2013)999http://hmf.icrar.org to calculate the halo mass functions, applying the same cosmological parameters and virial overdensities as those adopted in the simulations. We take the fitting functions from Behroozi et al. (2013b) in HMFcalc, which is a modified Tinker et al. (2008) halo mass function. It matches well with the halo mass functions directly extracted from our large-volume dark matter-only cosmological boxes in the redshift range we consider here. We bin the simulated sample in the two-dimensional $\log{\rm M}_{\rm halo}$–$\log(1+z)$ space with bin widths $\Delta\log{\rm M}_{\rm halo}=0.4$ from ${\rm M}_{\rm halo}=10^{7.5}$–$10^{12}\,{\rm M}_{\sun}$ and $\Delta\log(1+z)=0.04$ from $z=5$–12. We have confirmed that our results are not sensitive to the bin widths we adopt. In each bin, we count the number of halos in the simulated catalog $N_{\rm sim}$ and calculate the number of halos expected in the Universe $N_{\rm expect}$ from the halo mass function. All halos in the same bin are assigned the same weight $w=N_{\rm expect}/N_{\rm sim}$. In other words, by summing $w$ over all simulated halos in certain halo mass and redshift intervals and dividing $\sum_{i}w_{i}$ by the corresponding comoving volume, one should recover the halo number densities given by the input halo mass functions. When necessary, each halo snapshot in the simulated catalog is weighted by its $w$. This is important when we consider statistical properties of simulated galaxies at a fixed stellar mass or magnitude, where not all galaxies have equal halo mass (e.g. Sections 3.3 and 3.4). The weights will also be used to construct stellar mass functions and luminosity functions in Sections 4 and 5.

Strictly speaking, this approach is valid only if the halos in our samples are not strongly biased. However, after an extensive check, we find no significant difference between halos in various environments (see also e.g. O’Shea et al., 2015) and at different resolution in our simulations regarding their properties studied in this paper. Although our sample is still possibly biased due to complex selection criteria – for example, all halos below $10^{9}\,{\rm M}_{\sun}$ by $z=5$ in our sample live within a few virial radii of a more massive halo (i.e., the target halo in the zoom-in region) and we lack isolated halos at such low masses down to $z=5$, our conclusions in this paper are likely robust.

3 Galaxies in the reionization era

3.1 Morphology

In Figure 2, we show the stellar $u/g/r$-composite images at $z=5$ (left) for the central galaxy in the most massive halo in each zoom-in region. The stellar masses and halo masses are listed in Table 1. We use the BPASSv2.0 binary models to determine the stellar SEDs, assuming a Kroupa (2002) IMF from 0.1–$100\,{\rm M}_{\sun}$. Note that we only consider intrinsic stellar continuum emission, and ignore dust extinction and nebular line emission at this point101010Full spectral modeling of high-redshift galaxies in cosmological simulations has been developed recently by other groups (e.g. Wilkins et al., 2016; Barrow et al., 2017; Zackrisson et al., 2017).. The right panel shows the noise-free mock images as observed by the Near Infrared Camera (NIRCam) on JWST at F277W band (rest-frame 4600 Å), applying a Gaussian point spread function (PSF) with full width half maximum (FWHM) of two pixels with pixel size 0.065 arcsec (0.42 kpc in physical length)111111The PSF and pixel sizes are adopted from the NIRCam pocket guide from https://jwst.stsci.edu/instrumentation/nircam.. The three images in the same row are shown using the same color scale, which spans eight magnitudes in surface brightness, but the depth increases from $m_{\rm AB}=29.5$ mag arcsec*-2* in the top row to 31.5 mag arcsec*-2* in the bottom row (pixels below these limits are shown as black).

In Figure 3, we show the stellar $u/g/r$-composite images and noise-free mock JWST NIRCam F444W-band images (rest-frame 4000 Å) at $z=10$ for the most massive galaxy in each zoom-in simulation. These images are rearranged in place to ensure a descending order in halo mass from the top-left panel to the bottom-right panel. The mock JWST images have a pixel size 0.065 arcsec and 0.28 kpc in physical length. Again, the color scale in each image spans eight magnitudes in surface brightness, but the depth increases from $m_{\rm AB}=32$ mag arcsec*-2* in the top row to 36 mag arcsec*-2* in the bottom row.

Almost all of the simulated galaxies at $z\geq 5$ show clumpy, irregular morphologies even in rest-frame optical bands, possibly due to high merger rates and clumpy, gas-rich star formation at these redshifts. This is in contrast to galaxies at low and intermediate redshifts, which show a mix of late-type, early-type, and irregular morphologies at similar masses (e.g. Feldmann et al., 2017; El-Badry et al., 2018). Galaxies with similar stellar mass may have a variety of sizes and surface brightness, so their detectability can differ significantly. Therefore, our high-resolution simulations provide a useful database for understanding future multi-band, spatially resolved observations of $z\gtrsim 5$ galaxies, as well as determining the completeness of a flux-limited galaxy survey at these redshifts.

3.2 The stellar mass–halo mass relation

Figure 4 shows the stellar mass–halo mass relation (top panels) and the stellar baryon fraction–halo mass relation (bottom panels) for central galaxies at $z=6$, 8, and 10. The stellar baryon fraction is defined as ${\rm M}_{\ast}/(f_{b}{\rm M}_{\rm halo})$, where $f_{b}=\Omega_{b}/\Omega_{m}$ is the cosmic baryonic fraction. The two-dimensional histograms represent the number of halo snapshots in the simulated catalog (as defined in Section 2.3) within $\Delta z=0.5$ (e.g. from $z=7.5$–8.5 in the $z=8$ panels) in each pixel in logarithmic scale (as shown by the color scale). We remind the readers that we re-sample each halo multiple times to account for time-variability in where galaxies lie on this relation, but we refer to Figure 1 for the number of independent halos in our sample at these redshifts (see Section 2.3 for details). We also show the empirical relations at $z=6$ and 8 from Behroozi et al. (2013b, the cyan lines). At all redshifts, there is a tight correlation between stellar mass and halo mass, with the scatter increasing at the low-mass end121212We caution that our simulated sample have more independent galaxies at low masses than at higher masses, so we may underestimate the scatter at the high-mass end. Nonetheless, a halo mass-dependent scatter in the stellar mass–halo mass relation does exist at low redshift in both observations and FIRE-2 simulations (e.g. Garrison-Kimmel et al., 2014; Fitts et al., 2017). The current simulations are consistent with increased scatter at low masses.. We also examine the relation for satellite galaxies (not shown), which tend to have systematically higher stellar mass and larger scatter than central galaxies at a given halo mass, due to the fact that their halos are usually stripped. However, we note that the halo mass of a satellite depends strongly on which halo finder one uses. We find a smaller offset between central and satellite galaxies using the rockstar subhalo catalog than using the ahf catalog. Because satellite galaxies contribute no more than a few per cent of the total galaxy population at a given mass, we do not further quantify the difference in this paper.

We find little evolution in the stellar mass–halo mass relation at these redshifts, in line with recent empirical constraints (e.g. Rodríguez-Puebla et al. 2017; however, see Behroozi & Silk 2015). We will show it explicitly below. Using all halo snapshots in the simulated catalog at $z=5$–12, we calculate the median and $1\sigma$ dispersion in $\log{\rm M}_{\ast}$ at every 0.5 dex in $\log{\rm M}_{\rm halo}$ from ${\rm M}_{\rm halo}=10^{7.5}$–$10^{12}\,{\rm M}_{\sun}$. We assume a simple power-law relation between ${\rm M}_{\ast}$ and ${\rm M}_{\rm halo}$

[TABLE]

and a halo mass-dependent scatter

[TABLE]

and fit the median stellar mass–halo mass relation and $1\sigma$ dispersion obtained from the simulated sample as described above. We obtain the best-fit parameters as

[TABLE]

These results are nearly identical to those obtained from the FIRE-1 simulations (e.g. Hopkins et al., 2014; Ma et al., 2015; Feldmann et al., 2017) at similar halo mass and redshift, despite the subtle differences in numerical details and resolution between these simulations. Our predictions also broadly agree with those in the literature (e.g. Ceverino et al., 2017). We show our best-fit $1\sigma$ stellar mass–halo mass relation in every panel in Figure 4 (the red dashed lines). Visual inspection implies that equations 1 and 2 describe our simulated sample reasonably well at any redshift. We also confirmed that the median relation obtained from a subsample at a given redshift does not deviate from Equation 1 by more than 0.1 dex at most halo masses we consider here. It is an intriguing question why the ${\rm M}_{\ast}$–${\rm M}_{\rm halo}$ relation does not evolve at these redshifts. We speculate that this is probably due to feedback regulating star formation to zeroth order. A detailed analysis of the relation between halo growth rate, gas accretion rate, and stellar mass growth rate to understand the weak-evolution of the ${\rm M}_{\ast}$–${\rm M}_{\rm halo}$ relation is worth future investigation.

How do galaxies evolve on the ${\rm M}_{\ast}$–${\rm M}_{\rm halo}$ relation? All of our simulated galaxies experience bursty star formation because of stellar feedback. The stellar mass can grow by a factor of 2 or more during a short time period at the peak of a starburst, while it can remain almost unchanged during the troughs of its star formation history (see Section 3.3 and Figure 5 for examples). In contrast, the halo mass grows relatively smoothly via dark matter accretion, which is less affected by feedback. As a consequence, a galaxy moves vertically on the ${\rm M}_{\ast}$–${\rm M}_{\rm halo}$ plane during the peak of a starburst and reaches some point above the median ${\rm M}_{\ast}$–${\rm M}_{\rm halo}$ relation, while it then moves horizontally during a trough in its star formation history and reaches somewhere below the median relation until the next starburst episode. We confirmed in our simulations that the scatter in the ${\rm M}_{\ast}$–${\rm M}_{\rm halo}$ relation caused by bursty star formation is a physical effect.

There are several caveats in this analysis. First, at the high-mass end (${\rm M}_{\rm halo}\geq 10^{11}\,{\rm M}_{\sun}$), our approach only captures the scatter due to bursty star formation, but the sample does not contain sufficient numbers of independent halos to account for halo-to-halo variance. Therefore, we may underestimate the scatter at these halo masses. At lower masses, our simulations include considerably larger numbers of independent halos and the scatter is reliably measured. Second, our simulations do not include more massive halos above ${\rm M}_{\rm halo}=10^{12}\,{\rm M}_{\sun}$. At these masses, slowly-cooling hot halos and feedback from supermassive black holes may play an important role. Studying early galaxy formation in such high-redshift massive halos is beyond the scope of this paper. It may lead to a turnover in the ${\rm M}_{\ast}$–${\rm M}_{\rm halo}$ relation at these redshifts similar to what is seen at lower redshifts (e.g. Behroozi et al., 2013b). We caution that our best-fit ${\rm M}_{\ast}$–${\rm M}_{\rm halo}$ relation may break down at ${\rm M}_{\rm halo}>10^{12}\,{\rm M}_{\sun}$. Lastly, the best-fit ${\rm M}_{\ast}$–${\rm M}_{\rm halo}$ relation does not apply to halos with ${\rm M}_{\rm halo}\lesssim 10^{8}\,{\rm M}_{\sun}$ below $z\sim 6$. We will show in Section 3.3 that star formation in these low-mass halos is suppressed by the ionizing background near the end of reionization.

3.3 Star formation histories

Figure 5 shows the star formation histories (top panels) and the stellar mass growth histories (bottom panels) for four example galaxies in the $z=5$ halo mass range ${\rm M}_{\rm halo}=10^{8}$–$4\times 10^{11}\,{\rm M}_{\sun}$ (as labeled in each panel). The black solid lines and red dashed lines in the top panels show the star formation rates (SFRs) averaged over 10 Myr and 100 Myr, respectively. These are proxies for the H$\alpha$- and UV-inferred SFRs observationally (e.g. Sparre et al., 2017). All simulated galaxies show significant ‘bursty’ star formation histories, with starbursts occurring on time-scales of 50–100 Myr. This feature is also seen in other cosmological zoom-in simulations with comparably high resolution and detailed physics despite different numerical methods (e.g. Kimm & Cen, 2014; Ceverino et al., 2018). As discussed in Section 3.2, the stellar mass can grow almost instantaneously by a factor of 2 or more at the peak of a burst, while it remains nearly constant when the SFR is low. In the least massive halo (${\rm M}_{\rm halo}=10^{8}\,{\rm M}_{\sun}$ at $z=5$), the feedback from a starburst at $z\sim 8$ expels most of its gas. At later times, gas accretion and cooling becomes inefficient as heating from the ionizing background becomes significant for such low mass halos (e.g. Efstathiou, 1992; Thoul & Weinberg, 1996; Quinn et al., 1996; Okamoto et al., 2008; Faucher-Giguère et al., 2011; Noh & McQuinn, 2014; Sawala et al., 2016). Star formation in these halos is thus suppressed at later times. More massive halos are able to maintain star formation until the end of our simulation at $z=5$.

In Figure 6, we show the SFR (left panel) and specific SFR (sSFR, right panel) averaged over the past 100 Myr as a function of stellar mass for all galaxy snapshots in the simulated sample, color coded by their redshift. Note that the sharp upper limits at a given stellar mass in both relations are due to the fact that some galaxies form essentially all of their stellar mass in a starburst during the past 100 Myr. We also compare our results with the observed relation for $z=6$–9 galaxies from McLure et al. (2011) (black dashed lines). Our simulations agree well with observations at ${\rm M}_{\ast}\geq 10^{8}\,{\rm M}_{\sun}$. At lower masses (where there are no observations), the scatter is larger due to stronger burstiness in their star formation histories, as illustrated by the solid lines in the right panel. Moreover, at fixed stellar mass, galaxies at lower redshifts have lower star formation rates on average than galaxies at higher redshifts. This trend is expected because the stellar mass growth time-scale (the ratio of stellar mass to star formation rate) of galaxies at a given redshift should be comparable to the Hubble time at that redshift and has also been found in previous studies (e.g. Behroozi & Silk, 2015; Wilkins et al., 2017).

We now derive the weighted average SFR as a function of halo mass (or stellar mass) and redshift at every $\Delta\log{\rm M}_{\rm halo}=0.5$ (or $\Delta\log{\rm M}_{\ast}=0.5$) and $\Delta\log(1+z)=0.04$ for the simulated sample. All halo snapshots are included in this calculation. We then fit the results using two-dimensional linear functions

[TABLE]

and

[TABLE]

and obtain the best-fit parameters $(\alpha,\,\gamma,\,\delta)=(1.58,\,2.20,\,-1.58)$ and $(\alpha^{\prime},\,\gamma^{\prime},\,\delta^{\prime})=(1.03,\,2.01,\,1.36)$. In Figure 7, we show the ${\rm SFR}({\rm M}_{\rm halo},\,z)$ (left panel) and ${\rm SFR}({\rm M}_{\ast},\,z)$ (right panel) relations derived from the simulated sample. The colors represent the weighted average SFR relative to the $z=5$ best-fit relation (in logarithmic scale, as labeled at the colorbars). This eliminates the wide dynamic range shown in the left panel of Figure 6 and allows us to see the evolution more clearly. Note that the white regions show the parameter space where our simulated sample contains no galaxies. At fixed halo mass or stellar mass, the average SFR increases by $\sim 0.7\,{\rm dex}$ from $z=5$–12, in broad agreement with the qualitative trend shown in Figure 6.

Furthermore, after $z\sim 6$, the average SFRs in low-mass galaxies (halo mass below ${\rm M}_{\rm halo}\sim 10^{8}\,{\rm M}_{\sun}$ or stellar mass below ${\rm M}_{\ast}\sim 10^{4}\,{\rm M}_{\sun}$, the dark blue region in Figure 7) are significantly lower (by about 0.6 dex) than those inferred from the fitting functions. This is because, as mentioned above, at the end of the reionization era, the ionizing background can heat the gas in low-mass halos efficiently and prevent it from cooling and forming stars. The star formation in these galaxies is thus suppressed. We note that halos of similar masses at higher redshifts ($z\gtrsim 7$) or more massive halos (${\rm M}_{\rm halo}\gtrsim 10^{8.5}\,{\rm M}_{\sun}$) at any redshift are not affected and continue normal star formation.

In Figure 5, we show one example of such low-mass galaxies (the left most panel), where star formation is suppressed at lower redshifts. There are also halos at these masses which are completely ‘dark’ (containing no stars). The dark halo fraction is negligible for halos above ${\rm M}_{\rm halo}=10^{8.5}\,{\rm M}_{\sun}$ at any redshift, whereas at ${\rm M}_{\rm halo}\sim 10^{8}\,{\rm M}_{\sun}$, the dark fraction increases from less than 1% at $z=12$ to approximate 50% at $z=5$.131313Note that the increasing dark fraction at ${\rm M}_{\rm halo}\sim 10^{8}\,{\rm M}_{\sun}$ and below $z\sim 6$ indicates that the suppression of star formation in these halos is not purely due to stellar feedback but rather points to the importance of reionization. We will show later in Sections 4 and 5 that this effect leaves an imprint in the stellar mass function and luminosity functions at $z\sim 6$.

Our findings are broadly in line with other simulations in the literature. For example, Wise et al. (2014) found no dark halos above ${\rm M}_{\rm halo}\sim 10^{8}\,{\rm M}_{\sun}$ at $z>8$. Sawala et al. (2016) found that the dark fraction decreases sharply from nearly 90% at ${\rm M}_{\rm halo}\sim 10^{8}\,{\rm M}_{\sun}$ to 0% at ${\rm M}_{\rm halo}\sim 10^{8.5}\,{\rm M}_{\sun}$ at $z\sim 10$. They also find an increasing dark fraction with decreasing redshift at a fixed halo mass. The subtle differences are likely due to different models of the ionizing background adopted in these studies, as well as to different star formation and stellar feedback physics. Wise et al. (2014) modeled the ionizing fields more self-consistently using radiative-hydrodynamic methods, while Sawala et al. (2016) adopted the uniform Haardt & Madau (2001) ionizing background at these redshifts. In addition, a dark halo only means that the expected stellar mass is lower than the mass of a few star particles. This further complicates the comparison between these results obtained at different resolutions. The effects of the ionizing radiation fields prior to complete reionization on low-mass galaxies merits future investigation.

3.4 Broad-band photometry

We use the BPASSv2.0 stellar population synthesis models to calculate the broad-band luminosities and magnitudes for the simulated galaxies, using the binary models with a Kroupa (2002) IMF from 0.1–$100\,{\rm M}_{\sun}$ as our default model. We only consider intrinsic stellar continuum here, and ignore dust extinction and strong nebular line emission in the rest-frame UV and optical, as well as dust re-emission in the infrared (IR). We will explore the effect of dust attenuation in Section 6.2. In Figure 8, we show the magnitude–stellar mass relation (top panels) for our simulated sample with ${\rm M}_{\ast}>10^{3.5}\,{\rm M}_{\sun}$ for three example redshift and band combinations. We also show the inverse relation, the stellar mass–magnitude relations for the same combinations in the bottom panels. Only galaxies brighter than ${\rm M_{AB}}=-10$ are shown to ensure that our simulations are complete. The two-dimensional histogram represents the total weight (as defined in Section 2.4) of all galaxy snapshots in a pixel in logarithmic scale. This reflects the ‘correct’ relative number of galaxies in the Universe between pixels141414We remind the reader that we include all snapshots in the analysis to account for time variability of galaxy properties, which is important for UV luminosities, but we may underestimate the scatter for halos above ${\rm M}_{\rm halo}>10^{11}\,{\rm M}_{\sun}$ where our sample does not contain large number of independent halos (see Section 2.3 for details).. The red solid and dashed lines illustrate the best-fit median relation and $1\sigma$-scatter (16–84 per cent) as obtained below.

At fixed stellar mass, the distribution of magnitudes in a specific band tends to be asymmetric, with a broader spread at the bright end. The asymmetry is driven by the evolution of stellar populations: the luminosity of a stellar population declines rapidly as the most massive stars die (in about 3–30 Myr). Therefore, the luminosity of a galaxy depends not only on its total stellar mass but also on its recent star formation history. This feature is more prominent in low-mass galaxies which have significant bursty star formation histories. Figure 8 also shows that this effect is strongest in the rest-frame UV where young stars overwhelmingly dominate the starlight and becomes weaker at longer wavelengths, as rest-frame optical B-band relation has smaller scatter than that of rest-frame 1500 Å.

Furthermore, galaxies at higher redshifts appear brighter on average than those of similar stellar masses at lower redshifts, simply due to the fact that high-redshift galaxies have younger stellar populations and higher ongoing SFRs. We parametrize the magnitude–stellar mass relation with a linear function

[TABLE]

where we assume a fixed slope $a$ at any redshift but a redshift-dependent normalization to capture this feature. We fit the weighed median, $1\sigma$ lower- and upper-bound relations (above ${\rm M}_{\ast}=10^{3.5}\,{\rm M}_{\sun}$) obtained from eight subsamples in different redshift intervals from $z=5$–12 all together to determine the parameters for a given band. The top panels of Figure 8 illustrate three examples of this relationship, and we list the best-fit parameters for rest-frame 1500 Å, B band, and J band in the top half of Table 2.

We similarly assume a linear function for the stellar mass–magnitude relation

[TABLE]

and fit the weighted median and $1\sigma$ relations for galaxies brighter than ${\rm M_{AB,\,band}}=-10$ to obtain the parameters. Some examples are shown in Figure 8 and the best-fit parameters for rest-frame 1500 Å, B band, and J band are listed in Table 2 (the bottom block). We emphasize that the two relations are fundamentally different from each other – the distribution of stellar mass at fixed magnitude is biased toward low-mass galaxies, simply because they have much higher number densities in the Universe than more massive galaxies (see also Section 4). The stellar mass–magnitude relation is also redshift-dependent, with the median stellar mass decreasing by about 1 dex from $z=12$ to 5 at a given magnitude.

4 Stellar mass functions in the early universe

Now we calculate the stellar mass function using two distinct approaches. First, we utilize the weights constructed in Section 2.4: at a certain redshift, we collect all halo snapshots within $\Delta z=\pm 0.5$ from our simulated catalog. We then add the weights of galaxies in stellar mass bins and divide $\sum_{i}w_{i}$ by the comoving volume corresponding to the $\Delta z=\pm 0.5$ redshift interval. Only halos above ${\rm M}_{\rm halo}=10^{7.5}\,{\rm M}_{\sun}$ are taken into account. In Figure 9, we show the results in number $\rm dex^{-1}\,Mpc^{-3}$ above ${\rm M}_{\ast}=10^{3.5}\,{\rm M}_{\sun}$ at $z=6$, 8, 10, and 12 with the open symbols. The data are tabulated in Tables 3 and 4 for readers to use.

Alternatively, we can model the stellar mass function by directly convolving the stellar mass–halo mass relation derived in Section 3.2 with the halo mass function at a given redshift. We use a Monte Carlo method: we generate a large number of mock halos more massive than ${\rm M}_{\rm halo}=10^{7.5}\,{\rm M}_{\sun}$ with number densities following the halo mass function, assign each halo a stellar mass as described below, and derive the stellar mass function from the mock catalog. We assume a) every halo hosts one galaxy (considering only central galaxies) and b) the stellar mass follows a lognormal distribution at a given halo mass, with the median and $1\sigma$ dispersion following Equations 1 and 2. In this calculation, we use Equations 1 and 2 at all redshifts and all halo masses, but we caution that uncertainties may arise at the high-mass end (see discussion in Section 3.2). The results are shown by dashed lines in Figure 9. For simplicity, we refer them as model stellar mass functions thereafter. We also make these results publicly available (see Appendix C for details). We will compare our predictions with observations later in Section 6.3.

We highlight the following features shown in Figure 9. (1) The low-mass end of the stellar mass function (asymptotic form $\phi\,{\rm d}\log{\rm M}_{\ast}\sim{\rm M}_{\ast}^{\alpha+1}\,{\rm d}\log{\rm M}_{\ast}$) steepens with increasing redshift, with the slope decreasing from $\alpha=-1.80\pm 0.02$ at $z=6$ to $\alpha=-2.13\pm 0.12$ at $z=12$. The evolution of the slope is robust, although the exact slope at a given redshift may vary according to how it is computed151515The slopes quoted here are obtained by fitting the stellar mass functions derived from the simulated catalog with a Schechter (1976) function. We also experiment with fitting the model mass functions in different dynamic ranges or using a double power-law function. The slope obtained at a given redshift varies systematically with method by about 0.2.. This trend is consistent with the observed stellar mass functions (e.g. Song et al., 2016) and has been widely reproduced in cosmological simulations (e.g. Wilkins et al., 2016; Ceverino et al., 2017). Such a feature is directly inherited from the halo mass functions, which also steepen with increasing redshift at the low-mass end (e.g. Reed et al., 2003). (2) The model stellar mass functions agree well with those derived from the simulated catalog for a broad range of stellar mass and redshift. This demonstrates that the fitting functions in Section 3.2 describe the stellar mass–halo mass relation for the simulated sample reasonably well. (3) The discrepancies between the two stellar mass functions in the highest-mass bin is due to small numbers of galaxies in the simulated sample at the high-mass end. (4) The apparent flattening of the model stellar mass functions at ${\rm M}_{\ast}\sim 10^{3.5}\,{\rm M}_{\sun}$ is an expected artifact because we exclude all halos below ${\rm M}_{\rm halo}=10^{7.5}\,{\rm M}_{\sun}$.

More importantly, the $z=6$ stellar mass function derived from the simulated sample shows a flattening below ${\rm M}_{\ast}\sim 10^{4.5}\,{\rm M}_{\sun}$ and falls below the model mass function by a factor of 2. This is caused by the 50% fraction of dark halos at ${\rm M}_{\rm halo}\sim 10^{8}\,{\rm M}_{\sun}$ at $z\sim 6$. In other words, the assumption we adopted in the model that every halo hosts one galaxy breaks down at ${\rm M}_{\rm halo}\sim 10^{8}\,{\rm M}_{\sun}$. Note that if we ignore all the ‘dark halos’ in the simulated catalog and repeat the exercise, the two $z=6$ stellar mass functions agree well with each other. The large dark fraction in low-mass halos at lower redshifts is because of the suppression of star formation by the ionizing background near the end of reionization (see Section 3.3). The stellar mass functions at higher redshift do not show such flattening. This effect may relieve the tension between the number of low-mass galaxies in the Local Group and that needed for cosmic reionization (e.g. Bullock et al., 2000; Somerville, 2002; Boylan-Kolchin et al., 2014).

5 Multi-band luminosity functions

We calculate the luminosity functions at several bands following the same method described in Section 4. Again, only halos above ${\rm M}_{\rm halo}=10^{7.5}\,{\rm M}_{\sun}$ are taken into account. In Figure 10, we show the results in number $\rm mag^{-1}\,Mpc^{-3}$ at $z=6$, 8, 10, and 12 (the open symbols). The data are also provided in Appendix C. We also model the luminosity functions using a Monte Carlo method by convolving the magnitude–stellar mass relation derived in Section 3.4 (Equation 6) with the model stellar mass functions in Section 4 in a similar way described above. The model luminosity functions are shown by the dashed lines in Figure 10. Note that we only consider intrinsic stellar continuum emission here, but we will explore the effect of dust extinction in Section 6.2 and compare with observations in Section 6.3. We do not model nebular line emission in this paper.

The luminosity functions derived from the simulations agree well with models for a broad range of magnitude and redshift. Again, the faint-end slope (asymptotic form $\phi\,{\rm d}{\rm M_{AB,\,band}}\sim 10^{-0.4(\alpha+1){\rm M_{AB,\,band}}}{\rm d}{\rm M_{AB,\,band}}$) steepens with increasing redshift (e.g. from $\alpha=-1.85\pm 0.06$ at $z=6$ to $\alpha=-2.17\pm 0.10$ at $z=12$ at 1500 Å). The trend is in good agreement with observations (e.g. Bouwens et al., 2015; Finkelstein et al., 2015a), semi-analytic galaxy formation models (e.g. Clay et al., 2015; Cowley et al., 2018), and other simulations (e.g. Gnedin, 2016; Ceverino et al., 2017; Wilkins et al., 2017) at these redshifts, but the exact slopes depend largely on the fitting method. The flattening at the faintest bin at any redshift is due to the incompleteness of halos below ${\rm M}_{\ast}=10^{7.5}\,{\rm M}_{\sun}$. Similarly, the luminosity functions show a flattening below $\rm M_{1500\,\text{\AA}}\sim-12$, $\rm M_{B}\sim-12$, and $\rm M_{J}\sim-12$ at $z=6$, as seen from the fact that the luminosity functions derived from the simulated catalog fall below the model luminosity functions roughly by a factor of 2. This is caused by the large fraction of dark halos and the rapid drop in SFR at $z\sim 6$ below halo mass ${\rm M}_{\rm halo}\sim 10^{8}\,{\rm M}_{\sun}$, where star formation is suppressed by the ionizing background (Figure 7).

6 Discussion

6.1 SFR and stellar mass densities at $z\geq 5$

We now derive the cosmic SFR density (SFRD) from $z=5$–12 by convolving the SFR–${\rm M}_{\rm halo}$ relation given by Equation 4 (the average SFR at a given halo mass) and the halo mass function at the same redshift. In Figure 11, we show the results obtained by integrating over the halo mass range above ${\rm M}_{\rm halo}=10^{7.5}\,{\rm M}_{\sun}$ (dashed), ${\rm M}_{\rm halo}=10^{8.5}\,{\rm M}_{\sun}$ (solid), and ${\rm M}_{\rm halo}=10^{9.5}\,{\rm M}_{\sun}$ (dotted) to ${\rm M}_{\rm halo}=10^{12}\,{\rm M}_{\sun}$. The contributions from more massive halos are negligible at these redshifts (because they are extremely rare). In Figure 11, we also show observationally inferred SFRDs from Ellis et al. (2013), Oesch et al. (2013), Oesch et al. (2014), Bouwens et al. (2015), Finkelstein et al. (2015a), McLeod et al. (2016), and CLASH detections (Zheng et al., 2012; Coe et al., 2013; Bouwens et al., 2014). Data corrected (uncorrected) for dust attenuation are shown by open (filled) symbols161616We note that the conversion between rest-frame UV luminosity and SFR and the amount of dust correction are still very uncertain.. At $z\lesssim 8$, our predictions broadly agree with data within observational uncertainties. The SFRD at $z\gtrsim 9$ are still poorly constrained observationally. These results are derived by integrating the best-fit UV luminosity functions brighter than $\rm M_{UV}\sim-17$. This limit does not correspond to a unique halo mass, but is roughly consistent with what we obtain by integrating down to ${\rm M}_{\rm halo}=10^{9.5}\,{\rm M}_{\sun}$ at these redshifts (cf. Figures 4 and 8). Note that the number of galaxies in the observed $z\gtrsim 9$ sample is small, and some works are based on single galaxy detections. Our simulations suggest that the majority of star formation takes place in halos below ${\rm M}_{\rm halo}=10^{9.5}\,{\rm M}_{\sun}$ at $z\gtrsim 9$, but these low-mass galaxies are too faint to be detectable with current observational facilities. This may account for the apparent rapid decline in SFRD at these redshifts (e.g. Oesch et al., 2014). Future deep surveys by JWST at these redshifts are expected to put strong constraints on the $z\geq 9$ SFRD.

We also calculate the stellar mass density from $z=5$–12 by integrating the model stellar mass functions in Section 4 in certain stellar mass intervals. The three lines in Figure 12 show the results for ${\rm M}_{\ast}>10^{6}\,{\rm M}_{\sun}$ (solid), ${\rm M}_{\ast}>10^{4}\,{\rm M}_{\sun}$ (dashed), and ${\rm M}_{\ast}>10^{8}\,{\rm M}_{\sun}$ (dotted). At these redshifts, high-mass galaxies (${\rm M}_{\ast}>10^{10}\,{\rm M}_{\sun}$) only contribute a negligible fraction (less than 0.05 dex) of the total stellar mass due to their low number densities, so the total stellar mass density is insensitive to our uncertainties in the high-mass end of the stellar mass functions. In Figure 12, we also compare our predictions with observationally inferred results in the literature (symbols with errorbars, including González et al., 2011; Duncan et al., 2014; Oesch et al., 2014; Grazian et al., 2015; Song et al., 2016). Note that these results are derived by integrating the best-fit stellar mass functions above ${\rm M}_{\ast}>10^{8}\,{\rm M}_{\sun}$. Our predictions (the dotted line, which uses the same limit) broadly agree with more recent studies (Oesch et al., 2014; Grazian et al., 2015; Song et al., 2016). We note that although some groups report consistent stellar mass densities at these redshifts, their stellar mass functions do not usually agree with each other (see Section 6.3 or figure 9 in Song et al. 2016).

6.2 Dust extinction in rest-frame UV

So far we have only focused on intrinsic luminosity of our simulated galaxies, while dust obscuration can be very important in relatively massive galaxies (e.g. Cen & Kimm, 2014; Cullen et al., 2017; Wilkins et al., 2017). In this section, we estimate the amount of dust extinction in the rest-frame UV band. We follow the method from Hopkins et al. (2005) and calculate the extinction by ray-tracing the emission from star particles including dust attenuation self-consistently from the dust and metals in the simulation (see also Hopkins et al., 2014; Feldmann et al., 2016, 2017). We assume a canonical dust-to-metal ratio of 0.4 (Dwek, 1998). Following a Small Magellanic Cloud-like extinction curve from Pei (1992), we obtain a dust opacity of $2.06\times 10^{3}\,{\rm cm^{2}\,g^{-1}}$ at 1500 Å at solar metallicity. For each simulated galaxy, we include all gas particles within $\frac{1}{2}R_{\rm max}$ (about 1.5 times the size of the stellar component, see Section 2.3) and calculate the extinction by ray-tracing from every star particle to a hypothetical observer along ten different sightlines. Note that the gas in these high-redshift galaxies is clumpy, so the extinction can differ by several magnitudes between sightlines. In Figure 13, we show the relation between extinction A${}_{1500}=-2.5\log(F_{1500}/F_{1500,\,0})$ and intrinsic UV magnitude M1500 for all sightlines and all simulated galaxies from $z=5$–8. We do not find significant redshift dependence in our sample, but we caution that this may be due to the small sample size at the massive end. We fit the results with a parabolic function

[TABLE]

and quote a uniformly distributed scatter with half-width $\Delta{\rm A_{1500}}=-0.125\,({\rm M_{1500}}+15)$ at ${\rm M_{1500}}<-15$. A full radiative transfer calculation of dust extinction, scatter, and re-emission is beyond the scope of the current paper, but will be the subject of a future study.

6.3 Comparison with observations

In this section, we compare our predicted stellar mass functions and luminosity functions with observations. In the top panel in Figure 14, we show the $z=6$ stellar mass function derived from the simulated catalog (open squares) and from direct convolution between the stellar mass–halo mass relation and the halo mass function (dashed lines). They are identical to those in Figure 9, but we only show ${\rm M}_{\ast}\geq 10^{7}\,{\rm M}_{\sun}$ where the observational results are available. We also show a compilation of observations from González et al. (2011), Duncan et al. (2014), Grazian et al. (2015), Song et al. (2016), and Stefanon et al. (2017) (symbols with errorbars). Above ${\rm M}_{\ast}\sim 10^{9}\,{\rm M}_{\sun}$, our model agrees well with Song et al. (2016) and Stefanon et al. (2017), but falls below some other results. Below ${\rm M}_{\ast}\sim 10^{8}\,{\rm M}_{\sun}$, we predict slightly higher number densities than Song et al. (2016).

In the bottom panel, we show the $z=6$ UV luminosity function (rest-frame 1500 Å) from our predictions (open squares and the thin dashed line, identical to those in Figure 10) and from observations in Atek et al. (2015), Bowler et al. (2015), Bouwens et al. (2015), Finkelstein et al. (2015a), Bouwens et al. (2017b), and Livermore et al. (2017) (symbols with errorbars). First, we only consider the intrinsic stellar luminosities without accounting for dust extinction (the thin dashed line), which results in the fact that our model predicts higher number densities than observed at the bright end.

To quantify the effect of dust attenuation, we use a Monte Carlo method to apply the dust attenuation determined in Section 6.2 to the model UV luminosity function. We adopt the median attenuation from Equation 8, with a magnitude-dependent scatter following a uniform distribution with half-width $\Delta{\rm A_{1500}}=-0.09375\,({\rm M_{1500}}+15)$ at ${\rm M_{1500}}<-15$. The model UV luminosity function after dust extinction is shown by the thick dashed line in Figure 14, which agrees surprisingly well with observations at the bright end. This result suggest that the bright-end of the UV luminosity function is mostly set by dust obscuration, in line with predictions from semi-analytic models (Somerville et al., 2012) and cosmological simulations (Wilkins et al., 2017). We find that dust extinction becomes significant for galaxies with intrinsic UV magnitude brighter than $\rm M_{1500}\sim-20$. The star formation in these galaxies cannot be fully probed in the rest-frame UV. Approximately 37% of the UV light from galaxies brighter than $\rm M_{1500}=-13$ at $z=6$ is obscured by dust according to our model. The obscured fraction is 54% if only galaxies brighter than $\rm M_{1500}=-17$ are considered. These numbers are broadly in line with observational estimates of the dust obscured fraction of star formation at these redshifts (see e.g. Finkelstein et al., 2015a; Bouwens et al., 2015).

Our predicted UV luminosity function (after dust attenuation) is in good agreement with current data in a broad range of magnitudes171717Note that the sample in Livermore et al. (2017) has only one galaxy in the faintest bin, and no galaxy in the next two bins (upper limits)., but the predicted stellar mass function shows considerable discrepancies with observational measurements. We note that the stellar mass functions from different groups also do not in general agree perfectly with each other. We discuss several systematic uncertainties that might be important in these measurements. First, a non-negligible fraction of the light from galaxies will be missed due to the finite surface brightness depth of an observational campaign. Therefore, the stellar mass of a galaxy is possibly underestimated. This effect becomes much stronger at lower masses (e.g. Ma et al., 2017a). Second, the incompleteness correction at the low-mass end for a flux-limited sample is sensitive to the a priori distribution of magnitude at a given stellar mass. We show in Section 3.4 that this distribution is biased toward the faint end (top panels in Figure 8). One could underestimate the incompleteness if this bias is not properly accounted for. Third, measurement uncertainties in stellar mass will introduce contamination from low-mass galaxies in a given mass bin, and thereby lead to an overestimate of their number density, especially at the high-mass end where the stellar mass function is steep (e.g. Davidzon et al., 2017). In addition, cosmic variance may also lead to discrepancies at the massive end.

Nevertheless, our simulations do not include halos more massive than ${\rm M}_{\rm halo}\sim 10^{12}\,{\rm M}_{\sun}$ and only include a small number of independent halos above ${\rm M}_{\rm halo}\sim 10^{11}\,{\rm M}_{\sun}$. We may underestimate the scatter of certain galaxy properties at these masses. Moreover, we do not consider primordial chemistry or the ionizing background fluctuation prior to reionization, which may have important effects on halos below ${\rm M}_{\rm halo}\sim 10^{8}\,{\rm M}_{\sun}$. Our predictions should be tested by future observations to better understand the uncertainties in the current model.

6.4 Differences between stellar population models

In this paper, we use the BPASSv2.0 binary model with a Kroupa (2002) IMF from 0.1–$100\,{\rm M}_{\sun}$ as our default stellar population model for post-processing. To illustrate the difference between stellar population models, we show the synthetic spectrum for simulation z5m12a at $z=5$ (stellar mass ${\rm M}_{\ast}=3\times 10^{9}\,{\rm M}_{\sun}$) in Figure 15 using two models in BPASSv2.0: the default binary model (black) and the single-star model with the same IMF (red). Again, the spectra only include intrinsic star light without accounting for dust attenuation and line emission. The binary models produce slightly higher luminosities in the rest-frame UV at wavelengths bluer than the Balmer break at 3648 Å, but about 0.2–0.8 mag weaker emission in the rest-frame optical and IR than single-star models. This is because of consequences of binary interaction: (a) the ‘effective’ IMF is changed and (b) a large fraction of red supergiants are removed and replaced with hot stripped stars (J. J. Eldridge, private communication). These effects are particularly important in stellar populations younger than 1 Gyr, which are dominant in galaxies at $z\geq 5$ when the age of the Universe is comparable and even younger. If single-star models are used, the predicted B-band and J-band magnitudes will be brighter by 0.2 dex and 0.5 dex, respectively. Future observations of high-redshift galaxies at the rest-frame optical bands will provide more hints of the stellar populations in these galaxies. Another important difference is that binary models tend to produce more ionizing photons (wavelengths shorter than 912 Å). The production of these photons even extends to 30 Myr after the formation of a stellar population (as opposed to 10 Myr in single-star models), so these photons are more likely to escape the galaxy and play an important role in cosmic reionization (Ma et al., 2016b; Götberg et al., 2017). However, ionizing photons only contribute a small fraction (less than 10%) of the bolometric luminosity. Note that the differences between different stellar evolution calculations for non-rotating, single-star models (e.g. starburst99 and BPASS, using the same IMF and stellar atmosphere models) are much smaller than the effects of binaries.

7 Conclusion

In this paper, we present a suite of cosmological zoom-in simulations at $z\geq 5$ covering the $z=5$ halo mass range ${\rm M}_{\rm halo}\sim 10^{8}$–$10^{12}\,{\rm M}_{\sun}$. These are high-resolution simulations (100–$7000\,{\rm M}_{\sun}$ baryonic mass resolution) using physically motivated models of the multi-phase ISM, star formation, and stellar feedback from the FIRE-2 simulation suite of the FIRE project (Hopkins et al., 2017). These simulations provide useful guidance for future observations with JWST and next-generation ground-based telescopes. Our simulations are complementary to simulations of the first stars and low-mass galaxies at $z\gtrsim 15$ with sophisticated primordial chemistry as well as large-volume simulations using empirically calibrated star formation and feedback models at much poorer resolution.

By utilizing all properly resolved halos in each zoom-in region, we obtain a simulated sample containing hundreds of independent halos at a given redshift (Figure 1). We also include all snapshots (separated by about 20 Myr in time) in our analysis to account for time variability in galaxy properties. At low halo masses (e.g. ${\rm M}_{\rm halo}\leq 10^{11}\,{\rm M}_{\sun}$ at $z\sim 6$ or ${\rm M}_{\rm halo}\leq 10^{10.5}\,{\rm M}_{\sun}$ at $z\sim 10$), our sample includes a large number of independent halos to account for halo-to-halo variance. At higher halo masses, our sample is small, so we may underestimate the scatter of galaxy properties due to halo-to-halo variance (cf. Sections 2.3 and 3.2).

We use the BPASSv2.0 binary stellar population models with a Kroupa (2002) IMF to compute the broad-band photometry from starlight in these galaxies. We also use analytic halo mass functions to assign each simulated galaxy a proper number density that reflects its relative abundance in the Universe. In this paper, we study the stellar mass–halo mass relation, star formation histories, the relation between broad-band magnitude and stellar mass, and stellar mass function and luminosity functions. Our main results include the following:

(i) The stellar mass–halo mass relation shows little evolution at redshift $z=5$–12 (Figure 4). The best-fit median relation and $1\sigma$-scatter are $\log{\rm M}_{\ast}=1.58\,(\log{\rm M}_{\rm halo}-10)+7.10$ and $\Delta\log{\rm M}_{\ast}=\exp\,[-0.14\,(\log{\rm M}_{\rm halo}-10)-1.10]$ in the halo mass range ${\rm M}_{\rm halo}=10^{7.5}$–$10^{12}\,{\rm M}_{\sun}$ (Section 3.2). The ${\rm M}_{\ast}$–${\rm M}_{\rm halo}$ relation may bend at ${\rm M}_{\rm halo}>10^{12}\,{\rm M}_{\sun}$ (as is inferred at lower redshifts), but this regime is not probed by our simulations.

(ii) The relation between SFR, halo mass (stellar mass), and redshift can be best described by $\log{\rm SFR}=1.58\,(\log{\rm M}_{\rm halo}-10)+2.20\,\log\left(\frac{1+z}{6}\right)-1.58$ and $\log{\rm SFR}=1.03\,(\log{\rm M}_{\ast}-10)+2.01\,\log\left(\frac{1+z}{6}\right)+1.36$. The slopes of the SFR–${\rm M}_{\rm halo}$ and SFR–${\rm M}_{\ast}$ relations do not depend on redshift, but the average SFR at fixed halo mass (stellar mass) increases with increasing redshift by approximate 0.7 dex from $z=5$ to $z=12$ (Figure 7).

(iii) The mean SFR for galaxies below ${\rm M}_{\rm halo}\sim 10^{8}\,{\rm M}_{\sun}$ or ${\rm M}_{\ast}\sim 10^{4}\,{\rm M}_{\sun}$ below $z\sim 6$ drops significantly (Figure 7), because star formation is suppressed in low-mass galaxies by the ionizing background near the end of reionization (see also the left most panel in Figure 5). About 50% of the halos at ${\rm M}_{\rm halo}\sim 10^{8}\,{\rm M}_{\sun}$ at $z\sim 6$ are dark halos that contain no stars. Halos of similar masses above $z\sim 7$ or halos more massive than ${\rm M}_{\rm halo}\sim 10^{8.5}\,{\rm M}_{\sun}$ at any redshift continue normal star formation.

(iv) We provide the median and $1\sigma$-scatter for the magnitude–stellar mass relation and stellar mass–magnitude relation at rest-frame 1500 Å, B band, and J band (Table 2 and Figure 8). Both relations have large scatter. We emphasize that the two relations are fundamentally different from each other. At fixed stellar mass, the distribution of magnitudes is set by the range of recent star formation histories. At fixed magnitude, the distribution of stellar mass is biased toward the low-mass end, due to the higher abundance of low-mass galaxies in the Universe (Section 3.4).

(v) We predict the stellar mass function and luminosity functions at rest-frame 1500 Å, B band, and J band from $z=5$–12 (Figures 9 and 10). Our results are broadly consistent with current observations (Figures 14) and can be tested by future observations with JWST and next-generation ground-based telescopes. We make our predictions public for future use (see Appendix C for details).

(vi) Both the stellar mass function and luminosity functions show steepening low-mass-end or faint-end slopes with increasing redshift (from $\alpha=-1.85$ at $z\sim 6$ to $\alpha=-2.18$ at $z\sim 12$, Figures 9 and 10), as inherited from the steepening of the low-mass-end slope of the halo mass function.

(vii) The stellar mass function slightly flattens below ${\rm M}_{\ast}\sim 10^{4.5}\,{\rm M}_{\sun}$ at $z\sim 6$. This results from the high dark halo fraction at ${\rm M}_{\rm halo}\sim 10^{8}\,{\rm M}_{\sun}$, due to star formation being suppressed by the ionizing background at these redshifts. Similarly, the $z=6$ luminosity functions also show a flattening at magnitudes fainter than $\rm M_{1500}\sim-12$, $\rm M_{B}\sim-12$ and $\rm M_{J}\sim-12$ (Section 4). There is no such flattening at higher redshifts.

(viii) We derive the star formation rate and stellar mass density at $z=5$–12 (Figures 11 and 12). Our results are in good agreement with current observational constraints at $z\leq 8$. At higher redshifts, both are dominated by low-mass galaxies. Future JWST observations can put more robust constraints on the mass assembly histories at these redshifts by measuring galaxy number densities below $\rm M_{UV}\sim 15$ or ${\rm M}_{\ast}\sim 10^{8}\,{\rm M}_{\sun}$.

(ix) Dust attenuation in the rest-frame UV becomes important for galaxies with intrinsic 1500 Å-magnitude brighter than $\rm M_{1500}\sim-20$ (Figure 13). In our analysis, the bright-end shape of the UV luminosity function is primarily set by dust attenuation (Figure 14). Approximately 37% (54%) of the UV luminosity from galaxies brighter than $\rm M_{1500}=-13$ ($\rm M_{1500}=-17$) is obscured by dust at $z\sim 6$.

We note the caveat that our simulations do not include primordial chemistry and $\rm H_{2}$ formation and dissociation, nor try to model Pop III star formation. These are important in understanding the cooling in primordial gas and metal enrichment at very high redshifts ($z\geq 15$), which may affect the star formation efficiency in halos below ${\rm M}_{\rm halo}\sim 10^{8}\,{\rm M}_{\sun}$ (e.g. Wise et al., 2014; Chen et al., 2014). Furthermore, we do not model cosmic reionization self-consistently; instead, we only apply a spatially uniform, redshift-dependent ionizing background. This ignores the fact that reionization is highly inhomogeneous (e.g. Barkana & Loeb, 2001; Furlanetto et al., 2004; Iliev et al., 2006; McQuinn et al., 2007) and that even after reionization the ionizing background has large spatial fluctuations (e.g. Becker et al., 2015; Davies & Furlanetto, 2016; D’Aloisio et al., 2018). Preliminary results indicate that increasing the ionizing background strength by a factor of 10–100 may lower the stellar mass in halos of $z=5$ mass ${\rm M}_{\rm halo}\sim 10^{9}\,{\rm M}_{\sun}$ by a factor of 2. This may lead to larger scatter in the stellar mass–halo mass relation at the low-mass end that we do not capture in the current study. These questions are worth further exploration.

These simulations have many applications. In the future, we will use them to study the size evolution of high-redshift galaxies, dust attenuation and IR luminosity functions, nebular line emissions, the escape fraction of ionizing photons, [C ii] and CO luminosity functions, metal-enriched absorbers in the circum-galactic medium, Lyman-$\alpha$ radiative transfer, globular cluster formation, and more. We will also expand the simulation suite to lower and higher masses and more extreme environments at these redshifts.

Acknowledgments

We thank J. J. Eldridge, Steven Finkelstein, Renyue Cen, Frank van den Bosch, Priya Natarajan, Avi Loeb, and Rychard Bouwens for helpful discussions. The simulations used in this paper were run on XSEDE computational resources (allocations TG-AST120025, TG-AST130039, TG-AST140023, and TG-AST140064). The analysis was performed on the Caltech compute cluster “Zwicky” (NSF MRI award #PHY-0960291). Support for PFH was provided by an Alfred P. Sloan Research Fellowship, NASA ATP Grant NNX14AH35G, and NSF Collaborative Research Grant #1411920 and CAREER grant #1455342. Support for SGK was provided by NASA through Einstein Postdoctoral Fellowship grant number PF5-160136 awarded by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS8-03060. CAFG was supported by NSF through grants AST-1412836 and AST-1517491, by NASA through grant NNX15AB22G, and by STScI through grant HST-AR-14562.001. EQ was supported by NASA ATP grant 12-APT12-0183, a Simons Investigator award from the Simons Foundation, and the David and Lucile Packard Foundation. MBK was also partially supported by NASA through HST theory grants (programs AR-12836, AR-13888, AR-13896, and AR-14282) awarded by the Space Telescope Science Institute (STScI), which is operated by the Association of Universities for Research in Astronomy (AURA), Inc., under NASA contract NAS5-26555. RF is supported by the Swiss National Science Foundation (grant No. 157591). DK was supported by NSF grant AST-1412153, funds from the University of California, San Diego, and a Cottrell Scholar Award from the Research Corporation for Science Advancement.

Appendix A The weighting method

In Section 2.4, we introduce a weighting method by assigning each halo ‘snapshot’ a weight $w$ according to its halo mass and redshift to reflect its real abundance in the Universe. First, we bin our simulated catalog in $\log{\rm M}_{\rm halo}$–$\log(1+z)$ space with bin sizes $\Delta\log{\rm M}_{\rm halo}=0.4$ and $\Delta\log(1+z)=0.04$ and count the number of halo snapshots in each bin $N_{\rm sim}$ as shown in the left panel of Figure 16. Next, we compute the expected number of halos in the Universe in each bin $N_{\rm expect}=\phi\,\Delta\log{\rm M}_{\rm halo}\Delta V_{\rm com}$, where $\phi$ is the halo mass function obtained from HMFcalc (Murray et al., 2013, number of halos per dex per comoving volume) and $\Delta V_{\rm com}$ is the comoving volume corresponding to the redshift range of each bin and $\rm 1\,arcmin^{2}$ area on the sky (this is to avoid $w$ being too large or too small). Each halo in the same bin will then be given the same weight $w=N_{\rm expect}/N_{\rm sim}$. Therefore, summing over the weights of halos in a given bin leads to the expected number of halos in the Universe (on $\rm 1\,arcmin^{2}$ area of the sky). We show $N_{\rm expect}$ and $w$ in the middle and right panel of Figure 16.

Appendix B Resolution tests

Every zoom-in simulation presented in the paper has been run at several resolution levels. The main text shows results only from the highest-resolution runs available. In Figure 17, we show the stellar mass–halo mass relation in the $z=5$ snapshots for simulations using different mass resolution (shown by different symbols). The large symbols represent the most massive halo in each simulation and smaller symbols show less massive isolated halos in the zoom-in regions with more than $10^{4}$ particles and zero contamination. Note than we show the total stellar mass in the halo instead of the central galaxy stellar mass defined in Section 2.3, to reduce the effects of stochastic fluctuation in galaxy mass induced by mergers. Simulations at mass resolution $m_{b}\sim 5.6\times 10^{4}\,{\rm M}_{\sun}$ tend to systematically over-predict stellar mass by about a factor of two. This is also found in our previous work using ultra-high-resolution dwarf galaxy and Milky Way-mass galaxy simulations run with the same code down to $z=0$ (see Wetzel et al., 2016; Hopkins et al., 2017). At resolution $m_{b}\sim 7\times 10^{3}\,{\rm M}_{\sun}$ and better, we do not find significant systematic differences in the stellar mass–halo mass relation for a fairly large sample of galaxies (as we show for ${\rm M}_{\rm halo}<10^{11}\,{\rm M}_{\sun}$). The difference in the stellar mass of individual galaxy is mainly due to stochastic effects: when and where a star particle forms and a SN occurs are stochastically sampled from the SFR and SNe rates. Any perturbations may affect the final stellar mass of each galaxy, but the statistics in the stellar mass–halo mass relation is unchanged. This is the way we define convergence for our simulations. Therefore, we adopt mass resolution $m_{b}\sim 7\times 10^{3}\,{\rm M}_{\sun}$ for halos above ${\rm M}_{\rm halo}=10^{11}\,{\rm M}_{\sun}$ and even better resolution for our lower mass systems for final production runs, to ensure reasonable convergence and computational costs. For more extensive mass and spatial resolution tests, and other numerical details, see Hopkins et al. (2017).

Appendix C Stellar mass functions and luminosity functions

In Sections 4 and 5, we describe the method used to compute the stellar mass functions and luminosity functions from the simulated sample. Here in Tables 3 and 4, we provide these results at $z=5$–12. The first two columns give the stellar mass functions above ${\rm M}_{\ast}=10^{3.5}\,{\rm M}_{\sun}$, and the remaining columns give the luminosity functions brighter than $\rm M_{AB}=-8$ at rest-frame 1500 Å, B, and J band, respectfully. In addition, we also make our model stellar mass function and luminosity functions public. A machine-readable version of these results is available at http://www.tapir.caltech.edu/~xchma/data/hiz_smf_lf.zip. Those derived from the simulated catalog are tabulated in files SMF_sim_zxx.txt, LF_UV_sim_zxx.txt, LF_B_sim_zxx.txt, and LF_J_sim_zxx.txt (these are identical to the data in Tables 3 and 4). The model stellar mass functions and luminosity functions are tabulated in files SMF_model_zxx.txt, LF_UV_model_zxx.txt, LF_B_model_zxx.txt, and LF_J_model_zxx.txt (these are shown with the dashed lines in Figures 9 and 10). The UV luminosity functions after accounting for dust attenuation are tabulated in LF_UV_red_zxx.txt. The two digits xx in all file names represent the redshift. We encourage readers to use our results and confront them with future observations and other model predictions.