DRAGONS – Galaxy formation and the EoR

Dark-ages reionization and galaxy formation simulation - III. Modelling galaxy formation and the epoch of reionization


We introduce Meraxes, a new, purpose-built semi-analytic galaxy formation model designed for studying galaxy growth during reionization. Meraxes is the first model of its type to include a temporally and spatially coupled treatment of reionization and is built upon a custom (100 Mpc) N-body simulation with high temporal and mass resolution, allowing us to resolve the galaxy and star formation physics relevant to early galaxy formation. Our fiducial model with supernova feedback reproduces the observed optical depth to electron scattering and evolution of the galaxy stellar mass function between =5 and 7, predicting that a broad range of halo masses contribute to reionization. Using a constant escape fraction and global recombination rate, our model is unable to simultaneously match the observed ionizing emissivity at . However, the use of an evolving escape fraction of 0.05–0.1 at , increasing towards higher redshift, is able to satisfy these three constraints. We also demonstrate that photoionization suppression of low mass galaxy formation during reionization has only a small effect on the ionization history of the inter-galactic medium. This lack of ‘self-regulation’ arises due to the already efficient quenching of star formation by supernova feedback. It is only in models with gas supply-limited star formation that reionization feedback is effective at regulating galaxy growth. We similarly find that reionization has only a small effect on the stellar mass function, with no observationally detectable imprint at . However, patchy reionization has significant effects on individual galaxy masses, with variations of factors of 2–3 at =5 that correlate with environment.

galaxies: formation – galaxies: high redshift – dark ages, reionization, first stars


1 Introduction

There are several key observational areas in which substantial progress will be made in the study of the first galaxies during the coming decade. Of particular importance will be forthcoming programmes searching for galaxies beyond the current redshift frontier using the Hubble Space Telescope and, in the future, the James Webb Space Telescope (e.g. Bouwens et al., 2011; McLure et al., 2013; Schenker et al., 2013). However, even next generation surveys will not extend to the faint luminosities of the faintest galaxies thought to drive the reionization of inter-galactic neutral hydrogen in the early Universe (Robertson et al., 2013; Duffy et al., 2014). Thus, alongside new probes provided by high redshift gamma-ray bursts (e.g. Trenti et al., 2015) and metal pollution of the inter-galactic medium (IGM; e.g. Díaz et al., 2014), an important new observational window for study of the first galaxies will be provided by experiments to measure the redshifted radio signal (Furlanetto, 2006; Morales & Wyithe, 2010). These observations will both provide the first direct probe of the neutral hydrogen content in the high redshift Universe and, through modelling, provide a route to study the early dwarf galaxies thought to exist during reionization alongside their more massive counterparts whose star formation can be directly detected.

Within this context, the development of theoretical models that include a self-consistent treatment of the physics of galaxy formation and intergalactic hydrogen will play a key role. Traditional approaches to the study of galaxies and their effects on the IGM utilize either numerical simulation or analytic modelling. The latter allows investigation of average behaviours on large scales but the calculations are inherently linear, meaning that complex feedback processes cannot be addressed (e.g. Furlanetto et al., 2004; Wyithe & Loeb, 2004, 2013). Numerical simulations, on the other hand, include non-linear effects but at the expense of computational cost. To achieve a volume sufficiently large to study ionized structure, a popular and effective approach to simulating reionization is to begin with a collisionless N-body simulation (e.g. Ciardi et al., 2003; Sokasian et al., 2003; Iliev et al., 2007, 2008; Trac & Cen, 2007; Zahn et al., 2007; Shin et al., 2008; Trac et al., 2008) and use a simple prescription to relate halo mass to ionizing luminosity. A radiative transfer method (for example ray-tracing algorithms) can then be used to calculate the ionization structure on large scales. In recent years, new hybrid, or semi-numerical models (Mesinger & Furlanetto, 2007; Geil & Wyithe, 2008; Kim et al., 2013a) have been developed that combine N-body simulations with analytical methods to enable the calculation of reionization structure in very large volumes with high efficiency. These methods have elucidated the primary features of the ionization structure during reionization, but do not capture the physics of galaxy formation.

Therefore, to better understand the physics of galaxy formation, many authors have performed hydrodynamic simulations of galaxy formation (Finlator et al., 2011; Salvaterra et al., 2011; Jaacks et al., 2012, e.g.) which are able to directly model the growth of stellar mass in high-redshift galaxies when coupled with sub-grid models for processes including metal enrichment and feedback. These simulations are able to broadly reproduce the luminosity function of galaxies at high redshift, however, computational expense limits their ability to self-consistently model reionization in volumes large enough to statistically describe the spatial evolution of this process. Instead, a common approach is to impose a simple parametrized model to approximate the average ionizing background as a function of redshift, independent of the properties of the ionizing source population (e.g. Feng et al., 2016) or their spatial distribution (e.g. Genel et al., 2014; Schaye et al., 2015). Recently, hydrodynamical simulations of galaxy formation with coupled radiative transfer have been used to compute the effects of reionization on galaxy formation self-consistently for the first time (e.g. So et al., 2014; Norman et al., 2015; Ocvirk et al., 2015; Pawlik et al., 2015). However, the extreme computational expense of these simulations limit their size to relatively small volumes and/or few variations on galaxy formation physics or reionization scenarios that can be explored. In addition, the modelling of sub-grid physical processes remains uncertain, requiring systematic studies of the available parameter space in order to draw robust conclusions. Such studies represent an extreme computational challenge which has yet to be overcome.

Another approach to the realistic modelling of high redshift galaxies has been through the use of semi-analytic galaxy formation models (Benson et al., 2006; Lacey et al., 2011; Raičević et al., 2011; Zhou et al., 2013). While large volumes are available to such models, until now they have not been fully coupled to an accurate description of reionization. This is in part due to the structure of most existing semi-analytic models, which utilize so-called ‘vertical’ halo merger trees (e.g. Springel et al., 2005; Bower et al., 2006; Harker et al., 2006; De Lucia & Blaizot, 2007) in which galaxies belonging to each tree branch are evolved independently from the rest of the simulation volume. Since galaxies drive the process of reionization, which in turn affects their subsequent evolution, galaxies spatially separated by tens of Mpc cannot be considered and evolved independently as has traditionally been the case (Wyithe & Loeb, 2004). Self-consistently studying reionization instead requires a semi-analytic model designed to run on ‘horizontal’ merger trees where all haloes at each snapshot of the parent N-body simulation are processed simultaneously. Additionally, the reduced dynamical time of dark matter haloes at high redshift requires snapshots with a much higher cadence than is needed to model galaxy formation at lower redshifts.

This is the third paper in a series describing the Dark-ages Reionization And Galaxy Observables from Numerical Simulations (DRAGONS) project3, which integrates detailed semi-analytic models constructed specifically to study galaxy formation at high redshift, with semi-numerical models of the galaxy–reionization process interaction. In this work, we introduce Meraxes, the new semi-analytic model of galaxy formation developed for DRAGONS, integrating the 21cmFAST semi-numerical model for ionization structure described in Mesinger & Furlanetto (2007). Meraxes is implemented within the large-volume, high-resolution, and high-cadence Tiamat N-body simulation described in Poole et al. (2016, hereafter Paper I) and Angel et al. (2016, Paper II). In subsequent papers we will use Meraxes to carry out a range of studies including the investigation of the high redshift galaxy luminosity function (Liu et al., 2015, Paper IV), and the ionization structure of the IGM (Geil et al., 2015, Paper V). Complimentary, high resolution hydrodynamic simulations are described in Duffy et al. (2014).

The outline of the paper is as follows. In Section 2 we provide a full description of Meraxes including an overview of the Tiamat simulation and associated merger trees which act as input for our semi-analytic model, the physical prescriptions employed, and the methodology of our reionization coupling and the integration of 21cmFAST. In Section 3 we then go on to describe the calibration of the model’s free parameters against the evolution of the high-redshift galaxy stellar mass function. In Section 4 we investigate a range of different extreme reionization and galaxy physics modifications in order to elucidate the roles of reionization suppression and galactic feedback processes in the build up of stellar mass and the evolution of the global neutral hydrogen fraction. We also highlight the important consequences of utilizing a patchy, self-consistent reionization model compared to more commonly employed, parametrized descriptions of reionization. Finally, in Section 5 we summarize our study and conclusions. Throughout this work, we employ a standard, spatially flat cold dark matter cosmology with the most up-to-date cosmological parameters as determined by Planck Collaboration (2015): (, , , , , )=(0.678, 0.308, 0.0484, 0.692, 0.815, 0.968).

2 Meraxes

Modern semi-analytic galaxy formation models are capable of providing statistically accurate representations of the global properties of galaxies across a broad range of redshifts (e.g. Baugh, 2006; Mutch et al., 2013; Henriques et al., 2015), and are therefore able to describe the distribution and evolution of the ionizing photons which drive the process of cosmic reionization. These photons generate regions of ionized hydrogen (H ii) with characteristic sizes of tens of during reionization (Wyithe & Loeb, 2004). Thus, in order to take advantage of this information and to self-consistently model the effect of these photons on the growth of galaxies, one must consider the contributions of galaxies separated by similar scales.

Traditionally, semi-analytic models have therefore used parametrized descriptions to include the average effect of reionization and the associated photo-suppression of baryonic infall on the growth of galaxies (Benson et al., 2006; Croton et al., 2006; Somerville et al., 2008). These parametrizations are typically calibrated using radiative transfer simulations and are provided as a function of redshift and halo mass alone (e.g. Gnedin, 2000). Whilst it is computationally efficient to include reionization in this manner, there are a number of important drawbacks. First, the progression of reionization is not self-consistently modified by the growth of the galaxies which are driving it. Therefore it is impossible to investigate how different galaxy physics affect the ionization state of the IGM or to quantify the potential back-reaction on galaxy evolution. Secondly, these simple reionization prescriptions miss the potentially important effects of spatially dependent self-regulation (Iliev et al., 2007; Sobacchi & Mesinger, 2013), whereby massive galaxies located at the peaks in the density distribution can reionize their surroundings, delaying or preventing the onset of star formation in nearby lower mass haloes.

Our new semi-analytic galaxy formation model, Meraxes, has been written from the ground up to facilitate these modelling requirements. Its key features include the ‘horizontal’ processing of merger trees constructed from a purpose run N-body simulation (Tiamat; see Section 2.1 below) and the incorporation of the semi-numerical reionization algorithm, 21cmFAST, as a core component. When combined, these features allow Meraxes to efficiently couple the growth of galaxies to the process of reionization, both temporally and spatially. It can therefore be used to investigate the potentially complex effects of various reionization models on the properties of high- galaxies, as well as to test for observational discriminants of different galaxy physics in the distribution and evolution of inter-galactic neutral hydrogen.

In order to develop confidence in our newly developed framework, as well as provide a solid foundation for future additions and improvements, our initial implementation of the baryonic physics processes in Meraxes is heavily based on the well-studied L-Galaxies semi-analytic model (Kauffmann, 1996; De Lucia & Blaizot, 2007; Guo et al., 2013; Henriques et al., 2015), in particular the version described in Croton et al. (2006) and extended in Guo et al. (2011). However, as well as our improved treatment of reionization, the excellent temporal resolution provided to us by the Tiamat merger trees has also necessitated the development of a number of important updates to the treatment of supernova feedback and stellar mass recycling.

In the following sub-sections we describe Meraxes in full, including its input data set in the form of halo merger trees extracted from the Tiamat N-body simulation, the details of the implemented galaxy physics prescriptions, and our methodology for integrating 21cmFAST to self-consistently model reionization.

2.1 Input – the Tiamat N-body simulations

The Tiamat collisionless N-body simulation has been designed for the DRAGONS study of high-redshift galaxy formation and the epoch of reionization (EoR). It contains dark matter particles within a (comoving) periodic box and was run using a modified version of the GADGET-2 N-body code and the latest Planck 2015 (Planck Collaboration, 2015) cosmology. The volume of Tiamat allows for the investigation of the statistical signatures of reionization and its observational signal, whilst the resulting particle mass of provides the necessary resolution to identify the low-mass sources thought to be driving this process. Furthermore, Tiamat provides high temporal resolution in the form of 100 output snapshots evenly spaced in cosmic time between and 5, resulting in a cadence of per snapshot. This level of temporal resolution is a unique feature of Tiamat which allows our semi-analytic model to accurately simulate the stochastic nature of star formation in a regime where the dynamical time of a typical galactic disc is shorter than the lifetime of the least massive Type II supernova progenitor ().

In addition to the main Tiamat volume, a suite of smaller, higher mass resolution N-body simulations have been run as part of the DRAGONS programme (Paper-I). For this work we make particular use of the Tiny Tiamat and Medi Tiamat volumes in order to quantify the effect of resolution on our results (see Section 4.1 and Appendix A). Tiny Tiamat is the highest resolution simulation of the DRAGONS suite, with a particle mass of in a small box of side length , whilst Medi Tiamat bridges the resolution gap with the main simulation by providing a particle mass of in a box. Both simulations maintain the same snapshot cadence as the main Tiamat volume and are described in detail in Paper-I.

Halo identification in all simulations used in this work was carried out using the Subfind (Springel et al., 2001) real-space halo finder down to a minimum mass of 32 particles (corresponding to , , for Tiamat, Medi Tiamat, and Tiny Tiamat respectively). The resulting halo catalogues comprise of friends-of-friends (FoF) groups of gravitationally bound particles which themselves are made up of a single mass dominant ‘central’ subhalo along with zero or more sub-dominant ‘satellite’ subhaloes. For further details, interested readers are referred to Papers I and II.

Merger trees

The formation history of subhaloes, in the form of hierarchical merger trees, acts as the raw input to Meraxes and is used to define the positions and growth of galaxies. Many traditional semi-analytic models, such as the L-Galaxies (e.g. De Lucia & Blaizot, 2007; Guo et al., 2011; Henriques et al., 2015) and GALFORM (e.g. Bower et al., 2006; Lagos et al., 2012; Kim et al., 2013b) variants, process such trees in a depth-first (or ‘vertical’) order, whereby small collections of directly interacting dark matter haloes are processed one after the other from high to low redshift and independently of each other. Whilst computationally efficient in terms of minimizing the memory overhead required to process the simulation, the inherent assumption is that haloes (and by extension galaxies) which do not directly interact do not affect each other’s evolution. This assumption breaks down when considering the process of reionization during which ionizing photons from galaxies tens of Mpc away can heat the IGM, raising the local Jeans mass and altering the accretion rate of baryons (Dijkstra et al., 2004). Meraxes instead processes trees breadth-first (or ‘horizontally’). In this method all of the haloes in the entire volume are loaded into memory and the associated galaxies evolved for each snapshot sequentially. This allows Meraxes to more efficiently model reionization than previous comparable works (e.g. Kim et al., 2013a). More detailed information, including the precise details of our merger tree construction technique, can be found in Poole et al. (in preparation).

2.2 Baryonic infall

We begin by making the standard assumption that as FoF groups grow, any freshly accreted mass, always carries with it the universal baryon fraction, , in the form of pristine primordial gas. However, the fraction of these infalling baryons which will remain bound to the FoF group and participate in galaxy formation may be reduced by a number of factors. In particular, ionizing ultraviolet background (UVB) radiation from both local and external sources can heat the IGM, increasing the local Jeans mass and leading to a non-negligible reduction in the amount of baryons successfully captured by low mass systems (Dijkstra et al., 2004). We parametrize this reduction in terms of a baryon fraction modifier, , which represents the attenuation of the total baryon mass that could have ever been successfully captured by an FoF group in its lifetime:


where is the infalling baryonic mass, , and is the number of galaxies in the FoF group. The baryonic reservoirs , , , and are described in the following sections along with the physical prescriptions which govern their evolution. If the mass of the FoF group or the value of decreases then it is possible for to become negative. In this case baryons are stripped from the system as described in Section 2.8. An accurate spatially and temporally dependent calculation of the value of is a key feature of Meraxes and a subject which we return to in detail in Section 2.11.

Any baryons which are successfully captured are assumed to be shocked to the virial temperature of the host FoF group and added to a diffuse hydrostatic hot reservoir where they mix with any already present hot gas.

2.3 Cooling

At each time step in the simulation some fraction of the hydrostatic hot reservoir may cool and condense down into the central regions of the group where it can then participate in galaxy formation. In order to calculate the rate at which this occurs we follow the commonly employed methodology outlined in White & Frenk (1991). In this model, the cooling time of a quasi-static isothermal hot halo is given by the ratio of the specific thermal energy to cooling rate per unit volume:


where is the mean particle mass ( for a fully ionized gas), is the Boltzmann constant, is the cooling function (Sutherland & Dopita, 1993), is the temperature of the gas and is its density profile. As mentioned above, we assume that the hot gas is shocked to the virial temperature of the FoF group, therefore we set . For simplicity, we also assume that the hot gas follows a singular isothermal sphere density profile:


With knowledge of the cooling time, we can define an appropriate cooling radius, , within which there is enough time for the material to lose pressure support and condense to the system centre. Following Croton et al. (2006), we take this to be the radius at which is equal to the dynamical time of the host FoF group, . As discussed by White & Frenk (1991), this model for cooling naturally leads to three distinct regimes.

  1. When , any infalling gas will cool so rapidly that there will be no time for a stable shock to form and thus for the gas to reach hydrostatic equilibrium. In this case we assume that the infalling material flows directly into the central regions of the halo over a dynamical (free-fall) time, .

  2. When the cooling time will be sufficiently long that a quasi-static hot atmosphere will form. The cooling rate from this atmosphere can then be calculated from a simple continuity equation for the mass flux across the evolving cooling radius:

  3. When , we set and no cooling occurs. In the standard model of galaxy formation, haloes with this temperature represent the lowest mass scale for galaxy formation. Below this, the primary mechanism for gas cooling is via molecular hydrogen which is easily photo-dissociated by trace amounts of star formation, making it an inefficient pathway for Pop II star formation. Above this temperature, atomic line cooling provides an efficient mechanism to dissipate energy and remove pressure support (Barkana & Loeb, 2001). The mass resolution of our input N-body simulation, Tiamat, was chosen such that the minimum halo mass at is close to the atomic cooling mass threshold of . Although earlier Pop III and Pop II star formation is possible in smaller haloes, the level of contribution of these objects to reionization remains unclear, being heavily dependent on the masses of the first supernovæ which could potentially delay future star formation by tens to hundreds of Myr (e.g. Chen et al., 2014; Jeon et al., 2014). We therefore do not include these objects in our current model.

All material which successfully cools into the central regions of the FoF group is assumed to be deposited directly into the cold gas reservoir of the galaxy hosted by the central halo. This assumption is commonly employed by a number of other semi-analytic models (e.g. Bower et al., 2006; De Lucia & Blaizot, 2007; Somerville et al., 2008; Lu et al., 2011a) which utilize halo catalogues created by the Subfind halo finder. It is also warranted in the vast majority of systems where the central halo dominates the mass of the FoF group and ensures a physically meaningful match between the galaxy formation physics of Meraxes and the substructure hierarchy produced by the halo finder employed for this programme (see Paper I where this issue is raised).

2.4 Star formation

As discussed in the previous section, gas which cools from the FoF group hot reservoir is assumed to be deposited into the galaxy hosted by the central halo of the group. Here we assume that it settles into a rotationally supported cold gas disc with an exponential surface density profile. Under the simplifying assumption of full conservation of specific angular momentum, the scale radius of the disc can be approximated from the spin of the host dark matter halo, , to be , where we use the definition of provided by Bullock et al. (2001).

Based on the well-established observational work of Kennicutt (1998), the star formation rate of local spiral galaxies can be related to the surface density of cold gas above a given threshold. The value of this threshold can be understood in terms of the gravitational instability required to form massive star-forming clouds (Kennicutt, 1989). Assuming a constant gas velocity dispersion and a flat rotation curve with a circular velocity equal to that of the host dark matter halo, , Kauffmann (1996) demonstrated that this stability criterion can be expressed as


Kauffmann (1996) originally assumed a thin isothermal disc with a gas velocity dispersion appropriate for low-redshift spiral galaxies of , resulting in . However, both observations and simulations of high-redshift star forming galaxies indicate that they typically possess highly turbulent, clumpy discs (e.g. Wisnioski et al., 2011; Glazebrook, 2013; Bournaud et al., 2014). This suggests the need for a modified value. Given the uncertainty in what value this should take, we choose to leave it as a free parameter in our model. We note that Henriques et al. (2015) also advocate allowing freedom in the choice of ; however, they instead motivate this by the observation that star formation is more closely linked to molecular, rather than total, gas density (e.g. Leroy et al., 2008). This suggests that the surface density of total gas required for star formation could plausibly be lower than the Kauffmann (1996) value.

Equation (5) can be converted to a total critical mass, , by integrating out to the disc radius, :


Following Croton et al. (2006), we assume . The factor of 3 was chosen by Croton et al. (2006) based on the properties of the Milky Way and therefore may not be representative of the high-redshift galaxies which we consider in this work. However, instead of considering this to be another free parameter of the model we note that , meaning any such freedom can already be considered to be included in . In future work, we will compare our predicted disc sizes to high- observations and investigate the success of these simple scaling relations in detail.

If the total amount of cold gas in the disc is greater than the critical mass, the star formation rate is assumed to be given by


where is the dynamical time of the disc and is a free parameter describing the efficiency of star formation in the form of the formation time-scale in units of the dynamical time.

In summary, the star formation prescription we employ in Meraxes is almost identical to that of Croton et al. (2006) with the addition of as an extra free parameter (as previously advocated by Henriques et al., 2015).

2.5 Supernova feedback

The radiative and mechanical energy liberated by supernovæ can have a profound impact on galaxy evolution, potentially heating significant amounts of gas and even ejecting it from a galaxy or host dark matter halo entirely. This is especially so at high redshift where haloes are on average less massive than in the local Universe, and possess correspondingly shallower potential wells. Supernovæ also enrich the inter stellar medium (ISM), altering the chemical composition of future stellar generations and changing the efficiency with which gas can cool (c.f. Equation 2).

Delayed supernova feedback

Many semi-analytic models make the simplifying assumption that all supernova feedback energy is released instantaneously, during the same snapshot in which the relevant stars formed. This approximation is motivated by the reasonable further assumption that the majority of supernova feedback energy is released by massive stars () which have short (; Portinari et al., 1998) lifetimes, ending in violent SN-II. In the cases where the time span between each simulation snapshot is large (e.g. 250 Myr in the case of the Millennium Simulation; Springel et al., 2005), the approximation of the instantaneous deposition of all supernova energy into the ISM is valid. However, motivated by the short dynamical time of systems at high redshift, the separation between snapshots in our input simulation is approximately 11.1 Myr. It therefore takes at least three snapshots after a single coeval star formation episode for all stars more massive than to have gone supernova. In order to accommodate this matching of time-scales in Meraxes, we have implemented a simple delayed supernova feedback scheme which we outline in this section.

We begin with the standard assumption that all supernova feedback energy is released by SN-II which are the end result of the evolution of stars with initial masses greater than . The basic methodology of our delayed feedback scheme is then to calculate the total amount of energy which should be injected into the ISM by a single star formation episode, and to release this energy gradually over time in proportion to the fraction of SN-II which will have occurred.

We assume a standard Salpeter (1955) initial mass function (IMF) with upper and lower mass limits of and respectively:


where, by definition


and thus = 0.1706. With this choice of IMF, the number fraction of stars that will end their lives as type II supernovæ () is given by:


If we further assume that each supernova produced injects a constant of energy into the ISM then, for a burst of mass , the total amount of energy deposited into the ISM () is


where is a free parameter describing the efficiency with which the supernova energy couples to the surrounding gas. As is common practice, we model the mass of gas which is reheated by this energy deposition () as


where is a free parameter commonly referred to as the mass loading factor.

Croton et al. (2006) used constant values of and for all galaxies. However, we find that we are unable to replicate the observed shallow low-mass slope of the stellar mass function at without adopting a value for these parameters that scale with mass. We therefore follow Guo et al. (2013) who encountered a similar issue (although at lower redshifts) leading them to adopt the following parametrizations for :


and similarly for :


where , , , , , and are all free parameters. Since corresponds to an efficiency, we enforce that it must take a value in the range 0–1 at all times. We have also imposed the additional constraint that the mass loading factor, , cannot exceed an upper limit which we nominally set to be based on reasonable expectations for typical high- galaxies (e.g. Martin, 1999; Uhlig et al., 2012). For a standard energy-driven wind, (Murray et al., 2005). However, the value of is far less certain and depends on the poorly understood efficiency with which injected supernova energy is thermalized (Murray et al., 2005) and potential variations in the IMF of stars.

As discussed above, it takes approximately 40Myr for an star to go supernova (Portinari et al., 1998). As a result, the total amount of supernova energy released by a galaxy at any given snapshot, , will be dependent on the mass of stars formed both in the current and previous snapshots. We therefore explicitly track the total mass of stars formed in each galaxy (and all of its progenitors) for the last snapshots4. The value of is dependent on the input N-body simulation and is chosen such that at least the last of star formation is recorded at all times. For Tiamat this corresponds to . At snapshot , the value of is then


Similarly, the amount of cold gas reheated by this energy is


The term in the two equations above denotes the fraction of stars formed during snapshot , that go supernova during snapshot . This can be calculated by integrating the stellar IMF () between suitably chosen mass limits:


The values of and are set by the range of stellar masses formed during snapshot which will have had time to expend their fuel and go nova during the time spanned by snapshot . To calculate these we use a functional form fit to the H and He core burning lifetimes tabulated by Portinari et al. (1998), under the assumption that all stars go supernova immediately upon expending their H and He cores,


where and is the time since the stars formed. This fit is accurate to within 6% for all values of appropriate for this work. For simplicity, we also approximate the star formation which occurs during any given snapshot by a single coeval burst at the middle of that snapshot. Although crude, this approximation results in an error of at most in which rapidly becomes negligible as the time since the burst increases (i.e. increases). Finally we note that since we are assuming that all supernova feedback energy is produced by SN-II, we enforce a minimum value of 8 when evaluating Equations 15 & 16 above.

The eventual fate of the reheated material depends on both its mass, , and the amount of energy injected, . If we assume the gas to be adiabatically heated to the virial temperature of its host halo, the associated change in thermal energy is given by:


If then there is more energy injected into the reheated gas than is required to raise it to the virial temperature of the halo. We therefore assume the gas to be added to the hot halo of the host FoF group. Any excess energy is assumed to then go into ejecting some fraction of the FoF group hot reservoir from the system entirely:


where is the virial velocity of the FoF group. If instead , only an energetically feasible fraction of the total reheated mass is added to the FoF group hot reservoir:


with the rest raining back down on to the galaxy in a galactic fountain.

Any gas and metals which are successfully expelled from the system entirely are placed into a separate ‘ejected’ reservoir. Here they are assumed to play no further role in the evolution of the galaxies in the host FoF group until the group falls into a more massive system. At this point, the ejected material is assumed to be re-accreted into the new group and is added to its hot halo component.

Delayed versus contemporaneous feedback

In practice, we apply our supernova scheme in two phases. First, the amount of mass reheated and ejected due to past star formation episodes is calculated as described above. After the masses of the various baryonic reserves have been appropriately updated, the amount of new star formation in the current snapshot is then determined (cf. Section 2.4). In this way, ongoing energy injection from past star formation episodes is able to prevent new stars from forming in the current time step altogether.

After calculating the mass of stars formed, the corresponding reheated and ejected masses due to any stars with short enough lifetimes to go nova in the current time step are also calculated. If the total amount of cold gas removed from the galaxy due to both star formation and the corresponding contemporaneous supernova feedback exceeds that which is available, the mass of stars formed in the current time step is reduced until consistency is achieved.

2.6 Metal enrichment

As is common in semi-analytic models (e.g. De Lucia et al., 2004; Somerville et al., 2008; Guo et al., 2011), we implement a simple metal enrichment scheme whereby a fixed yield, , of metals is released into the ISM per unit mass of stars formed. Again, we assume that these metals are released predominantly by massive stars which end their lives as SN-II and we gradually release them over time as these supernovæ occur. However, since a more massive star will generally release more metals during its lifetime than a less massive counterpart, we release these metals in proportion to the mass fraction of SN-II (as opposed to the number fraction as was used above). In other words


where is the mass of metals released during snapshot , is the total fraction of stars with initial masses greater than 8, and is the fraction of stars formed during snapshot that go nova during snapshot .

Analogous to Equations 10 and 17 above, and are given by


where and are as defined in Section 2.5.

The metals released into the ISM in this manner are assumed to be uniformly mixed with the cold gas of the galaxy. From here they can be further distributed to the hot halo or ejected from the system entirely via supernova feedback. Metals which do end up in the hot gas reservoir can then enhance the cooling rate of gas on to the galaxy through metal line emission (see Equation 2).

2.7 Stellar mass recycling

A common assumption of many semi-analytic models is the so-called instantaneous recycling approximation (IRA), in which some fixed fraction of the stellar mass formed during each time step is instantaneously recycled back into the ISM. The precise value of this fraction varies from model to model and is often left as a free parameter; however, most works employ a value of approximately 30–40% (e.g. Cole et al., 2000; Croton et al., 2006; Somerville et al., 2008; Henriques et al., 2015).

Given our choice of IMF (see Section 2.5), a recycle fraction of 40% corresponds to all stars more massive than approximately instantaneously going supernova. However, the lifetime of a star is close to the current age of the Universe (e.g. Portinari et al., 1998) and hence the IRA can only be considered valid for galaxies around (i.e. well past the peak of the Universal star formation rate density). At , where the majority of galaxies have stellar populations dominated by recent star formation, this approximation becomes invalid and we are forced to consider a more realistic alternative.

Our stellar mass recycling prescription is divided into two parts:

  1. When implementing our delayed supernova feedback prescription, we assume that the initial stellar mass of all supernovæ is returned to the cold gas reservoir of the galaxy (thus ignoring any mass that may be locked up in long lived remnants such as neutron stars and black holes).

  2. As discussed in Section 2.5, we explicitly track the star formation history of each galaxy for the last snapshots. In order to calculate the recycled mass from older stars we approximate these as having formed in a single coeval burst which occurred at a time defined by their mass-weighted age. Equation 24 then allows us to calculate the relevant mass of stars which would have gone nova and this is again assumed to be returned the cold ISM in its entirety. Although crude, the approximation of a single coeval burst for all older stars provides the correct stellar mass loss to within less than 5% error at all times5 given the snapshot cadence of Tiamat and our fiducial value of .

2.8 Halo infall and gas stripping

As haloes inspiral towards more massive systems, tidal forces experienced during repeated pericentric passages can lead to the stripping of loosely bound material from the outer regions. In Meraxes, if the mass of an FoF group drops, then a pro rata fraction of the ejected and/or hot baryonic content of the halo is also removed. In practice, the amount of mass which must be removed is given by the value of (as defined in Equation 1) which will be negative in such systems. This material is taken first from the ejected reservoir, with further mass being removed from the hot halo component if required. No baryons are ever taken from the cold gas or stellar mass reservoirs as these are assumed to be protected from such tidal losses by their position deep in the central potential well of their haloes.

Further to these long-range tidal forces, galaxies infalling into groups or clusters are observed to be subjected to a number of dynamical processes which remove gas from the outskirts of the system, including ram-pressure stripping, strangulation, and harassment (e.g. van den Bosch et al., 2008; Peng et al., 2015). We model the combined effects of these processes by assuming that all FoF groups are instantly stripped of their entire hot and ejected gas reservoirs upon infall into a more massive structure, with their combined mass and metals being added to the hot component of the new parent. Although such a rapid stripping represents the most extreme scenario possible, we note that this approximation has been made in a number of previous semi-analytic models and defer an improved and more realistic treatment to future works.

2.9 Mergers

Mergers play an important role in the build up of galaxy stellar mass, both through hierarchical mass assembly and induced star formation. This is particularly so at high- where their prevalence is enhanced (cf. Paper I). In Meraxes (as in almost all semi-analytic models) these galaxy merger events are triggered by the merging of the corresponding host dark matter haloes. Following Croton et al. (2006), when a dark matter halo is marked as having merged, we utilize dynamical friction arguments to approximate the time taken for the orbit of the incoming galaxy to decay and the corresponding galaxy–galaxy merger to occur (Binney & Tremaine, 2008):


where it is standard to take , is the distance between the most-bound particle of the parent and the infalling halo, is the total mass of the infalling galaxy, and and are the virial properties of the parent. In Meraxes, all of these quantities are evaluated at the last time the infalling halo was successfully identified in the trees.

The value of is based on the assumption that is calculated at the moment the infalling halo crosses the virial radius of the parent. However, we instead calculate at the time at which the infalling halo can no longer be identified in the Tiamat simulation and is thus marked as having merged in our input merger trees. This results in an overestimate of the merger time-scale which worsens the longer the infalling halo remains identified after crossing the virial radius of the parent (e.g. Hopkins et al., 2010). Even in dense environments, the accurate merger trees produced from Tiamat results in haloes being identified for extended periods before the merger event occurs. Furthermore, previous authors have found that changes to the value of have been necessary in order to match observational constraints on the luminous end of the galaxy luminosity function (De Lucia & Blaizot, 2007) and idealized N-body halo merger simulations (e.g. Boylan-Kolchin et al., 2008). On average, we find that infalling haloes are successfully tracked in our input merger trees until , with a weak trend to be identified to smaller fractional radii with increasing redshift. Noting that , we therefore choose to fix . We also note that if, after starting the merger clock, the parent galaxy itself experiences a merger, then we assume that all of its infalling galaxies also undergo a merger with the same target.

Galaxy mergers can drive strong shocks and turbulence in any participating cold gas, driving this material towards the inner regions of the parent galaxy and resulting in an efficient burst of star formation. We model the fraction of cold gas consumed by such a burst, , using the prescription introduced by Somerville et al. (2001):


where and are the corresponding baryonic masses (i.e. cold gas + stellar mass), and we follow Croton et al. (2006) by setting the parameters and . This relation agrees well with the results of numerical simulations of mergers with baryonic mass ratios in the range 0.1–1.0 (Cox et al., 2004). For merger events where the mass ratio is less than 0.1, we suppress any merger-driven star formation.

For simplicity, we assume that all of the stars formed in a merger-driven burst do so within a single snapshot. At , the median dynamical time of a galaxy disc in Meraxes is 60% of the time between two consecutive snapshots of Tiamat (11.2 Myr). Hence, this approximation is roughly equivalent to the assumption that the merger-driven burst occurs on a time-scale approximately less than one disc dynamical time for the majority of galaxies. Although the disc dynamical time does increase with decreasing redshift, by the median is still only equal to one snapshot and hence this approximation remains valid for the majority of galaxies.

2.10 Ghost galaxy evolution

A ghost dark matter halo is one which is temporarily unresolved in our input merger trees. This can be due to a number of reasons, but is most commonly a result of a smaller halo passing through or nearby a much more massive structure. The Tiamat merger trees used in this work are carefully constructed to identify these artefacts, resulting in the skipping of a potentially large number of snapshots between haloes and their descendants. In many semi-analytic models, the galaxies hosted by such haloes are simply ignored until their halo is later re-identified. In some cases, re-identification fails or is not even attempted, resulting in spurious galaxy merger and creation events. However, we must ensure that we correctly include these objects at all snapshots in order to account for their ionizing photon contribution.

Due to the lack of knowledge of the properties of a ghost’s host dark matter halo, we are unable to implement many of the physics prescriptions outlined above. We therefore simply allow these galaxies to passively evolve during the time over which they are identified as ghosts, forming no new stars but experiencing the delayed supernova feedback from previously formed generations. When the host halo is eventually re-identified, we assume that any associated star formation occurred in a single coeval burst at , where is the time since the halo was last identified. We then back-fill the stellar mass history appropriately for use in our delayed supernova feedback scheme.

2.11 Reionization

A key goal of DRAGONS is to connect the evolution of the reionization structure to the formation of the source galaxy population. As such, Meraxes has been developed to be the first semi-analytic galaxy formation model to fully and self-consistently couple the process of reionization (in particular, the presence of a photo dissociating UVB) to the evolution of galaxies both temporally and spatially. To achieve this we have embedded a specially modified version of the semi-numerical reionization code, 21cmFAST (Mesinger et al., 2011), that includes the calculation of the local ionizing UVB described by Sobacchi & Mesinger (2013) and, importantly, makes full use of the realistic galaxy properties provided by Meraxes.

Self-consistent reionization with 21cmFAST

The basic methodology of 21cmFAST is to use an excursion set formalism in order to identify ionized bubbles where the integrated number of ionizing photons is greater than the number of absorbing atoms and associated recombinations:


Here is the integrated number of atoms being ionized within a sphere of radius , is the number of stellar baryons in the same volume, is the mean number of ionizing photons produced per stellar baryon, is the escape fraction of these photons, and is the mean number of recombinations per baryon. If we assume that helium is singly ionized at the same rate as cosmic hydrogen, then expanding Equation 27 in terms of the integrated stellar () and total () masses within gives


where is the helium mass fraction, is the proton mass, and the term corresponds to the combined number of hydrogen and helium atoms per baryon.

It is common for Equation 28 to be re-written in terms of an H ii ionizing efficiency, :




and we have excluded the term based on studies of the high-redshift Lyman- forest which suggest in the diffuse IGM (e.g. Bolton & Haehnelt, 2007; McQuinn et al., 2011). Despite this simplifying assumption of , we note that we implicitly include a mean-free path of ionizing photons through the IGM in our calculation by starting our excursion set calculation at an appropriate scale (Sobacchi & Mesinger, 2013). The right-hand side of this equation includes our fiducial values for each of the physical variables. The number of ionizing photons per stellar baryon, , is set by the assumed stellar IMF whilst both and are well constrained by cosmology. Of all of these terms, only is poorly known for high-redshift galaxies (Wise & Cen, 2009; Raičević et al., 2011; Kuhlen & Faucher-Giguère, 2012, e.g.). Our fiducial value of 0.2 is primarily chosen to provide a reionization history which is consistent with the latest Planck 2015 electron scattering optical depth measurements (Planck Collaboration, 2015, see Section 3 below).

By applying our integrated 21cmFAST algorithm to grids of stellar and total mass within Meraxes we can use Equation 29 to produce a neutral hydrogen fraction () grid for the entire simulation volume. In order to then determine how this spatially and temporally evolving ionization structure affects the baryon fraction modifier, , of individual FoF groups we utilize the UVB feedback model of Sobacchi & Mesinger (2013). Using idealized 1D hydrodynamical simulations of a static, uniform ionizing UVB impinging on collapsing dark matter haloes, Sobacchi & Mesinger (2013) found that was well described by


where is the mass of the halo and is the ‘filtering mass’ representing the mass at which :


Here is the redshift at which the collapsing halo was first exposed to the UVB and the parameters were found by the authors to provide the best fit to their results. The term represents the local UVB intensity:


where for a stellar-driven UV spectrum (Thoul & Weinberg, 1996).

In order to calculate using this formalism, we need to know both the redshift at which the IGM surrounding each halo was first ionized, , and the local ionizing background intensity at this time, . The average UVB intensity which a galaxy is exposed to within an ionized region, , is given by


where is the comoving mean-free path of ionizing photons (which is assumed here to be equal to the radius of the ionized bubble, ) and is the Planck constant. The term is introduced to account for the effect of galaxy clustering on boosting the ionizing emissivity at halo locations relative to the spatial average (Mesinger & Dijkstra, 2008). The term represents the ionizing emissivity, the time-averaged number of ionizing photons emitted into the IGM per unit time, per unit comoving volume. Approximating the average rate of stellar mass growth of all galaxies within this region as , where is the gross stellar mass formed within (i.e. without any decrement due to stellar evolution) and is the Hubble time, can be expressed as


Our utilization of , instead of the true instantaneous star formation rate predicted by Meraxes, is motivated by the need to smooth out the often bursty star formation histories of our galaxies. The filtering mass formula of Equation 32 assumes that the ionizing background intensity within a cell remains constant with time. This is a reasonable approximation due to the weak sensitivity of on the UVB intensity (). Sobacchi & Mesinger (2013) also find that the remains approximately constant within H ii regions, further validating this approximation. However, fixing the UVB intensity at the snapshot of ionization results in smaller bubbles having artificially high or low values depending on the current star forming state of the source galaxies. This introduces artificial scatter into the calculation which can become important when comparing the effect of reionization on the evolution of individual galaxies (see Section 4.5).

In practice, the coupling between galaxy evolution and reionization within Meraxes is implemented as follows.

  1. For a single time step in the simulation, Meraxes evolves all of the galaxies in the entire Tiamat volume.

  2. Our 21cmFAST algorithm constructs and processes halo mass, stellar mass, and averaged star formation rate grids, along with pre-computed total matter density grids. For this work, and all of the results herein, we use a grid resolution of . After applying the excursion set formalism and equations outlined above, a grid of and values is generated.

  3. Using the grid, Meraxes keeps track of the redshift at which each cell first became ionized () and the corresponding . It then calculates the baryon fraction modifier of each grid cell following Equation 31.

  4. In order to calculate the amount of freshly infalling baryonic material it should accrete, each FoF group uses the baryon fraction modifier of the grid cell in which it is located in the following simulation time step (see Section 2.2 above).

  5. This process is then repeated for each of the time steps (of which there are 100 for our input simulation, Tiamat, between ).

Through this procedure, the evolution of galaxies and reionization in the simulation are self-consistently coupled both temporally and spatially.

Homogeneous model

For comparison, we also employ a ‘homogeneous’ reionization prescription using values that depend on redshift alone (i.e. no information about the spatial distribution of the IGM ionization state is required). Prescriptions such as this (e.g. Gnedin, 2000; Kravtsov et al., 2004) have been commonly employed by semi-analytic models for many years (e.g. Benson et al., 2002; De Lucia & Blaizot, 2007; Somerville et al., 2008).

Figure 1: The mean filtering mass, , of FoF groups in the fiducial patchy reionization model as a function of redshift. For our homogeneous comparison model, these filtering mass values are applied to every FoF group in the simulation, regardless of their environment or local reionization history, thus allowing us to directly assess the effect of ignoring this information as is commonly done in traditional semi-analytic galaxy formation models. The blue shaded region indicates the 68% confidence intervals calculated using the spatial variation in the mean from 125 non-overlapping subvolumes comprising the full simulation.

At each time step in the simulation, we use our fiducial patchy reionization prescription described above to calculate the number-weighted mean filtering mass of all FoF groups in the volume. We then re-run Meraxes, applying this redshift dependent value to all FoF groups when calculating their baryon fraction modifier as per Equation 31. By generating from our fiducial model in this way, we are able to use the homogeneous model result as a baseline with which to assess the detailed effects of a self-consistent, spatially dependent reionization prescription on the evolution of the source galaxy population, whilst simultaneously providing a useful filtering mass relation for , that can be applied quickly and easily. In Fig. 1, we show the evolution of . The blue shaded region indicates the 68% confidence intervals calculated using the spatial variation of the mean in 125 subvolumes. The evolution in the filtering mass is well approximated by the following functional form:


where , and . The fitted parameter values were calculated using Markov chain Monte Carlo (MCMC) methods with a standard chi-squared likelihood, and provide a fit which is accurate to within of the mean model result across all redshifts. As we will demonstrate in Section 4, our homogeneous prescription does a reasonable job of reproducing the mean evolution of the stellar mass functions and global neutral fractions predicted by the full, patchy reionization implementation.

3 Model calibration

Parameter Prescription Equation Description Munich model
values value
Star formation (§2.4) 5 Critical cold gas surface density normalization    0.26 0.38 0.2
7 Star formation efficiency 0.01 0.07 0.03
Supernova feedback (§2.5) 13 Energy coupling efficiency normalization 0.18 0.7 0.5
13 Coupling efficiency scaling 0 3.5 2.0
13 Coupling efficiency normalization 70 336 70.0
14 Mass loading normalization 2.1 10.3 6.0
14 Mass loading scaling 0 3.5 0.0
14 Mass loading normalization 70 430 70.0
14 Maximum mass loading value 10.0
Metal enrichment (§2.6) 22 Mass of metals per unit mass of SN 0.03 0.047 0.03
Reionization (§2.11) 30 Ionizing photon escape fraction 0.2
Table 1: The fiducial parameter values used throughout this work. Values were constrained to visually reproduce the observed evolution in the galaxy stellar mass function between and 7 (see Fig. 2). The quoted Munich model values represent the range of fiducial values utilized in the following works (where appropriate): Croton et al. (2006), Guo et al. (2011, 2013), Mutch et al. (2013), and Henriques et al. (2013, 2015). All of the parameters in these works were calibrated against observations. However, they are presented here as a rough guide to the range of plausible values.

The free parameters of the model were manually calibrated (by hand) to replicate the observed evolution of the galaxy stellar mass function between redshifts 5 and 7, as well as the integrated free electron Thomson scattering optical depth measurements. The evolution of the stellar mass function has been shown by previous statistical investigations of semi-analytic models to provide a tight constraint on both the star formation efficiency and supernova feedback parameters (e.g. Henriques et al., 2013; Mutch et al., 2013). By combining this with the Thomson scattering optical observations, we can additionally put constraints on the escape fraction of ionizing photons (), and thus all of the free parameters of our model, as listed in Table 1.

It could be argued that the luminosity function would provide a more fundamental constraint on the model, rather the stellar mass function. However, whilst it is true that converting observed galaxy luminosities to stellar masses involves a number of assumptions and potentially unreliable conversions, the same is also true for the inverse procedure of converting model stellar masses to luminosities. Stellar masses are intrinsic predictions of semi-analytic models, whilst luminosities require an extra layer of modelling. For example, changing the IMF has a relatively small impact on the stellar masses predicted by the model, but can have a significant impact on the resulting luminosity function. In order to accurately model luminosities in various bands, one must typically calculate a full spectral energy distribution (SED) for every object, applying model-dependent Lyman absorption, sample selection (e.g. in colour–colour space), and poorly understood dust corrections. Furthermore, doing this for all galaxies at multiple redshifts can take a significant amount of time and memory which prohibits its usefulness when running a model many times for calibration purposes.

The redshift range corresponds to the highest redshifts for which reliable observed stellar mass functions are available. In particular, we make use of the mass functions estimated by González et al. (2011), Duncan et al. (2014), Grazian et al. (2015), and Song et al. (2016). Both Duncan et al. (2014) and Grazian et al. (2015) utilize data collected from the Cosmic Assembly Near-infrared Deep ExtragaLactic Survey (CANDELS; Grogin et al., 2011; Koekemoer et al., 2011) GOODS South field with stellar masses directly obtained from SED fitting of combined optical and near-infrared space-based observations, and include the effects of both nebular line and continuum emission. In addition, Grazian et al. (2015) include a detailed treatment of the effects of Eddington bias (Eddington, 1913) on the normalization and slope of their derived mass functions. Whilst Song et al. (2016) also utilize CANDELS infrared data, they instead carry out a hybrid SED stacking technique to derive a redshift dependent stellar mass–UV luminosity relation which is then combined with measured UV luminosity functions to estimate the galaxy stellar mass function. Similarly, González et al. (2011) utilized data combined from Hubble Space Telescope and Spitzer observations, but with stellar masses obtained via mass–UV luminosity relations calibrated at alone.

In addition to the stellar mass function which constrains the integrated amount of star formation driving reionization, the corresponding ionizing photon budget is also set by the escape fraction. One of the primary observational constraints on the timing and duration of reionization comes from the measured integrated optical depth to Thomson scattering of cosmic microwave background photons by free electrons, :


where is the Thomson scattering cross-section, is the mass-weighted global ionized fraction of species , and and are the average comoving density of hydrogen and helium, respectively (Wyithe & Loeb, 2003). For this work, we have assumed that helium is singly ionized at the same rate as hydrogen (i.e. ) and only becomes doubly ionized at (i.e. or 0 for greater than or less than 4 respectively; e.g. Kuhlen & Faucher-Giguère, 2012). As we shall demonstrate in the following sections, primarily constrains the escape fraction of ionizing photons, , in the model.

Figure 2: The evolution of the galaxy stellar mass function from 5–7. Coupled with the measured Thomson scattering optical depth (see Fig. 4), these are the only observational constraints applied to Meraxes throughout this work unless explicitly stated. Data points show the observations of Duncan et al. (2014), González et al. (2011) and Song et al. (2016). The purple dashed lines show the best-fitting Schechter functions of Grazian et al. (2015). Solid blue lines show the self-consistent patchy reionization model using our fiducial parameter values. Best fit low-mass slopes are provided in Table 2. All observations have been corrected to a Salpeter IMF where necessary.

The resulting parameter values for our fiducial model constrained against both the high- stellar mass function evolution and integrated electron scattering optical depth are presented in Table 1. It is important to note that, although these parameter values provide a good match to the constraining observations and are broadly consistent with comparable works where appropriate (see second-to-last column), they may not be the only possible solution. However, through testing extreme (to the point of being physically implausible) parameter combinations we are able to ascertain that supernova feedback is the only feedback mechanism in our model capable of producing a stellar mass function with a slope consistent with observations. Regardless, we impress upon the reader that all of the results in this work must be interpreted within the context of these particular chosen parameter values alone. In future work we will carry out a full MCMC analysis to accurately constrain the free model parameters against a wider range of observational quantities as well as explore any degeneracies which may exist between them (e.g. Lu et al., 2011b; Mutch et al., 2013).

In Fig. 2 we present the fiducial model stellar mass functions (blue solid lines) along with the constraining observations. All observations have been converted to a Salpeter IMF and where necessary. The error bars on observed data points include contributions from Poisson noise and uncertainties in photometric redshift determinations. However, they neglect the systematic uncertainties associated with the estimation of stellar masses from photometric data (e.g. stellar population synthesis model variations, and photometric uncertainties).

We are able to achieve an excellent match to the normalization, shape, and evolution of the observed mass function across all plotted redshifts. At , where there is the largest divergence between different observational data sets, we chose parameter values which provided a reasonable compromise between each. However, at the low-mass end we have chosen to follow the observations of Duncan et al. (2014) as they use a large data set with stellar masses obtained from SED fitting and provide actual data points rather than a Schechter fit. The quality of the agreement between our model and the observational data gives us faith that our implemented physical prescriptions are both reasonable and applicable at the high redshifts of interest in this work.

Although typically producing fewer ionizing photons than their more massive counterparts, low mass galaxies with are expected to contribute a large fraction of the overall ionizing photon budget due to their high number density. For this reason, the low-mass slope of the stellar mass function, , is of particular importance to reionization. In Table 2, we provide the Meraxes best-fit parameters at each redshift, obtained by fitting a standard Schechter function to the model results using MCMC methods6 and are in good statistical agreement with the corresponding values found by Duncan et al. (2014) of , and for redshifts 5, 6 and 7 respectively.

Figure 3: Upper: the distribution of total stellar mass as a function of FoF group virial mass for the fiducial model at . The solid and dashed black lines show the median and 68% confidence intervals of the distribution. The grey shaded region indicates halo masses below the atomic cooling mass threshold. Lower: the evolution of the fiducial model median relation between redshifts 10–5. There is no significant evolution in either the slope or normalization with redshift. The grey dashed line indicates the theoretically motivated slope of () suggested by Wyithe & Loeb (2013) for supernova regulated galaxy growth. The normalization of this line has been arbitrarily chosen to allow a comparison with the slope of 1.4 predicted by Meraxes.

In the top panel of Fig. 3 we present the distribution of total stellar mass (i.e. summed over all galaxies) in each FoF group as a function of the group virial mass in our fiducial model. The median relation (solid black line) is well described by a power law with a slope of . This shows good agreement with simple energy conservation arguments which suggest a slope of for supernova feedback-regulated galaxy growth and a fixed cold gas mass fraction (Wyithe & Loeb, 2013). However, at FoF group masses below there is a rapid increase in the spread of stellar mass values. This is due to a combination of supernova and reionization feedback effects, as well as a low star formation efficiency in these small, often diffuse haloes. In the lower panel of Fig. 3 we show the evolution of the median FoF group relation as a function of redshift. Interestingly, there is no statistically significant evolution in either the slope or normalization of the relation with redshift. However, the same simple energy conservation arguments which provided a good agreement for the slope of the relation suggest that the normalization should evolve with redshift (Wyithe & Loeb, 2003). Despite this, our model indicates that in order to reproduce the observed evolution of the galaxy stellar mass function over the redshifts considered in this work, the efficiency of galaxy formation and the associated feedback processes must conspire to provide a constant star formation efficieny in haloes of a fixed mass. This agrees with similar findings from subhalo abundance-matching (SHAM) studies at lower redshifts (e.g. Behroozi et al., 2013a).

Figure 4: The integrated free electron scattering optical depth, , as a function of redshift. The grey horizontal line and shaded region indicate the constraints on to from the Planck 2015 data release (Planck Collaboration, 2015). The blue solid line shows the fiducial model which is constrained to reproduce the Planck result. The homogeneous model result is obscured by the fiducial line.

In Fig. 4, we present the electron scattering optical depth of our fiducial model (blue line) against the current best observational measurements provided by the Planck satellite (Planck Collaboration, 2015). The other model variations shown in this plot will be discussed in detail below. However, as can be seen, our fiducial model provides an integrated optical depth which is in excellent agreement with the Planck results.

An important consideration when assessing the results of any cosmological simulation is the potential loss of stellar mass (and therefore ionizing photon contribution) due to finite mass resolution. This issue is explored in detail for Meraxes and the Tiamat suite of simulations in Appendix A. In summary, we find that the halo mass function of the full Tiamat simulation is complete down to approximately above the atomic cooling mass threshold at , whilst the lower volume but higher mass resolution Medi Tiamat is fully complete down to this mass limit. Using Medi Tiamat to quantify the fraction of total stellar mass missing from the full simulation, we find that at we miss approximately 25% of stellar mass from low mass, unresolved systems. By this fraction falls to .

Figure 5: The ionization structure of a thick slab produced by Meraxes. The left- and right-hand columns correspond to a volume-weighted global neutral fraction of and 0.3, respectively. Top row: the ionization state of the slab. Orange indicates regions of the volume which are neutral, whilst the green structure shows the underlying matter distribution inside the ionized bubbles surrounding the densest structures. At (right), there is a large ‘lane’ of ionized IGM approximately wide extending from top to bottom where multiple bubbles overlap. Middle row: the same ionization structure with a random 1/50th of the ionizing galaxy population (selected at ) overlaid as points scaled by stellar mass. Bottom row: the same galaxy population as above but with point size scaled by instantaneous star formation rate. There are fewer galaxies in this row, compared to above, due to the bursty nature of star formation resulting in some galaxies forming no new stars in the current snapshot.

Finally, in Fig. 5, we present a deep slab extracted from the fiducial Meraxes model (see Section 4 below) at a volume-averaged neutral fraction of 0.8 (left) and 0.3 (right). The orange regions indicate the neutral portions of the simulation volume, whilst the green structure shows the underlying matter distribution within the ionized bubbles. At early times, these ionized bubbles surround the peaks in the density field where the first galaxies form. As reionization progresses the bubbles begin to overlap, and by a large ‘lane’ of ionized hydrogen, extending the entire length of the simulation volume is formed. The ability to investigate the distribution and evolution of bubble morphologies and the associated observable power spectra is a key feature of our Meraxes framework and is investigated in detail in Paper V. In the middle and bottom rows of Fig. 5, we overplot the positions of a random 1/50th of the full galaxy population in the slab, selected at . The galaxies plotted in the panels are the main progenitors of this subsampled population. In the middle row, the points are scaled by the stellar mass of each galaxy, whilst in the bottom row they are scaled by their instantaneous star formation rates. The fact that fewer points are present in the star formation rate panels is a result of the bursty nature of star formation in the model meaning some fraction of galaxies have no ongoing star formation at any particular redshift.

4 The interplay between reionization and galaxy growth

In this section, we use Meraxes to investigate the relative importance of reionization feedback for regulating both the growth of galaxy stellar mass and the timing and duration of the EoR. We also quantify how simple variations to the physics of reionization affect the galaxy stellar mass function, before finally elucidating the importance of local environment for determining the stellar mass of galaxies affected by photoionization suppression. We again note that the results of this work should be interpreted within the context of our chosen parameter values, constrained to match the evolution of the galaxy stellar mass function and the optical depth to electron scattering as described in Section 3.

Fiducial 0.03 0.2 7.96 1.47
No reionization feedback 0.03 0.2 8.08 1.42
No SN feedback 0.03 0.2 11.41 1.52
No feedback 0.03 0.2 11.50 1.47
Recalibrated no SN feedback 0.00106 0.2 7.75 1.29
Recalibrated no feedback 0.00106 0.2 7.8 1.26
Half 0.03 0.1 6.89 1.38
Double 0.03 0.4 9.04 1.50
Homogeneous 0.03 0.2 7.94 1.50
Table 2: Summary table of the different model runs explored in this work (see Section 4). indicates the redshift spanned between a global neutral fraction of 80% and 20%. The and columns indicate the star formation efficiency and ionizing photon escape fraction of each run respectively. The , 6, and 7 SMF parameters (last three columns) refer to the low mass slope of the corresponding stellar mass functions fit with a standard Schechter function. For reference, the parameters measured by Duncan et al. (2014) are , and for redshifts 5, 6, and 7, respectively.

Throughout we focus on the following model variations.

Fiducial ():

both the spatial and temporal evolution of the reionization structure and UVB intensity are fully coupled to the growth of the source galaxy population. This is the model calibrated in Section 3 above to match the evolution of stellar mass function and Thompson scattering optical depth measurements.

No feedback ():

the fiducial model with no supernova feedback or reionization feedback included, resulting in runaway star formation.

No supernova feedback ():

the fiducial model with supernova feedback turned off. All remaining physical processes, including reionization feedback, remain unchanged.

No reionization feedback ():

the fiducial model with reionization feedback removed by setting for all galaxies at all times. All remaining physical processes, including supernova feedback, remain unchanged.

Recalibrated no supernova feedback ():

the no supernova feedback model with a lower star formation efficiency, , chosen to replicate the total =5 stellar mass density of the fiducial model.

Recalibrated no feedback ():

the no feedback model with the same lower star formation efficiency ().

Half ():

the fully coupled fiducial model but with a lower escape fraction of .

Double ():

the fully coupled fiducial model but with a higher escape fraction of .

Homogeneous ():

the evolution of reionization is decoupled from that of the growth of galaxies. The baryon fraction modifier of each halo, , is determined using the halo number-weighted average as a function of redshift calculated using the fiducial model (see Section 2.11.2 above). All remaining physical processes remain unchanged from the fiducial case.

4.1 The relative importance of reionization feedback

In Section 3 above, we demonstrated that our fiducial model successfully reproduces the evolution of the high- galaxy stellar mass function, as well as the most recent electron scattering optical depth measurements. In this section, we utilize the resulting realistic population of galaxies to investigate how important photoionization suppression of baryonic infall is for regulating the stellar mass content of dark matter haloes when compared to galactic feedback processes such as supernova feedback.

The stellar mass function

Figure 6: The and 8 galaxy stellar mass functions. The solid blue line indicates the result of the fiducial patchy reionization model which has been calibrated to reproduce the observed =5–7 mass functions. Grey points in the left-hand panel show the relevant observational data. The small difference introduced by omitting reionization feedback (red dash–dotted) demonstrates the minor role which this mechanism plays in regulating stellar mass growth. Conversely, the removal of supernova feedback (gold solid) produces a much larger effect. However, comparison between the no feedback (gold dashed) and no supernova feedback (grey solid) models shows that, in the absence of the dominant supernova feedback, reionization does suppress the number density of galaxies if star formation is gas supply limited in a large fraction of haloes.

In Fig. 6 we show the (left) and (right) galaxy stellar mass functions from each of our model variations. The latter redshift value corresponds to a volume-averaged neutral hydrogen fraction of 50% in the fiducial model (; cf. Fig. 9 and Table 2). The bottom panels also indicate the fractional differences of a subset of the models with respect to the fiducial result. Immediately apparent is that models without supernova feedback (, , ,) produce the most significant change to the stellar mass function and are the only models which are not consistent with the observational data (grey points). The gold solid line () shows the predicted stellar mass content of haloes in the no supernova feedback model. Here, reionization feedback is still included using the fiducial escape fraction of ; however, it is unable to counter the runaway star formation which occurs in the absence of supernova feedback. The net result is a large boost in the number densities of galaxies with respect to the fiducial model () for . At higher masses, supernova feedback becomes inefficient and the mass function converges to the fiducial result. This is because these large galaxies preferentially reside in the most massive haloes where supernovæ are unable to provide the required energy to heat gas to/beyond the virial temperature, thus preventing it from being used for further star formation (cf. Section 2.5 and equations therein).

The red dash–dotted line () in Fig. 6 shows the results of our no reionization feedback model. Here we use the fiducial escape fraction of in order to calculate the ionization state of the IGM; however, we decouple the effect of photoionization suppression from the infall of baryons into each FoF group by setting (cf. Equation 1) for all groups at all times. Hence, in this model variation, galaxy evolution proceeds independently of the ionization state of the IGM. By turning off reionization feedback in this manner whilst still leaving the strong supernova feedback required to reproduce the observed high- stellar mass functions in place, we find little change to the stellar mass function during the EoR (right panel) with respect to the fiducial model (). However, the effects of reionization are cumulative over time (see Equation 32) and by (left panel) we find the space density of galaxies with stellar masses less than is increased by up to relative to the fiducial result. However, this effect is small compared to that found in the absence of supernova feedback. Furthermore, as we move to higher masses, the relative differences rapidly decrease. This is a reflection of the fact that reionization feedback is only effective in low-mass haloes (hosting typically low-mass galaxies) with shallow potential wells susceptible to accretion suppression from the UVB (cf. Section 4.1).

By comparing the no feedback () and no supernova feedback () model lines, we can investigate the isolated impact of reionization on the growth of galaxy stellar mass. In the former model, there is no reionization or supernova feedback included; however, in the latter we turn on the photoionization suppression of baryon accretion. This results in a clear decrement in the number density of low mass galaxies with at . Higher mass galaxies are largely unaffected since, by the time their massive host haloes were exposed to the UVB, they already provided a potential well deep enough to accrete gas despite the presence of the photoionizing UVB. During reionization (right-hand panel), we see that the effects of reionization feedback are more modest due to a smaller fraction of galaxies being exposed to the UVB as well as the typically shorter exposure times for those which have.

Simply turning off supernova feedback whilst leaving the remaining model parameters constant (no supernova feedback; ) leads to an over-prediction of the total =5 stellar mass density by a factor of 5 relative to the fiducial model as well as a mass function slope which is too steep to be consistent with observations. This enhanced star formation may lead to an under-estimate in the importance of reionization feedback. In order to test this hypothesis we have run a recalibrated no supernova feedback model () with a reduced star formation efficiency parameter, , chosen to provide the same total =5 stellar mass density as the fiducial model. The recalibrated no feedback model () additionally shows the result of this altered with reionization feedback also omitted. Again, these models predict a stellar mass function which is too steep to be consistent with observations. However, we find that there is no combination of remaining parameters in our model which can reproduce the slope of the observed stellar mass function in the absence of supernova feedback.

The small relative difference between the recalibrated no supernova feedback and recalibrated no feedback models in Fig. 6 indicates that reionization feedback is even less effective at regulating galaxy growth than was the case with the original no supernova feedback and no feedback variations. This is because our fiducial star formation efficiency parameter results in stellar mass growth in the majority of small haloes being gas supply limited. By reducing the star formation efficiency parameter in the recalibrated models, this is no longer the case and there is more gas available in these systems than can be converted into stars in a single time step. Reionization-driven photosuppression of accretion into these gas-rich systems therefore has little impact on the growth of stellar mass, resulting in a quantitatively similar change to the stellar mass function as is found in the presence of supernova feedback (i.e. comparing the fiducial and no reionization feedback model results).

Together, these results highlight an over-estimate of the importance of reionization feedback for regulating star formation in models which do not employ the galactic feedback processes necessary to reproduce the observed stellar mass functions at high-. Alternatively, at the very least, a star formation efficiency low enough to result in a gas-rich star formation scenario is needed.

The stellar mass content of haloes

Figure 7: The median fraction of baryons in the form of stars, , as a function of FoF group mass, at (left) and (right). Thick lines show the results of Meraxes when run on the full Tiamat simulation merger trees, whilst thin lines indicate the results of running on the higher resolution Medi Tiamat trees. The blue solid lines and surrounding shaded regions show the median result of the fiducial (Tiamat) model and associated 68 and 95% confidence intervals. The magnitudes of these statistical uncertainties are representative of those of all of the models shown in each panel. The grey shaded region at the left of each panel denotes halo masses below the atomic cooling mass threshold, , corresponding to a virial temperature of . The green data points in the panel display the SHAM results of Behroozi et al. (2013b) which are in excellent agreement with our fiducial model. Comparison between the no feedback (gold dashed) and no supernova feedback (grey solid) lines demonstrates that reionization feedback is most effective in low mass haloes. However, supernova feedback of the level required to reproduce the observed high- stellar mass functions dominates the suppression of star formation across all halo masses. The strong halo mass dependence of contrasts the constant value assumed by the majority of reionization structure studies.

In Section 4.1.1, we demonstrated that the impact of reionization on regulating the growth of galaxies and the production of ionizing photons in Meraxes is minimal owing to the importance of supernova feedback. In this section, we explore this topic further by investigating the stellar mass content of haloes as predicted by a subset of our model variations, both subsequent to and during reionization. In Fig. a, we present the fraction of baryons in the form of stars, , as a function of FoF group virial mass, at and 8 (the latter redshift corresponding to in the fiducial model). Thick lines indicate the results from running Meraxes on the full Tiamat simulation. Thin lines show a subset of the same models run on the higher resolution Medi Tiamat trees which are complete down to the atomic cooling mass threshold at both redshifts shown (cf. Appendix A).

At low halo masses, there are minor differences between the fiducial model results of each simulation. These are predominantly driven by an increased prevalence of merger-driven halo mass growth in Medi Tiamat which is unresolved in the lower resolution Tiamat trees. These merger events occur less frequently than in situ star formation episodes, but they are more efficient, giving rise to more energetic supernova feedback episodes capable of ejecting significant amounts of material from haloes and temporarily halting star formation. The net result is a reduction in the stellar mass growth of low mass haloes in the higher resolution trees. During reionization (right-hand panel), the lower atomic cooling mass threshold results in Tiamat missing a larger fraction of the lowest mass haloes than is the case at . This can be seen by comparing the no feedback () model lines from each simulation. However, when feedback is included, the stellar mass content of these haloes is greatly reduced and the results of the different simulations again come to a good agreement. We also note that despite the minor discrepancies at low masses, there is excellent overall agreement between the Tiamat and Medi Tiamat results across all other masses in both panels. Hence, Fig. a demonstrates that our full Tiamat volume has sufficient mass resolution to correctly model the growth of galaxies across the range of masses and redshifts relevant for this work.

Our fiducial model () predicts a strong decline in the stellar mass content of haloes with . For comparison, a constant stellar baryon fraction of is commonly employed by previous studies utilizing 21cmFAST (e.g. Mesinger et al., 2011; Sobacchi & Mesinger, 2013), as well as previous N-body-based radiative transfer calculations (e.g. Iliev et al., 2007). Whilst it is important to remember that, for the purposes of reionization, the precise value of is degenerate with other quantities such as the escape fraction of ionizing photons, it is clear that the approximation of a constant value with halo mass is poor. Correctly predicting and self-consistently utilizing this relation is a key feature of Meraxes, as is investigating and predicting further potential contributing variables such as environmental density (see Section 4.5 below) and ionizing escape fraction.

For comparison, we have also plotted the relation and 1- scatter found by the SHAM study of Behroozi et al. (2013b, green error bars). Our fiducial model () shows excellent agreement with these results.7 This is perhaps not unexpected given that both studies have utilized N-body simulations with free model parameters constrained to match observed high- stellar mass functions. However, the agreement is still noteworthy given that we have utilized different N-body simulations, halo finders, observational data sets and methodologies. Less obvious is the high level of agreement between the scatter in values at a fixed which is an unconstrained prediction of our Meraxes results. We also highlight that our model provides predictions down to masses at least two orders of magnitude lower than can be directly probed by current observations and SHAM studies.

A notable feature of the no feedback model () is the sharp turn-over in the