The Low-redshift Lyα Forest as a Constraint for Models of AGN Feedback

The Lyα forest is a key diagnostic for the state of diffuse baryons in the intergalactic medium (IGM) and for fundamental cosmological parameters (Gunn & Peterson 1965; Palanque-Delabrouille et al. 2013; McQuinn 2016). At high redshift (i.e., z ≥ 2), the Lyα forest is observed in optical wavelengths from ground-based observatories and has been used to constrain small-scale cosmic structure (Croft et al. 1998; McDonald et al. 2005; Lidz et al. 2010), the IGM gas temperature (Becker et al. 2011; Boera et al. 2014), and the evolution of the ultraviolet ionizing background radiation (Haardt & Madau 1996; Faucher-Giguère et al. 2008, 2009; Haardt & Madau 2012; Faucher-Giguère 2020). In particular, a comparison between cosmological simulations of the Lyα forest and observations has shed light on dark matter and baryon interactions, IGM heating and cooling, and phase mixing in the IGM (Cen et al. 1994; Zhang et al. 1995; Miralda-Escudé et al. 1996; Hernquist et al. 1996; Rauch et al. 1997). At high redshifts (z = 2.0−5.4), the commonly studied statistical diagnostics of the Lyα forest such as the H i column density distribution (Altay et al. 2011; Rahmati et al. 2013), the Doppler widths (b parameter) of the Lyα absorption lines (Hiss et al. 2018), and the distribution and power spectrum of the transmitted flux (Walther et al. 2019; Chabanier et al. 2020) can be matched between simulations and observations to allow only a factor of 2 uncertainty in the amplitude of the metagalactic UV background (UVB).

Early theoretical studies of the low-redshift Lyα forest (e.g., Davé et al. 1999) demonstrated its potential for constraining the UVB. However, observing the Lyα forest at z ≤ 2 is substantially more challenging than at higher redshifts as it requires state-of-the-art spectrographs above Earth's atmosphere to reach the rest-frame far-ultraviolet (FUV) band (Danforth et al. 2016; Gurvich et al. 2017; Viel et al. 2017; Khaire et al. 2019; Christiansen et al. 2020). Recent observations with the Cosmic Origins Spectrograph (COS) aboard the Hubble Space Telescope have enabled studies of the low-redshift Lyα forest in great statistical detail and with high sensitivity (Meiring et al. 2011; Khaire et al. 2019; Danforth et al. 2016; Kim et al. 2020). In particular, Danforth et al. (2016, henceforth D16) built on the heritage of many previous low-z IGM absorber catalogs from HST/FOS (Bahcall et al. 1996; Jannuzi et al. 1998; Weymann et al. 1998), HST/GHRS (Penton et al. 2000), FUSE (Danforth & Shull 2005), and HST/STIS (Penton et al. 2004; Lehner et al. 2007; Tripp et al. 2008; Danforth et al. 2010; Tilton et al. 2012). The D16 survey represents the largest sample of low-z absorbers to date, comprising 2611 distinct redshift systems and 5138 intergalactic absorption lines. In addition to the statistical significance of this survey, the exquisite sensitivity of the COS instrument has allowed for the detection of very low-column-density absorbers in the local universe, down to column densities of N_{H i} ≈ 10^12.5 cm⁻².

The COS observations described in D16 have challenged our theoretical understanding of the Lyα forest in ways that have clearly called for significant further numerical study. In particular, forward-modeling cosmological simulations without galactic feedback or with current UVB models produces a simulated Lyα forest that does not match the low-redshift Lyα forest observed by COS. Kollmeier et al. (2014) reported the first discrepancy between the observed column density distribution function (CDD) of absorbers in D16 and the CDD derived using synthetic Lyα spectra from a simulation run with the smoothed particle hydrodynamics (SPH) code GADGET. Kollmeier et al. (2014) found that the number of absorbers in simulations at low column density is too high by a normalizing factor of 3.3 (equivalent to a factor of 5× fewer ionizing UV photons) and termed the difference the "photon underproduction crisis."

Gurvich et al. (2017, hereafter G17) first suggested that feedback from active galactic nuclei (AGNs) has an effect on the low-redshift Lyα forest in cosmological simulations, finding that the simulated CDD was sensitive to the overall strength of AGN feedback, in addition to the UVB. When sufficiently strong AGN feedback is included, the rate of ionizing UV photons need only be increased by a factor of 1.8× rather than 5×. These photons can be accounted for in the uncertainty between different UVB models, e.g., the difference at z = 0 between the UVB models of Faucher-Giguère et al. (2009) and Haardt & Madau (2012), which differ in their assumptions of the escape fraction of UV photons from star-forming galaxies. Thus, correctly modeling both the UVB and the AGN feedback in simulations is critical to matching the observed Lyα forest. Later studies using different hydrodynamic simulations further confirmed that black hole feedback can heat the IGM and is most impactful at z < 2 (Viel et al. 2017; Christiansen et al. 2020), while other studies have focused on further constraining the amplitude of the UV background in the empirical models (Khaire & Srianand 2015; Shull et al. 2015; Puchwein et al. 2018; Faucher-Giguère 2020; Chang et al. 2012).

Like the CDD, the distribution of absorption line Doppler widths (i.e., the b-parameter distribution) is consistently found to be mismatched between hydrodynamic simulations and observations at low redshift (Viel et al. 2017; Bolton et al. 2022). Viel et al. (2017) first explored the b distribution for the Illustris, Sherwood, and GADGET-3 simulations and found that all three simulations underpredict the number of absorbers at b = 25 to −45 km s⁻¹. At the same time, the simulations overpredict absorbers below 20 km s⁻¹. Using a similar analysis, Bolton et al. (2022) found that the simulated b distribution in the IllustrisTNG simulations is also systematically lower than what is observed.

It is now clear that, from the perspective of numerical simulations, modeling the low-redshift Lyα forest correctly is challenging because it requires not only (1) a precise understanding of hydrodynamics, heating, and cooling, but also (2) accurate subgrid models of the highly nonlinear physics that governs the interactions between the gas and galactic physics (e.g., AGN feedback). That the statistics of the Lyα forest at z < 2 have been shown to be more sensitive to gas on the outskirts of collapsed objects such as filaments and galaxies further suggests that galactic physics can impact observations of the diffuse IGM (Tonnesen et al. 2017; Martizzi et al. 2019). In addition to the integrated strength, the nuances of the AGN subgrid model matter substantially for the simulated Lyα forest properties. For example, the Sherwood cosmological simulations have generally found a limited impact of AGN feedback on the Lyα forest line widths, with most hot gas from AGN too hot and highly ionized for Lyα absorption (Viel et al. 2017; Nasir et al. 2017; Bolton et al. 2022). Their proposed solution is additional sources of turbulence in the IGM and heating from IGM dust. On the other hand, G17 and Christiansen et al. (2020) have found that warm, highly ionized gas produced by the jet feedback model in the SIMBA (Davé et al. 2019) and Illustris simulations (Springel et al. 2018) significantly improves agreement with the mean Lyα forest transmission at z < 0.5. Thus, the details of the subgrid model for AGN feedback can substantially alter the phase structure of the low-redshift IGM in cosmological simulations, even when controlling for differences in UVB models and hydrodynamic solver.

Here, we explore the low-redshift Lyα forest in the IllustrisTNG cosmological simulation suite compared to the original Illustris simulations in the context of the observed D16 COS data set in order to track the impact of different AGN feedback subgrid models given the same UVB, hydrodynamic solvers, and nearly identical cosmology.

The paper is organized as follows: In Section 2 we briefly describe the IllustrisTNG cosmological simulation suite, with a focus on the changes made to the AGN feedback subgrid model compared to the original Illustris simulations. In Section 3 we present our results comparing the statistics of the Lyα forest in IllustrisTNG to those of Illustris and the D16 data set. Specifically, we compare the column density distribution, distribution of the b parameter (i.e., the line widths), and flux probability distribution function and flux power spectrum. In Section 4 we discuss our results and summarize and conclude in Section 5.

In this section, we provide an overview of the Illustris and IllustrisTNG simulations used in this study, both of which are performed with the AREPO code (Springel 2010; Weinberger et al. 2020).

Illustris is a cosmological hydrodynamic simulation in a box of length 75 h⁻¹ Mpc. Gravitational interactions evolve via the TreePM algorithm (Springel et al. 2005). Radiative cooling is implemented using the network described in Katz et al. (1996) and includes line cooling, free–free emission, and inverse Compton cooling. Illustris assumes ionization equilibrium and accounts for on-the-fly hydrogen column density shielding from the radiation background (Rahmati et al. 2013). Metals and metal-line cooling are included (Vogelsberger et al. 2012, 2013), and star formation is implemented using the Springel & Hernquist (2003) subgrid model. Illustris uses Ω_m = 0.2726, Ω_Λ = 0.7274, Ω_b = 0.0456, σ₈ = 0.809, and H₀ = 100 h km s⁻¹ Mpc⁻¹ with h = 0.704. These parameters are consistent with the Wilkinson Microwave Anisotropy Probe 9 measurements (Hinshaw et al. 2013).

While the Illustris simulation suite has been remarkably successful in reproducing a number of observed galaxy properties (Vogelsberger et al. 2014a, 2014b; Sijacki et al. 2015), as well as properties of high-redshift neutral gas (Bird et al. 2014), its subgrid models for black hole feedback are in tension with a number of other observed galaxy properties (e.g., the low gas fraction of groups of galaxies as shown in Genel et al. 2014; Nelson et al. 2015). This tension partially motivated the next generation of the Illustris simulations, IllustrisTNG (TNG), which included an entirely updated subgrid model for AGN feedback developed in Weinberger et al. (2017), which we briefly describe below.

The TNG suite consists of 18 magnetohydrodynamic cosmological simulations, which vary in mass resolution, volume, and complexity of the physics included (Pillepich et al. 2018a; Marinacci et al. 2018; Naiman et al. 2018; Springel et al. 2018; Nelson et al. 2018, 2019a; Pillepich et al. 2019; Nelson et al. 2019b). Three box sizes are employed with cubic volumes of 51.7, 110.7, and 302.6 Mpc side length, referred to as TNG50, TNG100, and TNG300, respectively. TNG's default cosmological parameters are consistent with the Planck Collaboration et al. (2016) results, with matter density Ω_m = Ω_dm + Ω_b = 0.3089, baryonic density Ω_b = 0.0486, cosmological constant Ω_Λ = 0.6911, Hubble constant H₀ = 100 h km s⁻¹ Mpc⁻¹ with h = 0.6774 and normalization σ₈ = 0.8159. The largest box simulation (TNG300) allows for the study of rare objects/events and galaxy clustering and provides the largest galaxy sample while the smallest physical volume simulation (TNG50) allows for the study of more detailed structural properties of galaxies and their kinematics as well as the opportunity to do resolution convergence testing. TNG100 (also referred to in this work as simply TNG) uses the same initial conditions as the original Illustris simulation, allowing for direct comparison with the original Illustris and testing of the new physics included in the TNG runs. Each of the simulation boxes has been run at three resolution levels, allowing for further convergence studies.

The main differences between TNG and Illustris are (I) a new implementation of stellar-driven galactic winds with modified velocity, thermal content, and energy scalings (Pillepich et al. 2018b), (II) a new implementation of black-hole-driven kinetic feedback at low accretion rates, and (III) the inclusion of magnetohydrodynamics. G17 have already shown that galactic winds (as implemented in these models) have essentially no impact on the low-redshift Lyα forest properties, as their effects cannot travel to substantial distances outside of halos into the IGM, so we do not expect these changes in the physics model to affect the interpretation of our results.

Of great importance for the statistics of the Lyα forest is the implementation of the ultraviolet background (UVB) and the related photoheating rate. Illustris and TNG both use the UVB model presented in Faucher-Giguère et al. (2009). Changes to the UVB model and their effects on the IGM at z < 0.5 have been explored in other works (e.g., Shull et al. 2015; Gaikwad et al. 2017; Khaire et al. 2019; Bolton et al. 2022), therefore, we focus on studying the differences between the AGN feedback subgrid models. This direct comparison is made possible by the fact that TNG keeps the simulation properties known to affect the Lyα forest the same as they were in Illustris (e.g., the UVB model) except for the redesign of the AGN feedback model at low black hole accretion rates.

In the original Illustris, the AGN feedback subgrid model is split into a two-accretion-mode scenario in which feedback energy is injected by a radio mode at low accretion rates (relative to the Eddington limit) and a radiative mode at high accretion rates. The radio mode inflates a hot explosive bubble in the circumgalactic medium whose radius scales with the energy and density of the gas; however, this resulted in too-low gas fractions in galaxy groups and low-mass clusters due to an over-efficient expulsion of the gas (Genel et al. 2014). Additionally, the mode generated too high stellar masses in the central galaxies. Further details of the Illustris AGN model are described in Genel et al. (2014) and Vogelsberger et al. (2013). By contrast, the AGN feedback model employed in TNG assumes an updated two-accretion-mode scenario in which the feedback energy is injected through (1) a relatively efficient kinetic mode at low accretion rates, and (2) a less efficient thermal mode at high accretion rates in which thermal energy cools rapidly, which is similar to that in Illustris. The kinetic mode in TNG imparts a randomly oriented momentum boost (i.e., isotropic when integrated over time) to the gas cells in the few-kiloparsec feedback region around the black hole (with no mediation by wind particles) in a series of discrete injection events, which occur once the accumulated energy output from the black hole has exceeded a threshold. For more details on the TNG AGN model, see Pillepich et al. (2018a) and Zinger et al. (2020).

In Figure 1 (leftmost panels), we show a visual comparison of the z = 0.1 IGM column density (left column) and temperature (right column) in Illustris (top row) and TNG (bottom row). The differences seen are primarily due to the changes in the AGN jet feedback model. Hot ionized medium (HIM) gas above 10⁷ K (traced in X-ray wavelengths; red color in the temperature color bar) dominates in Illustris while in TNG the IGM volume has more Lyα -emitting gas (gas around 10⁴ K, blue color in the temperature color bar). In both simulations, the Warm-Hot Ionized Medium (WHIM; 10⁵ K < T < 10⁷ K, n_H < 10⁻⁴(1+z) cm⁻³) e.g., Davé et al. 2001) gas fills the volume many megaparsecs away from the galaxies. However, this gas is too hot to effectively absorb in Lyα. Interestingly, at z = 2.0, the IGM column density and temperature are much more similar between Illustris and TNG, as shown in Figure 1 (i.e., rightmost panels). As opposed to the situation at z = 0.1, much of the volume away from the central galaxies is cold (T ≲ 10⁴K) and the thermal state of the IGM is controlled by the UVB. A comparison of these figures not only demonstrates the ability of AGN feedback to impact gas out to many megaparsecs (see also Gurvich et al. 2017; Borrow et al. 2020; Christiansen et al. 2020) but specifically that the strength and/or manner in which the AGN feedback is implemented matters substantially in determining the phase structure of the IGM. In particular, the bubble model in Illustris produces copious amounts of WHIM. TNG instead has the updated kinetic AGN feedback model at low black hole accretion rates, which produces less extremely hot gas in the IGM.

Zoom In Zoom Out Reset image size

In this section, we present a comparison of the column density distribution (CDD), the distribution of the Doppler b parameter, the transmitted flux probability distribution function (PDF), and the power spectrum between the Hubble COS spectra, TNG, and Illustris.

In order to generate realistic spectra from the simulations, we use the publicly available code outlined in Bird et al. (2015) and Bird (2017). ⁷ Similar to G17, we estimate the CDD from the column density along 5000 simulated sight lines randomly oriented throughout the simulation box. Column densities are computed by interpolating the neutral hydrogen mass in each gas element to the sight line using an SPH kernel. For the column density, the sight line is split into absorbers every 50 kms⁻¹ (we confirm that our results are not sensitive to this value). In addition to this "direct" CDD method, we also compute the Voigt profile fits to the optical depth of our simulated spectra (using a fitting algorithm included in the fake_spectra package used to generate the simulated spectra) to obtain best-fit column densities and b parameters. G17 showed that there is no substantial difference between the Voigt and "direct" CDD calculations, and therefore, we will omit such a comparison here.

3.1. The Column Density Distribution

We define the CDD function:

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}N}{\mathrm{dlog}(N){\rm{d}}z}=\displaystyle \frac{F(N)}{{\rm{\Delta }}\mathrm{log}(N)}{\rm{\Delta }}z\end{eqnarray} \tag{ 1 }$$

F(N) is the number of absorbers per sight line with column density in the interval [N_{H i} ;N_{H i} + dN_{H i} ], over the redshift interval Δz contained in our simulated spectra and ${\rm{\Delta }}\mathrm{log}(N)$ is the logarithmic bin width.

Figure 2 shows a comparison of the z = 0.1 Lyα forest CDD of Illustris (i.e., the results of G17; solid red line), new measurements of the same from TNG (solid black line), and the D16 data set (blue points). We also include Illustris modified with a UVB field strength 0.4 times the original (dashed red line). We find that neither Illustris nor the updated TNG simulations can reproduce the COS CDD of the low-redshift Lyα forest. As reported by G17, we find that the Illustris CDD (solid red line) matches the D16 data well in the column density range of N_{H i} = 10¹²–10¹⁴cm⁻². TNG (solid black line) has too many absorbers in this column density range to match the observations, indicated by a higher amplitude in the CDD. Ultimately, we find that the CDD slope of both TNG and Illustris does not match the D16 data. In particular, the slope of both TNG and Illustris is too steep to match the COS data.

Zoom In Zoom Out Reset image size

We find that AGN feedback can shift the amplitude of the CDD, in a similar way to when altering the UVB. We demonstrate this by modifying the Illustris Lyα flux with UVB field strength 0.4 times the original Illustris simulation UVB (dashed red line). This is done in postprocessing by recalculating the neutral hydrogen fractions using scaled UVB values. We choose the value of 0.4 in order to produce an Illustris CDD, which matches well with TNG in the column range of $\mathrm{log}\ {N}_{\mathrm{HI}}\lt 15$. There remains a slight difference in the slope of TNG and Illustris CDD at $\mathrm{log}\ {N}_{\mathrm{HI}}\gt 15$; however, in this saturated regime the determination of column densities is more sensitive to details of the line-fitting method and spectral resolution, and it is subject to greater statistical uncertainties. Viel et al. (2017) investigated altering the UVB to force a match between Illustris (and other simulations such as the Sherwood simulations) and the COS data, finding that doing so requires an aggressive change in the photoionization rate of a factor of 1.5–3 times higher than that of Haardt & Madau (2012). However, the increase in required photons for TNG is not as extreme as that for models that lack AGNs, e.g., the simulations used in Kollmeier et al. (2014) required a factor of 5 more UV photons than those in Haardt & Madau (2012). Including AGN feedback is a possible solution to requiring a severe increase in the ionization rate as it contributes to heating the gas and therefore reduces the temperature dependent recombination rate coefficient until the gas is no longer able to absorb in Lyα. Therefore, it is likely that the full resolution of the discrepancy between the simulated and observed CDDs will require both an alteration in UVB and a reworked AGN feedback model that is tuned to both galaxy halo properties and properties of the IGM.

3.2. The Doppler b Parameter

We show the Doppler b parameter distribution for Illustris (red lines) and TNG (black lines) in Figure 3. In the bottom-right panel, we show the entire range of column density values we resolve (12.5 $\lt \ \mathrm{log}\ {N}_{\mathrm{HI}}\lt 15$ cm⁻²). Note that fitting Voigt profiles at column densities lower than ${\text{}}{\rm{log}}\ {N}_{\mathrm{HI}}\lt 12.5$ is difficult to perform accurately, while column densities higher than $\mathrm{log}\ {N}_{\mathrm{HI}}\gt 15$ become increasingly rare and suffer from low number statistics. When considered over the entire column density range, Illustris and TNG show almost identical b distributions peaked at b ≈ 20 km s⁻¹ with a long tail to large values. For reference, b = 22 km s⁻¹ is equivalent to thermally broadened lines with T = 10^4.5 K. The similarity in the two b distributions suggests that, when considered over the entire column density range, the b distribution for the Lyα forest is largely insensitive to the subgrid AGN feedback model.

Zoom In Zoom Out Reset image size

Other panels in Figure 3 consider low-, intermediate-, and high-column-density subsets of the full column density range to determine if there is a regime in which the effects of different subgrid AGN feedback models can be distinguished. At the lowest column density range considered (top left; 12.5 $\lt \ \mathrm{log}\ {N}_{\mathrm{HI}}\lt 13$ cm⁻²), Illustris tends to have slightly higher b values than TNG but both are still peaked at 20 km s⁻¹. Higher b values in Illustris are expected because its AGN radio feedback model produces an IGM with more volume-filling hot bubbles than TNG (see Figure 1). At the opposite end of the spectrum at higher column density values (bottom left; N_{H i} = 10¹⁴–10¹⁵cm⁻²), both Illustris and TNG are peaked at 30 km s⁻¹ though the b distribution for Illustris is slightly less skewed than TNG. The skew in TNG is likely due to shock heating from the AGN kinetic mode, which can affect these higher IGM column densities (see Figure 1).

In the intermediate column density range (top right; N_{H i} = 10¹³–10¹⁴ cm⁻²) Illustris and TNG are identical. As a point of comparison, we also overplot the COS observations analyzed in D16 (blue points) in the same column density range. It is clear that, as in other studies, both Illustris and TNG produce b distributions much lower than what is observed (Gaikwad et al. 2017; Nasir et al. 2017; Bolton et al. 2022). Though the simulated spectra have not been corrected to match the spectral resolution of COS (FWHM = 15 km s⁻¹), the differences seen between simulation and observation in Figure 3 are too large for this to be the sole solution. One possible reason for the difference in b values between the simulations and observations is that the simulations are underresolving the amount of turbulence present in the IGM or are missing sources of turbulence, which we discuss further in Section 4.2.

Overall, though there are small differences at the lowest and highest column densities we resolve, we find that the overall b distribution is not affected by differences in the subgrid AGN feedback model employed. We note that the Illustris b distribution is largely unchanged when altering the UVB.

3.3. The IGM Density–Temperature Relation

To further understand the effects of the subgrid AGN feedback model on the IGM temperature and column density, we plot the 2D temperature–density relation (otherwise known as the equation of state) for IGM gas in Figure 4. In simulations of the diffuse IGM, the temperature and density are known to satisfy a power-law relation over two decades in density, T(Δ) = T₀Δ^γ−1, where Δ = ρ/ρ₀ is the overdensity factor, γ is the power-law index, and T₀ is the temperature at the mean density (Katz et al. 1996; Hernquist et al. 1996; Davé et al. 1999). This power-law relation quantifies the thermal state of the IGM (Hui & Gnedin 1997; McQuinn 2016).

Zoom In Zoom Out Reset image size

In Figure 4 we show the phase diagrams for TNG (far left), Illustris (far right), and the ratio of the two (center). The ratio plot is colored by the gas mass ratio of TNG to Illustris, i.e., red identifies regions of temperature–density phase space more populated in TNG and blue those more populated in Illustris. White color indicates the two simulations are the same in that part of phase space. The power-law region of the diffuse IGM (10³ K ≤ T ≤ 10^4.5 K; 10⁻⁸ cm⁻³ ≤ n ≤ 10⁻⁵ cm⁻³) is nearly identical between the two simulations. This should be expected as the b-value distributions are very similar between the two simulations. However, looking at the middle panel of Figure 4 around the power law, the IGM gas appears to be warmer in Illustris than in TNG. Given that these two simulations employ the same UVB, another source of heating is needed to explain the temperature shift near the power-law relation. The additional heating is likely due to the hot AGN bubble model in Illustris.

Looking at higher temperatures, more gas in Illustris is pushed into the HIM/WHIM phase state as compared to TNG: Martizzi et al. (2019) found that the mass of the WHIM in Illustris is ∼26% larger than in TNG, consistent with this excess. However, this excess is concentrated at either high density (i.e., in the circumgalactic gas) or at the very lowest densities. In Figure 4, intermediate densities ($-4\leqslant \mathrm{log}(n)\leqslant -7$ cm⁻³) are white in color, suggesting little to no difference between TNG and Illustris here. Illustris has far more HIM than TNG. However, hot gas above 10⁵K does not absorb in Lyα and is instead probed better by ionized gas tracers beyond the scope of this study (e.g., ionized metal lines and X-rays). We note an excess of T = 10⁴K gas at high density in Illustris that is absent in TNG. This gas is halo/ISM gas that is not relevant for the IGM and is discussed further in Martizzi et al. (2019).

3.4. The Flux PDF and Power Spectrum

The amplitude of the Lyα forest power spectrum is sensitive to the UV background and can be used to measure the H i photoionization rate Γ_{H i} . Recently, the power spectrum of the low-redshift Lyα forest was estimated from the COS data of Danforth et al. (2016) and then used to estimate the photoheating rate (Khaire et al. 2019). Based on the power spectrum, they found a Γ_{H i} that agreed with previous estimates (Khaire & Srianand 2019; Puchwein et al. 2018; Shull et al. 2015). Their analysis suggested that the low-z UV background is dominated by ionizing photons emitted by quasars and does not require any significant contribution from galaxies. In this section, we investigate if the flux statistics contain dependencies on the AGN feedback model.

We present the flux power spectrum of the z = 0.1 Lyα forest in TNG and Illustris in Figure 5. We overplot the COS data presented in Khaire et al. (2019) in blue. Compared to the observations, TNG is in much better agreement than Illustris for 8.8 × 10⁻⁴ < k < 0.07 s km⁻¹. Both simulations and observations show evidence of a thermal cutoff in the power at small scales (k > 0.02 s km⁻¹), which is a result of pressure smoothing of the IGM and thermal broadening of absorption lines. This cutoff is an important feature that probes the thermal state of the IGM (Zaldarriaga et al. 2001; Walther et al. 2019). The cutoff is shifted toward smaller scales in the simulations as compared to the data, reflecting the lower gas temperatures in the simulation. In particular, TNG produces more power at k > 0.07 s km⁻¹ when compared to COS. TNG and Illustris are different in both shape and normalization. The amplitude of the flux power spectrum in Illustris is significantly lower than that in TNG. However, these differences can be reduced by using Illustris with a rescaled UVB selected to match the CDD (i.e., the dashed red line in Figure 5). However, Illustris produces excess power on the largest scales probed, k < 2 × 10⁻⁴ s km⁻¹, which could be due to large-scale correlations induced by AGN heating. From this, we can conclude that the difference manifested by the altered feedback models in these simulations is not solely in the temperature distribution of the gas, but also in the fact that the physical scale of AGN feedback in Illustris is large enough to change the power spectrum significantly, even when the mean absorption is matched. The mismatch of Illustris with the COS data (blue points) suggests that the AGN feedback in Illustris is unrealistically strong. This conclusion was also previously discussed based on a mismatch of galaxy properties (e.g., low gas fraction of groups of galaxies as shown in Genel et al. 2014; Nelson et al. 2015), and our results suggest less extreme AGN feedback models, such as TNG, could also be tested using the Lyα flux power spectrum.

Zoom In Zoom Out Reset image size

We emphasize that one cannot attribute the amplitude and shape differences in the z = 0.1 power spectrum between Illustris and TNG to the photoheating and photoionization rates as these rates are identical (both simulations use the Faucher-Giguère et al. 2009 UVB). This suggests that differences in the subgrid AGN or stellar feedback models are responsible for the differences seen in the power spectrum. We note that G17 finds that the Lyα forest is largely insensitive to stellar feedback. Furthermore, our results suggest N-body simulations that only explore changes in photoionization rates and do not include AGN feedback processes do not capture all the subtleties that can affect the flux power spectrum and simply rescaling the power spectrum amplitude using a different UVB may not capture all the necessary IGM physics relevant at low redshift.

In addition to the power spectrum, we present the flux PDF of Illustris (red solid line) and TNG (black solid line) in Figure 6. We also include the Illustris UVBx0.4 comparison, which matches the CDD of TNG (dashed red line). As expected from the CDD, the flux PDF amplitude in Illustris is lower than TNG until the optically thin regime. Furthermore, we see that simply altering the UVB is not able to produce matching flux PDFs from the two simulations, as it is sensitive to absorbers with column densities too small to be Voigt fit, which are thus not included in the column density function. As was the case with the power spectrum, the noticeable differences between the two simulations suggest that the flux PDF may be a useful tool to constrain subgrid models of AGN feedback.

Zoom In Zoom Out Reset image size

4.1. The Importance of AGN Feedback on the Low-z Lyα Forest

4.2. Missing Turbulence in the IGM

In this paper, we have studied the low-redshift Lyα forest statistics in the Illustris and TNG cosmological simulations in comparison to the COS UV data. The TNG simulation has an identical UVB prescription and very similar cosmological parameters as Illustris but an entirely reworked AGN feedback prescription. We find that:

• 1.  
  Due to the AGN radio-mode model, the original Illustris simulations have a factor of 2–3 fewer absorbers than TNG at column densities N_{H i} < 10^15.5 cm⁻². Neither TNG nor Illustris can match the slope of the COS CDD.
• 2.  
  The b distributions in Illustris and TNG are similar, suggesting it is a poor barometer for evaluating the effects of subgrid AGN feedback models.
• 3.  
  In agreement with Bolton et al. (2022), we find that both Illustris and TNG produce smaller b values than the COS observations, suggesting the need for additional sources of nonthermal broadening in the low-redshift IGM.
• 4.  
  The amplitude and shape differences in the flux power spectrum and flux PDF between Illustris and TNG are significant (Figures 5 and 6) and cannot be due solely to changes in the UVB model. We find that the flux statistics are sensitive diagnostic constraints for subgrid AGN feedback models.
• 5.  
  The new kinetic mode AGN feedback model in TNG substantially improves agreement with the observed flux power spectrum; however, TNG still overpredicts the high-k power observed in COS.

Ultimately, resolving the discrepancies between the cosmological simulations and the HST COS data will lead to the development of an improved UVB and AGN feedback model that is based on both the properties of galaxies (e.g., Weinberger et al. 2017) and the properties of the Lyα forest. Developing an updated AGN feedback model based on a match of simulated and observed IGM properties and deepening our understanding of how AGN feedback alters the thermal and turbulent properties of neutral gas in the IGM will require a large campaign of numerical work that spans a range of parameters and codes. Fortunately, such a simulation campaign has already been carried out as part of the CAMELS project: Cosmology and Astrophysics with MachinE Learning Simulations (Villaescusa-Navarro et al. 2020). CAMELS includes a suite of 2184 state-of-the-art magnetohydrodynamic cosmological simulations performed in a volume of 25 Mpc. The simulations are run either with the AREPO or GIZMO codes and employ the same baryonic subgrid physics as the TNG and SIMBA simulations (Davé et al. 2019). CAMELS contains thousands of different cosmological and astrophysical models by way of varying cosmological parameters and, most importantly for Lyα forest studies at z = 0.1, the parameters controlling stellar and AGN feedback and the UVB. We will use these large AGN feedback parameter suites in future work constraining the nature of the effects of AGN feedback on the IGM (M. Tillman et al. 2022, in preparation).

We investigate how resolution and simulation volume size affects the Lyα forest statistics presented above.

Figure 7 presents the Lyα statistics of TNG50-1 (red lines), TNG100-1 (black lines), and TNG300-1 (blue lines). Each of these simulations is the highest-resolution realization but they have different simulation box sizes, with TNG50-1 being the smallest volume run and TNG300-1 being the largest. The left-most plot of Figure 7 presents the CDD. Below column densities of 10¹⁵ cm⁻² we do not find a substantial difference between the different volume runs of TNG. However, higher column densities do show differences in the CDD. As expected, these high column densities represent rare overdensities around halos and TNG50 departs most significantly from TNG100 and TNG300. Our study primarily focuses on the IGM regime of column densities (i.e., below 10¹⁶ cm⁻²). Differences between the CDDs in the column density range of interest are negligible and due mostly to sample variance, demonstrating that our results are insensitive to changes in both the box size and in resolution down to 2 × 10³ resolution elements.

Zoom In Zoom Out Reset image size

The second and third panels of Figure 7 present the Doppler parameter distributions and flux PDFs respectively. Neither of these statistics shows significant variation, implying simulation volume is not important here. The rightmost panel presents the flux power spectrum. On all scales, sensitivity of the flux power spectrum to resolution and volume is substantially less than both observational errors and the effect of the subgrid feedback model. Large scales (k ≲ 10⁻³ s/km) do however exhibit effects from the finite volume of the TNG50 simulation, although they are reasonably well converged for the TNG100 simulation we use for our main results.

Figure 8 shows the effects of changing the number of resolution elements in TNG100. TNG100-1 (black solid lines) matches the original Illustris box size of 75 Mpc h⁻¹ and 2 × 1820³ resolution elements. TNG100-2 (black dashed lines) has the same volume size but decreases the resolution elements to 2 × 910³. TNG100-3 (dotted black lines) has the same volume as TNG100-1 and TNG100-2 but decreases the resolution elements to 2 × 455³.

Zoom In Zoom Out Reset image size

The leftmost plot of Figure 8 shows the CDD for varying resolutions. Again, there are variations in the CDD between different resolution runs at column densities greater than log 15 cm⁻² as expected for these rare absorbers. TNG100-2 and TNG100-1 are fairly close, demonstrating that there is convergence in this statistic with the explored resolutions. TNG100-3 produces noticeable differences in the CDD across the full range of column densities, implying that the lowest resolution is not sufficient for analyzing the CDD. TNG100-3 produces fewer absorbers at all column densities, similar to the results for the Illustris simulation suite found in G17. Convergence is similarly observed for the b-value distribution, the flux PDF, and the flux power spectrum.

While all the TNG50 and TNG100 low redshift Lyα statistics appear to converge, the TNG300 runs statistics do not converge. Figure 9 is the same as Figure 8 but for TNG300 runs instead. It is clear that the resolutions explored in the TNG300 runs are not sufficient to resolve the low redshift Lyα forest statistics. The Doppler parameter distribution seems particularly sensitive to resolution.

Zoom In Zoom Out Reset image size