The Impact of the WHIM on the IGM Thermal State Determined from the Low-z Lyman α Forest

MNRAS 527, 11338–11359 (2024) https://doi.org/10.1093/mnras/stad3846 Advance Access publication 2023 December 14 The impact of the WHIM on the IGM thermal state determined from the lowz Lyman αforest Teng Hu , 1 ‹Vikram Khaire , 1 , 2 Joseph F. Hennawi , 1 , 3 Jose O ˜ norbe , 4 Michael Walther , 5 , 6 Zarija Lukic 7 and Frederick Davies 8 1 Physics Department, Broida Hall, University of California, Santa Barbara, CA 93106-9530, USA 2 Indian Institute of Space Science & Technology, Thiruvananthapuram, Kerala 695547, India 3 Leiden Observatory, Leiden University, PO Box 9513, NL-2300 RA Leiden, the Netherlands 4 Facultad de F ´ ısica, Universidad de Sevilla, Avda. Reina Mercedes s/n, Campus de Reina Mercedes, E-41012 Sevilla, Spain 5 Faculty of Physics, Ludwig-Maximilians-Universit ¨ at M ¨ unchen, University Observatory, Scheinerstr. 1, D-81679 Munich, Germany 6 Excellence Cluster ORIGINS, Boltzmannstr. 2, D-85748 Garching, Germany 7 Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA 8 Max-Planck-Institut f ¨ ur Astronomie, K ¨ onigstuhl 17, D-69117 Heidelberg, Germany Accepted 2023 December 8. Received 2023 December 7; in original form 2023 September 18 A B S T R A C T At z ࣠ 1, shock heating caused by large-scale velocity flows and possibly violent feedback from galaxy formation, converts a significant fraction of the cool gas ( T ∼10 4 K) in the intergalactic medium (IGM) into warm–hot phase (WHIM) with T > 10 5 K, resulting in a significant deviation from the previously tight power-law IGM temperature–density relationship, T = T 0 ( ρ/ ¯ρ) γ−1 . This study explores the impact of the WHIM on measurements of the lowz IGM thermal state, [ T 0 , γ], based on the b –N H I distribution of the Ly αforest. Exploiting a machine learning-enabled simulation-based inference method trained on Nyx hydrodynamical simulations, we demonstrate that [ T 0 , γ] can still be reliably measured from the b–N H I distribution at z = 0.1, notwithstanding the substantial WHIM in the IGM. To investigate the effects of different feedback, we apply this inference methodology to mock spectra derived from the IllustrisTNG and Illustris simulations at z = 0.1. The results suggest that the underlying [ T 0 , γ] of both simulations can be reco v ered with biases as low as |  log ( T 0 /K) | ࣠ 0.05 dex, | γ | ࣠ 0.1, smaller than the precision of a typical measurement. Given the large differences in the volume-weighted WHIM fractions between the three simulations (Illustris 38 per cent, IllustrisTNG 10 per cent, and Nyx 4 per cent), we conclude that the b –N H I distribution is not sensitive to the WHIM under realistic conditions. Finally, we investigate the physical properties of the detectable Ly α absorbers, and disco v er that although their T and  distributions remain mostly unaffected by feedback, they are correlated with the photoionization rate used in the simulation. Key words: galaxy: formation –intergalactic medium –quasars: absorption lines. 1 INTRODUCTION Being the largest reservoir of baryons in the Universe, the intergalactic medium (IGM) plays a crucial role in the evolution of the Universe and the formation of structures. Based on the canonical cosmological model constrained by many observational studies, the thermal evolution of the IGM is dominated by two major phase transition events of the Universe. The first phase transition is hydrogen reionization caused by the first generation of galaxies at redshift 6 < z < 20 (Becker et al. 2001 ; Fan et al. 2006 ; McGreer, Mesinger & D’Odorico 2015 ; Robertson et al. 2015 ). The second one is the double reionization of helium (He II → He III ) driven by quasistellar objects (QSOs; see e.g. Madau & Meiksin 1994 ; McQuinn et al. 2009 ; Khaire 2017 ), which is believed to occur at z ∼3 (Worseck et al. 2011 , 2018 ; Syphers & Shull 2014 ), where the quasar luminosity function reaches its peak (see e.g. Hopkins, Richards &  E-mail: [email protected] Hernquist 2007 ; Khaire & Srianand 2015 ; Kulkarni, Worseck & Hennawi 2019 ). These two-phase transition events heat up the IGM dramatically to a maximum of 15 000 K while ionizing the IGM. After the completion of hydrogen reionization ( z ∼1–2), the IGM thermal state is shaped by the quasi-equilibrium balance between photoionization heating from the extragalactic UV background (Haardt & Madau 2012 ; Khaire & Srianand 2019 ; Faucher-Gigu ` ere 2020 ) and various cooling processes including recombinations, excitation, cooling due to Hubble expansion, and inverse Compton scattering of electrons off of the cosmic microwave background (CMB; see e.g. McQuinn 2016 ). All these processes together drive the IGM to follow a power-law temperature–density ( T – ) relation: T (  ) = T 0  γ−1 , (1) where  = ρ/ ¯ρis the o v erdensity, T 0 is the temperature at mean density, and γis the adiabatic index (Hui & Gnedin 1997 ; Mc- Quinn & Upton Sanderbeck 2016 ). These two parameters [ T 0 , γ] thus characterize the thermal state of the IGM, making it feasible to constrain the IGM thermal history (Lidz et al. 2010 ; Becker ©The Author(s) 2023. Published by Oxford University Press on behalf of Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License ( https:// creativecommons.org/ licenses/ by/ 4.0/ ), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. Downloaded from https://academic.oup.com/mnras/article/527/4/11338/7473705 by Biblioteca Universidad de Sevilla user on 15 May 2024 Impact of the WHIM on IGM thermal state 11339 MNRAS 527, 11338–11359 (2024) Figure 1. Examples of simulation skewers for IllustrisTNG (left) and Illustris (right) simulations, probing the structure generated by the same initial condition, while the two simulations are post-processed to share the same UV backgrounds photoionization rate,  H I . The flux is plotted in black on the top panel, while the temperature T , o v erdensity  , and velocity along LOS v los are shown in black in the second, third, and bottom panels consecutively. et al. 2011 ; Rorai et al. 2017 ; Hiss et al. 2018 ; Walther et al. 2019 ; Gaikwad et al. 2021 ) by measuring [ T 0 , γ] at different epochs. These measurements impro v e our knowledge of the IGM thermal evolution and shed light on the underlying heating and cooling processes involved. Ne vertheless, the aforementioned po wer-law T – relationship for the IGM could potentially break down at z ࣠ 1, where shock heating caused by large-scale velocity flows (Nath & Silk 2001 ; Cen & Ostriker 2006 ) and various feedback mechanisms become more common (Scannapieco, Silk & Bouwens 2005 ; Khaire et al. 2024 ). Specifically, shock heating at lowz converts a notable fraction of the cool IGM into warm hot intergalactic medium (WHIM) with T > 10 5 K (Shull et al. 2012 ), causing a substantial dispersion in the IGM T – distrib ution (Da v ´ e et al. 2001 ; Cen & Ostriker 2006 ). As a result of such dispersion, the IGM T – distribution can no longer be fully described by the typical power-law relationship (see Fig. 1 ), which introduces additional complexities in the measurement of the IGM thermal state (Hu et al. 2022 , hereafter Hu22 ). The imperative question is whether the significant shock heating at lowz influences the observable, i.e. the Ly αforest, which serves as the primary probe of the IGM, and if it does, how might such impacts affect measurements of the IGM thermal state [ T 0 , γ]? In practice, the IGM thermal state can be measured through various statistical properties of the Ly αforest. Particularly, at z ࣠ 3, the Ly αforest is amenable to Voigt profile decomposition (see Hiss et al. 2018 ), where each line can be fit by three parameters: redshift z abs , Doppler broadening b , and neutral hydrogen column density N H I . The IGM thermal state at these redshifts can thus be measured using the 2D joint b –N H I distribution (Schaye et al. 1999 , 2000 ; Bolton et al. 2014 ; Hiss et al. 2018 ; Rorai et al. 2018 ). Hu22 introduced a new inference method to measure the thermal state [ T 0 , γ] and the photoionization rate  H I of the IGM based on the b –N H I distribution and Ly αline density, d N /d z, of the Ly αforest. Such a method performs Bayesian inference with the help of neural networks and Gaussian emulators, trained on a suite of Nyx simulations (Almgren et al. 2013 ; Luki ´ c et al. 2015 ), making it possible to measure the thermal state of the IGM to high precision for realistic mock data sets. Moreo v er, the thermal state of the IGM at z < 1.7 remains poorly constrained, since the Ly αtransition below such redshift lies below the atmospheric cutoff ( λ∼3300 Å), requiring UV observations from space with Hubble Space Telescope (HST). After He II reionization ( z < 3), the thermal state of the IGM is considered to be dominated by adiabatic cooling from Hubble expansion, which leads to an IGM thermal state with T 0 ∼5000 K and γ∼1.6 at the current epoch z = 0 (McQuinn & Upton Sanderbeck 2016 ). Ho we ver, such a prediction of low temperatures has not yet been confirmed observationally. Meanwhile, recent studies have suggested that the Ly αlines appear broader than predicted by numerical simulations at z < 0.5 (Gaikwad et al. 2017 ; Nasir et al. 2017 ; Viel et al. 2017 ; Burkhart et al. 2022 ; Bolton et al. 2022a ). This observation is based on the b parameters acquired from the HST Cosmic Origins Spectrograph (COS) spectra (Danforth et al. 2016 , referred hereafter as D16 ) data set. While it has been argued that such a mismatch might be resolved by additional sources of turbulence, an alternativ e e xplanation would be that the IGM is actually hotter than previously presumed, with T 0 conceivably approaching 10 000 K, implying the existence of unexpected sources of heating (Bolton et al. 2022a , b ), which, if true, would change our understanding of the IGM physics thoroughly. In this paper, we adopt the Hu22 inference method to investigate the impact of the WHIM on measurements of the IGM thermal state, [ T 0 , γ], based on the b –N H I distribution of the Ly αforest. First, we assess the ef fecti veness of [ T 0 , γ] as IGM parameters at lowz by comparing its performance as neural network training labels against the photoheating labels [ A , B ] (see Section 2.3 ). These latter labels are photoheating rate rescaling factors used to generate the Nyx simulation suite with various thermal histories (see e.g Becker et al. 2011 ). Since our emulators are trained on these Nyx simulations generated by varying [ A , B ], the inference method is naturally inclined to retrieve these photoheating labels. On the Downloaded from https://academic.oup.com/mnras/article/527/4/11338/7473705 by Biblioteca Universidad de Sevilla user on 15 May 2024 11340 T. Hu et al. MNRAS 527, 11338–11359 (2024) Figure 2. Volume weighted T – distribution for all three simulations at z = 0.1. The log T for each bin are plotted as solid dots with 1–σT error bars. The best-fitting power-law relationship is shown as dashed lines. The Nyx (left) model is the default model which has log ( T 0 /K) = 3.612, γ= 1.588; and IllustrisTNG (middle) yields log ( T 0 /K) = 3.627, γ= 1.593; whereas Illustris (right) has log ( T 0 /K) = 3.633, γ= 1.577. The gas phase fractions are shown in the annotation. other hand, if shock heating at lowz causes the T – distribution of the Ly αabsorbers to deviate from the power-law relationship, the ef fecti veness of [ T 0 , γ] as labels could be compromised. Thus, our comparison between these two sets of labels provides insight into the robustness of [ T 0 , γ] as IGM parameters at lowz, in the presence of substantial shock heating. Afterwards, we explore the potential effects of different feedback mechanisms, which are associated with WHIM, on measurements of the IGM thermal state, [ T 0 , γ]. In terms of our inference methodology, the question becomes: what would happen if we used a simulation grid without feedback to interpret a Universe that includes feedback? Would this lead to unbiased [ T 0 , γ]? To answer these questions, we apply the Hu22 inference methodology to mock data drawn from the Illustris (Genel et al. 2014 ) and IllustrisTNG (Weinberger et al. 2017 ) simulations at z = 0.1. These two simulations incorporate galaxy formation models and feedback mechanisms that are not included in the Nyx simulation, which heat up the IGM substantially at lowz, and transform the cool diffuse Ly αgas into WHIM more ef fecti vely compared with Nyx simulation (see Fig. 2 ). We examine the inference results based on these two simulations and explore whether feedback biases the measurement of the thermal state [ T 0 , γ]. To further investigate this problem, we explore the specific impacts of shock heating and other astrophysical processes, such as active galactic nucleus (AGN) feedback and UV background photoionization, on the physical properties of the Ly αforest at z = 0.1. Within the three aforementioned simulations, we identify simulated Ly αabsorbers in the simulations and establish a direct correlation between the physical properties of these absorbers (including temperature T , o v erdensity  , and line-of-sight velocity v los ) and the observed Ly αline parameters ( b , N H I ) derived from the absorption lines detected in corresponding mock spectra. We then examine the distributions of  and T of these simulated Ly αabsorbers across the three aforementioned simulations to study the detailed effects of the feedback and UV background photoionization rate,  H I , on the Ly α forest. This paper is organized as follows: In Section 2 , we outline the simulations and associated processes applied to generate synthetic Ly αforest. It includes post-processing, forward-modelling, and Voigt profile fitting. The inference framework and results for all three simulations are then presented in Section 3 . Section 4 is dedicated to the investigation of the physical characteristics of low-redshift Ly αforest absorbers in all three simulations. Finally, in Section 5 , we present a summary and discussion of our findings. For the sake of brevity, we use log as a shorthand to denote log 10 throughout the paper. 2 SIMULATIONS In this paper, we utilize the inference framework described in Hu22 , which employs the b –N H I distribution emulator built on neural networks trained on a set of Nyx simulations. We also use galaxy formation simulations IllustrisTNG and Illustris to investigate the lowz Ly αforest under different feedback mechanisms. Since this work focuses on the low redshift Ly αforest, we use z = 0.1 simulation snapshots for all three simulations. In this section, we first provide a description of the simulations and the implemented physical models, followed by the (mock) data processing procedures employed in our study. This includes the generation of simulated line-of-sight (LOS) of Ly αforest (hereafter referred to as skewers for simplicity), forward modelling, and the Voigt profile fitting of Ly αlines. The cosmological parameters and thermal states of the three simulations are summarized in Table 1 . 2.1 Nyx Nyx is an adaptive mesh, massively parallel, cosmological simulation code primarily developed to simulate the IGM (Almgren et al. 2013 ; Luki ´ c et al. 2015 ). Nyx simulates the dark matter evolution by treating the dark matter as self-gravitating Lagrangian particles, while it models baryons as an ideal gas on a uniform Cartesian grid following an Eulerian approach. The Eulerian hydrodynamics equations are Downloaded from https://academic.oup.com/mnras/article/527/4/11338/7473705 by Biblioteca Universidad de Sevilla user on 15 May 2024 Impact of the WHIM on IGM thermal state 11341 MNRAS 527, 11338–11359 (2024) Table 1. Parameters of cosmology and T – relation (at z = 0.1). Parameters Nyx IllustrisTNG Illustris m 0.3192 0.3089 0.2726  0.6808 0.6911 0.7274 b 0.0496 0.0486 0.0456 h 0.670 0.677 0.704 σ8 0.8288 0.8159 0.809 n s 0.96 0.97 0.963 T 0 4093 K 4241 K 4292 K γ1.588 1.593 1.577 solved using a second-order piece-wise parabolic method, which is capable of accurately capturing shocks. Nyx includes the major astrophysical processes rele v ant to the evolution of the the Ly αforest. First, gas in the Nyx simulation is treated as having a primordial composition with a hydrogen mass fraction of 0.76 and helium mass fraction of 0.24, and zero metallicity. Nyx takes into account the process of inverse Compton cooling off the microwave background and tracks the total thermal energy loss due to atomic collisional processes. Nyx also implements recombination, collisional ionization, dielectric recombination, and cooling following the prescription given in Luki ´ c et al. ( 2015 ). Ionizing radiation in Nyx is modelled by a spatially uniform but timevarying ultraviolet background radiation field of Haardt & Madau ( 2012 ), while assuming all cells in the simulation are optically thin. Furthermore, following standard practice, we allow the UV background photoionization rate,  H I , to be a free parameter in postprocessing while generating mock Ly αskewers. Lastly, Nyx does not implement any galaxy formation or feedback, which simplification reduces the required computational resources significantly, allowing us to run a large ensemble of simulations with different thermal histories (see Section 2.3 ), which is required for accurate statistical inference. Each Nyx simulation model used in this study was initialized with the same initial condition at z = 159 and evolved down to z = 0.03 in a L box = 20 cMpc h −1 simulation box with N cell = 1024 3 Eulerian cells and 1024 3 dark matter particles. The box size is a compromise between computational cost and the need for convergence at least to < 10 per cent on small scales (large k ). In short, such choices of box size and resolution should not affect the line parameters of the Ly α forest significantly. More discussion of the resolution, box size, and convergence issues can be found in Luki ´ c et al. ( 2015 ) and Hu et al. ( 2022 ). 2.2 IllustrisTNG and Illustris To e v aluate the ef fecti veness of the IGM thermal state [ T 0 , γ] as the IGM parameters and test the efficacy of our inference framework on the realistic IGM, which can be affected by astrophysical processes that are not included in Nyx simulation such as galaxy formation and AGN feedback, we employ Illustris (Genel et al. 2014 ) and IllustrisTNG (Weinberger et al. 2017 ; Nelson et al. 2019 ) simulations, and use them as mock observational data in our inference method. The IllustrisTNG and Illustris are cosmological hydrodynamic simulations powered by the AREPO code (Springel 2010 ). This code employs a moving mesh approach to solve hydrodynamics through the Euler equations, and it computes gravitational forces on a quasi- Lagrangian moving Voronoi mesh via the tree-PM method. Both simulations incorporate a wide range of astrophysical processes for galaxy formation, such as star formation, stellar and AGN feedback, galactic winds, and chemical enrichment (Marinacci et al. 2018 ; Naiman et al. 2018 ; Springel et al. 2018 ). They utilize the UV background detailed in Faucher-Gigu ` ere et al. ( 2009 ) for photoionization heating and cooling. Other processes for modelling the Ly αforest, like collisional ionization and inverse Compton cooling from the cosmic microwave background, are also taken into account. The primary distinction between IllustrisTNG and Illustris lies in their AGN feedback mechanisms, especially regarding AGN feedback. Both simulations implement AGN feedback in two modes based on the gas accretion rate on to the central supermassive black hole: the ‘quasar-mode’ at high accretion rates (Springel 2005 ; Hopkins et al. 2008 ; Debuhr, Quataert & Ma 2011 ) and the ‘radiomode’ at low rates (Bower et al. 2006 ; Croton et al. 2006 ; Sijacki et al. 2007 ). While both use continuous thermal feedback in ‘quasarmode’, their ‘radio-mode’ implementations differ. Illustris employs a bubble model for radio-mode feedback, accumulating substantial feedback energy for e xplosiv e release, often ejecting e xcessiv e hot gas (Genel et al. 2014 ). Conversely, IllustrisTNG models this feedback as a kinetic wind, injecting momentum into neighbouring regions from the central black hole. This approach better replicates astrophysical properties like star formation rates and galaxy colour distributions (Nelson et al. 2018 ; Pillepich et al. 2018a , b ). Both the IllustrisTNG and Illustris simulations we used in this study have box sizes of 75 cMpc h −1 and 1820 3 baryon and dark matter particles. Since AREPO is a moving mesh code, we convert the Voronoi mesh outputs to 1820 3 cartesian grids by dumping the smoothed quantities such as temperature. 1 , density, and velocities on grids to generate Ly αforest skewers. A Gaussian kernel with a size equal to 2.5 times the radius of each Voronoi cell is applied for the smoothing, assuming each Voronoi cell is spherical. We then generate skewers for IllustrisTNG and Illustris simulations following the approach discussed in Section 2.4 . In Fig. 1 , we plot two simulation skewers for IllustrisTNG and Illustris respectively, while the two simulations are post-processed to share the same UV backgrounds photoionization rate,  H I (see Section 2.5 for more discussion). The flux ( e −τ) is plotted in the top panel, and the temperature, o v erdensity, and LOS velocity profiles are shown in the second, third, and bottom panels consecutively. It is worth mentioning that the two skewers probe the structure generated by the same initial condition, suggesting that the differences in T ,  , and v los are primarily caused by different feedback strengths. Specifically, the Illustris exhibits higher temperatures due to its stronger feedback, which results in weaker absorption features given the same UV backgrounds. More discussion on the differences between Ly αforest in IllustrisTNG and Illustris simulations can be found in Khaire et al. ( 2024 ) and Khaire et al. ( 2023 ). 2.3 IGM thermal state and parameter grid Following Hu22 , we make use of the Thermal History and Evolution in Reionization Models of Absorption Lines (THERMAL 2 ) suite of Nyx simulations (see Hiss et al. 2018 ; Walther et al. 2019 ) to model the IGM with various thermal histories. The suite consists of 48 models with varying thermal histories, each generated by changing 1 As presented in Appendix A of Martizzi et al. ( 2019 ), there exists a bug that affects the IGM temperature of the IllustrisTNG simulation. Ho we ver, its effect on the Ly αforest is minimal because the bug predominantly impacts the gas with the lowest density . Consequently , we continue to use the uncorrected temperature for the IllustrisTNG simulation. 2 For details of THERMAL suite, see http://thermal.joseonorbe.com Downloaded from https://academic.oup.com/mnras/article/527/4/11338/7473705 by Biblioteca Universidad de Sevilla user on 15 May 2024 11342 T. Hu et al. MNRAS 527, 11338–11359 (2024) Figure 3. Parameters grid (solid circles) from snapshots of Nyx simulations from the THERMAL suite at z = 0.1, parametrized by the thermal state [ T 0 , γ]. The left-hand panel shows the γ–T 0 grid, whose shape is determined by the photoheating labels [ A , B ] (see Fig. A3 ) and the evolution of the thermal state of the IGM. The right-hand panel is γ– H I grid, showing the 13  H I values for each point on the γ–T 0 grid. the photoheating rate of the simulation following the prescription described in Becker et al. ( 2011 ), in which the photoheating rate,  is assumed as a function of o v erdensity, i.e. = HM12 ( z) A B , (2) where HM12 ( z) stands for the time-varying photoheating rate per H II ion tabulated in Haardt & Madau ( 2012 ), and the constants A and B are free photoheating parameters that are varied in the different Nyx runs to achieve different thermal histories, which results in different thermal states at z = 0.1. The distribution of parameters in our thermal grid, i.e. the dif ferent v alues of [ T 0 , γ] are illustrated in Fig. 3 , and the corresponding values of [ A , B ] are presented in Fig. A3 (see Appendix A for more discussion). Conventionally, the thermal parameters [ T 0 , γ] are obtained by fitting a power law to the T – relationships (see equation 1 ). Such a fitting procedure is straightforward at higher redshift ( z  2) where the T – distributions of the IGM are tight. Ho we ver, in lo wz, the distributions of the IGM temperatures are noticeably broader due to the e xtensiv e shock heating, which heats up the IGM, resulting in more WHIM. The T – distributions for all three simulations (Nyx default model with A = 1, B = 0, and IllustrisTNG and Illustris) are shown in Fig. 2 . For each simulation, the gas is divided into four phases depending on the temperature and density, namely the WHIM, Diffuse Ly α, hot halo gas, and condensed, where the cutoffs are set to be T = 10 5 K and  = 120 3 . It can be seen that there exist significant dispersion in the T – distributions of the lowz IGM, i.e. the shock-heated WHIM, for all three simulations, and the fractions of the WHIM are directly proportional to the strength of the feedback. Specifically, f WHIM, Illustris > f WHIM, IllustrisTNG > f WHIM, Nyx , while Illustris implements extreme feedback, IllustrisTNG employs mild feedback and Nyx has no feedback. In order to fit the power-law relationship in the presence of dispersion in the IGM T – distribution, we utilize the fitting procedure presented in Hu22 , which fits the power-law T – relationship by 3 Here, we adopt the cutoff T = 10 5 K, and  = 120 for different gas phases following Dav ´ e et al. ( 2010 ), and more discussion about the different cutoff can be found in Gaikwad et al. ( 2017 ). binning the Diffuse Ly αgas ( T < 10 5 K and  < 120) into 20 bins based on log  , and applying a least-squares linear fit to the mean temperatures of the gas in each bin. Here, we modify the fitting range to −0.5 < log  < 1.5. 4 , which provides a more accurate representation of the  range of the Ly αabsorbers at z ∼0.1, which is the principal subject of this paper. Such a fitting procedure is applied to all simulations used in this study, including all Nyx models and IllustrisTNG and Illustris simulations. The best-fitting power-law relationship based on [ T 0 , γ] and the T – distributions are illustrated in Fig. 2 . The figure shows that although the three simulations yield very different overall T – distributions, their thermal state T 0 and γare ho we ver similar. Furthermore, as described in Hu22 , we vary the UV background photoionization rate,  H I , of the Nyx simulations in post-processing when the simulation skewers are generated, extending the parameter grid to [log T 0 , γ, log  H I ]. The value of  H I we used in this study spans from log (  H I / s −1 ) = −13.834 to −12.932 in logarithmic steps of 0.075 dex, which gives 13 values in total (see the right-hand panel of Fig. 3 ). The range of  H I used here co v ers more than twice the range obtained by UV background models of Khaire & Srianand ( 2019 ) at z = 0.1, achieved by varying the spectral energy distribution of quasars. Note that the range also co v es more than 2 σuncertainty in the  H I measurements (Gaikwad et al. 2017 ; Khaire et al. 2019 ). In total, the 3D thermal grid consists of 48 ×13 = 624 Nyx models. As mentioned earlier, T 0 and γcharacterize the IGM thermal state at z  2, where the IGM is dominated by the power-law T –  relationship. Ho we ver, their ef ficacy as parameters for the IGM thermal state remains uncertain at z ࣠ 1, where a significant fraction of the gas deviates from the power-law T – relationship due to shock heating and feedback. In this paper, we e v aluate the ef fecti veness of the thermal state [ T 0 , γ] as IGM parametrization at lowz using the inference framework presented in Hu22 , and we make use of the photoheating parameters [ A , B ] as an alternative set of labels as 4 Such a choice of fitting range of the power law T – relationship leads to slightly different thermal states [ T 0 , γ] for the three simulations compared with those presented in previous works ( Hu22 ; Khaire et al. 2024 ), but the difference is minor. Downloaded from https://academic.oup.com/mnras/article/527/4/11338/7473705 by Biblioteca Universidad de Sevilla user on 15 May 2024 Impact of the WHIM on IGM thermal state 11343 MNRAS 527, 11338–11359 (2024) Figure 4. One of the forward-modelled mock spectra. The simulated spectrum is shown in black, where the model spectrum determined from VPFIT is shown in blue, and the noise vector is plotted in red. a comparison. These labels are particularly rele v ant since all Nyx models used in the training procedure of our neural network, which is the major component of our inference method, are generated by varying [ A , B ]. This suggests that our inference framework should be capable of reco v ering the values of [ A , B ] efficiently. Therefore, [ A , B ] are particularly useful in the e v aluation of the [ T 0 , γ]. More information about the photoheating labels [ A , B ] is presented in Appendix. A . 2.4 Mock spectra, forward-modelling, and VPFIT We follow the procedure described in Hu22 and generate mock spectra by calculating the Ly αoptical depth ( τ) array along the mock LOS. For each simulation, including all Nyx models and IllustrisTNG and Illustris, an ensemble of 20 000 skewers is created. In this study, we include observational noise and instrumental effects to conduct our analysis under realistic conditions. We generate mock data sets with properties consistent with the D16 compilation of low redshift Ly αforest spectra, which comprises 82 unique quasar spectra observed with the COS on the HST. Among them, 34 spectra co v er our targeted redshift range of 0.06 < z < 0.16 and have signal-to-noise ratio (SNR)s greater than 5 per pixel. These spectra segments contribute to a total observational pathlength of z ob = 2.136. We follow the forward-modelling procedure described in Hu22 . For each mock spectrum, we select one of the 34 spectra segments randomly, stitching skewers to match the wavelength grid of the selected segment 5 . We then convolve these combined skewers with the COS 130M LSF, and add Gaussian noise based on the noise vectors associated with the chosen spectrum segment. This process 5 The length of the Nyx skewers is 25 Mpc h −1 , corresponding to rough z = 0.01. Roughly 10 Nyx skewers are needed for one forwarded-modelled mock spectrum. Similarly, for Illustris and IllustrisTNG, which have box sizes of 75Mpc h −1 , about 3 skewers are used to generate one mock spectrum. ensures that our mock data sets closely replicate the observational characteristics of the COS spectra, allowing for a more accurate and realistic analysis. For each simulation model, including all Nyx simulation models as well as both IllustrisTNG and Illustris simulations, we generated 1000 forward-modelled mock spectra. The total pathlength for each simulation is approximately z tot ∼60, which ensures that our training set for the b–N H I distribution emulator (see Section 3.1 ) is large enough and the resulting b–N H I distribution is not biased by the choice of (mock) spectra. We then use VPFIT (Carswell & Webb 2014 ) 6 to fit the Ly α lines in our simulated spectra to obtain a set of { b, N H I } pairs for all of the mock data sets, following the prescription given in Hu22 . In this paper, as is the convention in lowz Ly αforest analysis, we apply a filter for both b and N H I , and uses only b–N H I pairs with 12 . 5 ≤log ( N H I / cm −2 ) ≤14 . 5 and 0.5 ≤log ( b /kms −1 ) ≤2.5 in our analysis (Schaye et al. 2000 ; Rudie, Steidel & Pettini 2012 ; Hiss et al. 2018 ). A segment of one of the forward-modelled mock spectra is shown in Fig. 4 . The simulated spectrum is shown in grey, where the model spectrum determined from VPFIT is shown in blue, and the noise vector is plotted in red. The top panels of Fig. 5 display 1D histograms of both b (left) and N H I (right) for all three simulations, and the bottom panels illustrate their relati ve dif ferences when compared to the Nyx simulation (as discussed in Section 2.5 , the three simulations used here are d N /d z matched). The median value for log b and log N H I are indicated by dashed vertical lines for each simulation. Notably, while the median values of both b and N H I are comparable across the three simulations, there are distinct differences in the distributions of both parameters across the three simulations. We also notice that the differences in the b parameters are more significant across the three simulations compared with N H I . 6 VPFIT: http:// www.ast.cam.ac.uk/ ∼rfc/ vpfit.html Downloaded from https://academic.oup.com/mnras/article/527/4/11338/7473705 by Biblioteca Universidad de Sevilla user on 15 May 2024 11344 T. Hu et al. MNRAS 527, 11338–11359 (2024) Figur e 5. Top: Mar ginalized 1D b (left) and N H I (right) distributions for all three simulations. For each simulation, the { b, N H I } data set is obtained by VP-fitting an ensemble of 1000 forward-modelled mock spectra. The median values for log b and log N H I are indicated by dashed vertical lines. bottom: The relative difference compared with Nyx simulation. The three simulations used here are d N /d z matched. 2.5 Photoioniztion rate  H I and d N /d z It is noteworthy that the three simulations used in this study by default have different UV background photoionization rates  H I (for Nyx, here we are referring to the default model with log ( T 0 /K) = 3.612 and γ= 1.588.) This is because the photoionization rate  H I are tuned in post-processing across all three simulations to ensure they exhibit the same absorber density d N /d z as the one we measured from D16 data set at z = 0.1. Specifically, we apply the aforementioned VP-fitting procedure to D16 spectra (segments) with 0.06 < z < 0.16, and obtain d N /d z = 167.3 for absorbers within the limits 12 . 5 ≤log ( N H I / cm −2 ) ≤14 . 5 and 0.5 ≤log ( b /kms −1 ) ≤2.5. Such matching of d N /d z is analogous to the matching of the mean flux of simulations at highz for optically thin absorbers (Luki ´ c et al. 2015 ). To match this d N /d z, we tune the photoionization rate, following the prescription described in Section 2.3 , and set log (  H I , / s −1 ) = −13.093, −13.021, −13.414 for Nyx, IllustrisTNG, and Illustris, respectively (see Fig. 6 ). Such mismatch in  H I is caused by the de generac y between the photoionization rate and different feedback recipes used in the simulations. Since both the UV background and feedback suppress the formation of Ly αabsorbers (Khaire et al. 2024 ). More specifically, the feedback heats up the IGM, converting a significant amount of the diffuse Ly αgas into WHIM, which reduces the Ly αtransmission caused by the neutral hydrogen H I in the cool diffuse Ly αgas. To this end, simulations with stronger feedback exhibit lower d N /d z under the same  H I . We measure the d N /d z for the three simulations, including all Nyx simulation models and IllustrisTNG and Illustris, each based on its respective set of 1000 forward-modelled mock spectra. The relationships between UV background photoionization rate and d N /d z for all three models are shown in Fig. 6 , where the d N /d z for Nyx is plotted in blue, IllustrisTNG in green, and Illustris in red, while the d N /d z for the D16 data at z = 0.1 is shown as the horizontal dashed–dotted grey line. Fig. 6 demonstrates that while Illustis has the strongest feedback, which causes more gas to be collisionally Figure 6. d N /d z versus  H I for all three simulations at z = 0.1. Nyx (default model) is shown in blue, IllustrisTNG in green, and Illustris in red, while the observed d N /d z calculated from D16 data set for the corresponding redshift are shown in the horizontal grey dashed–dotted line. The  H I values used for each simulation to match the observed d N /d z are indicated by vertical dashed lines with the corresponding colour. ionized, reducing the Ly αabsorption, it requires the lowest  H I to match the d N /d z to the observed value, and IllustrisTNG, with mild feedback, has higher for the same UV background. Interestingly, whereas Nyx employs no feedback mechanism, it requires slightly lower  H I compared with IllustrisTNG, which implements feedback, to achieve the same d N /d z. Such a trend, which is opposite to the correlation between d N /d z and feedback strength (as seen between Illustris and IllustrisTNG), is caused by the small difference in Downloaded from https://academic.oup.com/mnras/article/527/4/11338/7473705 by Biblioteca Universidad de Sevilla user on 15 May 2024 Impact of the WHIM on IGM thermal state 11345 MNRAS 527, 11338–11359 (2024) the  , T distributions in Nyx and IllustrisTNG. More specifically, IllustrisTNG exhibits a slightly higher mass fraction of the diffuse Ly αgas in the particular  , T range that is probed by the Ly α forest 7 . It is possible that while the mild feedback in the IllustrisTNG simulation results in a slightly higher WHIM fraction compared with Nyx, it also produces more gas with T ∼10 4.5 K and  ∼10, which is the  , T range probed by the Ly αforest at z = 0.1 (see Section 4 ). Ho we ver, the detailed astrophysical mechanism leading to this specific  –T distribution in IllustrisTNG at z = 0.1 is still unclear, and we leave it to our future work. If not otherwise specified, the three simulations used in this study, including Illustris, IllustrisTNG, and Nyx default model, are tuned to have the same Ly αline densities, with d N /d z = 167.3, which is the same value we measured from the D16 data set. 3 INFERENCE METHOD 3.1 Emulating the { b, N H I } distribution In this work, we make use of the inference framework following Hu22 , which measures the thermal state and the photoionization rate  H I of the low redshift IGM using its b –N H I distribution and absorber line density d N /d z. The b –N H I distribution emulator is built on density-estimation likelihood-free inference (DELFI), which turns inference into a density estimation task by learning the distribution of a data set as a function of the labels or parameters (Papamakarios & Murray 2016 ; Alsing, Wandelt & Feeney 2018 ; Lueckmann et al. 2018 ; Papamakarios, Sterratt & Murray 2018 ; Alsing et al. 2019 ). Following Hu22 , we make use of pydelfi , the publicly available python implementation of DELFI, 8 which makes use of neural density estimation (NDE) to learn the sampling conditional probability distribution P ( d | θ) of the data summaries d , as a function of labels/parameters θ, from a training set of simulated data. Here, the data summaries d are [ log N H I , log b ], and our two sets of label parameters θare the thermal state [log T 0 , γ, log  H I ] and photoheating labels [ A , B , log  H I ]. The  H I grids are identical for the two sets of labels. We generate two training data sets by labelling the { b, N H I } pairs obtained from our simulated spectra with the two sets of labels respectively. We train the neural network on the summaryparameter pairs for each training data set separately. Our b–N H I distribution emulator learns the conditional probability distribution P ( b , N H I | T 0 , γ, log  H I ) and P ( b , N H I | A, B, log  H I ) from the corresponding training data set. These conditional b–N H I distributions are then used in our inference algorithm, where we try to find the best-fitting model given the observational/mock data set, which is described in the following section. 3.1.1 Likelihood function In Bayesian inference, a likelihood L = P ( data | model ) is used to describe the probability of observing the data for an y giv en model. We adopt the likelihood formalism introduced in Hu22 , which is 7 This is different from the mass-weighted or volume-weighted diffuse Ly α fraction, f Ly α, which is defined to include all gas with T ≤10 5 K and  ≤ 120 following Dav ´ e et al. ( 2010 ). Instead, the gas probed by the observed HST COS Ly αForest has a narrow range of  and T which also depends on the  H I value used in the simulation (see Section 4 for more details). 8 See https:// github.com/ justinalsing/ pydelfi summarized as follows: ln L = n  i= 1 ln ( μi ) −d N d z model z data , (3) where μi is the Poisson rate of an absorber occupying a cell in the b– N HI plane with area  N H I i ×bi, i.e. μi = d N d z model P ( b i , N H I ,i | θ) N H I b z data . (4) The P ( bi, N H I i | θ) in the equation is the probability distribution function at the point ( bi, N H I , i ) for any given model parameters θe v aluated by the DELFI b –N H I distribution emulator described in Section 3.1 . The z data is the total redshift pathlength co v ered by the quasar spectra from which we obtain our { b, N H I } data set, and ( d N / d z ) model is the absorber density which is e v aluated for any given set of parameters using a Gaussian process emulator (based on George , see Ambikasaran et al. 2016 ), which is also trained on our training data sets obtained from the Nyx simulation suite. To perform our analysis under realistic conditions, all tests performed in this paper are based on mock data sets consisting of 34 forward-modelled spectra, each corresponding to one of the 34 D16 quasar spectra, which gives these data sets the same pathlength as the observation data set with z ob = 2.136. Each of the mock data sets is constructed by randomly selecting 34 spectra from the 1000 forward-modelled spectra, while making sure that each of the 34 D16 quasar spectra is represented exactly once, thereby maintaining the integrity and representativeness of our mock data sets. An example of the MCMC posterior obtained based on the aforementioned likelihood function is given in Fig. 7 . The inference is conducted using the labels [ T 0 , γ, log  H I ]. The posterior appears compact, with the medians of the marginalized posteriors landing close to the true parameters for all three parameters, i.e. within 1 σerrors for marginalized 1D distributions. The b –N H I distribution reco v ered from the mock data set is presented in Fig. 8 , which is emulated by our b –N H I distribution emulator, trained on [ T 0 , γ, log  H I ], based on the inferred parameters, i.e. median values of the marginalized 1D MCMC posterior. The plot exhibits a good match between the mock data set (black dots) and the reco v ered b –N H I distribution (colour map). As a comparison to the IGM parametrization based on the thermal state, [ T 0 , γ, log  H I ], the inference result derived from the same mock data set using the photoheating labels [ A , B , log  H I ] is given in Appendix. A . 3.2 Inference test An inference test is an ef fecti ve method to e v aluate the robustness of a given inference algorithm, which usually consists of approximations and emulation/interpolation procedures that might induce additional uncertainties, altering the error budget. In practice, an inference test can be conducted by performing a set of realizations of the inference method using mock data sets and e v aluating the robustness of the resulting posterior probability distributions, which can be quantified by the co v erage probability P cov (Prangle et al. 2014 ; Ziegel & Gneiting 2014 ; Morrison & Simon 2018 ; Sellentin & Starck 2019 ), the proportion of the time that the true parameters used to generate a mock data set are contained within the posterior contour corresponding to a certain probability level P inf . Such calculations can be performed for many different probability levels, resulting in a series of co v erage probabilities. F or perfect inference, this co v erage probability P cov is al w ays equal to the probability level of the chosen posterior contour P inf (shown as the black dashed line in Fig. 9 ). Downloaded from https://academic.oup.com/mnras/article/527/4/11338/7473705 by Biblioteca Universidad de Sevilla user on 15 May 2024 11346 T. Hu et al. MNRAS 527, 11338–11359 (2024) Figure 7. An example of posterior obtained by our inference method based on inference labels [ T 0 , γ, log  H I ]. Projections of the thermal grid used for generating models are shown as blue dots, while the true model is shown as red dots. The inner (outer) black contour represents the projected 2D 1(2) σinterval. Red lines in the marginal distributions indicate the parameters of true models, while the dashed black lines indicate the 16, 50, and 84 percentile values of the marginalized 1D posterior. The true parameters are: log ( T 0 /K) = 3.612 and γ= 1.588, while log (  H I / s −1 ) = −13.093. In this study, we make use of the inference test described in Wolfson et al. ( 2022 ), which calculates the co v erage probability based on the MCMC posteriors. Compared with the one used in Hu22 , this inference test algorithm is more precise and automatically returns full co v erage probabilities from 0 to 1 rather than co v erage probabilities at only a few specific probability levels. To e v aluate the ef fecti veness of [ T 0 , γ, log  H I ] as IGM parameters for inference at lowz, where the IGM T -  distribution is no longer characterized by the power-law relationship, we perform inference tests based on different sets of labels. We compare the result of the inference test based on labels [ T 0 , γ, log  H I ] with the one based on the photoheating labels [ A , B , log  H I ]. As discussed in Section 2.3 , the comparison between these two sets of labels sheds light on the efficacy of [ T 0 , γ] as IGM parameters at lowz, where the perv asi ve shock heating causes significant dispersion in the IGM T – distribution. For each set of labels, we ran 300 realizations of our inference method, each based on a model randomly chosen from the grid. We exclude models that are close to the boundaries to mitigate boundary effects caused by the hard cutoff of the inference prior, which leads to the truncation of the posteriors. For [ T 0 , γ, log  H I ] gird, we specify 3.3 < log ( T 0 /K) < 3.9, 1.0 < γ< 2.3, −13 . 75 < log (  H I / s −1 ) < −13 . 0. We then calculate the full co v erage probabilities based on all 300 MCMC posteriors. The results of the inference tests are shown in Fig. 9 , where the x -axis stands for inferred probability P inf , and the y -axis shows the co v erage probability, P cov . The shaded re gions indicate the 1- σerror for P cov , which is calculated based on the binomial distribution. The y = x black-dashed line represents a perfect inference test. It can be seen that for Nyx simulations, our inference method is mildly o v erconfident, and the thermal state [ T 0 , γ] (blue) performs slightly better than the photoheating labels [ A , B ] (red), i.e. P cov / P inf is closer to unity . Quantitatively , for inference based on the thermal state [ T 0 , γ], the 68 per cent contour contains the true parameters 61.2 ±2.8% of the time, and the 95 per cent contour contains the true parameters 90.4 ±1.6% of the time. The results show that the [ T 0 , γ] are still robust inference labels for the IGM at lowz, although the shock heating alters the T -  distribution. This further suggests that shock heating alone does not significantly change our understanding in determining the thermal state of the IGM using the Ly αforest. Lastly, while the general efficacy of the inference framework remains robust, we attribute its imperfections to two primary sources: deficiencies within the neural network used in our inference algorithm, and the boundary Downloaded from https://academic.oup.com/mnras/article/527/4/11338/7473705 by Biblioteca Universidad de Sevilla user on 15 May 2024 Impact of the WHIM on IGM thermal state 11353 MNRAS 527, 11338–11359 (2024) Figure 16. Distributions of  and T of the simulated Ly αabsorbers in the Nyx simulation, compared with the full simulation. The ensemble consists of 10 510 absorbers, all obtained from the 1000 spectra discussed in Section 4.1 . The top panels show the 1D distributions of T (right) and  (left) for the whole simulation (green) compared with simulated Ly αabsorbers (blue). The medians of the T and  for the simulated absorbers are indicated by dashed black lines. The bottom panels plot the 2d T – distributions for the Ly αabsorbers(left) and for the whole simulation(right), while the contours for 1 σ, 3 σ(68 and 99.7 per cent) of the T – distribution of the absorbers are shown in both panels. The volume-weighted gas phases for absorbers and the whole simulation are given in the left-hand panel and the right-hand panel, respectively. The best-fitting power-law T – relationships are given in the bottom panels as comparisons. summarize the physical properties of the simulated absorbers. It is worth mentioning that here the IllustrisTNG and Illustris simulations are tuned to have identical d N /d z, which requires different  H I values (see Section 2.5 ). We plot the marginalized 1D distributions of the  and T for both IllustrisTNG (top) and Illustris (bottom) in Fig. 17 . The plots show that the o v erall distributions of T and  for the two simulations are e vidently dif ferent due to their different feedback recipes. F or instance, the e xtreme feedback in Illustris simulation results in much more WHIM compared with IllustrisTNG, causing a secondary peak in its T distribution. Ho we ver, the distributions of T and  for the absorbers in both simulations are comparable, with log ( T med /K) = 4.33, log  med = 0.97 for IllustrisTNG, and log ( T med /K) = 4.16, log  med = 0.72 for Illustris. Moreo v er, we disco v er that the differences in log ( T med /K) and log  med for the three simulations are actually caused by the dif ferent  H I v alues used for the three simulations. The rele v ant discussion is presented in Section 4.3 . The (volume-weighted) 2D T – distributions for simulated Ly α absorbers in both IllustrisTNG (top) and Illustris (bottom) simulations are shown in the left column of Fig. 18 , whereas the (volumeweighted) 2D T – distributions for the whole simulations are given in the right column as comparisons. The volume-weighted gas fractions are given in the annotation for simulated absorbers and full simulation in the left-hand and the right-hand panels, respectively. For the simulated Ly αabsorbers, 12.2 per cent (10.7 per cent) of the Ly αabsorbers arise from the WHIM for IllustrisTNG (Illustris), while the value for Nyx simulation is approximately 7 per cent. The 1 σand 3 σ(68 and 99.7 per cent) contours for the T – distributions for the simulated Ly αabsorbers are also given in the Figure, showing that their T -  distributions are more scattered compared with these in Nyx simulation, especially for the WHIM phase absorbers. These differences are caused by stronger shock heating in IllustrisTNG and Illustris simulations compared with Nyx simulation, caused by their feedback mechanisms. Ho we ver, while the (volume-weighted) WHIM fractions for the two simulations are remarkably different, Downloaded from https://academic.oup.com/mnras/article/527/4/11338/7473705 by Biblioteca Universidad de Sevilla user on 15 May 2024 11354 T. Hu et al. MNRAS 527, 11338–11359 (2024) Figur e 17. Mar ginalized 1D  and T distributions of the simulated Ly αabsorbers in the IllustrisTNG (top) and Illustris (bottom) simulation. The medians of the T and  for the simulated absorbers are indicated by dashed black lines. The o v erall  and T distributions of the full simulations are plotted as comparisons. The two simulations share the same d N /d z, which is the same value observed in the D16 data set. 9.8 per cent for IllustrisTNG and 38.0 per cent for Illustris, the WHIM fractions for the Ly αabsorbers are similar, both around 11 per cent. Furthermore, in Section 4.3 we demonstrate that the small difference in WHIM fractions for simulated absorbers in the two simulations actually arises from the different  H I values used in the two simulations. Such a fact implies that the lowz Ly αforest does not probe the WHIM ef fecti vely under realistic conditions, which is consistent with the conclusion drawn by Khaire et al. ( 2024 ). 4.3 Simulations under the same  H I Considering that the calculation of the Ly αoptical depth τLy α involves  H I , and given that the observed absorption feature (i.e. the Ly αforest) consistently probes regions with τLy α∼1, it follows that the T and  of these regions probed by the Ly α forest, are influenced by the  H I values. Such an argument can be qualitatively demonstrated by the fluctuating Gunn–Peterson approximation (FGPA; see Weinberg et al. 1997 ) τLy α∝ x HI n H ∝ n 2 H T −0 . 7  HI ∝  2 . 7 −γ  HI ∝ T 2 / ( γ−1) −0 . 7  HI , (6) where the τLy αdenotes the Ly αoptical depth and the n H is the hydrogen number density. Since the Ly αforest al w ays probes the region with τLy α∼1, the last two terms in equation ( 6 ) suggest that the  H I is in positive correlation with  and T, respecti vely, gi ven the γ∼1.6 at z = 0.1. In our analysis, the three simulations are tuned to match d N /d z. Ho we ver, due to the degeneracy between  H I and feedback mechanisms, each simulation ends up with a distinct  H I value (refer to Section 2.5 ). As a result, the T and  distributions of the simulated Ly αabsorbers in these simulations are influenced not just by the feedback but also by the varying  H I values. To isolate and examine solely the impact of feedback, we post-process the IllustrisTNG and Illustris simulations to align with the  H I value used in Nyx, set at log (  H I /s −1 ) = −13 . 093. With this consistent  H I across the three simulations, we re-perform the analysis from the prior section and summarize the results below. It is worth mentioning that the o v erall T – distributions of simulations are determined by the cooling and heating processes during their evolution and are not altered by the post-processing of the  H I . We plot the marginalized  and T distributions and their median values for Ly αabsorbers in Nyx, IllustrisTNG, and Illustris simulations with the same d N /d z in Fig. 19 . Interestingly, for simulations Downloaded from https://academic.oup.com/mnras/article/527/4/11338/7473705 by Biblioteca Universidad de Sevilla user on 15 May 2024 Impact of the WHIM on IGM thermal state 11355 MNRAS 527, 11338–11359 (2024) Figure 18. T – distributions of the Ly αabsorbers in the IllustrisTNG (top) and Illustirs (bottom) simulations, compared with the T – distributions of full simulations (right). The contours for 1 σ, 3 σ(68 and 99.7 per cent) of the T – distribution of the absorbers are shown in both panels. The volume-weighted gas phases for absorbers and the whole simulation are given in the left-hand panel and the right-hand panel, respectively. The best-fitting power-law T – relationships are given as comparisons. The two simulations are d N /d z matched. with the same d N /d z, the T and  for absorbers are correlated with its  H I . More specifically, with  H I , Illustris <  H I , Nyx <  H I , IllustrisTNG (see Fig. 6 ), we obtain T med, Illustris < T med, Nyx < T med, IllustrisTNG and  med, Illustris <  med, Nyx <  med, IllustrisTNG . We plot the marginalized  and T distributions and their median values for Ly αabsorbers in the three simulations under the same  H I in Fig. 20 . Under the same  H I , the T and  distributions for simulated absorbers in all three simulations become almost identical, having nearly the same median values for T and  , respectively. Such a result suggests that while feedback e vidently af fects the o v erall T – distributions of the lowz IGM (see Fig. 2 ), their impacts on the physical properties of the lowz Ly αforest (i.e. the T and  distributions) are not distinguishable under realistic conditions. In Fig. 21 , we plot the 2D T – distributions of the simulated Ly αabsorbers in IllustrisTNG (left) and Illustris (right), under the same  H I . While the o v erall IGM T – distributions for the two simulations are evidently different (see Fig. 2 ), the T – distributions of the simulated Ly αabsorbers in these two simulations are similar, and the gas phase fractions for absorbers in both simulations are almost identical, suggesting that the small difference in the WHIM fractions of the simulated absorbers shown in Fig. 18 are caused by different  H I values. Such results indicate that the  H I has a much stronger impact on the Ly αforest compared with the feedback mechanisms implemented in IllustrisTNG and Illustris simulations. 5 SUMMARY AND DISCUSSION In this paper, we explore the effects of the WHIM, which causes significant dispersion in the IGM T – distribution, on the lowz Ly αforest and the IGM thermal state [ T 0 , γ] measured from it. We first e v aluate the ef fecti veness of [ T 0 , γ] as IGM parameters under the inference framework presented in Hu22 , and compare its performance with the photoheating labels [ A , B ]. We disco v er that the thermal state [ T 0 , γ] still parametrizes the IGM ef fecti vely in spite of the dispersion in the IGM T – distribution. We further apply the inference method to IllustrisTNG and Illustris simulations which implement dif ferent v ariants of feedback, potentially making them Downloaded from https://academic.oup.com/mnras/article/527/4/11338/7473705 by Biblioteca Universidad de Sevilla user on 15 May 2024 11356 T. Hu et al. MNRAS 527, 11338–11359 (2024) Figur e 19. Mar ginalized 1d distribution of T (left), and  (right) of the simulated Ly αabsorbers identified in Nyx (blue), IllustrisTNG (green), and Illustris (yellow) simulations. The medians of log T and log  are indicated by vertical dashed lines. The three simulations are turned to have identical d N /d z. Figur e 20. Mar ginalized 1d distribution of T (left), and  (right) of the simulated Ly αabsorbers identified in Nyx (blue), IllustrisTNG (green), and Illustris (yellow) simulations. The medians of log T and log  are indicated by vertical dashed lines. The three simulations used here are post-processed to have the same UV background photoionization rate, with  H I = −13.093. better approximations to the real Universe. We discover that the [ T 0 , γ] of these two simulations can be reco v ered using the inference method within reasonable offsets. Considering the inference results and the huge difference across the three simulations in the IGM WHIM fractions, we conclude that the Ly αforest does not probe the WHIM ef fecti vely under realistic conditions, and the IGM thermal state [ T 0 , γ] is not affected by the shock heating caused by AGN feedback and other astrophysical processes significantly at z = 0.1. To further confirm our conclusion, we identified the Ly α absorbers in all three simulations at z = 0.1, and pair them to the corresponding absorption lines identified in the mock spectra. The physical properties of the simulated Ly αabsorbers support our conclusion that the observable Ly αforest are not affected by the substantial WHIM in the lowz, and the thermal state [ T 0 , γ] measured from the Ly αforest remains solid. In this section, we summarize our paper and present our discussion as follows. (i) We compare the performance of [ T 0 , γ] as neural network training labels against the photoheating labels [ A , B ], i.e. the photoheating rate rescaling factors used to generate the Nyx simulation suite with various thermal histories. Given that the [ A , B ] parameters were actually used to generate the simulation outputs, one might e xpect that the y would serv e as a better set of labels than [ T 0 , γ]. Ho we ver, the inference results show the efficacy of these two sets of labels are comparable, suggesting that the [ T 0 , γ], which parametrize the power law T – relationship, still ef fecti vely characterize the Ly αobserv ables at lo wz, notwithstanding the dispersion in the T – distribution induced by shock heating at lowz. (ii) We explored the degree to which the presence of feedback can influence or bias the inference of the IGM thermal state parameters from the b–N H I distribution. In the context of our inference framework, this question becomes: what would happen if we used a simulation grid without feedback to infer the thermal state of Downloaded from https://academic.oup.com/mnras/article/527/4/11338/7473705 by Biblioteca Universidad de Sevilla user on 15 May 2024 Impact of the WHIM on IGM thermal state 11357 MNRAS 527, 11338–11359 (2024) Figure 21. T – distributions of the simulated Ly αabsorbers in the IllustrisTNG (left) and Illustirs (right) simulations under the same  H I . The contours for 1 σ, 3 σ(68 and 99.7 per cent) of the T – distribution of the absorbers are shown in both panels. The volume-weighted gas phase fractions for absorbers in both simulations are given, respectively. a Universe that has strong feedback? Would the feedback lead to biased inference? To address these questions, we apply our inference procedure trained on Nyx simulations without feedback to mock data sets from the IllustrisTNG and Illustris simulations which include feedback, whereby the latter serve as potential proxies for the real Universe. We find that the [ T 0 , γ] of IllustrisTNG and Illustris can be reco v ered within small offset, where  log ( T 0 /K) = −0.022 dex, γ = 0.064 for IllustrisTNG and  log ( T 0 /K) = 0.047 dex, γ = −0.094 for Illustris. These offsets are smaller than the typical precision afforded by a realistic data set, i.e.  log T 0 ࣠ 0.5 σlog T 0 , and γ ࣠ σγ. (iii) We developed a method to identify regions in the simulation responsible for the Ly αabsorption lines identified via Voigtprofile fitting, allowing us to determine their temperature T and o v erdensity  from the simulation skewers. For the Nyx simulations, the simulated Ly αabsorbers have a median density log  median = 0.92, a median temperature T median = 10 4.27 K, and about 7 per cent of the simulated Ly αabsorbers have T > 10 5 , making them outliers from the power-law T – relationship. This low fraction is consistent with the previous study of Tepper-Garc ´ ıa et al. ( 2012 ) on the lowz BLAs. (iv) As pointed out in previous work (Bolton et al. 2022a ; Khaire et al. 2024 ; Tillman et al. 2023a ), the Ly αforest is affected by the UV background, which impacts the d N /d z. Nevertheless, we observe that the temperature and o v erdensity of the region probed by the lowz Ly αforest are also affected by the UV background photoionization rate  H I used in the simulation. For d N /d z matched simulations, the T and  of the simulated Ly αabsorbers are correlated with its  H I respectively . Specifically , the Ly αforest probes regions with higher T and  given a higher  H I . This is because for Ly αabsorbers with τLy α∼1, the fluctuating Gunn–Peterson approximation implies that  H I ∝  2 . 7 −γ∝ T 2 / ( γ−1) −0 . 7 , where γ∼1.6. (v) We post-processing the three simulations to share the same  H I , allowing us to explore the effects of different mechanisms. Under the same  H I , the T and  distributions of the simulated Ly αabsorbers across all three simulations become almost indistinguishable, converging to nearly identical median values, while the o v erall IGM T – distributions remain different among the simulations, due to their distinct feedback mechanisms. For the WHIM fractions, the volume-weighted WHIM fractions for IllustrisTNG and Illustris stand at 9.8 and 38.0 per cent, respectively, but the WHIM fractions for the simulated Ly αabsorbers in both simulations are nearly identical, averaging around 11.6 per cent. This suggests that while feedback significantly alters the lowz IGM T – distribution, especially the WHIM phase gas, their impacts on the lowz Ly αforest is indistinguishable under realistic conditions. Such a conclusion aligns with the results derived from various statistics of the lowz Ly αforest by Khaire et al. ( 2024 ). We have thus far demonstrated the robustness of the thermal state [ T 0 , γ] as IGM parameters at lowz, in spite of the dispersion in the T – distribution induced by shock heating. We also pro v ed that the Hu22 inference framework can effectively measure the thermal state [ T 0 , γ] notwithstanding the feedback mechanisms implemented in the IllustrisTNG and Illustris simulations. Looking ahead, we plan to apply the Hu22 inference methodology to simulations with more flexible and sophisticated feedback mechanisms, including EAGLE (Schaye et al. 2015 ) and CAMELS suite (Villaescusa-Navarro et al. 2021 ). The outcomes will provide us with a deeper understanding of the impact of various feedback processes on lowz IGM. Moreo v er, by applying our methodology on archi v al HST COS and space telescope imaging spectrograph (STIS) data sets, we expect precise measurements of the lowz IGM thermal state. These results will pinpoint the onset of the discrepancy in the b parameters of the lowz Ly αforest between current simulations and observations, which is essential for unravelling the underlying physics and acquiring a comprehensive picture of the IGM thermal evolution at lowz after the epoch of helium reionization. ACKNOWLEDGEMENTS We thank the members of the ENIGMA 10 , Siang Peng Oh, Timothy Brandt, and K.G. Lee for helpful discussions and suggestions. Calculations presented in this paper used the hydra and draco clusters of the Max Planck Computing and Data Facility (MPCDF, formerly known as RZG). MPCDF is a competence center of the Max Planck Society located in Garching (Germany). This research 10 http:// enigma.physics.ucsb.edu/ Downloaded from https://academic.oup.com/mnras/article/527/4/11338/7473705 by Biblioteca Universidad de Sevilla user on 15 May 2024 11358 T. Hu et al. MNRAS 527, 11338–11359 (2024) also used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility located at La wrence Berkele y National Laboratory, operated under Contract No. DE-AC02-05CH11231. In addition, we acknowledge PRACE for awarding us access to JUWELS hosted by JSC, Germany. JO acknowledges support from grants BEAGAL18/00057 and CNS2022-135878 from the Spanish Ministerio de Ciencia y Tecnologia. DATA AVAILABILITY The simulation data and analysis code underlying this article will be shared on reasonable request to the corresponding author. REFERENCES Almgren A. S. , Bell J. B., Lijewski M. J., Luki ´ c Z., Van Andel E., 2013, ApJ , 765, 39 Alsing J. , Wandelt B., Feeney S., 2018, MNRAS , 477, 2874 Alsing J. , Charnock T., Feeney S., Wandelt B., 2019, MNRAS , 488, 4440 Ambikasaran S. , F oreman-Macke y D., Greengard L., Hogg D. W., O’Neil M., 2016, IEEE Trans. Pattern Anal. Mach. Intell. , 38, 252 Becker R. H. et al., 2001, AJ , 122, 2850 Becker G. D. , Bolton J. S., Haehnelt M. G., Sargent W. L. W., 2011, MNRAS , 410, 1096 Bolton J. S. , Becker G. D., Haehnelt M. G., Viel M., 2014, MNRAS , 438, 2499 Bolton J. S. , Caputo A., Liu H., Viel M., 2022a, Phys. Rev. Lett. , 129, 211102 Bolton J. S. , Gaikwad P., Haehnelt M. G., Kim T.-S., Nasir F., Puchwein E., Viel M., Wakker B. P., 2022b, MNRAS , 513, 864 Bower R. G. , Benson A. J., Malbon R., Helly J. C., Frenk C. S., Baugh C. M., Cole S., Lacey C. G., 2006, MNRAS , 370, 645 Burkhart B. , Tillman M., Gurvich A. B., Bird S., Tonnesen S., Bryan G. L., Hernquist L., Somerville R. S., 2022, ApJ , 933, L46 Carswell R. F . , W ebb J. K., 2014, Astrophysics Source Code Library, record ascl:1408.015 Cen R. , Ostriker J. P., 2006, ApJ , 650, 560 Croton D. J. et al., 2006, MNRAS , 365, 11 Danforth C. W. et al., 2016, ApJ, 817, 111 Dav ´ e R. et al., 2001, ApJ , 552, 473 Dav ´ e R. , Oppenheimer B. D., Katz N., Kollmeier J. A., Weinberg D. H., 2010, MNRAS , 408, 2051 Debuhr J. , Quataert E., Ma C.-P., 2011, MNRAS , 412, 1341 Fan X. et al., 2006, AJ , 132, 117 Faucher-Gigu ` ere C.-A. , 2020, MNRAS , 493, 1614 Faucher-Gigu ` ere C.-A. , Lidz A., Zaldarriaga M., Hernquist L., 2009, ApJ , 703, 1416 Gaikwad P. , Srianand R., Choudhury T. R., Khaire V., 2017, MNRAS , 467, 3172 Gaikwad P. , Srianand R., Haehnelt M. G., Choudhury T. R., 2021, MNRAS , 506, 4389 Garzilli A. , Theuns T., Schaye J., 2015, MNRAS , 450, 1465 Genel S. et al., 2014, MNRAS , 445, 175 Haardt F. , Madau P., 2012, ApJ , 746, 125 Hiss H. , Walther M., Hennawi J. F., O ˜ n orbe J., O’Meara J. M., Rorai A., Luki ´ c Z., 2018, ApJ , 865, 42 Hopkins P. F. , Richards G. T., Hernquist L., 2007, ApJ , 654, 731 Hopkins P. F. , Hernquist L., Cox T. J., Kere ˇ s D., 2008, ApJS , 175, 356 Hu T. et al., 2022, MNRAS , 515, 2188 (Hu22) Hui L. , Gnedin N. Y., 1997, MNRAS , 292, 27 Khaire V. , 2017, MNRAS , 471, 255 Khaire V. , Srianand R., 2015, MNRAS , 451, L30 Khaire V. , Srianand R., 2019, MNRAS , 484, 4174 Khaire V. et al., 2019, MNRAS , 486, 769 Khaire V. , Hu T., Hennawi J. F ., W alther M., Davies F ., 2024, MNRAS , 527, 4545 Khaire V. , Hu T., Hennawi J. F., Burchett J. N., Walther M., Davies F., 2023, preprint ( arXiv:2311.08470 ) Kulkarni G. , Worseck G., Hennawi J. F., 2019, MNRAS , 488, 1035 Lidz A. , Faucher-Gigu ` ere C.-A., Dall’Aglio A., McQuinn M., Fechner C., Zaldarriaga M., Hernquist L., Dutta S., 2010, ApJ , 718, 199 Lueckmann J.-M. , Bassetto G., Karaletsos T., Macke J. H., 2018, preprint ( arXiv:1805.09294 ) Luki ´ c Z. , Stark C. W., Nugent P., White M., Meiksin A. A., Almgren A., 2015, MNRAS , 446, 3697 Madau P. , Meiksin A., 1994, ApJ , 433, L53 Marinacci F. et al., 2018, MNRAS , 480, 5113 Martizzi D. et al., 2019, MNRAS , 486, 3766 McGreer I. D. , Mesinger A., D’Odorico V., 2015, MNRAS , 447, 499 McQuinn M. , 2016, ARA&A , 54, 313 McQuinn M. , Upton Sanderbeck P. R., 2016, MNRAS , 456, 47 McQuinn M. , Lidz A., Zaldarriaga M., Hernquist L., Hopkins P. F., Dutta S., Faucher-Gigu ` ere C.-A., 2009, ApJ , 694, 842 Morrison J. , Simon N., 2018, J. Comput. Graph. Stat. , 27, 648 Naiman J. P. et al., 2018, MNRAS , 477, 1206 Nasir F. , Bolton J. S., Viel M., Kim T.-S., Haehnelt M. G., Puchwein E., Sijacki D., 2017, MNRAS , 471, 1056 Nath B. B. , Silk J., 2001, MNRAS , 327, L5 Nelson D. et al., 2018, MNRAS , 475, 624 Nelson D. et al., 2019, Computa. Astrophys. Cosmol. , 6, 2 Papamakarios G. , Murray I., 2016, Advances in Neural Information Processing Systems, The MIT Press, Cambridge, MA, p. 1028 Papamakarios G. , Sterratt D. C., Murray I., 2018, preprint ( arXiv:1805.07226 ) Pillepich A. et al., 2018a, MNRAS , 473, 4077 Pillepich A. et al., 2018b, MNRAS , 475, 648 Prangle D. , Blum M. G. B., Popovic G., Sisson S. A., 2014, Aust. N. Z. J. Stat. , 56, 309 Robertson B. E. , Ellis R. S., Furlanetto S. R., Dunlop J. S., 2015, ApJL , 802, L19 Rorai A. et al., 2017, Science , 356, 418 Rorai A. , Carswell R. F., Haehnelt M. G., Becker G. D., Bolton J. S., Murphy M. T., 2018, MNRAS , 474, 2871 Rudie G. C. , Steidel C. C., Pettini M., 2012, ApJ , 757, L30 Scannapieco E. , Silk J., Bouwens R., 2005, ApJ , 635, L13 Schaye J. , Theuns T., Leonard A., Efstathiou G., 1999, MNRAS , 310, 57 Schaye J. , Theuns T., Rauch M., Efstathiou G., Sargent W. L. W., 2000, MNRAS , 318, 817 Schaye J. et al., 2015, MNRAS , 446, 521 Sellentin E. , Starck J.-L., 2019, J. Cosmol. Astropart. Phys. , 2019, 021 Shull J. M. , Harness A., Trenti M., Smith B. D., 2012, ApJ , 747, 100 Sijacki D. , Springel V., Di Matteo T., Hernquist L., 2007, MNRAS , 380, 877 Springel V. , 2005, MNRAS , 364, 1105 Springel V. , 2010, MNRAS , 401, 791 Springel V. et al., 2018, MNRAS , 475, 676 Syphers D. , Shull J. M., 2014, ApJ , 784, 42 Tepper-Garc ´ ıa T. , Richter P., Schaye J., Booth C. M., Dalla Vecchia C., Theuns T., 2012, MNRAS , 425, 1640 Tillman M. T. et al., 2023a, AJ ,166,228 Tillman M. T. , Burkhart B., Tonnesen S., Bird S., Bryan G. L., Angl ´ es-Alc ´ azar D., Dav ´ e R., Genel S., 2023b, ApJ , 945, L17 Viel M. , Haehnelt M. G., Bolton J. S., Kim T.-S., Puchwein E., Nasir F., Wakker B. P., 2017, MNRAS , 467, L86 Villaescusa-Navarro F. et al., 2021, ApJ , 915, 71 Walther M. , O ˜ norbe J., Hennawi J. F., Luki ´ c Z., 2019, ApJ , 872, 13 Weinberg D. H. , Hernsquit L., Katz N., Croft R., Miralda-Escud ´ e J., 1997, in Petitjean P., Charlot S.eds, Structure and Evolution of the Intergalactic Medium from QSO Absorption Line System. p. 133 Paris: Editions Frontieres Paris, France Weinberger R. et al., 2017, MNRAS , 465, 3291 Wolfson M. , Hennawi J. F., Davies F. B., O ˜ norbe J., 2022, MNRAS,521,4056 Worseck G. et al., 2011, ApJ , 733, L24 Worseck G. , Davies F. B., Hennawi J. F., Prochaska J. X., 2018, APJ,875,111 Ziegel J. F. , Gneiting T., 2014, Electron. J. Stat. , 8, 2619 Downloaded from https://academic.oup.com/mnras/article/527/4/11338/7473705 by Biblioteca Universidad de Sevilla user on 15 May 2024 Impact of the WHIM on IGM thermal state 11359 MNRAS 527, 11338–11359 (2024) Figure A1. An example of posterior obtained by our inference method based on [ T 0 , γ,  H I () ]. Projections of the thermal grid used for generating models are shown as blue dots, while the true model is shown as red dots. The inner (outer) black contour represents the projected 2D 1(2) σinterval. The parameters of true models are indicated by red lines in the marginal distributions, while the dashed black lines indicate the 16, 50, and 84 percentile values of the posterior. The true parameters are: A = 1.0 and B = 0.0, log (  H I / s −1 ) = −13.093. APPENDIX A: INFERENCE BASED ON THE PHOTOHEATING LABELS [ A , B ] In this section, we present our inference results using the framework where different Nyx models are labelled by the photoheating parameters [ A , B ] instead of the thermal state [ T 0 , γ], and the inference method returns [ A, B, log  H I ]. The inference is conducted following the procedures described in Section 3 , based on the DELFI b–N H I distribution emulator trained on training data set labelled by [ A , B , log  H I ], which returns P ( b , N H I | A, B, log  H I ). The simulation grid, parametrized by the photoheating labels [ A , B ], is given in Fig. A3 . An example of the MCMC posterior obtained based on the aforementioned likelihood function is given in Fig. A1 . The inference method returns A = 1.321 (1.0), B = −0.190 (0.0),  H I = −13 . 160 ( −13.093), whereas the true values are given in the parentheses. The posterior appears compact, with the medians of the marginalized posteriors landing within 1 σerrors for all three parameters. The b–N H I distribution reco v ered from the mock data set is presented in Fig. A2 , which is emulated by our DELFI b–N H I distribution emulator based on the inferred parameters. We perform an inference test following the Section 3.2 , in which we also exclude models that are too close to the parameter boundaries to a v oid the truncation of the resulting posteriors. Specifically, we only use models with 3.3 < log ( T 0 /K) < 3.9, 1.0 < γ< 2.3, −13 . 75 < log (  H I / s −1 ) < −13 . 0. The result of the inference test is shown in Fig. 9 . The performance looks comparable to the one based on the thermal state [ T 0 , γ], suggesting that [ T 0 , γ] are still ef fecti ve IGM parameters at lowz, notwithstanding the substantial dispersion in the IGM T -  distribution induced by perv asi ve shock heating at this redshift. Figure A2. The colour map is the full b–N H I distribution reco v ered from the Nyx mock data set, which is emulated by our DELFI emulator based on the best-fitting parameters (median values of the marginalized MCMC posterior), where A = 3.695 (1.0) and B = 1.507 (0.0) and log (  H I / s −1 ) = −13.237 ( −13.093), the true parameters are given in parentheses. Black dots are the mock data sets we used in the inference. The contours correspond to cumulative probabilities of 68, 95, and 99.7 per cent. For illustration purposes, the values of pdf are multiplied by 100 in the colour bar. Figure A3. Parameters grid (blue circles) from snapshots of hydrodynamic simulations of the THERMAL suite at z = 0.1 parametrized by the thermal state [ A , B ]. This paper has been typeset from a T E X/L A T E X file prepared by the author. ©The Author(s) 2023. Published by Oxford University Press on behalf of Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License ( https://cr eativecommons.or g/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. Downloaded from https://academic.oup.com/mnras/article/527/4/11338/7473705 by Biblioteca Universidad de Sevilla user on 15 May 2024