Simulating the effect of photoheating feedback during reionization

TL;DR

This study uses advanced radiation hydrodynamic simulations to assess how photoheating during reionization influences galaxy formation, finding it less impactful than previously thought and not leaving clear observable signatures.

Contribution

It provides a self-consistent simulation framework showing that photoheating feedback is weaker and less observable than earlier models suggested.

Findings

Photoheating suppresses star formation mainly in halos below ~10^{8.4} M_sun.

Uniform UV background causes earlier star formation suppression than the fiducial model.

Photoheating alone cannot significantly quench star formation in halos above ~10^9 M_sun.

Abstract

We present self-consistent radiation hydrodynamic simulations of hydrogen reionization performed with Arepo-RT complemented by a state-of-the-art galaxy formation model. We examine how photoheating feedback, due to reionization, shapes the galaxies properties. Our fiducial model completes reionization by $z\approx6$ and matches observations of the Ly$\alpha$ forest, the CMB electron scattering optical depth, the high-redshift UV luminosity function, and stellar mass function. Contrary to previous works, photoheating suppresses star formation rates by more than $50\%$ only in halos less massive than $\sim10^{8.4}\ M_\odot$ ($\sim10^{8.8}\ M_\odot$) at $z=6$ $(z=5)$, suggesting inefficient photoheating feedback from photons within galaxies. The use of a uniform UV background that heats up the gas at $z\approx10.7$ generates an earlier onset of suppression of star formation compared to our fiducial model. This discrepancy can be mitigated by adopting a UV background model with a more realistic reionization history. In the absence of stellar feedback, photoheating alone is only able to quench halos less massive than $\sim10^9\ M_\odot$ at $z\gtrsim5$, implying that photoheating feedback is sub-dominant in regulating star formation. In addition, stellar feedback, implemented as a non-local galactic wind scheme in the simulations, weakens the strength of photoheating feedback by reducing the amount of stellar sources. Most importantly, photoheating does not leave observable imprints in the UV luminosity function, stellar mass function, or the cosmic star formation rate density. The feasibility of using these observables to detect imprints of reionization therefore requires further investigation.

Tables

Table 1. Table 1: Mean cross sections and photon energies above the ionization thresholds of each species used in the simulations. The energy intervals of the three frequency bins traced in the simulations are indicated in eV in the first column. The other six columns show σ ¯ X subscript ¯ 𝜎 𝑋 \bar{\sigma}_{X} and ϵ ¯ X subscript ¯ italic-ϵ 𝑋 \bar{\epsilon}_{X} for each species X ∈ [ H  i , He  i , He  ii ] 𝑋 H  i He  i He  ii X\in[{\rm\text{H\,{i}},\text{He\,{i}},\text{He\,{ii}}}] , which are the mean cross section and latent heat per ionization calculated using the zero-age zero-metallicity spectrum of the Bruzual & Charlot ( 2003 ) model.

Energy Range	${\overline{\sigma }}_{\text{H }\text{i}}$	${\overline{ϵ}}_{\text{H }\text{i}}$	${\overline{\sigma }}_{\text{He }\text{i}}$	${\overline{ϵ}}_{\text{He }\text{i}}$	${\overline{\sigma }}_{\text{He }\text{ii}}$	${\overline{ϵ}}_{\text{He }\text{ii}}$
$\left[\mathrm{eV}\right]$	$\left[{\mathrm{cm}}^2\right]$	$\left[\mathrm{eV}\right]$	$\left[{\mathrm{cm}}^2\right]$	$\left[\mathrm{eV}\right]$	$\left[{\mathrm{cm}}^2\right]$	$\left[\mathrm{eV}\right]$
$13.6-24.6$	$3.2\times {10}^{-18}$	3.4	0	0	0	0
$23.6-54.4$	$6.0\times {10}^{-19}$	17.0	$3.9\times {10}^{-18}$	6.1	0	0
$54.4-100.0$	$1.0\times {10}^{-19}$	43.9	$7.2\times {10}^{-19}$	33.0	$1.3\times {10}^{-18}$	3.2

Table 2. Table 2: The simulations used in this work. The table lists the volume side length, choice of stellar feedback, dark matter particle mass, mean gas cell mass in the initial conditions, gravitational softening length, and source of photoheating of each simulation.

Name	$L_{\mathrm{box}}$	Winds	$m_{\mathrm{DM}}$	$m_{\mathrm{target}}$	Softening Length	Source of
	$\left[\mathrm{cMpc}/h\right]$	(Stellar Feedback)	$\left[M_{\odot }\right]$	$\left[M_{\odot }\right]$	$\left[\mathrm{ckpc}/h\right]$	Photoheating
L6n256 fiducial-noRT	6	Yes	$1.4\times {10}^6$	$2.2\times {10}^5$	0.24	None
L6n256 fiducial-UVB	6	Yes	$1.4\times {10}^6$	$2.2\times {10}^5$	0.24	FG09 UVB
L6n256 fiducial-RT	6	Yes	$1.4\times {10}^6$	$2.2\times {10}^5$	0.24	RT
L6n256 NW-noRT	6	No	$1.4\times {10}^6$	$2.2\times {10}^5$	0.24	None
L6n256 NW-UVB	6	No	$1.4\times {10}^6$	$2.2\times {10}^5$	0.24	FG09 UVB
L6n256 NW-RT	6	No	$1.4\times {10}^6$	$2.2\times {10}^5$	0.24	RT
L3n256 fiducial-noRT	3	Yes	$1.8\times {10}^5$	$2.8\times {10}^4$	0.12	None
L3n256 fiducial-UVB	3	Yes	$1.8\times {10}^5$	$2.8\times {10}^4$	0.12	FG09 UVB
L3n256 fiducial-RT	3	Yes	$1.8\times {10}^5$	$2.8\times {10}^4$	0.12	RT

Equations

$$(1-f){\left[1+{\left(\frac{n_{\mathrm{H}}}{n_0}\right)}^{\beta }\right]}^{{\alpha }_1}+f{\left[1+\frac{n_{\mathrm{H}}}{n_0}\right]}^{{\alpha }_2}.$$

$$u_{\mathrm{eff}}=(1-x)u_{\mathrm{h}}+x{u}_{\mathrm{c}},$$

$$t_{\ast }=2.2\sqrt{\frac{n_{\mathrm{th}}}{n}}\mathrm{Gyr},$$

$$v_{\mathrm{w}}={\kappa }_{\mathrm{w}}{\sigma }_{\mathrm{DM}}^{1D},$$

$${\eta }_{\mathrm{w}}=\frac{\mathrm{eg}{\mathrm{y}}_{\mathrm{w}}}{v_{\mathrm{w}}^2},$$

$$p=\frac{M}{M_{\ast }}\left[1-\exp \left(-\frac{(1+{\eta }_{\mathrm{w}})\Delta t}{t_{\mathrm{SF}}}\right)\right]$$

$$\begin{gathered} {\tilde{N}}_i=a^3{N}_i, \\ {\tilde{\mathbf{F}}}_i=a^3{\mathbf{F}}_i, \end{gathered}$$

$$\frac{\partial {\tilde{E}}_i}{\partial t}+\frac{1}{a}\nabla \cdot {\tilde{\mathbf{F}}}_i=0,$$

$$\frac{1}{c}\frac{\partial {\tilde{\mathbf{F}}}_i}{\partial t}+\frac{c}{a}\nabla \cdot {\tilde{\mathbb{P}}}_i=0,$$

$${\sigma }_{ij}=\frac{{\int }_{{\nu }_{i,1}}^{{\nu }_{i,2}}\frac{4\pi J_{\nu }}{h\nu }{\sigma }_{\nu j}\mathrm{d}\nu }{{\int }_{{\nu }_{i,1}}^{{\nu }_{i,2}}\frac{4\pi J_{\nu }}{h\nu }\mathrm{d}\nu },$$

$${ϵ}_{ij}=\frac{{\int }_{{\nu }_{i,1}}^{{\nu }_{i,2}}\frac{4\pi J_{\nu }}{h\nu }{\sigma }_{\nu j}(h\nu -h{\nu }_{tj})\mathrm{d}\nu }{{\int }_{{\nu }_{i,1}}^{{\nu }_{i,2}}\frac{4\pi J_{\nu }}{h\nu }{\sigma }_{\nu j}\mathrm{d}\nu },$$

$${\Gamma }_{\text{Hi}}=\sum _{i=0}^3\tilde{c}{N}_i{\sigma }_{i,\text{Hi}}$$

$$\tau =c{\sigma }_{\mathrm{Th}}{\int }_{z_0}^0{n}_{\mathrm{e}}(z)\frac{\mathrm{d}t}{\mathrm{d}z}\mathrm{d}z,$$

$$s_{\mathrm{w}}(M_{\mathrm{vir}})=\frac{\mathrm{SFR}(\mathrm{fiducial}\text{-}\mathrm{noRT})}{\mathrm{SFR}(\mathrm{NW}\text{-}\mathrm{noRT})},$$

$$s_{\mathrm{h}}(M_{\mathrm{vir}})=\frac{\mathrm{SFR}(\mathrm{NW}\text{-}\mathrm{RT})}{\mathrm{SFR}(\mathrm{NW}\text{-}\mathrm{noRT})}.$$

$$s_{\mathrm{wh}}(M_{\mathrm{vir}})=\frac{\mathrm{SFR}(\mathrm{fiducial}\text{-}\mathrm{RT})}{\mathrm{SFR}(\mathrm{NW}\text{-}\mathrm{noRT})}.$$

Full text

††pubyear: 2019††pagerange: Simulating the effect of photoheating feedback during reionization–C

1 Introduction

Cosmological hydrodynamical simulations are among the most powerful tools to study the process of galaxy formation. One of the major challenges for galaxy formation models is to include realistic feedback mechanisms that can regulate gas cooling and star formation. These feedback processes are crucial for simulating realistic galaxy populations (e.g. Vogelsberger et al., 2014a; Schaye et al., 2015). Among them, stellar feedback in the form of galactic winds driven by supernovae (SNe) is a key ingredient in reducing star formation across a large range of halo masses and is particularly effective in suppressing the faint end slope of the galaxy luminosity function (e.g. Benson et al., 2003). Stellar feedback has been invoked in simulations to match the observed galaxy abundances and scaling relations between stars, gas and metals (e.g. Davé et al., 2006, 2011a, 2011b; Vogelsberger et al., 2013).

Stellar feedback has been extensively studied, producing a comprehensive understanding of its role in galaxy formation. However, at high redshifts ($z\gtrsim 5$), photoionization heating due to the reionization process provides another form of feedback, which particularly affects low-mass halos. The epoch of reionization is the era when radiation from the first stars and galaxies reionized the intergalactic medium (IGM), turning it from a cold and neutral medium to a hot and highly ionized one with temperatures of $\sim 20,000-30,000$ K (e.g. Miralda-Escudé & Rees, 1994; McQuinn, 2012). During this period, the virial temperatures of halos less massive than $\sim 10^{9}\ M_{\odot}$ became lower than the mean IGM temperature, leading to suppression of gas accretion onto these objects (Thoul & Weinberg, 1996; Gnedin, 2000; Hoeft et al., 2006; Okamoto et al., 2008; Noh & McQuinn, 2014; Katz et al., 2019). Moreover, such halos gradually lose their baryon content because their shallow potentials can no longer hold the photoheated gas, leading to an overall reduction of the star formation rate (SFR) of these halos (Petkova & Springel, 2011; Finlator et al., 2011; Hasegawa & Semelin, 2013; Ocvirk et al., 2016; Finkelstein et al., 2019). Since halos less massive than $\sim 10^{8}\ M_{\odot}$ can be easily disrupted by a single SN explosion (Finlator et al., 2011), halos in the mass range $\sim 10^{8}-10^{9}\ M_{\odot}$ are the most sensitive to the effects of photoheating feedback (for a recent review on the back-reaction of reionization on galaxy formation, see Dayal & Ferrara, 2018).

Some earlier works studied the mass loss and the suppression of star formation of low-mass halos due to photoheating using a spatially uniform UV background (UVB; e.g. Hoeft et al., 2006; Okamoto et al., 2008; Pawlik & Schaye, 2009). Large-scale cosmological simulations of galaxy formation (e.g. Vogelsberger et al., 2014a; Dubois et al., 2014; Okamoto et al., 2014; Schaye et al., 2015; Davé et al., 2016; Weinberger et al., 2017; Pillepich et al., 2018a, b; Nelson et al., 2018; Springel et al., 2018; Naiman et al., 2018; Marinacci et al., 2018) also often adopt a homogeneous time-varying UVB as an approximation of reionization (e.g. Haardt & Madau, 2012; Faucher-Giguère et al., 2009). However, reionization is usually believed to be a spatially inhomogeneous and temporally extended process. It is unclear whether a uniform UVB and a patchy reionization produce the same amount of suppression of star formation in low-mass halos. To better address this issue, radiative transfer (RT), or more precisely, radiation hydrodynamics (RHD) simulations, are needed to model the growth of ionized bubbles in a self-consistent manner.

The suppression of star formation in low-mass halos by photoheating feedback has been argued to be observable. For instance, the cosmic star formation rate density (SFRD) may experience a drop during the epoch of reionization (Barkana & Loeb, 2000). The faint-end slope of the galaxy UV luminosity function (UVLF) may also be sensitive to the reionization history (Gardner et al., 2006). However, if stellar feedback dominates the regulation of star formation, imprints of photoheating feedback on these observables can become less evident (e.g. Mutch et al., 2016). The contribution of low-mass galaxies to reionization may also be reduced (e.g. Wyithe & Loeb, 2013). A thorough understanding of the interplay between stellar and photoheating feedback requires self-consistent RHD simulations.

A number of previous works have explored the above-mentioned problems using RHD simulations. For instance, Finlator et al. (2011) and Ocvirk et al. (2016) both found that at the end of reionization, there is a sharp decrease in SFR of halos less massive than $\sim 10^{9}\ M_{\odot}$ due to photoheating feedback. In Ocvirk et al. (2016), even halos of $10^{10}-10^{11}\ M_{\odot}$ show a factor of $\sim 2$ difference in their SFR when compared to simulations without RT. Finlator et al. (2011) also illustrated that for low-mass halos ($\lesssim 10^{9}\ M_{\odot}$), stellar feedback weakens the strength of photoheating feedback on suppressing star formation by reducing the radiation field produced by stellar sources. Contrarily, Pawlik et al. (2015) found an amplification of the effect of stellar feedback by photoheating feedback (for a larger suite of simulations, see Pawlik et al., 2017). However, different simulations seem to agree that photoheating is sub-dominant in regulating star formation compared to stellar feedback, and that the latter plays the major role in shaping the galaxy properties (e.g. Rosdahl et al., 2018). The effects of photoheating feedback on observables are also unclear. While Gnedin & Kaurov (2014) did not see a drop in the cosmic SFRD or a significant change in the faint end slope of the UVLF in the reionization simulations of Gnedin (2014), Finlator et al. (2018) found a small change in the UVLF at ${\rm M_{1500}}>-14$ mag.

In this work we present RHD simulations run with arepo-rt coupled to the Illustris galaxy formation model (Vogelsberger et al., 2013, 2014a, 2014b; Genel et al., 2014; Sijacki et al., 2015; Nelson et al., 2015) in order to explore the aforementioned open questions. The Illustris model has been shown to be able to reproduce a number of observed properties of galaxies (e.g. Vogelsberger et al., 2014a; Genel et al., 2014) and the IGM (e.g. Vogelsberger et al., 2014a; Bird et al., 2014) at various redshifts, making it a state-of-the-art model for galaxy formation studies. By post-processing the Illustris simulation with RT, Bauer et al. (2015) showed that the Illustris star formation history is able to generate a realistic hydrogen reionization history assuming rather low escape fractions ($\lesssim 20\%$). In this paper we show how reionization proceeds in the Illustris model when the radiation field is evolved self-consistently with hydrodynamics. We explore how photoheating feedback due to reionization suppresses star formation in halos of different masses. We also analyze the relative importance of stellar feedback and photoheating feedback by comparing simulations with and without stellar feedback. In addition, we assess the differences in the amount of suppression in star formation by performing RHD versus using the uniform Faucher-Giguère et al. (2009) UVB (hereafter FG09, updated in 2011111https://galaxies.northwestern.edu/uvb/). We briefly examine the feasibility of using the high redshift UVLF and the cosmic SFRD to detect imprints of reionization.

This paper is organized as follows. In Sec. 2 we present the galaxy formation model and the RHD scheme. In Sec. 3 we show the reionization histories in the simulations and discuss the baryon depletion and suppression of star formation in low mass halos due to photoheating feedback. We also explore the implications for observables, including the UVLF, stellar mass function, and cosmic SFRD. We give detailed discussions about the relative importance of photoheating feedback and stellar feedback in Sec. 4, and summarize our work in Sec. 5.

2 Methods

We use the arepo-rt code (Kannan et al., 2019) to solve the coupled equations of gravity, hydrodynamics, and radiative transfer. arepo-rt is a RHD extension of the moving-mesh cosmological hydrodynamic code arepo (Springel, 2010) which uses an unstructured Voronoi tessellation of the computational domain. The mesh-generating points are allowed to move freely, offering significant flexibility for representing the geometry of the flow. The mesh is then used to solve the equations of ideal hydrodynamics using a second-order unsplit Godunov scheme with an exact Riemann solver. arepo has been shown to surpass traditional smoothed particle hydrodynamics (SPH) and adaptive mesh refinement (AMR) codes in terms of its accuracy (Kereš et al., 2012; Sijacki et al., 2012; Torrey et al., 2012; Vogelsberger et al., 2012; Genel et al., 2013; Nelson et al., 2013). Gravitational forces are computed using a Tree-PM scheme (Xu, 1995), where short-range and long-range forces are calculated using a hierarchical octree algorithm (Barnes & Hut, 1986) and a Fourier particle-mesh method respectively. We briefly describe the galaxy formation model and the RHD scheme in detail below, which are the key modules of the code needed to perform the simulations presented in this work.

2.1 Galaxy Formation Model

We adopt the galaxy formation model outlined in Vogelsberger et al. (2013). Briefly, gas cells change their internal energy via radiative cooling and heating processes including collisional excitation, collisional ionization, recombination, dielectric recombination, bremsstrahlung, Compton cooling off the cosmic microwave background (CMB), photoionization and photoheating (Katz et al., 1996). In the original Illustris implementation, gas is assumed to be in ionization equilibrium with the spatially uniform and time-dependent FG09 UVB. The FG09 UVB includes contributions from quasars and star-forming galaxies, the latter dominating at $z\gtrsim 3$. It was calibrated to satisfy the observed mean transmission of the Ly$\alpha$ forest at $z=2-4.2$ (Faucher-Giguère et al., 2008a, b), have He ii reionization by $z\sim 3$ (McQuinn et al., 2009), and complete hydrogen reionization by $z=6$. Gas self-shielding is taken into account at $z<6$ by suppressing the photoionization and photoheating rates by a factor of (see Rahmati et al., 2013, equation A1):

[TABLE]

where $n_{\rm H}$ is the physical hydrogen number density of the cell. The parameters $(\alpha_{1},\alpha_{2},\beta,f,n_{0})$ are linearly interpolated in redshift with the values given in table A1 of Rahmati et al. (2013). In the simulations that adopt the FG09 UVB, we continue to use this setup for the treatment of radiative cooling. The RHD implementation of gas cooling adopts a non-equilibrium hydrogen and helium thermochemistry network, which will be presented in Sec. 2.2.

Gas cooling triggers star formation. We follow the scheme of Springel & Hernquist (2003, to which we refer the reader for more details) and model the star-forming interstellar medium (ISM) gas using an effective equation of state (eEOS). Specifically, we describe the star-forming ISM as a fluid composed of dense cold clouds in pressure equilibrium with an ambient hot gas. Assuming equilibrium, it can be shown that the effective internal energy per unit mass of the two-phase gas is given by

[TABLE]

where $u_{\rm h}$ and $u_{\rm c}$ are the internal energy per unit mass of the hot and cold phases respectively, and $x$ is the mass fraction of the cold gas (computed by the model as a function of gas density). Equation 2 defines the effective equation of state for the star-forming gas.

We consider a gas cell to be star-forming when it exceeds the physical number density threshold222This is also the density above which the eEOS is imposed. of $n_{\rm th}\simeq 0.13$ cm*-3*. Following Springel et al. (2005), we determine the temperature of a star-forming gas cell via a weighted mean between the full Springel & Hernquist (2003) eEOS value and an isothermal EOS at $10^{4}$ K. In computing the mean, we assign a weight $q=0.3$ to the eEOS value and, correspondingly, a weight $1-q$ to the isothermal EOS. For a star-forming gas cell, its star formation timescale is given by

[TABLE]

where $n$ indicates the physical gas number density. We calculate the SFR as the ratio of the cold gas mass and $t_{*}$.

We assume that each stellar particle represents a co-eval, single metallicity stellar population which follows a Chabrier (2003) initial mass function (IMF). We calculate mass and metal return due to stellar evolution by integrating over this IMF the time evolution of stellar particles and using information from stellar evolution calculations on the expected main-sequence lifetime, mass return fraction, and heavy element production for a wide range of initial stellar masses and metallicities. We track nine chemical elements – H, He, C, N, O, Ne, Mg, Si, Fe, and the total gas phase metallicity.

The stellar feedback implementation adopts a non-local energy-driven wind model. In this model winds are directly launched from the star-forming ISM gas in the form of wind particles. After being created, these particles are decoupled from hydrodynamic forces, but not the gravitational forces, until they travel to a region with density below $0.1$ times the density threshold for star formation, or a maximum travel time of $50$ Myr has elapsed. When either of these two criteria is satisfied, we recouple the wind particle and deposit its mass, momentum, thermal energy, and metals into the gas cell where it is currently located. The initial wind velocity $v_{\rm w}$ is

[TABLE]

where $\kappa_{\rm w}=3.7$ is a dimensionless model parameter, and $\sigma^{\rm 1D}_{\rm Dm}$ is the local one-dimensional dark matter velocity dispersion at the current position of the gas cell. We determine the mass carried by galactic winds by computing the mass loading factor $\eta_{\rm w}$, which specifies the ratio of the wind mass flux to the star formation rate:

[TABLE]

Here ${\rm egy_{w}}=1.89\times 10^{49}$ erg ${\rm M}^{-1}_{\odot}$ is the specific energy available for wind generation, i.e. the available Type II SN energy per formed stellar mass. We assume the newly created wind particle has a metallicity that is 0.4 times that of the ambient ISM. The direction of the ejection velocity of wind particles is randomly drawn.

We probabilistically select star-forming gas cells to be converted either into stellar or wind particles, according to the values of star formation rate and wind mass loading factor computed by the model. At each time step $\Delta t$, we draw a random number $x$ from a uniform distribution $U(0,1)$. If $x<1/(1+\eta_{\rm w})$, we treat the spawning of star particles. Otherwise we consider launching winds. The probability of spawning a star or wind particle of mass $M_{*}$ from a gas cell of mass $M$ is given by

[TABLE]

where $t_{\rm SF}$ is the ratio of the cell mass and the cell SFR. The star/wind particle mass $M_{*}$ is set as follows: if $M<2m_{\rm target}$, then $M_{*}=M$ and the full gas cell is converted into a star or wind particle. Otherwise the cell only spawns a star (wind) particle of mass $m_{\rm target}$. $m_{\rm target}$ is the mean gas cell mass in the initial conditions (see Table 2 for values). We employ a (de-)refinement scheme that keeps the cell masses close (within a factor of 2) to $m_{\rm target}$.

Our simulations are only run until $z=5$, by which redshift hydrogen reionization has completed. Therefore, we do not include metal-line cooling or black hole formation and feedback in our galaxy formation model. This is done for simplicity, but these processes are not expected to have a significant impact at $z\gtrsim 5$ (e.g. Ocvirk et al., 2016). Nevertheless, we use stellar metallicities to calculate the luminosity of star particles, which determines how many photons a star particle should emit and deposit into its surrounding gas cells per unit time (see Sec. 2.2).

2.2 Radiative Transfer

The RT implementation solves the moment-based radiative transfer equations using the M1 closure relation (Levermore, 1984) on a moving mesh (see Kannan et al., 2019, for a detailed description of this scheme). We divide the UV continuum into three frequency bins relevant for hydrogen and helium photoionization: $[13.6,24.6]$ eV, $[24.6,54.4]$ eV, $[54.4,100]$ eV. For each frequency bin $i$, we evolve the comoving photon number density $\tilde{N}_{i}$ and photon flux $\tilde{\mathbf{F}}_{i}$, which are related to the physical photon density $N_{i}$ and photon flux $F_{i}$ via

[TABLE]

where the scale factor $a$ is adopted to account for the loss of photon energy due to cosmological expansion (e.g. Rosdahl et al., 2013). Assuming that the Universe does not expand significantly before a photon is absorbed, the transport equations take the form

[TABLE]

where $\tilde{\mathbb{P}}_{i}$ is the radiation pressure tensor and is related to $\tilde{E}_{i}$ via the Eddington tensor.

We solve photon transport with an explicit scheme, which constrains the simulation time step by the Courant condition. To lower the computational cost we use the reduced speed of light approximation (Gnedin & Abel, 2001) with $\tilde{c}=0.1c$, where $\tilde{c}$ and $c$ are the reduced and actual speed of light respectively. To further reduce the computing time we perform 32 RT sub-cycles for each hydro time step. In each RT step, we advect radiation by solving the Riemann problem at each cell interface and computing the flux using Godunov’s approach (Godunov, 1959). We adopt a Global-Lax-Friedrich flux function (Rusanov, 1961), and achieve second order accuracy by replacing the piecewise constant approximation of Godunov’s scheme with a slope limited linear spatial extrapolation and a first order time prediction step to obtain the values of the primitive variables on both sides of the cell interface. We perform the spatial extrapolations using a local least square fit gradient estimate (Pakmor et al., 2016).

During an RT time step, besides advecting photons, we track the non-equilibrium hydrogen and helium thermochemistry. To do so we use an implicit scheme that takes into account the same radiative processes described in Sec. 2.1. We adopt the on-the-spot approximation, assuming that recombination emission is absorbed within the same cell. We trace the ionization fractions of hydrogen and singly and doubly ionized helium for each gas cell. In each hydro time step, we advect these ionization fractions as passive scalars along with the gas.

In our simulations, star particles are the only source of the radiation. We compute the number of photons a star particle emits based on its spectral energy distribution (SED), which is a function of both its age and metallicity, as given by Bruzual & Charlot (2003). We integrate the SED in each frequency bin to calculate the number of photons to deposit into the surrounding gas cells. Each neighboring gas cell receives a fraction of the total emitted photons of the star particle, proportional to their volume weighted by the evaluation of an SPH cubic spline kernel. The kernel smoothing length is defined as the "standard" SPH smoothing length, which is the length enclosing a predefined number of effective neighbours (in our case $64$). To take into account absorption of photons on unresolved scales, we assume an escape fraction $f_{\rm esc}$ for each star particle, which represents the escape fraction from the birth cloud. We adopt $f_{\rm esc}=0.7$ for all star particles in all the RT simulations. This choice ensures that reionization completes at $z\approx 6$ in the simulation with the fiducial stellar feedback model, i.e. the volume-average H i fraction drops to $\sim 10^{-4}$ at $z=6$.

In order to calculate the photoionization and photoheating rates for a gas cell, for each species $j\in[{\rm\text{H\,{i}},\text{He\,{i}},\text{He\,{ii}}}]$ we compute the mean ionization cross section in each frequency bin $i$ that runs from frequency $\nu_{i,1}$ to $\nu_{i,2}$:

[TABLE]

where $J_{\nu}$ is the mean specific intensity. We also calculate the latent heat per photoionization event of species $j$ as:

[TABLE]

where $h\nu_{tj}$ is the ionization potential of the ionic species $j$. In principle $\sigma_{ij}$ and $\epsilon_{ij}$ vary among gas cells due to the different shapes of the spectrum received by the gas cells, which the current code is unable to track. Rosdahl et al. (2013) circumvents this by assuming the same $\sigma_{ij}$ and $\epsilon_{ij}$ for all gas cells and updating them every 10 coarse time steps from the luminosity-weighted averages of the spectra of all star particles in the simulation volume, making $\sigma_{ij}$ and $\epsilon_{ij}$ representative of the average photon population. We note that for the Bruzual & Charlot (2003) spectra, star particles emit most of their photons during the first $\sim 5$ Myr of their lifetime, when the calculated $\sigma_{ij}$ and $\epsilon_{ij}$ stay roughly constant and do not change much with metallicity (Fig. B2 of Rosdahl et al., 2013). We therefore calculate $\sigma_{ij}$ and $\epsilon_{ij}$ using the zero-age zero-metallicity spectrum of the Bruzual & Charlot (2003) model and adopt the same values of $\sigma_{ij}$ and $\epsilon_{ij}$ for all gas cells333We note that spectral hardening during photon propagation or using a stellar SED model harder than Bruzual & Charlot (2003) can heat the gas to higher temperatures. Appendix B will show that a $\sim 10,000$ K difference in the halo gas temperature does not have a strong impact on the suppression of halo SFR. Moreover, the IGM temperature is only weakly dependent on the spectral slope (D’Aloisio et al., 2018). We therefore expect harder spectra to have minor effects on our results in Sec. 3.. The resulting $\sigma_{ij}$ and $\epsilon_{ij}$ are tabulated in Table 1.

2.3 Simulations

Table 2 summarizes the key features of our simulations. The fiducial simulations have a volume of $(6\ {\rm cMpc}/h)^{3}$ with $256^{3}$ dark matter particles and an initial number of $256^{3}$ gas cells (denoted as L6n256). We run two sets of L6n256 simulations, one using the fiducial stellar feedback model outlined in Sec. 2.1, and the other without stellar feedback (tagged as "NW"). For each set of these simulations we run three variations, one without RT or UVB, one with the FG09 UVB, and another with RT, named with postfixes "-noRT", "-UVB", "-RT" respectively. In order to check the convergence of our results, we run three additional simulations with a $3\ {\rm cMpc}/h$ side length and $2\times 256^{3}$ resolution elements (L3n256) using the fiducial stellar feedback model, adopting no RT, FG09 UVB, and RT respectively. Appendix A discusses results of the convergence tests.

We adopt a Planck 2016 cosmology with $\Omega_{m}=0.3089$, $\Omega_{\Lambda}=0.6911$, $\Omega_{b}=0.0486$, $h=0.6774$, and $\sigma_{8}=0.8159$ (Planck Collaboration et al., 2016). Hence the L6n256 and L3n256 simulations have dark matter particle masses of $1.4\times 10^{6}\ M_{\odot}$ and $1.8\times 10^{5}\ M_{\odot}$ respectively. The minimum gravitational softening lengths are $0.24$ ckpc$/h$ and $0.12$ ckpc$/h$ in L6n256 and L3n256 respectively. Gas cells use an adaptive softening length tied to the cell radius, limited by a minimum value of $0.03$ ckpc$/h$ in L6n256 and $0.015$ ckpc$/h$ in L3n256 respectively. Our fiducial L6n256 runs are thus able to resolve halos of $10^{8}\ M_{\odot}$ with $\sim 100$ dark matter particles.

We identify dark matter halos by a friends-of-friends (FOF) algorithm with a minimum particle number of 32 (Davis et al., 1985) and a linking length of 0.2 times the mean particle separation. Stellar particles and gas cells are attached to these FOF primaries in a secondary linking stage (Dolag et al., 2009). We then use the Subfind algorithm to identify gravitationally bound structures (Springel et al., 2001; Dolag et al., 2009). In the following sections, we quote halo mass as the halo virial mass $M_{\rm vir}$, defined as the mass contained in a spherical region with average density that is $200$ times the critical density of the Universe at that time.

3 Results

3.1 Reionization History

We first present how reionization proceeds in our simulations and compare the simulation results with different observational constraints. Fig. 1 shows maps of the H i fraction, gas temperature, ionizing flux density, and gas density in the fiducial-RT run (from the left to right columns), obtained by projecting a slice of the simulation with a dimension of $6\times 6\times 0.5\ ({\rm cMpc}/h)^{3}$. From top to bottom, the maps are taken at $z=10,8.5,7.5,6.8,6.0$, respectively. At $z=10$ the volume-averaged H i fraction is only about $5\%$. Ionized bubbles form around the early galaxies, which lie on the peaks of the cosmological density field. At $z=8.5$ the global ionized fraction reaches $\sim 20\%$, and ionized bubbles are still isolated from each other. By $z=7.5$ the bubbles begin to overlap, and the volume becomes $\sim 40\%$ ionized. When the simulated box is $\sim 70\%$ ionized at $z=6.8$, only two large neutral islands remain in the slice. Finally complete overlap of the ionized bubbles happens at $z\approx 6$. The gas temperature and ionizing flux density evolve in a similar manner, as the ionization fronts sweep through the IGM. Some regions reach temperatures as high as $\sim 25000$ K at $z=6.6$. These maps clearly show how the ionized bubbles grow and overlap with each other, illustrating the patchiness of the reionization process.

Fig. 2 illustrates the volume-averaged H i fraction $\langle x_{\rm\text{H\,{i}}}\rangle^{V}$ as a function of redshift in both the fiducial-RT run (blue line) and the NW-RT run (magenta line), compared with the observations of Fan et al. (2006), Mason et al. (2018), and Hoag et al. (2019) (black triangles, circles, and stars, respectively). In the fiducial-RT run, the reionization process is $50\%$ complete at $z\approx 7$ and finishes at $z\approx 6$, when $\langle x_{\rm\text{H\,{i}}}\rangle^{V}$ drops to $\sim 10^{-4}$. $\langle x_{\rm\text{H\,{i}}}\rangle^{V}$ at $z=7-8$ roughly matches the observations, though slightly lower. The post-reionization $\langle x_{\rm\text{H\,{i}}}\rangle^{V}$ matches the observational data well. We caution, however, that this good match is the consequence of our choices for the values of the escape fraction ($f_{\rm esc}=0.7$) and reduced speed of light ($\tilde{c}=0.1c$). Adopting the actual speed of light with the same $f_{\rm esc}$ would reduce the post-reionization $\langle x_{\rm\text{H\,{i}}}\rangle^{V}$ by a factor of $\sim 10$ (Ocvirk et al., 2018a; Deparis et al., 2019). The NW-RT run has a much earlier-ending reionization and a lower post-reionization $\langle x_{\rm\text{H\,{i}}}\rangle^{V}$. The increased SFR in the NW-RT simulation (about a factor of $10$ higher than the SFR in the fiducial-RT run, see Fig. 6 and Fig. 9) pushes the time of overlap of ionized bubbles to as early as $z\approx 9.5$ and decreases the post-reionization $\langle x_{\rm\text{H\,{i}}}\rangle^{V}$ by about two orders of magnitude.

Fig. 3 presents the simulated hydrogen photoionization rates ($\Gamma_{\rm\text{H\,{i}}}$) compared to the predictions from the FG09 UVB model (yellow line) and the observations of Wyithe & Bolton (2011) and Calverley et al. (2011) (black circles and squares, respectively). $\Gamma_{\rm\text{H\,{i}}}$is the volume-weighted average of the hydrogen photoionization rate in ionized gas, defined as any gas cell having a hydrogen ionized fraction larger than $50\%$. The photoionization rate for each gas cell is given by

[TABLE]

where $\tilde{c}$ is the reduced speed of light, $N_{i}$ is the photon number density of frequency bin $i$, and $\sigma_{i,{\rm\text{H\,{i}}}}$ is the hydrogen photoionization cross section of this frequency bin. Results from the fiducial-RT simulation match the observational data and the FG09 background at $z\sim 5-6$. However, similar to the behavior of $\langle x_{\rm\text{H\,{i}}}\rangle^{V}$, the use of the actual speed of light with the same $f_{\rm esc}$ would raise the post-reionization UVB amplitude by a factor of $\sim 10$ (Ocvirk et al., 2018a). The NW-RT run, due to the enhanced star formation, has a post-reionization UV background that is nearly two orders of magnitude higher than the fiducial-RT run.

The cumulative optical depth to Thomson scattering is another key observable that constrains reionization models. It quantifies the probability of CMB photons scattering off of the free electrons after the epoch of recombination. This optical depth at a redshift $z_{0}$ is calculated as

[TABLE]

where $\sigma_{\rm Th}$ is the Thomson scattering cross section and $n_{\rm e}$ is the number density of free electrons. For our calculations, at $z\geq 5$, $n_{\rm e}$ takes the volume-averaged value obtained from the simulations. From $z=5$ to $3$, $n_{\rm e}/n_{\rm H}=1.08$, since hydrogen reionization is complete and helium is singly ionized. After $z=3$, the time when He ii reionization is usually thought to happen, $n_{\rm e}/n_{\rm H}=1.158$, assuming full ionization of hydrogen and helium. Fig. 4 shows $\tau$ as a function of redshift in the fiducial-RT and NW-RT runs compared to the Planck Collaboration et al. (2018) observations (black lines). A relatively rapid and late-ending reionization in the fiducial-RT run leads to a good match to the observations of Planck Collaboration et al. (2018). A much earlier and/or a much more extended reionization process increases $\tau$, as seen in the NW-RT run. More importantly, $\tau$ is relatively insensitive to the choice of the reduced speed of light value because the evolution of $\langle x_{\rm\text{H\,{i}}}\rangle^{V}$ at $\langle x_{\rm\text{H\,{i}}}\rangle^{V}\gtrsim 0.01$ is independent of the speed of light (see Appendix B, and also Deparis et al., 2019; Ocvirk et al., 2018a). Our fiducial model therefore is able to roughly match various observational constraints on the hydrogen reionization process.

3.2 Halo Properties

We now turn our attention to the effects of photoheating feedback on the properties of low mass ($\lesssim 10^{9}\ M_{\odot}$) halos. Fig. 5 shows the evolution of the median baryon mass fraction as a function of halo mass, at $z=7,6,5$ (in panels from left to right). The baryon fraction is computed as the ratio of the total baryon (gas and stars) mass to the total mass within the halo virial radius, normalized to the global value of $\Omega_{\rm b}/\Omega_{\rm m}\sim 0.157$. Results from the fiducial runs and the NW runs are illustrated in the top and bottom panels, respectively. The black, blue, and red lines show results from the no RT, FG09 UVB, and RT simulations, respectively. Errorbars represent the $1\sigma$ scatter, which are only shown for the fiducial runs to avoid clutter. In both the fiducial-noRT and NW-noRT simulations, the halos retain roughly all their baryons except the $\sim 10^{8}\ M_{\odot}$ ones because the fiducial stellar feedback model has outflow velocities that are lower than the escape velocity of the halo (Genel et al., 2014; Suresh et al., 2015). However, both the UVB and RT simulations show a gradual decrease in the baryon fraction of the $\lesssim 10^{9}\ M_{\odot}$ halos with time due to photoheating. In the fiducial-RT run at $z=7$ when the simulation volume is only $\sim 50\%$ ionized, the baryon fraction–halo mass relation is about the same as that of the fiducial-noRT run. Interestingly, at $z=6$ when ionized bubbles have completely overlapped, there is less than $\lesssim 0.1$ dex reduction in the baryon content of the $10^{8}\ M_{\odot}$ halos in the fiducial-RT run compared to the fiducial-noRT run. At $z=5$, suppression of the baryon fraction can be seen in halos less massive than $\sim 10^{9}\ M_{\odot}$, with the $10^{8}\ M_{\odot}$ halos having a $\sim 80\%$ depletion in their baryons. This delayed response of the halo gas reservoir to the reionization process indicates that the internal photoheating feedback due to photons in the same halo is not efficient at evaporating gas in our simulations, contrary to the findings of Hasegawa & Semelin (2013). It also implies that the low-mass halos likely start to be exposed to the ionized bubbles at late stages of reionization, so external photoheating feedback due to photons from other galaxies takes effect late. This external photoheating feedback can also be delayed by gas self-shielding. Indeed, the gas self-shielding threshold at these high redshifts is about 20 times above the cosmic mean (Chardin et al., 2018), indicating that the higher density halo gas is exposed to a less intense UV background than the IGM. The halo gas may therefore require more time to be heated up.

Compared to the fiducial-RT run, there is more suppression of the halo baryon content in the fiducial-UVB run across all redshifts, lowering the baryon fraction of $10^{8}\ M_{\odot}$ halos by $75\%$, $88\%$, and $92\%$ at $z=7,6,5$, respectively. After the UVB is turned on at $z\approx 10.7$, it quickly heats up all the gas in the simulation volume and thus acts as an early reionization model in terms of its effect on low-mass halos. Using a UVB model that completes reionization at $z\approx 6$ with a more realistic thermal history of the IGM (e.g. Oñorbe et al., 2017; Puchwein et al., 2019) will lead to similar amounts of suppression of the halo baryon content as the fiducial-RT run, which will be demonstrated in Sec. 4.

Turning off stellar feedback leads to photoheating being able to generate much more suppression of the baryon content of the $M_{\rm vir}\lesssim 10^{9}\ M_{\odot}$ halos. At $z=7,6,5$, the $10^{8}\ M_{\odot}$ halos in the NW-RT run undergo a $85\%,91\%,93\%$ depletion in baryon fraction respectively. On the other hand, the amount of suppression of baryon fraction in the NW-UVB run is the same as in the fiducial-UVB run because the strength of the UVB is independent of the stellar feedback model. Photoevaporation due to feedback from internal UV photons is thus regulated by the intensity of the UV field. A higher photoionization rate leads to a lower H i fraction, resulting in less cooling in the dense halo gas because the H i fraction controls cooling at temperatures of $10^{4}-10^{5.5}$ K (Ocvirk et al., 2018a). This leads to the possibility that using the actual speed of light may strengthen the suppression of the baryon fraction. However, the photoionization rate is not affected by the reduced speed of light before the overlap of ionized bubbles (Ocvirk et al., 2018a), implying that the choice of the speed of light plays a minor role in determining the suppression of baryon content. We will discuss effects of the reduced speed of light approximation in detail in Sec. 4. On the other hand, feedback from the external radiation field is less affected by the UVB amplitude, because the IGM thermal evolution is independent of it (e.g. Hui & Gnedin, 1997; McQuinn & Upton Sanderbeck, 2016).

The suppression of halo baryon content translates to a similar trend in the suppression of the instantaneous SFR. Fig. 6 presents the median SFR–halo mass relations at $z=7,6,5$ in simulations without RT (black lines), with FG09 UVB (blue lines), and with RT (red lines), performed with both the fiducial (top panels) and NW (bottom panels) models. Errorbars, representing the $1\sigma$ scatter, are only shown for the fiducial-RT run. At $z=7$, no suppression of the halo SFR is seen in the fiducial-RT run. At $z=6$, the SFR–halo mass relation of the fiducial-RT run starts deviating from that of the fiducial-noRT run at about $10^{8.4}\ M_{\odot}$. In contrast, suppression of SFR at $z=7$ and $6$ is seen in the fiducial-UVB run in halos as massive as $\sim 10^{9}\ M_{\odot}$. The fiducial-RT run thus only generates a small decrease in the SFR–halo mass relation444In fact the convergence test in Appendix A will show that the L3n256 RT run likely shows no suppression at all at $z=6$.. At $z=5$ the suppression of SFR in low-mass halos in the fiducial-RT and fiducial-UVB run are comparable. $50\%$ suppression of SFR happens at $M_{\rm vir}\sim 10^{8.8}\ M_{\odot}$ in both simulations, consistent with the findings of previous works (e.g. Finlator et al., 2011; Ocvirk et al., 2016). However, the lack of suppression of SFR in the fiducial-RT run at the end of reionization ($z=6$) is in tension with these previous studies. Since the strength of photoheating feedback relies on the total amount of radiation sources, the extent of reduction in star formation by photoheating feedback is SFR dependent. The no RT simulation in Finlator et al. (2011) produces an SFR of about $10^{-2.6}\ M_{\odot}/{\rm yr}$ in $10^{8.5}\ M_{\odot}$ halos at $z=6$, which is an order of magnitude higher than our fiducial-noRT simulation ($10^{-3.6}\ M_{\odot}/{\rm yr}$). The higher SFR in their simulations leads to a larger radiation field intensity, allowing photoheating feedback to be more effective at an earlier time. In our simulations stellar feedback is strong enough to suppress star formation efficiently, producing a lower impact of photoheating feedback. Moreover, the recently updated simulations of Ocvirk et al. (2018b) show much less quenching in $10^{8}-10^{9}\ M_{\odot}$ halos compared to the simulations of Ocvirk et al. (2016) after a recalibration of the star formation sub-grid model, confirming our analysis.

Additionally, both Ocvirk et al. (2016) and Finlator et al. (2011) may suffer from insufficient resolution due to large grid sizes. While the $10^{8}-10^{9}\ M_{\odot}$ halos have virial radii of $10-20$ ckpc$/h$, the gas cells of the Ocvirk et al. (2016) simulations are of $\sim 16$ ckpc$/h$ in width. The finest RT grid in the simulations of Finlator et al. (2011) is even coarser, having a side length of $\sim 93$ ckpc$/h$, so one grid cell could cover an entire low mass halo. A degraded RT grid can smooth out the inhomogeneity of the ionizing background. Thus their RT simulations can possibly mimic the uniform UVB on suppressing halo SFR, making the effect of photoheating feedback seemingly stronger. This shows the importance of resolving these low mass halos both spatially and in mass.

Turning off stellar feedback results in a much larger suppression of SFR in low mass halos due to photoheating feedback, similar to the findings in the baryon fraction–halo mass relation. At $z=7$, the SFR–halo mass relation in the NW-RT run starts deviating from that of the NW-noRT run at $M_{\rm vir}\lesssim 10^{8.7}\ M_{\odot}$. At $z=6$ and $z=5$ the majority of the $M_{\rm vir}\lesssim 10^{9}\ M_{\odot}$ halos are completely quenched. The NW-RT run also produces more suppression of SFR than the NW-UVB run at all times because of the increased ionizing radiation intensity. Our findings imply that the strength of photoheating feedback is suppressed by stellar feedback, consistent with Finlator et al. (2011) but in contrast with Pawlik et al. (2015). We will come back to the interplay between stellar feedback and photoheating feedback in Sec. 4 and discuss how this depends on the implementation of the galactic wind.

3.3 UV Luminosity Function (UVLF), Stellar Mass Function, Cosmic SFR Density

Since the baryon mass fraction–halo mass relation and SFR–halo mass relation are not directly observable, we evaluate how photoheating feedback shapes the more directly observable quantities, including the UVLF, stellar mass function, and cosmic SFR density. The simulated UVLFs are calculated by the following procedure. For each star particle in the Subfind subhalos, we compute its rest-frame spectrum by interpolating the Flexible Stellar Population Synthesis (FSPS) library with nebular emission (Conroy et al., 2009; Conroy & Gunn, 2010) based on its age and metallicity. We do not include dust extinction since it has negligible impact on the UVLF for UV magnitudes $\gtrsim-18$ mag (e.g. Tacchella et al., 2018). The rest-frame spectrum of each galaxy is then the summation of the spectra of its star particles. The rest-frame 1500 Å luminosity is obtained by convolving the galaxy’s spectrum with a top-hat filter centered at 1500 Å with 400 Å in width.

Fig. 7 illustrates the simulated UVLFs in the no RT (black lines), FG09 UVB (blue lines), and RT (red lines) runs for both the fiducial (top panels) and NW (bottom panels) models at $z=7,6,5$, compared with the observations of Bouwens et al. (2015) (green crosses). In the top panels, observational data from Bouwens et al. (2017) ($z=6$, magenta triangles), Livermore et al. (2017) ($z=7$ and $z=6$, yellow squares), and Atek et al. (2018) ($z=6$, gray stars) are also shown. Errorbars, representing the $1\sigma$ scatter, are only plotted for the fiducial-RT run. We cut off the UVLFs at $-15$ mag, because for higher magnitudes the UVLFs will exhibit a turnover caused by lack of resolution (see Appendix A for a demonstration). The UVLFs in the fiducial simulations roughly match the observational data at these redshifts, proving the ability of the Illustris galaxy formation model to reproduce high redshift observations. There is no observable bend in the faint end slope of the UVLF in the fiducial-RT and fiducial-UVB simulations at ${\rm M_{1500}}<-15$ mag compared to the fiducial-noRT run. Since the most massive galaxy in a $10^{9}\ M_{\odot}$ halo in the fiducial simulations is about $10^{6}\ M_{\odot}$ in stellar mass, which has a UV magnitude of about $-14$ mag, an observable flattening of the faint end slope of the UVLF is more likely to be seen at ${\rm M_{1500}}\gtrsim-14$ mag. In principle we can combine the UVLFs from our L3n256 run to push the simulated UVLFs to lower luminosities, but we refrain to do so because the 3 cMpc$/h$ boxes suffer more from cosmic variance. We therefore conclude that the addition of photoheating feedback from reionization does not induce an observable difference in the faint end slope of the $z>5$ UVLF at ${\rm M_{1500}}<-15$ mag.

In the NW simulations, the higher level of star formation raises the UV luminosity of a halo of a given mass. Suppression of SFR in the $\lesssim 10^{9}\ M_{\odot}$ in the NW-RT run is therefore reflected in a suppression of the faint end slope of the UVLF at ${\rm M_{1500}}\gtrsim-16.5$ mag at $z=6$ and $5$ compared to the NW-noRT run, which does not exist in the fiducial runs.

Fig. 8 presents the simulated stellar mass functions at $z=7,6,5$ compared with the measurements of Song et al. (2016) (green circles), which the fiducial runs roughly reproduce555We note that at high redshifts, the stellar mass function estimate from observations is very uncertain, due to limited sample size and systematic uncertainties in the modeling of galaxy SEDs (see e.g. Tacchella et al., 2018, for a detailed discussion).. The fiducial-RT and fiducial-UVB runs do not generate any observable suppression of the abundance of low-mass galaxies down to $\sim 10^{6}\ M_{\odot}$ compared to the fiducial-noRT run. As pointed out earlier, the most massive galaxy in a $10^{9}\ M_{\odot}$ halo has about $M_{*}=10^{6}\ M_{\odot}$ in the fiducial simulations. Hence, the simulated stellar mass function is not expected to show much change at $\gtrsim 10^{6}\ M_{\odot}$ when photoheating is included. Contrarily, due to the higher stellar mass of galaxies in a given halo in the NW simulations, a suppression of the number of low-mass galaxies is seen at $10^{7}\ M_{\odot}$ at $z=7$ in the NW-RT and NW-UVB runs compared to the NW-noRT run, which gets stronger with time. Thus, in the dynamic range that we are able to probe in the simulations, no observable difference is seen in the faint end of the UVLF or the low-mass end of the stellar mass function unless stellar feedback is turned off.

Finally, Fig. 9 illustrates the cosmic SFRD as a function of redshift in the fiducial (top panel) and NW (bottom panel) simulations. The cosmic SFRD in the fiducial runs does not match the observations of Bouwens et al. (2015) (green triangles), mainly because of the lack of bright sources in the small simulation volume that lowers the cosmic SFRD at late times. Interestingly, there is no observable dip in the cosmic SFRD in either the fiducial-UVB or the fiducial-RT run. Even in the NW-UVB and NW-RT simulations, where the $\lesssim 10^{9}\ M_{\odot}$ halos are largely quenched due to photoheating feedback, there is no drop in the cosmic SFRD during or after reionization, in contrast with the prediction of Barkana & Loeb (2000). The cosmic SFRD may experience a fall-off if it is dominated by halos of masses $\lesssim 10^{9}\ M_{\odot}$, which indicates that the reduction of SFR in these halos is not reflected in the cosmic SFRD. The magenta dashed lines in both panels of Fig. 9 represent the star formation histories of three halo mass bins in the fiducial-RT and NW-RT runs, respectively: $10^{8}-10^{9}\ M_{\odot}$, $10^{9}-10^{10}\ M_{\odot}$, and $>10^{10}\ M_{\odot}$. Regardless of whether stellar feedback is included, the $10^{8}-10^{9}\ M_{\odot}$ halos dominate the cosmic SFRD before $z\approx 9.5$. At $z\sim 9.5-6.5$ and $z\lesssim 6.5$, the major contribution to the cosmic SFRD comes from the intermediate mass halos and the most massive halos, respectively. Therefore, the dominance of the cosmic SFRD by halos that are not affected by photoheating feedback during and after reionization compensates for the suppressed star formation in the low-mass halos. Our results suggest that it is unlikely that reionization can be probed by an observable dip in the evolution of the cosmic SFRD.

4 Discussions

4.1 Interplay Between Photoheating And Stellar Feedback

Results in Sec. 3 indicate that stellar feedback is able to suppress the strength of photoheating feedback, because the latter generates much more suppression of halo baryon content and SFR when the former is turned off. In this section we examine in detail the non-linear coupling between the two feedback processes.

To better understand the efficacy of the two feedback mechanisms, we define the amplitude of suppression of SFR due to stellar feedback only as

[TABLE]

and suppression amplitude caused only by photoheating as

[TABLE]

The inclusion of both photoheating and stellar feedback gives a suppression amplitude of

[TABLE]

These definitions are the inverse of those in Pawlik & Schaye (2009), but we calculate the suppression amplitudes in this way to avoid division by zero. Fig. 10 shows $s_{\rm w}$ (green), $s_{\rm h}$ (red), $s_{\rm wh}$ (black) as a function of halo mass at $z=7,6,5$ in panels from left to right, respectively. The suppression amplitude $s_{\rm h}$ drops below $\sim 0.01$ at halo masses $\lesssim 10^{8.4}\ M_{\odot}$, $\lesssim 10^{8.6}\ M_{\odot}$, and $\lesssim 10^{8.8}\ M_{\odot}$ at $z=7,6,5$, respectively. For higher mass halos ($\gtrsim 10^{9}\ M_{\odot}$), the power of photoheating feedback quickly fades away, with the value of $s_{\rm h}$ rapidly increasing to $1$. Photoheating feedback therefore only suppresses star formation in halos of $10^{8}-10^{9}\ M_{\odot}$666We will show in Appendix A that this sharp drop in SFR is indeed caused by photoheating feedback, not due to lack of resolution.. Stellar feedback, in contrast, is able to reduce star formation across the entire halo mass range of $10^{8}-10^{11}\ M_{\odot}$ by a factor of $\sim 20$, making it the dominant mechanism in regulating star formation.

Stellar feedback and photoheating feedback do not seem to amplify the effect of each other when coupled together, contrary to the findings of Pawlik & Schaye (2009) and Pawlik et al. (2015). While halos less massive than $\sim 10^{8.4}\ M_{\odot}$ and $\sim 10^{8.6}\ M_{\odot}$ are quenched in the NW-RT run at $z=7$ and $6$ respectively, halos in the same mass range in the fiducial-RT run are still forming stars. The value of $s_{\rm wh}$ in these halo mass ranges at $z=7$ and $6$ is at most $\sim 0.1$ dex lower than $s_{\rm w}$. If photoheating and stellar feedback boost the power of each other, we should get $s_{\rm wh}<s_{\rm w}s_{\rm h}$, which is not seen. This demonstrates that the strength of photoheating feedback is weakened by stellar feedback, as found in Sec. 3. The driving force of this effect is the large difference in the suppression amplitudes $s_{\rm w}$ and $s_{\rm h}$. As stellar feedback dominates the regulation of star formation, it reduces the strength of the radiation field, thus suppressing the impact of photoheating feedback. This effect mainly concerns the internal photoheating feedback, because the external photoheating feedback is less affected by changes in the UVB intensity (see Sec. 3). This likely causes the lack of internal photoheating feedback found in Sec. 3.

The major reason why we see a different interplay of the two feedback mechanisms from Pawlik & Schaye (2009) and Pawlik et al. (2015) lies in the galactic wind scheme. Ionizing radiation from new born stars heats up the surrounding medium and decreases its density, hence reducing the thermal losses that the wind undergoes after the SNe go off (Stinson et al., 2013; Hopkins et al., 2014; Rosdahl et al., 2015; Kannan et al., 2018). In local feedback implementations where the SN thermal energy is released into the adjacent gas cells of the star particle, SN feedback works more efficiently when photoheating feedback is included (e.g. Hasegawa & Semelin, 2013). By decoupling the wind particles from hydrodynamic forces at n\,\lower 3.01385pt\hbox{;\stackrel{{\scriptstyle\textstyle>}}{{\sim}};}0.01 cm*-3*, the wind no longer suffers from thermal losses in the high density ISM gas. Effects of photoheating feedback and SN feedback are thus “decoupled” in the high density ISM, preventing the boost that photoheating might have on the stength of SN feedback. This also contributes to the insufficient suppression of halo baryon content and SFR at the end of reionization in the fiducial-RT run. Another reason for the disagreement with Pawlik & Schaye (2009) is that they used the same UVB (Haardt & Madau, 2001) in simulations with and without stellar feedback, while the UV radiation field can get stronger in the absence of stellar feedback. A larger UVB intensity heats up the dense halo gas more, leading to larger suppression of star formation. We therefore conclude that how stellar feedback and photoheating feedback affect each other is model-dependent.

However, the sub-dominance of photoheating feedback in regulating star formation compared to stellar feedback is in agreement with many other works using different galaxy formation prescriptions. Our findings are consistent with those of Pawlik et al. (2015) that stellar feedback plays the dominant role in shaping the galaxies properties, and that photoheating does not leave detectable imprints on the UVLF. Moreover, Rosdahl et al. (2018) showed that there is little change in the SFR–halo mass relation when switching from using the single star SED model to binary star SED model. The former failed to complete reionization by $z=6$ in their work, while the latter did. This supports the idea that radiation feedback is sub-dominant in suppressing star formation compared to stellar feedback. Our results are also in agreement with the semi-analytic models of Wyithe & Loeb (2013) and Mutch et al. (2016), that stellar feedback plays a greater regulatory role than photoheating. Our predictions on the effects of photoheating on the observables are thus robust. The simulated IGM clumping factor is also relatively robust because the non-local wind scheme does not affect this low-density regime (see Appendix C).

4.2 The Reduced Speed of Light Approximation

As discussed in Sec. 3, using the actual speed of light should boost the post-reionization photoionization rate by a factor of $\sim 10$, which possibly leads to more suppression of star formation after reionization. Indeed, Appendix B will show that the post-reionization temperatures of gas cells with overdensities of $10-1000$ are lowered by $5,000-10,000$ K when using a reduced speed of light of $0.1c$. By comparing the fiducial-UVB run to a simulation using a scaled version of the FG09 UVB where the photoionization and photoheating rates are raised by a factor of $10$, we found that a $\sim 10,000$ K difference in the halo gas temperature results in a $0.1\sim 0.2$ dex change in the halo mass threshold where SFR suppression begins to show up. This implies that the post-reionization SFR suppression is not much affected by the adoption of reduced speed of light approximation. On the other hand, the gas temperatures are well-converged when $\langle x_{\rm\text{H\,{i}}}\rangle^{V}$ is above $\sim 0.01$, so the SFR suppression before the overlap of ionized bubbles is also not influenced by the reduced speed of light approximation. Our findings on the SFR suppression are therefore robust to the choice of the reduced speed of light.

4.3 Other UVB Models

We explore whether using a UVB model with a more realistic reionization history can generate similar trends in the baryon fraction and SFR suppression as the fiducial-RT simulation. For instance, the Puchwein et al. (2019) UVB model is designed to complete reionization at $z\approx 6$ and generates an IGM thermal history with a peak at $z\approx 6$. We hence perform an additional fiducial simulation with their “equivalent-equilibrium” photoionization and photoheating rates777https://arxiv.org/src/1801.04931v1/anc/TREECOOL. Fig. 11 presents a comparison of the baryon fraction–halo mass relation (top panel) and SFR–halo mass relation (bottom panel) at $z=7,6,5$ (colors from light to dark), with the fiducial-RT and Puchwein et al. (2019) UVB simulations shown in red and teal, respectively. This UVB model generates similar suppression of the halo baryon fraction and SFR as the fiducial-RT simulation at all redshifts, thus providing comparable external photoheating feedback as the fiducial-RT run.

5 Conclusions

In this work we present a suite of state-of-the-art cosmological radiation hydrodynamic simulations with arepo-rt using the Illustris galaxy formation model to simulate the process of reionization. We examined the effects of photoheating feedback due to reionization on galaxy properties and compared the impact of photoheating feedback with that of stellar feedback. Our main results are listed as follows.

• •

Reionization completes at $z\approx 6$ in the fiducial-RT run. The simulation is able to match the observed volume-averaged neutral hydrogen fraction at $z=5-6$ (Fan et al., 2006), the intensity of the post-reionization ionizing background (Calverley et al., 2011; Wyithe & Bolton, 2011), the cumulative optical depth to Thomson scattering (Planck Collaboration et al., 2018), the UVLFs (Bouwens et al., 2015) and stellar mass functions (Song et al., 2016) at $z=5-7$, although for the first two there is a dependence on the choice of the reduced speed of light at fixed escape fraction. This demonstrates the ability of our RHD scheme to simulate a realistic reionization process, and the capability of the Illustris galaxy formation model to reproduce high redshift observations.

• •

At $z=6$ ($z=5$), suppression of the baryon content and SFR of low-mass halos ($10^{8}-10^{9}\ M_{\odot}$) due to photoheating feedback only begins to be seen at $\lesssim 10^{8.4}\ M_{\odot}$ ($10^{8.8}\ M_{\odot}$) in our fiducial-RT run, indicating insufficient internal photoheating feedback from photons in the same halo. However, turning off stellar feedback leads to quenching of these low-mass halos at $z\leq 7$. The FG09 UVB acts as an early reionization model and begins suppressing star formation earlier. This discrepancy can be mitigated by using a UVB model with a more realistic reionization history (e.g. Puchwein et al., 2019).

• •

Photoheating does not generate any observable flattening in the faint-end slope of the UVLFs up to ${\rm M_{1500}}=-15$ mag in the fiducial simulations, or of the low-mass end of the stellar mass functions down to $10^{6}\ M_{\odot}$. However, we point out that there may be an observable difference in the faint-end slope of the UVLF if one can probe down to ${\rm M_{1500}}\gtrsim-14$ mag. We also did not see any dip in the cosmic SFRD during or after reionization, because the SFRD is dominated by halos more massive than $10^{10}\ M_{\odot}$ which are not affected by photoheating at $z\lesssim 6.5$ near the end of reionization.

• •

Photoheating quenches star formation in low-mass halos with masses $\lesssim 10^{9}\ M_{\odot}$ at $z\gtrsim 5$ without the presence of stellar feedback. Its effect on higher mass halos is negligible. On the contrary, stellar feedback is able to reduce star formation across the entire sampled halo mass range by a factor $\sim 20$. When coupled together, stellar feedback suppresses the strength of photoheating feedback by reducing the amount of radiation sources. This interplay between the two feedback mechanisms is a result of the non-local galactic wind scheme, but the dominance of stellar feedback in regulating star formation is consistent with other works using different galaxy formation models.

In addition to the impact of photoheating being weakened by stellar feedback, another likely cause of the lack of suppression of halo SFR at the end of reionization in the fiducial-RT run is the diversity in the reionization times of halos of different masses. If low-mass halos are exposed to the bulk of the ionized bubbles later than the most massive ones, there may be a delay in the response of their SFR to the reionization process. The evolution of the halo baryon fraction already hints upon this hypothesis (see Sec. 3). However, this cannot be checked in our current simulations because we did not include tracer particles that track the reionization time of each gas cell. We defer this analysis and test this scenario in future work.

Future observational facilities, especially the James Webb Space Telescope (JWST), will be able to observe a number of high redshift galaxies, thus offering new insights into the sources that reionized the Universe. Some of the deepest JWST surveys in the first two years can provide a complete sample of galaxies with ${\rm M_{1500}}\lesssim-17$ mag (Williams et al., 2018), indicating the need for deeper surveys to explore the $z\gtrsim 5$ UVLFs at ${\rm M_{1500}}\gtrsim-14$ mag. However, if photoheating only introduces a $\sim 0.05$ change in the faint end slope of the UVLF (Gnedin & Kaurov, 2014), detecting imprints of the reionization process from the faint end slope of the UVLF seems questionable. Moreover, the dominance of $\gtrsim 10^{10}\ M_{\odot}$ halos on the cosmic SFRD at $z\lesssim 6.5$ also makes it unlikely to detect an observable dip in the cosmic SFRD during reionization. More careful investigation is therefore needed to explore the feasibility of using these observables to explore reionization.

Acknowledgements

We thank Romeel Davé and Hui Li for valuable discussions regarding this work, and Sandro Tacchella for generous feedback on the draft of this paper. We also thank Daisuke Nagai, Nick Gnedin, and Kristian Finlator for useful comments on some of the key results of this work. We thank Joop Schaye, Pierre Ocvirk, Stuart Wyithe, Ewald Puchwein, Pratika Dayal, and Steve Finkelstein for important feedback on the first version of this paper. RK acknowledges support from NASA through Einstein Postdoctoral Fellowship grant number PF7-180163 awarded by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS8-03060. FM is supported by the program “Rita Levi Montalcini” of the Italian MIUR. MV acknowledges support through an MIT RSC award, a Kavli Research Investment Fund, NASA ATP grant NNX17AG29G, and NSF grants AST-1814053 and AST-1814259. The simulations were performed on the Harvard computing cluster supported by the Faculty of Arts and Sciences.

Appendix A Numerical convergence of the halo SFR and UV luminosity function

In this section we discuss the numerical convergence of the of our simulations by using the L6n256 L3n256 boxes, both run with the fiducial model.

We generated three different initial conditions (ICs) for the L3n256 runs using three different random number seeds. An escape fraction of 0.7 is adopted for running all these ICs. Due to the larger cosmic variance of a 3 cMpc$/h$ box, only one of the ICs is able to generate a hydrogen reionization history that experiences complete overlap of ionized bubbles at $z\approx 6$. Fig. 12 shows the resulting volume-averaged hydrogen ionization fraction as a function of redshift of this IC (red line) compared with that of the L6n256 fiducial-RT run (blue). We will focus the analysis of the remaining of this section on results from this particular IC. The other two ICs give reionization histories that end at $z<5$, because in general L3n256 requires a higher escape fraction to finish reionization at the same redshift as L6n256, due to the lack of bright sources responsible for completing reionization (Katz et al., 2018) and the ability to resolve more small scale clumping.

Fig. 13 shows the median SFR as a function of halo mass extracted from L3n256 (dashed lines) and L6n256 (solid lines) at $z=7,6,5$ (from left to right panels). Black, blue, and red lines represent no RT, with FG09 UVB, and with RT respectively. The L6n256-noRT SFR–halo mass relations start turning downwards from the L3n256-noRT ones at $\sim 10^{8.2}\ M_{\odot}$, $\sim 10^{8.5}\ M_{\odot}$, and $\sim 10^{8.7}\ M_{\odot}$ at $z=7,6,5$, respectively. The deviations of the L6n256-UVB run from L3n256-UVB happen at larger halo masses because of the quenching by photoheating. The L3n356-RT run at $z=6$ basically shows no suppression of SFR compared to the L3n256-noRT run, due to the slightly later overlap of ionized bubbles than L6n256-RT. Based on this convergence study, we are more inclined to conclude that at $z=6$ there is little or no suppression of SFR due to photoheating by RT. Nevertheless, at $z=5$ we find good agreement between our L3n256 and L6n256 RT and UVB runs, in terms of the position of 50% suppression of SFR (halo mass $\sim 10^{8.8}\ M_{\odot}$). Therefore although the low-mass halos may not be completely quenched in the L3n256-RT and L3n256-UVB runs, the suppression of SFR in the L6n256 simulations is indeed a photoheating effect, and it is not caused by lack of resolution.

Fig. 14 compares the UVLFs from the L6n256-RT (blue) and L3n256-RT (red) runs at $z=7,6,5$ (from left to right). Comparisons between the no RT or UVB runs are similar. Since the L3n256 UVLFs suffer more from stochasticity, we do not try to combine the UVLFs from the two simulation volumes to get a large dynamical range in UV luminosity. At ${\rm{\rm M_{1500}}}\gtrsim-15$ mag, the L6n256 UVLFs experience a turnover compared to the L3n256 ones due to lack of resolution. Increasing the mass resolution by a factor of 8 generates more star particles to sample the star formation history, raising the faint end of the UVLF. We thus cut off the UVLFs at -15 mag in Fig. 7.

Appendix B Numerical convergence of the reduced speed of light approximation

To assess the possible effects of using the reduced speed of light on the suppression of halo baryon fraction and SFR, we have performed simulations of $25$ cMpc$/h$ box size with $2\times 256^{3}$ resolution elements (L25n256) with $0.1c$, $0.3c$, and $1.0c$. All three simulations are run with $f_{\rm esc}=1$. The mass resolution of L25n256 does not allow us to probe the suppression of SFR in $\lesssim 10^{9}\ M_{\odot}$ halos directly, but these simulations provide information about how the gas temperature changes with the adopted speed of light. The top panel of Fig. 15 shows the resulting $\langle x_{\rm\text{H\,{i}}}\rangle^{V}$ evolution with $z$. Simulations using $0.1c$, $0.3c$, and $1.0c$ are represented by black, blue, and red lines, respectively. The reionization histories are well-converged before $\langle x_{\rm\text{H\,{i}}}\rangle^{V}$ drops to $\sim 0.01$. Using $0.3c$ gives good convergence in terms of the time of reaching $\langle x_{\rm\text{H\,{i}}}\rangle^{V}=10^{-4}$, but adopting $0.1c$ delays this redshift by $\sim 0.5$. The post-reionization $\langle x_{\rm\text{H\,{i}}}\rangle^{V}$ scales as the inverse of the value of the reduced speed of light. These results are consistent with the findings of Deparis et al. (2019) and Ocvirk et al. (2018a).

The bottom panel of Fig. 15 illustrates the temperature evolution of gas with overdensities of $1$ (solid), $10$ (dashed), $100$ (dot-dashed), and $1000$ (dotted). The temperature of the IGM does not depend on the amplitude of the post-reionization UVB, so $T(\bar{\rho})$ is well-converged. $T(1000\bar{\rho})$ is also converged at all redshifts, consistent with the findings of Ocvirk et al. (2018a) that $x_{\rm\text{H\,{i}}}$ is converged at overdensities $\gtrsim 1000$ after reionization. Before $\langle x_{\rm\text{H\,{i}}}\rangle^{V}$ drops to $\sim 0.01$, gas temperature at all overdensities are converged. After that point, using $0.1c$ can lead to gas temperatures at overdensities of $10$ and $100$ being underestimated by $5,000-10,000$ K.

In order to understand how this $5,000-10,000$ K underestimation in the dense gas temperature can affect the suppression of star formation, we ran an additional L6n256 simulation with the FG09 UVB, but with the photoionization and photoheating rates scaled by a factor of $10$ (denoted by FG09x10). This mimics the effects of using the actual speed of light, especially after reionization. Fig. 16 shows the temperature–density diagrams at $z=6$ in the original FG09 simulation (left) and the FG09x10 simulation (right). The IGM temperature is unchanged, as expected. Gas with overdensities larger than $\sim 10$ have $\sim 10,000$ K higher temperatures in FG09x10 than FG09. Fig. 17 illustrates the baryon fraction–halo mass relation (top panels) and SFR–halo mass relation (bottom panels) at $z=6$ in the FG09 (blue lines) and FG09x10 (green lines) simulations. The no RT simulation results are shown in black. Despite the $\sim 10,000$ K difference in the halo gas temperature, the $10^{8}\ M_{\odot}$ halos in FG09x10 only experience a $\sim 0.2$ dex more decrease in the baryon content than FG09. The halo mass range of SFR suppression is enlarged by $0.1\sim 0.2$ dex in FG09x10, implying the effect of $\sim 10,000$ K temperature difference is not strong. We therefore conclude that our results on the suppression of halo SFR is relatively robust with the choice of the reduced speed of light.

Appendix C IGM Clumping

We investigate how photoheating and galactic wind reduce the IGM recombination rate by computing the clumping factor $C=\langle\rho^{2}\rangle/\langle\rho\rangle^{2}$. Here we focus on $C_{100}$, which parametrizes the average recombination rate of gas with overdensities $\leq 100$. Fig. 18 shows the evolution of $C_{100}$ with redshift in the fiducial (solid lines) and NW (dashed lines) simulations. Black, red, and blue represent no RT, RT, and UVB runs, respectively. Photoheating strongly decreases the clumping factor by increasing the Jeans mass of the ionized gas. Galactic wind raises the clumping factor because it blows gas out of galaxies into the IGM. The fiducial-UVB run also generates much lower $C_{100}$ than the fiducial-RT run because the FG09 UVB turns on at a high redshift. These results are consistent with the findings of Pawlik et al. (2015), although the NW-RT run produces a much lower $C_{100}$ than their corresponding simulation because of a much earlier reionization. The effects of photoheating are also qualitatively similar to what was found by Finlator et al. (2012), but their galactic wind model does not seem to move gas out of galaxies as efficiently as ours. We therefore conclude that the simulated effects of photoheating on the IGM properties are relatively robust.