© The Authors 2024

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⁻⁴ cm⁻³, 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 XIII 2016; Planck Collaboration VI 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 techniques, such as X-ray observations and Sunyaev–Zeldovich (SZ) effects, remains one of the primary investigation goals in modern astrophysics (Driver 2021).

Currently, absorption studies by using OVI lines in the far ultraviolet bands can successfully trace the 10^5.5 K low-temperature phases WHIM (Danforth et al. 2016). For X-ray absorption lines that trace ∼10⁶ K phases WHIM, while there are many cases of plausible detections (Nicastro et al. 2005, 2018; Fang et al. 2010), 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 et 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.

Wide-area optical spectroscopic surveys have mapped the large-scale galaxy distributions up to redshifts of 0.7 (Reid et al. 2016), thus offering 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. 2020a, 2022), demonstrating the potential of X-ray stacking analysis to trace the emission in cosmic filaments. Similarly, mild positive detections have also been reported by stacking the SZ effect signal that traces the pressure of the intergalactic medium (de Graaff et al. 2019; Tanimura et al. 2019, 2020b). 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⁷ 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.

In this work, we assume a fiducial Lambda cold dark matter cosmology with parameters: H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_m = 0.3, Ω_Λ = 0.7. We adopt the ratio between baryon density parameter to matter density parameter from the Planck18 cosmology results, namely: Ω_b/Ω_m = 0.158 (Planck Collaboration VI 2020). At z = 0, the mean baryon density is ${\overline{\rho }}_{\mathrm{b}}(0)=4.51\times {10}^{-31}$ g cm⁻³. The baryon density contrast is defined as ${\mathrm{\Delta }}_{\mathrm{b}}\equiv {\rho }_{\mathrm{b}}/{\overline{\rho }}_{\mathrm{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 improvements 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.

We 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+CMASS footprint (Reid et al. 2016) and the western Galactic hemisphere of the eROSITA All-Sky Survey determines a ∼3200 deg² 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². 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).

Fig. 1. Smoothed 0.3–1.2 keV eROSITA count rate map of the analysis footprint. The missing pixels show source-masked regions.

Fig. 2. Distribution of the selected filaments in the redshift and physical length space.

2.1. Source masking

We masked out all eRASS:4 X-ray sources in the 0.2–2.3 keV eSASS detection catalog¹. 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, we did not apply selection cuts to the contaminating X-ray emitting sources.

We used a background level-dependent radius to exclude detected sources as follows. For each source, the surface brightness profile can be approximately described as

$$\begin{array}{c}\hfill S(r)=S_0{(1+\frac{r^2}{r_{\mathrm{PSF}}^2+r_{\mathrm{EXT}}^2})}^{-1.5},\end{array}$$(1)

where 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.

For 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₅₀₀. Additionally, we adopted an optical cluster catalog compiled by the redMaPPer algorithm (Rykoff 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⁻¹ 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.

To 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₅₀₀. If a source displays 1.5 × r₅₀₀ < 400 kpc, we used 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 large-scale structure (see Malavasi et al. 2020, for details). Details on all masked sources are listed in Table 1.

We 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 suffer from significantly high contamination due to projection effects (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² (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 = ∑_jC_ij, where C_ij 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 = (∑_jt_ij × Ω_ij)/Ω_HP, i, where ∑_jt_ij and Ω_ij 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 𝒫_i, the j th filament ℱ_j and the k th distance bin ℬ_k, we assigned a Boolean coefficient b_ijk, whose value is 1 (or 0) if 𝒫_i is in (or not in) the k th distance bin of ℱ_j. Therefore, for each filament, ℱ_j, the count number in the kth transverse distance bin is $C_{\mathit{jk}}={\displaystyle \sum _i^{n_{\mathrm{HP}}}}{b}_{\mathit{ijk}}\times C_i$, and the exposure time in the kth transverse distance bin is $t_{\mathit{jk}}={\displaystyle \sum _i^{n_{\mathrm{HP}}}}{b}_{\mathit{ijk}}\times t_i$. We used exposure-weighted averages to obtain the stacked surface brightness profile. The stacked profile,

$$\begin{array}{c}\hfill S_{\mathrm{X}}(r_k)=\frac{1}{{\sum }_j^{n_{\mathrm{fil}}}{t}_{\mathit{jk}}}{\displaystyle \sum _j^{n_{\mathrm{fil}}}}{t}_{\mathit{jk}}\times \frac{C_{\mathit{jk}}}{t_{\mathit{jk}}}=\frac{{\sum }_j^{n_{\mathrm{fil}}}{C}_{\mathit{jk}}}{{\sum }_j^{n_{\mathrm{fil}}}{t}_{\mathit{jk}}},\end{array}$$(2)

is the quotient between the sum of the count profile and the sum of the exposure profile, where r_k 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

$$\begin{array}{c}\hfill S_{\mathrm{X},\mathrm{net}}(r_k)=S_{\mathrm{X}}(r_k)-S_{\mathrm{X},\mathrm{bkg}}.\end{array}$$(3)

Due to the projection effects, the distance bins of different filaments at different redshifts could be simultaneously located at a single position. Therefore, the surface brightness in different 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 lth resample of the pixel set with the replacement is

$$\begin{array}{c}\hfill {\mathcal{P}}_l^{\prime }=({\mathcal{P}}_{l1}^{\prime }{\mathcal{P}}_{l2}^{\prime }\dots {\mathcal{P}}_{l{n}_{\mathrm{HP}}}^{\prime }),\mathrm{with}{\mathcal{P}}_{\mathit{li}}^{\prime }\in \mathcal{P}.\end{array}$$(4)

Following Eqs. (2) and (3), for each resampled pixel set [BOLD]𝒫_l′, we obtained a stacked surface brightness profile ${\mathit{S}}_{\mathrm{X},\mathrm{net},l}\prime (r_k)$. After drawing n_bs random resamples, the bootstrapping estimate of the net surface brightness profile is

$$\begin{array}{c}\hfill {\widehat{S}}_{\mathrm{X},\mathrm{net}}(r_k)=\frac{1}{n_{\mathrm{bs}}-1}{\displaystyle \sum _l^{n_{\mathrm{bs}}}}{S}_{\mathrm{X},\mathrm{net},l}^{\prime }(r_k),\end{array}$$(5)

and the bootstrapping estimate of the covariance between the k₁ th and k₂ th bin is

$$\begin{array}{cc}\hfill {\mathcal{C}}_{k_1{k}_2}=& \frac{1}{n_{\mathrm{bs}}-1}\hfill \\ \hfill & \times {\displaystyle \sum _l^{n_{\mathrm{bs}}}}[S_{\mathrm{X},\mathrm{net},l}^{\prime }(k_1)-{\widehat{S}}_{\mathrm{X},\mathrm{net}}(k_1)][S_{\mathrm{X},\mathrm{net},l}^{\prime }(k_2)-{\widehat{S}}_{\mathrm{X},\mathrm{net}}(k_2)].\hfill \end{array}$$(6)

The uncertainty of ${\widehat{S}}_{\mathrm{X},\mathrm{net},k}$ can be taken from the diagonal element of the covariance matrix ${\widehat{\sigma }}_{S_{\mathrm{X}},k}^2={\mathcal{C}}_{\mathit{kk}}$.

3.1. Stacking the surface brightness profile

We 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. (2020a); Tanimura et al. (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 effects 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 X-ray 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.

Fig. 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 random-positioned 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.

To 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 different 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.

Furthermore, we investigated the impact of different 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.05S​_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 different redshifts could be overlapped with each other due to projection effects. 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, respectively. There are two different cases of overlapping. In the first case, the two overlapped filaments have different 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.

Fig. 4. Sketches of two types of overlapped filaments due to the projection effects. Deep and shallow regions are the source and local background regions, respectively. In the first case (left), the two filaments have different 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.

To 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. We simulated 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 effect 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 level, our stacking method can successfully recover the emission profile.

Fig. 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.

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 binaries (XRBs). The steps required to model the contribution from these sources in the stacked total signal are given below.

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.

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¹¹ 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 different stellar masses and redshifts are used to correct the predicted count rate of each galaxy.

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 − 2 keV} − M_* relation for unmasked X-ray halos in different redshift bins by applying the same source selection criteria.

The galaxy sample used in Anderson et al. (2015) was selected from an SDSS-DR7 based New York University Value Added 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 difference on the estimated unmasked X-ray 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 cross-match 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 − 2 keV} − M_* relation. The modified X-ray halo L_X − M_* relation is plotted in the top left panel of Fig. 6).

Fig. 6. Rest-frame galaxy averaged 0.5–2 keV luminosity – stellar mass scaling relations for unmasked X-ray halos (top-left), AGNs (top-right), and XRBs (bottom-left) at different 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.

To convert the rest-frame 0.5–2 keV luminosity to the 0.3–1.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).

With 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 differential 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:

$$\begin{array}{c}\hfill F(E)={\displaystyle \sum _i}\frac{l_i\times F(E,z_i)}{D_{\mathrm{L}}{(z_i)}^2},\end{array}$$(7)

where l_i is the physical length of the i-th filament, z_i is the redshift of the i-th filament, F(E, z_i) is the multi-temperature differential emission measure spectrum at the redshift of z_i, 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 high-resolution 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).

We applied a selection criteria that removes AGNs with F_{0.5 − 2 keV} > 10⁻¹⁴ erg s⁻¹ cm⁻², which would be detected in eRASS:4, to select sources in the mock catalog. In each redshift-stellar 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 obtain 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 different redshifts are shown in the top right panel of Fig. 6.

We 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-ray-excess, 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 different 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 different 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¹¹ M_⊙ (the result of the redshift evolution of star-forming rate and star-forming galaxy fraction).

We adopted a power-law distribution with intrinsic absorption spectral model for flux conversion, with parameters of Γ = 1.7 and n_H = 2 × 10²¹ cm⁻², 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 different stellar masses as the observed relation and the slight luminosity difference 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 contribution 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%.

Fig. 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.

Because 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.

We 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

$$\begin{array}{c}\hfill f=\frac{D_{\mathrm{L}}^2(z_{\mathrm{true}})}{n_{\mathrm{MC}}}{\displaystyle \sum _i}\frac{w_i}{D_{\mathrm{L}}^2(z_i)},\end{array}$$(8)

where n_MC is the number of realizations, and z_i is the i th randomized redshift value, while the value of w_i is 1 if 0.15 < z_i < 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%.

Fig. 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.

The ∼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 modeling, which is based on our current knowledge of different X-ray emitters in galaxy systems at different scales. A better understanding of all the different types of X-ray sources in the Universe will offer 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 different 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 Galactic n_H value (HI4PI Collaboration 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).

We used tools addspec, addarf, and addrmf in FTOOLS³ (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 different definitions from those in other X-ray satellite data products⁴. 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:

$$\begin{array}{c}\hfill N_{\mathrm{stacked},\mathrm{src}}(E)={\displaystyle \sum _i}{N}_{\mathrm{src},i}(E).\end{array}$$(9)

We then scale each local background spectrum to match the extraction solid angle of the corresponding source spectrum to get the stacked background spectrum,

$$\begin{array}{c}\hfill N_{\mathrm{stacked},\mathrm{bkg}}(E)={\displaystyle \sum _i}\frac{{\mathrm{\Omega }}_{\mathrm{REG},\mathrm{src},i}}{{\mathrm{\Omega }}_{\mathrm{REG},\mathrm{bkg},i}}\times N_{\mathrm{bkg},i}(E).\end{array}$$(10)

For the final stacked source spectrum, we assign the exposure time as a solid angle weighted averaged exposure ${\overline{t}}_{\mathrm{ave}}=({\displaystyle \sum _i}{\mathrm{\Omega }}_{\mathrm{REG},i}\times t_{\mathrm{EXP},i})/{\displaystyle \sum _i}{\mathrm{\Omega }}_{\mathrm{REG},i}$. To match this exposure, the stacked ARF is defined as

$$\begin{array}{c}\hfill A_{\mathrm{stacked}}(E)=\frac{1}{{\overline{t}}_{\mathrm{ave}}}{\displaystyle \sum _i}{t}_{\mathrm{EXP},i}{\mathrm{\Omega }}_{\mathrm{REG},i}{A}_i(E),\end{array}$$(11)

which takes each extraction solid angle into account. The modeled source flux using the exposure and stacked ARF is scaled per deg². Similarly, the stacked RMF is weighted by Ω_REG × t_ave, which is

$$\begin{array}{c}\hfill {\mathcal{R}}_{\mathrm{stacked}}=\frac{1}{{\sum }_i{\mathrm{\Omega }}_{\mathrm{BS},i}{t}_{\mathrm{EXP},i}}\times {\displaystyle \sum _i}{\mathrm{\Omega }}_{\mathrm{BS},i}{t}_{\mathrm{EXP},i}{\mathcal{R}}_i.\end{array}$$(12)

The 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.

One clear advantage of the spectral stacking method in the rest frame is reducing the effect 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 effect 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.

Fig. 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 signal-to-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.

We used XSPEC version 12.13.0 and the BXA⁵ 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.5–6.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.

We 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; Biffi et al. 2017) and > 0.1 Z_⊙ in z < 1 filaments from numerical simulations (Martizzi et 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.

Recently, 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 best-fit temperature using a 1-T model does not agree with the majority of the gas phases around 10⁶ K predicted 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.8–0.9 keV energy range, dominated by several Fe XVII lines emitted by 10^6.8 K gas.

We 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. Detected phases of WHIM

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³ to 10⁷ K in temperature and from 10⁻⁶ to 10⁻¹ cm⁻³ 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

$$\begin{array}{c}\hfill {\overline{ϵ}}_{\det }(r)=1.20\times n_{\mathrm{H}}{(r)}^2{f}_{\mathrm{vol}}(r),\end{array}$$(13)

where 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,

$$\begin{array}{cc}& n_{\mathrm{H}}(r)=n_{\mathrm{H},0}{[1+{(r/r_{\mathrm{c},\mathrm{n}})}^2]}^{{\alpha }_{\mathrm{n}}};\hfill \end{array}$$(14)

$$\begin{array}{cc}& f_{\mathrm{vol}}(r)=f_{\mathrm{vol},0}{[1+{(r/r_{\mathrm{c},\mathrm{fvol}})}^2]}^{{\alpha }_{\mathrm{fvol}}},\hfill \end{array}$$(15)

where 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.

In 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,

$$\begin{array}{c}\hfill n_{\mathrm{H},0}(z)=\frac{{\mathrm{\Delta }}_{\mathrm{b},0}{\rho }_{\mathrm{crit}}{\mathrm{\Omega }}_{\mathrm{b}}}{1.43\times m_{\mathrm{H}}}\times {(1+z)}^3,\end{array}$$(16)

where ρ_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

$$\begin{array}{c}\hfill S_{\mathrm{X},\mathrm{model},\mathrm{conv}}(r_{\mathrm{proj}},z)=G({\sigma }_{\mathrm{sm}})\ast S_{\mathrm{X},\mathrm{model}}(r_{\mathrm{proj}},z),\end{array}$$(17)

where G(σ_sm) is a Gaussian convolving kernel that accounts for the smoothing effect of the observed filament width. The surface brightness is

$$\begin{array}{cc}\hfill S_{\mathrm{X},\mathrm{model}}(r_{\mathrm{proj}},z)=& A_{\mathrm{CF}}^{-1}{K}_{\mathrm{ph}}{(1+z)}^{-3}\hfill \\ \hfill & \times {\int }_{-l_{\mathrm{out}}}^{l_{\mathrm{out}}}{\overline{ϵ}}_{\det }(\sqrt{l^2+r_{\mathrm{proj}}^2},z)\mathrm{\Lambda }\mathrm{d}l,\hfill \end{array}$$(18)

where 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)⁻³ 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.

When 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

$$\begin{array}{c}\hfill S_{\mathrm{X},\mathrm{model},\mathrm{stack}}=\frac{1}{N}{\displaystyle \sum _i^N}{\mathrm{\Omega }}_i\times S_{\mathrm{X},\mathrm{model},\mathrm{conv},i}.\end{array}$$(19)

The fitting statistics is defined as

$$\begin{array}{c}\hfill {\chi }^2={{\mathit{S}}_{\mathrm{X},\mathrm{diff}}}^{⊺}{\mathcal{C}}^{-1}{\mathit{S}}_{\mathrm{X},\mathrm{diff}},\end{array}$$(20)

where S_X, diff = S_{X, net, corr} − S_{X, unmasked, net} − S_{X, model, stack} × 1.13, is the difference 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 [BOLD]𝒞 is the covariance matrix estimated from the bootstrapping method.

We 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 statistics χ², the likelihood function includes an additional constraint from SZ stacking results (Tanimura et al. 2020b), which is necessary to break up the degeneracy between Δ_b, 0 and f_vol, 0. The predicted Compton-y strength is

$$\begin{array}{c}\hfill y_{\mathrm{pd}}(r_{\mathrm{proj}})=\left[\frac{{\sigma }_{\mathrm{T}}}{m_{\mathrm{e}}{c}^2}{\displaystyle {\int }_{-l_{\mathrm{out}}}^{l_{\mathrm{out}}}}kT{n}_{\mathrm{e}}\left(\sqrt{l^2+r_{\mathrm{proj}}^2}\right)f_{\mathrm{vol}}\left(\sqrt{l^2+r_{\mathrm{proj}}^2}\right)\mathrm{d}l\right],\end{array}$$(21)

where σ_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. (2020b) show similar broadening effects as our stacked X-ray profile. Therefore, we also convolved the predicted Compton-y profile with the same Gaussian kernel to create a smoothed Compton-y profile as

$$\begin{array}{c}\hfill y_{\mathrm{pd},\mathrm{conv}}=G({\sigma }_{\mathrm{sm}})\ast y_{\mathrm{pd}}.\end{array}$$(22)

Given 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. (2020b). We used a likelihood term P(log y_pd, conv(0)) = 𝒩(μ_{log y}, σ_{log y}) as the addition constraint. We do not know how much of the total stacked Compton-y is contributed by the X-ray-emitting gas. Nevertheless, we adopt μ_{log y} = −8.3 and σ_{log y} = 0.5 in our analysis. The combined likelihood function is

$$\begin{array}{c}\hfill ln\mathcal{L}=\{\begin{array}{cc}-\frac{1}{2}{\chi }^2-\frac{1}{2}\frac{{[y_{\mathrm{pd},\mathrm{conv}}(0)-{\mu }_{\log y}]}^2}{{\sigma }_{\log y}^2}\hfill & \mathrm{when}log{y}_{\mathrm{pd},\mathrm{conv}}(0)\leqq -7.5,\hfill \\ -\infty \hfill & \mathrm{elsewhere}.\hfill \end{array}\end{array}$$(23)

We used a Markov chain Monte Carlo (MCMC) method and the python package Zeus⁶ (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{\overline{\mathrm{\Delta }}}_{\mathrm{b}}=1.88\pm 0.18$, which corresponds to a hydrogen density $4.4_{-1.5}^{+2.2}\times {10}^{-5}$ cm⁻³ 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.

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).

To 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⁷. 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.

We 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⁷ K is not as prominent as those in the z = 0 snapshot (Galárraga-Espinosa 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).

Fig. 10. llustrisTNG 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.

To 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 represents the traceable population of SZ observations (e.g. de Graaff et al. 2019; Tanimura et al. 2019, 2020b), 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 differ. SZ observations trace 10^5.5 − 10^6.8 K temperature and a Δ_b = 10^0.3 − 10² 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 multi-temperature nature.

6.3. Impact of photoionization

It is known that the WHIM is also partially ionized by the cosmic ultraviolet background (UVB), which alters the ionization balance and exhibits different 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.

The ionization parameter is defined as

$$\begin{array}{c}\hfill \xi \equiv \frac{4\pi F_{1-1000\mathrm{Ry}}}{n_{\mathrm{H}}},\end{array}$$(24)

where F_{1 − 1000 Ry} 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 energy distribution for the calculation. We obtained a relation of log ξ[erg s⁻¹ cm²] = 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 different ξ 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.

The top panel of the Fig. 11 shows the 0.5–2 keV photon luminosity per unit of emission measure of the UVB photonionized WHIM at different temperature and density contrast. When Δ_b > 2 (log ξ[erg s⁻¹ cm²]≲0), the radiative cooling reaches the CIE limit and simply becomes a function of temperature. When Δ_b ≲ 1 (log ξ[erg s⁻¹ cm²]≳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⁶ 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.

Fig. 11. Impact of the cosmic UV background photoionization on the X-ray detectable phases. Top: Photon luminosity per unit of emission measure in the 0.5–2 keV band of the photoionized WHIM at different density (x-axis) and temperature (y-axis) by assuming Z = 0.2Z⊙. 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.

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. 2020a, 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:

1. We stacked 0.3–1.2 keV surface brightness profiles of 7817 SDSS optical filaments with physical lengths of 20 − 100 Mpc and redshifts of 0.2–0.6 in a 2275 deg² footprint. We detected an X-ray emission excess over the background at a 9σ significance.

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% difference, which could be the result of a possible boosting from projectively overlapped filaments in the same orientations.

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 WHIM in the filaments and the WHIM detection significance is 5.4σ.

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.5–0.6 keV energy range, suggesting the presence of the OVII abundant temperature phases.

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.

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.

This 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.

Acknowledgments

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 Rafforzamento 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 opinions 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

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. (2020a); Tanimura et al. (2022), we observe significant X-ray emission in the < 0.6 keV energy band.

We also split the filament sample into 20–40 Mpc and 40–100 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).

Fig. A.1. Same as Fig. 3, but in the 0.3–0.7 keV (left) and 0.7–1.2 keV (right) bands, respectively.

Fig. 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.

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 different cluster mask radii and point source mask radii on the results. They discovered that the increase in the cluster mask radii from r₅₀₀ to 3r₅₀₀ and the increase in the point source mask radii from 30″ to 45″ do not affect the results (within the uncertainties). In this section, we present our independent check of this systematic uncertainty in our data.

For cluster masks, instead of the 1.5r₅₀₀ for the eRASS1 X-ray selected sample and 1.5r_λ 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 3r₅₀₀ and 3r_λ 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).

Similarly, for the eRASS:4 X-ray sources, we tested on the masking radii that are calculated as S(r_cut) = 0.05S​_bkg and plot the comparison in the right panel of Fig. B.1. The profiles using two different mask radii are identical within the uncertainties.

Fig. B.1. Impact of the source masking radius on the stacked profile. Left: Comparison between stacked net profile using 1.5r500 (red) and 3r500 (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.1S​bkg (red) and r|S = 0.05S​bkg (blue). A more conservative choice of the source mask radius does not change the result within uncertainties.

Fig. B.2. Same as Fig. 1, but with 3r500 (3rλ) mask radii for X-ray selected (optically identified) clusters. The analysis footprint without source masks is significantly reduced.

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.

We 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 C_l, rel up to a maximum multipole moment l = 3000 using the function anafast in the HEALPY package. 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.

Fig. 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.

Appendix D: Spectral models for undetected galactic sources

In Fig. D.1, we plot the emission measure distribution of unmasked halos at different redshifts (left panel), the composite spectra of unmasked halos at different redshifts (middle panel), and the averaged spectrum (right panel). In Fig. D.2, we plot the spectral models in different redshift bins, stellar mass bins used in the mock catalog, and the median spectrum.

Fig. D.1. Left: Multi-temperature differential emission measure distribution of undetected X-ray halos at different redshifts after applying the selection. Middle: Corresponding composite spectra of each differential emission measure distribution at different redshifts. Right: Synthesis rest-frame spectral model of undetected X-ray halos by taking the filament sample length and redshift distributions into account.

Fig. D.2. Individual undetected AGN spectral models in different 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.

All Tables

All Figures

	Fig. 1. Smoothed 0.3–1.2 keV eROSITA count rate map of the analysis footprint. The missing pixels show source-masked regions.
In the text

	Fig. 2. Distribution of the selected filaments in the redshift and physical length space.
In the text

Fig. 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 random-positioned 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.

In the text

Fig. 4. Sketches of two types of overlapped filaments due to the projection effects. Deep and shallow regions are the source and local background regions, respectively. In the first case (left), the two filaments have different 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.

In the text

Fig. 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.

In the text

	Fig. 6. Rest-frame galaxy averaged 0.5–2 keV luminosity – stellar mass scaling relations for unmasked X-ray halos (top-left), AGNs (top-right), and XRBs (bottom-left) at different 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.
In the text

	Fig. 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.
In the text

Fig. 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.

In the text

Fig. 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 signal-to-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.

In the text

Fig. 10. llustrisTNG 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.

In the text

Fig. 11. Impact of the cosmic UV background photoionization on the X-ray detectable phases. Top: Photon luminosity per unit of emission measure in the 0.5–2 keV band of the photoionized WHIM at different density (x-axis) and temperature (y-axis) by assuming Z = 0.2Z⊙. 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.

In the text

	Fig. A.1. Same as Fig. 3, but in the 0.3–0.7 keV (left) and 0.7–1.2 keV (right) bands, respectively.
In the text

	Fig. 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.
In the text

Fig. B.1. Impact of the source masking radius on the stacked profile. Left: Comparison between stacked net profile using 1.5r500 (red) and 3r500 (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.1S​bkg (red) and r|S = 0.05S​bkg (blue). A more conservative choice of the source mask radius does not change the result within uncertainties.

In the text

	Fig. B.2. Same as Fig. 1, but with 3r500 (3rλ) mask radii for X-ray selected (optically identified) clusters. The analysis footprint without source masks is significantly reduced.
In the text

	Fig. 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.
In the text

Fig. D.1. Left: Multi-temperature differential emission measure distribution of undetected X-ray halos at different redshifts after applying the selection. Middle: Corresponding composite spectra of each differential emission measure distribution at different redshifts. Right: Synthesis rest-frame spectral model of undetected X-ray halos by taking the filament sample length and redshift distributions into account.

In the text

	Fig. D.2. Individual undetected AGN spectral models in different 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.
In the text