An Exploration of AGN and Stellar Feedback Effects in the Intergalactic Medium via the Low-redshift Lyα Forest

We explore the role of galactic feedback on the low-redshift Lyα (Lyα) forest (z ≲ 2) statistics and its potential to alter the thermal state of the intergalactic medium. Using the Cosmology and Astrophysics with Machine Learning Simulations (CAMELS) suite, we explore variations of the AGN and stellar feedback models in the IllustrisTNG and Simba subgrid models. We find that both AGN and stellar feedback in Simba play a role in setting the Lyα forest column density distribution function (CDD) and the Doppler width (b-value) distribution. The Simba AGN jet feedback mode is able to efficiently transport energy out to the diffuse IGM, causing changes in the shape and normalization of the CDD and a broadening of the b-value distribution. We find that stellar feedback plays a prominent role in regulating supermassive black hole growth and feedback, highlighting the importance of constraining stellar and AGN feedback simultaneously. In IllustrisTNG, the AGN feedback variations explored in CAMELS do not affect the Lyα forest, but varying the stellar feedback model does produce subtle changes. Our results imply that the low-z Lyα forest can be sensitive to changes in the ultraviolet background, stellar and black hole feedback, and that AGN jet feedback in particular can have a strong effect on the thermal state of the IGM.


INTRODUCTION
Over the past 50 years, the Lyα forest has been used at high redshift (z > 2) to constrain the amplitude of cosmological fluctuations, the thermal properties of the intergalactic medium (IGM), and the process of reionization (see McQuinn 2016, for a review).The majority of matter in the universe resides within the IGM, the space between galaxies, making the absorption lines that make up the forest a powerful tool for probing otherwise unobservable gas.An understanding of radiative processes reveals from the Lyα forest spectra the column density of the absorbing structure as well as its thermal properties and kinetics properties.The observed forest has been used to study small-scale cosmic structure (Croft et al. 1998;McDonald et al. 2005), the temperature of dark matter (Viel et al. 2009), gas temperature (Becker et al. 2011;Boera et al. 2014), the evolution of the metagalactic ionizing background (Fan et al. 2006;Faucher-Giguère et al. 2009;Haardt & Madau 2012), and the mapping of neutral hydrogen (Lee & White 2016;Ozbek et al. 2016).Additionally, simulations of the forest have enabled constraints on baryons, dark matter interactions, and photoionization heating and adiabatic cooling in the IGM (Cen et al. 1994;Zhang et al. 1995;Miralda-Escudé et al. 1996;Hernquist et al. 1996;Rauch et al. 1997).In more recent years, efforts with the high-z Lyα forest data have led to further constraints on ionizing background models and the termperature of the IGM (Gaikwad et al. 2017a;Oñorbe et al. 2017;Chardin et al. 2017;D'Aloisio et al. 2017;Walther et al. 2018;Puchwein et al. 2019;Faucher-Giguère 2020), properties of dark matter (Iršič et al. 2017;Rogers & Peiris 2021), the time frame and processes of the epoch of HI and HeII reionization (Zhu et al. 2021;Gaikwad et al. 2021;Villasenor et al. 2022;Wang et al. 2022;Yang et al. 2023), and the 3D mapping of cosmic structures (Lee et al. 2018;Horowitz et al. 2022;Qezlou et al. 2022a;Newman et al. 2022).
When utilizing the high redshift IGM to constrain cosmology, the role of astrophysical processes can impact predictions on the linear matter power spectrum (Pritchard et al. 2007;Wyithe & Dijkstra 2011), and the use of eBOSS data is sufficiently constraining to break most of these degeneracies (Bird et al. 2011).Additionally, localized fluctuations of the ionizing background or inhomogeneous reionization can impact the observed Lyα forest flux power spectrum leading to bias in predictions of the matter power spectrum (McQuinn 2016).These challenges require studies to be mindful of the impact that astrophysical processes, including feedback, can have on predictions that are based on the Lyα forest.However, radiation backgrounds are often approximated as uniform for large cosmological simulations so the effects of photoionization has largely been included via models for the ultraviolet background (UVB).These uniform ionizing background models include ionizing photons from galactic feedback by assuming some fraction of the photons come from feedback by AGN and supernovae.
Despite the robustness of high-z IGM studies, the low redshift (z < 2) IGM has received considerably less attention, with one early exception being the work of Davé et al. (1999).At low-z the Lyα forest is observed in the far-ultraviolet (FUV) band, requiring the use of spacebased telescopes.Only a few instruments such as the Hubble Space Telescope's (HST) Cosmic Origins Spectrograph (COS) and the Far Ultraviolet Spectroscopic Explorer (FUSE) satellite can observe in the UV band and collect Lyα data (Danforth & Shull 2005).In particular, HST COS can observe Lyα from z ≈ 0 − 0.5 and the Danforth et al. (2016;henceforth D16) catalog has built on previous surveys to produce the largest col-lection of low-z Lyα forest absorber data to date.This allows for the exploration of the low-z Lyα forest in a way previously unattainable.
Since the D16 study, many theoretical efforts have explored the low-z Lyα forest due to its relatively untapped potential.In particular, the role of feedback in the low redshift Lyα forest remains relatively unclear.Kollmeier et al. (2014) reported a disconnect between the observed and simulated Lyα forest column density distribution function (CDD), which required either an HI photoionizing rate a factor of 5 times higher than UVB models predicted at the time or the existence of nontraditional heating sources.While the factor of 5 may be partly attributed to poor resolution, it was later shown that a large increase in the UVB magnitude from Haardt & Madau (2012) (the UVB model studied in Kollmeier et al. 2014) is possible (Khaire & Srianand 2015a;Puchwein et al. 2019;Bolton et al. 2022).As for additional heating sources, several studies explore the potential impact of AGN feedback in heating the IGM, with some finding that different AGN feedback models can have a dramatic impact on the low-z Lyα forest statistics (Gurvich et al. 2017;Nasir et al. 2017;Christiansen et al. 2020;Burkhart et al. 2022;Tillman et al. 2023;Khaire et al. 2023).Fortunately, modern day simulations have put a greater emphasis on incorporating feedback mechanisms through sophisticated modeling based on observations (e.g., Vogelsberger et al. 2014a,b;Hopkins et al. 2014Hopkins et al. , 2018;;McAlpine et al. 2016;Weinberger et al. 2017;Kaviraj et al. 2017;Davé et al. 2019;Bird et al. 2022, andVogelsberger et al. 2020 for a review) making the study of feedback effects on the low-z forest possible through simulated data.
While the effects of some AGN feedback models on the statistics of the forest are degenerate with a re-scaling of the assumed UVB (e.g.see Burkhart et al. 2022;Khaire et al. 2023;Mallik et al. 2023), other AGN models had different, idiosyncratic, effects, highlighting the importance of including feedback sub-grid models within cosmological simulations (Tillman et al. 2023).The unique effects seen on the Lyα forest column density distribution function (CDD) due to AGN jet feedback (Tillman et al. 2023) are primarily attributed to the far-reaching heating effects that the jet feedback has on the largescale environment (see Tonnesen et al. 2017, for discussion of the dependence of the CDD on the large-scale environment).These results emphasize the need for an improved knowledge of the physical conditions of the IGM and how they vary on large scales.This not only includes a better understanding of UV ionizing photon sources, but also any astrophysical effects such as feed-back that might affect the HI within the IGM via heating or gas transport.
The Lyα forest Doppler width, or b-value, distribution can provide information about the thermal and turbulent state of absorbers.However, matching the simulated b-value distribution with the observed distribution has proved more difficult than matching the CDD.The simulations consistently under-predict the number of high b absorbers as compared to observations.Several studies have posed the idea of introducing nontraditional heating sources, such as galactic feedback, or extra line of sight turbulent velocity components in simulations to alleviate the tension (Viel et al. 2017;Gaikwad et al. 2017b;Bolton et al. 2022;Burkhart et al. 2022).Bolton et al. (2022) and Viel et al. (2017) both discuss the difficulties of using AGN feedback as a solution to correcting the b-value distribution.The AGN feedback models explored in these studies (which were heavily based on the Illustris and IllustrisTNG models) failed to reproduce the observed b-value distribution, however, whether this statistic might be used to constrain other feedback models is unclear.It is likely that the combination of unresolved turbulence and heating from galactic feedback that does not over-ionize HI in the IGM is required to reproduce the observed b-value distribution.
To further understand the role of AGN feedback in the low-z Lyα CDD and b-value distributions, a systematic exploration of multiple different AGN feedback models is required.The CAMELS simulations provide an excellent opportunity for this type of analysis (Villaescusa-Navarro et al. 2021).By varying multiple parameters for both AGN and stellar feedback we can not only explore the effects of AGN feedback models in different simulations, but we can also explore the variation of feedback parameters within each model.
In this study, we explore variations on the Simba and IllutrisTNG stellar and AGN feedback models using the CAMELS project simulations.We analyze both the Lyα forest CDD and b-value distributions within the context of these simulations.In Section 2 we describe the CAMELS project and provide a brief description of the relevant numerical models of the simulation suites we utilize in this study (Section 2.1).We also discuss how we generate our synthetic spectra from the simulations (Section 2.2), how we calculate the Lyα forest CDD and b-value distribution (Section 2.3), and the different supermassive black hole (SMBH) statistics that we analyze (Section 2.4).In Section 3 we present our results, followed by a discussion in Section 4. Finally we conclude and discuss implications and future work in Section 5.

The CAMELS Simulations
The CAMELS1 project consists of thousands of Nbody and (magneto-)hydrodynamical simulations run with the AREPO, GIZMO, and GADGET codes (Springel 2005;Villaescusa-Navarro et al. 2021).The CAMELS simulations are run with the same sub-grid models as pre-existing simulations while varying different cosmological and astrophysical parameters.Currently the CAMELS project has made simulations run with the IllustrisTNG (Weinberger et al. 2017;Pillepich et al. 2018), Simba (Davé et al. 2019), and Astrid (Bird et al. 2022) sub-grid models publicly available (Villaescusa-Navarro et al. 2023;Ni et al. 2023).However, in this work we will focus on IllustrisTNG and Simba as the Astrid AGN feedback closely resembles that of IllustrisTNG.
For the IllustrisTNG and Simba suites, 6 parameters are varied for each sub-grid model.The 2 cosmological parameters explored are Ω m , the matter density, and σ 8 , the variance of the linear field on 8 Mpc/h scales at z = 0.The 4 astrophysical parameters explored vary the different feedback models with 2 parameters assigned to stellar feedback and 2 for AGN feedback.The exact physical processes that the feedback parameters control vary between the different simulation sub-grid models as the feedback models differ dramatically between them.The extensive range of cosmological and astrophysical parameters explored makes the CAMELS simulations a unique data set to constrain astrophysical and cosmological models.
Each CAMELS simulation has 256 3 gas resolution elements in a periodic comoving volume with a side length of 25 Mpc/h.The mass resolution of the CAMELS simulations is similar to the original IllustrisTNG300-1 simulation.Each simulation also utilizes cosmological parameters Ω b = 0.049, h = 0.6711, n s = 0.9624, M ν = 0.0 eV, w = −1, and Ω K = 0.In this study we do not explore simulations that vary Ω m and σ 8 , and therefore these values are set to the fiducial Ω m = 0.3 and σ 8 = 0.8 in every case.For each suite (i.e.subgrid model for IllustrisTNG or Simba) there are 3 main "sets" of simulations: the 1P set (61 simulations), the CV set (27 simulations), and the LH set (1000 simulations).1P stands for 1-Parameter and the simulations in this set vary one parameter at a time.CV stands for Cosmic Variance and these simulations vary only the initial random seed.LH stands for Latin Hypercube and these simulations vary all 6 cosmological and astrophys- Wind speed (Eq.2) IllustrisTNG Energy per unit BHAR (Eq.9) Burstiness (Eq.8) Energy per unit SFR (Eq.5) Wind speed (Eq.6) ical parameters and the random seed simultaneously.
The LH set is ideal for machine learning applications.
A short summary of the suite and set terminology utilized in this study as well as a brief description of each feedback parameter is presented in Table 1.This study primarily focuses on the 1P set simulations that vary the stellar and AGN feedback parameters individually over the range shown in Table 1 while all other feedback parameters are held at their fiducial value of 1.Additionally, we only explore the Simba and IllustrisTNG suites as the Astrid feedback models (especially the AGN feedback model) resemble those of Il-lustrisTNG.The feedback parameters (labeled A AGN1 , A AGN2 , A SN1 , and A SN2 ) correspond to different aspects of the feedback models in IllustrisTNG vs. Simba as the models are dramatically different.The parameters A AGN1 and A AGN2 vary aspects of the AGN feedback models while A SN1 and A SN2 vary aspects of the stellar feedback models.Mathematical descriptions of the feedback parameters are provided in the following subsections along with an outline of the different feedback models.

SIMBA
The CAMELS Simba suite follows the same subgrid model as the original Simba simulations and is fully described in Davé et al. (2019).Simba utilizes GIZMO (Hopkins 2015), an N-body hydrodynamics code, in its Meshless Finite Mass hydrodynamics mode.Radiative cooling and photoionizing heating are modeled with Grackle3.0(Smith et al. 2017) which includes non-equilibrium evolution of primordial elements, a partially-uniform ionizing background (Haardt & Madau 2012), hydrogen self-shielding (Rahmati et al. 2013), and metal cooling.growth by accretion of metals, and destruction by thermal sputtering and SNe are also tracked (Li et al. 2019).

SIMBA STELLAR FEEDBACK
The stellar feedback drives galactic winds kinetically with hydrodynamically-decoupled two-phase metalenriched winds.30% of those wind particles are heated

Low
Fiducial High based on the SNe energy and wind kinetic energy.The gas elements are ejected stochastically depending on the mass loading factor η * ≡ Ṁwind /SFR (Star Formation Rate) and the wind velocity v w .For SIMBA, A SN1 controls the normalization of η * such that increases in A SN1 results in a proportional increase in the mass ejected per SFR: where M 0 = 5.2 × 10 9 M ⊙ .Both the mass loading factor and the SNe wind speed prescriptions are based on the FIRE zoom-in simulations (Hopkins et al. 2014;Muratov et al. 2015;Anglés-Alcázar et al. 2017b).The SN wind speed at which mass is ejected is proportional to A SN2 plus some velocity corresponding to the potential difference between position of the SN and 0.25 times the virial radius (R vir ).This leaves A SN2 to control the normalization of the SN wind speed as follows (2) In summary, the parameter A SN1 changes the mass loading factor of SN while A SN2 varies the SN wind speed.These parameters are varied within a range of A SN1 ∈ [0.25, 4] and A SN2 ∈ [0.5, 2], and have a fiducial value of 1 to match the original Simba simulations.

SIMBA SUPERMASSIVE BLACK HOLE FEEDBACK
In the Simba suite, the black hole module is based on the gravitational torque and Bondi accretion models combined with the kinetic feedback sub-grid model in GIZMO (Anglés-Alcázar et al. 2017a).SMBHs are seeded with mass M seed = 10 4 M ⊙ h −1 in galaxies with M * ≥ 10 9.5 M ⊙ .SMBHs are re-positioned to the potential minimum of their host group if it is within 4 × R 0 .R 0 is the size of the BH kernel which encloses the nearest 256 gas elements and any BHs within R 0 of each other are merged so long as their relative velocity is less than 3 times the mutual escape velocity.
SMBH growth occurs via a two-phase model with cold gas accreted at a rate controlled through the transport of angular momentum by gravitational torques from stars (Hopkins & Quataert 2011;Anglés-Alcázar et al. 2017a), and hot gas is accreted via spherical Bondi accretion (Bondi 1952).For both modes, a radiative efficiency of 0.1 is assumed for accretion.The cold-gas torquebased accretion is capped at three times the Eddington limit and the hot gas Bondi accretion mode has a strict maximum of the Eddington limit.
AGN feedback follows a two-mode model: a high Eddington ratio (η ≡ ṀBH / ṀEdd where ṀEdd = L Edd /ϵ r c 2 ) radiative mode with high mass loading outflows, and a jet mode with faster outflows but lower mass-loading at low Eddington ratios.For all AGN feedback, gas is ejected in a bipolar fashion parallel and anti-parallel to the angular momentum vector of the gas within the SMBH kernel.The total momentum flux of AGN feedback follows: where L bol = ϵ r ṀBH c 2 is the bolometic luminosity, ϵ r = 0.1 is the radiative efficiency, and c is the speed of light.
In the radiative feedback mode, the outflow velocity scales with M BH like v rad = 500 + 500(log 10 (M BH ) − 6)/3 km s −1 .When the SMBH has a mass M BH < 10 7.5 M ⊙ or a high Eddington ratio of η > 0.2 it produces feedback in the radiative mode.Otherwise, the feedback is produced following the jet mode prescription.
The jet mode feedback gains an additional velocity kick of v jet = 7000 × min[1, log 10 (0.2/η)] km s −1 such that full jet speeds are reached when η < 0.02.Based on both the radiative and jet feedback mode prescriptions, the outflow of any SMBH producing feedback follows: otherwise.
(4) In this formulation A AGN2 controls the maximum jet velocity achieved in the jet feedback mode.At the maximum jet speed and when M gas /M * < 0.2 in the galaxy, SMBHs will also inject energy into the gas immediately surrounding it in the form of X-rays.The particles ejected in the jet mode are hydrodynamically de-coupled from the rest of the gas in the box for a time that scales with the Hubble time at the moment of ejection.This results in the jets traveling a distance of up to ∼ 10 kpc and losing a maximum of 1000 km/s in speed due to gravitational effects before re-coupling to the surroundings.The jet particles are heated to the virial temperature of the host halo (so long as v jet > 2000 km/s which is almost always the case for jet mode AGN at z < 2.0) and are ejected in a bipolar fashion parallel to the angular momentum vector of the disk with a zero degree opening angle.
The total momentum flux formulation in Equation 3, where Ṗout is a set value which can be varied by A AGN1 , applies to both the radiative and jet modes in Simba.However, the way it scales with the outflow velocity (v out ) means that significantly less mass is ejected via AGN feedback in the jet mode than the radiative mode.Changing A AGN1 will change the amount of ejected mass in both modes but will have a greater impact on the radiative mode.Changing A AGN2 will also change the amount of mass ejected but only for the jet mode.For example, increasing only A AGN2 will cause less mass to be ejected at higher velocities in the jet mode since the momentum flux is set, but it would not affect the amount of mass ejected in the radiative mode.
In summary, in the Simba suite A AGN1 controls the momentum flux and A AGN2 controls the additional jet velocity kick.With how v rad is calculated, it is not affected by A AGN1 but the amount of mass ejected ( Ṁout ) is.Conversely, the jet speed (v jet ) is directly proportional to A AGN2 allowing the maximum jet boost to reach up to 14,000 km s −1 in the most extreme case (A AGN2 =2).The AGN feedback parameters are varied with the ranges A AGN1 ∈ [0.25, 4] and A AGN2 ∈ [0.5, 2.0], with fiducial values of 1.
For additional information, see the CAMELS documentation (Villaescusa-Navarro et al. 2021) and the original Simba documentation (Davé et al. 2019) and references therein.

IllustrisTNG
The IllustrisTNG (TNG) suite consists of magnetohydrodynamic cosmological simulations that utilize the same sub-grid models as the original TNG simulation set (Pillepich et al. 2018).The simulations are run utilizing the AREPO code (Springel 2010;Weinberger et al. 2020) and gravitational interactions evolve through the TreePM algorithm (Springel et al. 2005).Radiative cooling from hydrogen and helium is implemented using the Katz et al. (1996) network and includes cooling via line free-free emission, inverse Compton, and line cooling.Metals and metal-line cooling are included as described in (Vogelsberger et al. 2012(Vogelsberger et al. , 2013)).TNG assumes ionization equilibrium with a Faucher-Giguère et al. ( 2009) UV background and accounts for on-the-fly hydrogen column density shielding from the radiation background (Rahmati et al. 2013).The star formation and interstellar medium sub-grid model is from Springel & Hernquist (2003).

TNG STELLAR FEEDBACK
TNG models stellar feedback from SNe and tracks chemical enrichment from 9 elements (H, He, C, N, O, Ne, Mg, Si, Fe).Enrichment is modeled from Type Ia SNe, Type II SNe, the Asymptotic Giant Branch, and neutron star-neutron star mergers.In TNG, stellar feedback from SNe drives kinetically implemented (with a thermal energy sub-component) galactic winds that are modeled as temporarily hydrodynamically decoupled particles.These winds are stochastically and isotropically ejected from star-forming gas.The total energy injected per unit star-formation rate and stellar feedback driven galactic outflow speed are modified by the parameters A SN1 and A SN2 , respectively, via the formulations below.The energy per unit SFR is: where Z is the gas metallicity.The wind speed of galactic outflows is:  is then η w ≡ Ṁwind /SFR = 1.8v −2 w e w ∝ A SN1 /A 2 SN2 .This results in both A SN1 and A SN2 playing a role in calculating the mass loading factor for SNe feedback.
In summary, the energy injection is proportional to A SN1 with a fiducial value of 1 and the range of variation A SN1 ∈ [0.25, 4].The SN wind speed is proportional to A SN2 with a fiducial value of 1 and a range of variation of A SN2 ∈ [0.5, 2].The fiducial values correspond to the original TNG runs.

TNG SUPERMASSIVE BLACK HOLE FEEDBACK
In the TNG suite, SMBH particles are seeded in halos with mass M halo > 5 × 10 10 M ⊙ h −1 with mass M seed = 8 × 10 5 M ⊙ h −1 .BH accretion in TNG uses the Bondi (1952) accretion prescription capped at the Eddington limit.A radiative efficiency of 0.2 is assumed for all BH growth.
The SMBH feedback in TNG includes three modes: thermal, kinetic, and radiative.The radiative feedback mode is always on, adding the radiation flux of the SMBH to the cosmic ionizing background heating the gas in and around the host halo.In addition to the radiative mode, either the kinetic or thermal mode is active.Which mode the SMBH produces feedback in depends on the Eddington ratio.The efficiency fraction (fraction of accreted mass converted into energy) for the thermal mode is a constant 0.02.The efficiency fraction in the kinetic mode is calculated as min[0.2,ρ/(0.05ρSFthresh ]), where ρ is the density of the gas around the SMBH and ρ SFthresh is the star formation threshold density.
In both the kinetic and thermal feedback modes, energy is injected in the 'feedback sphere' of the SMBH.The feedback sphere has a size that scales with reso- baryon .The size of the sphere is roughly constant within a simulation but varies slightly with the particles neighboring the SMBH (Weinberger et al. 2017;Pillepich et al. 2018).Unlike the stellar feedback, there is no de-coupling of the cells within the feedback sphere from hydrodynamical forces or radiative cooling.The transition from thermal to kinetic feedback mode happens at Eddington ratios (η) lower than some threshold value (χ) with (7) Due to this definition, the transition from high-accretion thermal feedback mode to low-accretion kinetic feedback mode occurs at M BH ∼ 10 8 M ⊙ .In CAMELS, the highaccretion rate thermal mode is set as an injection of thermal energy into the defined sphere around the SMBH.The efficiency fraction of the mass-to-energy conversion for energy injection is Ėhigh = 0.02× ṀBH c 2 where ṀBH is the instantaneous SMBH accretion rate and c is the speed of light.
The CAMELS simulations explores variations of the low-accretion rate kinetic mode.For the kinetic mode, energy is accumulated until a certain threshold is reached (since the last event) after which the energy is injected into the feedback sphere in a random direction (averaging over multiple events the injections become isotropic).The energy threshold at which injection occurs is: where σ2 DM is a one-dimensional dark matter velocity dispersion around the SMBH and m enc is the gas mass in the feedback sphere.The amount of energy produced per accretion event is proportional to both A AGN1 and the gas density around the SMBH (ρ) with: where ρ SFthresh is the density threshold for starformation.
Given these relationships, when A AGN1 is doubled so is the efficiency at which mass is converted into energy per accretion event.When A AGN2 is doubled so is the energy threshold in which injection happens, therefore larger values for A AGN2 results in less frequent but stronger AGN feedback events.In this sense, for the TNG SMBH feedback model in CAMELS, A AGN1 controls the amount of energy produced per SMBH accretion event while A AGN2 controls the burstiness and strength of the SMBH feedback.
The parameters controlling AGN feedback in the CAMELS TNG suite are intertwined such that increasing A AGN1 alone does not increase the strength of any individual AGN feedback injection but only increases the frequency at which these events occurs.If A AGN2 is also increased, the feedback events grow stronger but at the expense of the frequency at which the events occur.The range explored for these parameters are A AGN1 ∈ [0.25, 4] and A AGN2 ∈ [0.5, 2.0] with a fiducial value of 1 corresponding to the original TNG runs.
For additional information, see the CAMELS documentation (Villaescusa-Navarro et al. 2021) and the original TNG documentation (Pillepich et al. 2018) and references therein.

Synthetic Spectra
We generate synthetic spectra from the CAMELS simulations using the publicly available fake-spectra 2 code outlined in Bird et al. (2015); Bird (2017) with MPI support from Qezlou et al. (2022b).From simulation snapshots the code generates and analyzes mock spectra.The fake-spectra package is fast, parallel, and is written in C++ and Python 3 with the user interface being Python-based.
Column densities (CDs, N HI ) are computed by interpolating the neutral hydrogen mass in each gas element to the sightline using an SPH (smoothed particle hydrodynamics) kernel.The method used is based on what is appropriate for the corresponding simulation; for the CAMELS-TNG simulations a tophat (or uniform) kernel is used while for CAMELS-Simba a cubic spline kernel is used.The CDs have units of neutral hydrogen atoms (HI) per cm −2 .
We generate 5,000 sightlines randomly placed in each simulation box; a number found to be sufficient for avoiding variations due to sampling (Tillman et al. 2023).We do not add noise to the spectra generated from the simulation box.Adding noise to the spectra will not affect the results of this study as it would largely alter the lowest CD values, where observational errors are large, and the b-values predicted, which we do not compare to observations herein.In this study, we focus primarily on CDs of 10 12 < N HI < 10 15 cm −2 .Below 10 12 cm −2 the lines are too faint to detect and characterize at current observational sensitivities, while above 10 15 cm −2 lines are saturated and it becomes difficult to accurately determine CDs from flux.Additionally, absorbers at the highest CDs are rarer, which makes them difficult to study both observationally, and in the smallbox (25 Mpc h −1 ) 3 CAMELS simulations.Tillman et al. (2023) found that at least for TNG, the CAMELS boxes produce a converged CDD as compared to the original TNG100-1 simulation, within observational error bars.

Lyα Statistics
The CCD (f (N HI )) is defined as where F (N HI ) is the fraction of absorbers with column densities in the range [N HI , N HI + ∆N HI ], and ∆z is the redshift distance of the sightline.The CDD describes the number of absorbers within a logarithmic column density bin width and redshift distance.
To calculate the CDD we utilize a direct integration method as described in Tillman et al. (2023).The particles in 525 kpc/h slices are amalgamated into absorbers that are then used to calculate column densities for the CDD.As discussed in Tillman et al. (2023) this method and the defined size of the absorbers results in a well converged CDD.Variations in the chosen absorber size (reasonable sizes ranging from 300 kpc/h to 800 kpc/h) results in a less than ∼ 15% effect on the CDD at CD values lower than 10 12.5 cm −2 , a less than 1% effect at higher CD values, and all differences are well within 1σ of observational error bars.It has been previously found that the direct integration method as compared to Voigt profile fitting provides a CDD that is converged within 1 sigma with respect to observational errors (Gurvich et al. 2017).At the lowest column densities (∼ 10 12 cm −2 ) the CDD produced via Voigt profile fitting diverges from the direct integration significantly likely due to the difficulty in fitting those absorbers with an automated procedure.However, the divergence is still within 1 sigma of the observational error bars where observational data exists and does not influence the main results of this paper.
To calculate the b-value (Doppler width) distribution, Voigt profile fitting is necessary.To conduct our fits we utilize the Voigt fitting algorithm included in the fakespectra package3 .The algorithm is closely based on that of AUTOVP (Davé et al. 1997).Peaks are found, fit, and iteratively removed after which the peaks are refit to the spectrum all together.The algorithm starts with the largest peaks and continues until adding another peak to the fit no longer improves the fit by some chosen significance value.Another stopping condition is when the largest peak remaining is less than 10 −4 ×max[1, max(τ )] where τ is optical depth.These conditions are in place since, as the fit continues, smaller and smaller peaks are more likely to become fitting errors rather than an actual absorber.It is due to this that automatic Voigt profile fitting is so difficult, and this is why the lower end of the CDD diverges when using the Voigt fitting method.
The original Voigt profile fitting algorithm in the fakespectra package minimizes the squared difference between the fit and the original spectrum using the Nelder-Mead algorithm.For our work we use a modified version of the Voigt profile fitter that minimizes the earth mover's distance using a simple bounds limited-memory Broyden-Fletcher-Goldfarb-Shanno algorithm.We found that this method had an easier time fitting smaller peaks in the spectra and ran into fewer errors during computation, but this choice did not change the overall results of our study.
The b-value is affected by both the temperature and turbulence of an absorber.This manifests mathematically as: where λ 0 is the wavelength of the Lyα transition, m H is the mass of a hydrogen atom, k is the Boltzman constant, T is the temperature, and σ v is the root meansquare turbulent velocities of the absorber.The b-value can also have components from peculiar and Hubble velocities.At least at high redshift (z ∼ 2 to 4) the b-value of absorbers tends to be dominated by the Hubble flow (corresponding to the physical width of the absorber along the line of sight) but is also affected by thermal broadening and pressure support along the Jean's length, but the narrowest absorption features are dominated by thermal broadening (Bryan & Machacek 2000;Peeples et al. 2010a).At lower redshifts, such as those explore in this work, the Hubble flow becomes less important in absorber broadening.
The formulation above also puts a lower limit on bvalues we examine in this work.Assuming an absorber is non-turbulent and recognizing that the cooling floor of the simulations is T min = 10 4 K, we should assume a minimum b-value of b min ≈ 13 km s −1 .We limit the upper bounds of our b-value analysis to 100 km s −1 as this is where observational data is available.However, this bound does not affect the results of our study because there are so few absorbers at such high values.
We do not compare to observational values or the original simulations since the b-value distribution is more sensitive to mass resolution than the CDD (see Appendix in Burkhart et al. 2022).Lower numerical resolution causes a nonphysical broadening of the b-value distribution shifting to higher b-values, making the comparison to observations challenging (artificial broadening due to numerical resolution limits as seen in Peeples et al. 2010a).Due to this, the b-value distributions generated from the CAMELS simulations are most robustly interpreted when comparing between simulations using the same numerical resolution, as done in this paper.See Appendix A for additional discussion on Lyα forest statistics convergence.
We explored introducing a Gaussian line spread function (LSF) with a full width half maximum of 6.5 to our data to determine how noise due to the HST-COS instrumentation might affect our results.While this is not the exact LSF utilized when fitting COS data, it represents a good approximation and should be sufficient in determining the general instrumentation effects on this study.The LSF broadens the absorption features, which leads to larger b-values and smaller CDs.As previously mentioned, noise in the spectra largely affects lower b and CD values.The effects of introducing this LSF on the b-value and CD distributions are small and our main results remain unchanged therefore the results presented herein include no instrumentation or random noise.
For both the CDD and the b-value distributions, we will explore how these statistics vary for the various feedback parameters explored in CAMELS.To help visualize the effects on the various statistics we also generate projection plots for the various CAMELS simulations of both absorber CDs and average temperatures.The CD projections help interpret the CDD while the temperature projections can help interpret both the CDD and the b-value distribution.These projections can be seen in Figures 1, 2 and 3.

Black Hole Statistics
We also analyze certain BH statistics to determine the extent to which stellar feedback can suppress SMBH growth and thus AGN feedback, and to see how efficient AGN feedback is at self-suppressing SMBH growth.We analyze the accretion rate density for SMBHs in the CAMELS simulations which is the total accretion rate from all SMBHs in the box divided by the box volume ( ρBH = i ṀBH,i /25 Mpc 3 for all SMBHs i).The energy emitted in feedback scales with accretion rate, thus ρBH gives a sense of the total energy in the box due to SMBH feedback.
For CAMELS-TNG, the thermal AGN mode produces feedback as thermal energy injected into the area immediately surrounding the BH whereas the kinetic mode drives matter and energy out to larger distances.Considering this, we expect the kinetic mode to be more likely to affect the Lyα forest as opposed to the thermal mode in TNG.However previous work has shown, by completely removing the kinetic feedback, that this feedback mode does not affect the CDD in any meaningful way (Tillman et al. 2023).In this work, we instead find the radiative AGN feedback mode to be the most impactful on the forest statistics in TNG.The contribution from the radiative mode scales with the amount of matter being accreted, therefore, for TNG, the contribution to the UVB from AGN feedback scales with ρBH .
For Simba, the jet mode feedback velocity boost depends on η such that a lower η produces faster jet speeds up to a maximum jet speed boost of A AGN2 ×7000 km s −1 for η ≤ 0.02.As seen in Christiansen et al. (2020); Tillman et al. (2023), the jet mode in Simba can have a dramatic effect on various IGM statistics particularly due to the ability for jets to reach far into the IGM.Analyzing the AGN jet mode through the SMBH accretion rate will aid in interpreting the effects of varying the feedback parameters.Lower values for ṀBH in the jet mode correspond to higher jet speeds but less energetic events overall.
For both TNG and Simba, we supplement the accretion rate statistics with the number of SMBHs in the different feedback modes.We will analyze these BH statistics for variations of the CAMELS-TNG and CAMELS-Simba feedback parameters A SN1 , A SN2 , A AGN1 , and A AGN2 .

RESULTS
In this study we explore statistics from the Lyα forest and from the SMBHs themselves to determine what influences the feedback has on the neutral hydrogen in the IGM and why those effects manifest.We analyze the mass-weighted mean temperature and CD projections in Section 3.1.We also explore the CDD and bvalue posterior distribution function in Sections 3.2 and 3.3 respectively.Finally, we look at SMBH statistics regarding the Eddington ratio, accretion rate, number of SMBHs in the simulation, and the number of SMBHs accreting in each feedback mode in Section 3.4.

Temperature and Column Density Projections
The projection plots in Figures 1, 2 and 3 clearly show how varying the different feedback modes in TNG and Simba affects the IGM.The projections show a 525 kpc/h slice of the simulations corresponding with the typical size of an absorber found in the Lyα forest as defined in this study.For most definitions, used in previous works, the absorber length is found to be less than 500 kpc/h (Tonnesen et al. 2017) but for temperatures as low as 10 4 K an absorber length of 800 kpc/h could be expected (Peeples et al. 2010b).However, for the analysis methods used herein, it was found that variations of the absorber length within these ranges had less than a 15% impact on the resulting CDD (Tillman et al. 2023).The projection plots are good visual indicators of the overall effect of varying the stellar and AGN feedback parameters in CAMELS.
Figures 1 and 2 show mass-weighted mean temperature projections and column density projections for the Simba suite when varying the feedback parameters.Figure 1 shows projections for variations in the AGN feedback parameters (A AGN1 and A AGN2 respectively), and Figure 2 shows variations in the stellar feedback parameters (A SN1 and A SN2 respectively).For both figures, the top two rows display temperature projections and the bottom two rows are column density projections.The left column is for the lowest value of the parameter explored, the middle column indicates fiducial results, and the right column displays the highest value, as indicated at the top of each figure.
Varying AGN momentum flux (A AGN1 ) and AGN jet speed (A AGN2 ) have clear consequences for both the temperature in the box and the amount of neutral hydrogen in the IGM.With increases in either AGN parameter value, we see higher temperatures and less neutral hydrogen populating the IGM.However, much of the change in HI abundances is seen in CDs outside of the range of interest for the Lyα forest.Changes in A AGN2 demonstrate these effects more dramatically than does A AGN1 .
Varying the SN mass loading factor (A SN1 ) and the SN wind speed (A SN2 ) in Simba also has a clear impact on both the temperature and amount of neutral hydrogen in the box.For increases in the SN mass loading factor we see increases in temperature and decreases in the amount of neutral hydrogen.The extent of this effect appears similar to that of the AGN feedback parameters.However, increasing the SN wind speed has the opposite effect, with higher wind speed decreasing the temperature and increasing the fraction of neutral hydrogen.
The SN wind speed has the most dramatic effect on the IGM out of the four feedback parameters explored in CAMELS-Simba.Stellar feedback as implemented in these simulations does not, by itself, cause IGM scale effects.If heating from SN directly influenced the forest, then increases in SN wind speed should see a decrease in Lyα absorbers, however the opposite is true.As we will explore in more detail in this section and the discussion (Section 4), the stellar feedback's effect on the IGM comes indirectly through its impact on SMBH growth and AGN feedback.
The complex interaction between stellar feedback and SMBH growth in simulations has been explored in many studies (Booth & Schaye 2013;Dubois et al. 2015;Habouzit et al. 2017;Anglés-Alcázar et al. 2017;van Daalen et al. 2019;Lapiner et al. 2021;Tillman et al. 2022;Byrne et al. 2023;Delgado et al. 2023).Star formation is able to regulate the growth of galaxies and their central SMBHs, and stronger stellar feedback can both reduce the need for AGN feedback to regulate star formation as well as suppress SMBH growth.Weaker AGN feedback will then diminish the overall impact on the IGM.We explore SMBH statistics for variations in the feedback parameters to help illustrate this interplay between AGN and stellar feedback.
For the TNG suite we only display temperature projections because the column density projections exhibit no visible changes when varying feedback parameter values.Figure 3 resembles closely that of Figures 1  and 2, but shows only the temperature projections for CAMELS-TNG.AGN feedback parameter variations are on the left and stellar feedback variations are on the right.The temperature of the IGM in TNG is overall much lower than that in Simba.Increases in the AGN feedback parameters show minimal increases in temperature and those changes appear largely confined to host halos and the surrounding area.Decreases in the stellar feedback parameters show more substantial temperature changes but the effects still reside within ∼1-2 Mpc of host halos.Similar to the Simba suite, the stellar feedback parameter variations in TNG also high-light the importance of AGN-stellar feedback interplay.However, the impacts are not as pronounced as they are for the Simba suite as the TNG AGN feedback model has minimal effect on the IGM.

Column Density Distribution Functions
Next, we examine the CDD to quantify the impacts of the parameter variations seen in the previous section.The two leftmost plots of Figure 4 show the CDD for different values of A AGN1 and A AGN2 .Only the Simba suite is displayed in these plots as the TNG results for AGN feedback showed no discernible differences, as evident in the projection plots.
For Simba, as implied in the projection plots, decreases in the AGN feedback parameters show increases in the amount of neutral hydrogen in the Lyα forest.However, increases in A AGN1 , relative to the fiducial value, show no observable effect on the CDD (with respect to the observational error bars) while increases in A AGN2 result in less neutral hydrogen overall.Increases in the A AGN1 parameter affect the temperature and HI fraction of gas not associated with the Lyα forest since differences are clearly visible in the projection plots but not in the CDDs.The gas affected is too hot and at too low column densities to be associated with the forest.
The right four plots of Figure 4 show the same results for the stellar feedback parameters.For the Simba suite, decreasing A SN1 shows an increase in neutral hydrogen absorbers as implied by the results in Figure 2.However, increasing A SN2 shows an increase in HI absorbers.For the remainder of the parameters (A SN2 for Simba, and A SN1 and A SN2 for TNG) larger values result in more absorption while smaller values result in less.The effect is the most dramatic for the SN wind speed in Simba and minimal for the stellar feedback parameters in TNG.However, both of these results imply the suppression of SMBH growth by stellar feedback.We discuss the relationship between stellar and AGN feedback further in Section 3.4.
Figure 9 shows the redshift evolution of the TNG and Simba suites from z = 2.0 to z = 0.In these plots, the TNG CDD is corrected to utilize the Haardt & Madau (2012) UV ionizing background to match the one used in Simba.Removing the difference in the assumed UVB model allows for easier comparison of AGN and stellar feedback effects between the two simulations.In TNG, the AGN radiative feedback mode adds to the assumed ionizing background which effectively results in a slightly stronger UVB than the assumed Faucher-Giguère et al. ( 2009 the magnitude of the difference made by the radiative feedback mode in TNG is well visualized at z = 2.0.If we instead rescale both of the simulations at z = 2.0 to have the same Lyα mean flux, the predicted CDDs of TNG and Simba lie on top of one another.At z = 2 the shape of the CDD is nearly identical between the two simulations, after which they diverge around z ∼ 0.5.

The Doppler Width (b-value)
Similar to the CDD analysis, we explore the b-value distributions for variations in the CAMELS feedback parameters for both the TNG and Simba suites.Figure 5 shows how the b-value probability density function (PDF) varies with the feedback parameters.Overall Simba has a broader distribution of b-values due to having a hotter IGM than TNG.However, the peak of the distribution is around 20 km/s for both simulations.
The left two plots of Figure 5 show the effects of AGN feedback in Simba on the b-value distribution.Higher bvalues are observed for stronger AGN feedback, with the AGN jet speed (A AGN2 ) demonstrating a stronger effect than the mass loading factor (A AGN1 ).Increases in the AGN feedback parameters, in Simba, show an overall increase in the IGM temperature (seen in Figure 1) so increased b-values, seen in Figure 5, are expected.As with the CDD results, we do not include the b-value distributions for variation in the TNG AGN feedback parameters as no discernible difference is seen.
The right four plots of Figure 5 show variations in the stellar feedback parameters for both Simba and TNG.For Simba, decreases in the stellar feedback parameters show a shift to higher b-values with the SN wind speed (A SN2 ) having a more dramatic effect than the SN mass loading factor (A SN1 ).For TNG, variations in the stellar feedback parameters have only a small effect with marginal shifts to higher b-values for lower SN energy per unit SFR (A SN1 ) and higher SN wind speed (A SN2 ) but these differences are subtle and would not be measurable in observations.The results from Figure 5 parallel the effects seen on the CDD, but perhaps to a smaller degree than was seen for the CDD.This can be explained by the fact that both the CDD and the b-value distribution have a temperature scaling, although the b-value tends to have a stronger temperature correlation.These temperature effects are due to the AGN jets in Simba and the radiative mode in TNG.

Supermassive Black Hole Statistics
Now that we have demonstrated the impact of varying CAMELS parameters on the Lyα forest statistics, we explore its cause.When varying the different feedback parameters available in the CAMELS simulation suite, Simba's feedback exhibited clear effects on the Lyα forest CDD and the b-value PDF, while TNG's feedback produced no effects that could be observationally measured, at least for the parameter variations explored.These results are consistent with previous work comparing the AGN feedback models in Simba and TNG and their effects on the Lyα forest CDD (Tillman et al. 2023).
However, the stellar feedback parameters show a subtle effect on the CDD and b-value distribution for TNG.This may be initially surprising, as we do not expect the effects of stellar feedback to directly manifest in low CD Lyα forest absorbers, which can be greater than a Mpc away from any SNe.The scale of a single SN (which individually drives winds on the scale of pcs) is not relevant when thinking about the larger scale impacts of galactic winds contributed by multiple SNe.However, we find in this work, for both Simba and TNG, that changes in Lyα statistics due to varying stellar feedback are not due to galactic wind impact but rather due to stellar feedback suppressing SMBH growth.Therefore, these results imply that the AGN feedback in TNG has at least a small impact on the Lyα forest and that the parameters chosen to be varied for TNG's AGN feedback are not representative of said impact.To explore the interplay between stellar and AGN feedback, we analyze SMBH statistics over the redshift range of z = 2 to 0. The results are displayed in Figures 6 and 7. Exploring how SMBH accretion and seeding in the simulations vary with stellar feedback will help illustrate how SMBHs are affected.
The top row of Figure 6 shows the accretion rate density ρBH for the Simba AGN jet feedback mode.Note that all of these SMBHs have masses M BH > 10 7.5 M ⊙ and Eddington ratios η ≤ 0.2 in order to be in the jet mode.The second row shows the number of SMBHs producing jet mode feedback.Each column of plots represents a different feedback parameter that is being varied (labeled at the top), with the red lines corresponding to the highest values for the feedback parameters, the black lines corresponding to the fiducial results, and the blue lines correspond to the lowest values.
Figure 7 is similar to Figure 6 in that it shows ρBH for the TNG thermal feedback mode in the top two row and the number of SMBHs in the bottom row.The radiative mode is always on and scales with the overall accretion rate but the kinetic mode is radiatively inef-ficient so we are mostly interested in the thermal mode (Weinberger et al. 2017).Variations in ρBH due to the AGN parameter variations are not shown since no effect is seen.
Figures 7 and 6 clearly show that the stellar feedback parameters in both simulation suites have an impact on SMBH growth in the simulations.We discuss the interplay of stellar and AGN feedback for each simulation suite individually.

SIMBA AGN FEEDBACK PARAMETERS
By focusing on the first row of Figure 6, we find that when increasing A AGN1 , the BHs in the box have overall larger AGN jet velocities (via Equation 4 i.e. smaller ṀBH means smaller η) starting at higher redshifts.However, these BHs are less massive overall and are accreting less mass (similar to what was found in Anglés-Alcázar et al. 2017a, when increasing AGN momentum flux).These results imply that BH fueling is hindered by increases in A AGN1 .This form of selfregulation from the BH feedback results in a decrease in energy propagated out to IGM scales.Additionally, a higher value for A AGN1 results in fewer AGN producing feedback in the jet mode (bottom panel).Despite AGN jets reaching maximum velocities starting at higher redshifts and AGN ejecting more mass in feedback events, the fact that less AGN jet feedback occurs with overall less energy available (lower ρBH ) results in no impact on the CDD.
Decreasing the AGN momentum flux increases the number of BHs producing feedback in the jet mode, increases the overall mass of these BHs, and increases the amount of mass these BHs are accreting.However, the ṀBH values of these BHs tend to be larger at higher redshifts which results in a delay of when maximum jet speeds are reached, and the lower value for A AGN1 means less mass is ejected from the accreting SMBH overall.Therefore, less energy overall is being propagated far out to IGM scales.The greater amount of energy available (larger ρBH ) is not able to overcome this fact resulting in a null effect.
The effect of varying the AGN jet speed (A AGN2 ) is more intuitive, as ρBH does not vary significantly.The energy per feedback event does not change for variations in A AGN2 but the fraction of energy going into the outflow velocity does.This allows for jets to travel further before hydrodynamically re-coupling to the gas in the box and depositing kinetic and thermal energy.Increasing jet speed also decreases the number of BHs producing jet feedback at low redshifts.This is likely due to faster jet speeds suppressing stellar growth in halos which prevents SMBHs from being seeded and from growing.However, reducing the number of SMBHs seeded this way does not reduce the impact that the AGN jets have on the Lyα statistics.

SIMBA STELLAR FEEDBACK PARAMETERS
Figure 6 shows that increasing the mass loading, from the fiducial value, increases ρBH of the SMBHs, at z < 1.0, in the jet mode while simultaneously decreasing the number of SMBHs producing feedback in the jet mode, at z ≳ 0.5.The suppression of AGN jet feedback at higher redshifts, due to increases in the SN mass loading factor hinders the ability of the jet feedback to remove HI from the Lyα forest by z = 0.1.This effect is most likely dominated by the slower jet speeds from higher ṀBH rather than the smaller number of SMBHs considering results from other parameter variations.
Mass loading factors smaller than the fiducial value results in lower ρBH and fewer SMBHs producing jet feedback by z = 0.1.Despite AGN jet feedback in the box reaching maximum jet velocities at higher redshifts, it appears the effect of an overall reduced SMBH population producing jets at z < 0.8 dominates.Moreover, lower ρBH means less energy from feedback is in the box.However, the increase in HI absorbers due to larger A SN1 is not as dramatic as the increase due to smaller A SN1 .
More intuitive than the mass loading parameter is the SN wind speed parameter (A SN2 ).Values for ρBH remain quite similar by z = 0 but the number of SMBHs producing jet feedback decreases with increasing SN wind speed and ρBH is dramatically lower at higher z.Due to the large decrease in the number of SMBHs producing jet feedback and the amount of energy from feedback when increasing SN wind speed, the overall HI surviving in the IGM increases and the temperature of that HI tends to be cooler.Strong stellar feedback in the form of faster SN wind speeds efficiently suppresses SMBH growth and feedback in the simulation box.Since the BH seeding mechanism in Simba is tied to the stellar mass of the galaxy, it is clear the lack of SMBHs results from suppressed stellar growth within galaxies due to the fast SN wind speeds.

TNG SMBH Statistics
As previously shown in Tillman et al. (2023), the kinetic mode in the TNG AGN feedback model has no discernible effect on the low-z Lyα forest.From Figure 3 a subtle heating effect can be seen when varying the AGN feedback parameters explored in CAMELS and more prominent heating is seen for variations in stellar feedback.From here, it is clear the feedback parameters affect the temperature in the box to some degree, but on Mpc scales this is small.Despite there being twice   as many SMBHs in TNG as in Simba no effect is seen for the Lyα forest statistics explored herein when varying the CAMELS-TNG AGN feedback parameters.No changes were seen in the SMBH statistics when varying the AGN feedback parameters, either, implying that the AGN feedback parameters in CAMELS-TNG play no role in self-suppression of SMBH growth.This makes sense as both A AGN1 and A AGN2 modify the kinetic feedback mode whereas the thermal feedback mode selfregulates the SMBH growth the most (Weinberger et al. 2017).
Increasing the stellar feedback parameters tends to decrease the value of ρBH and in the case of A SN1 (energy per unit SFR) the number of BHs in the simulation is decreased as well.These results imply that the AGN feedback in TNG affects the forest to some degree since small variations in the Lyα statistics are seen when stronger stellar feedback suppresses BH growth.As discussed previously, the radiative feedback mode is the AGN feedback mode that affects the IGM in TNG.
The radiative mode adds flux directly to the cosmic ionizing background and heats the gas around the host halo.The flux and heat added will scale with the accretion rate of the SMBHs in the box and both of these effects have implications for the gas in the IGM.The increased energy in the ionizing background will be most impactful on the CDD.In the forest, CD is proportional to the inverse of the photoionizing rate but for temperature the relation is weaker, following N HI ∝ T −0.7 .We not only see the temperature scaling affecting the CDD but we also see a slight difference in the b-value distributions since b scales with temperature as in Equation 11.
The AGN thermal feedback mode is more radiatively efficient so changes to ρBH from the thermal mode are expected to have the most impact on Lyα statistics.From Figure 7, increases in both A SN1 and A SN2 relative to the fiducial values show decreases in ρBH (but less so for A SN2 ).This results in less flux added to the UVB and thus more HI absorption in the forest explaining the shift in the TNG CDD seen in Figure 4.There is also a reduction in the number of SMBHs in total for increased A SN1 and a reduction in thermal mode SMBHs for increased A SN2 .These changes to the number of SMBHs seem to have minimal impact most likely due to the change being relatively small when compared to the total number in the box (note the y-axis in the bottom panels only extends from 1000 to 1600).
On the other hand, decreases in A SN1 , relative to the fiducial value, show an increase in the total number of SMBHs and in ρBH .This should result in more flux being added to the UVB and thus reduce the number of absorbers in the forest, which is exactly what is seen from the CDD.Decreases in A SN2 result in a small increase in ρBH and a small decrease in the number of SMBHs in the box.Larger ρBH should result in a stronger UVB and thus less HI absorption but since these changes in ρBH for A SN2 are less dramatic than for A SN1 the effect will be less discernible.The changes to the CDD due to smaller A SN2 are much smaller than changes due to smaller A SN1 which is exactly as expected.
Since the AGN feedback parameters explored thus far have shown no impact on the Lyα forest statistics we also analyze results from additional parameters explored by CAMELS.The CAMELS IllustrisTNG extended 1P set is similar to the original 1P set but explores 22 additional parameters.Variations of all these parameters simultaneously composes the SB28 set (the original 6 parameters plus 22 more results in 28 total parameters Ni et al. 2023).The additional SMBH sub-grid model parameters explored in this set are: The Bondi rate multiplier, the high-accretion mode feedback efficiency, the Eddington rate multiplier for the BH accretion rate limiter, the BH seed mass, the BH radiative efficiency (fraction of rest mass energy released in feedback), the Eddington ratio for the transition between the AGN feedback modes, and the steepness of the mass transition between the AGN feedback modes.As with the original AGN feedback parameters explored herein, we found no observable difference, relative to the observational error bars, in the Lyα forest CDD for independent variations of any of the parameters listed above.At least, no difference was found that showed more dramatic results than the stellar feedback parameters for TNG.It appears that apart from the radiative mode in TNG, the TNG AGN feedback model has essentially no effect on the low-z Lyα forest.

DISCUSSION
In the last section we found that the galactic feedback in Simba, both from AGN and SNe, can have a dramatic impact on the predicted Lyα forest.AGN momentum flux, SN mass loading, and SN wind speed all have the ability to regulate the growth of SMBHs and variations in these parameters can manifest in different predictions for the Lyα forest CDD and b-value distribution -specifically, feedback that lowers the number of SMBHs or decreases the AGN feedback flux will result in more absorption.The AGN jet speed in Simba has the largest direct effect on the predicted Lyα forest statistics explored herein, with faster jets heating a larger swath of the IGM and decreasing the number of absorbers.Indeed, at late times, in fiducial Simba, the AGN heating from jets becomes almost volume filling and the low CD absorbers are particularly vulnerable while the higher CD absorbers are less susceptible to strong shocks on re-coupling (Christiansen et al. 2020;Tillman et al. 2023).This results in an impact on the CDD preferentially at the low CD end, flattening the distributions, as seen in Figure 4.
While the kinetic and thermal AGN feedback in TNG appears not to have an impact on the predicted Lyα forest, it is clear that the AGN radiative feedback mode has at least a small effect.Variations in the stellar feedback parameters in TNG, that can regulate the growth of SMBHs in the box, help illustrate the consequences that AGN have on the predicted Lyα forest in TNG.
Multiple previous studies have posed the use of the low-z Lyα forest as a tool for constraining galactic feedback models and the results herein further motivate that idea.How to move forward with this idea relies on a careful analysis of other processes that have an effect on the predicted Lyα forest not limited to but including the assumed UVB.We also must acknowledge the limitations of analyzing large-scale parameter spaces such as the CAMELS simulations.We discuss these ideas and more in the following sub-sections.

Degeneracy in Lyα Forest Statistics
In Tillman et al. (2023), the degeneracy between the UVB and AGN feedback models was analyzed.That study found that the impact of AGN jet feedback on the CDD resulted in both a simple translation of the distribution function that was degenerate with changes in the UVB, as well as a change in the slope of the CDD that could not be replicated with a change in the UVB.With regards to scaling the UVB, assuming that the Lyα forest is in ionizing equilibrium and optically thin, column density scales with the photoionization rate of hydrogen as N HI ∝ Γ −1 HI .This results in a normalization shift as changing the value of Γ HI for the UVB applies to all the gas in the simulation box, and the effects are largely equivalent for gas in CDs associated with the forest.Temperature changes also affect the CDD due to the aforementioned assumptions leading to the relationship N HI ∝ T −0.7 .Since the heating due to AGN jets is able to reach distances on the order of Mpcs into the diffuse IGM, low density gas is efficiently heated.This means certain CDs of the forest are affected more than others causing a slope change in the CDD.This is the result seen in Tillman et al. (2023) where AGN jet feedback in Simba flattens the CDD.
Understanding other impacts on the forest has important implications for how one might constrain AGN feedback in the future.For example, a better fit to observational data might be achieved with a weaker AGN feedback model accompanied by a slightly stronger UVB.This scenario could be likely in Simba since the fiducial AGN jet mode aggressively blows out gas which can negatively affect group statistics (Robson & Davé 2020, 2021;Oppenheimer et al. 2021;Lovisari et al. 2021;Yang et al. 2022).But as we saw when varying the jet speed, the amount of gas ejected is less important than the ability for that ejected gas to transport heat out to the diffuse IGM.
We also recognize an additional important factor with a degenerate effect on Lyα forest statistics in Simba.The stellar feedback parameters, especially the SN wind speed, can be efficient at suppressing supermassive black hole growth and thus affect the forest.For Simba, increases in the SN wind speed amount to an effect similar to that of decreasing the AGN jet speed.In fact, smaller changes in SN wind speed are required for the same change in the CDD or b-value distribution when varying AGN jet speed.Instead of decreasing the AGN jet speed, one can slightly increase SN wind speed and achieve the same result.

Redshift effects
Figure 9 shows the CDD redshift evolution of TNG vs. Simba.In the plots, the TNG results have been corrected to use the Haardt & Madau (2012) UVB (the same UVB as in Simba) to make the results more directly comparable.We check if the difference between TNG and Simba at z = 2.0 could stem from a different mean flux since the simulations not only have different UVBs but TNG also includes a radiative AGN feedback mode which adds flux to the assumed UVB.Rescaling to the observed effective optical depth, as seen in Kim et al. (2007), results in CDDs at z = 2.0 that are converged dN/dz Simba N HI < 10 14 cm 2 IllustrisTNG N HI < 10 14 cm 2 Simba N HI > 10 14 cm 2 IllustrisTNG N HI > 10 14 cm 2 Figure 10.The redshift evolution of the total number of absorbers per redshift distance for CD bins of 10 12 < NHI < 10 14 cm −2 and 10 14 < NHI < 10 17 cm −2 .Results from CAMELS-Simba and CAMELS-TNG are shown.The dotted grey line marks z = 0.7, the threshold at which the shape of the CDDs diverge significantly between Simba and TNG.Like in Figure 9, the CAMELS-TNG data is corrected in post-processing to utilize the Haardt & Madau (2012) UVB for comparison to the CAMELS-Simba data.
The CDDs of Simba vs TNG at z = 2.0 are almost identical in that they produce similar shapes and abundances.As the universe evolves to lower redshift the differences between the simulations reveal themselves.The slope of the Simba CDD remains largely unchanged while a steepening in the TNG CDD appears above N HI = 10 14 cm −2 becoming more dramatic at lower z.The differences seen between z = 2 and z = 0 could be partially attributed to the density of gas the forest probes at these redshifts and the extent to which feedback might affect those densities.An absorber of CD N HI = 10 14 cm −2 at z = 2 is likely to be probing higher density gas than the same CD absorber at z = 0 (Davé et al. 1999).If we instead looked at similar density gas in the forest at z = 2 we may see a larger difference than we do in Figure 9.
To further explore the difference in the shapes of the evolving CDDs, we analyze the total number of absorbers per redshift distance dN/dz in two different CD bins in Figure 10.The 10 12 < N HI < 10 14 cm −2 CD range probes weak absorbers while the 10 12 < N HI < 10 14 cm −2 range probes strong absorbers.From z = 2 to z ∼ 0.7 the abundance of weak absorbers is about the same between Simba and TNG, but, as seen in their CDDs, the two simulations diverge at lower redshift and that divergence is mainly driven by the loss of absorbers in Simba.AGN heating in Simba likely causes the drop in absorber abundance seen while in TNG the number of absorbers in each bin appears to plateau after z ∼ 0.7.
The difference between the CDDs predicted by Simba and TNG as the universe evolves highlights additional degenerate effects not explored herein.The differences are not due to cosmic variance or cosmology as the CAMELS 1P set simulations use identical initial conditions and cosmological parameters.In Tillman et al. (2023) removing the AGN jet feedback in Simba did not produce a slope change dramatic enough to match the TNG results.However, the simulations explored in Tillman et al. ( 2023) did not have the same box size, initial conditions, and had slightly different cosmological values.Regardless, those results imply that part of the difference seen between the two simulations (starting at z ≲ 0.7) might arise from differences in the temperature and density distribution of the IGM in Simba and TNG that are caused by some factor in addition to AGN jet feedback.A closer look at individual absorbers in the simulations via a cross-correlation method is likely necessary to fully explain these differences.

The effect of cosmic variance.
Figure 11 shows the CDD for Simba and TNG when allowing the initial random seed of the simulation to vary.The maximum variation of the CDD in the CV set due to initial conditions appears within ∼2-2.5 sigma (relative to the smallest observational error bars) and the interquartile range of variation is well within 1 sigma.For TNG, the median of the CV set is well converged to the original TNG300-1 simulation (which has comparable initial mass resolution as CAMELS).For Simba, the median of the CV set is not converged to the original Simba simulations for N HI < 10 14.5 cm −2 and diverges more for lower N HI .However, the divergence between CAMELS and original Simba is on par with the maximum variation from the CV set.This variation comes from the fact that the lowest Lyα forest CD statistics in Simba are extremely sensitive to the AGN jet feedback and the CAMELS-Simba simulations are not volume-converged due to the long-range impact of jets (seen in Borrow et al. 2020, and Gebhardt et al. submitted).
As shown in Figure 4, the lowest column densities are particularly susceptible to variations in the AGN jet feedback mode.Since the AGN jets can reach the outskirts of galaxies before recoupling to the heating and cooling of the simulation (distances on the order of ∼ 10 kpc) they are able to carry most of their energy beyond galaxy scales.This energy, in many cases, is then enough to reach beyond the CGM and halo scales to reach the IGM.That heating will be more effective in lower density gas which could cause the divergence we see approaching lower N HI .In fact, the fiducial CAMELS Simba simulation has an overall hotter IGM than that of original Simba.For the lowest values of n HI , corresponding to more diffuse absorbers, this temperature difference is hotter by a factor of ∼ 1.7.Since in the diffuse IGM N HI ∝ T −0.7 this temperature difference leads to N HI values in original Simba that are ∼ 1.5 times larger than that in CAMELS Simba.Thus the temperature difference is enough to explain the offset between CAMELS and original Simba at the lowest column densities.
The origin of the temperature difference is almost certainly from the number of BHs producing jet feedback in the simulation box.This number will be affected by the random seed due to different halo mass functions being produced in each CAMELS box.A few more massive halos in the box due to a different random seed will result in an increased number of SMBHs large enough to produce impactful long-range jet feedback.A greater number of massive halos with strong AGN jets can be compounded by the small CAMELS box.The periodic boundary conditions in combination with extremely far impacting jets could also result in the energy of these jets being trapped within the box from the point of view of the host halo rather than being dispersed at some distance as it would in reality.
This idea is further motivated by the variance of the TNG CDD due to initial conditions.The maximum variance of the CDD due to different random seeds seen in the CV set is significantly lower than the variance seen in Simba.We saw repeatedly that modifying the AGN feedback parameters explored in the CAMELS TNG suite does not affect the Lyα forest statistics implying that the impact of TNG's AGN feedback might have minimal impact on the forest (apart from the radiative mode which mimics UVB effects).Different AGN feedback parameters may affect the forest or the same AGN parameters explored herein may affect the forest in a larger simulation box.
The minimal impact seen from the AGN feedback implies that the variance of the TNG CDD due to the random seed is likely to be a better indicator of the effect of the random seed on just the population of Lyα absorbers (i.e.whether or not the box size used herein is conducive to a representative population of Lyα absorbers).Since the variance of TNG at low N HI is so much lower than that of Simba, it is likely that another factor (e.g.AGN feedback) is affecting the Simba variance.The large variance seen at higher N HI is similar in both simulations and is due to the rarity of those absorbers.Different random seeds will produce more high CD absorbers than others and this will be especially variable in the small box sizes of CAMELS.See Appendix A for additional discussion on cosmic variance and convergence.

Broader Picture
In Burkhart et al. (2022) the differences between the Illustris and TNG Lyα forest statistics due to different AGN feedback models was found to be difficult to disentangle from effects of the UVB.The CDD in Illustris was able to match that of TNG by a re-scaling of the assumed UVB and the b-value distribution showed no difference between the two simulations.Those results implied that current observations of the CDD and b-values are unable to constrain the different AGN feedback models in Illustris and TNG.However, Burkhart et al. (2022) found that the Lyα flux statistics such as the flux PDF and the 1D flux power spectrum might show observable differences between the two simulations.
Another recent study, Khaire et al. (2023), comparing Illustris and TNG found similar results with the CDD and b-value distribution showing no observable differences but the Lyα flux power spectrum potentially showing observable differences at the highest k values (on small physical scales).Further exploration of the flux power spectrum (particularly around z = 2 to 3) could reveal a dependence on AGN feedback that may become important for certain studies.For example, in the PRIYA simulation suite the effects of AGN feedback were very small (< 0.1%) in the 120 Mpc box and slightly larger in 25 Mpc box but only due to cosmic variance (Bird et al. 2023).However, the AGN feedback model in those simulations resembles closely that of TNG and it is unclear if the Simba model would have more of an effect.
Despite the Illustris and TNG simulation results implying the observed Lyα forest CDD cannot be used to constrain AGN feedback, further work done in Tillman et al. (2023) found that AGN feedback in Simba might be observable.The changes in the Simba simulation CDD when removing AGN jet feedback showed potentially observable differences in the intermediate CD range (N HI ∼ 10 13 cm −2 ) even when allowing the UVB strength to vary to obtain the best fit.Conducting a similar analysis to the one seen in Tillman et al. (2023) on the simulations in this work, we analyze what parameter variations lead to observable changes that can be disentangled from degenerate effects of the UVB.We conduct fits in the CD range of 10 13 < N HI < 10 14.5 cm −2 as this is the range that observation and simulation data is most robust.Reasonable variations of the chosen fitting range does not change the main results of this work.We utilize the observational error bars as weights in our fits and we assume Gaussian distributed random variables.
Figure 8 is the same as Figure 4 but now allowing the strength of the UVB to vary in order to find the best fit to the observed data.Each fit has two free parameters, the value of the feedback parameter being varied (A SN or A AGN ) and the factor the assumed UVB model is multiplied by (Haardt & Madau (2012) for Simba and Faucher-Giguère et al. ( 2009) for TNG).The factor by which the UVB is multiplied and the reduced χ 2 value for a given parameter variation can be found in Table 2 in Appendix B. The most noticeable differences (relative to the observational error bars) in the CDD, when accounting for UVB degeneracy, can be seen for variations in the Simba AGN jet speed A AGN2 and the Simba SN wind speed A SN2 .A noticeable flattening of the CDD can be seen for increases in A AGN2 and decreases in A SN2 .
In particular, the most dramatic differences can be seen when increasing A SN2 from the fiducial value.We have already confirmed that the effects on the forest due to A SN2 are due to the link between stellar feedback and the suppression of SMBH seeding and growth.The dramatic reduction in the number of AGN producing jet feedback results in a much steeper CDD that resembles that of TNG.This effect resembles that of turning off the AGN jets in Simba seen in Tillman et al. (2023), but manifests more prominently due to the fact that the AGN feedback has a larger impact on the forest in the CAMELS-Simba box than in the original Simba run (an idea we discussed to explain the CV set results).The results from Figure 8 imply that the effects from the AGN feedback in Simba could be unique enough to disentangle from the assumed UVB model and that observations of the Lyα CDD could be used to constrain said feedback.In future work, we plan to explore the Lyα forest flux statistics in the context of Simba's AGN feedback to determine if, similar to TNG vs. Illustris, the flux power spectrum might be used to constrain AGN jets.
Finally, we will briefly discuss the discrepancy seen between the observed b-value distribution and the distribution predicted from simulations.Simulations have consistently under-predicted the number of high b absorbers (at low-z) with an observed mean lying around 30 km/s and the simulations predicting a mean around 20 km/s.While we cannot directly compare the CAMELS predicted b-value distribution to observations due to numerical resolution limitations, we can acknowledge the role that faster AGN jets and slower SN wind speeds play roles in broadening the b-value distribution.The heating causing the changes in Figure 5 originates from AGN jet shocks and dispersion of the heat the jets carried into the IGM when said jets re-couple hydrodynamically to the gas in the box.While this heating alone cannot explain the discrepancy between the observed and simulated b-value distribution, it could act as a partial solution.Then the amount of missing turbulence required to fully resolve the statistic would be slightly less than what has been predicted by previous studies (e.g. by Viel et al. 2017;Gaikwad et al. 2017b;Bolton et al. 2022).

CONCLUSION
In this study we use the CAMELS 1P set of simulations to explore effects on Lyα forest statistics due to parameter variation in the Simba and IllustrisTNG feedback models.We find that, in Simba, all four feedback parameters that were varied -the AGN momentum flux, AGN jet speed, SN mass loading factor, and SN wind speed -have clear effects on the Lyα forest CDD and bvalue distribution when ignoring degeneracy due to the assumed UVB.When accounting for the plausible effects of the UVB on the CDD, we found that the AGN jet speed and the SN wind speed showed noticeable differences in the CDD when varied.
For TNG, none of the AGN feedback parameter variations explored in the CAMELS simulation suite affected the Lyα forest statistics in a discernible way.The only effect AGN feedback in TNG appeared to have on the Lyα forest in this study was due to the radiative mode which strengthens the UVB.Varying stellar feedback parameters in TNG (SN energy per unit SFR and SN wind speed) showed minimal effects on both the CDD and the b-value distribution.
In both Simba and TNG we found that adjusting stellar feedback parameters had an indirect effect on the Lyα forest.Stellar feedback in these simulations has the power to limit BH seeding, growth, and feedback to the point of having an impact on the forest statistics.The Simba simulation stellar feedback parameters had a significantly more dramatic effect than that seen in TNG, which follows from the fact that AGN feedback in Simba has a dramatic affect on the IGM whereas TNG AGN does not.Additionally, the Simba SMBH seeding prescription depends on stellar mass rather than halo mass which results in the stellar feedback in Simba having a large effect on the number of SMBHs in the box.We explored the extent to which stellar and AGN feedback affected the BHs in the simulation by looking at Eddington ratios and SMBH accretion rates for relevant AGN feedback modes, the number of SMBHs in each mode, and the number of SMBHs in the box.The main conclusions of our work are as follows: • As found in previous work, Simba's AGN jet feedback plays a dominant role in the low-redshift Lyα forest as the jets heat gas well outside of halos and into the IGM.The distance the heat is transported is more important than the amount of gas ejected.The AGN jet feedback can change the shape of the CDD by flattening it.
• The b-value distribution in Simba is broadened by heating from stronger AGN feedback (faster jets and to a lesser degree higher momentum flux).
Heating from AGN jet feedback (clearly seen in Figure 1) may be a partial solution in resolving the discrepancy between the observed and simulated b-value distribution.
• In agreement with previous work, we find the TNG AGN feedback model has minimal effect on the IGM and thus the Lyα forest.The AGN radiative mode affects the IGM which results in small changes when varying stellar feedback, but the effects of the radiative mode are largely degenerate with that of the UVB.
• Stellar feedback plays a role in SMBH growth suppression and thus should be considered as a degenerate parameter along with the strength of the UVB and AGN feedback when analyzing the low-z Lyα forest statistics.
The CAMELS-Simba simulations for parameter variations closest to the fiducial values produced the best fits to the observed CDD overall when allowing the strength of the UVB to vary.The CAMELS-TNG simulations could not produce a CDD that matches observations for any parameter variation.Both CAMELS-TNG and CAMELS-Simba produced b-value distributions with mean values around 20 km/s which is too low relative to observations.Stronger AGN jet feedback in CAMELS-Simba (via faster jet speeds or lower SN wind speeds) broadened the b-value distribution implying heating from AGN jets could be a partial solution to resolving the predicted b-values from simulations to the observed b-value distribution.However, nonthermal broadening is likely necessary to fully resolve the b-values distribution via un-resolved or un-modeled forms of turbulence in simulations or via other instability mechanisms such as pressure from cosmic rays.
Understanding the interplay between degenerate factors that affect the neutral hydrogen in the IGM is the first step in constraining these mechanisms with observational data.The next step to determine the extent of the degeneracy between the UVB, stellar feedback, and AGN feedback in these sub-grid models is to employ machine learning techniques to the CAMELS project simulations to determine best fits.Many other studies have noted degeneracies between UVB and AGN feedback effects on the forest, and it is clear that constraining feedback models to additional observables will be necessary to unravel the relationship (Burkhart et al. 2022;Tillman et al. 2023;Mallik et al. 2023;Khaire et al. 2023).This will not only help simulations reporoduce a wider range of observables, but will aid in the construction of new more physical feedback models.
Finally it is important to acknowledge that additional factors not discussed in this study likely affect the largescale environment of the IGM.For example, the role of cosmic rays in determining the amount and distribution of HI in the IGM is relatively unconstrained (Lacki 2015;Leite et al. 2017;Butsky et al. 2023).Exploring other cosmological and astrophysical mechanisms that affect the IGM in addition to the UVB and galactic feedback will be essential in fully resolving the discrepancies between the observed and simulated Lyα forest.In future work we plan to do a more thorough analysis of individual absorbers between simulations to determine which factors apart from AGN feedback are defining the forest statistics.We also plan to analyze the Lyα flux power spectrum in the Simba simulation for variations in the AGN feedback models to determine the extent to which those effects might be observable or affect constraints are dark matter properties.Figure 11 shows the range of variation in the CDD from the CV set.Displayed is the median CDD from all the CV simulations for each suite.The shaded region represents the range of variation in the CDD for the CV set.Also displayed are the original Simba and IllustrisTNG CDDs.We compare to the TNG300-1 simulation as it has a comparable mass resolution to that of CAMELS.The CAMELS TNG results are converged to the original TNG results however the original Simba results diverge from the CAMELS Simba results at lower CDs.The Simba divergence is slightly larger than the allowed CV range.Previous work had already found that the TNG CDD is converged for the resolution explored herein (Burkhart et al. 2022).However, for Simba it is unclear if even the original simulation produces a converged CDD due to the lack of higher resolution simulations for comparison.Since the AGN feedback model in the Simba simulations has such a dramatic effect on the predicted CDD, it is possible that a higher resolution simulation may be necessary to properly constrain said feedback model.
Figure 12 shows the median b-value PDF from the CV set as well as the range of variation.The variation in the PDF is higher at lower b-values.This is likely a result of the fitting algorithm as low b-value absorbers are harder to fit via an automated method and are more likely to be noise from a previous fit rather than a unique absorber.The predicted b-value distributions from CAMELS are not converged for any CDs range, but diverge more particularly for lower CDs (N HI < 10 14 cm −2 ).This was determined prior to this work through a resolution study with the IllustrisTNG simulations (Burkhart et al. 2022).As mentioned throughout, we do not directly compare the CAMELS b-value distributions to observations due to the lack of convergence.We expect the CAMELS-TNG results to be close to convergence, but both CAMELS suites lack the mass resolution necessary for confident comparisons.However, we include the Danforth et al. (2016) observational data in Figure 12 for reference.The range of column densities shown in Figure 12 corresponds to that of the Danforth et al. (2016) observational data (10 13 < N HI < 10 14 cm −2 ) and the CAMELS b-value distribution in this range is closer to convergence than if we were to include the N HI < 10 13 cm −2 data.Furthermore, if we only look at absorbers with N HI > 10 14 cm −2 then the CAMELS predicted b-value distributions do converge to the original Simba and TNG simulations exemplifying the importance of mass resolution for low CD absorber statistics.Regardless of convergence, none of the simulations produce a distribution that looks like the observed data, in fact it appears the lower resolution (not converged) distributions actually produce a closer match than the higher resolution results.This is consistent with previous studies' findings that simulations consistently under predict b-values when compared to observations indicating a lack of heating or turbulence in the simulated low-z Lyα forest (Viel et al. 2017;Gaikwad et al. 2017b;Bolton et al. 2022;Burkhart et al. 2022).

B. UVB CORRECTIONS AND BEST FITS
We conduct a least squares fit to find the UVB correction factor required for the best fit of the CDDs predicted by the CAMELS 1P simulations explored herein to the D16 observational data.The fitting procedure is heavily based on that used in Tillman et al. (2023).For our fitting procedure, we assume Gaussian distributed random variables.We conduct this fit within a CD range of N HI = 10 13 cm −2 to 10 14.5 cm −2 as this is where both the simulation data and observations are most robust.Reasonable variations of this fitting range does not change the main results of this work.Excluding lower CD values from the fit due to observational scarcity will not affect the main results of these fits due to the high observational error bars below N HI ≈ 10 13 cm −2 .
We closely follow a UVB correction method as outlined in Kollmeier et al. (2014) which uses the approximation that N HI ∝ 1/Γ HI where Γ HI is the hydrogen photoionization rate.The relation works since the low redshift Lyα forest can be well approximated as an optically thin region in photoionization equilibrium.This method breaks down when absorbers are no longer optically thin but is well converged for CDs explored herein and can be applied in post-processing.
For the reduced χ 2 (χ 2 R ), the number of degrees of freedom is the number of observational points being fit with two variable parameters: the UVB correction factor and the feedback parameter varied in CAMELS.For several parameter variations of CAMELS-Simba, the value for χ 2 R is below 1 implying the CDD is over fitted (this could be partially fixed by removing the large observational error bar points at low CDs).However, many variations including the fiducial results produce χ 2 R values close to 1 exemplifying the remarkable fit to the observed data predicted by CAMELS-Simba.The CAMELS-Simba fits are better than the original Simba fits found in Tillman et al. (2023) due to a further flattening of the CDD in CAMELS-Simba that manifests due to box-size.This effect was discussed in Section 4.
Recent studies have found hydrogen photoionizing values at z = 0.1 that are ∼ 1. 77, 1.78, 2.56, and 1.74 (for Gaikwad et al. 2017a;Khaire & Srianand 2019;Puchwein et al. 2019;Faucher-Giguère 2020, respectively) times stronger than the Haardt & Madau (2012) values.These values can go as high as ∼ 5 times stronger when allowing the escape fraction of HI ionizing photons from galaxies to vary (Khaire & Srianand 2015b).While UVB correction factors lower than Haardt & Madau (2012) found for many of the CAMELS-Simba best fits are disfavored by these more recent UVB model studies, the UVB at low-z is still not well constrained which is why it is often corrected out of Lyα statistics in studies like this one.

Figure 1 .
Figure 1.Simba temperature and column density projections of a single absorber slice as defined in this study (525 kpc/h, see text in 2.3 and 3.1 for justification of this value).The top plots display the projected mean temperature weighted by mass.The bottom plots display the integrated HI column density of the slice.We show projections for varying AGN momentum flux (AAGN1) and AGN jet speed (AAGN2).The left plots are the lowest value for the parameter, the right plots are for the highest value, and the middle plots are the fiducial results.The projections help visualize the effect of varying the AGN feedback parameters.

Figure 2 .
Figure 2. Same as Figure 1 except for variations of the SN mass loading factor (ASN1), and SN wind speed (ASN2).

Figure 3 .
Figure 3.The same temperature projections as seen in Figures 1 and 2 but for the CAMELS IllustrisTNG 1P simulation set feedback parameters.AGN feedback variations are on the left and stellar feedback variations are on the right.The parameters varied are AGN energy per SMBH accretion rate (AAGN1), AGN Burstiness (AAGN2), SN energy per SFR (ASN1), and SN wind speed (ASN2).Column density projections show minimal to no changes and are thus not included.Changes in the temperature projection due to feedback remain close to or confined to the host halos.

Figure 4 .
Figure 4.The Simba and TNG CDDs at z = 0.1 for the CAMELS 1P set varying the AGN feedback parameters (left two plots) and stellar feedback parameters (right four plots).The blue scatter points are observational data from D16.The blue lines represent a decrease in the parameter value (labeled in corresponding color-bars) while the red lines represent an increase.The dashed black lines are the fiducial results.The panel below each plot displays the ratio, ∆, of the parameter variation result to the observed D16 data.The gray shaded region corresponds to the observational error bars.The TNG results when varying the AGN feedback parameters are not shown as no changes in the CDD are seen.

Figure 5 .
Figure 5. Same as Figure 4 but for the b-value distribution and ∆ is the ratio of the parameter variation results to the fiducial results.The distribution is normalized to integrate to unity.No change is seen in the b-value distribution for the TNG simulations when varying AGN feedback so those results are excluded.

Figure 6 .
Figure 6.BH accretion rate density in the jet mode and overall number of SMBHs in the jet mode when varying different feedback parameters in the CAMELS-Simba simulations.All plots: The black lines correspond to fiducial results, the blue lines correspond to the smallest value for the parameter varied, and the red lines correspond to the largest value for the parameter varied.Top row: The SMBH accretion rate density ρBH in the AGN jet feedback mode vs. redshift for different values of the parameters.Bottom row: The number of SMBHs producing jet mode feedback for parameter variations.

Figure 7 .
Figure 7. Similar to Figure 6 but instead for the CAMELS-TNG simulations.The top row display the SMBH accretion rate density for the thermal feedback mode.The bottom row displays the total number of SMBHs (solid lines) along with the number of SMBHs in the thermal mode (dashed lines).The results for variations in the AGN feedback parameters are not shown since those statistics are largely unaffected.

Figure 8 .
Figure 8. Same as Figure 4 but fitting the results of each parameter variation to the observational data by allowing the strength of the UVB to vary.A table including the UVB correction value and reduced χ 2 values for the fit can be found inTable 2 in Appendix B.

Figure 9 .
Figure 9.The redshift evolution of the SIMBA (black solid line) and TNG (red dashed lines) simulation column density distributions.The TNG CDD has been adjusted so as to use the same Haardt & Madau (2012) UVB model as used in Simba.This has been done to more accurately determine what differences are caused by either feedback or large-scale environment as opposed to the ionizing background.
VARIANCE AND CONVERGENCEWe examine the effect that different initial random seeds have on the CDD and b-value distribution.By understanding to what extent the results may change based on initial conditions we can better examine what is an effect of feedback vs. a sampling effect.

Figure 11 .
Figure11.The CDD for all the CAMELS CV set simulations for both the SIMBA and IllustrisTNG suites.The solid black lines are the median of the CAMELS CV CDDs, the shaded region is the 90th percentile range of the CV CDDs, and the dashed red lines are the original Simba and TNG300-1 simulations.Variations in the CDDF due to varying initial conditions results in a maximum shift in the overall normalization by 0.25 dex.Both the original Simba and TNG300-1 simulations are largely contained within the allowed variation of the CDD due to cosmic variance.

Figure 12 .
Figure 12.Same as Figure 11 but instead for the b-value distribution.The blue points are observational data from Danforth et al. (2016).The b-value distributions are calculated only for absorbers in the column density range corresponding to the observational data as written in the figure panels.The x-axis range ends at 100 km/s to avoid blends of multiple lines and continuum fitting errors.

a
Factor by which the Haardt & Madau (2012) UVB is multiplied by at z = 0.1.b Factor by which the Faucher-Giguère et al. (2009) UVB is multiplied by at z = 0.1.

Table 1 .
Summary of CAMELS Simulation Terminology.
a 1P Varies one parameter at a time.Variable parameters are Ωm, σ8, ASN1, ASN2, AAGN1, AAGN2.Each suite has 61 1P simulations.b CV uses fiducial values for all parameters but varies the initial random seed.Each suite has 27 CV simulations.
Pillepich et al. (2018) the local dark matter velocity dispersion and H is the Hubble constant.Additional parameters (ē w , f w,Z , Z w,ref , γ w,Z , N SNII , E SNII,51 , κ w , and v w,min ) are constants with values and descriptions in Table1ofPillepich et al. (2018).The wind mass loading factor All authors acknowledge support for this work by NASA ATP 80NSSC22K0823.MTT thanks the Simons Foundation and the Center for Computation Astrophysics at the Flatiron Institute for providing the computational resources used for this analysis.MTTthanks Mahdi Qezlou for helpful conversations regarding proper usage and trouble shooting of the fake-spectra package.MTT thanks Doug Rennehan for insightful conversations about the difference between the Simba and CAMELS CDD due to cosmic variance.BB is grateful for additional support by NSF grant 2009679, the Packard Fellowship, and the Sloan Fellowship.Part of this work was completed during the KITP workshop Galevo 23.DAA acknowledges support by NSF grants AST-2009687 and AST-2108944, CXO grant TM2-23006X, Simons Foundation Award CCA-1018464, and Cottrell Scholar Award CS-CSA-2023-028 by the Research Corporation for Science Advancement.GLB acknowledges support from the NSF (AST-2108470, XSEDE grant MCA06N030), NASA TCAN award 80NSSC21K1053, and the Simons Foundation (grant 822237) and the Simons Collaboration on Learning the Universe.SH acknowledges support for Program number HST-HF2-51507 provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, incorporated, under NASA contract NAS5-26555.SB acknowledges additional support from NSF AST-2107821.