The low redshift Lyman-$\alpha$ Forest as a constraint for models of AGN feedback

We study the sensitivity of the $z=0.1$ Lyman-$\alpha$ Forest observables, such as the column density distribution function (CDD), flux PDF, flux power spectrum, and line width distribution, to sub-grid models of active galactic nuclei (AGN) feedback using the Illustris and IllustrisTNG (TNG) cosmological simulations. The two simulations share an identical Ultraviolet Background (UVB) prescription and similar cosmological parameters, but TNG features an entirely reworked AGN feedback model. Due to changes in the AGN radio mode model, the original Illustris simulations have a factor of 2-3 fewer Lyman-$\alpha$ absorbers than TNG at column densities $N_{\rm HI}<10^{15.5}$ cm$^{-2}$. We compare the simulated forest statistics to UV data from the Cosmic Origins Spectrograph (COS) and find that neither simulation can reproduce the slope of the absorber distribution. Both Illustris and TNG also produce significantly smaller line width distributions than observed in the COS data. We show that TNG is in much better agreement with the observed $z=0.1$ flux power spectrum than Illustris. We explore which statistics can disentangle the effects of AGN feedback from alternative UVB models by rescaling the UVB of Illustris to produce a CDD match to TNG. While this UVB rescaling is degenerate with the effect of AGN feedback on the CDD, the amplitude and shape of the flux PDF and 1D flux power spectrum change in a way distinct from a scaling of the UVB. Our study suggests that the $z=0.1$ Lyman-$\alpha$ forest observables can be used as a diagnostic of AGN feedback models.


INTRODUCTION
The Lyman-α 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 Lyman-α 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, 2009Haardt & Madau 2012;Faucher-Giguère 2020). In particular, comparison between cosmological simulations of the Lyman-α forest and observations have 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 Lyman-α forest such as the HI column density distribution (Altay et al. 2011;Rahmati et al. 2013a), the Doppler widths (b-parameter) of the Lyman-α absorption lines (Hiss et al. 2018), and the distribution and power spectrum of the transmitted flux (Walther et al. 2019a;Chabanier et al. 2020) can be matched between simulations and observations to allow only a factor of two uncertainty in the amplitude of the metagalactic UV background (UVB).
Early theoretical studies of the low-redshift Lymanα forest (e.g. Davé et al. (1999)) demonstrated its potential for constraining the UVB. However, observing the Lyman-α forest at z ≤ 2 is substantially more challenging than at higher redshifts as it requires state-of-theart 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;Christiansen et al. 2020). Recent observations with the Cosmic Origins Spectrograph (COS) aboard the Hubble Space Telescope have enabled studies of the low-redshift Lyman-α forest in great statistical detail and with high sensitivity (Meiring et al. 2011;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 detection of very low column density absorbers in the local universe, down to column densities of N HI ≈ 10 12.5 cm −2 .
The COS observations described in D16 have challenged our theoretical understanding of the Lyman-α 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 Lyman-α forest that does not match the low redshift Lyman-α 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 Lyman-α 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 5x 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 (AGN) has an effect on the low-redshift Lyman-α 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.8x rather than 5x. 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 andMadau 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 Lyman-α forest. Later studies using different hydrodynamic simulations further confirmed that black hole feedback can heat the IGM, and is most impactful at z < 2 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 Bolton et al. 2021). 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 − 45km s −1 . At the same time, the simulations over predict absorbers below 20 km s −1 . Using a similar analysis, Bolton et al. (2021) found that the simulated b distribution in the IllustrisTNG simulations are also systematically lower than what is observed.
It is now clear that, from the perspective of numerical simulations, modeling the low-redshift Lyman-α forest correctly is challenging because it requires not only: 1) a precise understanding of hydrodynamics, heating, and cooling, but also 2) accurate sub-grid models of the highly non-linear physics that governs the interactions between the gas and galactic physics (e.g. AGN feedback). That the statistics of the Lyman-α 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 sub-grid model matter substantially for the simulated Lyman-α forest properties. For example, the Sherwood cosmological simulations have generally found a limited impact of AGN feedback on the Lyman-α forest line widths, with most hot gas from AGN too hot and highly ionized for Lyman-α absorption Nasir et al. 2017;Bolton et al. 2021). 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 ) significantly improves agreement with the mean Lyman-α forest transmission at z < 0.5. Thus the details of the sub-grid 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 Lyman-α forest in the IllustrisTNG cosmological simulation suite as com-pared to the original Illustris simulations in the context of the observed D16 COS dataset in order to track the impact of different AGN feedback sub-grid 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 sub-grid model compared to the original Illustris simulations. In Section 3 we present our results comparing the statistics of the Lyman-α 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.

NUMERICAL SIMULATIONS
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).
While the Illustris simulation suite has been remarkably successful in reproducing a number of observed galaxy properties (Vogelsberger et al. 2014a,c;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., 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 sub-grid 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. 2018Nelson et al. , 2019aPillepich 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, cosmolog-ical constant Ω Λ = 0.6911, Hubble constant H 0 = 100 h km s −1 Mpc −1 with h = 0.6774, normalisation σ 8 = 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 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 has already shown that galactic winds (as implemented in these models) have essentially no impact on the lowredshift Lyman-α 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 Lymanα 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;Bolton et al. 2021) therefore we focus on studying the differences between the AGN feedback sub-grid models. This direct comparison is made possible by the fact that TNG keeps the simulation properties known to affect the Lyman-α 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 sub-grid 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 energy and density of the gas, however this resulted in too-low gas fractions in galaxy groups and low-mass clusters due to an overefficient expulsion of the gas ). 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 twoaccretion-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 . At z=2.0, there are very minor differences between Illustris and TNG but by z=0.1 significant visual differences are obvious in the temperature and density of the IGM. The differences seen here are nearly entirely due to the changes in the AGN jet feedback model. Hot ionized medium (HIM) gas above 10 7 K (traced in x-rays wavelengths; red color in the temperature color-bar) dominates in Illustris while the IGM volume is more filled by warm Lyman-α emitting gas in TNG (gas around 10 4 K (blue color in the temperature color-bar).
gas cells in the few kpc 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); Zinger et al. (2020).
In Figure 1 (left most 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. HIM gas above 10 7 K (traced in x-rays wavelengths; red color in the temperature color-bar) dominates in Illustris while in TNG the IGM volume has more Lyman-α emitting gas (gas around 10 4 K, blue color in the temperature color-bar). In both simulations, the Warm-Hot Ionized Medium (WHIM; 10 5 K < T < 10 7 K, n H < 10 −4 (1+z) cm −3 ) e.g. Davé et al. (2001)) gas fills the volume many Mpcs away from the galaxies. However, this gas is too hot to effectively absorb in Lyman-α. 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., right most panels). As opposed to the situation at z = 0.1, much of the volume away from the central galaxies is cold (T 10 4 K) and the thermal state of the IGM is controlled by the UVB. Comparison of these figures not only demonstrates the ability of AGN feedback to impact gas out to many Mpcs (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.

RESULTS
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 power spectrum between 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  1 . Similar to G17, we estimate the CDD from the column density along 5,000 simulated sightlines randomly oriented throughout the simulation box. Column densities are computed by interpolating the neutral hydrogen mass in each gas element to the sightline using an SPH kernel. For the column density, the sightline is split into absorbers every 50 kms −1 (we confirm that our results are not sensitive to this value). In addition to this "direct" CDD method, we also compute 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.

The Column Density Distribution
We define the CDD function: redshift interval ∆z contained in our simulated spectra and ∆ log(N ) is the logarithmic bin width. Figure 2 shows a comparison of the z = 0.1 Lymanα 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 Lyman-α 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 HI = 10 12 -10 14 cm −2 . 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.
We find that AGN feedback can shift the amplitude of the CDD, in a similar way as when altering the UVB. We demonstrate this by modifying the Illustris Lymanα flux with UVB field strength 0.4 times the original Illustris simulation UVB (dashed red line). This is done in post processing 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 logN HI < 15. There remains a slight difference in the slope of TNG and Illustris CDD at logN HI > 15, however in this saturated regime the determination of column densities is more sensitive to details of the line fitting method, 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 Haardt & Madau (2012). However, the increase in required photons for TNG is not as extreme as models that lack AGN, e.g., the simulations used in Kollmeier et al. (2014) required a factor of 5 more UV photons than 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 Lyman-α. 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 < logN HI < 15 cm −2 ). Note that fitting Voigt profiles at column densities lower than logN HI < 12.5 is difficult to perform accurately, while column densities higher than logN HI > 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 −1 with a long tail to large values. For reference, b = 22 km s −1 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 Lyman-α Forest is largely insensitive to the sub-grid AGN feedback model.
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 sub-grid AGN feedback models can be distinguished. At the lowest column density range considered (top left; 12.5 < logN HI < 13 cm −2 ) Illustris tends to have slightly higher b values than TNG but both are still peaked at 20 km s −1 . 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 HI = 10 14 -10 15 cm −2 ), both Illustris and TNG are peaked at 30 km s −1 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 HI = 10 13 -10 14 cm −2 ) 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. 2021). 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 sub-grid AGN feedback model employed. We note that the Illustris b distribution is largely unchanged when altering the UVB.

The IGM Density-Temperature Relation
To further understand the effects of the sub-grid 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 0 ∆ γ−1 where ∆ = ρ/ρ 0 is the overdensity factor, γ is the power law index, and T 0 is the temperature at the mean density Hernquist et al. 1996;Davé et al. 1999). This power-law relation quantifies the thermal state of the IGM (Hui & Gnedin 1997;McQuinn 2016).
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 3 K ≤ T ≤ 10 4.5 K; 10 −8 cm −3 ≤ n ≤ 10 −5 cm −3 ) 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 ≤ log(n) ≤ −7 cm −3 ) 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 5 K does not absorb in Lymanα 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 4 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).

The Flux PDF and Power Spectrum
The amplitude of the Lyman-α forest power spectrum is sensitive to the UV background and can be used to measure the HI photoionization rate Γ HI . Recently, the power spectrum of the low redshift Lyman-α forest was estimated from the COS data of Danforth et al. (2016) and then used to estimate the photoheating rate . Based on the power spectrum, they found a Γ HI 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 Lyman-α forest in TNG and Illustris in Figure 5. We overplot the COS data presented in  in blue. Compared to the observations, TNG is in much better agreement than Illustris for 8.8 × 10 −4 < k < 0.07 s/km. Both simulations and observations show evidence of a thermal cut-off in the power at small scales (k > 0.02 skm −1 ), which is a result of pressure smoothing of the IGM and thermal broadening of absorption lines. This cut-off is an important feature that probes the thermal state of the IGM (Zaldarriaga et al. 2001;Walther et al. 2019b). The cutoff is shifted towards 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 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 −4 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 Lyman-α flux power spectrum.
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 sub-grid AGN or stellar feedback models is responsible for the differences seen in the power spectrum. We note that G17 finds that the Lyman-α 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 sub-grid models of AGN feedback.

The importance of AGN feedback on the low-z
Lyman-α forest AGN feedback is crucial for reproducing galaxy properties (Somerville & Davé 2015), but its effect on the physical state of the low redshift IGM is not well understood. The last decade has witnessed the development of various AGN feedback models for cosmological hydrodynamic simulations that are tuned to quench galaxies in agreement with observations (Vogelsberger et al. 2014b;Schaye et al. 2015;Weinberger et al. 2017;Henden et al. 2018). It has only recently become clear that AGN feedback can potentially carry matter and energy far from the central regions of galaxies and have a substantial effect on the thermal state of the IGM (Borrow et al. 2020;Christiansen et al. 2020). In this paper we have demonstrated that AGN heating is impactful at low redshift and the non-linear growth of structure driving the WHIM and the diffuse IGM are subject to the details of the subgrid feedback-physics implementation. TNG and Illustris employ nearly identical cosmological parameters and numerical schemes but differ dramatically in their stellar and black hole feedback prescriptions. The sensitivity of the low-redshift Lyman-α forest to hydrodynamical effects emphasizes the importance of AGN feedback physics and related effects, which makes the low-redshift Lyman-α forest an ideal test-bed for realistic hydrodynamical simulations. In addition to the CDD, the differences between Illustris and TNG in the flux power spectrum and PDF are significant and cannot be attributed to changes in the photoheating rate.
The gas in recent numerical simulations is either too hot to contribute to Lyman-α absorption or too cold (in terms of temperature or amplitude of turbulence) to produce the required linewidths, fluxes and absorber distribution. The radio mode (hot bubbles produced in the black hole low accretion mode) in Illustris produces an unphysically large amount of volume filling hot gas at z < 2 and produces too few Lyman-α absorbers for N HI > 10 13.5 cm −2 . The updated TNG kinetic mode used for low black hole accretion is more gentle in terms of heating and alters the IGM thermal state to a lesser degree. However, TNG still does not match the COS CDD well for N HI < 10 14 cm −2 . In particular, the match of TNG to the COS data is particularly poor due to the steep slope of the CDD. Increasing the UVB of TNG by about a factor of 2 produces a better fit for column densities logN HI < 14 because the reduced equilibrium ionization fraction shifts the CDDF horizontally. The mismatch between the predicted and observed slopes in the N HI = 10 13 -10 15 cm −2 range in Figure 2 is seemingly worse than in many other simulations (Kollmeier et al. 2014;Schaye et al. 2015;Bolton et al. 2021). We will investigate the origin of the slope mismatch in a future work.
Our study demonstrates again the key role that COS HST UV spectroscopy plays in our understanding of the low-z IGM. Given the time-sensitive state of HST's remaining mission, a focus on UV absorption now is critical to provide further constraints on the physical state of IGM gas covering the Lyman-α transition at 0.4 < z < 1.7, where the Lyman-α forest statistics transition from dark matter/UVB dominated to baryonic feedback/UVB dominated (See also, Kim et al. 2020). Bolton et al. (2021) recently suggested that in order for the simulated and observed diffuse IGM Lyman-α b distribution to match, unresolved turbulent broadening at the ratio of b turb /b therm = 0.73 is needed. They found that this translates to an upper limit of v turb = 8.5 km s −1 (N HI /10 13.5 cm −2 ). Observationally, non-thermal broadening can be measured via the Doppler parameters of species with different masses in the same gas phase, e.g. using OVI, CIV and HI absorbers (Tripp et al. 2008;Werk et al. 2016;Burkhart 2021).

Missing turbulence in the IGM
Like Bolton et al. (2021), we find that both Illustris and TNG produce smaller Doppler b values than the COS observations. 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. -The various Lyman-α statistics presented in this paper for the TNG50-1, TNG100-1 and TNG300-1 simulations. The highest resolution realizations, TNG50-1, TNG100-1 and TNG300-1, include 2×2160 3 , 2×1820 3 and 2×2500 3 resolution elements, respectively. The simulation box sizes go as 35, 75, and 205 Mpc h −1 for TNG50-1, TNG100-1, and TNG300-1 respectively.
However, simulations have generally found that the Doppler parameter distribution is insensitive to resolution (see our analysis in the Appendix and also Shull et al. 2015). It may be that even higher resolution than TNG50 is required to resolve turbulent motions in the IGM and that there exists a critical resolution threshold that has not yet been reached in order to begin to see true convergence (Burkhart et al. 2015). In addition to the issue of resolution, there may be unmodeled sources of turbulence in the IGM missing from the simulations (e.g., cosmic ray induced instabilities). In order to address these questions, including a sub-grid turbulence model for the IGM gas, similar to what has been done in disk galaxies (e.g. Semenov et al. 2016) and including cosmic rays will be topics of a future exploration.

CONCLUSIONS
In this paper we have studied the low redshift Lymanα forest statistics in the Illustris and TNG cosmological simulations in comparison to the Cosmic Origins Spectrograph (COS) UV data. The TNG simulation has an identical Ultraviolet background (UVB) prescription and very similar cosmological parameters as Illustris but an entirely reworked AGN feedback prescription. We find that: • 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 HI < 10 15.5 cm −2 . Neither TNG nor Illustris can match the slope of the COS CDD.
• The b distributions in Illustris and TNG are similar, suggesting it is a poor barometer for evaluating the effects of sub-grid AGN feedback models.
• In agreement with Bolton et al. (2021), we find that both Illustris and TNG produce smaller b values than the COS observations, suggesting the need for addition sources of non-thermal broadening in the low redshift IGM.
• The amplitude and shape differences in the flux power spectrum and flux PDF between Illustris and TNG are significant (Figures 5 and 6) and can not be due solely to changes in the UVB model. We find that the flux statistics are sensitive diagnostic constraints for sub-grid AGN feedback models.
• 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 Lyman-α 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 2,184 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 Lyman-α 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 (Tillman et al. in prep).

RESOLUTION EFFECTS ON THE CDD
We investigate how resolution and simulation volume size affects the Lyman-α forest statistics presented above. Figure 7 presents the Lyman-α 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 -The Lyman-α statistics for different resolution elements in TNG300. TNG300-1 (black solid line) has 2× 2500 3 resolution elements. TNG300-2 (black dashed line) decreases the resolution elements to 2×1250 3 . TNG300-3 (dotted black line) decreases the resolution elements to 2×625 3 . The box size of TNG300 is 205 Mpc h −1 .
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 15 cm −2 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 over-densities 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 16 cm −2 ). 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 3 resolution elements.
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 −3 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 −1 and 2x 1820 3 resolution elements. TNG100-2 (black dashed lines) has the same volume size but decreases the resolution elements to 2x910 3 . TNG100-3 (dotted black lines) has the same volume as TNG100-1 and TNG100-2 but decreases the resolution elements to 2x455 3 .
The left-most 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 −2 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 Lyman-α 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 Lyman-α forest statistics. The Doppler parameter distribution seems particularly sensitive to resolution.