THE METAGALACTIC IONIZING BACKGROUND: A CRISIS IN UV PHOTON PRODUCTION OR INCORRECT GALAXY ESCAPE FRACTIONS?

For an isotropic radiation field of specific intensity ${I}_{\nu }$, the normally incident photon flux per frequency is $(\pi {I}_{\nu }/h\nu )$ into an angle-averaged, forward-directed effective solid angle of π steradians. The isotropic photon flux striking an atom or ion is $4\pi ({I}_{\nu }/h\nu )$ and the hydrogen photoionization rate follows by integrating this photon flux times the photoionization cross-section over frequency from threshold (${\nu }_{0}$) to $\infty$

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{\rm{H}}}={\int }_{{\nu }_{0}}^{\infty }\displaystyle \frac{4\pi \;{I}_{\nu }}{h\nu }{\sigma }_{\nu }\;d\nu \approx \displaystyle \frac{4\pi {I}_{0}{\sigma }_{0}}{h(\alpha +3)}\ .\end{eqnarray} \tag{ 1 }$$

Here, we approximate the frequency dependence of the specific intensity and photoionization cross-section by power laws, ${I}_{\nu }={I}_{0}{(\nu /{\nu }_{0})}^{-\alpha }$ and ${\sigma }_{\nu }={\sigma }_{0}{(\nu /{\nu }_{0})}^{-3}$, where ${\sigma }_{0}=6.30\times {10}^{-18}$ cm² and $\alpha \approx 1.4$ for active galactic nuclei (AGNs; Shull et al. 2012c; Stevans et al. 2014). The integrated unidirectional flux of ionizing photons is then

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{0}={\displaystyle \int }_{{\nu }_{0}}^{\infty }\displaystyle \frac{\pi {I}_{\nu }}{h\nu }\ d\nu =\displaystyle \frac{\pi {I}_{0}}{h\alpha },\end{eqnarray} \tag{ 2 }$$

which is related to the density of hydrogen-ionizing photons by ${{\rm{\Phi }}}_{0}={n}_{\gamma }(c/4)$ for an isotropic radiation field. We then have the relations among parameters:

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{\rm{H}}}=4{\sigma }_{0}{{\rm{\Phi }}}_{0}(\displaystyle \frac{\alpha }{\alpha +3})=(8.06\times {10}^{-14}\;{{\rm{s}}}^{-1}){{\rm{\Phi }}}_{4}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{\rm{H}}}=\displaystyle \frac{4\pi {I}_{0}{\sigma }_{0}}{h(\alpha +3)}=(2.71\times {10}^{-14}\;{{\rm{s}}}^{-1}){I}_{-23},\end{eqnarray} \tag{ 4 }$$

where we normalize the incident flux of ionizing photons and specific intensity to characteristic values, ${{\rm{\Phi }}}_{0}=({10}^{4}\;{\mathrm{cm}}^{-2}\;{{\rm{s}}}^{-1}){{\rm{\Phi }}}_{4}$ and ${I}_{0}=({10}^{-23}\ \mathrm{erg}{\mathrm{cm}}^{-2}\;{{\rm{s}}}^{-1}{\mathrm{Hz}}^{-1}{\mathrm{sr}}^{-1}){I}_{-23}$ at the hydrogen Lyman limit ($h{\nu }_{0}=13.60\;\mathrm{eV}$). In an HST/COS survey of IGM metallicity at $z\leqslant 0.4$, Shull et al. (2014) found that ${{\rm{\Phi }}}_{4}\approx 1$ and ${I}_{-23}\approx 3$ gave reasonable fits to the observed ratios of adjacent ionization states of carbon and silicon, (Si iii/Si iv) = ${0.67}_{-0.19}^{+0.35}$, (C iii/C iv) = ${0.70}_{-0.20}^{+0.43}$, and their sum, $({{\rm{\Omega }}}_{{\rm{C}}\;{\rm{III}}}+{{\rm{\Omega }}}_{{\rm{C}}\;{\rm{IV}}})/({{\rm{\Omega }}}_{\mathrm{Si}\;{\rm{III}}}+{{\rm{\Omega }}}_{\mathrm{Si}\;{\rm{IV}}})={4.9}_{-1.1}^{+2.2}$.

Over the past 20 years, numerous papers have estimated the ionizing background and photoionization rate. Table 1 lists ${{\rm{\Gamma }}}_{{\rm{H}}}$ for several models, with values at z = 0 ranging from $(2.28-13.5)\times {10}^{-14}$ s⁻¹. The HM12 rate, ${{\rm{\Gamma }}}_{{\rm{H}}}=2.28\times {10}^{-14}\;{{\rm{s}}}^{-1}$, corresponds to a one-sided ionizing flux ${{\rm{\Phi }}}_{0}=[{{\rm{\Gamma }}}_{{\rm{H}}}(\alpha +3)/4{\sigma }_{0}\;\alpha ]\approx 2630\;{\mathrm{cm}}^{-2}\;{{\rm{s}}}^{-1}$ and specific intensity ${I}_{-23}\approx 0.823$ for the radio-quiet AGN spectral index ($\alpha =1.57$) assumed by HM12. These fluxes are lower by a factor of three compared to the z = 0 metagalactic radiation fields from AGN and galaxies calculated by Shull et al. (1999), ${I}_{\mathrm{AGN}}={1.3}_{-0.5}^{+0.8}\times {10}^{-23}$ and ${I}_{\mathrm{Gal}}={1.1}_{-0.7}^{+1.5}\times {10}^{-23}$, respectively. Adding these two values with propagated errors gives a total intensity and hydrogen ionization rate of ${I}_{\mathrm{tot}}={2.4}_{-0.9}^{+1.7}\times {10}^{-23}\ \mathrm{erg}{\mathrm{cm}}^{-2}\;{{\rm{s}}}^{-1}{\mathrm{Hz}}^{-1}{\mathrm{sr}}^{-1}$ and ${{\rm{\Gamma }}}_{{\rm{H}}}\ ={6.0}_{-2.1}^{+4.2}\times {10}^{-14}\ {{\rm{s}}}^{-1}$. The difference between these radiation fields appears to be the contribution from galaxies. The HM12 background adopts a negligible UVB from galaxies, whereas the background models of Shull et al. (1999), HM01, and HM05 have comparable intensities from galaxies, owing to higher assumed LyC escape fractions.

Previous ultraviolet spectroscopic surveys of low-redshift Lyα absorbers estimated their contribution to the baryon census through the distribution of H i column densities in the diffuse Lyα forest (e.g., Penton et al. 2000, 2004; Lehner et al. 2007; Danforth & Shull 2008). Our recent HST surveys of the low-redshift IGM were more extensive, obtaining 746 Lyα absorbers with STIS (Tilton et al. 2012) and 2577 Lyα absorbers with COS (Danforth et al. 2015). The COS survey used the medium-resolution far-UV gratings (Green et al. 2012) with coverage between 1135–1460 Å (G130M) and 1390–1795 Å (G160M); a few spectra extended slightly outside these boundaries. Our COS survey probed 82 AGN sight lines with cumulative pathlength ${\rm{\Delta }}z=21.7$, with H i column densities ${N}_{{\rm{H}}\;{\rm{I}}}$ (in cm⁻²) between $12.5\lt \mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}\lt 17$. For this paper, we analyze a "uniform redshift-limited sample" of Lyα absorbers at $z\leqslant 0.47$. As discussed in these survey papers, the Lyα statistics are consistent with 24%–30% of the baryons residing in the Lyα forest and partial Lyman limit systems.

Table 2 summarizes the data selected from the STIS survey, with 613 Lyα absorbers at $z\leqslant 0.4$ over the range $12.9\leqslant \mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}\leqslant 14.7$, reprocessed as described in Danforth et al. (2015). Table 3 gives similar results from the COS survey at $z\leqslant 0.47$, with 2074 absorbers between $12.6\;\leqslant \mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}\leqslant 15.2$. For each bin in $\mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}$, we list the number of absorbers (${{\mathcal{N}}}_{\mathrm{abs}}$) and effective redshift pathlength (${\rm{\Delta }}{z}_{\mathrm{eff}}$) over which our Lyα survey is sensitive. From these data, we derive the bivariate distribution, ${d}^{2}{\mathcal{N}}/d(\mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}){dz}$, of absorbers in H i column density and redshift by dividing ${{\mathcal{N}}}_{\mathrm{abs}}$ by ${\rm{\Delta }}{z}_{\mathrm{eff}}$ and ${\rm{\Delta }}\mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}=0.2$. An accurate calculation of this distribution function requires both good absorber counting statistics and knowledge of the pathlength ${\rm{\Delta }}{z}_{\mathrm{eff}}$ covered in the survey for a given column density. As in Danforth & Shull (2008), we compute ${\rm{\Delta }}{z}_{\mathrm{eff}}$ corresponding to the $4\sigma$ minimum Lyα equivalent width as a function of wavelength in each spectrum. Asymmetric error bars arise from statistical uncertainties in ${N}_{\mathrm{abs}}$ and ${\rm{\Delta }}{z}_{\mathrm{eff}}$, computed for each bin using the formalism of Gehrels (1986). At low column densities, the errors are dominated by ${\rm{\Delta }}{z}_{\mathrm{eff}}$, while at high column densities the errors are dominated by the small values of ${N}_{\mathrm{abs}}$ in bins at $\mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}\geqslant 14.6$.

Table 2.  Column Density Distribution^a (STIS Survey)

$\langle \mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}\rangle$	Range in $\mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}$	${{\mathcal{N}}}_{\mathrm{abs}}$	${\rm{\Delta }}{z}_{\mathrm{eff}}$	$f({N}_{{\rm{H}}\;{\rm{I}}},z)$
13.0	(12.9–13.1)	115	4.637	${124}_{-13}^{+13}$
13.2	(13.1–13.3)	115	5.003	${115}_{-11}^{+11}$
13.4	(13.3–13.5)	106	5.216	${102}_{-10}^{+10}$
13.6	(13.5–13.7)	75	5.327	${70}_{-8}^{+9}$
13.8	(13.7–13.9)	73	5.341	${68}_{-8}^{+9}$
14.0	(13.9–14.1)	50	5.360	${47}_{-7}^{+8}$
14.2	(14.1–14.3)	35	5.379	${33}_{-5}^{+6}$
14.4	(14.3–14.5)	26	5.382	${24}_{-5}^{+6}$
14.6	(14.5–14.7)	18	5.382	${17}_{-4}^{+5}$

Note.

^aFor 613 low-redshift Lyα absorbers taken from the HST/STIS survey (Tilton et al. 2012), the columns show: (1) the mean column density (${N}_{{\rm{H}}\;{\rm{I}}}$ in ${\mathrm{cm}}^{-2})$ with bin width ${\rm{\Delta }}\mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}=0.2$; (2) the bin range in $\mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}$; (3) the number of Lyα absorbers (${{\mathcal{N}}}_{\mathrm{abs}}$) in bin; (4) the total redshift pathlength ${\rm{\Delta }}{z}_{\mathrm{eff}}$ at each ${N}_{{\rm{H}}\;{\rm{I}}}$; and (5) the bivariate distribution of absorbers in column density and redshift, $f({N}_{{\rm{H}}\;{\rm{I}}},z)\equiv {d}^{2}{\mathcal{N}}/d(\mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}){dz}$, computed as $({{\mathcal{N}}}_{\mathrm{abs}}/0.2\;{\rm{\Delta }}{z}_{\mathrm{eff}})$ and plotted in Figures 2–7.

Download table as:  ASCIITypeset image

Table 3.  Column Density Distribution^a (COS Survey)

$\langle \mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}\rangle$	Range in $\mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}$	${{\mathcal{N}}}_{\mathrm{abs}}$	${\rm{\Delta }}{z}_{\mathrm{eff}}$	$f({N}_{{\rm{H}}\;{\rm{I}}},z)$
12.7	(12.6–12.8)	194	3.67	${260}_{-120}^{+1000}$
12.9	(12.8–13.0)	292	9.79	${150}_{-40}^{+80}$
13.1	(13.0–13.2)	383	16.46	120 ± 10
13.3	(13.2–13.4)	368	18.78	98 ± 5
13.5	(13.4–13.6)	267	19.23	${69}_{-12}^{+17}$
13.7	(13.6–13.8)	146	12.03	61 ± 5
13.9	(13.8–14.0)	123	12.25	50 ± 5
14.1	(14.0–14.2)	96	12.96	37 ± 4
14.3	(14.2–14.4)	75	13.75	${27}_{-3}^{+4}$
14.5	(14.4–14.6)	47	14.17	${17}_{-2}^{+3}$
14.7	(14.6–14.8)	45	14.24	${16}_{-2}^{+3}$
14.9	(14.8–15.0)	19	14.25	${6.7}_{-1.5}^{+1.9}$
15.1	(15.0–15.2)	19	14.25	${6.7}_{-1.5}^{+1.9}$

Note.

^aStatistics of low-redshift Lyα absorbers ($12.6\leqslant \mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}\lt 15.2$) taken from the HST/COS survey (Danforth et al. 2015, updated 2015 July) using the "uniform sub-sample" of 2074 H i absorbers. The columns show: (1) the mean column density (${N}_{{\rm{H}}\;{\rm{I}}}$ in cm⁻²) with bin width ${\rm{\Delta }}\mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}=0.2$; (2) the bin range in $\mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}$; (3) number of Lyα absorbers (${{\mathcal{N}}}_{\mathrm{abs}}$) in bin; (4) total redshift pathlength ${\rm{\Delta }}{z}_{\mathrm{eff}}$ at each ${N}_{{\rm{H}}\;{\rm{I}}}$; and (5) the bivariate distribution of absorbers in column density and redshift, $f({N}_{{\rm{H}}\;{\rm{I}}},z)\equiv {d}^{2}{\mathcal{N}}/d(\mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}){dz}$, computed as $({{\mathcal{N}}}_{\mathrm{abs}}/0.2\;{\rm{\Delta }}{z}_{\mathrm{eff}})$ and plotted in Figures 2–7.

Download table as:  ASCIITypeset image

An important effect of the IGM is the flux decrement, D_A, produced by Lyα forest line blanketing of the continuum flux. Previous measurements of D_A were reported by Kirkman et al. (2007) using low-resolution UV observations with the HST Faint Object Spectrograph (FOS) toward 74 AGNs. Their measured decrements were fitted to a power law, ${D}_{A}(z)=(0.016){(1+z)}^{1.01}$, over the range $0\lt z\lt 1.6$ with considerable scatter. Based on our moderate-resolution HST/COS survey, Figure 1 shows the fraction of light removed by Lyα absorption from the continuum for $0\lt z\lt 0.47$. These statistics were found by summing the observed-frame equivalent widths of all identified ($\gt 4\sigma$) Lyα absorbers in the Danforth et al. (2015) catalog within a given redshift range, normalized by the effective redshift pathlength probed by the survey for absorbers of that strength. The uncertainty, ${\sigma }_{{D}_{A}}$, is the quadratic sum of measured equivalent width uncertainties normalized in the same manner. The errors are small at lower redshifts ($z\lesssim 0.35$), where most of the surveyed Lyα absorbers are found. However, they remain small compared to D_A even at higher redshifts where ${\rm{\Delta }}{z}_{\mathrm{eff}}$ is smaller. Figure 1 shows our results for bins of width ${\rm{\Delta }}z=0.01$ and ${\rm{\Delta }}z=0.05$ to illustrate the variance on small scales. The scatter in D_A is significantly larger than ${\sigma }_{{D}_{A}}$ and cosmic variance is significant, particularly at $z\gt 0.35$. These fluctuations are large in the ${\rm{\Delta }}z=0.01$ bins, which are comparable to the mean absorber separations for the observed Lyα absorption-line frequencies $d{\mathcal{N}}/{dz}\approx 100\mbox{--}200$ for $\mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}=12.7\mbox{--}13.3$.

Zoom In Zoom Out Reset image size

Figure 1 shows a small dip in D_A between $0.2\lt z\lt 0.3$. This redshift band corresponds to Lyα absorbers at 1459–1580 Å, observed primarily by the COS/G160M grating (1390–1795 Å) since G130M ends at 1460 Å. The dip appears in several ranges of column density, $13\lt \mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}\lt 14$ and $14\lt \mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}\lt 15$. The first range is expected to dominate line blanketing in the Lyα forest for line profiles with Doppler velocity parameters $b\approx 25\mbox{--}30$ km s⁻¹. The dip also remains when we offset the redshift bins. Thus, we believe the dip to be real, although we do not have a plausible physical explanation for its presence.

We fitted the mean values of the HST observations with ${\rm{\Delta }}z=0.05$ bins to the power-law form ${D}_{A}(z)=(0.014\pm 0.001){(1+z)}^{2.2\pm 0.2}$. For $z\leqslant 0.2$, these observations are in reasonable agreement with the lower-resolution HST/FOS results of Kirkman et al. (2007), although we find steeper redshift evolution. By post-processing the H i absorption-line profiles in our simulations (Sections 2.3 and 2.4) over the range $0\lt z\lt 0.2$, we find mean decrements at $\langle z\rangle =0.1$ of ${D}_{A}=0.0116$, 0.0097, and 0.0332 for the HM01, HM05, and HM12 radiation fields, respectively. At $\langle z\rangle =0.1$, these UVBs correspond to hydrogen photoionization rates ${{\rm{\Gamma }}}_{{\rm{H}}}=3.54\times {10}^{-14}\;{{\rm{s}}}^{-1}$ (HM12), $1.38\times {10}^{-13}\;{{\rm{s}}}^{-1}$ (HM01), and $1.86\times {10}^{-13}\;{{\rm{s}}}^{-1}$ (HM05). Our simulations find that ${D}_{A}\propto {{\rm{\Gamma }}}_{{\rm{H}}}^{-0.7}$. In their SPH simulations, K14 found decrements at $z\approx 0.1$ of ${D}_{A}=0.024$ (HM01) and ${D}_{A}=0.050$ (HM12), higher than our simulations by factors of 1.9 and 2.9, respectively. Similar offsets are seen in the distribution of H i column densities.

Our observations at $z\leqslant 0.47$ with STIS (Tilton et al. 2012) and COS (Danforth et al. 2015) yield consistent values for $f({N}_{{\rm{H}}\;{\rm{I}}},z)$, the bivariate distribution of Lyα absorbers in column density and redshift. A least-squares fit of the COS data over the range of $12.7\leqslant \mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}\leqslant 15.2$ had the form, $f({N}_{{\rm{H}}\;{\rm{I}}},z)\equiv {d}^{2}{\mathcal{N}}/d(\mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}){dz}\approx (167){N}_{13}^{-0.65\pm 0.02}$, where we define the dimensionless column density ${N}_{13}=[{N}_{{\rm{H}}\;{\rm{I}}}/{10}^{13}\;{\mathrm{cm}}^{-2}]$. We can use this line frequency to make an analytic estimate of the Lyα line blanketing and flux decrement,

$$\begin{eqnarray}{D}_{A} & = & {\displaystyle \int }_{{N}_{1}}^{{N}_{2}}(\displaystyle \frac{{W}_{\lambda }}{\lambda })f(N,z)\;d(\mathrm{log}N)\\ & = & (\displaystyle \frac{\pi {e}^{2}f\lambda }{{m}_{{\rm{e}}}{c}^{2}})(\displaystyle \frac{167\times {10}^{13}\;{\mathrm{cm}}^{-2}}{2.303}){\displaystyle \int }_{{N}_{1}}^{{N}_{2}}{N}_{13}^{-0.65}\;{{dN}}_{13}.\end{eqnarray} \tag{ 5 }$$

In the last expression of Equation (5), we assume unsaturated Lyα lines, with equivalent widths ${W}_{\lambda }/\lambda =(\pi {e}^{2}/{m}_{{\rm{e}}}c)({Nf}\lambda /c)$, where ${W}_{\lambda }$ and $\lambda =1215.67$ Å are defined in the rest-frame and $f=0.4164$ is the Lyα oscillator strength. We adopt the limits $\mathrm{log}{N}_{1}=12.5$ and $\mathrm{log}{N}_{2}=14.0$ for the integral, which we denote as $I({N}_{1},{N}_{2})\approx 4.49$. The estimated decrement is ${D}_{A}=(0.00323)I({N}_{1},{N}_{2})\approx 0.0145$, close to our observations at $z\approx 0$. Because the line frequency is $f({N}_{{\rm{H}}\;{\rm{I}}},z)\propto {(1+z)}^{2.2}$ out to $z\approx 0.5$ (Danforth et al. 2015), we expect D_A to have the same dependence, as seen in Figure 1. This calculation shows that for the observed steep distribution of H i column densities, most of the line blanketing occurs from moderate-strength Lyα absorbers at ${N}_{{\rm{H}}\;{\rm{I}}}\approx {N}_{2}$, beginning to saturate and appear on the flat portion of the curve of growth. Owing to line saturation, the stronger but rarer Lyα absorbers add a small amount to this estimate.

Our simulations of the H i absorber distributions were run using the Eulerian N-body + hydrodynamics code Enzo (Bryan et al. 2014). The N-body dynamics of the dark matter particles in our simulations were calculated using a particle-mesh solver (Hockney & Eastwood 1988), and the hydrodynamic equations were solved using a direct-Eulerian piecewise parabolic method (Colella & Woodward 1984; Bryan et al. 1995) with a Harten–Lax–van Leer-Contact Riemann solver (Toro et al. 1994). Our modifications of this code and its applications to the IGM are discussed in Smith et al. (2011), where we conducted tests of convergence, examined various feedback schemes, and compared the results to IGM thermal phases, O vi absorbers, and star formation histories using box sizes of $25{h}^{-1}\;\mathrm{Mpc}$ and $50{h}^{-1}\;\mathrm{Mpc}$ and grids of 256³, 384³, 512³, 768³, and 1024³ cells. We note that the Smith et al. (2011) simulations used the HM96 background, rather than HM01 as stated in that paper.³ This code has a substantial IGM heritage, with applications to O vi absorbers (Smith et al. 2011), the baryon census at $z\lt 0.4$ (Shull et al. 2012b), clumping factors and critical star formation rates during reionization (Shull et al. 2012a), and synthetic absorption-line spectra for a comparison of simulations with HST observations (Egan et al. 2014).

Our current simulations were initialized at a redshift z = 99 and run to z = 0, using the WMAP-9 maximum likelihood concordance values (Hinshaw et al. 2013) with ${{\rm{\Omega }}}_{m}=0.282$, ${{\rm{\Omega }}}_{{\rm{\Lambda }}}=0.718$, ${{\rm{\Omega }}}_{{\rm{b}}}=0.046$, ${H}_{0}=69.7$ km s⁻¹ Mpc⁻¹, ${\sigma }_{8}=0.817$, and ${n}_{s}=0.965$ to create the initial conditions. The dark matter density power spectrum used the Eisenstein & Hu (1999) transfer function. Three independent realizations of the initial conditions were created in order to constrain the effects of cosmic variance due to the finite box size. Two of the realizations were created using Music (Hahn & Abel 2011), with second-order Lagrangian perturbation theory. The third realization was created using the initial condition generator packaged with Enzo. The basic simulations discussed in this paper were run on static, uniform grids with a box size of $50{h}^{-1}$ comoving Mpc and 512³ and 768³ cells. We also analyzed a previous simulation with 1536³ cells, produced by Britton Smith on the XSEDE supercomputer and discussed in our study of IGM clumping factors (Shull et al. 2012a), as well as the (${1024}^{3},50{h}^{-1}\;\mathrm{Mpc}$) simulations of Smith et al. (2011). Because those 1024³ models were run with an older UVB from Haardt & Madau (1996), we chose not to compare them to our 512³, 768³ runs, which used more recent backgrounds (HM01, HM05, HM12).

In view of the different predictions of our code for the distribution of low-redshift Lyα absorbers, it is appropriate to compare the resolution⁴ of our grid-code simulations to the SPH simulations run by K14. With a comoving box size of $50{h}^{-1}\;\mathrm{Mpc}$, our models have spatial resolutions (cell sizes) of $97.6{h}^{-1}$, $65.1{h}^{-1}$, and $32.6{h}^{-1}$ kpc for grids of 512³, 768³, and 1536³, respectively. The corresponding dark matter mass resolutions are ${m}_{\mathrm{dm}}=6.1\times {10}^{7}\;{h}^{-1}$ ${M}_{\odot }$ (512³), $1.8\times {10}^{7}\;{h}^{-1}$ ${M}_{\odot }$ (768³), and $2.2\times {10}^{6}\;{h}^{-1}$ ${M}_{\odot }$ (1536³). For comparison, the simulations analyzed by K14 used the SPH code of Davé et al. (2010), also in a $50{h}^{-1}\;\mathrm{Mpc}$ box with 576³ dark matter (and baryon) particles. The characteristic spatial resolution in their SPH method is the interparticle distance, $50{h}^{-1}\;\mathrm{Mpc}/576\approx 86.8{h}^{-1}$ kpc. Under gravitational evolution, the SPH particles are concentrated in regions with baryon overdensity ${{\rm{\Delta }}}_{{\rm{b}}}={\rho }_{{\rm{b}}}/{\bar{\rho }}_{{\rm{b}}}\gt 1$, and the interparticle distance is reduced by a factor of ${{\rm{\Delta }}}_{{\rm{b}}}^{1/3}$. Although K14 do not quote a mass resolution, their simulation ($50{h}^{-1}\;\mathrm{Mpc}$, 576³) would have ${m}_{\mathrm{dm}}\approx 5\times {10}^{7}{h}^{-1}$ ${M}_{\odot }$, scaling from values quoted in the (512³, $100{h}^{-1}\;\mathrm{Mpc}$) SPH calculation of Tepper-Garcia et al. (2012). Therefore, the mass resolution in our 512³ runs is comparable to that of K14, and our 768³ and 1536³ resolutions are superior. Resolution can be important when applied to the large (100–200 kpc) Lyα absorption systems observed in the low-redshift IGM.

Star particles were formed in the simulations using the sub-grid physics prescription of Cen & Ostriker (1992). Star formation occurs in cells where the baryon density is greater than 100 times the critical density, the divergence of the velocity is negative, and the cooling time is less than the dynamical time. As described by Smith et al. (2011), we adopted a star formation efficiency of 10% in these cells. Feedback from star formation returns gas, matter, and energy from star particles to the ISM and IGM. Although star particles are formed instantaneously, feedback occurs gradually with an exponential decay over time. Star particles produce 90% of their feedback within four dynamical times of their creation. Of the star particles' total mass, 25% is returned as gas, with 10% of that taking the form of metals and ${10}^{-5}$ of the rest-mass energy returned as thermal feedback. The feedback is distributed evenly over 27 grid cells centered on the star particle. This scheme avoids the overcooling and subsequent runaway star formation that occurs when feedback is returned to a single cell. We analyzed this "distributed feedback" method in our previous study (Smith et al. 2011) and found that it produces better convergence and removes the overcooling. These simulations were able to reproduce both the observed star formation history and the number density of O vi absorbers per unit redshift over the range $0\lt z\lt 0.4$. Many SPH studies of the IGM (Davé et al. 2010; Tepper-Garcia et al. 2012) suppress the overcooling by suspending the cooling or turning off the hydrodynamics for a period of time, allowing the injected energy and metals to diffuse into the surrounding gas.

The ionization states of H and He in the IGM were calculated with the non-equilibrium chemistry module in Enzo (Abel et al. 1997; Anninos et al. 1997). Metal cooling of the gas was computed using precompiled Cloudy ⁵ (Ferland et al. 2013) tables that assume ionization equilibrium for the metals and are coupled to the chemistry solver (Smith et al. 2008). Photoionization and radiative heating of the gas come from a spatially uniform metagalactic UVB. For each realization of the initial conditions, we ran three simulations that were identical except for the form of the UVB. Two simulations used analytic fits to the HM01 and HM05 backgrounds, initialized at $z=8.9$ and reaching full strength by z = 8. A third background used the HM12 table implemented in the Grackle ⁶ chemistry and cooling library (Bryan et al. 2014; Kim et al. 2014) that begins at $z=15.13$.

In order to analyze the distribution of Lyα absorbers in the simulations, we created synthetic quasar sight lines using the YT⁷ analysis package (Turk et al. 2011). For each simulation box, we created 200 synthetic sight lines, each constructed from simulation outputs spanning the redshift range $0.0\leqslant z\leqslant 0.4$. A randomly oriented ray was chosen through each output, representing the portion of the sightline beginning at the output redshift and with sufficient pathlength to cover the ${\rm{\Delta }}z$ to the next output file. These rays were then combined to create a sightline spanning the entire redshift range. We used identical random seeds for simulations sharing the same initial conditions; differences in the absorber distributions therefore directly reflect the effects of different UVBs. Each sightline is composed of many line elements, from the portions of the sightline passing through a single grid cell. We identified Lyα absorbers as sets of contiguous line elements with neutral hydrogen density ${n}_{{\rm{H}}\;{\rm{I}}}\geqslant {10}^{-12.5}\;{\mathrm{cm}}^{-3}$. Egan et al. (2014) compared this method of identifying absorbers to a more sophisticated method involving synthetic spectra. In their Figure 11, they showed that the resulting absorber column density distributions are in good agreement. We also analyzed a large (${1536}^{3},50{h}^{-1}\;\mathrm{Mpc}$) simulation run by Britton Smith and discussed in Shull et al. (2012a). The light rays from this simulation were generated by Devin Silvia, and we analyzed them in the same manner.

The flux decrement D_A is an imperfect metric for estimating the UVB, since measurements of the observed flux decrement depend on spectral resolution. In addition, both measurements and simulations exhibit significant variance in D_A. A more robust IGM diagnostic comes from the distribution of H i column densities. Figure 2 illustrates the observed STIS and COS distributions, $f({N}_{{\rm{H}}\;{\rm{I}}},z)$, over the range $12.5\leqslant \mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}\leqslant 14.5$. The STIS and COS bins are offset by ${\rm{\Delta }}\mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}=0.1$, using data from Tables 2 and 3. The COS survey is more extensive in both absorber numbers and column density, and its statistical accuracy is superior to that of the STIS survey. Over the overlapping range in $\mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}$, the agreement is good. We overplot values of $f({N}_{{\rm{H}}\;{\rm{I}}},z)$ from our new IGM simulations ($50{h}^{-1}\;\mathrm{Mpc}$ box and 768³ grid cells) using three different ionizing backgrounds (HM01, HM05, HM12). As analyzed further in Section 2.5, the inferred baryon density requires ionization corrections for the neutral fraction, ${n}_{{\rm{H}}\;{\rm{I}}}/{n}_{{\rm{H}}}$. The Appendix suggests that the amplitude of the Lyα absorber distribution, $f({N}_{{\rm{H}}\;{\rm{I}}},z)$, is proportional to ${{\rm{\Gamma }}}_{{\rm{H}}}^{-1/2}{T}_{4}^{-0.363}$ for fixed baryon content (${{\rm{\Omega }}}_{{\rm{b}}}$) and temperature $T=({10}^{4}\;{\rm{K}}){T}_{4}$, and it therefore provides a constraint on the UVB. In our simulations, we find that $f({N}_{{\rm{H}}\;{\rm{I}}},z)\propto {{\rm{\Gamma }}}_{{\rm{H}}}^{-0.773}$ for absorbers in the range $13\lt \mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}\lt 14$.

Zoom In Zoom Out Reset image size

Figures 2–7 illustrate the differential absorber distribution, ${d}^{2}\;{\mathcal{N}}/d(\mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}){dz}$, with tests of convergence, cosmic variance, box size, feedback, and redshift evolution. The simulated distributions are compared to observed distributions from HST surveys of Tilton et al. (2012) and Danforth et al. (2015). Figure 2 presents the primary results of this paper: the H i column-density distributions for three radiation fields (HM01, HM05, HM12) computed with 768³ grids and a $50{h}^{-1}\;\mathrm{Mpc}$ box. The lower HM12 background produces significantly more absorbers than earlier UVBs, ranging from a factor of two for absorbers between $12.5\leqslant \mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}\leqslant 12.7$ to a factor of seven for $14.3\leqslant \mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}\leqslant 14.5$. The observed distributions fall between the results for the HM01 and HM12 backgrounds, approaching the HM12 simulations for $\mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}\gt 14.0$. This suggests that there may be less tension between the HM12 background and observations than implied by the K14 results. Figure 3 explores the effects of cosmic variance, comparing three models with $50{h}^{-1}\;\mathrm{Mpc}$ boxes and 512³ grids. For each radiation field (HM01, HM05, HM12) the simulations labeled with subscripts 1 and 2 are those with the Music initial conditions, and those labeled 3 use the Enzo-packaged initial conditions. For a fixed UVB the scatter among simulations is at the $15\%$ level or less, indicating that variance is unlikely to play an important role. Figure 4 explores the convergence of our simulations, comparing models with $50{h}^{-1}\;\mathrm{Mpc}$ boxes on grids of 512³, 768³, and 1536³. The differences between these models are relatively small, at the 10% level and within the expected differences arising from sample variance. Figure 5 investigates the potential effects of box size, which could suppress structure formation down to low redshift if the box was too small. We ran two simulations ($100{h}^{-1}\;\mathrm{Mpc}$, 512³ and $50{h}^{-1}\;\mathrm{Mpc}$, 256³), both using the same HM01 ionizing background. The grid cell sizes are identical ($195{h}^{-1}$ kpc) and the results are essentially the same for $\mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}\lt 14$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 6 explores the potential influence of our feedback methods. In our 512³ simulations, feedback has little effect on the distribution of H i column densities. We find little difference if we inject mass, metals, and energy within a single cell, in adjacent cells, or eliminating feedback entirely. Clumpiness of the IGM on small scales can be affected by the feedback method, even though the large-scale structural parameters, ($T,{\rho }_{{\rm{b}}}$) and (${{\rm{\Delta }}}_{{\rm{b}}},{N}_{{\rm{H}}\;{\rm{I}}}$), remain nearly the same. Figure 7 investigates the possible effects of redshift evolution of the Lyα forest. We split the simulated results into two redshift intervals ($0\lt z\lt 0.2$ and $0.2\lt z\lt 0.4$) in our 768³ simulations run with the HM05 and HM12 radiation fields. Differences between the simulated distributions are quite small and consistent with COS observations (Danforth et al. 2015) in which the weak Lyα absorbers ($\mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}\lt 14$) exhibit little redshift evolution between z = 0 and $z=0.4$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To compare our grid-code simulations with previous SPH models, we examined two structural measures of the IGM. The first measure is the thermal-phase diagram of $(T,{\rho }_{{\rm{b}}}),$ the temperature versus baryon density in the low-redshift Lyα forest. For our grid models, these diagrams were shown as Figure 19 in Smith et al. (2011) and as Figure 1 in Shull et al. (2012b). For the SPH models, we examined Figure 8 of Davé et al. (2010) and Figure 7 of Tepper-Garcia et al. (2012). All of these diagrams are consistent with a relation $T=(5000\;{\rm{K}}){{\rm{\Delta }}}_{{\rm{b}}}^{0.6}$ where ${{\rm{\Delta }}}_{{\rm{b}}}$ is the baryon overdensity. A second measure is the correlation of baryon overdensity and column density $({{\rm{\Delta }}}_{{\rm{b}}},{N}_{{\rm{H}}\;{\rm{I}}})$, typically expressed as a power law, ${{\rm{\Delta }}}_{{\rm{b}}}={{\rm{\Delta }}}_{0}{N}_{14}^{\alpha }$, where ${{\rm{\Delta }}}_{0}$ is the normalization at a fiducial H i column density ${N}_{{\rm{H}}\;{\rm{I}}}=({10}^{14}\;{\mathrm{cm}}^{-2}){N}_{14}$. These correlations vary with redshift, but typically are evaluated at $z=0.25$. Davé et al. (2010) used $({384}^{3},48{h}^{-1}\;\mathrm{Mpc})$ simulations with the HM01 background to derive a fit for absorbers with $T\lt {10}^{4.5}$ K,

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{{\rm{b}}}=(35.5\pm 0.03){N}_{14}^{0.741\pm 0.003}\;{f}_{\tau }^{-0.741}\;{10}^{-0.365z}.\end{eqnarray} \tag{ 6 }$$

Here ${f}_{\tau }\approx 2/3$ is a renormalization factor of optical depths introduced to multiply the simulated optical depths to match the mean flux decrement D_A. Thus, at $z=0.25$ and with ${f}_{\tau }\approx 2/3$, they found ${{\rm{\Delta }}}_{{\rm{b}}}=(38.9){N}_{14}^{0.741}$. From $({512}^{3},100{h}^{-1}\;\mathrm{Mpc})$ simulations with the HM01 background, Tepper-Garcia et al. (2012) found ${{\rm{\Delta }}}_{{\rm{b}}}=(48.3){N}_{14}^{0.786\pm 0.010}$, with no renormalization factor ${f}_{\tau }$. They noted a theoretical expectation that ${{\rm{\Delta }}}_{{\rm{b}}}\propto {N}_{{\rm{H}}\;{\rm{I}}}^{0.738}$ for absorbing gas in hydrostatic and photoionization equilibrium (Schaye 2001). In our $({768}^{3},50{h}^{-1}\;\mathrm{Mpc}$) simulations, also using the HM01 radiation field, we find ${{\rm{\Delta }}}_{{\rm{b}}}=(36.9){N}_{14}^{0.65}$ over the range $0.2\lt z\lt 0.3$ for column densities $12.5\lt \mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}\lt 14.5$. The overdensity is calculated as a weighted average over the line elements that make up the absorber, using the column density as the weight field. We find similar normalizations, ${{\rm{\Delta }}}_{0}=(36.9,38.9,48.3)$ for the three simulations (Davé et al. 2010; Tepper-Garcia et al. 2012; the current paper) evaluated at $z=0.25$ with the HM01 background.

Evidently, the structural measures of the low-z Lyα forest give similar results. The thermal-phase diagrams are essentially identical, while the differences in the $({{\rm{\Delta }}}_{{\rm{b}}},{N}_{{\rm{H}}\;{\rm{I}}})$ correlation likely arise from different methods of identifying and characterizing H i absorbers in column density and overdensity. These comparisons are complicated by the fact that Davé et al. (2010) corrected their simulated optical depths by a factor (${f}_{\tau }$), and Tepper-Garcia et al. (2012) corrected for the temperature–density correlation of gas in hydrodynamic and photoionization equilibrium. The simulations also used different cosmological parameters. Tepper-Garcia et al. (2012) adopted density fractions (${{\rm{\Omega }}}_{m}=0.238$ and ${{\rm{\Omega }}}_{{\rm{b}}}=0.0418$) from WMAP-3, which are lower by 15.6% and 9.1% than our WMAP-9 values (${{\rm{\Omega }}}_{m}=0.282$ and ${{\rm{\Omega }}}_{{\rm{b}}}=0.046$). Davé et al. (2010) adopted WMAP-7 parameters (${{\rm{\Omega }}}_{m}=0.28$ and ${{\rm{\Omega }}}_{{\rm{b}}}=0.046$) similar to our values. Kollmeier et al. (2014) used ${{\rm{\Omega }}}_{m}=0.25$ (10% lower than our value) and ${{\rm{\Omega }}}_{{\rm{b}}}=0.044$ (4.3% lower). The simulations also have different resolutions: our standard simulations use (768³, $50{h}^{-1}\;\mathrm{Mpc}$) while Tepper-Garcia use (512³, $100{h}^{-1}\;\mathrm{Mpc}$). Given the scatter in these correlations, their redshift dependence, and the correction factors, we believe  that the modest differences do not warrant concern.