A Cosmic UV/X-ray Background Model Update

TL;DR

This paper updates the cosmic UV/X-ray background model to better match recent empirical data, revealing the dominant sources of ionization at different redshifts and implications for galaxy escape fractions.

Contribution

The new model incorporates recent observational constraints and improves upon previous models, providing a more accurate depiction of the cosmic ionizing background across redshifts.

Findings

AGN dominate the HI ionizing background at z<~3.

Star-forming galaxies dominate at higher redshifts.

An escape fraction of ~1% from galaxies suffices for reionization at z=3.

Abstract

We present an updated model of the cosmic ionizing background from the UV to the X-rays. Relative to our previous model (Faucher-Giguere et al. 2009), the new model provides a better match to a large number of up-to-date empirical constraints, including: 1) new galaxy and AGN luminosity functions; 2) stellar spectra including binary stars; 3) obscured and unobscured AGN; 4) a measurement of the non-ionizing UV background; 5) measurements of the intergalactic HI and HeII photoionization rates at z~0-6; 6) the local X-ray background; and 7) improved measurements of the intergalactic opacity. In this model, AGN dominate the HI ionizing background at z<~3 and star-forming galaxies dominate it at higher redshifts. Combined with the steeply declining AGN luminosity function beyond z~2, the slow evolution of the HI ionization rate inferred from the high-redshift HI Lya forest requires an escape fraction from star-forming galaxies that increases with redshift (a population-averaged escape fraction of ~1% suffices to ionize the intergalactic medium at z=3 when including the contribution from AGN). We provide effective photoionization and photoheating rates calibrated to match the Planck 2018 reionization optical depth and recent constraints from the HeII Lya forest in hydrodynamic simulations.

Figures

Equations

$$\left(\frac{\partial }{\partial t}-\nu H\frac{\partial }{\partial \nu }\right)J_{\nu }=-3H{J}_{\nu }-c{\alpha }_{\nu }{J}_{\nu }+\frac{c}{4\pi }{ϵ}_{\nu },$$

$${ϵ}_{\nu }={ϵ}_{\nu }^{\star }+{ϵ}_{\nu }^{\mathrm{AGN}}+{ϵ}_{\nu }^{\mathrm{rec}}.$$

$$J_{{\nu }_0}(z_0)=\frac{1}{4\pi }{\int }_{z_0}^{\infty }dz\frac{dl}{dz}\frac{(1+z_0{)}^3}{(1+z{)}^3}{ϵ}_{\nu }(z)\exp [-\bar{\tau }({\nu }_0,z_0,z)],$$

$$\bar{\tau }({\nu }_0,z_0,z)={\int }_{z_0}^{z}d{z}^{\prime }{\int }_0^{\infty }d{N}_{\mathrm{HI}}\frac{{\partial }^{2}N}{\partial N_{\mathrm{HI}}\partial z^{\prime }}(1-e^{-{\tau }_{\nu }}),$$

$${\tau }_{\nu }=N_{\mathrm{HI}}{\sigma }_{\mathrm{HI}}(\nu )+N_{\mathrm{HeI}}{\sigma }_{\mathrm{HeI}}(\nu )+N_{\mathrm{HeII}}{\sigma }_{\mathrm{HeII}}(\nu ),$$

$$f(N_{\mathrm{HI}})\equiv \frac{{\partial }^{2}N}{\partial z\partial N_{\mathrm{HI}}}=A(N_{\mathrm{HI}},z)N_{\mathrm{HI}}^{-\beta (N_{\mathrm{HI}},z)}(1+z{)}^{\gamma (N_{\mathrm{HI}},z)},$$

$${\Gamma }_i=4\pi {\int }_{{\nu }_i}^{\infty }\frac{d\nu }{h\nu }{J}_{\nu }{\sigma }_i(\nu ),$$

$${\dot{q}}_i=4\pi {\int }_{{\nu }_i}^{\infty }\frac{d\nu }{h\nu }{J}_{\nu }{\sigma }_i(\nu )(h\nu -h{\nu }_i).$$

$${\log }_{10}{ϵ}_{\nu 1500}^{\star ,\mathrm{com}}=(25.62,25.79,25.92,26.04,26.15,26.26,$$

$$26.32,26.42,26.48,26.54,26.61,26.60,$$

$$26.60,26.52,26.30,26.11,25.90,25.67)$$

$$S_{\nu }^{\mathrm{dust}}=S_{\nu }^{\mathrm{intr}}10^{-0.4E(B-V)k(\lambda )},$$

$$S_{\nu}=\left\{\begin{array}[]{ll}S_{\nu}^{\rm dust}(1-f_{\rm esc}^{\star})+S_{\nu}^{\rm intr}f_{\rm esc}^{\star}&\nu<{\rm 1~{}Ry}\\ S_{\nu}^{\rm intr}f_{\rm esc}^{\star}&\nu\geq{\rm 1~{}Ry}\end{array}\right..$$

$$f_{\mathrm{esc}}^{\star }=\min \left(f_{\mathrm{esc}}^{\star ,z=3}{\left(\frac{1+z}{4}\right)}^{2.5},0.2\right).$$

$${ϵ}_{\nu 912}^{\mathrm{AGN},\mathrm{com}}=10^{24.11}(1+z{)}^{5.87}\frac{\exp (0.73z)}{\exp (3.06z)+15.60}.$$

$$F_E^{\mathrm{tot}}=(1-f_{\mathrm{obsc}})F_E^{\mathrm{int}}+f_{\mathrm{obsc}}{F}_E^{\mathrm{obsc}},$$

$$F_{\mathrm{E}}^{\mathrm{obsc}}=\frac{F_{\mathrm{E}}^{\mathrm{tot}}-(1-f_{\mathrm{obsc}}{F}_{\mathrm{E}}^{\mathrm{int}})}{f_{\mathrm{obsc}}}.$$

$$F_{E}^{\rm tot}=A\left\{\begin{array}[]{ll}E^{-\alpha}e^{-E/E_{1}}&2~{}{\rm keV}\leq E<E_{0}\\ B(1+kE^{\beta-\gamma})E^{-\beta}&E\geq E_{0}\equiv(\beta-\alpha)E_{1}.\end{array}\right.$$

$${\Gamma }_{\mathrm{HI}}{n}_{\mathrm{HI}}=C_{\mathrm{HI}}{\alpha }_{\mathrm{HII}}(T_{IGM0,\mathrm{HI}})n_{\mathrm{HII}}{n}_{\mathrm{e}},$$

$${\Gamma }_{\mathrm{HI}}^{\mathrm{rei}}(z)=\frac{x_{\mathrm{II}}^2(z)}{1-x_{\mathrm{II}}(z)}(1+\chi )C_{\mathrm{HI}}{\alpha }_{\mathrm{HII}}(T_{IGM0,\mathrm{HI}})n_{\mathrm{H}}(z),$$

$${\Gamma }_{\mathrm{HeI}}^{\mathrm{rei}}(z)=\frac{{\alpha }_{\mathrm{HeII}}(T_{IGM0,\mathrm{HI}})}{{\alpha }_{\mathrm{HII}}(T_{IGM0,\mathrm{HI}})}{\Gamma }_{\mathrm{HI}}^{\mathrm{rei}}(z).$$

$${\Gamma }_{\mathrm{HeII}}{n}_{\mathrm{HeII}}=C_{\mathrm{HeII}}{\alpha }_{\mathrm{HeIII}}(T_{IGM0,\mathrm{HeII}})n_{\mathrm{HeIII}}{n}_{\mathrm{e}},$$

$${\Gamma }_{\mathrm{HeII}}^{\mathrm{rei}}(z)=\frac{y_{\mathrm{III}}(z)}{1-y_{\mathrm{III}}(z)}$$

$$\times C_{\mathrm{HeII}}{\alpha }_{\mathrm{HeIII}}(T_{IGM0,\mathrm{HeII}})n_{\mathrm{H}}(z).$$

$$\Theta (z;z_{\mathrm{rei},\mathrm{i}},\Delta z_{\mathrm{rei},\mathrm{i}})=1+\sin \left[\frac{(z_{\mathrm{rei},\mathrm{i}}-z-\Delta z_{\mathrm{rei},\mathrm{i}})\pi }{4\Delta z_{\mathrm{rei},\mathrm{i}}}\right],$$

$$x_{\mathrm{II}}(z)=\Theta (z;z_{\mathrm{rei},\mathrm{HI}},\Delta z_{\mathrm{rei},\mathrm{HI}})x_{\mathrm{II}}^{\mathrm{eq}}(z_{\mathrm{rei},\mathrm{HI}}-\Delta z_{\mathrm{rei},\mathrm{HI}})$$

$$y_{\mathrm{III}}(z)=\Theta (z;z_{\mathrm{rei},\mathrm{HeII}},\Delta z_{\mathrm{rei},\mathrm{HeII}})y_{\mathrm{III}}^{\mathrm{eq}}(z_{\mathrm{rei},\mathrm{HeII}}-\Delta z_{\mathrm{rei},\mathrm{HeII}}),$$

$$\displaystyle\Gamma_{\rm i}^{\rm eff}(z)=\left\{\begin{array}[]{ll}\Gamma_{\rm i}^{\rm eq}(z)&z\leq z_{\rm rei,i}-\Delta z_{\rm rei,i}\\ \Gamma_{\rm i}^{\rm rei,i}(z)&z>z_{\rm rei,i}-\Delta z_{\rm rei,i}.\end{array}\right.$$

$${\frac{dT}{dt}|}_{\mathrm{HI}\mathrm{rei}}=\Delta T_{\mathrm{HI}}\frac{d{x}_{\mathrm{II}}}{dt},$$

$${\dot{q}}_{\mathrm{HI}}{n}_{\mathrm{HI}}$$

$${{\dot{q}}_{\mathrm{tot}}|}_{\mathrm{HI}\mathrm{rei}}{n}_{\mathrm{HI}}$$

$${{\dot{q}}_{\mathrm{tot}}|}_{\mathrm{HI}\mathrm{rei}}$$

Full text

1 Introduction

Galaxies and active galactic nuclei (AGN) produce a diffuse background of ultra-violet (UV) radiation that permeates the intergalactic medium (IGM). This cosmic UV background (UVB) keeps the IGM ionized following the epoch of reionization. UVB models are widely used for two purposes. The first is ionization corrections, which are necessary to convert measured abundances of ions to total masses in different elements (e.g., Stern et al. 2016; Prochaska et al. 2017; Chen et al. 2017). The second is to model the ionization and cooling of cosmic gas in cosmological hydrodynamic simulations (e.g., Cen & Ostriker 1992; Katz et al. 1996; Wiersma et al. 2009; Gnedin & Hollon 2012). AGN also produce a background of X-ray radiation which can ionize heavy elements and further heat the IGM through inverse Compton scattering (Madau & Efstathiou 1999). The cosmic X-ray background (CXB) is directly measured at $z=0$ by X-ray telescopes (e.g., Giacconi et al. 1962; Hickox & Markevitch 2007; Ajello et al. 2008). For simplicity, we will generally refer to the UVB in this paper, but we also include X-rays from AGN.

Several authors have modeled the spectrum and redshift evolution of the UVB (e.g., Miralda-Escude & Ostriker 1990; Haardt & Madau 1996; Giroux & Shapiro 1996; Fardal et al. 1998; Faucher-Giguère et al. 2009). Of these, the series of models by Haardt & Madau (Haardt & Madau 1996, 2001, 2012) have been the most widely used. We refer to the Haardt & Madau (2012) model as HM12. In the last decade, the Faucher-Giguère et al. (2009, hereafter FG09) model has also been used in a number of studies, e.g. in the Illustris and FIRE galaxy formation simulations (Vogelsberger et al. 2014; Pillepich et al. 2018; Hopkins et al. 2014; Hopkins et al. 2018). Recently, Khaire & Srianand (2019, hereafter KS19) produced new synthesis models of the extragalactic background light from the far infrared to the TeV $\gamma-$rays and Puchwein et al. (2019, hereafter P19) updated the HM12 model with an emphasis on the effects of reionization.

Despite the fact that a UVB model is assumed in most cosmological hydrodynamic simulations and that many observational inferences from quasar absorption data require an ionization correction, there is a relative paucity of available models which take into account the most up-to-date empirical constraints on the background intensity and its sources. This not only means that widely-used models may not be consistent with all the latest observational data, but also that researchers are often unable to assess uncertainties in inferences based on these models. For example, a comparison of the ionization rate required by the low-redshift Ly$\alpha$ forest to the HM12 model suggested that known astrophysical sources could not explain the observationally-inferred intergalactic ionization rate, implying a “photon underproduction crisis” (Kollmeier et al. 2014). However, as subsequent studies pointed out, the tension between the Ly$\alpha$ forest observations and UVB models could be alleviated by considering the FG09 model instead, since the latter predicted a higher UVB amplitude at $z\sim 0-0.5$ than HM12 (e.g., Shull et al. 2015; Gurvich et al. 2017, see also Khaire & Srianand 2015). Our new UVB synthesis model, which provides a good match to the most recent measurements of the low-redshift ionization rate, confirms that there is at present no crisis, since most or all of the low-redshift ionization rate can be explained by AGN (see §2.5).

In this paper, we present an update of the FG09 UVB model. Our methods are similar to those described in FG09, but we incorporate several more more up-to-date empirical constraints. Specifically, our new model is informed by:

Recent measurements of the galaxy UV luminosity function out to the epoch of HI reionization and of the total cosmic UV emissivity at low redshifts. 2. 2.

A stellar spectral template including binary stars, as well as recent observational constraints on dust attenuation and the escape of ionizing photons from galaxies. 3. 3.

New measurements of the AGN luminosity function. 4. 4.

An AGN spectral template including both obscured and unobscured sources and constrained to match the local X-ray background. 5. 5.

Improved measurements of the IGM opacity. 6. 6.

Updated constraints from the Ly$\alpha$ forest on the integrated HI photoionization rate, especially at low redshift ($z<0.5$) and approaching the epoch of HI reionization ($z\sim 5-6$). 7. 7.

The Planck 2018 constraint on reionization from the optical depth to the surface of last scattering. 8. 8.

The latest observational constraints on HeII reionization.

Another important motivation for our new calculations is to provide “effective” UVB models designed to address an important issue with the use of spatially homogeneous UVB models in cosmological simulations. Namely, most current simulation codes calculate ionization states under the assumption of optically thin gas in photoionization equilibrium with an homogeneous ionizing background. While this approximation is valid in the IGM when the mean free path of ionizing photons is large, it breaks down before and during reionization events. During reionization, IGM patches are not in photoionization equilibrium but rather transition rapidly from neutral to ionized as ionization fronts propagate. As a result, simulations that use a homogeneous UVB model do not in general correctly model the timing of reionization events and the photoheating that accompanies reionization (e.g., Puchwein et al. 2015).

This issue was briefly discussed in FG09, and was partially addressed by the 2011 release of an updated version of the FG09 UVB model modified to produce a HI reionization redshift $z_{\rm rei,HI}\sim 10$ consistent with the WMAP 7-year optical depth (Komatsu et al. 2011). This model however produces a reionization redshift substantially earlier than the best fit $z_{\rm rei,HI}=7.82\pm 0.71$ from the Planck 2018 results. Relative to earlier Planck results, the 2018 measurement of the electron scattering optical depth to the surface of last scattering ($\tau_{\rm e}=0.054\pm 0.007$) benefited from improved measurements of the large-scale polarization of the cosmic microwave background (CMB; Planck Collaboration et al. 2018).

Recently, Oñorbe et al. (2017) substantially developed the concepts of effective photoionization and photoheating rates and showed explicitly how to modify UVB models so that they produce a correct mean reionization history and more accurate photoheating when used in standard simulation codes. Puchwein et al. (2019) published a modified version of the HM12 UVB model incorporating similar effective rates, but calibrated to produce a higher optical depth $\tau_{\rm e}=0.065$, closer to the Planck 2015 results (Planck Collaboration et al. 2016).

The plan of this paper is as follows. We review our UVB modeling methodology, model ingredients, and compare our new spectral synthesis results to observational constraints in §2. We then derive effective photoionization and photoheating rates calibrated to match desired reionization histories in §3. We discuss our results and conclude in §4. A series of appendices provide additional details, including on how the results depend on model assumptions.

Throughout, we assume a standard flat $\Lambda$CDM cosmology with $\Omega_{\rm m}=0.32$, $\Omega_{\Lambda}=1-\Omega_{\rm m}$, $\Omega_{\rm b}=0.049$, and $H_{0}=67$ km s*-1* Mpc*-1* (e.g., Planck Collaboration et al. 2018). We use $X$ to denote the the hydrogen mass fraction and $Y$ for the helium mass fraction. Lower-case $x$ and $y$ are used to denote ionized fractions of hydrogen and helium, respectively. For example, $x_{\rm II}$ is the hydrogen mass fraction in HII, and $y_{\rm II}$ and $y_{\rm III}$ are the helium mass fractions in HeII and HeIII.

2 UV/X-ray background modeling

2.1 Radiative transfer equations

The angle-averaged specific intensity of the homogeneous background is denoted by $J_{\nu}$ and satisfies the cosmological radiative transfer equation,

[TABLE]

where $H(t)$ is the Hubble parameter, $c$ is the speed of light, $\alpha_{\nu}$ is the proper absorption coefficient per unit length, and $\epsilon_{\nu}$ is the proper emissivity (e.g., Haardt & Madau 1996). We consider three components to the total emissivity, corresponding to star-forming galaxies (superscript $\star$), AGN (superscript AGN), and recombination emission (superscript rec):

[TABLE]

We discuss our assumptions for the different emissivity components in §2.2-2.4.

Integrating equation (1) and expressing the result in terms of redshift gives

[TABLE]

where $\nu=\nu_{0}(1+z)/(1+z_{0})$, the proper line element $dl/dz=c/[(1+z)H(z)]$, and the “effective optical depth” $\bar{\tau}$ quantifies the attenuation due to absorption in the IGM of photons of frequency $\nu_{0}$ at redshift $z_{0}$ that were emitted at redshift $z$ through $e^{\bar{\tau}}\equiv\langle e^{-\tau}\rangle$ (the usual optical depth is related to the absorption coefficient via $d\tau_{\nu}=\alpha_{\nu}dl$). For Poisson-distributed absorbers, each of column density $N_{\rm HI}$,

[TABLE]

where $\partial^{2}N/\partial N_{\rm HI}\partial z^{\prime}$ is the column density distribution of intergalactic absorbers versus redshift (Paresce et al. 1980).

The optical depth $\tau_{\nu}$ shortward of the Lyman limit is dominated by the photoelectric opacity of hydrogen and helium,

[TABLE]

where the $N_{i}$ and $\sigma_{i}$ are the column densities and photoionization cross sections of ion $i$. While the HI column density distribution is reasonably well determined over large redshift and column density intervals, the incidence of HeI and HeII absorbers is not accurately measured. As in previous work, we therefore use a model for the HeI and HeII column densities given the HI column.

We parameterize the HI column density distribution using broken power laws in $N_{\rm HI}$ and $z$:

[TABLE]

where the power-law parameters $A$, $\beta$, and $\gamma$ are functions of the HI column and redshift, and are constrained so that the distribution is everywhere continuous with respect to both $N_{\rm HI}$ and $z$. We adopt a modified version of the broken power-law parameters in Table 1 of Puchwein et al. (2019). The Puchwein et al. (2019) column density distribution is itself is a modified version of the broken power laws used by HM12. This model consists of five column density regimes ($\log{(N_{\rm HI}/{\rm cm^{-2}})}=11-16,~{}16-18,~{}18-19.5,~{}19.5-20.3,$ and $20.3-21.55$) and two redshift regimes ($z\leq 1.56$ and $z>1.56$).

We use the same parameters as in Puchwein et al. (2019) except for the redshift scaling at $z<1.56$. While $f(N_{\rm HI})\propto(1+z)^{0.16}$ at $z\leq 1.56$ in HM12 and Puchwein et al. (2019), we assume a stronger redshift scaling $f(N_{\rm HI})\propto(1+z)$ in this regime (i.e., $\gamma(z\leq 1.56)=1$). The motivation for the stronger redshift is that it is more consistent with the observed redshift evolution of Lyman limit systems (Ribaudo et al. 2011), which most directly determine the mean free path of ionizing photons. To maintain continuity at $z=1.56$, we renormalize the low-redshift $A$ coefficients from Puchwein et al. (2019) by a factor $(1+1.56)^{0.16}/(1+1.56)=0.454$. Appendix A shows that the HI ionizing mean free path implied by this column density distribution is in excellent agreement with observational constraints.

To obtain HeI and HeII columns from $N_{\rm HI}$, we define the ratios $\eta=N_{\rm HeII}/N_{\rm HI}$ and $\zeta=N_{\rm HeI}/N_{\rm HI}$, which are functions of $N_{\rm HI}$ and the photoionization rates. FG09 used radiative transfer calculations to compute $\eta$ and $\zeta$ over a representative range of parameters and improved analytic approximations from Fardal et al. (1998). HM12 further improved on these results by using similar radiative transfer calculations but deriving fitting functions for $\eta$ and $\zeta$ that are more accurate in the limit of optically thick absorbers. In this work, we use the improved fitting functions from HM12. We have verified, however, that we obtain nearly identical radiative transfer results using the $\eta$ and $\zeta$ fitting functions from FG09 instead. This is because, in the limit of optically thick absorbers, the precise values of $N_{\rm HeI}$ and $N_{\rm HeII}$ do not significantly affect the radiative transfer.

For $i\in\{{\rm HI,~{}HeI,~{}HeII}\}$, the photoionization rates are defined as

[TABLE]

where $\sigma_{i}$ is the photoionization cross section for the species of interest and $\nu_{i}$ is the frequency at the photoionization edge. The photoheating rates are defined similarly:

[TABLE]

Finally, we note that we have rewritten our radiative transfer solver from F09 to improve both speed and accuracy. In F09, we solved the radiative transfer equation on a fixed grid of redshifts and frequencies. Because photons redshift as the universe expands, this approach required a very large number of grid points in order to resolve fine structure in rest-frame frequency (e.g., the HeII Ly-series sawtooth). We now instead integrate the solution to the radiative transfer equation (3) along light cones using a recursive method. This method allows us to accurately compute $J_{\nu_{0}}$ at any frequency without requiring a large grid that resolves the full spectrum including all the relevant fine-scale structure at higher redshifts.

2.2 Star-forming galaxies

For the evolution of the emissivity from star-forming galaxies, we assume that the rest-UV comoving emissivity at rest wavelength $\approx 1500$ Å is given by

[TABLE]

at $z=(0,0.2,0.4,0.6,0.8,1,1.2,1.4,1.6,1.8,2,2.1,3,3.8,$ $4.9,5.9,6.8,7.9)$, and interpolate linearly between these redshifts (emissivity in erg s*-1* Hz*-1* Mpc*-3*). At $z\leq 2$, these values correspond to the total UV emissivity as function of redshift inferred by Chiang et al. (2019) by cross-correlating GALEX UV imaging with sources with spectroscopic redshifts in SDSS. The total emissivity at these redshifts is consistent with the emissivity predicted by integrating the galaxy luminosity function down to very faint magnitudes (e.g., Driver et al. 2016) and our results below also indicate that AGN contribute only a small fraction of the total 1500 Å emissivity at these redshifts. At $z\geq 2.1$, the above luminosity densities are inferred from measurements of the rest-UV galaxy luminosity function integrated down to -13.0 AB mag (Bouwens et al., in prep.). Relative to the values reported by Bouwens et al., we have adjusted the luminosity density at some redshifts (within uncertainties) to avoid non-monotonic behavior with redshift.

Figure 1 summarizes the comoving emissivities from star-forming galaxies and AGN adopted in our UVB synthesis model as a function of redshift, at rest-frame wavelengths of 1500 Å and 912 Å. We tie our ionizing background models to UV luminosity densities inferred before any correction for dust obscuration because ionizing photons are also absorbed by dust. Thus, the ionizing emissivity relevant for modeling the intergalactic background should track the emissivity of UV photons that reach the IGM. As Bouwens et al. (2015) show, the dust-corrected and dust-uncorrected UV luminosity densities diverge increasingly from high to low redshift.

Galaxies contribute most importantly to the HI photoionization rate at high redshifts $z\gtrsim 2$, past which the quasar luminosity function begins to drop much more steeply than the nearly flat photoionization rate (e.g., Bolton et al. 2005; Faucher-Giguère et al. 2008a, b; McQuinn et al. 2009). We therefore base our fiducial stellar spectral template on recent results from spectroscopic surveys of $z\sim 3$ Lyman break galaxies (LBGs). While star-forming populations are observed out to much higher redshifts (e.g., Bouwens et al. 2017; Livermore et al. 2017), the stellar populations of epoch-of-reionization galaxies have not yet been studied in as much detail as those at $z\sim 3$. Moreover, absorption by the intergalactic medium makes the direct detection of escaping ionizing radiation from galaxies nearly impossible past $z\sim 3.5$.

For the spectra of star-forming galaxies, we use templates produced by the BPASS stellar population synthesis model, which includes binary stars (v2.2.1; Eldridge et al. 2017). The distributions of binary parameters used in BPASS are as in table 13 of Moe & Di Stefano (2017). Our treatment of stellar binaries improves on previous UVB models, which used either simple power-law approximations for spectra (e.g., FG09) or stellar population synthesis models including single stars only (e.g., HM12, P19, KS19). Instead of a detailed convolution of the spectral model with a distribution of stellar ages and metallicities, we use fixed templates representative of high-redshift star-forming galaxies, where stars are most important for the UVB. Detailed observations of $z\sim 2-3.5$ LBGs have recently shown that their nebular emission lines are better modeled by stellar populations including binaries (Bordoloi et al. 2014; Steidel et al. 2016; Strom et al. 2017), whose effective temperatures are higher (for an alternative interpretation of observed line ratios based on chemical abundance ratios see, e.g., Masters et al. 2014; Shapley et al. 2015). When binaries are included, massive hot stars can persist significantly past the few Myrs predicted by single-star models. Support for the importance of stellar binaries at high redshift has also recently emerged from galaxy formation simulations, which find that longer-lived hot stars help boost the effective escape fraction to values high enough to explain HI reionization (e.g., Ma et al. 2015, 2016; Rosdahl et al. 2018). We use the default ‘imf135_300’ initial mass fraction (IMF) in BPASS, corresponding to a power-law slope $\alpha_{1}=1.30$ between 0.1 $M_{\odot}$ and 0.5 $M_{\odot}$ and a slope $\alpha_{2}=2.45$ between 0.5 $M_{\odot}$ and 300 $M_{\odot}$ ($dN(<M)/dM\propto M^{\alpha_{i}}$). Steidel et al. (2018) showed that this choice of IMF yields good fits to LBG spectra. We show in Appendix B that our results would not change significantly if we instead assumed a Chabrier (2003) IMF.

Another recent development regarding high-redshift galaxies is the realization that their stellar and gas-phase metallicities can differ substantially. Steidel et al. (2016) showed that the stellar and nebular spectra of $z\sim 2.4$ LBGs can be reconciled if the massive stars have a metallicity $Z_{\star}\approx 0.1$ Z*⊙* but the nebular gas has a higher metallicity, $Z_{\rm neb}\approx 0.5$ Z*⊙. This apparent discrepancy can be explained by a $\sim 5\times$ supersolar O/Fe abundance for the nebular gas, which is expected for enrichment dominated by core-collapse SNe. While O dominates the physics of the nebular gas observed in emission lines, Fe dominates the extreme and far-UV opacity and controls the mass loss rate from stars. Thus, the Fe abundance most strongly affects the properties of the massive stars which produce most of the ionizing photons. For our fiducial stellar template, we therefore assume a low stellar metallicity $Z_{\star}=0.1$ Z⊙* (corresponding to a metal mass fraction $Z=0.002$).

For dust attenuation, we use the Reddy et al. (2016) attenuation curve calibrated to $z\sim 3$ observations:

[TABLE]

where $k(\lambda)$ is given in Reddy et al. (2015) for $0.15\leq\lambda\leq 2.85$ $\mu$m and in Reddy et al. (2016) for $\lambda<0.15$ $\mu$m. The Reddy et al. (2016) attenuation curve provides a better match to observations of high-redshift galaxies than the Calzetti et al. (2000) attenuation curve, and implies a lower attenuation in the far-UV for a given $E(B-V)$. $S_{\nu}^{\rm intr}$ is the intrinsic spectrum of the stellar population before any dust correction is applied and $S_{\nu}^{\rm dust}$ is the same spectrum but attenuated by dust.

For the escape fraction of ionizing photons ($f_{\rm esc}^{\star}$), we use the “holes” model in which $f_{\rm esc}^{\star}$ is determined by the fraction of sight lines from hot stars along which photons can escape unimpeded, while ionizing photons are completely absorbed along other sight lines. We furthermore assume no dust obscuration along holes. A holes model is supported by high-resolution galaxy formation simulations (e.g., Kim et al. 2013; Cen & Kimm 2015; Ma et al. 2015) as well as observations (e.g., Steidel et al. 2018; Fletcher et al. 2019). The net angle-averaged emergent spectrum from galaxies is thus

[TABLE]

Figure 2 shows our stellar spectral template for fiducial $E(B-V)$ and $f_{\rm esc}^{\star}$ values. We fix $E(B-V)=0.129$ in our calculations, a value which Steidel et al. (2018) found to provide a good fit to a composite $z\sim 3$ LBG spectrum. The stellar spectra assume a continuous star formation history with an age of 300 Myr, which is appropriate for LBGs (e.g., Reddy et al. 2012). The exact age of the stellar population is not critical as the stellar UV flux becomes approximately constant after $\sim 100$ Myr of continuous star formation. The normalization of the spectral template in equation (11) is arbitrary; at any given redshift, the spectral shape is used to evaluate the emissivity at any frequency $\nu$ given the prescribed 1500 Å emissivity (eq 9).

The escape fraction is a major uncertainty in models of the ionizing background. In detail, it likely depends on galaxy properties, including redshift and mass (e.g., Gnedin et al. 2008; Xu et al. 2016). For our purpose, we require an “effective” escape fraction weighted by the intrinsic emissivity of the sources.

By combining constraints on the ionization of the IGM at $z\sim 3$ from the Ly$\alpha$ forest and the observed UV galaxy luminosity function, Kuhlen & Faucher-Giguère (2012) found that reconciling the galaxy UV luminosity function with the reionization history inferred from the WMAP-7 optical depth could be achieved by assuming an escape fraction that increases with increasing redshift (other options are to invoke a large population of galaxies below detection limits, or a substantially higher ionizing efficiency at very high redshift). A similar conclusion is supported by more recent data. For example, Robertson et al. (2015) assume a fiducial escape fraction of 20% in $z\gtrsim 6$ reionization models constrained by HST and Planck observations. Furthermore, the most recent UVB synthesis models all assume a galaxy escape fraction that increases with redshift, albeit with different functional forms (Haardt & Madau 2012; Puchwein et al. 2019; Khaire & Srianand 2019).

Motivated by these results, we assume an effective escape fraction that increases continuously with redshift but constrained to stay under 20%:111This expression implies that the modeled galaxy escape fraction can be very low at $z\sim 0$. Although the constraints on the escape fraction remain relatively poor in this regime (direct Lyman continuum detection requires space-based observations), this is consistent with the fact that many previous attempts to detect escaping ionizing radiation at low redshift have primarily yielded upper limits (e.g., Leitherer et al. 1995; Heckman et al. 2001). Significant escaping ionizing radiation has been detected from some low-redshift galaxies (see the compilation of HST/COS detections in Chisholm et al. 2018), but the existing detections tend to be biased toward extreme starburst galaxies with properties similar to LBGs. Heckman et al. (2011) proposed that extreme feedback by a powerful starburst may be required to enable the escape of ionizing radiation from galaxies. High star formation rates and powerful galactic winds are common at high redshift but much rarer in the nearby Universe, which could explain strong redshift evolution in the escape fraction.

[TABLE]

For our assumed emissivity and IGM properties, we find that $f_{\rm esc}^{\star,z=3}=0.01$ provides a good match to the HI photoionization rate inferred from the Ly$\alpha$ forest at $z\gtrsim 3$ (see §2.5). This value is similar to the absolute escape fraction inferred by other studies by comparing the $\approx 1500$ Å (non-ionizing) galaxy emissivity to the total HI photoionization rate at $z\gtrsim 3$. In particular, after accounting for differences in definitions, Faucher-Giguère et al. (2008b), Haardt & Madau (2012), Puchwein et al. (2019), and Khaire & Srianand (2019) all find that the integrated observational constraints can be matched assuming $f_{\rm esc}^{\star,z=3}\approx 0.01-0.02$.

It is interesting, however, that the absolute escape fraction inferred from integral constraints at $z=3$ is substantially lower than indicated by recent “direct” measurements of the escape fraction from $z\approx 3$ star-forming galaxies. Steidel et al. (2018) carried out a comprehensive analysis of a sample of $z\approx 3$ LBGs and found an average absolute escape fraction of $9\pm 1$%. In an HST survey $z\approx 3.1$ galaxies, Fletcher et al. (2019) found an average escape fraction (accounting for non-detections) of $\approx 5$%. Appendix B discusses some effects that may contribute to the difference, but the resolution is not clear. One possibility is that the higher escape fractions inferred from direct observations are not representative of the UV emissivity-weighted escape fraction. This is plausible since, for example, the state-of-the-art Steidel et al. (2018) and Fletcher et al. (2019) studies only detect escaping ionizing photons from relatively luminous galaxies with $L_{\rm UV}\gtrsim 0.25L^{\star}_{\rm UV}$. For a Schechter (1976) luminosity function with a faint-end slope $\alpha\approx-1.7$, such galaxies account for only $\approx 30\%$ of the total UV emissivity.

At face value, this “incompleteness” effect goes in the direction opposite to what is needed to explain the discrepancy, since observations indicate that $f_{\rm esc}^{\star}$ correlates positively with the Ly$\alpha$ equivalent width $W({\rm Ly\alpha})$, while $W({\rm Ly\alpha})$ tends to increase with decreasing $L_{\rm UV}$ (Steidel et al. 2018). This suggests that, if anything, on average fainter galaxies should have higher escape fractions. Nonetheless, accounting for this $W({\rm Ly\alpha})$ dependence, Steidel et al. (2018) estimate that $L_{\rm UV}\gtrsim 0.25L^{\star}_{\rm UV}$ galaxies alone produce a comoving ionizing emissivity $\epsilon_{\nu 912}^{\rm\star,com}\approx 5.9\times 10^{24}$ erg s*-1* Hz*-1* Mpc*-3* at $z\approx 3$, which is almost exactly equal to the total ionizing emissivity needed to match the Ly$\alpha$ forest-inferred $\Gamma_{\rm HI}$ at $z=3$ in our UVB synthesis model (see Fig. 1). There is thus, at face value, no room for $L_{\rm UV}\lesssim 0.25L^{\star}_{\rm UV}$ galaxies to produce a significant amount of escaping ionizing photons. If the ionizing emissivity predicted using the Shen et al. (2020) AGN luminosity function in the next section is accurate at $z=3$, AGN contribute about half of the total ionizing emissivity at this redshift, and a galaxy escape fraction $\approx 9\%$ may overproduce $\Gamma_{\rm HI}(z=3)$ even if restricted to $L_{\rm UV}\gtrsim 0.25L^{\star}_{\rm UV}$ galaxies. We note, however, that Bian & Fan (2020) did not detect any escaping ionizing photons in a recent study of faint ($L_{\rm UV}\approx 0.1L^{\star}_{\rm UV}$), $z\approx 3.1$ Ly$\alpha$ emitters. Bian & Fan (2020) conclude that this result is in tension with an extrapolation of Steidel et al. (2018)’s relation between $f_{\rm esc}^{\star}$ and $W({\rm Ly\alpha})$ to their sample of high-$W({\rm Ly\alpha})$ galaxies. Thus, it is possible that faint galaxies have lower than escape fractions than suggested by trends measured in samples of brighter galaxies.

2.3 Active galactic nuclei

In FG09, we used the bolometric luminosity function model from Hopkins et al. (2007, hereafter HRH07), which was based on a large compilation of obervational data available at the time, to evaluate the rest-frame AGN emissivity at 4400 Å ($B$ band). We then adopted simple power-law approximations for the AGN spectrum from the optical to the extreme UV. In this paper, we update our AGN treatment to account for new luminosity function measurements (especially at high redshift, where the HRH07 luminosity function relied heavily on extrapolations) and with an improved spectral model including both unobscured and obscured sources.

We use the new analysis of the AGN bolometric luminosity function from Shen et al. (2020). This new analysis is an update of the HRH07 bolometric luminosity function which includes a large number of new AGN measurements that have become available since HRH07’s study. We start with the AGN ionizing emissivity derived by Shen et al. by integrating their luminosity function (at the rest wavelength 912 Å). For the comoving emissivity in erg s*-1* Hz*-1* Mpc*-3*, this gives:

[TABLE]

A spectral template for the AGN population is then used to evaluate the emissivity at other wavelengths.

We derive and use an AGN spectral template that accounts for both unobscured and obscured sources. For unobscured (intrinsic) AGN spectra, our template closely follows that used by Shen et al. (2020) in their analysis of the bolometric luminosity function. Briefly:

• •

Optical-UV ($\lambda\geq 912$ Å): The Krawczyk et al. (2013) template for the mean intrinsic quasar spectrum based on $\sim 100,000$ spectra at $0.064<z<5.46$ (the template extends into the infrared, but we do not use that part in this work).

• •

Ionizing UV ($600<\lambda\leq 912$ Å): A power-law $F_{\nu}\propto\nu^{\alpha}$, where $\alpha=-1.7$ based on a stack of $z\sim 2.4$ quasars corrected for intervening Ly$\alpha$ forest and Lyman continuum absorption (Lusso et al. 2015).

• •

X-rays ($>0.5$ keV): An exponentially-truncated power-law template with photon index $\Gamma$ and cut-off energy $E_{\rm c}$: $F_{E}/E\propto E^{-\Gamma}\exp{(-E/E_{\rm c})}$. We set $\Gamma=1.8$ and $E_{\rm c}=300$ keV, which are representative of the best fits from several X-ray studies (e.g., Dadina 2008; Ueda et al. 2014; Aird et al. 2015). To this, we add a Compton reflection component using the PEXRAV model (Magdziarz & Zdziarski 1995), assuming a reflection strength of unity, an inclination angle of $60^{\circ}$, and solar abundances.

• •

Extreme UV to X-rays (600 Å to 0.5 keV): We assume a single power law in the energy ‘gap’ between the UV and X-ray components defined above. The power law connects the template fluxes at 600 Å and 0.5 keV, where the relative amplitude of the X-ray component relative to the ionizing UV component is set to match observational measurements of the X-ray-to-optical ratio, $\alpha_{\rm OX}$, given the optical-to-UV and X-ray spectral shapes. Let $L_{\nu}(2~{}{\rm keV})$ and $L_{\nu}(2500~{}{\rm\AA})$ be specific luminosities in units of erg s*-1* Hz*-1*. Observations have found a correlation of the form $\log{L_{\nu}(2~{}{\rm keV})}=B\log{L_{\nu}(2500~{}{\rm\AA})}+C$, which implies a luminosity-dependent $\alpha_{\rm OX}=$ $\frac{\log{L_{\nu}(2~{}{\rm keV})}-\log{L_{\nu}(2500~{}{\rm\AA})}}{\log{\nu({\rm 2~{}keV})}-\log{\nu(2500~{}{\rm\AA})}}$ $=0.384\log{\left(\frac{L_{\nu}(2~{}{\rm keV})}{L_{\nu}(2500~{}{\rm\AA})}\right)}$. Since $\alpha_{\rm OX}$ is luminosity-dependent, we must choose a representative AGN luminosity to set the relative amplitude of the UV and X-ray components. We use the value of $L_{\nu}(2500~{}{\rm\AA})$ corresponding to a $B$-band optical luminosity $L_{\rm B}=\nu L_{\nu}|_{\rm 4400~{}\AA}=10^{43.3}$ erg s*-1*. This is roughly the optical luminosity of AGN that contribute most to the AGN UV emissivity at $z=0.25$ (see §2.3).

One approach to obtain a spectral template for obscured AGN is to convolve the intrinsic AGN spectrum with a distribution of obscuring columns, $N_{\rm H}$, making assumptions about the dust content of obscuring material. This is, for example, what is done to model obscured AGN in the bolometric luminosity function analyses of HRH07 and Shen et al. (2020). We follow a different approach, similar to what was done by Sazonov et al. (2004, hereafter SOS04) to construct an older AGN spectral template. The idea is to assume that the CXB (at $z=0$) above some minimum energy is dominated by a sum over all AGN (obscured and non-obscured), while the optical-UV background is dominated by non-obscured AGN. For a given AGN luminosity function, we can then derive what X-ray spectral shape is needed to match the observed CXB.

More specifically, we define the following SED for the total AGN population averaged over unobscured and obscured sources:

[TABLE]

where $f_{\rm obsc}$ is the obscured fraction, $F_{E}^{\rm int}$ is the mean intrinsic AGN SED defined above, and $F_{E}^{\rm obsc}$ is a template for the mean SED of obscured AGN. For simplicity, we first neglect Compton thick AGN ($N_{\rm H}>10^{24}$ cm*-2*) and normalize the templates such that $F_{E}^{\rm int}({\rm 50keV})=F_{E}^{\rm obsc}({\rm 50keV})=F_{E}^{\rm tot}({\rm 50keV})$ (i.e., we neglect AGN that are significantly absorbed at 50 keV). We can then solve for the obscured AGN template:

[TABLE]

We use the fact that obscured AGN contribute non-negligibly only to X-rays and higher energies, and that both obscured and non-obscured sources contribute to the CXB. Furthermore, we assume (as observations indicate, e.g., Hickox & Markevitch 2006; Cappelluti et al. 2017) that most of the $\gtrsim 2$ keV CXB is produced by AGN. As in SOS04, we thus have the requirement that the total AGN spectral template $F_{E}^{\rm tot}$ should produce the measured CXB after integrating over the AGN luminosity function and redshifting. We use a functional form based on SOS04 for the total AGN template in the X-ray regime:

[TABLE]

SOS04 found that the following parameters provided a good fit to a variety of observational constraints on AGN spectra, including the CXB: $A=2^{\alpha}e^{2/E_{1}}$, $B=E_{0}^{\beta-\alpha}e^{-(\beta-\alpha)/(1+kE_{0}^{\beta-\alpha})}$, $\alpha=0.24$, $\beta=1.6$, $\gamma=1.06$, $E_{1}=83$ keV, $k=4.1\times 10^{-3}$. For our UVB synthesis model, we adopt the same parameters as SOS04, except for $E_{1}$ which we set to 67 keV. The different value is due to the fact that we assume a different redshift evolution for the AGN luminosity function in this work, so different spectral parameters are necessary to produce a good fit to the CXB (we compare our model results to the observed CXB in §2.5). Figure 3 shows the resulting AGN spectral template.

We find that $f_{\rm obsc}=0.75$ provides an excellent match to the relative intensities of the UV and X-ray backgrounds. Since we normalized the obscured template such that it matches the intrinsic template at 50 keV, $f_{\rm obsc}$ may be roughly interpreted as the fraction of AGN that are obscured in the UV (or optical) but not at 50 keV. Let us call these the Type II AGN. According to this definition, Type II AGN exclude the most heavily obscured, Compton-thick AGN ($N_{\rm H}>10^{24}$ cm*-2*), which are needed to explain the shape of the CXB (e.g., Gilli et al. 2007). Although some CXB models imply that Compton-thick AGN are roughly as abundant as Type II AGN (e.g., Ueda et al. 2014), which would perhaps imply an uncomfortably high total obscured fraction, the inferred abundance of Compton-thick AGN from CXB modeling is degenerate with the assumed strength of the X-ray reflection component from the total AGN population. Recent studies based on the direct detection and modeling of heavily obscured AGN in hard X-rays using the NuSTAR satellite however indicate a much lower fraction of Compton-thick AGN, corresponding to $\sim 10-20\%$ of the Type II population (Masini et al. 2018; Georgantopoulos & Akylas 2019).

We also considered the AGN UV luminosity function recently derived by Kulkarni et al. (2019) from a compilation of different surveys. However, unlike the Shen et al. (2020) model, the Kulkarni et al. (2019) luminosity function has a faint-end slope so steep that the total emissivity is sensitive to the limiting magnitude of the integral (e.g., $M_{\rm lim}=-18$ vs. $M_{\rm lim}=-21$), i.e. the abundance of very faint AGN for which direct observational constraints are poor. We compare the ionizing emissivities implied by the HRH07, Kulkarni et al. (2019), and Shen et al. (2020) AGN luminosity functions and expand on the issue of completeness in Appendix C.

2.4 Recombination emission

For the recombination emissivity, we include Lyman continuum (LyC; 1 Ry) and Ly$\alpha$ (0.75 Ry) recombination emission from HI; and LyC (4 Ry), Balmer continuum (BalC; 1 Ry), and Ly$\alpha$ recombination emission from HeII (3 Ry). For HI LyC and HeII LyC, BalC, and Ly$\alpha$, we use the approximations developed in FG09.

For HI Ly$\alpha$, we assume that 0.68 HI Ly$\alpha$ photon is produced for each ionizing photon absorbed in the ISM of a star-forming galaxy ($\propto(1-f_{\rm esc}^{\star}$)), corresponding to case B recombination. We then assume that a fraction of these Ly$\alpha$ photons are destroyed by dust before escaping the galaxy by applying the same dust attenuation that we use to attenuate the UV continuum of galaxies. At any redshift, the emissivity from this process is modeled as a spatially-homogeneous $\delta$-function in frequency. Note that this is only a crude approximation for HI Ly$\alpha$ as it neglects resonant scattering effects, which can strongly affect the radiative transfer in this line (e.g., Dijkstra et al. 2006; Verhamme et al. 2006; Faucher-Giguère et al. 2010; Smith et al. 2019). Moreover, we neglect Ly$\alpha$ photons produced by quasars and by recombinations from intergalactic clouds.

Relative to FG09, we have improved our radiative transfer code to include “sawtooth” absorption by HeII Lyman series (Lys) lines between 3 and 4 Ry and the corresponding recombination emission following Madau & Haardt (2009).

2.5 Integral constraints and synthesis results

The most stringent constraints on our UVB model are integral constraints that do not require detailed assumptions about the radiation sources. In this section, we compare our new synthesis model with several integral constraints as well as to the previous FG09, HM12, and fiducial KS19 model (their Q18 model). In §3, we will also compare our new model to results from Oñorbe et al. (2017) and Puchwein et al. (2019) on modeling the effects of reionization.

Figure 4 compares the frequency-resolved total UV emissivity at $z\leq 1.5$ inferred by Chiang et al. (2019) to the model emissivities for galaxies and AGN based on the assumptions detailed in §2.2-2.3. Interestingly, the spectral shape of the total emissivity is in excellent agreement with the BPASS stellar template used in our UVB model (including an apparent but subtle change in spectral slope around 1216 Å), providing support for this stellar template.

Figures 5 and 6 compile $\Gamma_{\rm HI}$ (Becker & Bolton 2013; Kollmeier et al. 2014; Gaikwad et al. 2017; Viel et al. 2017; D’Aloisio et al. 2018; Khaire et al. 2019) and $\Gamma_{\rm HeII}$ (Khaire 2017; Worseck et al. 2018) measurements from quasar absorption spectra. These constraints include contributions from all ionizing sources, even galaxies or AGN too faint to be individually detected. We also include an estimate of $\Gamma_{\rm HI}$ at $z=0$ based on modeling H$\alpha$ fluorescence emission from a nearby HI disk (Fumagalli et al. 2017).

Figure 7 shows the spectrally-resolved model background intensity at selected redshifts ($z=0,~{}1,~{}2$, and $3$) and compares to the previous FG09, HM12, and KS19 UVB models. In Figure 8, we compare the models with measurements of the local X-ray background from Cappelluti et al. (2017, using Chandra), Churazov et al. (2007, using INTEGRAL), and Ajello et al. (2008, using Swift BAT).

Overall, our new UVB model provides an excellent simultaneous match to the empirical constraints considered. This is in large part because we calibrated some model parameters to satisfy these constraints, but it is noteworthy that we were able to so by including only radiation from galaxy and AGN populations consistent with observed luminosity functions, and assuming well-motivated spectral templates for each.

One apparent exception is $\Gamma_{\rm HI}$ at $z\sim 4.5-5$, which our UVB model appears to underestimate relative to the inference from the Ly$\alpha$ forest by Becker & Bolton (2013). However, our synthesis model agrees well with the measurements from D’Aloisio et al. (2018), which overlap at $z\sim 5$ with the redshift range covered by Becker & Bolton (2013). Since the discrepancy suggests that either the Becker & Bolton (2013) or the D’Aloisio et al. (2018) data points are affected by systematic effects at $z\sim 5$, we do not require our model to match the discrepant points from Becker & Bolton (2013).

It is also noteworthy that our model provides a good match to most $\Gamma_{\rm HI}$ observational constraints at $z<0.5$. In this regime, the UVB in our model is dominated by AGN whose contribution is calculated using relatively standard assumptions: an AGN luminosity function and spectral template calibrated to the latest data (§2.3), and an unity escape fraction for AGN. Thus, our UVB model does not appear to require non-standard ionizing sources or modeling assumptions to explain the low-redshift HI ionizing background (cf. Kollmeier et al. 2014).

Kollmeier et al. (2014) reached a different conclusion by comparing their $\Gamma_{\rm HI}$ measurement plotted in Figure 6 to the HM12 synthesis model. Our results suggest that the “photon underproduction crisis” noted by Kollmeier et al. (2014) can be resolved by a combination of two factors. First, our new synthesis model predicts a higher total $\Gamma_{\rm HI}$ than HM12 at low redshift, by factor $\approx 2$ at $z=0$. Second, the Kollmeier et al. (2014) measurement is a factor $\approx 3$ outlier in the opposite direction relative to subsequent $\Gamma_{\rm HI}$ measurements around the same redshift (also shown in Fig. 6). Together, these effects account for the factor of $\approx 5$ discrepancy reported by Kollmeier et al. (2014) between the HM12 synthesis model and their $\Gamma_{\rm HI}$ measurement. Shull et al. (2015) previously noted that assuming a higher escape fraction from star-forming galaxies could also reconcile observed galaxies and AGN with the low-redshift Ly$\alpha$ forest. Although this is possible, our analysis suggests that this is not required since AGN can account for most of the low-redshift ionizing background.

The HM12 synthesis model also significantly under-predicts the X-ray background measured by Cappelluti et al. (2017) using Chandra data at $\sim 2-7$ keV, by up to $\approx 30$% (see also Hickox & Markevitch 2006). However, HM12 show that their model are in better agreement with the HEAO-1 measurement from Gruber et al. (1999) in this energy range. Figure 7 shows that the difference in the X-ray spectrum between our synthesis model and HM12 increases in magnitude with increasing redshift. This is thus a regime where the two synthesis models can produce different results. The KS19 synthesis model is in better agreement with our new model at $\gtrsim 1$ keV. The FG09 model did not include X-rays from obscured AGN and therefore under-predicts the X-ray background at all redshifts.

Another regime where current UVB models differ significantly is the energy range between 4 Ry (the HeII photoionization edge) and 1 keV. Figures 5, 6, and 7 in particular show that our new UVB model predicts a specific intensity $J_{\nu}$ at energies near (but above) 4 Ry lower by a factor $\approx 2$ than the HM12 and KS19 models. Differences in this energy range arise because observational constraints on intrinsic AGN spectra and the IGM opacity in this regime are relatively poor, so the results are affected by different assumptions consistent with available data. In our model, the AGN spectral shape in this regime is set by connecting observational measurements in the UV ($\approx 600$ Å) and in the soft X-rays ($\sim 1$ keV), while the HeII opacity is not directly constrained by observations but rather predicted by radiative transfer calculations (§2.1). As can be seen in Figure 3, our assumption that the AGN spectral template is a simple power-law from $600$ Å ($\approx 21$ eV) to 1 keV (distinct from the independently-constrained power-law immediately blueward of 912 Å) implies a softening of the ionizing continuum at $600$ Å. For models calibrated to match the same HI ionization rate, this produces a lower and softer AGN flux near 4 Ry, relative to a single power-law extrapolation of the $F_{\nu}\propto\nu^{-1.7}$ continuum immediately blueward of 912 Å all the way to X-rays. This lower and softer HeII ionizing flux explains the systematic offsets relative to the HM12 and KS19 models, which assumed single power-law spectra from 912 Å to the soft X-rays.

3 Effective reionization history

We now derive “effective” photoionization and photoheating rates designed to produce reionization histories with specified parameters. While reionization is expected to be highly inhomogeneous (e.g., Barkana & Loeb 2001; Furlanetto et al. 2006; Trac & Gnedin 2011), most simulation codes do not follow the time-dependent radiative transfer necessary to model the inhomogeneities. Instead, most codes assume that the IGM is in photoionization equilibrium with a prescribed uniform homogeneous UVB. Our goal here is to provide homogeneous photoioionization rates constructed to produce a prescribed volume-averaged ionized fraction versus redshift under the assumption of photoionization equilibrium. We follow an approach similar to Oñorbe et al. (2017), who developed some ideas introduced in FG09.

3.1 Effective photoionization rates

We begin by describing how we modify photoionization rates. In photoionization equilibrium (neglecting collisional ionization),

[TABLE]

where $\alpha_{\rm HII}$ is the recombination coefficient of HII into HI. The ionic number densities and IGM temperature $T_{\rm IGM0,HI}$ in equation (17) refer to values at the mean density of the Universe. The effects of clumping of intergalactic gas on photoionization balance are modeled through the clumping factor $C_{\rm HI}$. The HI subscripts here do not refer to the HI gas, but rather to the IGM temperature and clumping factor to use during HI reionization (recombinations require ionized gas). We assume a constant IGM temperature in constructing the effective photoionization rates, although in reality the IGM temperature changes during reionization events, e.g. owing to photoheating (§3.2). We use the case A recombination coefficient because this is the limit normally used in cosmological simulation codes.

Given a redshift-dependent hydrogen ionized fraction $x_{\rm II}(z)$, we can solve for the effective $\Gamma_{\rm HI}^{\rm rei}(z)$ that will produce the desired reionization history in photoionization equilibrium:

[TABLE]

where $\chi=Y/4X$. This result assumes that HeI is reionized simultaneously to HI (these two atoms have similar ionization potentials and can be simultaneously ionized by star-forming galaxies), but that HeII reionization is delayed until the rise of the AGN luminosity function. Using an expression analogous to equation (17) for photoionization equilibrium of HeI and assuming the same clumping factor, the assumption that hydrogen and helium are first ionized in concert ($x_{\rm II}=y_{\rm II}$) implies an effective HeI photoionization rate during reionization

[TABLE]

Similarly, photoionization equilibrium for HeII in the two-state approximation that all He is either in HeII or HeIII, while H is fully ionized, implies

[TABLE]

where $T_{\rm IGM0,HeII}$ is the assumed IGM temperature at mean density during HeII reionization. The effective HeII photonization rate is therefore

[TABLE]

We parametrize reionization events using a redshift $z_{\rm rei,i}$ and a redshift width $\Delta z_{\rm rei,i}$. We define a smoothly varying function which tends to one for $z\ll z_{\rm rei,i}$ and to zero for $z\gg z_{\rm rei,i}$. We choose a function with compact support, which has the advantage of unambiguously defining redshifts at which reionization starts and ends:

[TABLE]

for $z\in[z_{\rm rei,i}-\Delta z_{\rm rei,i},~{}z_{\rm rei,i}+\Delta z_{\rm rei,i}]$. Note that this corresponds to one quarter of a sine period, which we define this way so the ionization fraction evolves more slowly at the beginning of reionization and more rapidly toward the end as ionized bubbles overlap. Because of this asymmetry, the redshift of 50% ionization is slightly below $z_{\rm rei,i}$. We use this function to model the evolving ionization fractions during HI and HeII reionization as

[TABLE]

and

[TABLE]

where $x_{\rm II}^{\rm eq}(z)$ and $y_{\rm III}^{\rm eq}(z)$ are the ionization fractions assuming photoionization equilibrium (eqs (17) and (20)), solved for by using the post-reionization HI and HeII photoionization rates calculated in the homogeneous background approximation (to be explicit, we denote these post-reionization, equilibrium rates $\Gamma_{\rm i}^{\rm eq}$ below). This ensures continuity in the final effective photoionization rates:

[TABLE]

For our effective reionization models, we use $\Delta z_{\rm rei,i}$ values that can be smaller than the actual duration of reionization events. This is because the equilibrium ionization fraction for a spatially-homogeneous photoionization rate depends on local density. If we used a large $\Delta z_{\rm rei,i}$, then low-density regions would be reionized significantly before high-density regions. This would introduce large differences in reionization time tied to local density, which is not an accurate model for the propagation of large-scale ionizing fronts during reionization (e.g., Furlanetto et al. 2004; McQuinn et al. 2007). It would also introduce large differences between volume-weighted and mass-weighted ionization histories. For these reasons, we do not attempt to match the reionization history implied by our effective photoionization rates to the detailed redshift evolution of the HI ionized fraction constrained by observations other than the CMB optical depth (e.g., Davies et al. 2018; Greig et al. 2019; Mason et al. 2019). Instead, we force different gas densities to be reionized at roughly the same redshift by using relatively small $\Delta z_{\rm rei,i}$ values.

3.2 Effective photoheating rates

We can also improve the accuracy of the reionization treatment in simulations that assume a homogeneous UVB by defining effective photoheating rates during reionization. This is needed because the homogeneous approximation neglects optical depth effects which significantly affect photoheating during reionization (e.g. Abel & Haehnelt 1999).

Let us first focus on HI reionization. Assuming that photoheating due to reionization is perfectly coupled to the reionization process, we postulate that the rate of change of the IGM temperature during reionization can be approximated as

[TABLE]

where $\Delta T_{\rm HI}$ is the total temperature increment due to reionization heating. Such an approach was suggested by Faucher-Giguère et al. (2009) and implemented in more detail by Oñorbe et al. (2017). The HI photoheating rate is defined such that

[TABLE]

and analogous expressions define the HeI and HeII photoheating rates.

As for the effective photoionization rates in the previous section, we can use this to define an effective total photoheating rate during HI reionization (arising as the sum of photoheating from the photoionization of HI and HeI):

[TABLE]

where $n_{\rm free}$ is the total number of free particles that share the thermal energy. During HI reionization, $\Delta n_{\rm free}/n_{\rm free}\ll\Delta T/T$ since $n_{\rm free}$ changes only by order unity while the IGM temperature increases from $T\sim 10$ K to $T\sim 10^{4}$ K. This implies that $n_{\rm free}(dT/dt)\gg T(dn_{\rm free}/dt)$ and therefore

[TABLE]

where we have used equation (28) in the second step, and the same assumptions as before regarding the simultaneous reionization of HI and HeI in the third step.

Since we assume that HeI is reionized simultaneously to HI, we partition this total heating between the photionization of HI and HeI according to the relative number densities of H and He nuclei:

[TABLE]

Using analogous approximations as for HI reionization but now applied to HeII reionization, which we assume proceed after HI/HeI reionization has completed, we can derive a similar effective photoheating rate:

[TABLE]

3.3 Reionization history parameters

The parameters of the reionization history are the subject of active observational research. It is also useful for theoretical studies to explore the effects of different reionization histories, e.g. on the properties of dwarf galaxies in cosmological simulations (e.g. Benítez-Llambay et al. 2015; Garrison-Kimmel et al., in prep.). Our approach is therefore to provide a UVB model with a fiducial reionization history motivated by current constraints and to also make available, in electronic form, data for other reionization parameters.

An important calibration for our UVB models is that they should produce a reionization history consistent with the electron scattering optical depth to the surface of last scattering measured from CMB observations. This optical depth is given by

[TABLE]

where $\sigma_{\rm T}$ is the Thomson cross section.

For our fiducial reionization model, we assume the following parameters for HI reionization:

[TABLE]

and the following parameters for HeII reionization:

[TABLE]

Figure 9 shows the resulting ionization fractions of HI and HeII, assuming photoionization equilibrium. The electron scattering optical depth for this model is $\tau_{\rm e}=0.054$, equal to the best fit to the Planck 2018 CMB data.

The other reionization parameters were chosen as follows. The redshift of HeII reionization, $z_{\rm rei,HeII}=3.5$, is consistent with simulations of HeII reionization that match the observationally-inferred HeII photoionization rate at $z\sim 2.5-3$ (e.g., McQuinn et al. 2009; Worseck et al. 2018) and is also broadly consistent with inferences from the evolution of the IGM temperature (e.g., Lidz et al. 2010; Becker et al. 2011; Walther et al. 2019). We use approximate clumping factors motivated by a combination of numerical simulations and analytic arguments for the gas in which IGM recombinations take place (e.g., McQuinn et al. 2011; Finlator et al. 2012; Jeeson-Daniel et al. 2014; Kaurov & Gnedin 2014, 2015). The mean-density IGM temperatures to use here should be representative of the redshift when most recombinations take place during reionization, i.e. around the reionization mid-point. We use a rough value of $T_{\rm IGM0,HI}=T_{\rm IGM0,HeII}=10^{4}$ K; since the recombination coefficients scale as $\propto T^{-0.6}$ in the relevant temperature range, the exact choice has only a modest impact on the results. For HI reionization heating, we use a value $\Delta T_{\rm HI}=20,000$ K representative of standard galaxy-driven reionization (e.g. McQuinn 2012) while for HeII reionization heating we use a value $\Delta T_{\rm HeII}=15,000$ K appropriate for reionization by quasars with an extreme-UV spectrum $f_{\nu}\propto\nu^{-1.7}$ (e.g., McQuinn et al. 2019; FG09).

The effective photoionization and photoheating rates for this fiducial reionization model are plotted in Figure 10 and tabulated in Appendix D. The figure also shows effective reionization rates from Oñorbe et al. (2017) and Puchwein et al. (2019).

4 Discussion and conclusions

In this paper, we updated the Faucher-Giguère et al. (2009) cosmic UVB model. Our goal was to incorporate several new observational constraints on the galaxy and AGN populations that putatively dominate the background, an improved model of the IGM opacity, as well as some new modeling elements. Among the new modeling elements, we used the BPASS spectral model for star-forming galaxies including binary stars, a combination of obscured and unobscured AGN (necessary to simultaneously match the empirical UV and X-ray backgrounds), and the HeII Lyman-series sawtooth feature (Madau & Haardt 2009). We also computed effective photoionization and photoheating rates to produce a reionization history consistent with the latest electron scattering optical depth from Planck (Planck Collaboration et al. 2018) and constraints on HeII reionization from quasar absorption spectra in simulations that assume photoionization equilibrium with a homogeneous background.

For our new model, which provides an overall good match to a range of observational constraints, we find the following:

The HI ionizing background, quantified by the photoionization rate $\Gamma_{\rm HI}$, is dominated by AGN at $z\lesssim 3$ and by star-forming galaxies at $z\gtrsim 3$. This is the case even though the non-ionizing background at rest-frame wavelength UV $1500$ Å is dominated by star-forming galaxies at all redshifts. 2. 2.

AGN can explain the entire HI ionizing background at $z<0.5$ inferred by recent modeling of the Ly$\alpha$ forest, assuming a standard spectral energy distribution template and unity escape fraction for AGN. Our UVB model thus suggests that there is no “photon underproduction crisis” (cf., Kollmeier et al. 2014). 3. 3.

The fact that AGN can explain the low-redshift HI ionizing background but fall increasingly short of explaining the Ly$\alpha$ forest transmission at $z>3$ due to the steeply declining AGN luminosity function suggests that the population-averaged escape fraction of ionizing photons from star-forming galaxies increases strongly from $z\sim 0$ to $z\gtrsim 6$. This could be due to powerful stellar feedback clearing large “holes” in the ISM of early galaxies. 4. 4.

Our UVB synthesis model matches the total HI photoionization rate inferred from the Ly$\alpha$ forest for an effective (UV emissivity-weighted) absolute escape fraction of 1% for $z=3$ galaxies. This is substantially lower than the escape fractions $\approx 5-10$% implied by recent direct measurements of escaping Lyman continuum photons (e.g., Steidel et al. 2018; Fletcher et al. 2019). This may indicate that escape fraction trends observed in current direct measurement samples cannot be extrapolated to fainter galaxies. 5. 5.

The non-ionizing part of the low-redshift UVB spectrum inferred from cross-correlating GALEX observations with SDSS spectroscopic objects (Chiang et al. 2019) is in excellent agreement with the UVB spectrum predicted using the BPASS stellar template. 6. 6.

The low-redshift HI ionizing background and the local X-ray background can be simultaneously explained by an AGN population that includes both obscured and unobscured sources and are consistent with a large AGN obscured fraction $f_{\rm obsc}\approx 0.75$. 7. 7.

As in previous models, the HeII ionizing background is dominated by AGN at all redshifts due to the strong spectral break at the HeII photoionization edge in stellar spectra. While star-forming galaxies drive HI reionization, AGN drive HeII reionization at lower redshifts.

One motivation for updating our UVB model was to better understand the uncertainties in available UVB models. We find that different recent synthesis models (e.g., HM12, KS19, P19, and this work) are qualitatively similar in their predicted spectra and redshift evolution. However, important quantitative differences remain. Relative to HM12, our new model provides a better match to recent measurements of the low-redshift Ly$\alpha$ forest (higher $\Gamma_{\rm HI}$) as well as to the local X-ray background determined by Chandra. Furthermore, the differences in X-ray predictions increase with increasing redshift. Different models also differ notably in their predictions in the energy range between the HeII photoionization edge (4 Ry) and $\approx 1$ keV, due to relatively poor observational constraints on the intrinsic AGN spectrum and IGM opacity in this regime. As for the modeling of reionization, to our knowledge this work presents the first full UVB spectral synthesis model to include effective photoionization and photoheating rates that match the best-fit Planck 2018 electron scattering optical depth, corresponding to a relatively late redshift of HI reionization $z_{\rm rei,HI}\sim 7.8$.

We plan to release model updates as improved empirical constraints become available.222See http://galaxies.northwestern.edu/uvb for the electronic data.

Acknowledgments

We acknowledge José Oñorbe, Jonathan Stern, and Luke Zoltan Kelley for useful discussions and comments. We thank Andrew Wetzel, Nick Frontiere and JD Emberson for testing our new model UVB in hydrodynamic simulations. We are grateful to Rychard Bouwens for providing galaxy luminosity function measurements ahead of publication, and to Jacob Shen for providing updated AGN luminosity function and spectral models also before publication. Detailed comments by an anonymous referee substantially improved the paper. This work was supported by NSF through grants AST-1517491, AST-1715216, and CAREER award AST-1652522; by NASA through grants NNX15AB22G and 17-ATP17-0067; by STScI through grants HST-GO-14681.011, HST-GO-14268.022-A, and HST-AR-14293.001-A; and by a Cottrell Scholar Award from the Research Corporation for Science Advancement.

Appendix A Ionizing mean free paths

UVB models depend on assumptions about source populations as well as the transfer of radiation through intervening gas. A key parameter is the mean free path (m.f.p.) of ionizing photons. This is most clearly seen by considering the “local source approximation,” which is valid at redshifts and frequencies such that the m.f.p. in the IGM, $\Delta l$, is much shorter than the Hubble length. In this limit, the UVB intensity at any location can be thought of as arising from the emissivity integrated over all sources within a radius $\Delta l$, i.e.

[TABLE]

(see equation D1 in FG09).

Since photoionization rates scale linearly with $J_{\nu}$ we have, for example,

[TABLE]

in the high-redshift limit in which star-forming galaxies dominate the HI ionizing background. Equation (38) highlights the degeneracy between the m.f.p. model and the escape fraction inferred by comparing the emissivity obtained by integrating over the galaxy luminosity function ($\epsilon_{\nu 1500}^{\star}$) and the total IGM photoionization rate measured using the Ly$\alpha$ forest ($\Gamma_{\rm HI}$).

The m.f.p. corresponds to the length over which the increment of effective optical depth $\bar{\tau}$ is unity. It can therefore be obtained by evaluating the following derivative:

[TABLE]

where $\bar{\tau}$ is given by equation (4). The left hand side of Figure 11 shows the HI ionizing m.f.p. (at $\nu_{\rm HI}$) corresponding to the HI column density distribution adopted in our UVB model (see §2.1). This model m.f.p. is compared to observational measurements from O’Meara et al. (2013), Fumagalli et al. (2013), Prochaska et al. (2009), and Worseck et al. (2014). We also show on this panel the mean distance between HI absorbers with optical depth $\tau=2$ at the Lyman limit from Ribaudo et al. (2011). This is an upper limit to the HI ionizing m.f.p. because it neglects opacity from absorbers with $\tau<2$. Overall, our model is in excellent agreement with observational constraints on the HI ionizing m.f.p.. The dotted curve shows the Hubble length $c/H(z)$. Comparing the m.f.p. with the Hubble length indicates that the IGM photoionization rate is limited by the Hubble volume at $z\lesssim 1$, and therefore that the local source approximation is only valid at $z\gtrsim 1$.

The HeII ionizing m.f.p. (at $\nu_{\rm HeII}$) cannot be computed simply by integrating over the assumed HI column density distribution since it depends on the ratio $\eta=N_{\rm HeII}/N_{\rm HI}$, which itself is a function of the UVB. The HeII ionizing m.f.p. can however been inferred from a full calculation of the UVB spectrum vs. redshift; it is shown on the right hand side of Figure 11 for the fiducial UVB model in this paper.

Appendix B Stellar spectra and escape fraction

Assuming a simple power-law spectrum $J_{\nu}=J_{\nu_{\rm HI}}(\nu/\nu_{\rm HI})^{-\alpha_{\rm HI}}$ blueward of the Lyman limit, the local source approximation implies

[TABLE]

(eq. 25 in Faucher-Giguère et al. 2008b). Consider a high redshift where star-forming galaxies dominate the HI ionizing background,

[TABLE]

where $S_{\nu}$ is the spectral template for stellar radiation escaping galaxies into the IGM (eq. 11). In the limit $f_{\rm esc}^{\star}\ll 1$ and $f_{\rm esc}^{\star}\ll f_{\rm dust}^{1500}$, where $f_{\rm dust}({\rm 1500~{}\AA})\equiv 10^{-0.4E(B-V)k(1500~{}{\rm\AA})}$, the holes model implies the simple scaling

[TABLE]

where we have introduced the notation $(f_{900}/f_{1500})_{\rm out}$ to connect to observational studies which measure the ratio of the flux escaping galaxies at $\approx 900$ Å to the flux escaping at $\approx 1500$ Å. This ratio is more directly constrained by observations than the absolute fraction of ionizing photons that escape galaxies, since the estimating the latter requires a model for the intrinsic ionizing flux of a stellar population (for a discussion of different definitions of the escape fraction, see Steidel et al. 2018). In our numerical evaluations of the fluxes from BPASS spectra, we average $S_{\nu}$ over the wavelength windows 1460-1510 Å and 890-910 Å to avoid local fluctuations due to emission and absorption lines (including an SiII absorption line at 1526 Å). The right-hand side expresses this ratio in terms of the ratio of intrinsic specific fluxes at these rest wavelengths, and the effects of the escape fraction of ionizing photons and dust attenuation. Only the dust attenuation at 1500 Å appears in this expression because escaping ionizing photons are not affected by dust in the holes model (the holes are clear of both ionizing and dust opacity). Combining equations (40)-(42), we find:

[TABLE]

This is the more complete expression corresponding to the scaling in equation (38) above.

Equation (43) shows that, for assumed $\alpha_{\rm HI}$ and $\Delta l(\nu_{\rm HI})$, the combination of a Ly$\alpha$ forest $\Gamma_{\rm HI}$ measurement and the UV luminosity function (which yields $\epsilon_{\nu 1500}^{\rm\star,prop}$) determines what the effective $(f_{900}/f_{1500})_{\rm out}$ must be. In Figure 12, we plot $(f_{900}/f_{1500})_{\rm out}$ vs. $E(B-V)$ for different assumed IMFs, stellar metallicities $Z_{\star}$, and absolute escape fractions $f_{\rm esc}^{\star}$. These curves can be used to gauge how the inferred absolute escape fraction depends on assumed parameters. All curves assume BPASS v2.2.1. We find relatively small differences ($\sim 10\%$) when either the stellar IMF or metallicity is varied. Going from the fiducial $E(B-V)=0.129$ to no dust extinction can reduce $(f_{900}/f_{1500})_{\rm out}$ at given $f_{\rm esc}^{\star}$ by a factor $\sim 2$, i.e. imply a $\sim 2\times$ larger $f_{\rm esc}^{\star}$ for a given $(f_{900}/f_{1500})_{\rm out}$. Unless our understanding of how ionizing photons propagate in the IGM is affected by much larger errors than is presently appreciated (e.g., Prochaska et al. 2014), we do not expect that the $\Delta l(\nu_{\rm HI})$ term or the prefactor involving $\alpha_{\rm HI}$ in equation (43) to be off by more than tens of percent. Thus, overall, it appears difficult to reconcile the effective $f_{\rm esc}^{\star,z=3}\approx 1$% absolute escape fraction implied by our UVB synthesis model (and others; Haardt & Madau 2012; Puchwein et al. 2019; Khaire & Srianand 2019) with recent measurements from the direct detection of escaping Lyman continuum photons implying $\approx 5-10$% escape fractions at $z\approx 3$ (Steidel et al. 2018; Fletcher et al. 2019), other than if faint galaxies not represented in direct studies have much smaller escape fractions (see §2.2).

Appendix C AGN luminosity functions

In Figure 13 we compare the ionizing emissivities vs. redshift implied by different AGN luminosity functions (Hopkins et al. 2007 [used in FG09]; Kulkarni et al. 2019; and Shen et al. 2020 [used in this work]). The emissivities implied by the Hopkins et al. and Shen et al. luminosity functions converge well with limiting magnitude (see Fig. 14), so the curves shown correspond to total emissivities. For the Kulkarni et al. luminosity function, we show the results of integrating town to limiting rest 1450 Å UV magnitudes of -18 and -21, respectively; the curves are significantly different owing to the faint-end slope of the Kulkarni et al. luminosity function. We assume an escape fraction of unity for ionizing photons in each case.

A significant difference between the Hopkins et al. (2007) luminosity function and the more recent Kulkarni et al. (2019) and Shen et al. (2020) determinations is the emissivity implied at $z<0.5$. Namely, the more recent luminosity functions imply emissivities that decline much more rapidly from $z\approx 0.5$ to $z=0$, resulting in a factor $\sim 3$ difference by $z=0$. This difference is important because this is a redshift regime in which our fiducial model implies that AGN dominate the ionizing emissivity over star-forming galaxies, so the AGN ionizing emissivity primarily determines whether we can explain the ionization rate implied by the Ly$\alpha$ at low redshift (see Figure 6). Our results indicate that AGN can indeed provide the majority of the ionizing photons required by the low-redshift Ly$\alpha$ forest, but only if the AGN that contribute significantly the ionizing background have an escape fraction near unity, as we have assumed.

At high redshift ($z\gtrsim 3$), the different luminosity functions also predict integrated ionizing emissivities that can differ by a factor $\sim 3$. The largest difference is between the Shen et al. (2020) luminosity function and the Kulkarni et al. (2019) luminosity function integrated down to a limiting magnitude of -18. On the other hand, the Kulkarni et al. (2019) luminosity function integrated to a limiting magnitude of -21, and the Hopkins et al. (2007) and Shen et al. (2020) luminosity functions integrated over all luminosities predict emissivities within a factor of 2 of each other from $z\sim 0.5$ to $z\approx 6$. Since star-forming galaxies increasingly dominate the HI ionizing background at high redshift, the relatively large differences between AGN luminosity functions have only modest effects on the total HI ionizing background, although the effects can be larger at higher energies.

We now show how the shape of the AGN luminosity function affects which AGN contribute the most to the ionizing emissivity. We consider an arbitrary waveband and drop waveband-specific subscripts to simplify the notation. We define $d\phi/dL$ as the number density of sources with luminosity in the interval $dL$ and $d\phi/dM$ as the number density of sources with absolute magnitude in the interval $dM$. It is common to use a double power-law function to characterize the AGN luminosity function. In terms of magnitudes, we define:

[TABLE]

where $\phi_{\star}$ is the amplitude, $M_{\star}$ is the break magnitude, $\alpha$ is the bright-end slope, and $\beta$ is the faint-end slope. Using the usual conversion between AB magnitude and specific luminosity,

[TABLE]

(Oke & Gunn 1983), the AGN luminosity function can also be expressed as

[TABLE]

where $L_{\star}$ is the break luminosity corresponding to $M_{\star}$.

Since the emissivity $\epsilon=\int dLL(d\phi/dL)=\int d(\ln{L})L^{2}(d\phi/dL)$, each logarithmic bin of luminosity contributes $d(\ln{L})L^{2}(d\phi/dL)$ to the total emissivity. Defining $x=L/L_{\star}$ and expressing in terms of luminosity decades, the double power-law AGN luminosity function implies a fractional contribution to the total emissivity per luminosity decade

[TABLE]

Figure 14 shows how different AGN luminosity bins contribute to the total emissivity. On the left, we show the results for double power-law models with varying faint-end slopes while on the right we show the results for the Hopkins et al. (2007), Kulkarni et al. (2019), and Shen et al. (2020) luminosity functions at $z=0.25$, as a function of UV luminosity at rest-frame 1450 Å. We note that the Hopkins et al. and Shen et al. UV luminosity functions do not follow exact double power-law models, since they are computed from (observationally calibrated) models for the AGN bolometric luminosity function, taking into account luminosity-dependent obscuration and scatter in the bolometric corrections. The figure shows that the Shen et al. luminosity function used in this work converges more rapidly with decreasing limiting luminosity than the Kulkarni et al. luminosity function.

Appendix D Photoionization/heating rates

Table D lists the HI, HeI, and HeIII photoionization and photoheating rates for our UVB model and the fiducial reionization history parameters in §3.3. Before reionization events, these rates are effective homogeneous values appropriate for use in simulation codes that assume photoionization equilibrium with a homogenous UVB.

\ctable

[ caption=Effective Photoionization/heating Rates for the Fiducial Planck 2018 Reionization History,center,star ]lcccccccl \tnote[ ]Reionization history parameters (defined in §3). HI reionization: $z_{\rm rei,HI}=7.8$, $\Delta z_{\rm rei,HI}=0.5$, $C_{\rm HI}=3$, $T_{\rm IGM0,HI}=10,000$ K, $\Delta T_{\rm HI}=20,000$ K. HeII reionization: $z_{\rm rei,HeII}=3.5$, $\Delta z_{\rm rei,HeII}=0.5$, $C_{\rm HeII}=3$, $T_{\rm IGM0,HeII}=10,000$ K, $\Delta T_{\rm HeII}=15,000$ K.

$z$
$\Gamma_{\rm HI}^{\rm eff}$
$\Gamma_{\rm HeI}^{\rm eff}$
$\Gamma_{\rm HeII}^{\rm eff}$
$\dot{q}_{\rm HI}^{\rm eff}$
$\dot{q}_{\rm HeI}^{\rm eff}$
$\dot{q}_{\rm HeII}^{\rm eff}$

10*-12* s*-1*
10*-12* s*-1*
10*-12* s*-1*
10*-24* erg s*-1*
10*-24* erg s*-1*
10*-24* erg s*-1*

$0.0$ $3.62\times 10^{-2}$ $1.19\times 10^{-2}$ $2.84\times 10^{-4}$ $0.203$ $0.127$ $8.77\times 10^{-3}$

$0.2$ $8.72\times 10^{-2}$ $3.50\times 10^{-2}$ $6.71\times 10^{-4}$ $0.516$ $0.359$ $1.93\times 10^{-2}$

$0.4$ $0.177$ $8.04\times 10^{-2}$ $1.40\times 10^{-3}$ $1.075$ $0.837$ $3.81\times 10^{-2}$

$0.6$ $0.312$ $0.152$ $2.58\times 10^{-3}$ $1.926$ $1.635$ $6.72\times 10^{-2}$

$0.8$ $0.488$ $0.246$ $4.20\times 10^{-3}$ $3.034$ $2.732$ $0.105$

$1.0$ $0.682$ $0.352$ $6.00\times 10^{-3}$ $4.260$ $3.995$ $0.146$

$1.2$ $0.862$ $0.452$ $7.64\times 10^{-3}$ $5.414$ $5.218$ $0.182$

$1.4$ $1.003$ $0.533$ $8.74\times 10^{-3}$ $6.338$ $6.229$ $0.203$

$1.6$ $1.091$ $0.590$ $9.12\times 10^{-3}$ $6.939$ $6.950$ $0.207$

$1.8$ $1.136$ $0.622$ $9.10\times 10^{-3}$ $7.264$ $7.386$ $0.198$

$2.0$ $1.136$ $0.631$ $8.58\times 10^{-3}$ $7.316$ $7.518$ $0.180$

$2.2$ $1.109$ $0.627$ $7.47\times 10^{-3}$ $7.191$ $7.432$ $0.153$

$2.4$ $1.069$ $0.611$ $6.23\times 10^{-3}$ $6.960$ $7.196$ $0.125$

$2.6$ $1.021$ $0.589$ $5.02\times 10^{-3}$ $6.678$ $6.874$ $9.76\times 10^{-2}$

$2.8$ $0.968$ $0.564$ $3.91\times 10^{-3}$ $6.359$ $6.508$ $7.33\times 10^{-2}$

$3.0$ $0.915$ $0.537$ $2.14\times 10^{-3}$ $6.027$ $6.112$ $5.19\times 10^{-2}$

$3.2$ $0.857$ $0.509$ $2.30\times 10^{-4}$ $5.671$ $5.716$ $1.65\times 10^{-2}$

$3.4$ $0.805$ $0.482$ $8.23\times 10^{-5}$ $5.349$ $5.340$ $8.34\times 10^{-3}$

$3.6$ $0.757$ $0.457$ $3.16\times 10^{-5}$ $5.035$ $4.995$ $4.91\times 10^{-3}$

$3.8$ $0.711$ $0.431$ $8.08\times 10^{-6}$ $4.744$ $4.679$ $2.42\times 10^{-3}$

$4.0$ $0.670$ $0.410$ $2.81\times 10^{-7}$ $4.465$ $4.393$ $2.80\times 10^{-4}$

$4.2$ $0.631$ $0.388$ $1.00\times 10^{-28}$ $4.203$ $4.136$ $1.00\times 10^{-16}$

$4.4$ $0.592$ $0.365$ $1.00\times 10^{-28}$ $3.946$ $3.899$ $1.00\times 10^{-16}$

$4.6$ $0.558$ $0.348$ $1.00\times 10^{-28}$ $3.718$ $3.706$ $1.00\times 10^{-16}$

$4.8$ $0.525$ $0.331$ $1.00\times 10^{-28}$ $3.502$ $3.526$ $1.00\times 10^{-16}$

$5.0$ $0.493$ $0.315$ $1.00\times 10^{-28}$ $3.300$ $3.369$ $1.00\times 10^{-16}$

$5.2$ $0.463$ $0.301$ $1.00\times 10^{-28}$ $3.116$ $3.233$ $1.00\times 10^{-16}$

$5.4$ $0.436$ $0.287$ $1.00\times 10^{-28}$ $2.947$ $3.101$ $1.00\times 10^{-16}$

$5.6$ $0.409$ $0.273$ $1.00\times 10^{-28}$ $2.779$ $2.971$ $1.00\times 10^{-16}$

$5.8$ $0.384$ $0.259$ $1.00\times 10^{-28}$ $2.615$ $2.835$ $1.00\times 10^{-16}$

$6.0$ $0.356$ $0.245$ $1.00\times 10^{-28}$ $2.441$ $2.704$ $1.00\times 10^{-16}$

$6.2$ $0.329$ $0.230$ $1.00\times 10^{-28}$ $2.268$ $2.564$ $1.00\times 10^{-16}$

$6.4$ $0.304$ $0.216$ $1.00\times 10^{-28}$ $2.102$ $2.418$ $1.00\times 10^{-16}$

$6.6$ $0.279$ $0.202$ $1.00\times 10^{-28}$ $1.943$ $2.282$ $1.00\times 10^{-16}$

$6.8$ $0.257$ $0.189$ $1.00\times 10^{-28}$ $1.801$ $2.160$ $1.00\times 10^{-16}$

$7.0$ $0.237$ $0.177$ $1.00\times 10^{-28}$ $1.670$ $2.040$ $1.00\times 10^{-16}$

$7.2$ $0.167$ $0.125$ $1.00\times 10^{-28}$ $1.188$ $1.446$ $1.00\times 10^{-16}$

$7.4$ $1.25\times 10^{-2}$ $8.03\times 10^{-3}$ $1.00\times 10^{-28}$ $0.107$ $9.42\times 10^{-2}$ $1.00\times 10^{-16}$

$7.6$ $1.44\times 10^{-4}$ $9.96\times 10^{-6}$ $1.00\times 10^{-28}$ $5.77\times 10^{-3}$ $4.81\times 10^{-4}$ $1.00\times 10^{-16}$

$7.8$ $2.80\times 10^{-5}$ $1.93\times 10^{-6}$ $1.00\times 10^{-28}$ $2.43\times 10^{-3}$ $2.03\times 10^{-4}$ $1.00\times 10^{-16}$

$8.0$ $4.26\times 10^{-6}$ $2.94\times 10^{-7}$ $1.00\times 10^{-28}$ $1.11\times 10^{-3}$ $9.27\times 10^{-5}$ $1.00\times 10^{-16}$

$8.2$ $2.81\times 10^{-7}$ $1.94\times 10^{-8}$ $1.00\times 10^{-28}$ $3.56\times 10^{-4}$ $2.97\times 10^{-5}$ $1.00\times 10^{-16}$

$8.4$ $1.24\times 10^{-9}$ $8.62\times 10^{-11}$ $1.00\times 10^{-28}$ $3.13\times 10^{-5}$ $2.61\times 10^{-6}$ $1.00\times 10^{-16}$

$8.6$ $1.00\times 10^{-28}$ $1.00\times 10^{-28}$ $1.00\times 10^{-28}$ $1.00\times 10^{-16}$ $1.00\times 10^{-16}$ $1.00\times 10^{-16}$

$8.8$ $1.00\times 10^{-28}$ $1.00\times 10^{-28}$ $1.00\times 10^{-28}$ $1.00\times 10^{-16}$ $1.00\times 10^{-16}$ $1.00\times 10^{-16}$