11institutetext: Max-Planck-Institut fur Extraterrestrische Physik,
Giessenbachstrasse, 85748, Garching, Germany
11email<EMAIL_ADDRESS>22institutetext: Hamburger Sternwarte,
Gojenbergsweg 112, 21029, Hamburg, Germany 33institutetext: University
Observatory Munich, Faculty of Physics, Ludwig-Maximilians-Universität,
Scheinerstr. 1, 81679, Munich, Germany 44institutetext: Universität Innsbruck,
Institut für Astro- und Teilchenphysik, Technikerstrasse 25/8, 6020,
Innsbruck, Austria 55institutetext: Arnold Sommerfeld Center for Theoretical
Physics, LMU Munich, Theresienstr. 37, 80333 Munich, Germany 66institutetext:
INAF, Osservatorio Astronomico di Brera, Via E. Bianchi 46, 23807, Merate,
(LC), Italy 77institutetext: Argelander Institute for Astronomy, University of
Bonn, Auf dem Hügel 71, 53121, Bonn, Germany

# The SRG/eROSITA all-sky survey.

X-ray emission from the warm-hot phase gas in long cosmic filaments

X. Zhang The SRG/eROSITA all-sky survey.The SRG/eROSITA all-sky survey. E.
Bulbul The SRG/eROSITA all-sky survey.The SRG/eROSITA all-sky survey. N.
Malavasi The SRG/eROSITA all-sky survey.The SRG/eROSITA all-sky survey. V.
Ghirardini The SRG/eROSITA all-sky survey.The SRG/eROSITA all-sky survey. J.
Comparat The SRG/eROSITA all-sky survey.The SRG/eROSITA all-sky survey. M.
Kluge The SRG/eROSITA all-sky survey.The SRG/eROSITA all-sky survey. A. Liu
The SRG/eROSITA all-sky survey.The SRG/eROSITA all-sky survey. A. Merloni The
SRG/eROSITA all-sky survey.The SRG/eROSITA all-sky survey. Y. Zhang The
SRG/eROSITA all-sky survey.The SRG/eROSITA all-sky survey. Y. E. Bahar The
SRG/eROSITA all-sky survey.The SRG/eROSITA all-sky survey. E. Artis The
SRG/eROSITA all-sky survey.The SRG/eROSITA all-sky survey. J. Sanders The
SRG/eROSITA all-sky survey.The SRG/eROSITA all-sky survey. C. Garrel The
SRG/eROSITA all-sky survey.The SRG/eROSITA all-sky survey. F. Balzer The
SRG/eROSITA all-sky survey.The SRG/eROSITA all-sky survey. M. Brüggen The
SRG/eROSITA all-sky survey.The SRG/eROSITA all-sky survey. M. Freyberg The
SRG/eROSITA all-sky survey.The SRG/eROSITA all-sky survey. E. Gatuzz The
SRG/eROSITA all-sky survey.The SRG/eROSITA all-sky survey. S. Grandis The
SRG/eROSITA all-sky survey.The SRG/eROSITA all-sky survey. S. Krippendorf The
SRG/eROSITA all-sky survey.The SRG/eROSITA all-sky survey.The SRG/eROSITA all-
sky survey.The SRG/eROSITA all-sky survey. K. Nandra The SRG/eROSITA all-sky
survey.The SRG/eROSITA all-sky survey. G. Ponti The SRG/eROSITA all-sky
survey.The SRG/eROSITA all-sky survey.The SRG/eROSITA all-sky survey.The
SRG/eROSITA all-sky survey. M. Ramos-Ceja The SRG/eROSITA all-sky survey.The
SRG/eROSITA all-sky survey. P. Predehl The SRG/eROSITA all-sky survey.The
SRG/eROSITA all-sky survey. T. H. Reiprich The SRG/eROSITA all-sky survey.The
SRG/eROSITA all-sky survey. A. Veronica The SRG/eROSITA all-sky survey.The
SRG/eROSITA all-sky survey. M. C. H. Yeung The SRG/eROSITA all-sky survey.The
SRG/eROSITA all-sky survey. S. Zelmer The SRG/eROSITA all-sky survey.The
SRG/eROSITA all-sky survey.

The properties of the warm-hot intergalactic medium (WHIM) in the cosmic
filaments are among the least quantified in modern astrophysics. eROSITA All
Sky Survey (eRASS) provides us with 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
optical galaxy samples. We obtained a $9\sigma$ significant detection of the
total X-ray signal from filaments in the 0.3–1.2 keV band. We introduce a
novel method to estimate the contamination fraction from unmasked X-ray halos,
active galactic nuclei, and X-ray binaries that are associated with galaxies
in filaments. We found an approximately 40% contamination fraction by these
unmasked sources, suggesting that the remaining 60% of the signal could be
from the WHIM and a $5.4\sigma$ 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/\mathrm{K})=6.84\pm 0.07$, obtained by using a single
temperature model is marginally higher than the simulation results, which
could be due to the fitting of a single temperature model on a multi-
temperature spectrum. The best-fit baryon density contrast
$\log\Delta_{\mathrm{b}}=1.88\pm 0.18$ generally agrees with the X-ray
emitting phases in the IllustrisTNG simulation, suggesting that the broadband
X-ray emission traces the high temperature and density end of the entire WHIM
population.

###### Key Words.:

large-scale structure of Universe – intergalactic medium – X-rays: diffuse
background

##  1 Introduction

The warm-hot intergalactic medium (WHIM) permeating the cosmic filaments with
temperature of $10^{5-7}$ K and density $n_{\mathrm{H}}<10^{-4}$ cm-3 is 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 that trace the total line-of-sight free electron
column density quantify the cosmic baryon density in the low redshift Universe
in an 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 (BBN) (Pitrou et al. 2018),
which conclude that baryons are not missing in the low redshift Universe.
However, the FRB dispersion measure technique is unable to provide additional
detailed information about the baryons such as density and temperature.
Exploring these properties by using other observation techniques such as X-ray
observations and Sunyaev–Zeldovich (SZ) effects is still one of the ongoing
prime investigations of 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 $\sim 10^{6}$ K phases
WHIM, though 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; 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 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), 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\sigma$ and $3.8\sigma$) of stacked X-ray surface brightness
emission from cosmic filaments using ROSAT All-Sky Survey (RASS) and  Final
Equatorial-Depth Survey (eFEDS) data (Tanimura et al. 2020b, 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, 2020a).
However, in the previous work, the detection in the very soft X-ray band
$<0.6$ keV at the observed frame, which is needed to constrain the temperature
of the redshifted $T<10^{7}$ K gas, is 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  All-Sky Survey (eRASS) data
that covers 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_{0}=70$ km s-1 Mpc-1, $\Omega_{\mathrm{m}}=0.3$,
$\Omega_{\Lambda}=0.7$, and adopt the ratio between baryon density parameter
to matter density parameter from _Planck_ 18 cosmology results, which is
$\Omega_{\mathrm{b}}/\Omega_{\mathrm{m}}=0.158$ (Planck Collaboration et al.
2020). At $z=0$, the mean baryon density is
$\overline{\rho}_{\mathrm{b}}(0)=4.51\times 10^{-31}$ g cm-3. The baryon
density contrast is defined as
$\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. Without further indication, all
sizes/scales/lengths are 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 use. 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. Sect. 5 presents the spectral
stacking method and temperature measurements. The estimation of the detected
X-ray emitting phases and the comparison with numerical simulations will be
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 eROSITA All-Sky Survey (hereafter eRASS:4)
data were collected from 12th December 2019 to 19th December 2021. The data
was 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 select the events from telescope
modules 1,2,3,4,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 adopt an optical filament catalog (Malavasi et al. 2020) compiled using the
DisPerSE software (Sousbie 2011; Sousbie et al. 2011) 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\sigma$ 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 convert 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 filaments is 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 $\sim 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^{\circ}$ are excluded from the analysis
for simplicity. The total analysis footprint before source masking has a solid
angle shown in Fig. 1, amounting to 2275 deg2. In this footprint, we select
7817 filaments in a redshift range $0.2<z<0.6$ and a physical length range 20
Mpc $<l_{\mathrm{fil}}<$ 100 Mpc for stacking analysis (see Fig. 2).

Figure 1: Smoothed 0.3–1.2 keV eROSITA count rate map of the analysis
footprint. The missing pixels show source-masked regions. Figure 2:
Distribution of the selected filaments in the redshift and physical length
space.

###  2.1 Source masking

We mask out all eRASS:4 X-ray sources in the 0.2–2.3 keV eSASS detection
catalog111Proprietary eROSITA data. Version
all_s4_SourceCat1B_221031_poscorr_mpe_photom. The detection pipeline used in
this work is the same as the one version used for the first release of the
All-Sky Survey data (Merloni et al. 2024). The catalog contains point-like
sources and extended sources, such as supernova remnants, galaxies, galaxy
groups, and clusters. To be conservative, we do not apply selection cuts on
the contaminating X-ray emitting sources.

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

$S(r)=S_{0}\left(1+\frac{r^{2}}{r_{\mathrm{PSF}}^{2}+r_{\mathrm{EXT}}^{2}}\right)^{-1.5},$
(1)

where $r_{\mathrm{PSF}}=9^{\prime\prime}$ is the survey averaged point spread
function (PSF) core radius parameter obtained from point source stacking and
$r_{\mathrm{EXT}}$ is the parameter EXT given by the catalog if the source is
extended. We adopt the radius at which the source profile reaches 0.1 (10%) of
the background level as the masking radius, i.e.,
$S(r_{\mathrm{cut}})=0.1\times S_{\mathrm{bkg}}$, where $S_{\mathrm{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 mask out
all sources in the eRASS1 galaxy cluster and group catalog with radii of
$1.5\times r_{500}$. Additionally, we adopt 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 exclude sources with richness $\lambda>20$ and redshift $0.2<z<0.8$
in radii of $1.5\times r_{\lambda}$, where $r_{\lambda}=(\lambda/100)^{1/5}$
$h^{-1}$ Mpc, is the richness scaled radius from the redMaPPer algorithm. This
cut safely removes clusters of galaxies with both 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 adopt a galaxy group catalog (Tinker 2021) based on
the SDSS Main Galaxy Sample (Strauss et al. 2002). We exclude sources with
$n_{\mathrm{sat}}>1$ in radii of $1.5\times r_{500}$. If a source has
$1.5\times r_{500}<400$ kpc, we use a fixed value of 400 kpc instead. In
addition to those cluster and group catalogs, we mask out critical points
identified by DisPerSE of type maximum and bifurcation in the filament catalog
in a redshift range $0.2<z<0.6$ with fixed radii of 2 Mpc as we consider these
to be the nodes of the large scale structure (see Malavasi et al. 2020, for
the details). The details of all masked sources are listed in table 1.

We note that we did not mask optically selected clusters with $\lambda<20$, in
which richness range the identified clusters or groups usually 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
$\lambda<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 estimate 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 in Sect. 4.

Table 1: Summary of masked sources. Catalog | Type | Selection | Masking Radius
---|---|---|---
eRASS:4 0.2–2.3 keV detection | X-ray | all | $r|_{S=0.1S_{\mathrm{bkg}}}$
eRASS1 cluster & group | X-ray | all | $1.5\times r_{\mathrm{500c}}$
redMaPPer on LS DR9 & DR10 | Optical |  $\lambda>20$, $0.2<z<0.8$ | $1.5\times r_{\lambda}$
SDSS MGS group | Optical | $n_{\mathrm{sat}}>1$ | $\max\\{400\mathrm{kpc},1.5\times r_{\mathrm{500c}}\\}$
Critical points in filament catalog | DisPerSE | maxima + bifurcations | 2 Mpc

###  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\sim 0.47$ (redshifted), assuming a $kT<1$ keV gas temperature. We
split the analysis footprint into 27 tiles, each size of 10 deg $\times$ 10
deg. For each tile, we merge the event list using evtool and obtain flare-free
time interval using flaregti with parameters pimin=5000, threshold=1.1,
source_size=150. We use 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^{\prime\prime}$. We use expmap to create
a corresponding vignetting corrected exposure map for each count image. After
masking sources, we bin the tile image pixels to NSIDE=1024
HEALPix222http://healpix.sourceforge.net (Górski et al. 2005; Zonca et al.
2019) map pixels. The NSIDE=1024 HEALPix has a $3.4^{\prime}$ pixel size
corresponding to 1.4 Mpc at $z=0.6$. For the count map, the count number of
the $i$ th pixel, $C_{\mathrm{HP},i}=\sum_{j}C_{ij}$, where $C_{ij}$ is the
count number of the $j$ th image pixel that is located in the $i$ th HEALPix
pixel. For the exposure map, the exposure time of the $i$ th pixel,
$t_{\mathrm{HP},i}=(\sum_{j}t_{ij}\times\Omega_{ij})/{\Omega_{\mathrm{HP},i}}$,
where $\sum_{j}t_{ij}$ and $\Omega_{ij}$ are the exposure and solid angle of
the $j$ th unmasked image pixel that is located in the $i$ th HEALPix pixel,
and $\Omega_{\mathrm{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_{\mathrm{HP}}$ HEALPix
pixels, $n_{\mathrm{fil}}$ filaments and $n_{\mathrm{bin}}$ distance bins. For
each combination of the $i$ th pixel $\mathcal{P}_{i}$, the $j$ th filament
$\mathcal{F}_{j}$ and the $k$ th distance bin $\mathcal{B}_{k}$, we assign a
Boolean coefficient $b_{ijk}$, whose value is 1 (or 0) if $\mathcal{P}_{i}$ is
in (or not in) the $k$ th distance bin of $\mathcal{F}_{j}$. Therefore, for
each filament $\mathcal{F}_{j}$, the count number in the $k$ th transverse
distance bin is $C_{jk}=\sum_{i}^{n_{\mathrm{HP}}}{b_{ijk}\times C_{i}}$, and
the exposure time in the $k$ th transverse distance bin is
$t_{jk}=\sum_{i}^{n_{\mathrm{HP}}}{b_{ijk}\times t_{i}}$. We use exposure-
weighted averages to obtain the stacked surface brightness profile. The
stacked profile

$S_{\mathrm{X}}(r_{k})=\frac{1}{\sum_{j}^{n_{\mathrm{fil}}}{t_{jk}}}\sum_{j}^{n_{\mathrm{fil}}}t_{jk}\times\frac{C_{jk}}{t_{jk}}=\frac{\sum_{j}^{n_{\mathrm{fil}}}C_{jk}}{\sum_{j}^{n_{\mathrm{fil}}}t_{jk}},$
(2)

is the quotient between the sum of the count profile and the sum of the
exposure profile. We use the averaged surface brightness in the outer part of
the stacked profile as the local background $S_{\mathrm{X,bkg}}$ and subtract
it from the entire profile to get a net surface brightness profile as

$S_{\mathrm{X,net}}(r_{k})=S_{\mathrm{X}}(r_{k})-S_{\mathrm{X,bkg}}.$ (3)

Due to the projection of filaments at different redshifts, on average, each
pixel is shared by 12 filaments if the profile is extracted up to 20 Mpc. It
means that, for a single pixel, distance bins of different filaments at
different redshifts could be simultaneously located at its position.
Therefore, the surface brightness in different bins is possibly correlated. We
apply a bootstrapping method to construct resampled pixel sets to obtain the
estimates of the mean and covariance of the stacked profile. We use $l$ to
denote the index of the randomized pixel sample. The $l$ th resample of the
pixel set with replacement is

$\boldsymbol{\mathcal{P}}_{l}^{{}^{\prime}}=\left(\mathcal{P}_{l1}^{{}^{\prime}}\quad\mathcal{P}_{l2}^{{}^{\prime}}\quad\dots\quad\mathcal{P}_{ln_{\mathrm{HP}}}^{{}^{\prime}}\right),\text{
with }\mathcal{P}_{li}^{{}^{\prime}}\in\boldsymbol{\mathcal{P}}.$ (4)

Following equations 2 and 3, for each resampled pixel set
$\boldsymbol{\mathcal{P}}_{l}^{{}^{\prime}}$, we obtain a stacked surface
brightness profile $\boldsymbol{S}_{\mathrm{X,net},l}^{{}^{\prime}}(r_{k})$.
After drawing $n_{\mathrm{bs}}$ random resamples, the bootstrapping estimate
of the net surface brightness profile is

$\hat{S}_{\mathrm{X,net}}(r_{k})=\frac{1}{n_{\mathrm{bs}}-1}\sum_{l}^{n_{\mathrm{bs}}}S_{\mathrm{X,net},l}^{{}^{\prime}}(r_{k}),$
(5)

and the bootstrapping estimate of the covariance between the $k_{1}$ th and
$k_{2}$ th bin is

$\displaystyle\mathcal{C}_{k_{1}k_{2}}=$
$\displaystyle\frac{1}{n_{\mathrm{bs}}-1}$
$\displaystyle\times\sum_{l}^{n_{\mathrm{bs}}}\left[S_{\mathrm{X,net},l}^{{}^{\prime}}(k_{1})-\hat{S}_{\mathrm{X,net}}(k_{1})\right]\left[S_{\mathrm{X,net},l}^{{}^{\prime}}(k_{2})-\hat{S}_{\mathrm{X,net}}(k_{2})\right].$
(6)

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

###  3.1 Stacking the surface brightness profile

We applied our stacking method to the HEALPix count and exposure maps created
in Sect. 2.2. Fig. 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 $\sim$10 Mpc
in the transverse distance. We use the 10–20 Mpc distance region to evaluate
the local background level and subtract 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
construct 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\sigma$ confidence level in the 0–10 Mpc bins.

To investigate the dependence of the significance of the stacked profile on
the selection of the energy band and filament lengths, we apply our analysis
of two sub-bands and two distinct filament populations with different lengths
and find consistent results with our analysis with default selection (see
Appendix A for details).

Figure 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\sigma$ uncertainty
estimated from a pixel-wise bootstrapping method, and that of the control
samples is the $1\sigma$ scatter. The correlation matrix of the 5 bins in the
0–10 Mpc range is plotted in the inset figure.

###  3.2 Stacking method validation

Filaments at different redshifts could be overlapped with each other due to
projection effects. 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), and
the emission in the overlapping source region will be double counted and
therefore increase the stacked signal. Meanwhile, our stacking method assumes
that systematic uncertainties from the spatial variations of multiple
components, e.g. the Milky Way foreground emission and absorption, the Cosmic
X-ray background (CXB) emission, as well as the inhomogeneity of exposure is
minimized by a large number of filaments randomly distributed in the
footprint.

To evaluate the impact of the second case overlapping and multiple possible
systematic uncertainties, we first simulate mock count maps and then apply 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 mock map creation are presented in Appendix C. We simulate 200 mock count
maps and apply the same profile stacking method to all. The stacked mock
profiles from the mock maps are plotted in Fig. 5, together with the injected
model profile. The averaged mock profile in general agrees with the injected
model within the scatter and is slightly higher than the injected model by
$13\%$. This $13\%$ enhanced signal could be due to the second case of
overlapping, and we correct this effect when modeling gas properties in the
subsequent sections. The validation also suggests although the averaged
filament signal at the spine is $0.5\%$ of the total background level, the
emission profile can be successfully recovered by our stacking method.

Figure 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, which is evaluated in Sect. 3.2. Figure 5: The results of 200 net
surface brightness profiles of the simulated maps 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 $\sim 0.5\%$ of the background value, the stacking method can
successfully recover the injected 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 $\lambda>20$. We note
that X-ray sources outside these selection criteria (i.e., X-ray sources 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
evaluate 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 luminous active galactic nuclei, and X-ray binaries. The steps
of modeling the contribution from these sources in the stacked total signal
are:

  1. 1.

We build up a redshift dependent $L_{\mathrm{X}}-M_{*}$ scaling relation of
each type of the sources in a $9.5<\log(M_{*}/M_{\odot})<12$ stellar mass
range. We describe the contribution from each type of unmasked source in
detail in Sect. 4.1.

  2. 2.

We adopt 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 select galaxies in the redshift range of [0.15, 0.7]
and clean up this catalog by removing sources with redshift uncertainty
$>20\%$. We predict 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
position of the galaxy. We further correct for imcompleteness of the observed
stellar mass function (SMF) in 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_{\sun}$ 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. 3.

For each selected galaxy, we add the predicted count rate to the HEALPix pixel
that the galaxy is located in. By this method, we create a predicted count
rate map of unmasked galaxy X-ray emission in the selected redshift range. We
apply the same source mask used in the profile extracting and stacking
pipeline to extract the profile of unmasked signals.

Figure 6: Rest-frame galaxy averaged 0.5–2 keV luminosity - stellar mass
scaling relations for unmasked X-ray halos (_top left_), AGNs (_top right_),
XRBs (_bottom left_) at different redshifts. The scaling relations of total
unmasked X-ray signals associated with galaxies are plotted in the _bottom
right_ panel, together with the recent eRASS stacking results of $z<0.2$
central galaxy sample in black.

###  4.1 Properties of contaminating unmasked sources

To estimate the flux distribution of the contaminating signal, we construct
$L_{\mathrm{X}}-M_{*}$ for three types of unmasked X-ray sources. Meanwhile,
we also investigate the spectral properties of these to be used in count rate
conversion in spectral fitting. The details of the modeling are given in the
subsections below for each component.

####  4.1.1 Unresolved haloes

We use the X-ray halo $L_{\mathrm{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 exclude
eRASS:4 detected X-ray sources and clusters detected in the optical band with
richness higher than 20 before stacking, to obtain the $L_{\mathrm{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), and the X-ray
properties of each halo are painted by a phenomenological modeling approach
(Comparat et al. 2020) and then scaled to agree the observed
$L_{\mathrm{X}}-M_{*}$ relation from RASS galaxy stacking results of Anderson
et al. (2015). By applying the same source selection criteria, we tabulate the
averaged rest-frame $L_{\mathrm{0.5-2keV}}-M_{*}$ relation for unmasked X-ray
halos in different redshift bins.

The galaxy sample used in Anderson et al. (2015) is 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_{\mathrm{X}}-M_{*}$ relation has a steep slope of
$\gtrsim 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 obtain the discrepancy between the two
stellar mass schemes. We use this stellar mass discrepancy to empirically
modify the tabulated $L_{\mathrm{0.5-2keV}}-M_{*}$ relation. The modified
X-ray halo $L_{\mathrm{X}}-M_{*}$ relation is plotted in the top left panel of
Fig. 6).

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 necessary. For this
purpose, we assume that each source can be represented by a single temperature
collisional ionization equilibrium (CIE) model and estimate the temperature
using the luminosity - temperature ($L_{\mathrm{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_{\mathrm{X}}-kT$ scaling relation, we also 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. 17). The spectral model
on each single temperature is calculated using Astrophysical Plasma Emission
Code (APEC) with $Z=0.4Z_{\sun}$. Then we use the redshift and length
distributions of the selected filament sample to create the synthesis spectral
model as,

$F(E)=\sum_{i}\frac{l_{i}\times F(E,z_{i})}{D_{\mathrm{L}}(z_{i})^{2}},$ (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}$,
$D_{\mathrm{L}}$ is the luminosity distance and is a function of redshift. The
constructed model spectrum is shown in the right panel of Fig. 17.

####  4.1.2 Unresolved AGN

We calculate the galaxy averaged $L_{\mathrm{X}}-M_{*}$ relation for
unresolved AGN 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 apply a selection criteria that removes AGN with
$F_{\mathrm{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 redshift-stellar mass
bin, we sum 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 will
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 investigate the stacked spectral properties of unresolved AGNs. We
adopt the same spectral model used for creating 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$\alpha$ emission line, a soft-X-ray-excess, and a
cold-reflection component, respectively. For each unmasked source, based on
the hard band luminosity and the intrinsic absorbing column density, we
predict a corresponding mock spectrum. The stacked spectra of different bins
as well as the median spectrum are plotted in Fig. 18. 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 use a stellar mass and star forming rate (SFR) dependent model of
integrated X-ray luminosity of XRB 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) to construct the averaged
$L_{\mathrm{X}}-M_{*}$ relations of XRB at different redshifts. The results
are plotted in the bottom left panel of Fig. 6, where the relations show
redshift evolution at $M_{*}<10^{11}$ $M_{\odot}$, which is the result of the
redshift evolution of star-forming rate and star-forming galaxy fraction.

We adopt a power-law distribution with intrinsic absorption spectral model for
flux conversion, with parameters $\Gamma=1.7$ and $n_{\mathrm{H}}=2\times
10^{21}$ cm-2 suggested by observations (Lehmer et al. 2019).

####  4.1.4 Comparison with observed $L_{\mathrm{X}}-M_{*}$ relations

Recently, Zhang et al. (2024) investigated the $L_{\mathrm{X}}-M_{*}$ scaling
relation of galaxies using eRASS:4 data. We compare our constructed total
unresolved galactic X-ray sources and their $L_{\mathrm{X}}-M_{*}$ relations
with their results from eROSITA observations. We note that the stellar mass
used in Zhang et al. (2024) is adopted from Tinker (2021), which was
originally estimated by Chen et al. (2012). Similar to the analysis in Sect.
4.1.1, we scale the observed $L_{\mathrm{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 $\sim 0.1$. Our constructed
relation at $z=0.2$ has the same slope at different stellar masses as the 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 calculate 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 divide the stellar mass range into two subranges of
$9.5<\log M_{*}(M_{\sun})<11$ and $11<\log M_{*}(M_{\sun})<12$, and repeat the
same analysis. The results of fractions are listed in Table 2. We find that 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 AGN and XRB have more
contribution. We use a bootstrapping method to estimate the statistical
uncertainties of the fraction. In short, we randomly select galaxies with
replacement to construct bootstrapping samples, and calculate 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\%$.

Because the uncertainty on photometric redshifts that we used in this analysis
is not negligible and is a function of magnitude (it can reach a few tens
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$, and that of the galaxies we used to model
unmasked source contribution is $0.15<z<0.7$. It means that for a galaxy whose
true redshift $z_{\mathrm{true}}$ is within the filament redshift range, but
the measured redshift $z_{\mathrm{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_{\mathrm{true}}$ outside of the
filament redshift range but $z_{\mathrm{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), $\Delta
z_{\mathrm{norm}}\equiv(z_{\mathrm{phot}}-z_{\mathrm{spec}})/(1+z_{\mathrm{spec}})$,
and assume that $z_{\mathrm{spec}}=z_{\mathrm{true}}$. The ratio between the
recovered flux and the true flux of a galaxy is

$f=\frac{D_{\mathrm{L}}^{2}(z_{\mathrm{true}})}{n_{\mathrm{MC}}}\sum_{i}{\frac{w_{i}}{D_{\mathrm{L}}^{2}(z_{i})}},$
(8)

where $n_{\mathrm{MC}}$ is the number of realizations, $z_{i}$ is the $i$ th
randomized redshift value, the value of $w_{i}$ is one if $0.15<z_{i}<0.7$ and
is zero elsewhere. 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-\Delta z$ relation shown in Fig. 5 of Zou et al.
(2019), we predict the recovered flux fractions for unmasked galaxies in the
$9.5<\log M_{*}(M_{\sun})<11$ and $11<\log M_{*}(M_{\sun})<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, which is due to the high quality of photometric redshift of bright
central galaxies. This fraction is $\sim 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
$\sim 40\%$.

The $\sim 40\%$ fraction of the unmasked source contribution suggests that the
remaining $\sim 60\%$ of the stacked signal is from the WHIM in the filaments,
which reach a $5.4\sigma$ detection significance. We note that 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. Future better understanding of the all different types of X-ray
sources in the Universe will bring us more accurate results of the
contamination fraction.

Figure 7: The 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 $\sim 60\%$ based on the unmasked signal estimation.  Table 2: Fraction of X-ray signal from unmasked sources in different stellar mass bins. $\log(M_{*}/M_{\sun})$ | Halo | AGN | XRB
---|---|---|---
9.5–12 | 0.20 | 0.12 | 0.05
9.5–11 | 0.01 | 0.09 | 0.04
11–12 | 0.19 | 0.03 | 0.01
Figure 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_{\sun})<12$ and $9.5<\log M_{*}(M_{\sun})<11$
stellar mass bins, respectively. For high stellar mass galaxies that
contribute most of the unmasked signals, the impact of photometric redshift
uncertainty is only a few percent.

##  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
in the data, we stack 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.
We study the spectral properties and constrain the temperature of the detected
WHIM using the stacked filament spectra through this analysis.

###  5.1 Stacking processes

For each filament, we first merge the event files of TM8 and shift the Pulse
Invariant of each event by a factor of $1+z$. Then, we use 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, and whose region masks for spectral extraction are created by
constructing images using and excluding sources in Table 1. For each ARF, we
first correct it by multiplying the foreground absorption curve given the
averaged Galactic $n_{\mathrm{H}}$ value (HI4PI Collaboration et al. 2016) in
the extraction region and the TBABS absorption model (Wilms et al. 2000), and
then shift it by a factor of $1+z$ to the source rest frame. The response
matrix file (RMF) is taken from the calibration database and then shifted to
the source rest frame by adopting a two-dimensional interpolation strategy on
both axes of induced photon energy and resulting spectral channels (see Bulbul
et al. 2014, for further details).

We use tools addspec, addarf, and addrmf in
FTOOLS333http://heasarc.gsfc.nasa.gov/ftools (Nasa High Energy Astrophysics
Science Archive Research Center (2014) Heasarc) 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
products444See
https://erosita.mpe.mpg.de/edr/DataAnalysis/srctool_doc.html#Output_files for
the details of the srctool products.. In short, the EXPOSURE value of a
spectrum $t_{\mathrm{EXP}}$ is the total time of an extraction region on a
detector, instead of the averaged exposure time in the region
$t_{\mathrm{ave}}$; BACKSCAL value $\Omega_{\mathrm{BS}}$ is the averaged
intersection of the extraction region and the detector area; the geometric
solid angle of the extraction region $\Omega_{\mathrm{REG}}$ has a keyword
REGAREA, and
$\Omega_{\mathrm{REG}}/\Omega_{\mathrm{BS}}=t_{\mathrm{EXP}}/t_{\mathrm{ave}}$.
Furthermore, the normalization of an ARF curve is scaled by a factor of
$t_{\mathrm{ave}}/t_{\mathrm{EXP}}$. Therefore, we add up all individual
source spectra to obtain a stacked source spectrum.

$N_{\mathrm{stacked,src}}(E)=\sum_{i}N_{\mathrm{src},i}(E),$ (9)

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

$N_{\mathrm{stacked,bkg}}(E)=\sum_{i}\frac{\Omega_{\mathrm{REG,src},i}}{\Omega_{\mathrm{REG,bkg},i}}\times
N_{\mathrm{bkg},i}(E).$ (10)

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

$A_{\mathrm{stacked}}(E)=\frac{1}{\overline{t}_{\mathrm{ave}}}\sum_{i}t_{\mathrm{EXP},i}\Omega_{\mathrm{REG},i}A_{i}(E),$
(11)

which 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 $\Omega_{\mathrm{REG}}\times t_{\mathrm{ave}}$,
which is

$\boldsymbol{\mathcal{R}}_{\mathrm{stacked}}=\frac{1}{\sum_{i}\Omega_{\mathrm{BS},i}t_{\mathrm{EXP},i}}\times\sum_{i}\Omega_{\mathrm{BS},i}t_{\mathrm{EXP},i}\boldsymbol{\mathcal{R}}_{i}.$
(12)

The blueshifted response files at the rest frame corresponding to each
filament are then later stacked using FTOOLS addarf and addrmf to obtain one
final stacked response file. These files are used in the modeling of the
stacked spectrum described in the section below.

One clear advantage of the spectral stacking method in the rest frame is the
reduction of 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/background features, which are frozen in the detector frame, are
smoothed over the redshift range covered by the sample, therefore reducing
their overall effect on the net spectrum. The spectral stacking method is
particularly 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 5 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.

We use XSPEC version 12.13.0 and the package
BXA555https://github.com/JohannesBuchner/BXA (Buchner et al. 2014) for
spectral analysis. The spectrum is first optimally binned (Kaastra & Bleeker
2016) using the tool ftgrouppha and then loaded into XSPEC. We use the 0.5–6.0
keV range for analysis and adopt 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 has been described in Sect. 4.1. We calculate 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 use a single temperature (1-T) APEC
model with fixed abundance $Z=0.2Z_{\sun}$. This value is based on the
reported $\lesssim 0.3Z_{\odot}$ abundance in cluster outskirts by both
observations and numerical simulations (Bulbul et al. 2016; Mernier et al.
2017; Biffi et al. 2017) and $>0.1Z_{\odot}$ 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\pm 0.10$ keV and
a fitting statistics $c_{\mathrm{stat}}/\mathrm{d.o.f.}=105/111$. We also
tried to fix the metalicity to $0.1Z_{\sun}$ and $0.3Z_{\sun}$ to test the
impact of metalicity assumption and found there is has no impact on the best-
fit temperature.

Recently, Vijayan & Li (2022) demonstrates that fitting a multi-temperature
nature spectrum using a 1-T model will lead to an overestimation of 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^{6}$ 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 $\sim 0.56$ keV O VII lines emitted by cooler
phases of gas around $10^{6.3}$ K ($\sim 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 try to fit the spectrum using a multi-temperature model. We adopt a
log-normal temperature distribution DEM model with a peak at
$\log(T/\mathrm{K})=5.9$ and $1\sigma$ width of 0.45 dex, which is suggested
by our analysis of numerical simulations (see Sect. 6.2 later). We use a
simplified assumption of CIE condition and $0.2Z_{\sun}$ 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. It could be due to that
the true temperature distribution is more complex than the log-normal
distribution we use, which creates discrepancy in some energy ranges.

Figure 9: _Left_ : 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. The spectrum is binned to a minimum signal-to-noise ratio of 2.5 for
illustration. 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\pm 0.10$ keV. _Right_ : Same as the left
panel, but fitted by a log-normal temperature distribution with
$\log(\mu/\mathrm{K})=5.9$ and $\sigma=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 $\sim
0.7$ keV.

##  6 The detected phases of WHIM

###  6.1 Surface brightness profile modeling

Table 3: Parameters and the corresponding priors for constraining the density profile of the observed phase WHIM. Parameter | Unit | Description | Range | Prior
---|---|---|---|---
$\log\Delta_{\mathrm{b,0}}$ |  | Central baryon density contrast in the logarithmic scale | [0.5,3.5] | $\mathcal{N}(1.5,0.5)$
$\log f_{\mathrm{vol,0}}$ |  | Central volume fraction in the logarithmic scale | [-3, -1] | $\mathcal{N}(-2,0.3)$
$\log r_{\mathrm{c,n}}$ | Mpc | Core radius of the gas density profile in the logarithmic scale | [-1, 0.8] | $\mathcal{N}(0.3,0.2)$
$\alpha_{\mathrm{n}}$ |  | Slope of the density profile in the outskirts | [-3,0] | $\mathcal{N}(-1,0.3)$
$\alpha_{\mathrm{fvol}}$ |  | Slope of the volume fraction profile in the outskirts | [-3,0] | $\mathcal{N}(-1,0.3)$
$\log\sigma_{\mathrm{sm}}$ | Mpc | Gaussian smoothing kernel width in logarithmic scale | [0,0.8] | $\mathcal{N}(0.3,0.2)$

With the best-fit temperature obtained in the previous section, we fit the
density profiles of the warm gas in cosmic filaments. In this section, we use
$r$ to denote the physical radius of gas cells/particles to the filament spine
in the 3-D space. Based on results from numerical simulations (Cui et al.
2019; Galárraga-Espinosa et al. 2021; Tuominen 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 introduce a parameter $f_{\mathrm{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

$\bar{\epsilon}_{\mathrm{det}}(r)=1.20\times
n_{\mathrm{H}}(r)^{2}f_{\mathrm{vol}}(r),$ (13)

where $n_{\mathrm{H}}$ is the hydrogen density, 1.20 is the
$n_{\mathrm{e}}/n_{\mathrm{H}}$ ratio for the plasma at the detected
temperature. We use a King-type profile (King 1962), which is widely used to
fit central flattened profiles, to model the radial change of both
$n_{\mathrm{H}}$ and $f_{\mathrm{vol}}$,

$\displaystyle n_{\mathrm{H}}(r)=$ $\displaystyle
n_{\mathrm{H,0}}\left[1+(r/r_{\mathrm{c,n}})^{2}\right]^{\alpha_{\mathrm{n}}};$
(14) $\displaystyle f_{\mathrm{vol}}(r)=$ $\displaystyle
f_{\mathrm{vol,0}}\left[1+(r/r_{\mathrm{c,fvol}})^{2}\right]^{\alpha_{\mathrm{fvol}}},$
(15)

where, $n_{\mathrm{H,0}}$ ($f_{\mathrm{vol,0}}$), $r_{\mathrm{c,n}}$
($r_{\mathrm{c,fvol}}$), and $\alpha_{\mathrm{n}}$ ($\alpha_{\mathrm{fvol}}$)
are the central normalization, core radius, and slope in the outskirts for
$n_{\mathrm{H}}$ (or $f_{\mathrm{vol}}$), respectively. Because the radial
dependencies of $n_{\mathrm{H}}$ and $f_{\mathrm{vol}}$ have a strong
degeneracy, we assume that they have the same core radii, i.e.,
$r_{\mathrm{c,n}}=r_{\mathrm{c,fvol}}$.

In addition, we make several assumptions. First, we assume a constant baryon
density contrast $\Delta_{\mathrm{b,0}}$ at the filament spine for filaments
at all redshifts. In this case,

$\displaystyle
n_{\mathrm{H,0}}(z)=\frac{\Delta_{\mathrm{b,0}}\rho_{\mathrm{crit}}\Omega_{\mathrm{b}}}{1.43\times
m_{\mathrm{H}}}\times\left(1+z\right)^{3},$ (16)

where $\rho_{\mathrm{crit}}$ is the critical density at $z=0$, the factor 1.43
is the ratio between total mass density and hydrogen mass density,
$m_{\mathrm{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_{\mathrm{out}}$, beyond which the gas density of our detected
phase is zero. In this work, we adopt 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_{\mathrm{proj}}$ is

$S_{\mathrm{X,model,conv}}(r_{\mathrm{proj}},z)=G(\sigma_{\mathrm{sm}})*S_{\mathrm{X,model}}(r_{\mathrm{proj}},z),$
(17)

where $G(\sigma_{\mathrm{sm}})$ is a Gaussian convolving kernel that accounts
for the smoothing effect of the observed filament width. The surface
brightness is

$\displaystyle S_{\mathrm{X,model}}(r_{\mathrm{proj}},z)=$ $\displaystyle
A_{\mathrm{CF}}^{-1}K_{\mathrm{ph}}(1+z)^{-3}$
$\displaystyle\times\int_{-l_{\mathrm{out}}}^{l_{\mathrm{out}}}\bar{\epsilon}_{\mathrm{det}}\left(\sqrt{l^{2}+r_{\mathrm{proj}}^{2}},z\right)\Lambda\mathrm{d}l,$
(18)

where $l$ is depth along the line of sight; $K_{\mathrm{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; $\Lambda$ is the
cooling function calculated using the APEC model given the temperature,
metallicity and energy band; $A_{\mathrm{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 construct 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

$S_{\mathrm{X,model,stack}}=\frac{1}{N}\sum_{i}^{N}\Omega_{i}\times
S_{\mathrm{X,model,conv},i}.$ (19)

The fitting statistics is defined as

$\chi^{2}={\boldsymbol{S}_{\mathrm{X,diff}}}^{\intercal}\boldsymbol{\mathcal{C}}^{-1}\boldsymbol{S}_{\mathrm{X,diff}},$
(20)

where
$\boldsymbol{S}_{\mathrm{X,diff}}=\boldsymbol{S}_{\mathrm{X,net,corr}}-\boldsymbol{S}_{\mathrm{X,unmasked,net}}-\boldsymbol{S}_{\mathrm{X,model,stack}}\times
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 $\boldsymbol{\mathcal{C}}$ is the covariance
matrix estimated from the bootstrapping method.

We use a Bayesian inference approach to obtain posterior distributions of the
six parameters, which are $\log\Delta_{\mathrm{b,0}}$, $\log
f_{\mathrm{vol,0}}$, $\log r_{\mathrm{c,n}}$, $\alpha_{\mathrm{n}}$,
$\alpha_{\mathrm{fvol}}$ and $\log\sigma_{\mathrm{sm}}$, as listed in Table 3.
In addition to the fitting statistics $\chi^{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
$\Delta_{\mathrm{b,0}}$ and $f_{\mathrm{vol,0}}$. The predicted Compton-$y$
strength is

$y_{\mathrm{pd}}(r_{\mathrm{proj}})=\left[\frac{\sigma_{\mathrm{T}}}{m_{\mathrm{e}}c^{2}}\int_{-l_{\mathrm{out}}}^{l_{\mathrm{out}}}kTn_{\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],$
(21)

where $\sigma_{\mathrm{T}}$ is the Thomson scattering cross section,
$m_{\mathrm{e}}$ is the electron mass, $c$ is the speed of light. The profiles
of Tanimura et al. (2020a) 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,

$y_{\mathrm{pd,conv}}=G(\sigma_{\mathrm{sm}})*y_{\mathrm{pd}}.$ (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_{\mathrm{pd,conv}}(0)$ to be -7.5, which is approximately the stacked SZ
signal at filament spine reported by (Tanimura et al. 2020a). We use a
likelihood term $P(\log
y_{\mathrm{pd,conv}}(0))=\mathcal{N}(\mu_{\mathrm{logy}},\sigma_{\mathrm{logy}})$
as the addition constraint. We do not exactly know how much of the total
stacked Compton-$y$ is contributed by the X-ray emitting gas, nevertheless, we
adopt $\mu_{\mathrm{logy}}=-8.3$ and $\sigma_{\mathrm{logy}}=0.5$ in our
analysis. The combined likelihood function is

$\ln\mathcal{L}=\cases{-}\frac{1}{2}\chi^{2}-\frac{1}{2}\frac{\left[y_{\mathrm{pd,conv}}(0)-\mu_{\mathrm{logy}}\right]^{2}}{\sigma_{\mathrm{logy}}^{2}}&\text{when\
}\log y_{\mathrm{pd,conv}}(0)<=-7.5,\\\ -\infty\text{elsewhere.}{}$ (23)

We use the Markov Chain Monte Carlo (MCMC) method and the python package
Zeus666https://zeus-mcmc.readthedocs.io/en/latest/ (Karamanis et al. 2021) to
perform Bayesian inference. The marginalized posterior distributions shows
that the $\log\Delta_{\mathrm{b,0}}=2.28\pm 0.22$ and $\log
f_{\mathrm{vol,0}}=-1.91\pm 0.28$. By using all MCMC chains for calculating
the averaged baryon density contrast, we obtain
$\log\bar{\Delta}_{\mathrm{b}}=1.88\pm 0.18$, which corresponds to a hydrogen
density $4.4^{+2.2}_{-1.5}\times 10^{-5}$ cm-3 at the median redshift of the
selected filament sample $z=0.47$.

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

###  6.2 X-ray emitting phases in numerical simulations

To compare the measured physical quantities with those in numerical
simulations, we use the public data (Nelson et al. 2019) of IllustrisTNG
simulations (Nelson et al. 2018; Pillepich et al. 2018; Naiman et al. 2018;
Springel et al. 2018; Marinacci et al. 2018) and select the $z=0.48$ snapshot
of the TNG300-1 box. We aim to create phase diagrams of the filament gas. We
use the DisPerSE algorithm to construct the filaments in that snapshot using
$9<\log(M_{*}/M_{\sun})<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_{\mathrm{fil}}$, to the nearest maxima or bifurcation critical point
$d_{\mathrm{CP}}$, and the nearest critical point but along the direction of
filament $d_{\mathrm{skel}}$ follow the definition in Malavasi et al. (2022).

To reduce the computational time, we randomly select 1/500 of the total gas
cells, which is enough to represent the phase diagram (Galárraga-Espinosa et
al. 2021). We select gas cells with $d_{\mathrm{fil}}<1$ Mpc and
$d_{\mathrm{CP}}>2$ Mpc, and remove 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 recipe777https://www.tng-project.org/data/docs/faq/#gen6. 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^{7}$ 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).

To investigate the traceable populations of the filament gas, we create 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, 2020a), and 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 are different. SZ observations trace $10^{5.5}-10^{6.8}$ K
temperature and a $\Delta_{\mathrm{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 $\Delta_{\mathrm{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 cosmic
ultraviolet background (UVB) (e.g., Nicastro et al. 1999; Krongold et al.
2003), which alters the ionization balance and exhibit different spectral
features from the CIE plasma (Khabibullin & Churazov 2019; Churazov et al.
2023). In this section, we combine the WHIM phase distributions in numerical
simulations (left panel of Fig. 10) and the photoionization emission model to
characterize the impact of including photoionization model on our conclusions.

The ionization parameter is defined as

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

where $F_{\mathrm{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 adopt
the Faucher-Giguère (2020) UVB spectral energy distribution for the
calculation. We obtain a relation of $\log\xi[\mathrm{erg\ s^{-1}\
cm^{2}}]=2.26-\log\Delta_{\mathrm{b}}$ at $z=0.48$. We use 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 $\xi$ and temperature. We note that by
setting the parameter tmod=1, the pion model in SPEX allows calculations 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.

Figure 11: _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_{\sun}$. Towards the high density
(low ionization parameter) end, the radiative cooling reaches the CIE limit
and is simply a function of temperature. When $\Delta_{\mathrm{b}}\lesssim 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.

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 $\Delta_{\mathrm{b}}>2$ ($\log\xi[\mathrm{erg\ s^{-1}\
cm^{2}}]\lesssim 0)$, the radiative cooling reaches the CIE limit and is
simply a function of temperature. When $\Delta_{\mathrm{b}}\lesssim 1$
($\log\xi[\mathrm{erg\ s^{-1}\ cm^{2}}]\gtrsim 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 keeps high for the WHIM
with 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 of 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 are in the highest temperature and density phases among the whole
filament gas population.

##  7 Summary

eRASS provides us with the deepest All-Sky X-ray survey in the soft band,
allowing a more precise investigation of X-ray emission in cosmic filaments
than previous pioneering works (Tanimura et al. 2020b, 2022). The aim of this
work was 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. Comparing to the pioneering works, we adopt a novel method to
estimate the fraction of non-WHIM emission in the stacked signal; improve the
spectral stacking method to minimize the systematics from the instrumental and
sky background features at the observer frame; combine the SZ stacking results
to quantify the volume fraction and density of the X-ray emitting WHIM
simultaneously. We summarize our main results as follows:

  1. 1.

We stacked 0.3–1.2 keV surface brightness profiles of 7817 SDSS optical
filaments with a physical lengths of $20-100$ Mpc and redshifts of 0.2–0.6 in
a 2275 $\deg^{2}$ footprint and detect an X-ray emission excess over the
background at a $9\sigma$ significance.

  2. 2.

We validated our profile stacking method by applying it to simulated count
maps that include 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 13% difference, which could be from the
possible boosting from projectively overlapped filaments in same orientations.

  3. 3.

We carefully estimate the fraction of stacked signal from unmasked sources
including X-ray halos, AGNs and XRBs that is associated with galaxies. We
first paint the observed galaxies in the analysis footprint with X-ray count
rate from $L_{\mathrm{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 types 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 $\sim 60\%$ of the signal is from
WHIM in the filaments and the WHIM detection significance is $5.4\sigma$.

  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 by using a 1-T model is $0.58\pm 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 is able to reduce the residual in
the 0.5–0.6 keV energy range, suggesting the presence of the OVII abundant
temperature phases.

  5. 5.

We modeled the baryon density contrast $\Delta_{\mathrm{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_{\mathrm{vol}}$, which is necessary because only a small fraction of WHIM
emits X-rays that can be observed by eROSITA. We use the SZ stacking results
in literature to break the degeneracy between the $\Delta_{\mathrm{b}}$ and
volume fraction $f_{\mathrm{vol}}$, and obtained an averaged
$\log\Delta_{\mathrm{b}}=1.88\pm 0.18$, suggesting that the observed emission
is from the densest population of the WHIM.

  6. 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 is emitted by phases of
$6.2<\log(T/\mathrm{K})<6.9$ and $1.2<\log\Delta_{\mathrm{b}}<2.2$, which
generally agree with our estimated values.

It is the first time we detect $>5\sigma$ X-ray emission from cosmic filaments
and apply 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 from the WHIM within them, will be significantly improved from
multiple aspects. The 4MOST Cosmology Redshift Survey (Richard et al. 2019)
will bring us a 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 Monte Carlo Physarum
Machine (Elek et al. 2022) will be worked on the new survey data to provide a
more accurately constructed filament catalog. Meanwhile, better understandings
of X-ray properties of galaxy groups, AGNs, and XRBs will help us to better
estimate the contamination fractions of unmasked sources. Moreover, high
spectral resolution X-ray missions such as HUBS (Cui et al. 2020) and LEM
(Kraft et al. 2022) will be able to explore a wider parameter space of the
WHIM properties.

###### Acknowledgements.

We acknowledge Jelle Kaastra for the instruction on using the SPEX Pion model.
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.

## References

  * Aird et al. (2015) Aird, J., Coil, A. L., Georgakakis, A., et al. 2015, MNRAS, 451, 1892
  * Akamatsu et al. (2017) Akamatsu, H., Fujita, Y., Akahori, T., et al. 2017, A&A, 606, A1
  * Anderson et al. (2015) Anderson, M. E., Gaspari, M., White, S. D. M., Wang, W., & Dai, X. 2015, MNRAS, 449, 3806
  * Arnaud (1996) 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
  * Bahar et al. (2022) Bahar, Y. E., Bulbul, E., Clerc, N., et al. 2022, A&A, 661, A7
  * Behroozi et al. (2019) Behroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. 2019, MNRAS, 488, 3143
  * Bernardi et al. (2010) Bernardi, M., Shankar, F., Hyde, J. B., et al. 2010, MNRAS, 404, 2087
  * Biffi et al. (2017) Biffi, V., Planelles, S., Borgani, S., et al. 2017, MNRAS, 468, 531
  * Blanton & Roweis (2007) Blanton, M. R. & Roweis, S. 2007, AJ, 133, 734
  * Blanton et al. (2005) Blanton, M. R., Schlegel, D. J., Strauss, M. A., et al. 2005, AJ, 129, 2562
  * Böhringer et al. (1994) Böhringer, H., Briel, U. G., Schwarz, R. A., et al. 1994, Nature, 368, 828
  * Brunner et al. (2022) Brunner, H., Liu, T., Lamer, G., et al. 2022, A&A, 661, A1
  * Buchner (2016) Buchner, J. 2016, Statistics and Computing, 26, 383
  * Buchner (2019) Buchner, J. 2019, PASP, 131, 108005
  * Buchner (2021) Buchner, J. 2021, The Journal of Open Source Software, 6, 3001
  * Buchner et al. (2014) Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125
  * Bulbul et al. (2024) Bulbul, E., Liu, A., Kluge, M., et al. 2024, arXiv e-prints, arXiv:2402.08452
  * Bulbul et al. (2014) Bulbul, E., Markevitch, M., Foster, A., et al. 2014, ApJ, 789, 13
  * Bulbul et al. (2016) Bulbul, E., Randall, S. W., Bayliss, M., et al. 2016, ApJ, 818, 131
  * Cen & Ostriker (1999) Cen, R. & Ostriker, J. P. 1999, ApJ, 514, 1
  * Chen et al. (2012) Chen, Y.-M., Kauffmann, G., Tremonti, C. A., et al. 2012, MNRAS, 421, 314
  * Churazov et al. (2023) Churazov, E., Khabibullin, I. I., Dolag, K., Lyskova, N., & Sunyaev, R. A. 2023, MNRAS, 523, 1209
  * Comparat et al. (2020) Comparat, J., Eckert, D., Finoguenov, A., et al. 2020, The Open Journal of Astrophysics, 3, 13
  * Comparat et al. (2023) Comparat, J., Luo, W., Merloni, A., et al. 2023, A&A, 673, A122
  * Comparat et al. (2019) Comparat, J., Merloni, A., Salvato, M., et al. 2019, MNRAS, 487, 2005
  * Costanzi et al. (2019) Costanzi, M., Rozo, E., Rykoff, E. S., et al. 2019, MNRAS, 482, 490
  * Cui et al. (2020) Cui, W., Chen, L. B., Gao, B., et al. 2020, Journal of Low Temperature Physics, 199, 502
  * Cui et al. (2019) Cui, W., Knebe, A., Libeskind, N. I., et al. 2019, MNRAS, 485, 2367
  * Danforth et al. (2016) Danforth, C. W., Keeney, B. A., Tilton, E. M., et al. 2016, ApJ, 817, 111
  * de Graaff et al. (2019) de Graaff, A., Cai, Y.-C., Heymans, C., & Peacock, J. A. 2019, A&A, 624, A48
  * Dey et al. (2019) Dey, A., Schlegel, D. J., Lang, D., et al. 2019, AJ, 157, 168
  * Driver (2021) Driver, S. 2021, Nature Astronomy, 5, 852
  * Eckert et al. (2015) Eckert, D., Jauzac, M., Shan, H., et al. 2015, Nature, 528, 105
  * Elek et al. (2022) Elek, O., Burchett, J. N., Prochaska, J. X., & Forbes, A. G. 2022, arXiv e-prints, arXiv:2204.01256
  * Fang et al. (2010) Fang, T., Buote, D. A., Humphrey, P. J., et al. 2010, ApJ, 714, 1715
  * Fang et al. (2022) Fang, T., Mathur, S., & Nicastro, F. 2022, Absorption Studies of the Most Diffuse Gas in the Large-Scale Structure, ed. C. Bambi & A. Santangelo (Singapore: Springer Nature Singapore), 1–43
  * Faucher-Giguère (2020) Faucher-Giguère, C.-A. 2020, MNRAS, 493, 1614
  * Fukugita et al. (1998) Fukugita, M., Hogan, C. J., & Peebles, P. J. E. 1998, ApJ, 503, 518
  * Galárraga-Espinosa et al. (2021) Galárraga-Espinosa, D., Aghanim, N., Langer, M., & Tanimura, H. 2021, A&A, 649, A117
  * Gatuzz et al. (2023) Gatuzz, E., García, J. A., Churazov, E., & Kallman, T. R. 2023, MNRAS, 521, 3098
  * Georgakakis et al. (2017) Georgakakis, A., Aird, J., Schulze, A., et al. 2017, MNRAS, 471, 1976
  * Górski et al. (2005) Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759
  * Grandis et al. (2021) Grandis, S., Mohr, J. J., Costanzi, M., et al. 2021, MNRAS, 504, 1253
  * HI4PI Collaboration et al. (2016) HI4PI Collaboration, Ben Bekhti, N., Flöer, L., et al. 2016, A&A, 594, A116
  * Ishiyama et al. (2021) Ishiyama, T., Prada, F., Klypin, A. A., et al. 2021, MNRAS, 506, 4210
  * Johnson et al. (2019) Johnson, S. D., Mulchaey, J. S., Chen, H.-W., et al. 2019, ApJ, 884, L31
  * Kaastra (2017) Kaastra, J. S. 2017, A&A, 605, A51
  * Kaastra & Bleeker (2016) Kaastra, J. S. & Bleeker, J. A. M. 2016, A&A, 587, A151
  * Kaastra et al. (1996) Kaastra, J. S., Mewe, R., & Nieuwenhuijzen, H. 1996, in UV and X-ray Spectroscopy of Astrophysical and Laboratory Plasmas, 411–414
  * Kaastra et al. (2024) Kaastra, J. S., Raassen, A. J. J., de Plaa, J., & Gu, L. 2024, SPEX X-ray spectral fitting package
  * Karamanis et al. (2021) Karamanis, M., Beutler, F., & Peacock, J. A. 2021, MNRAS, 508, 3589
  * Khabibullin & Churazov (2019) Khabibullin, I. & Churazov, E. 2019, MNRAS, 482, 4972
  * King (1962) King, I. 1962, AJ, 67, 471
  * Kluge et al. (2024) Kluge, M., Comparat, J., Liu, A., et al. 2024, arXiv e-prints, arXiv:2402.08453
  * Klypin et al. (2016) Klypin, A., Yepes, G., Gottlöber, S., Prada, F., & Heß, S. 2016, MNRAS, 457, 4340
  * Kraft et al. (2022) Kraft, R., Markevitch, M., Kilbourne, C., et al. 2022, arXiv e-prints, arXiv:2211.09827
  * Krolik et al. (1981) Krolik, J. H., McKee, C. F., & Tarter, C. B. 1981, ApJ, 249, 422
  * Krongold et al. (2003) Krongold, Y., Nicastro, F., Brickhouse, N. S., et al. 2003, ApJ, 597, 832
  * Lehmer et al. (2019) Lehmer, B. D., Eufrasio, R. T., Tzanavaris, P., et al. 2019, ApJS, 243, 3
  * Li & White (2009) Li, C. & White, S. D. M. 2009, MNRAS, 398, 2177
  * Liu et al. (2022) Liu, A., Bulbul, E., Ghirardini, V., et al. 2022, A&A, 661, A2
  * Liu et al. (2024) Liu, A., Bulbul, E., Kluge, M., et al. 2024, A&A, 683, A130
  * Liu et al. (2017) Liu, T., Tozzi, P., Wang, J.-X., et al. 2017, ApJS, 232, 8
  * Macquart et al. (2020) Macquart, J. P., Prochaska, J. X., McQuinn, M., et al. 2020, Nature, 581, 391
  * Malavasi et al. (2020) Malavasi, N., Aghanim, N., Douspis, M., Tanimura, H., & Bonjean, V. 2020, A&A, 642, A19
  * Malavasi et al. (2022) Malavasi, N., Langer, M., Aghanim, N., Galárraga-Espinosa, D., & Gouin, C. 2022, A&A, 658, A113
  * Marinacci et al. (2018) Marinacci, F., Vogelsberger, M., Pakmor, R., et al. 2018, MNRAS, 480, 5113
  * Martizzi et al. (2019) Martizzi, D., Vogelsberger, M., Artale, M. C., et al. 2019, MNRAS, 486, 3766
  * McCall et al. (2024) McCall, H., Reiprich, T. H., Veronica, A., et al. 2024, arXiv e-prints, arXiv:2401.17296
  * Mehdipour et al. (2016) Mehdipour, M., Kaastra, J. S., & Kallman, T. 2016, A&A, 596, A65
  * Merloni et al. (2024) Merloni, A., Lamer, G., Liu, T., et al. 2024, A&A, 682, A34
  * Mernier et al. (2017) Mernier, F., de Plaa, J., Kaastra, J. S., et al. 2017, A&A, 603, A80
  * Miller et al. (2015) Miller, J. M., Kaastra, J. S., Miller, M. C., et al. 2015, Nature, 526, 542
  * Moustakas et al. (2013) Moustakas, J., Coil, A. L., Aird, J., et al. 2013, ApJ, 767, 50
  * Naiman et al. (2018) Naiman, J. P., Pillepich, A., Springel, V., et al. 2018, MNRAS, 477, 1206
  * Nasa High Energy Astrophysics Science Archive Research Center (2014) (Heasarc) Nasa High Energy Astrophysics Science Archive Research Center (Heasarc). 2014, HEAsoft: Unified Release of FTOOLS and XANADU, Astrophysics Source Code Library, record ascl:1408.004
  * Nelson et al. (2018) Nelson, D., Pillepich, A., Springel, V., et al. 2018, MNRAS, 475, 624
  * Nelson et al. (2019) Nelson, D., Springel, V., Pillepich, A., et al. 2019, Computational Astrophysics and Cosmology, 6, 2
  * Nicastro et al. (1999) Nicastro, F., Fiore, F., & Matt, G. 1999, ApJ, 517, 108
  * Nicastro et al. (2018) Nicastro, F., Kaastra, J., Krongold, Y., et al. 2018, Nature, 558, 406
  * Nicastro et al. (2005) Nicastro, F., Mathur, S., Elvis, M., et al. 2005, Nature, 433, 495
  * Pillepich et al. (2018) Pillepich, A., Nelson, D., Hernquist, L., et al. 2018, MNRAS, 475, 648
  * Pitrou et al. (2018) Pitrou, C., Coc, A., Uzan, J.-P., & Vangioni, E. 2018, Phys. Rep, 754, 1
  * Planck Collaboration et al. (2016) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A13
  * Planck Collaboration et al. (2020) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6
  * Predehl et al. (2021) Predehl, P., Andritschke, R., Arefiev, V., et al. 2021, A&A, 647, A1
  * Predehl et al. (2020) Predehl, P., Sunyaev, R. A., Becker, W., et al. 2020, Nature, 588, 227
  * Reid et al. (2016) Reid, B., Ho, S., Padmanabhan, N., et al. 2016, MNRAS, 455, 1553
  * Reiprich et al. (2021) Reiprich, T. H., Veronica, A., Pacaud, F., et al. 2021, A&A, 647, A2
  * Ricci et al. (2017) Ricci, C., Trakhtenbrot, B., Koss, M. J., et al. 2017, Nature, 549, 488
  * Richard et al. (2019) Richard, J., Kneib, J. P., Blake, C., et al. 2019, The Messenger, 175, 50
  * Rykoff et al. (2014) Rykoff, E. S., Rozo, E., Busha, M. T., et al. 2014, ApJ, 785, 104
  * Shull et al. (2012) Shull, J. M., Smith, B. D., & Danforth, C. W. 2012, ApJ, 759, 23
  * Sousbie (2011) Sousbie, T. 2011, MNRAS, 414, 350
  * Sousbie et al. (2011) Sousbie, T., Pichon, C., & Kawahara, H. 2011, MNRAS, 414, 384
  * Springel et al. (2018) Springel, V., Pakmor, R., Pillepich, A., et al. 2018, MNRAS, 475, 676
  * Strauss et al. (2002) Strauss, M. A., Weinberg, D. H., Lupton, R. H., et al. 2002, AJ, 124, 1810
  * Tanimura et al. (2020a) Tanimura, H., Aghanim, N., Bonjean, V., Malavasi, N., & Douspis, M. 2020a, A&A, 637, A41
  * Tanimura et al. (2022) Tanimura, H., Aghanim, N., Douspis, M., & Malavasi, N. 2022, A&A, 667, A161
  * Tanimura et al. (2020b) Tanimura, H., Aghanim, N., Kolodzig, A., Douspis, M., & Malavasi, N. 2020b, A&A, 643, L2
  * Tanimura et al. (2019) Tanimura, H., Hinshaw, G., McCarthy, I. G., et al. 2019, MNRAS, 483, 223
  * Tinker (2021) Tinker, J. L. 2021, ApJ, 923, 154
  * Tuominen et al. (2021) Tuominen, T., Nevalainen, J., Tempel, E., et al. 2021, A&A, 646, A156
  * Vijayan & Li (2022) Vijayan, A. & Li, M. 2022, MNRAS, 510, 568
  * Wang & Wei (2023) Wang, B. & Wei, J.-J. 2023, ApJ, 944, 50
  * Werner et al. (2008) Werner, N., Finoguenov, A., Kaastra, J. S., et al. 2008, A&A, 482, L29
  * Wilms et al. (2000) Wilms, J., Allen, A., & McCray, R. 2000, ApJ, 542, 914
  * Yang et al. (2022) Yang, K. B., Wu, Q., & Wang, F. Y. 2022, ApJ, 940, L29
  * Zhang et al. (2024) Zhang, Y., Comparat, J., Ponti, G., et al. 2024, arXiv e-prints, arXiv:2401.17309
  * Zonca et al. (2019) Zonca, A., Singer, L., Lenz, D., et al. 2019, Journal of Open Source Software, 4, 1298
  * Zou et al. (2019) Zou, H., Gao, J., Zhou, X., & Kong, X. 2019, ApJS, 242, 8

##  Appendix A Surface brightness profiles in different bands and length
populations

In Fig. 12, 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\sigma$ and $6.5\sigma$ 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.

We also split the filament sample into 20–40 Mpc and 40–100 Mpc length
populations, we detect a marginal $2.5\sigma$ discrepancy between the two
populations (see Fig. 13). 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).

Figure 12: Same as Fig. 3, but in the 0.3–0.7 keV (_left_) and 0.7–1.2 keV
(_right_) bands, respectively.  Figure 13: 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 cluster masking radius

We investigate the impact of the cluster masking radius on the stacked
profile. Instead of the $1.5r_{\lambda}$ masking radii for the optically
identified cluster with $\lambda>20$ as described in Sect. 2.1, we use
$3r_{\lambda}$ masking radii for a test.

Following Sects. 2.1 and 2.2, we create new 0.3–1.2 keV count and exposure
Healpix maps with optically identified clusters masked using $3r_{\lambda}$
masks. Other source masks are the same as those listed in Table 1. We plot the
smoothed rate map with $3r_{\lambda}$ masking radii in Fig. 14. Comparing to
the count rate map in Fig. 1, the doubled cluster mask size significantly
reduced the unmasked footprint.

In Fig. 15, we plot stacked net profiles using the two masking radii. The
profile with 3$r_{\lambda}$ masks generally follows the one with
$1.5r_{\lambda}$. The integrated net signal within 10 Mpc transverse distance
is 90% of the integrated signal from profiles with $1.5r_{\lambda}$.

Figure 14: Same as Fig. 1, but with $3r_{\lambda}$ mask radii for optically
identified clusters with $\lambda>20$. Figure 15: Comparison between stacked
net profile using $1.5r_{\lambda}$ (red) and $3r_{\lambda}$ (blue) cluster
mask radii.

##  Appendix C Mock maps for profile stacking validation

We first 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 create the mock filament emission
rate map based on the projected geometric configurations of the 7817 filaments
(top left panel of Fig. 16). We note that the best-fit parameters in this
section are empirical, and we do 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 simulate mock background rate maps based on the angular power spectrum of
the observed rate map (top right panel of Fig. 16). 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^{\circ}$
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 compute the angular power spectrum of
the relative amplitude map $C_{l,\mathrm{rel}}$ up to a maximum multipole
moment $l=3000$ using the function anafast in the HEALPY package. We simulate
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 sum up the mock filament rate map and the mock background rate map to
get a mock sky rate map. We obtain a mock sky count map by multiplying the
exposure map from observations and applying Poisson randomization to the
expected count value in each pixel.

|
---|---
|
Figure 16: 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 in as a comparison.

##  Appendix D Spectral models for undetected galactic sources

In Fig. 17, 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. 18, we plot the spectral models in different redshift bins and stellar
mass bins used in the mock catalog, as well as the median spectrum.

Figure 17: _Left_ : Multi-temperature differential emission measure
distribution of undetected X-ray halos at different redshifts after applying
the selection. _Middle_ : The corresponding composite spectra of each
differential emission measure distribution at different redshifts. _Right_ :
The synthesis rest-frame spectral model of undetected X-ray halos by taking
the filament sample length and redshift distributions into account.  Figure
18: 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.

  *[WHIM]: warm-hot intergalactic medium
  *[FRB]: fast radio burst
  *[CMB]: cosmic microwave background
  *[BBN]: Big Bang Nucleosynthesis
  *[SZ]: Sunyaev–Zeldovich
  *[RASS]: ROSAT All-Sky Survey
  *[]:
  *[eFEDS]:  Final Equatorial-Depth Survey
  *[eRASS]:   () All-Sky Survey
  *[TMs]: telescope module
  *[SDSS]: Sloan Digital Sky Survey
  *[PSF]: point spread function
  *[CXB]: Cosmic X-ray background
  *[SMF]: stellar mass function
  *[SMFs]: stellar mass function
  *[AGNs]: active galactic nucleus
  *[XRBs]: X-ray binary
  *[NYU-VAGC]: New York University Value Added Galaxy Catalogue
  *[CIE]: collisional ionization equilibrium
  *[DEM]: differential emission measure
  *[APEC]: Astrophysical Plasma Emission Code
  *[AGN]: active galactic nucleus
  *[XRB]: X-ray binary
  *[SFR]: star forming rate
  *[ARF]: ancillary response file
  *[RMF]: response matrix file
  *[ARFs]: ancillary response file
  *[1-T]: single temperature
  *[MCMC]: Markov Chain Monte Carlo
  *[UVB]: cosmic ultraviolet background
  *[eROSITA]: extended ROentgen Survey with an Imaging Telescope Array
  *[FWHM]: full width at half maximum