The properties of the warmhot intergalactic medium WHIM in cosmic filaments are among the least quantified units in modern astrophysics. The Spectrum Roentgen GammaeROSITA All Sky Survey SRGeRASS provides a unique opportunity to study the Xray emission of the WHIM. We applied both imaging and spectroscopic stacking techniques to the data of the first four eRASS scans to inspect the Xray emissions from 7817 cosmic filaments identified from Sloan Digital Sky Survey SDSS optical galaxy samples. We obtained a 9 significant detection of the total Xray signal from filaments in the 0.31.2 keV band. Here we introduce a novel method to estimate the contamination fraction from unmasked Xray halos active galactic nuclei and Xray binaries associated with filament galaxies. We found an approximately 40 contamination fraction for these unmasked sources suggesting that the remaining 60 of the signal could be coming from the WHIM and a 5.4 detection significance of the WHIM. Moreover we modeled the temperature and baryon density contrast of the detected WHIM by fitting the stacked spectrum and surface brightness profile. The bestfit temperature logTK  6.84  0.07 obtained by using a single temperature model is marginally higher than in the simulation results. This could be due to the fitting of a single temperature model on a multitemperature spectrum. Assuming a 0.2 solar abundance the bestfit baryon density contrast logSUBbSUB  1.88  0.18 is in general agreement with the Xray emitting phases in the IllustrisTNG simulation. This result suggests that the broadband Xray emission traces the high end of the temperature and density values that characterize the entire WHIM population.
['intergalactic medium', 'large-scale structure of Universe', 'X-rays: diffuse background', 'Astrophysics - High Energy Astrophysical Phenomena', 'Astrophysics - Cosmology and Nongalactic Astrophysics']
The SRGeROSITA allsky survey Xray emission from the warmhot phase gas in long cosmic filaments
{'X-ray emission from the warm-hot phase gas in long cosmic filaments': 'X. Zhang 1 , E. Bulbul 1 , N. Malavasi 1 , V. Ghirardini 1 , 2 , J. Comparat 1 , M. Kluge 1 , A. Liu 1 , A. Merloni 1 , Y. Zhang 1 , Y. E. Bahar 1 , E. Artis 1 , J. S. Sanders 1 , C. Garrel 1 , F. Balzer 1 , M. Brüggen 3 , M. Freyberg 1 , E. Gatuzz 1 , S. Grandis 5 , S. Krippendorf 4 , 6 , K. Nandra 1 , G. Ponti 7 , 1 , M. Ramos-Ceja 1 , P. Predehl 1 , T. H. Reiprich 8 , A. Veronica 8 , M. C. H. Yeung 1 , and S. Zelmer 1  \n- 1 Max-Planck-Institut fur Extraterrestrische Physik, Giessenbachstrasse, 85748, Garching, Germany e-mail: xzhang@mpe.mpg.de\n- 2 INAF, Osservatorio di Astrofisica e Scienza dello Spazio, via Piero Gobetti 93 / 3, 40129 Bologna, Italy\n- 3 Hamburger Sternwarte, Gojenbergsweg 112, 21029, Hamburg, Germany\n- 4 University Observatory Munich, Faculty of Physics, Ludwig-Maximilians-Universität, Scheinerstr. 1, 81679, Munich, Germany\n- 5 Universität Innsbruck, Institut für Astro- und Teilchenphysik, Technikerstrasse 25 / 8, 6020, Innsbruck, Austria\n- 6 Arnold Sommerfeld Center for Theoretical Physics, LMU Munich, Theresienstr. 37, 80333 Munich, Germany\n- 7 INAF, Osservatorio Astronomico di Brera, Via E. Bianchi 46, 23807, Merate, (LC), Italy\n- 8 Argelander Institute for Astronomy, University of Bonn, Auf dem Hügel 71, 53121, Bonn, Germany', 'ABSTRACT': 'The properties of the warm-hot intergalactic medium (WHIM) in cosmic filaments are among the least quantified units in modern astrophysics. The Spectrum Roentgen Gamma / eROSITA All Sky Survey ((SRG / eRASS) provides a unique opportunity to study the X-ray emission of the WHIM. We applied both imaging and spectroscopic stacking techniques to the data of the first four eRASS scans to inspect the X-ray emissions from 7817 cosmic filaments identified from Sloan Digital Sky Survey (SDSS) optical galaxy samples. We obtained a 9 σ significant detection of the total X-ray signal from filaments in the 0.3-1.2 keV band. Here, we introduce a novel method to estimate the contamination fraction from unmasked X-ray halos, active galactic nuclei, and X-ray binaries associated with filament galaxies. We found an approximately 40% contamination fraction for these unmasked sources, suggesting that the remaining 60% of the signal could be coming from the WHIM and a 5 . 4 σ detection significance of the WHIM. Moreover, we modeled the temperature and baryon density contrast of the detected WHIM by fitting the stacked spectrum and surface brightness profile. The best-fit temperature log( T / K) = 6 . 84 ± 0 . 07, obtained by using a single temperature model, is marginally higher than in the simulation results. This could be due to the fitting of a single temperature model on a multi-temperature spectrum. Assuming a 0.2 solar abundance, the best-fit baryon density contrast log ∆ b = 1 . 88 ± 0 . 18 is in general agreement with the X-ray emitting phases in the IllustrisTNG simulation. This result suggests that the broadband X-ray emission traces the high end of the temperature and density values that characterize the entire WHIM population.  \nKey words. large-scale structure of Universe - intergalactic medium - X-rays: di ff use background', '1. Introduction': "The warm-hot intergalactic medium (WHIM) permeating cosmic filaments, with a temperature of 10 5 -7 K and density of n H < 10 -4 cm -3 , makes up the majority of the 'missing baryon' in the low-redshift Universe predicted by numerical simulations (e.g., Fukugita et al. 1998; Cen & Ostriker 1999; Shull et al. 2012; Tuominen et al. 2021). Recent measurements on the dispersions of localized fast radio bursts (FRBs) that trace the total line-of-sight free electron column density quantify the cosmic baryon density in the low-redshift Universe in a new way (Macquart et al. 2020; Yang et al. 2022; Wang & Wei 2023). These FRB cosmic baryon density measurements are consistent with the results from other probes, such as cosmic microwave background (CMB Planck Collaboration et al. 2016, 2020) and big bang nucleosynthesis (Pitrou et al. 2018), which conclude that baryons are not missing in the low-redshift Universe. However, the FRB dispersion measure technique cannot provide additional detailed information on the baryons, such as density and temperature. Exploring these properties using other observation tech-  \nniques, such as X-ray observations and Sunyaev-Zeldovich (SZ) e ff ects, remains one of the primary investigation goals in modern astrophysics (Driver 2021).  \nCurrently, absorption studies by using OVI lines in the far ultraviolet bands can successfully trace the 10 5 . 5 K lowtemperature phases WHIM (Danforth et al. 2016). For X-ray absorption lines that trace ∼ 10 6 K phases WHIM, while there are many cases of plausible detections (Nicastro et al. 2005; Fang et al. 2010; Nicastro et al. 2018), these constraints suffer from potential contamination from other absorption sources (e.g. Johnson et al. 2019; Dorigo Jones et al. 2022; Gatuzz et al. 2023). This method is also currently restricted to a handful of bright background blazars (see Fang et al. 2022 for a review). Meanwhile, X-ray emission studies have been mostly focused on the brightest part of the WHIM with a relatively high density contrast, which resides in bridges connecting galaxy cluster pairs or along the filaments (Werner et al. 2008; Bulbul et al. 2016; Akamatsu et al. 2017; Reiprich et al. 2021) or circumcluster space of massive galaxy clusters (Eckert et al. 2015; Bulbul  \net al. 2016; Reiprich et al. 2021). Similarly, the number of these systems is limited, and there is less information about the properties of the WHIM in the environments of long cosmic filaments.  \nWide-area optical spectroscopic surveys have mapped the large-scale galaxy distributions up to redshifts of 0.7 (Reid et al. 2016), thus o ff ering opportunities for a complete examination of the WHIM emission in the vast space of the cosmic filaments. Two pilot studies have shown mild positive detections (4 . 2 σ and 3 . 8 σ ) of stacked X-ray surface brightness emission from cosmic filaments using ROSAT All-Sky Survey (RASS) and extended ROentgen Survey with an Imaging Telescope Array (eROSITA) Final Equatorial-Depth Survey (eFEDS) data (Tanimura et al. 2020b, 2022), demonstrating the potential of Xray stacking analysis to trace the emission in cosmic filaments. Similarly, mild positive detections have also been reported by stacking the SZ e ff ect signal that traces the pressure of the intergalactic medium (de Graa ff et al. 2019; Tanimura et al. 2019, 2020a). However, in these previous works, detections in the very soft X-ray band < 0 . 6 keV at the observed frame, needed to constrain the temperature of the redshifted T < 10 7 K gas, are still lacking. Meanwhile, the energy resolution was limited in the observer frame due to the broad redshift distribution of the filament sample, limiting the ability to measure the plasma temperature in the soft band. Moreover, the pilot X-ray image stacking studies did not account for the signal from unmasked X-ray sources, which could contribute a significant fraction of the total emission. In this work, we improve upon the previous stacking analyses by extending the bandpass for spectral analysis, increasing the energy resolution by using a rest-frame spectral stacking scheme, accounting for the emission from unmasked low-mass halos and point sources and utilizing the new eROSITA All-Sky Survey (eRASS) data. This gives us access to an extensive area with an unprecedented depths, with excellent soft band sensitivity that is essential for detecting the faint X-ray emission from WHIM gas.  \nIn this work, we assume a fiducial Lambda cold dark matter cosmology with parameters: H 0 = 70 km s -1 Mpc -1 , Ω m = 0 . 3, ΩΛ = 0 . 7. We adopt the ratio between baryon density parameter to matter density parameter from the Planck 18 cosmology results, namely: Ω b / Ω m = 0 . 158 (Planck Collaboration et al. 2020). At z = 0, the mean baryon density is ρ b (0) = 4 . 51 × 10 -31 g cm -3 . The baryon density contrast is defined as ∆ b ≡ ρ b /ρ b . Throughout this article, log refers to the logarithm to the base 10, and ln refers to the natural logarithm. Unless otherwise indicated, all sizes, scales, and lengths are given in physical (proper) scale. This article is organized as follows. In Sect. 2, we describe the X-ray data products and the optical filament catalog we used. In Sect. 3, we present the surface brightness profile stacking method and results. In Sect. 4, we estimate the fraction of detected signal due to unmasked galactic sources. Section 5 presents the spectral stacking method and temperature measurements. The estimation of the detected X-ray emitting phases and the comparison with numerical simulations are presented in Sect. 6. We conclude our work in Sect. 7.", '2. Observations, optical filament catalog, and data preparation': 'The stacked first four scans of the eROSITA All-Sky Survey (hereafter eRASS:4) data were collected from December 12, 2019 to December 19, 2021. The data were processed with eSASS (Brunner et al. 2022) pipeline version 020, which is the same as the version 010 used for eRASS1 data release with  \nFig. 1. Smoothed 0.3-1.2 keV eROSITA count rate map of the analysis footprint. The missing pixels show source-masked regions.  \n<!-- image -->  \nimprovements on boresight correction, detector noise suppression, and subpixel position computation (Merloni et al. 2024). We only selected events from telescope modules (TMs) 1, 2, 3, 4, and 6 (hereafter, TM8) to avoid the systematic uncertainties caused by the optical light leak in TMs 5 and 7 (Predehl et al. 2021). We used the tools in the software package eSASS version eSASSusers\\_211214\\_0\\_4 to generate eRASS data products.  \nWe adopted an optical filament catalog (Malavasi et al. 2020), compiled using the DisPerSE software (Sousbie 2011; Sousbie et al. 2011) and based on the Sloan Digital Sky Survey (SDSS) LOWZ + CMASS galaxy catalogs (Reid et al. 2016). The filament catalog contains 63391 filaments after a 3 σ persistence cut without any smoothing applied on the galaxy density field. The original filament catalog is stored in the comoving frame and Cartesian coordinates. Therefore, we converted the coordinates back to (right ascension, declination, redshift) using the conversion functions and the cosmology used in Malavasi et al. (2020). The physical length of the filaments was calculated using the assumed cosmological parameters above. The cross-section of the SDSS LOWZ + CMASSfootprint (Reid et al. 2016) and the western Galactic hemisphere of the eROSITA AllSky Survey determines a ∼ 3200 deg 2 common analysis footprint. In addition, the bright foreground regions of the eROSITA Bubble (Predehl et al. 2020), the Virgo Cluster (Böhringer et al. 1994; McCall et al. 2024), and the Galactic plane | b | < 20 · were excluded from the analysis for simplicity. The total analysis footprint before source masking has a solid angle is shown in Fig. 1, amounting to 2275 deg 2 . In this footprint, we select 7817 filaments in a redshift range 0 . 2 < z < 0 . 6 and a physical length range 20 Mpc < l fil < 100 Mpc for the stacking analysis (see Fig. 2).', '2.1. Source masking': "We masked out all eRASS:4 X-ray sources in the 0.2-2.3 keV eSASS detection catalog 1 . The detection pipeline used in this work is the same as the version used for the first release of the All-Sky Survey data (Merloni et al. 2024). The catalog contains point-like and extended sources, such as supernova remnants, galaxies, galaxy groups, and clusters. To remain conservative,  \nFig. 2. Distribution of the selected filaments in the redshift and physical length space.  \n<!-- image -->  \nwe did not apply selection cuts to the contaminating X-ray emitting sources.  \nWe used a background level-dependent radius to exclude detected sources as follows. For each source, the surface brightness profile can be approximately described as  \nS ( r ) = S 0       1 + r 2 r 2 PSF + r 2 EXT       -1 . 5 , (1)  \nwhere r PSF = 9 '' is the survey averaged point spread function (PSF) core radius parameter obtained from point source stacking and r EXT is the parameter EXT given by the catalog if the source is extended. We adopted the radius at which the source profile reaches 0.1 (10%) of the background level as the masking radius, namely, S ( r cut) = 0 . 1 × S bkg, where S bkg is obtained from ML\\_BKG\\_1 from the detection pipeline.  \nFor extended sources with large angular sizes, the eSASS detection pipeline usually returns several detections known as split sources (Liu et al. 2022, 2024). These large angular size sources are already included in the eRASS1 galaxy cluster and group catalog (Bulbul et al. 2024). Therefore, we masked out all sources in the eRASS1 galaxy cluster and group catalog with radii of 1 . 5 × r 500. Additionally, we adopted an optical cluster catalog compiled by the redMaPPer algorithm (Ryko ff et al. 2014) on the DESI Legacy Survey (Dey et al. 2019) DR9 + DR10 grz images in a blind mode (Kluge et al. 2024). We excluded sources with a richness of λ > 20 and redshift of 0 . 2 < z < 0 . 8 at radii of 1 . 5 × r λ , where r λ = ( λ/ 100) 1 / 5 h -1 Mpc (this is the richness scaled radius from the redMaPPer algorithm). This cut safely removes clusters of galaxies with high completeness and high purity, especially for high-redshift systems whose X-ray luminosity is below the eRASS1 detection limit.  \nTo conservatively remove the emission from low-richness galaxy groups in the local Universe, we also adopted a galaxy group catalog (Tinker 2021) based on the SDSS Main Galaxy Sample (Strauss et al. 2002). We exclude sources with n sat > 1 in radii of 1 . 5 × r 500. If a source displays 1 . 5 × r 500 < 400 kpc, we  \nused a fixed value of 400 kpc instead. In addition to those cluster and group catalogs, we masked out critical points identified by DisPerSE of type maximum and bifurcation in the filament catalog in a redshift range of 0 . 2 < z < 0 . 6 with fixed radii of 2 Mpc. This range is what we consider as the nodes of the largescale structure (see Malavasi et al. 2020, for details). Details on all masked sources are listed in table 1.  \nWe note that we did not mask optically selected clusters with λ < 20, which is the richness range where the identified clusters or groups are usually seen to su ff er from significantly high contamination due to projection e ff ects (Costanzi et al. 2019; Grandis et al. 2021). Masking these sources out will significantly reduce the useful solid angle in the analysis footprint. Given that the λ < 20 halos and the unmasked sources below the X-ray detection limit also contribute to the stacked signal in this work (due to their associations with the large-scale structure of the Universe), we estimated the fraction of these unmasked X-ray sources and model it out from our final analysis when determining the significance of the emission from the WHIM (see Sect. 4).", '2.2. Healpix map creation': "We adopted the 0.3-1.2 keV band in the imaging analysis. This band is optimized for maximizing the WHIM emission from our selected filament sample with z ∼ 0 . 47 (redshifted), assuming a kT < 1 keV gas temperature. We split the analysis footprint into 27 tiles, each size of 10 deg × 10 deg. For each tile, we merged the event list using evtool and obtain flare-free time interval using flaregti with parameters pimin=5000 , threshold=1.1 , and source\\_size=150 . We used evtool and criteria PATTERN=15 , FLAG=0xE000F000 , and GTI=FLAREGTI to select events in the flare-free time interval and generate count images with a pixel size of 16 '' . We used expmap to create a corresponding vignetting corrected exposure map for each count image. After masking the sources, we binned the tile image pixels to NSIDE = 1024 HEALPix 2 (Górski et al. 2005; Zonca et al. 2019) map pixels. The NSIDE = 1024 HEALPix has a 3 . 4 ' pixel size corresponding to 1.4 Mpc at z = 0 . 6. For the count map, the count number of the i th pixel is C HP , i = P j Cij , where Cij is the count number of the j th image pixel that is located in the i th HEALPix pixel and is not masked. For the exposure map, the exposure time of the i th pixel, t HP , i = ( P j ti j × Ω i j ) / Ω HP , i , where P j ti j and Ω i j are the exposure and solid angle of the j th image pixel that is located in the i th HEALPix pixel and is not masked, and Ω HP , i is the solid angle of the i th HEALPix pixel.", '3. Surface brightness stacking': 'We use i , j , and k to denote the indices of pixel, filament, and radial transverse distance bin, respectively, and we have n HP HEALPix pixels, n fil filaments and n bin distance bins. For each combination of the i th pixel P i , the j th filament F j and the k th distance bin B k , we assigned a Boolean coe ffi cient bijk , whose value is 1 (or 0) if P i is in (or not in) the k th distance bin of F j . Therefore, for each filament, F j , the count number in the k th transverse distance bin is Cjk = P n HP i bijk × Ci , and the exposure time in the k th transverse distance bin is t jk = P n HP i bijk × ti . We used exposure-weighted averages to obtain the stacked surface brightness profile. The stacked profile,  \nTable 1. Summary of the masked sources.', '3.1. Stacking the surface brightness profile': "S X( rk ) = 1 P n fil j t jk n fil X j t jk × Cjk t jk = P n fil j Cjk P n fil j t jk , (2)  \nis the quotient between the sum of the count profile and the sum of the exposure profile, where rk is the radial transverse distance of the k th bin. We used the averaged surface brightness in the outer part of the stacked profile as the local background, S X , bkg, and subtracted it from the entire profile to get a net surface brightness profile as  \nS X , net( rk ) = S X( rk ) -S X , bkg . (3)  \nDue to the projection e ff ects, the distance bins of di ff erent filaments at di ff erent redshifts could be simultaneously located at a single position. Therefore, the surface brightness in di ff erent bins is possibly correlated. We applied a bootstrapping method to construct resampled pixel sets to obtain the estimates of the mean and covariance of the stacked profile. We used l to denote the index of the randomized pixel sample. The l th resample of the pixel set with the replacement is  \nP ' l = GLYPH<16> P ' l 1 P ' l 2 . . . P ' ln HP GLYPH<17> , with P ' li ∈ P . (4)  \nFollowing Eqs. 2 and 3, for each resampled pixel set P ' l , we obtained a stacked surface brightness profile S ' X , net , l ( rk ). After drawing n bs random resamples, the bootstrapping estimate of the net surface brightness profile is  \nˆ S X , net( rk ) = 1 n bs -1 n bs X l S ' X , net , l ( rk ) , (5)  \nand the bootstrapping estimate of the covariance between the k 1 th and k 2 th bin is  \nC k 1 k 2 = 1 n bs -1 × n bs X l h S ' X , net , l ( k 1) -ˆ S X , net( k 1) i h S ' X , net , l ( k 2) -ˆ S X , net( k 2) i . (6)  \nThe uncertainty of ˆ S X , net , k can be taken from the diagonal element of the covariance matrix ˆ σ 2 S X , k = C kk .  \nWe applied our stacking method to the HEALPix count and exposure maps we created (as described in Sect. 2.2). Figure 3 shows the stacked surface brightness against the transverse distance from the filament spine in the 0.3-1.2 keV band. The profile is similar to those in Tanimura et al. (2020b, 2022), exhibiting a centrally peaked signal, and the intensity gradually decreases to ∼ 10 Mpc in the transverse distance. We used the 10-20 Mpc distance region to evaluate the local background level and subtracted it from the profile to obtain the net surface brightness profile of the filaments. To evaluate the significance of detecting the source signal by accounting for the magnitude of fluctuations and systematic e ff ects due to varying background levels across the sky, we constructed 500 control samples of filaments. Each control sample contains 7817 filaments, but we moved them to random positions in the analysis footprint and then rotated them by random angles. The mean surface brightness profile and the corresponding standard deviation of these 500 control samples are also plotted in Fig. 3. The excess Xray emission from the source filaments over the control sample (i.e.,background in this case) is significant, reaching a 9 . 0 σ confidence level in the 0-10 Mpc bins.  \nTo investigate the dependence of the significance of the stacked profile on the selection of the energy band and filament lengths, we applied our analysis to two sub-bands and two distinct filament populations of di ff erent lengths and found consistent results with our analysis with the default selection (see Appendix A for details). The signals in the 0.3-0.7 keV and the 0.7-1.2 keV bands are significant. Meanwhile, the detected filament width of the 20-40 Mpc lengths is marginally broader than that of the 40-100 Mpc length.  \nFurthermore, we investigated the impact of di ff erent source mask radii on the stacked profile (see Appendix B for details). If we increase the cluster mask radii by a factor of two, the resulting profile agrees with the original one within the uncertainties and the integrated signal within the 0-10 Mpc transverse distance bin is only reduced by 7% (see the left panel of Fig. B.1). If we adopt a more conservative mask radius criteria for eRASS:4 X-ray sources, which is S ( r cut) = 0 . 05 S bkg, the resulting profile agrees with the original one within 1 σ uncertainties as well and the integrated signal is enhanced by 3% (see the right panel of Fig. B.1). Therefore, we conclude that increasing the source mask sizes has a trivial impact on the net profile normalization, which is in agreement with the inspections presented in Tanimura et al. (2022).", '3.2. Stacking method validation': "Filaments at di ff erent redshifts could be overlapped with each other due to projection e ff ects. The average number of projected filaments on the pixels within 4 Mpc transverse distance for our selected filament sample is 2.5. This number increases to 5.5 and 11 for pixels within transverse distances of 10 Mpc and 20 Mpc,  \nFig. 3. Stacked filament surface brightness profile in the 0.3-1.2 keV band at the observer's frame. The X-axis is the transverse distance from the filament spine. The orange and purple lines and shaded regions denote the profiles of the source and the averaged value of the randompositioned control samples, respectively. The average surface brightness in the 10-20 Mpc transverse distance range is used as the local background and is subtracted from the profile to obtain the surface brightness signal of the net source profile. The shaded area of the source profile is the 1 σ uncertainty estimated from a pixel-wise bootstrapping method, and that of the control samples is the 1 σ scatter. The correlation matrix of the five bins in the 0-10 Mpc range is plotted in the inset figure.  \n<!-- image -->  \nrespectively. There are two di ff erent cases of overlapping. In the first case, the two overlapped filaments have di ff erent orientations (left in Fig. 4). When extracting the profile of filament A, the contribution from filament B in the overlapped source region will be accounted for by the filament B emission in the local background region and vice versa. In the second case, the two overlapped filaments have similar orientations (right in Fig. 4), while the emission in the overlapping source region will be double counted; therefore increasing the stacked signal. Meanwhile, our stacking method assumes that systematic uncertainties from the spatial variations of multiple components, for instance, the Milky Way foreground emission and absorption,cosmic X-ray background (CXB) emission. In addition, the inhomogeneity of exposure is minimized by a large number of filaments randomly distributed in the footprint.  \nTo evaluate the impact of the second case overlapping and multiple possible systematic uncertainties, we first simulated the mock count maps and applied the same surface brightness profile stacking method to the mock maps to obtain the averaged stacked result and to compare it with the injected model. The details of the mock map creation are presented in Appendix C. Wesimulated 200 mock count maps and applied the same profile stacking method. The stacked mock profiles from the mock maps are plotted in Fig. 5, together with the injected model profile. The averaged mock profile is in general agreement with the injected model within the scatter; whereas, it is slightly higher than the injected model by 13%. This 13% enhanced signal could be due to the second case of overlapping and we corrected this e ff ect when modeling the gas properties (described in the subsequent sections). The validation also suggests that although the averaged filament signal at the spine is 0 . 5% of the total background  \n<!-- image -->  \nFig. 4. Sketches of two types of overlapped filaments due to the projection e ff ects. Deep and shallow regions are the source and local background regions, respectively. In the first case ( left ), the two filaments have di ff erent orientations, and the emission in the overlapping region will be accounted for by local background subtraction. In the second case ( right ), the two filaments have similar orientations and will lead to overestimated stacked results, evaluated in Sect. 3.2.  \n<!-- image -->  \nFig. 5. Results of 200 net surface brightness profiles of the simulated maps obtained by applying our stacking method (in gray). The red dashed curve is the injected surface brightness profile model. The curve in purple and the shaded region are the mean and standard deviation of the 200 profiles. Although the filament signal at the spine is only ∼ 0 . 5% of the background value, the stacking method can recover the injected profile successfully.  \n<!-- image -->  \nlevel, our stacking method can successfully recover the emission profile.", '4. Modeling the contamination from unresolved extragalactic sources': "The source masks applied in this work successfully mask out all X-ray-detected point sources and galaxy clusters and groups above the eRASS:4 limit, as well as optically selected galaxy clusters and groups with λ > 20. We note that X-ray sources outside these selection criteria (i.e., those below the eRASS:4 detection limit and low-richness clusters and groups) may also contribute to the stacked X-ray signal. To remove the contaminating signal, we evaluated the contribution from three types of unmasked X-ray signal associated with galaxies: gravitationally bound hot gas in galaxy or galaxy group halos, undetected low-luminosity active galactic nuclei (AGNs), and X-ray bina-  \nies (XRBs). The steps required to model the contribution from these sources in the stacked total signal are given below.  \n- 1. We built a redshift-dependent L X -M ∗ scaling relation of each type of the sources in a 9 . 5 < log( M ∗ / M ⊙ ) < 12 stellar mass range. We describe the contribution from each type of unmasked source in detail in Sect. 4.1.\n- 2. We adopted the DESI Legacy Survey DR9 galaxy photometric redshift and stellar mass catalog (Zou et al. 2019), which is relatively complete and covers our analysis footprint. We selected galaxies in the redshift range of [0.15, 0.7] and cleaned up this catalog by removing sources with a redshift uncertainty of > 20%. We predicted the expected count rate in the 0.3-1.2 keV band for each galaxy using the constructed scaling relations, given the stellar mass, redshift, spectral shape of the source, and the Galactic absorption at the galaxy's position. We further corrected for the incompleteness of the observed stellar mass function (SMF) from Zou et al. (2019) catalog to a fiducial empirical SMF model from the result of UniverseMachine (Behroozi et al. 2019), which is a project characterizing the link between galaxy growth and halo growth. In the M ∗ < 10 11 M ⊙ range, where the galaxies in Zou et al. (2019) catalog are not complete, the ratio between the SMF of the selected galaxies from Zou et al. (2019) and the fiducial SMFs at di ff erent stellar masses and redshifts are used to correct the predicted count rate of each galaxy.\n- 3. For each selected galaxy, we added the predicted count rate to the HEALPix pixel that the galaxy is located in. Thanks to this method, we were able to create a predicted count rate map of unmasked galaxy X-ray emission in the selected redshift range. We applied the same source mask used in the profile extracting and stacking pipeline to extract the profile of unmasked signals.", '4.1. Properties of contaminating unmasked sources': 'To estimate the flux distribution of the contaminating signal, we constructed L X -M ∗ for three types of unmasked X-ray sources. Meanwhile, we also investigated the spectral properties to be used in count rate conversion in spectral fitting. The details of the modeling for each component are given in the subsections below.', '4.1.1. Unresolved haloes': 'We used the X-ray halo L X -M ∗ relation from RASS galaxy stacking analysis (Anderson et al. 2015) to predict the X-ray luminosity of bound gas in unmasked galaxy (group and cluster) halos. Because we excluded eRASS:4-detected X-ray sources and clusters detected in the optical band with a richness higher than 20 before stacking, to obtain the L X -M ∗ relation that meets the same selection - we use our galaxy cluster mock catalog that is constructed for eRASS predictions. The catalog is based on HugeMDPL dark matter-only N-body simulation (Klypin et al. 2016), while the X-ray properties of each halo are illustrated by a phenomenological modeling approach (Comparat et al. 2020) and then scaled to agree with the observed L X -M ∗ relation from RASS galaxy stacking results of Anderson et al. (2015). We tabulated the averaged rest-frame L 0 . 5 -2keV -M ∗ relation for unmasked X-ray halos in di ff erent redshift bins by applying the same source selection criteria.  \nThe galaxy sample used in Anderson et al. (2015) was selected from an SDSS-DR7 based New York University Value  \nAdded Galaxy Catalogue (NYU-VAGC) (Blanton et al. 2005). The stellar mass in NYU-VAGC is calculated based on the method of Blanton & Roweis (2007). It is known that the Blanton & Roweis (2007) stellar mass is systematically lower than many other estimates at the high stellar mass end (e.g., Li & White 2009; Bernardi et al. 2010; Moustakas et al. 2013). Because the X-ray halo L X -M ∗ relation has a steep slope of ≳ 3, a small underestimate of the stellar mass at the high-mass end could introduce a huge di ff erence on the estimated unmasked Xray halo emission. Therefore, we need to properly modify the current scaling relation to match the Zou et al. (2019) stellar mass scheme. We selected a z < 0 . 2 subsample of Zou et al. (2019) galaxies with available spectroscopic redshift to crossmatch the NYU-VAGC catalog and obtained the discrepancy between the two stellar mass schemes. We used this stellar mass discrepancy to empirically modify the tabulated L 0 . 5 -2keV -M ∗ relation. The modified X-ray halo L X -M ∗ relation is plotted in the top left panel of Fig. 6).  \nTo convert the rest-frame 0.5-2 keV luminosity to the 0.31.2 keV observed count rate, the spectral properties of the sources are needed. For this purpose, we assumed that each source can be represented by a single temperature collisional ionization equilibrium (CIE) model and we estimated the temperature using the luminosity - temperature ( L X -kT ) scaling relation constructed using eFEDS galaxy cluster and groups (Bahar et al. 2022).  \nWith the estimated temperature of each source in the mock catalog by using the L X -kT scaling relation, we were also able to predict a synthesis spectral model of unmasked X-ray halos for spectral analysis. After applying the selection, we construct the multi-temperature di ff erential emission measure (DEM) models for each redshift bin (left panel of Fig. D.1). The spectral model on each single temperature is calculated using Astrophysical Plasma Emission Code (APEC) with Z = 0 . 4 Z ⊙ . Then, we used the redshift and length distributions of the selected filament sample to create the synthesis spectral model as:  \nF ( E ) = X i li × F ( E , zi ) D L( zi ) 2 , (7)  \nwhere li is the physical length of the i -th filament, zi is the redshift of the i -th filament, F ( E , zi ) is the multi-temperature di ff erential emission measure spectrum at the redshift of zi , and D L is the luminosity distance as a function of redshift. The constructed model spectrum is shown in the right panel of Fig. D.1.', '4.1.2. Unresolved AGN': 'We calculate the galaxy averaged L X -M ∗ relation for unresolved AGNs using an AGN mock catalog based on the highresolution Uchuu dark matter only N-body simulation (Ishiyama et al. 2021). The AGN properties are painted to dark matter halos based on an empirical method (see Comparat et al. 2019, for further details), where the AGN model is based on the observed AGN 2-10 keV X-ray luminosity function (Aird et al. 2015), duty cycle (Georgakakis et al. 2017), obscuration distribution (Ricci et al. 2017), spectral model (Liu et al. 2017), and parameters of AGN halo occupation distribution (Comparat et al. 2023).  \nWe applied a selection criteria that removes AGNs with F 0 . 5 -2keV > 10 -14 erg s -1 cm -2 , which would be detected in eRASS:4, to select sources in the mock catalog. In each redshiftstellar mass bin, we summed the rest-frame 0.5-2 keV luminosity and divide it by the total number of galaxies in the same bin to  \nFig. 6. Rest-frame galaxy averaged 0.5-2 keV luminosity - stellar mass scaling relations for unmasked X-ray halos ( top-left ), AGNs ( topright ), and XRBs ( bottom-left ) at di ff erent redshifts. The scaling relations of total unmasked X-ray signals associated with galaxies are plotted ( bottom right ) together with the recent eRASS stacking results of z < 0 . 2 central galaxy sample in black.  \n<!-- image -->  \nobtain the galaxy averaged unresolved AGN luminosity. In our method, a galaxy-averaged scaling relation would distribute the unmasked AGN flux to all galaxies. Therefore, we do not need to additionally predict whether a galaxy hosts an AGN. The scaling relations at di ff erent redshifts are shown in the top right panel of Fig. 6.  \nWe also investigated the stacked spectral properties of unresolved AGNs. We adopted the same spectral model used to create the mock catalog. The detailed components in XSPEC (Arnaud 1996) are plcabs + zgauss + constant * powerlaw + pexrav * constant , representing a power-law model with intrinsic absorption, the 6.4 keV Fe K α emission line, a soft-X-rayexcess, and a cold-reflection component, respectively. We predicted a corresponding mock spectrum based on the hard band luminosity and the intrinsic absorbing column density for each unmasked source. The stacked spectra of di ff erent bins and the median spectrum are plotted in Fig. D.2. They show a strong self-similarity in the soft band. Therefore, we use the median spectrum to convert the stacked rest-frame 0.5-2 keV luminosity to 0.3-1.2 keV count rate.', '4.1.3. XRB in galaxies': 'We used an integrated X-ray luminosity model that is a function of stellar mass and star forming rate (SFR) to construct the averaged L X -M ∗ relations of XRB at di ff erent redshifts. The model is constrained by Chandra observations of local galaxies (Lehmer et al. 2019), and the UniverseMachine model of redshift dependent SMF and SFR (Behroozi et al. 2019). The results are plotted in the bottom left panel of Fig. 6, where the relations show redshift evolution at M ∗ < 10 11 M ⊙ (the result  \nof the redshift evolution of star-forming rate and star-forming galaxy fraction).  \nWe adopted a power-law distribution with intrinsic absorption spectral model for flux conversion, with parameters of Γ = 1 . 7 and n H = 2 × 10 21 cm -2 , as suggested by observations (Lehmer et al. 2019).', '4.1.4. Comparison with observed L X -M ∗ relations': 'Zhang et al. (2024) investigated the L X -M ∗ scaling relation of galaxies using eRASS:4 data. We compared our constructed total unresolved galactic X-ray sources and their L X -M ∗ relations with their results from eROSITA observations. We note that the stellar mass used in Zhang et al. (2024) was adopted from Tinker (2021), which was originally estimated by Chen et al. (2012). Similarly to the analysis in Sect. 4.1.1, we scaled the observed L X -M ∗ scaling relation in Zhang et al. (2024) to match the stellar mass scheme of Zou et al. (2019) for comparison. In the bottom right panel of Fig. 6, we plot the total unresolved source scaling relation model together with that of a central galaxy sample from Zhang et al. (2024), whose median redshift is ∼ 0 . 1. Our constructed relation at z = 0 . 2 has the same slope at di ff erent stellar masses as the observed relation and the slight luminosity di ff erence meets the expected redshift revolution of AGN and XRB scaling relations.', '4.2. Fractions in stacked signals': 'Next, we calculated the integrated signal of unmasked halos, AGNs, and XRBs in the 0-10 Mpc transverse distance range to account for the contamination in our imaging analysis. The total  \ncontribution is 37%, with each of 20%, 12%, and 5% for unmasked halos, AGNs, and XRBs, respectively (also see Fig. 7). We further divided the stellar mass range into two subranges (9 . 5 < log M ∗ ( M ⊙ ) < 11 and 11 < log M ∗ ( M ⊙ ) < 12) and repeated the same analysis. The results of the fractions are listed in Table 2. In the high stellar mass bin, where the unmasked galaxies are mostly central galaxies of group-size halos, the unmasked signal is dominated by halos gas; whereas in the low stellar mass bin, unmasked AGNs and XRBs demonstrate more of a contribution. We used a bootstrapping method to estimate the statistical uncertainties of the fraction. In short, we randomly select galaxies with replacements to construct bootstrapping samples and calculated the standard deviations of the fractions estimated using the bootstrapping samples as the statistical uncertainties. The relative uncertainties are 6 . 9%, 1 . 8%, and 1 . 4% for unmasked halos, AGNs, and XRB, respectively. Given that the unmasked halos contribute about half of the unmasked signal, the averaged statistical uncertainty is below 5%.  \nBecause the uncertainty on photometric redshifts that we used in this analysis is not negligible and is a function of magnitude (it can reach up to a few tens of percent for faint galaxies Zou et al. 2019), it has a potential impact on the recovered unresolved contamination signal. The redshift range of the filament sample we select is 0 . 2 < z < 0 . 6, while the range for the galaxies we used to model unmasked source contribution is 0 . 15 < z < 0 . 7. This means that for a galaxy whose true redshift z true is within the filament redshift range, but the measured redshift z phot is outside the redshift range of the selected galaxies, it may not be taken into account in the unmasked source profile. On the contrary, for a galaxy with z true outside of the filament redshift range but z phot within that range, because it is not spatially correlated with the filament we selected, it may not contribute to the unmasked source profile either. Therefore, we note that the fraction of unmasked source signals could be underestimated in our analysis due to the galaxy photo-z uncertainties.  \nWe use a Monte Carlo method to calculate the recovered flux fraction of galaxies at a given redshift with given redshift uncertainties. For creating the randomized redshifts, we follow the definition of normalized redshift error used by Zou et al. (2019), ∆ z norm ≡ ( z phot -z spec) / (1 + z spec) and assume that z spec = z true. The ratio between the recovered flux and the true flux of a galaxy is  \nf = D 2 L ( z true) n MC X i wi D 2 L ( zi ) , (8)  \nwhere n MC is the number of realizations, and zi is the i th randomized redshift value, while the value of wi is 1 if 0 . 15 < zi < 0 . 7 and it is 0 otherwise. The results are shown as dashed lines in Fig. 8. Based on the 90% percentile of the r -band magnitude distribution of the galaxies we selected and the r -∆ z relation shown in Fig. 5 of Zou et al. (2019), we were able to predict the recovered flux fractions for unmasked galaxies in the 9 . 5 < log M ∗ ( M ⊙ ) < 11 and 11 < log M ∗ ( M ⊙ ) < 12 stellar mass ranges, which are shown in the crossed and triangle curve in Fig. 8, respectively. The recovered fraction of galaxies in the high stellar bin is around one due to the high quality of photometric redshift of bright central galaxies. This fraction is ∼ 80% in the low stellar mass bin when z > 0 . 3. It means that even if we take the impact of photo-z uncertainty into account, the total contribution from unresolved sources can be as high as ∼ 40%.  \nThe ∼ 40% fraction of the unmasked source contribution suggests that the remaining ∼ 60% of the stacked signal is from the WHIM in the filaments, reaching a 5 . 4 σ detection significance. This fraction relies on the galaxy X-ray contamination  \nFig. 7. Contribution of unmasked X-ray halos (green), AGNs (light blue), and XRBs (grey) to the total stacked signal (orange). To find the net WHIM signal, the contribution from these contaminating signals is modeled out in the analysis. The signal from WHIM is ∼ 60% based on the unmasked signal estimation.  \n<!-- image -->  \nTransverse distance to spine (Mpc)  \nTable 2. Fraction of X-ray signal from unmasked sources in di ff erent stellar mass bins.  \nmodeling, which is based on our current knowledge of di ff erent X-ray emitters in galaxy systems at di ff erent scales. A better understanding of all the di ff erent types of X-ray sources in the Universe will o ff er more accurate results of the contamination fraction in the future.', '5. Spectral stacking': 'In this section, we describe the X-ray spectral analysis of the filaments. Due to the faint nature of the cosmic filaments, to increase the signal-to-noise ratio (S / N) in the data, we stacked the spectra extracted from individual filaments by following the recipe provided in Bulbul et al. (2014). The method successfully produces the stacked rest-frame spectrum of filaments at di ff erent redshifts. Through this analysis, we study the spectral properties and constrain the temperature of the detected WHIM using the stacked filament spectra.', '5.1. Stacking processes': 'For each filament, we first merged the event files of TM8 and shifted the pulse-invariant of each event by a factor of 1 + z . Then, we used srctool to extract the source and local background spectra and the vignetting corrected source ancillary response file (ARF). The source and local background extraction regions are defined as 0-10 Mpc and 10-20 Mpc transverse distance to the filament spine, whose region masks for spectral extraction are created by constructing images using and excluding sources in Table 1. For each ARF, we first corrected it by multiplying the foreground absorption curve given the averaged  \nFig. 8. Fraction of recovered flux of unresolved sources as a function of true redshift (x-axis) and normalized redshift uncertainty (dashed lines in color). The triangle and crossed lines are the fractions for our selected galaxies in the 11 < log M ∗ ( M ⊙ ) < 12 and 9 . 5 < log M ∗ ( M ⊙ ) < 11 stellar mass bins, respectively. The impact of photometric redshift uncertainty is only a few percent for high stellar mass galaxies that contribute most of the unmasked signals.  \n<!-- image -->  \ntrue  \nGalactic n H value (HI4PI Collaboration et al. 2016) in the extraction region and the TBABS absorption model (Wilms et al. 2000). We then shifted it by a factor of 1 + z to the source rest frame. The response matrix file (RMF) was taken from the calibration database and then shifted to the source rest frame by adopting a two-dimensional (2D) interpolation strategy on both axes of induced photon energy and resulting spectral channels (see Bulbul et al. 2014, for further details).  \nWe used tools addspec , addarf , and addrmf in FTOOLS 3 (Nasa High Energy Astrophysics Science Archive Research Center (Heasarc) 2014) to stack the spectra, ARFs, and RMFs, respectively. We note that some keywords in the headers of srctool products have di ff erent definitions from those in other X-ray satellite data products 4 . In short, the EXPOSURE value of a spectrum t EXP is the total time of an extraction region on a detector, instead of the averaged exposure time in the region t ave; BACKSCAL value Ω BS is the averaged intersection of the extraction region and the detector area; the geometric solid angle of the extraction region Ω REG has a keyword of REGAREA, along with the ratio Ω REG / Ω BS = t EXP / t ave. Furthermore, the normalization of an ARF curve is scaled by a factor of t ave / t EXP. Therefore, we need to add up all individual source spectra to obtain a stacked source spectrum, as follows:  \nN stacked , src( E ) = X i N src , i ( E ) . (9)  \nWe then scale each local background spectrum to match the extraction solid angle of the corresponding source spectrum to get the stacked background spectrum,  \nN stacked , bkg( E ) = X i Ω REG , src , i Ω REG , bkg , i × N bkg , i ( E ) . (10)  \nFor the final stacked source spectrum, we assign the exposure time as a solid angle weighted averaged exposure t ave = GLYPH<0>P i Ω REG , i × t EXP , i GLYPH<1> / P i Ω REG , i . To match this exposure, the stacked ARF is defined as  \nA stacked( E ) = 1 t ave X i t EXP , i Ω REG , i Ai ( E ) , (11)  \nwhich takes each extraction solid angle into account. The modeled source flux using the exposure and stacked ARF is scaled per deg 2 . Similarly, the stacked RMF is weighted by Ω REG × t ave, which is  \nR stacked = 1 P i Ω BS , it EXP , i × X i Ω BS , it EXP , i R i . (12)  \nThe blueshifted response files at the rest frame corresponding to each filament are later stacked using FTOOLS addarf and addrmf to obtain one final stacked response file. These files are used to model the stacked spectrum described in the section below.  \nOne clear advantage of the spectral stacking method in the rest frame is reducing the e ff ect on the instrumental and background features on the final net background-subtracted source spectrum as described in (Bulbul et al. 2014). The energies of the instrumental and celestial X-ray foreground and background features, frozen in the detector frame, are smoothed over the redshift range covered by the sample, reducing their overall e ff ect on the net spectrum. The spectral stacking method is well suited for obtaining a spectrum of faint sources, such as cosmic filaments.', '5.2. Stacked spectrum and spectral fitting': 'We stacked the 0-10 Mpc transverse distance region spectra as the source spectrum and the 10-20 Mpc transverse distance region spectra as the local background region spectrum. The stacked background spectrum is smoothed using a Wiener filter with a kernel size of five channels for direct subtraction. Additionally, we compare the 6-9 keV count rate of the two spectra and store the ratio as the backscal in the header for further correction. The stacked net source spectrum is plotted in Fig. 9.  \nWe used XSPEC version 12.13.0 and the BXA 5 package (Buchner et al. 2014) for our spectral analysis. The spectrum was first optimally binned (Kaastra & Bleeker 2016) using the tool ftgrouppha and then loaded into XSPEC. We used the 0.56.0 keV range for analysis and adopted C-stat (Kaastra 2017) to derive the fitting likelihood. The BXA package uses the nested sampling algorithm UltraNest (Buchner 2016, 2019, 2021) for Bayesian inference.  \nWe fit the spectrum using four components that account for the detected WHIM, unmasked halo gas, unmasked AGN, and unmasked XRB. The spectral model of the three unresolved sources is described in Sect. 4.1. We calculated the expected normalization ratios among the three components based on the modeled unresolved source fractions in Sect. 4.2 and coupled the ratios during fitting. To model the WHIM component, we used a single temperature (1-T) APEC model with fixed abundance Z = 0 . 2 Z ⊙ . This value is based on the reported ≲ 0 . 3 Z ⊙ abundance in cluster outskirts by both observations and numerical simulations (Bulbul et al. 2016; Mernier et al. 2017; Bi ffi et al. 2017) and > 0 . 1 Z ⊙ in z < 1 filaments from numerical simulations (Martizzi  \n<!-- image -->  \nFig. 9. Rest frame stacked source spectrum of the 0-10 Mpc region from the filament spine in the transverse direction after the subtraction of the local background extracted in the 10-20 Mpc region surrounding the filament. For illustration, the spectrum is binned to a minimum signalto-noise ratio of 2.5. The left panel shows that the combination of a single temperature CIE (red dashed line), and the unmasked sources models (green, blue, and grey) can well fit most of the spectrum, except the significant excess at E < 0 . 6 keV. The best-fit temperature of the WHIM gas is 0 . 58 ± 0 . 10 keV. The right panel shows a fitting by a log-normal temperature distribution with log( µ/ K) = 5 . 9 and σ = 0 . 45 dex. The 0.5-0.6 keV energy range can be well fitted by this model, though there is a likely deficit at ∼ 0 . 7 keV.  \n<!-- image -->  \net al. 2019). This model combination can well fit the E > 0 . 7 keV spectrum range with a best-fit parameter kT = 0 . 58 ± 0 . 10 keV and a fitting statistics c stat / d . o . f . = 105 / 111. We tried to fix the metalicity to 0 . 1 Z ⊙ and 0 . 3 Z ⊙ to test the impact of the metalicity assumption and found there is no impact on the best-fit temperature. In addition, we also fitted a stacked 0-4 Mpc transverse distance region spectrum and obtained a similar 0 . 56 ± 0 . 09 keV temperature.  \nRecently, Vijayan & Li (2022) demonstrated that fitting a multi-temperature nature spectrum using a 1-T model will lead to an overestimation in the temperature. It explains that our bestfit temperature using a 1-T model does not agree with the majority of the gas phases around 10 6 Kpredicted by simulations (e.g., Martizzi et al. 2019; Galárraga-Espinosa et al. 2021). In our stacked spectrum, we observe an excess feature in 0.5-0.6 keV, which agrees with the energy range of ∼ 0 . 56 keV O VII lines emitted by cooler phases of gas around 10 6 . 3 K ( ∼ 0 . 17 keV), suggesting the presence of cooler phases. The gas in these temperature phases only contributes to lines below 0.7 keV, while the fitting using a 1-T model tends to maximize the likelihood around the peak of the observed spectrum, which is in the 0.80.9 keV energy range, dominated by several Fe XVII lines emitted by 10 6 . 8 K gas.  \nWe also tried to fit the spectrum using a multi-temperature model. We adopted a log-normal temperature distribution DEM model, with a peak at log( T / K) = 5 . 9 and 1 σ width of 0.45 dex (as suggested by our analysis of numerical simulations; see more in Sect. 6.2). We used a simplified assumption of CIE condition and 0 . 2 Z ⊙ abundance of the WHIM. The log-normal DEM model can better fit the stacked spectrum in the 0.5-0.6 keV energy range (see the bottom panel of Fig. 9), though this fitting shows the feature of a likely deficit around 0.7 keV. This could be because the true temperature distribution is more complex than the log-normal distribution we use here, which creates discrepancies in some energy ranges.', '6.1. Surface brightness profile modeling': 'With the best-fit temperature obtained in the previous section, we fit the density profiles of the gas in cosmic filaments. In this section, we use r to denote the physical radius of gas cells or particles to the filament spine in the 3D space. Based on results from numerical simulations (Cui et al. 2019; Galárraga-Espinosa et al. 2021; Tuominen et al. 2021; Zhu et al. 2021), the gas in filaments has a multi-phase nature, ranging from 10 3 to 10 7 K in temperature and from 10 -6 to 10 -1 cm -3 in hydrogen density. Therefore, we introduced the parameter, f vol, to represent the fraction of the volume of the detected gas to the total filament volume. The volume-averaged emissivity of the detected phase gas at a certain radius is  \n¯ ϵ det( r ) = 1 . 20 × n H( r ) 2 f vol( r ) , (13)  \nwhere n H is the hydrogen density, 1.20 is the n e / n H ratio for the plasma at the detected temperature. We used a King-type profile (King 1962), which is widely used to fit central flattened profiles, to model the radial change of both n H and f vol,  \nn H( r ) = n H , 0 h 1 + ( r / r c , n) 2 i α n ; (14)  \nf vol( r ) = f vol , 0 h 1 + ( r / r c , fvol) 2 i α fvol , (15)  \nwhere n H , 0 ( f vol , 0), r c , n ( r c , fvol), and α n ( α fvol) are the central normalization, core radius, and slope in the outskirts for n H (or f vol), respectively. Because the radial dependencies of n H and f vol have a strong degeneracy, we assume that they have the same core radii, namely, r c , n = r c , fvol.  \nIn addition, we must make several assumptions. First, we assume a constant baryon density contrast ∆ b , 0 at the filament spine for filaments at all redshifts. In this case,  \nn H , 0( z ) = ∆ b , 0 ρ crit Ω b 1 . 43 × m H × (1 + z ) 3 , (16)  \nTable 3. Parameters and corresponding priors to constrain the density profile of the observed phase WHIM.  \nwhere ρ crit is the critical density at z = 0, the factor 1.43 is the ratio between total mass density and hydrogen mass density, m H is the hydrogen mass, and the redshift evolution comes from the expansion of the physical scale. Second, there is an outer boundary of the filaments r out, beyond which the gas density of our detected phase is zero. In this work, we have adopted a value of 10 Mpc. Third, the observed profile is smoothed and broadened due to the position uncertainty of the skeleton constructed from incomplete galaxy samples (e.g., Fig. 11 of Galárraga-Espinosa et al. 2021). Including these assumptions, the projected surface brightness at a given projected radial distance, r proj, is then  \nS X , model , conv( r proj , z ) = G ( σ sm) ∗ S X , model( r proj , z ) , (17)  \nwhere G ( σ sm) is a Gaussian convolving kernel that accounts for the smoothing e ff ect of the observed filament width. The surface brightness is  \nS X , model( r proj , z ) = A -1 CF K ph(1 + z ) -3 × Z l out -l out ¯ ϵ det GLYPH<18> q l 2 + r 2 proj , z GLYPH<19> Λ d l , (18)  \nwhere l is depth along the line of sight; K ph is the K-correction of photon spectrum given the spectral shape, redshift, and energy band; the term (1 + z ) -3 is the photon dimming factor; Λ is the cooling function calculated using the APEC model given the temperature, the assumption of 0 . 2 Z ⊙ abundance and the energy band; A CF is the conversion factor from photon flux to count rate, which is calculated given the on-axis ARF, observation band, observed photon spectrum, and Galactic absorption curve.  \nWhen computing the fitting statistics, we constructed 7817 individual model surface brightness count rate profiles based on the redshift of each filament. We use the solid angle of each filament as the weight to obtain an averaged surface brightness profile as  \nS X , model , stack = 1 N N X i Ω i × S X , model , conv , i . (19)  \nThe fitting statistics is defined as  \nχ 2 = S X , di ff ⊺ C -1 S X , di ff , (20)  \nwhere S X , di ff = S X , net , corr -S X , unmasked , net -S X , model , stack × 1 . 13, is the di ff erence between the model profile and observed net profile with the unmasked sources contamination subtracted. The factor of 1.13 takes into account the possible enhancement of the signal due to the overlapping filaments (see Sect. 3.2). The C is the covariance matrix estimated from the bootstrapping method.  \nWe used a Bayesian inference approach to obtain posterior distributions of the six parameters: log ∆ b , 0, log f vol , 0, log r c , n, α n, α fvol, and log σ sm, as listed in Table 3. In addition to the fitting  \nstatistics χ 2 , the likelihood function includes an additional constraint from SZ stacking results (Tanimura et al. 2020a), which is necessary to break up the degeneracy between ∆ b , 0 and f vol , 0. The predicted Comptony strength is  \ny pd( r proj) = " σ T m e c 2 Z l out -l out kTn e GLYPH<18> q l 2 + r 2 proj GLYPH<19> f vol GLYPH<18> q l 2 + r 2 proj GLYPH<19> d l # , (21)  \nwhere σ T is the Thomson scattering cross section, m e is the electron mass, and c is the speed of light. The profiles of Tanimura et al. (2020a) show similar broadening e ff ects as our stacked X-ray profile. Therefore, we also convolved the predicted Comptony profile with the same Gaussian kernel to create a smoothed Comptony profile as  \ny pd , conv = G ( σ sm) ∗ y pd . (22)  \nGiven that the phases of the WHIM we detect only contribute a fraction of the total stacked SZ signal, we tentatively set the upper limit of the log y pd , conv(0) to be -7.5, which is approximately the stacked SZ signal at filament spine reported by (Tanimura et al. 2020a). We used a likelihood term P (log y pd , conv(0)) = N ( µ logy , σ logy) as the addition constraint. We do not know how much of the total stacked Comptony is contributed by the Xray-emitting gas. Nevertheless, we adopt µ logy = -8 . 3 and σ logy = 0 . 5 in our analysis. The combined likelihood function is  \nln L =          -1 2 χ 2 -1 2 [ y pd , conv(0) -µ logy] 2 σ 2 logy when log y pd , conv(0) < = -7 . 5 , -∞ elsewhere.  \n(23)  \nWe used a Markov chain Monte Carlo (MCMC) method and the python package Zeus 6 (Karamanis et al. 2021) to perform Bayesian inference. The marginalized posterior distributions shows that the log ∆ b , 0 = 2 . 28 ± 0 . 22 and log f vol , 0 = -1 . 91 ± 0 . 28. By using all MCMC chains for calculating the averaged baryon density contrast, we obtained log ¯ ∆ b = 1 . 88 ± 0 . 18, which corresponds to a hydrogen density 4 . 4 + 2 . 2 -1 . 5 × 10 -5 cm -3 at the median redshift of the selected filament sample z = 0 . 47. We note that we assumed a Z = 0 . 2 Z ⊙ metalicity for this inference. If the WHIM has a metalicity of 0 . 1 Z ⊙ (and 0 . 05 Z ⊙ ), at the temperature of 0.6 keV, the radiative cooling function in our analysis band is a factor of 1.5 (or 2.5) smaller than that of the 0 . 2 Z ⊙ assumption, which would introduce a 0.09 (and 0.2) dex increase in the inferred log ∆ b , 0.  \nFig. 10. IllustrisTNG phase diagrams of the filament gas at z = 0 . 48 weighted by mass ( left ), electron pressure ( middle ), and 0.5-2 keV photon luminosity folded by eROSITA response ( right ). Green and white contours are 50% enclosed phases in electron temperature and X-ray emission, respectively. The best-fit temperature using a 1-T model and baryon density contrast of this work are shown in yellow with error bars.  \n<!-- image -->', '6.2. X-ray emitting phases in numerical simulations': 'To compare the measured physical quantities with those in numerical simulations, we used the public data (Nelson et al. 2019) from IllustrisTNG simulations (Nelson et al. 2018; Pillepich et al. 2018; Naiman et al. 2018; Springel et al. 2018; Marinacci et al. 2018). We selected the z = 0 . 48 snapshot of the TNG300-1 box and aimed to create phase diagrams of the filament gas. We used the DisPerSE algorithm to construct the filaments in that snapshot using 9 < log( M ∗ / M ⊙ ) < 12 galaxies (Nelson et al. 2019). The detailed procedure follows Galárraga-Espinosa et al. (2021) and Malavasi et al. (2022). The calculations of the distances of each gas cell to the nearest filament, d fil, to the nearest maxima or bifurcation critical point, d CP, and the nearest critical point but along the direction of filament, d skel, follow the definition in Malavasi et al. (2022).  \nTo reduce the computational time, we randomly selected 1 / 500 of the total gas cells, which is enough to represent the phase diagram (Galárraga-Espinosa et al. 2021). We selected gas cells with d fil < 1 Mpc and d CP > 2 Mpc. Then, we removed cells in halos or within a 25 Mpc range to simulation box boundaries. The temperature of a gas cell is calculated by using the quantities InternalEnergy and ElectronAbundance based on the standard recipe 7 . The baryon density contrast of a cell is converted from the density by comparing it with the mean baryon density at the redshift of 0.48.  \nWe plot the mass-weighted phase diagrams of filament gas in the left panel of Fig. 10, where the majority of the gas is in 10 4 . 5 -10 6 . 5 K temperature and 1 -100 baryon density contrast phases. We note that the presence of hot phases with T > 10 7 K is not as prominent as those in the z = 0 snapshot (GalárragaEspinosa et al. 2021). It could be due to the redshift evolution of filament gas, which is being continuously heated during structure formation (Cui et al. 2019; Martizzi et al. 2019).  \nTo investigate the traceable populations of the filament gas, we created two phase diagrams weighted by electron pressure and 0 . 5 -2 keV photon luminosity folded by eROSITA response, respectively (middle and right panels in Fig. 10). The former one  \nrepresents the traceable population of SZ observations (e.g. de Graa ff et al. 2019; Tanimura et al. 2019, 2020a), while the latter one represents the traceable population of soft band X-ray emission studies. The two diagrams show that the traceable phases of SZ signal and soft band X-ray emission di ff er. SZ observations trace 10 5 . 5 -10 6 . 8 K temperature and a ∆ b = 10 0 . 3 -10 2 density phases WHIM in filaments.On the other hand, X-ray emission in the 0.5-2 keV band traces 10 6 . 2 -10 6 . 9 K temperature and ∆ b = 10 1 . 2 -10 2 . 2 density phases WHIM. It also shows that the X-ray emitting population only contributes a fraction of the total SZ signal, which agrees with the assumption we used in Sect. 6.1. Our averaged density contrast measurement agrees with the TNG numerical simulations. However, the best-fit temperature using a single component model is close to the upper boundary of the contour of X-ray emitting WHIM (see the yellow data point with error bar in the left panel of Fig.10). This discrepancy has been addressed in Sect. 5.2, i.e., the best-fit temperature value using a 1-T model will be biased to the high end of the temperature distribution when fitting a spectrum with multitemperature nature.', '6.3. Impact of photoionization': 'It is known that the WHIM is also partially ionized by the cosmic cosmic ultraviolet background (UVB), which alters the ionization balance and exhibits di ff erent spectral features from the CIE plasma (e.g., Nicastro et al. 1999; Krongold et al. 2003; Khabibullin & Churazov 2019; Churazov et al. 2023). In this section, we describe how we combined the WHIM phase distributions in numerical simulations (left panel of Fig. 10) and the photoionization emission model to characterize the impact of including the photoionization model on our results.  \nThe ionization parameter is defined as  \nξ ≡ 4 π F 1 -1000Ry n H , (24)  \nwhere F 1 -1000Ry is the local flux in the 1-1000 Ry energy band (1 Ry = 13.6 eV) following the convention of Krolik et al. (1981). We adopted the Faucher-Giguère (2020) UVB spectral  \n4  \n<!-- image -->  \nFig. 11. Impact of the cosmic UV background photoionization on the Xray detectable phases. Top : Photon luminosity per unit of emission measure in the 0.5-2 keV band of the photoionized WHIM at di ff erent density (x-axis) and temperature (y-axis) by assuming Z = 0 . 2 Z ⊙ . Towards the high density (low ionization parameter) end, the radiative cooling reaches the CIE limit and is simply a function of temperature. When ∆ b ≲ 1, the radiative cooling keeps high at low temperature. Bottom : Same as the right panel of Fig. 10 but calculated using the photoionization model. By taking the cosmic UV background photoionization into account, the size of the 50% enclosed phases contour is marginally larger than that of the CIE condition.  \n<!-- image -->  \nenergy distribution for the calculation. We obtained a relation of log ξ [erg s -1 cm 2 ] = 2 . 26 -log ∆ b at z = 0 . 48. We used the pion model (Miller et al. 2015; Mehdipour et al. 2016) in SPEX-3.08.00 (Kaastra et al. 1996, 2024) for calculating the photoioized emission spectra and energy loss rates due to radiation at di ff erent ξ and temperature. We note that by setting the parameter tmod=1 , the pion model in SPEX allows for calculations to be carried out at fixed temperatures and ionization parameters without solving the energy balance, which meets the physical condition of WHIM, whose temperature is higher than the photoionization equilibrium temperature and is maintained by external heating mechanisms such as shock heating.  \nThe top panel of the Fig. 11 shows the 0.5-2 keV photon luminosity per unit of emission measure of the UVB photonion-  \nized WHIM at di ff erent temperature and density contrast. When ∆ b > 2 (log ξ [erg s -1 cm 2 ] ≲ 0), the radiative cooling reaches the CIE limit and simply becomes a function of temperature. When ∆ b ≲ 1 (log ξ [erg s -1 cm 2 ] ≳ 1), the radiative cooling deviates from the CIE condition and is a function of both temperature and density. The radiative cooling in the 0.5-2 keV band remains high for the WHIM, with a temperature below 10 6 K. We use the results of the photon luminosity of the photoionized WHIM to update the phase diagram weighted by the 0.5-2 keV emission and plot it in the bottom panel of Fig. 11. The 50% enclosed phases contour is at the same position as that of the CIE condition with a slightly larger size. It means that our conclusion drawn from Sect. 6 does not change by taking photonionization into account. The traceable WHIM by broadband X-ray emission is seen to be in the highest temperature and density phases among the whole filament gas population.', '7. Summary': "Overall, eRASS has provided the deepest All-Sky X-ray survey in the soft band, allowing for a more precise investigation of X-ray emission in cosmic filaments than in previous pioneering works (Tanimura et al. 2020b, 2022). This work aims to use the eRASS:4 X-ray data and SDSS optical filament catalog to stack X-ray emission from filaments and study the properties of the X-ray emitting WHIM. We adopted a novel method to: 1) estimate the fraction of non-WHIM emission in the stacked signal; 2) improve the spectral stacking method to minimize the systematics from the instrumental and sky background features at the observer frame; and 3) combine the SZ stacking results to simultaneously quantify the volume fraction and density of the X-ray emitting WHIM. We summarize our main results as follows:  \n- 1. We stacked 0.3-1.2 keV surface brightness profiles of 7817 SDSSoptical filaments with physical lengths of 20 -100 Mpc and redshifts of 0.2-0.6 in a 2275 deg 2 footprint. We detected an X-ray emission excess over the background at a 9 σ significance.\n- 2. We validated our profile stacking method by applying it to simulated count maps, including emission from filaments with an injected model and randomized X-ray background based on the observed angular power spectrum. After 200 realizations, the recovered averaged filament profiles agree well with the injected model with only a 13% di ff erence, which could be the result of a possible boosting from projectively overlapped filaments in the same orientations.\n- 3. We carefully estimate the fraction of stacked signal from unmasked sources, including X-ray halos, AGNs, and XRBs that are associated with galaxies. We first paint the observed galaxies in the analysis footprint with X-ray count rate from L X -M ∗ scaling relations to create an expected X-ray emission map from galaxies, where the incompleteness and bias of observed galaxy SMF is corrected. Then we apply the same surface brightness profile stacking method to obtain the fractions of stacked signal from each type of unmasked sources, which is 20%, 12%, and 5% for X-ray halos, AGNs, and XRBs, respectively. The photometric redshift uncertainties of the galaxies we used could underestimate another 3% of the contamination fraction. Based on these fractions, we conclude that the remaining ∼ 60% of the signal is from WHIMinthe filaments and the WHIM detection significance is 5 . 4 σ .  \n- 4. We stacked broadband spectra and response files of the filaments in the rest frame. We included the components of unmasked sources in spectral fitting. The best-fit temperature using a 1-T model is 0 . 58 ± 0 . 10 keV, which could be biased to the high end of the WHIM temperature distribution given its multi-temperature nature suggested by numerical simulations. The stacked spectra can also be fitted by a numerical simulation-motivated model with a log-normal temperature distribution, which will reduce the residual in the 0.50.6 keV energy range, suggesting the presence of the OVII abundant temperature phases.\n- 5. We modeled the baryon density contrast ∆ b of the X-ray emitting WHIM by fitting the stacked surface brightness profile, where we additionally introduced a parameter of volume filling factor f vol, which is necessary because only a small fraction of WHIM emits X-rays that can be observed by eROSITA. We used the SZ stacking results in the literature to break the degeneracy between the ∆ b and volume fraction f vol and obtained an averaged log ∆ b = 1 . 88 ± 0 . 18, suggesting that the observed emission is from the densest population of the WHIM.\n- 6. We analyzed the z = 0 . 48 snapshot of the TNG300-1 simulation to investigate the X-ray emitting phases in simulations. We found that most of the photons in the 0.5-2 keV energy range are emitted by phases of 6 . 2 < log( T / K) < 6 . 9 and 1 . 2 < log ∆ b < 2 . 2, which is in general agreement with our estimated values.  \nThis is the first time we detected > 5 σ X-ray emission from cosmic filaments. We applied a thorough investigation of the properties of the X-ray emitting WHIM. In the coming decade, our understanding of X-ray emissions from filaments and the WHIM within them will significantly improve from multiple aspects. The 4-metre Multi-Object Spectroscopic Telescope (4MOST) Cosmology Redshift Survey (Richard et al. 2019) will bring deep and wide field galaxy spectroscopic redshift data in the southern hemisphere, allowing us to take full use of the eRASS data of the German consortium. Improved filament finders such as the Monte Carlo Physarum Machine (Elek et al. 2022) and 1-DREAM (Canducci et al. 2022) will be applied to new survey data and provide a more accurately constructed filament catalog. Meanwhile, a better understanding of X-ray properties of galaxy groups, AGNs, and XRBs will help us to obtain a better estimate of the contamination fractions of unmasked sources. Moreover, high spectral resolution X-ray missions such as Hot Universe Baryon Surveyor (Cui et al. 2020) and Line Emission Mapper (Kraft et al. 2022) will be able to explore a wider parameter space of the WHIM properties.  \nAcknowledgements. We acknowledge Jelle Kaastra for the instruction on using the SPEX Pion model. We also acknowledge the anonymous referee for the insightful comments that improved the manuscript. X.Z., E.B., V.G., A.L., C.G., and S.Z. acknowledge financial support from the European Research Council (ERC) Consolidator Grant under the European Union's Horizon 2020 research and innovation program (grant agreement CoG DarkQuest No 101002585). M.C.H.Y. and M.F. acknowledge support from the Deutsche Forschungsgemeinschaft through the grant FR 1691 / 2-1. Y.Z. and G.P. acknowledge funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 865637). G.P. acknowledges support from Bando per il Finanziamento della Ricerca Fondamentale 2022 dell'Istituto Nazionale di Astrofisica (INAF): GO Large program and from the Framework per l'Attrazione e il Ra ff orzamento delle Eccellenze (FARE) per la ricerca in Italia (R20L5S39T9). M.B. acknowledges funding by the Deutsche Forschungsgemeinschaft under Germany's Excellence Strategy EXC 2121 'Quantum Universe' - 390833306. A.V. acknowledges funding by the Deutsche Forschungsgemeinschaft - 450861021. N. M. acknowledges funding by the European Union through a Marie Skłodowska-Curie Action Postdoctoral Fellowship (Grant Agreement: 101061448, project: MEMORY). Views and  \nopinions expressed are however those of the author only and do not necessarily reflect those of the European Union or of the Research Executive Agency. Neither the European Union nor the granting authority can be held responsible for them. This work is based on data from eROSITA, the soft X-ray instrument aboard SRG, a joint Russian-German science mission supported by the Russian Space Agency (Roskosmos), in the interests of the Russian Academy of Sciences represented by its Space Research Institute (IKI), and the Deutsches Zentrum für Luft und Raumfahrt (DLR). The SRG spacecraft was built by Lavochkin Association (NPOL) and its subcontractors and is operated by NPOL with support from the Max Planck Institute for Extraterrestrial Physics (MPE). The development and construction of the eROSITA X-ray instrument were led by MPE, with contributions from the Dr. Karl Remeis Observatory Bamberg & ECAP (FAU Erlangen-Nuernberg), the University of Hamburg Observatory, the Leibniz Institute for Astrophysics Potsdam (AIP), and the Institute for Astronomy and Astrophysics of the University of Tübingen, with the support of DLR and the Max Planck Society. The Argelander Institute for Astronomy of the University of Bonn and the Ludwig Maximilians Universität Munich also participated in the science preparation for eROSITA. The eROSITA data shown here were processed using the eSASS / NRTA software system developed by the German eROSITA consortium.", 'References': 'Aird, J., Coil, A. L., Georgakakis, A., et al. 2015, MNRAS, 451, 1892 Akamatsu, H., Fujita, Y., Akahori, T., et al. 2017, A&A, 606, A1  \n- Anderson, M. E., Gaspari, M., White, S. D. M., Wang, W., & Dai, X. 2015, MNRAS, 449, 3806\n- Arnaud, K. A. 1996, in Astronomical Society of the Pacific Conference Series, Vol. 101, Astronomical Data Analysis Software and Systems V, ed. G. H. Jacoby & J. Barnes, 17\n- Bahar, Y. E., Bulbul, E., Clerc, N., et al. 2022, A&A, 661, A7\n- Behroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. 2019, MNRAS, 488, 3143\n- Bernardi, M., Shankar, F., Hyde, J. B., et al. 2010, MNRAS, 404, 2087\n- Bi ffi , V., Planelles, S., Borgani, S., et al. 2017, MNRAS, 468, 531\n- Blanton, M. R. & Roweis, S. 2007, AJ, 133, 734\n- Blanton, M. R., Schlegel, D. J., Strauss, M. A., et al. 2005, AJ, 129, 2562\n- Böhringer, H., Briel, U. G., Schwarz, R. A., et al. 1994, Nature, 368, 828\n- Brunner, H., Liu, T., Lamer, G., et al. 2022, A&A, 661, A1\n- Buchner, J. 2016, Statistics and Computing, 26, 383\n- Buchner, J. 2019, PASP, 131, 108005\n- Buchner, J. 2021, The Journal of Open Source Software, 6, 3001\n- Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125\n- Bulbul, E., Liu, A., Kluge, M., et al. 2024, A&A, 685, A106\n- Bulbul, E., Markevitch, M., Foster, A., et al. 2014, ApJ, 789, 13\n- Bulbul, E., Randall, S. W., Bayliss, M., et al. 2016, ApJ, 818, 131\n- Canducci, M., Awad, P., Taghribi, A., et al. 2022, Astronomy and Computing, 41, 100658\n- Cen, R. & Ostriker, J. P. 1999, ApJ, 514, 1\n- Chen, Y.-M., Kau ff mann, G., Tremonti, C. A., et al. 2012, MNRAS, 421, 314  \nChurazov, E., Khabibullin, I. I., Dolag, K., Lyskova, N., & Sunyaev, R. A. 2023,  \nMNRAS, 523, 1209  \n- Comparat, J., Eckert, D., Finoguenov, A., et al. 2020, The Open Journal of Astrophysics, 3, 13\n- Comparat, J., Luo, W., Merloni, A., et al. 2023, A&A, 673, A122\n- Comparat, J., Merloni, A., Salvato, M., et al. 2019, MNRAS, 487, 2005\n- Costanzi, M., Rozo, E., Ryko ff , E. S., et al. 2019, MNRAS, 482, 490\n- Cui, W., Chen, L. B., Gao, B., et al. 2020, Journal of Low Temperature Physics, 199, 502  \nCui, W., Knebe, A., Libeskind, N. I., et al. 2019, MNRAS, 485, 2367  \n- Danforth, C. W., Keeney, B. A., Tilton, E. M., et al. 2016, ApJ, 817, 111\n- de Graa ff , A., Cai, Y.-C., Heymans, C., & Peacock, J. A. 2019, A&A, 624, A48 Dey, A., Schlegel, D. J., Lang, D., et al. 2019, AJ, 157, 168\n- Dorigo Jones, J., Johnson, S. D., Muzahid, S., et al. 2022, MNRAS, 509, 4330 Driver, S. 2021, Nature Astronomy, 5, 852\n- Eckert, D., Jauzac, M., Shan, H., et al. 2015, Nature, 528, 105\n- Elek, O., Burchett, J. N., Prochaska, J. X., & Forbes, A. G. 2022, Artificial Life, 28, 22\n- Fang, T., Buote, D. A., Humphrey, P. J., et al. 2010, ApJ, 714, 1715\n- Fang, T., Mathur, S., & Nicastro, F. 2022, Absorption Studies of the Most Di ff use Gas in the Large-Scale Structure, ed. C. Bambi & A. Santangelo (Singapore: Springer Nature Singapore), 1-43\n- Faucher-Giguère, C.-A. 2020, MNRAS, 493, 1614\n- Fukugita, M., Hogan, C. J., & Peebles, P. J. E. 1998, ApJ, 503, 518\n- Galárraga-Espinosa, D., Aghanim, N., Langer, M., & Tanimura, H. 2021, A&A, 649, A117\n- Gatuzz, E., García, J. A., Churazov, E., & Kallman, T. R. 2023, MNRAS, 521, 3098', 'Appendix A: Surface brightness profiles in different bands and length populations': 'In Fig. A.1, we plot the stacked net profile in the 0.3-0.7 keV and 0.7-1.2 keV bands. The detections in the two bands are both significant, reaching 7 . 0 σ and 6 . 5 σ levels, respectively. In contrast to previous works of Tanimura et al. (2020b, 2022), we observe significant X-ray emission in the < 0.6 keV energy band.  \nWe also split the filament sample into 20-40 Mpc and 40100 Mpc length populations and we detect a marginal 2 . 5 σ discrepancy between the two populations (see Fig. A.2). The short filaments have a broader stacked emission profile than long filaments, consistent with the predictions in numerical simulations (Galárraga-Espinosa et al. 2021).  \nFig. A.1. Same as Fig. 3, but in the 0.3-0.7 keV ( left ) and 0.7-1.2 keV ( right ) bands, respectively.  \n<!-- image -->  \nFig. A.2. Local background-subtracted surface brightness profiles of the 20-40 Mpc length (red) and the 40-100 Mpc length (blue) filament samples. The dashed horizontal line denotes the zero level.  \n<!-- image -->', 'Appendix B: Impact of the source masking radius': "The choice of source masking radius could lead to a change in the stacked profile. Tanimura et al. (2022) has investigated the impact of di ff erent cluster mask radii and point source mask radii on the results. They discovered that the increase in the cluster mask radii from r 500 to 3 r 500 and the increase in the point source mask radii from 30 '' to 45 '' do not a ff ect the results (within the uncertainties). In this section, we present our independent check of this systematic uncertainty in our data.  \nFor cluster masks, instead of the 1 . 5 r 500 for the eRASS1 Xray selected sample and 1 . 5 r λ for the redMaPPer sample (as described in Sect. 2.1), we increased the masking radii by a factor of two for a test. Following Sects. 2.1 and 2.2, we create new 0.3-1.2 keV count and exposure Healpix maps with clusters masked using 3 r 500 and 3 r λ masks. Other source masks are the same as those listed in Table 1. In the left panel of Fig. B.1, we plot the stacked net profiles with the doubled cluster masking radii. The profile with doubled mask radii agrees with the original one within uncertainties. The integrated net signal within 10 Mpc transverse distance is 93% of the original integrated signal. We note that the uncertainty in each bin is increased by 13%, which is due to the significant reduction of the available footprint (see the comparison between Fig. 1 and Fig. B.2).  \nSimilarly, for the eRASS:4 X-ray sources, we tested on the masking radii that are calculated as S ( r cut) = 0 . 05 S bkg and plot the comparison in the right panel of Fig. B.1. The profiles using two di ff erent mask radii are identical within the uncertainties.  \n<!-- image -->  \nFig. B.1. Impact of the source masking radius on the stacked profile. Left : Comparison between stacked net profile using 1 . 5 r 500 (red) and 3 r 500 (blue) cluster mask radii. The stacked profile with large cluster masks agrees with the original profile within the uncertainties. Right : Comparison between stacked net profiles using X-ray source masks with r | S = 0 . 1 S bkg (red) and r | S = 0 . 05 S bkg (blue). A more conservative choice of the source mask radius does not change the result within uncertainties.  \n<!-- image -->  \nFig. B.2. Same as Fig. 1, but with 3 r 500 (3 r λ ) mask radii for X-ray selected (optically identified) clusters. The analysis footprint without source masks is significantly reduced.  \n<!-- image -->", 'Appendix C: Mock maps for profile stacking validation': 'First, we phenomenically fit the stacked net surface brightness profile using a King-profile model (see the section on modeling gas properties in the filaments). With the best-fit parameters, we created the mock filament emission rate map based on the projected geometric configurations of the 7817 filaments (top left panel of Fig. C.1). We note that the best-fit parameters in this section are empirical and we did not use them for physical interpretation. The goal is to test the methodology and algorithms used to extract stacked surface brightness profiles. The best-fit profile for injection is plotted as the dotted line in Fig. 5.  \nWe simulated mock background rate maps based on the angular power spectrum of the observed rate map (top right panel of Fig. C.1). By comparing the mock filament emission map and the observed sky map, we know that the filament emission contributes < 1% of the total signal. Therefore, we assume that the angular power spectrum of the observed rate map adequately represents that of the background component. We first smooth the rate map using a 1 · full width at half maximum (FWHM) Gaussian kernel to reduce the Poisson noise. Then, we construct a relative amplitude map based on the mean and standard deviation of the smoothed rate map. We computed the angular power spectrum of the relative amplitude map Cl , rel up to a maximum multipole moment l = 3000 using the function anafast in the HEALPYpackage. We simulated a mock relative amplitude map using this angular power spectrum and convert it to a mock rate map using the mean and standard deviation of the smoothed rate map. We summed up the mock filament rate map and the mock background rate map to get a mock sky rate map. We obtained a mock sky count map by multiplying the exposure map from observations and applying Poisson randomization to the expected count value in each pixel.  \nFig. C.1. Example of simulated maps in the 0.3-1.2 keV band for validating the imaging stacking method. Top left : Mock count rate map for filament emission. Top right : Mock foreground / background emission mock map. Bottom left : Mock count map of the total emission. Bottom right : Observed count map as a comparison.  \n<!-- image -->', 'Appendix D: Spectral models for undetected galactic sources': 'In Fig. D.1, we plot the emission measure distribution of unmasked halos at di ff erent redshifts (left panel), the composite spectra of unmasked halos at di ff erent redshifts (middle panel), and the averaged spectrum (right panel). In Fig. D.2, we plot the spectral models in di ff erent redshift bins, stellar mass bins used in the mock catalog, and the median spectrum.  \nFig. D.1. Left : Multi-temperature di ff erential emission measure distribution of undetected X-ray halos at di ff erent redshifts after applying the selection. Middle : Corresponding composite spectra of each di ff erential emission measure distribution at di ff erent redshifts. Right : Synthesis restframe spectral model of undetected X-ray halos by taking the filament sample length and redshift distributions into account.  \n<!-- image -->  \nFig. D.2. Individual undetected AGN spectral models in di ff erent redshift - stellar mass bins of the mock AGN catalog with the normalizations rescaled. The spectra show strong self-similarity in the soft band. The dotted is the median of all the spectra.  \n<!-- image -->'}