{'fraction_non_alphanumeric': 0.04003497743010422, 'fraction_numerical': 0.029482664902039563, 'mean_word_length': 4.398851674641148, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 14, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Simulating galaxies in the reionization era with FIRE-2: galaxy scaling relations, stellar mass functions, and luminosity functions

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.

impact on the inferred faint-end luminosity functions in the Hubble Frontier Fields. The faint-end slope of the UV luminosity function fainter than MUV = −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 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 restframe 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. Highresolution 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 Mhalo ∼ 10 9 M ) at relatively high redshifts (z 10). These types of simulations have been used to predict the scaling relations of high-redshift, lowmass 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 (MUV > −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 empiricallycalibrated 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 10 5 M . 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 ∼ 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 semianalytic models of galaxy formation (e.g. Lu et al. 2014b).
In this work, we introduce a new suite of cosmological 'zoomin' simulations at z 5 in the z = 5 halo mass range Mhalo ∼ 10 8 -10 12 M . 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 zoomin 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 5: previous works using a similar technique have studied the escape fraction of ionizing photons (e.g. , 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) project 1 . 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.

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) 2 . 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.

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 sidelength 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 ∆vir ≈ 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 Mhalo = 2 × 10 9 -10 12 M , requiring that there is no more massive halo within 5Rvir from the target halo. This selection excludes 1/3 of the halos in the box 3 . . Number of independent halos in the simulated catalog at several redshifts. The simulated sample includes considerably larger numbers of independent halos below M halo = 10 11 M at z ∼ 6 or below M halo = 10 10.5 M at z ∼ 10, where we are able to account for (at least some) halo-to-halo variance. At higher masses, our sample is limited. We can, however, study time-variability and its impact on scatter in observational properties of galaxies.
We identify zoom-in regions based on particles within ∼ 3-5Rvir of the targeted halo, and iterate to ensure zero mass contamination from low-resolution particles within 2Rvir and less than 1% contamination within 3Rvir 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.

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 finestructure and molecular cooling at low temperatures and hightemperature 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 , 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 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 halo (the target halo in the zoom-in region) by design. This will partially compensate the selection bias due to the isolation criteria above. (1) Name: Simulation designation.
(2) M halo : Halo mass of the target halo at z = 5 and its progenitor mass at z = 10.
(3) M * : Stellar mass of the central galaxy in the target halo at z = 5 and its progenitor mass z = 10 (see Section 2.3).
(4) m b : Initial baryonic (gas and star) particle mass in the high-resolution region. A star particle loses about 25% of its initial mass during its entire life due to mass return via supernovae and stellar winds.
(5) b : Minimum Plummer-equivalent force softening for gas particles. Force softening for gas particles is adaptive. The gas inter-particle separation defined in Hopkins et al. (2017) is about 1.4 b .
(6) star: Plummer-equivalent force softening for star particles (fixed in comoving units until z = 9 and in physical units thereafter).
(7) m DM : Dark matter particle mass in the high-resolution region.
(8) DM : Plummer-equivalent force softening for high-resolution dark matter particles (fixed in comoving units until z = 9 and in physical units thereafter).
Faucher- Giguère et al. (2009) 4 , 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 nth = 1000 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 nth 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 M . All feedback quantities are directly calculated from STARBURST99 (Leitherer et al. 1999).
Note that STARBURST99 is a single-star stellar population model 5 , which assumes each star evolves independently, but most 4 The 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 ∼ 7. A tabulated version of the background is available at http://galaxies.northwestern.edu/uvb/. 5 Note that the stellar population models used in the simulations do not including stellar rotation, which is another key ingredient in stellar popula-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 effects 6 . 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) 7 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 ∼ 2-3 galaxies (e.g. Steidel et al. 2016;Strom et al. 2017). In this paper, we mainly tion synthesis (e.g. Leitherer et al. 2014;Choi et al. 2017) and could have important implications for reionization (e.g. Topping & Shull 2015). 6 Binary 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). 7 http://bpass.auckland.ac.nz consider stellar continuum emission, while detailed modeling of dust extinction and nebular line emission will be the subject of future studies.

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 1 3 Rmax from the halo center, where Rmax is the radius at which the maximum circular velocity Vmax is reached (Rmax is already computed by AHF). We thus define a galaxy by including all star particles within 1 3 Rmax after removing contributions from subhalos outside 1 5 Rmax. 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 Rvir and contain more than 10 4 particles in total 8 . 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 Mhalo = 10 11 M at z ∼ 6 and below Mhalo = 10 10.5 M at z ∼ 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. Mhalo 10 11 M at z ∼ 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 (Mhalo 10 11 M ), 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.

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) 9 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 Mhalo-log(1 + z) space with bin widths ∆ log Mhalo = 0.4 from Mhalo = 10 7.5 -10 12 M and ∆ 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 Nsim and calculate the number of halos expected in the Universe Nexpect from the halo mass function. All halos in the same bin are assigned the same weight w = Nexpect/Nsim. In other words, by summing w over all simulated halos in certain halo mass and redshift intervals and dividing i wi 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 M 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.

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 zoomin region. The stellar masses and halo masses are listed in Table  1. We use the BPASSv2.0 binary models to determine the stellar u/g/r-composite images (z = 5) Mock JWST NIRCam F277W images (z = 5) Figure 2. Left: Stellar u/g/r-composite images for the central galaxies of the 'target' halos from Table 1 at z = 5. Right: Noise-free mock JWST NIRCam F277W-band images (rest-frame 4600 Å). The PSF is a Gaussian function with FWHM of 2 pixels. The pixel size is 0.065 arcsec and 0.42 kpc in physical length. The three images in the same row use the same color scale, which spans eight magnitudes in surface brightness, but the depth increases from m AB = 29.5 mag arcsec −2 for the most massive galaxies in the top row to 31.5 mag arcsec −2 for low-mass galaxies in the bottom row. We use the BPASSv2.0 binary models to determine the SED of each star particle from its age and metallicity, and then ray trace along the line-of-sight without dust attenuation. Nebular line emission is also ignored. The scale bar in each panel indicates 1 kpc (physical).
SEDs, assuming a Kroupa (2002) IMF from 0.1-100 M . Note that we only consider intrinsic stellar continuum emission, and ignore dust extinction and nebular line emission at this point 10 . The right panel shows the noise-free mock images as observed by the Near Infrared Camera (NIRCam) on JWST at F277W band (restframe 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) 11 . 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 10 Full 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). 11 The PSF and pixel sizes are adopted from the NIRCam pocket guide from https://jwst.stsci.edu/instrumentation/nircam. mAB = 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 (restframe 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 mAB = 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 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 red- Mock JWST NIRCam F444W images (z = 10) Figure 3. Left: Stellar u/g/r-composite images for all central galaxies at z = 10 (as Figure 2). Right: Noise-free mock JWST NIRCam F444W-band images (rest-frame 4000 Å). The PSF is a Gaussian function with FWHM of 2 pixels. The pixel size is 0.065 arcsec and 0.28 kpc in physical length. The three images in the same row use the same color scale, which spans eight magnitudes in surface brightness, but the depth increases from m AB = 32 mag arcsec −2 for the most massive galaxies in the top row to 36 mag arcsec −2 for low-mass galaxies in the bottom row. The scale bar in each panel indicates 500 pc (physical).
shifts, 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 5 galaxies, as well as determining the completeness of a flux-limited galaxy survey at these redshifts. 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 M * /( f b Mhalo), where f b = Ω b /Ω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 ∆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 lowmass end 12 . 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  The stellar baryon fraction-halo mass relation at the same redshifts. The two-dimensional histograms represent the number of simulated halo snapshots in each pixel in logarithmic scale (as shown by the color scale). All central galaxies that meet the selection criteria described in Section 2.3 are included. The cyan lines show the abundance matching results at z = 5-8 from Behroozi et al. (2013b). The red dashed lines show the best-fit 1σ region of all central galaxies at z = 5-12 (lines are identical in all three panels, see Section 3.2 for details). At each redshift, the stellar mass correlates tightly with halo mass, and there is no significant evolution in the stellar mass-halo mass relation from z = 5-12.

The stellar mass-halo mass relation
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σ dispersion in log M * at every 0.5 dex in log Mhalo from Mhalo = 10 7.5 -10 12 M . We assume a simple power-law relation between M * and Mhalo and a halo mass-dependent scatter and fit the median stellar mass-halo mass relation and 1σ dispersion obtained from the simulated sample as described above. We obtain the best-fit parameters as (α, β, γ, δ) = (1.58, 7.10, −0.14, −1.10).
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σ stellar masshalo 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 M * -Mhalo 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 M * -Mhalo relation is worth future investigation.
How do galaxies evolve on the M * -Mhalo 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 M * -Mhalo plane during the peak of a starburst and reaches some point above the median M * -Mhalo relation, The z = 5 halo masses are labeled in the upper-left corner. Bottom: Stellar mass growth histories for the same galaxies. All high-redshift galaxies show strong 'bursty' star formation histories. The least massive halo (M halo = 10 8 M , left panels) does not form any stars after z ∼ 8: the feedback following a starburst removes the majority of its gas at that time, and the ionizing background prevents fresh gas from accreting and cooling efficiently onto such low-mass halos at late times. More massive halos remain star forming until the end of the simulation at z = 5.
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 M * -Mhalo relation caused by bursty star formation is a physical effect.
There are several caveats in this analysis. First, at the highmass end (Mhalo 10 11 M ), 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-tohalo 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 Mhalo = 10 12 M . 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 M * -Mhalo relation at these redshifts similar to what is seen at lower redshifts (e.g. Behroozi et al. 2013b). We caution that our best-fit M * -Mhalo relation may break down at Mhalo > 10 12 M . Lastly, the best-fit M * -Mhalo relation does not apply to halos with Mhalo 10 8 M below z ∼ 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. 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 Mhalo = 10 8 -4 × 10 11 M (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α-and UVinferred 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. 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 (Mhalo = 10 8 M at z = 5), the feedback from a starburst at z ∼ 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.

Star formation histories
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 M * 10 8 M . 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 ∆ log Mhalo = 0.5 (or ∆ log M * = 0.5) and ∆ 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 Our simulations agree well with observations at the most massive end. At lower masses, the scatter is large as a result of strong bursty star formation histories. Galaxies at lower redshifts have lower SFRs and larger scatter on average than galaxies at higher redshifts at a given stellar mass. . Colors show the weighted average SFR for the simulated sample, relative to the best-fit relation at z = 5 (as labeled beside the colorbars). This emphasizes the dependence on redshift. The white regions represent parameter space where we have no simulations. At fixed halo mass or stellar mass, the mean SFR increases with increasing redshift, by about 0.7 dex from z = 5 to z = 12. Since z ∼ 6, the average SFRs for galaxies below M halo ∼ 10 8 M or M * ∼ 10 4 M decrease by approximately 0.6 dex, because star formation in low-mass galaxies is suppressed by the ionizing background near the end of reionization. and log SFR = α (log M * − 10) + γ log 1 + z 6 + δ , and obtain the best-fit parameters (α, γ, δ) = (1.58, 2.20, −1.58) and (α , γ , δ ) = (1.03, 2.01, 1.36). In Figure 7, we show the SFR(Mhalo, z) (left panel) and SFR(M * , 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 ∼ 0.7 dex from z = 5-12, in broad agreement with the qualitative trend shown in Figure 6.
Furthermore, after z ∼ 6, the average SFRs in low-mass galaxies (halo mass below Mhalo ∼ 10 8 M or stellar mass below M * ∼ 10 4 M , 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 7) or more massive halos (Mhalo 10 8.5 M ) 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 Mhalo = 10 8.5 M at any redshift, whereas at Mhalo ∼ The magnitude-stellar mass relation for rest-frame 1500 Å at z = 6 (left) and z = 10 (middle), and for rest-frame B band at z = 6 (right). The two-dimensional histogram represents the total weight of simulated galaxy snapshots in each pixel (in logarithmic scale), taking into account our simulation results and the halo mass function (see Section 2.4 for details). This reflects the true relative number of galaxies between pixels in the Universe. The red solid and dashed lines show the best-fit weighted median relation and 1σ scatter (see text for details). The cyan dotted line in the middle panel shows the median relation at z = 6 for reference. At fixed stellar mass, the distribution of magnitudes is asymmetric, with a broader spread at the bright end, and the median magnitude becomes more negative (galaxies being brighter) at higher redshifts. The scatter gets smaller from rest-frame UV to longer wavelength. Bottom: The stellar mass-magnitude relation for the same bands and redshifts. At fixed magnitude, the distribution of stellar mass is skewed toward low-mass galaxies, simply due to the fact that low-mass galaxies are more abundant in the Universe. The median stellar mass decreases by about 1 dex from z = 12 to 5. 10 8 M , the dark fraction increases from less than 1% at z = 12 to approximate 50% at z = 5. 13 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 ∼ 6. Our findings are broadly in line with other simulations in the literature. For example, Wise et al. (2014) found no dark halos above Mhalo ∼ 10 8 M at z > 8. Sawala et al. (2016) found that the dark fraction decreases sharply from nearly 90% at Mhalo ∼ 10 8 M to 0% at Mhalo ∼ 10 8.5 M at z ∼ 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 radiativehydrodynamic 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 13 Note that the increasing dark fraction at M halo ∼ 10 8 M and below z ∼ 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.
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.

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 M 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 magnitudestellar mass relation (top panels) for our simulated sample with M * > 10 3.5 M 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 MAB = −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 Note: All magnitudes are derived from intrinsic stellar luminosities without accounting for dust attenuation and nebular line emission. galaxies in the Universe between pixels 14 . The red solid and dashed lines illustrate the best-fit median relation and 1σ-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 MAB,band = a (log M * − 10) + c log 1 + z 6 + d, where we assume a fixed slope a at any redshift but a redshiftdependent normalization to capture this feature. We fit the weighed median, 1σ lower-and upper-bound relations (above M * = 14 We 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 M halo > 10 11 M where our sample does not contain large number of independent halos (see Section 2.3 for details).  convolving between the stellar masshalo mass relation from Section 3.2 and the halo mass function, assuming each halo contains one central galaxy. Both methods only account for halos more massive than M halo = 10 7.5 M . The two stellar mass functions agree well with each other for a broad range of mass and redshift. The low-massend slope steepens with increasing redshift (from α = −1.83 at z ∼ 6 to α = −2.18 at z ∼ 12). At z ∼ 6, the stellar mass function derived from simulations flattens and falls below the model stellar mass function by a factor of 2 below M * ∼ 10 4.5 M , owing to the 50% fraction of dark halos around M halo ∼ 10 8 M . A comparison with observations is shown later in Section 6.3 ( Figure 14). We make our predictions publicly available (see Appendix C for details). 10 3.5 M ) 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 massmagnitude relation and fit the weighted median and 1σ relations for galaxies brighter than MAB,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.  Figure 9, the faint-end slope steepens with increasing redshift. The stronger flattening of the z = 6 luminosity functions fainter than M AB ∼ −12 is due to the suppression of star formation in halos around M halo ∼ 10 8 M by the strong ionizing background. We will compare our predicted UV luminosity function with observations later in Section 6.3 and Figure 14. We make our predictions publicly available (see Appendix C for details).

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 ∆z = ±0.5 from our simulated catalog. We then add the weights of galaxies in stellar mass bins and divide i wi by the comoving volume corresponding to the ∆z = ±0.5 redshift interval. Only halos above Mhalo = 10 7.5 M are taken into account. In Figure 9, we show the results in number dex −1 Mpc −3 above M * = 10 3.5 M at z = 6, 8, 10, and 12 with the open symbols. The data are tabulated in Tables  C1 and C2 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 Mhalo = 10 7.5 M 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σ 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 φ d log M * ∼ M α+1 * d log M * ) steepens with increasing redshift, with the slope decreasing from α = −1.80 ± 0.02 at z = 6 to α = −2.13 ± 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 computed 15 . 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 M * ∼ 10 3.5 M is an expected artifact because we exclude all halos below Mhalo = 10 7.5 M .
More importantly, the z = 6 stellar mass function derived from the simulated sample shows a flattening below M * ∼ 10 4.5 M and falls below the model mass function by a factor of 2. This is caused by the 50% fraction of dark halos at Mhalo ∼ 10 8 M at z ∼ 6. In other words, the assumption we adopted in the model that every halo hosts one galaxy breaks down at Mhalo ∼ 10 8 M . 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 15 The 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. The z 9 SFRD is poorly constrained due to the small size of each observational sample. Our simulations suggest that low-mass halos dominate the SFRD, due to their rapidly increasing number densities at these redshifts. This is beyond the detection limits of current observational facilities (M UV ∼ −17, roughly corresponding to halo mass M halo ∼ 10 9.5 M in our simulations). Future deep surveys by JWST will be able to put stronger constrains on the z 9 SFRD. reionization (e.g. Bullock et al. 2000;Somerville 2002;Boylan-Kolchin et al. 2014).

MULTI-BAND LUMINOSITY FUNCTIONS
We calculate the luminosity functions at several bands following the same method described in Section 4. Again, only halos above Mhalo = 10 7.5 M are taken into account. In Figure 10, we show the results in number 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 φ dMAB,band ∼ 10 −0.4(α+1)M AB, band dMAB,band) steepens with increasing redshift (e.g. from α = −1.85 ± 0.06 at z = 6 to α = −2.17 ± 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 We now derive the cosmic SFR density (SFRD) from z = 5-12 by convolving the SFR-Mhalo 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 Mhalo = 10 7.5 M (dashed), Mhalo = 10 8.5 M (solid), and Mhalo = 10 9.5 M (dotted) to Mhalo = 10 12 M . 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 agree with data within observational uncertainties. The SFRD at z 9 are still poorly constrained observationally. These results are derived by integrating the best-fit UV luminosity functions brighter than MUV ∼ −17. This limit does not correspond to a unique halo mass, but is roughly consistent with what we obtain by integrating down to Mhalo = 10 9.5 M at these redshifts (cf. Figures 4 and  8). Note that the number of galaxies in the observed z 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 Mhalo = 10 9.5 M at z 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 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 M * > 10 6 M (solid), M * > 10 4 M (dashed), and M * > 10 8 M (dotted). At these redshifts, high-mass galaxies (M * > 10 10 M ) 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 M * > 10 8 M . 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).

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. 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 raytracing 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. 2016Feldmann et al. , 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 × 10 3 cm 2 g −1 at 1500 Å at solar metallicity. For each simulated galaxy, we include all gas particles within 1 2 Rmax (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 A1500 = −2.5 log(F1500/F1500,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  Figure 13. The relation between UV extinction and intrinsic UV magnitude in the simulations. The points show all simulated galaxies from z = 5-8, with ten different sightlines for each galaxy. We determine the dust extinction by ray-tracing and assuming a constant dust-to-metal ratio, using the distribution of gas and metals in the simulations. The blue dashed line shows the best-fit relation in Equation 8. end. We fit the results with a parabolic function A1500 = (0.0306±0.0002) (M1500 + 15) 2 , and quote a uniformly distributed scatter with half-width ∆A1500 = −0.125 (M1500 + 15) at M1500 < −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.

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 ). 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  Stefanon et al. (2017). Our results are broadly consistent with observations. The discrepancies might be due to systematic uncertainties in deriving stellar mass from single-band magnitude, incompleteness corrections at the low-mass end, and cosmic variance at the massive end. Bottom: The UV luminosity at z = 6. A compilation of observations from 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) are shown by symbols with errorbars. Using intrinsic luminosities, the model tends to predict higher number densities than observed at the bright end. The thick dashed line shows the luminosity function after accounting for the dust extinction in the simulations (see text for details). The good agreement with data suggests that the turnover at the bright-end of the UV luminosity function is largely due to dust extinction. Approximately 37% (54%) of the UV luminosity from galaxies brighter than M 1500 = −13 (−17) is obscured by dust at z ∼ 6. scatter following a uniform distribution with half-width ∆A1500 = −0.09375 (M1500 + 15) at M1500 < −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 M1500 ∼ −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 M1500 = −13 at z = 6 is obscured by dust according to our model. The obscured fraction is 54% if only galaxies brighter than M1500 = −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 magnitudes 17 , 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 Mhalo ∼ 10 12 M and only include a small number of independent halos above Mhalo ∼ 10 11 M . 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 Mhalo ∼ 10 8 M . Our predictions should be tested by future observations to better understand the uncertainties in the current model.

Differences between stellar population models
In this paper, we use the BPASSv2.0 binary model with a Kroupa (2002) IMF from 0.1-100 M 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 M * = 3 × 10 9 M ) in Figure Figure 15. Synthetic spectrum for a M * = 3 × 10 9 M galaxy (simulation z5m12a at z = 5) using different stellar population models. Binary models produce slightly higher rest-frame UV luminosities at 1500 Å, but approximately 0.2-0.8 mag lower rest-frame optical-to-IR luminosities than single-star models. Moreover, binary models produces more ionizing photons (wavelengths shorter than 912 Å), but these photons only contribute a small fraction of the bolometric luminosity. 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 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 restframe 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. STAR-BURST99 and BPASS, using the same IMF and stellar atmosphere models) are much smaller than the effects of binaries.

CONCLUSION
In this paper, we present a suite of cosmological zoom-in simulations at z 5 covering the z = 5 halo mass range Mhalo ∼ 10 8 -10 12 M . These are high-resolution simulations (100-7000 M 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 . 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 lowmass galaxies at z 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. Mhalo 10 11 M at z ∼ 6 or Mhalo 10 10.5 M at z ∼ 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σscatter are log M * = 1.58 (log Mhalo − 10) + 7.10 and ∆ log M * = exp [−0.14 (log Mhalo − 10) − 1.10] in the halo mass range Mhalo = 10 7.5 -10 12 M (Section 3.2). The M * -Mhalo relation may bend at Mhalo > 10 12 M (as is inferred at lower redshifts), but this regime is not probed by our simulations.
(iii) The mean SFR for galaxies below Mhalo ∼ 10 8 M or M * ∼ 10 4 M below z ∼ 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 Mhalo ∼ 10 8 M at z ∼ 6 are dark halos that contain no stars. Halos of similar masses above z ∼ 7 or halos more massive than Mhalo ∼ 10 8.5 M at any redshift continue normal star formation.
(iv) We provide the median and 1σ-scatter for the magnitudestellar mass relation and stellar mass-magnitude relation at restframe 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 α = −1.85 at z ∼ 6 to α = −2.18 at z ∼ 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 M * ∼ 10 4.5 M at z ∼ 6. This results from the high dark halo fraction at Mhalo ∼ 10 8 M , 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 M1500 ∼ −12, MB ∼ −12 and MJ ∼ −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 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 MUV ∼ 15 or M * ∼ 10 8 M .
(ix) Dust attenuation in the rest-frame UV becomes important for galaxies with intrinsic 1500 Å-magnitude brighter than M1500 ∼ −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 M1500 = −13 (M1500 = −17) is obscured by dust at z ∼ 6.
We note the caveat that our simulations do not include primordial chemistry and H2 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 15), which may affect the star formation efficiency in halos below Mhalo ∼ 10 8 M (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 Mhalo ∼ 10 9 M 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-α 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 Mhalo-log(1 + z) space with bin sizes ∆ log Mhalo = 0.4 and ∆ log(1 + z) = 0.04 and count the number of halo snapshots in each bin Nsim as shown in the left panel of Figure A1. Next, we compute the expected number of halos in the Universe in each bin Nexpect = φ ∆ log Mhalo∆Vcom, where φ is the halo mass function obtained from HMFcalc , number of halos per dex per comoving volume) and ∆Vcom is the comoving volume corresponding to the redshift range of each bin and 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 = Nexpect/Nsim. Therefore, summing over the weights of halos in a given bin leads to the expected number of halos in the Universe (on 1 arcmin 2 area of the sky). We show Nexpect and w in the middle and right panel of Figure A1.

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 B1, 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 zoomin 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 ∼ 5.6 × 10 4 M 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 ∼ 7 × 10 3 M 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 Mhalo < 10 11 M ). 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 ∼ 7 × 10 3 M for halos above Mhalo = 10 11 M 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 C1 and C2, we provide these results at z = 5-12. The first two columns give the stellar mass functions above M * = 10 3.5 M , and the remaining columns give the luminosity functions brighter than MAB = −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 C1 and  C2). 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 Figure A1. The number of halo snapshots in our simulated catalog (left), the expected number of galaxies in the Universe (on 1 arcmin 2 sky, center), and the weight assigned to each halo snapshot (right) as a function of halo mass and redshift.  Figure B1. The stellar mass-halo mass relation produced by simulations at different mass resolution. Simulations at mass resolution m b ∼ 5.6 × 10 4 M systematically produce two times more stars. Simulations at resolution m b ∼ 7 × 10 3 M and better do not show statistically-significant systematic differences over a fairly large number of galaxies (as we are able to show below M halo ∼ 10 11 M ). The difference in the stellar mass of individual galaxy is usually due to stochastic star formation and feedback. the redshift. We encourage readers to use our results and confront them with future observations and other model predictions. Notes: the magnitudes are intrinsic magnitude without dust attenuation. Notes: the magnitudes are intrinsic magnitude without dust attenuation.