Issue 
A&A
Volume 675, July 2023



Article Number  A63  
Number of page(s)  14  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/202245187  
Published online  03 July 2023 
WHIMhunting through crosscorrelations between optical and SZ effect data in the Virgo cluster filaments
^{1}
Argelander Institut für Astronomie, Universität Bonn, Auf dem Hügel 71, 53121 Bonn, Germany
email: cagri.erciyes@dlr.de
^{2}
Department of Astronomy and Space Science, Chungnam National University, 99 Daehakro, Daejeon, 305764
Korea
Received:
12
October
2022
Accepted:
24
April
2023
Context. The physical state of most of the baryonic matter in the local universe is unknown, commonly referred to as the “missing baryon problem”. It has been theorized that at least half of these missing baryons are in a warmhot, lowdensity phase, outside of the virialized darkmatter halos.
Aims. We attempted to find the signature of this warmhot intergalactic medium (WHIM) phase in the filaments of the nearby Virgo cluster by using optical and SunyaevZeldovich effect data.
Methods. Specifically, we used a filamentgalaxy catalog created from the HyperLeda database and an allsky Comptony map extracted from the Planck satellite data for twodimensional crosscorrelation analysis by applying a spherical harmonics transform. The significance test is based on the nulltest simulations, which exploits advanced cutsky analysis tools for a proper map reconstruction. To place upper limits on the WHIM density in the Virgo filaments, realistic baryonic density modeling within the cosmic filaments was done based on stateoftheart hydrosimulations, within the signalboosting routine.
Results. The crosscorrelation signal is found to be too dim compared to the noise level in the Planckymap. At a 3 σ confidence level, the upper limit on volumeaverage WHIM density turns out to be ⟨ n_{e} ⟩< 4 × 10^{−4} cm^{−3}, which is indeed consistent with the WHIM parameter space, as predicted from simulations.
Key words: largescale structure of Universe / galaxies: clusters: intracluster medium / intergalactic medium
© The Authors 2023
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1. Introduction
The first indication that the universe might be holding 40%−50% of its baryons in its current epoch within the vast filamentary network of dark matter known as the “cosmic web” came from numerical simulations (Cen & Ostriker 2006). This phase of baryonic matter was named the warmhot intergalactic medium (WHIM) following its expected temperature and density ranges: T_{e} ∈ [10^{5} − 10^{7}] K and n_{e} ∈ [10^{−5} − 10^{−8}] cm^{−3} (e.g., Bertone et al. 2008; Shull et al. 2012; Martizzi et al. 2019) and to distinguish it from the denser and hotter plasma that is bound within the potential well of galaxy clusters. The expected properties of WHIM have made it obvious that it would be very difficult to detect it via its direct thermal emission in the soft Xray bands. Hence, the UV and Xray absorption lines were identified as optimal tools for the WHIM search, leading to several detections and upper limits (e.g., Danforth & Shull 2008; Fresco et al. 2020). Nevertheless, a full accounting of the local baryons via line absorption method is incomplete, which has led to the socalled “missing baryon problem”.
The onset of widearea millimeter and submillimeter surveys for detecting the SunyaevZeldovich (SZ) effect in galaxy clusters has provided a new approach to search for the WHIM phase. The SZ effect (Sunyaev & Zeldovich 1970) is a small spectral distortion on the cosmic microwave background (CMB) caused by the scattering off the highenergy thermal electrons embedded within the cosmic structures and, crucially, the linear dependence on density for the SZ signal makes it particularly attractive for the WHIM search (Mroczkowski et al. 2019, for a review). Moreover, SZ and Xray measurements can be considered a more “direct” probe for the WHIM than line absorption measurements, as the former two directly constrain the plasma properties without any assumption of its metal content.
Taking advantage of the fullsky SZ effect map provided by the Planck satellite (Planck Collaboration XXII 2016; Planck Collaboration XXVII 2016), a series of papers have used the stacking technique to claim detection of the WHIM in the cosmic filaments with roughly 5σ significance (Tanimura et al. 2019b; de Graaff et al. 2019). The stacking technique is a common approach in astronomy where the signaltonoise ratio of a faint signal is improved by averaging many (often several 1000 s) individual sky patches, where the location of the signal is inferred from some other data sets, thereby lowering the contribution from stochastic instrumental noise. In these studies, the location of filaments was inferred from bright, elliptical galaxy locations in the SDSS optical data (Alam et al. 2015). More recently, this stacking technique has been improved and newer filament catalogs based on SDSS data have been used, resulting in an improved detection from the same Planck SunyaevZeldovich (SZ) map (Tanimura et al. 2020a). Also, detections have been made in the Xray bands, with widearea Xray maps from the ROSAT and eROSITA satellites (Tanimura et al. 2020b, 2022). These WHIM detection results, based on the stacking technique, provide a statistical understanding of the WHIM properties, as the emissions and SZ imprints from many filaments are averaged into a mean signal.
Correlation measurement is one of the most popular techniques (besides the stacking approach) to reveal hidden information for the case when the targeted signal is below the system noise level. Correlation of fields measures the expectation value of the products of these fields at different spatial points; in other words, it provides the strength of the link between physical quantities in a specified scale. Crosscorrelation techniques can be used in a multipurposed manner; such as redshift estimation (Van Den Busch et al. 2020), revealing source characteristics over different observations (Dunkley et al. 2013; Cao et al. 2007) and exploring the nature of unresolved phenomena (Hashimoto et al. 2019). There are also studies that used the crosscorrelation technique for searching thermal SZ signature induced by hot gas (bound or unbound to clusters) by using CMB data (HernándezMonteagudo et al. 2004; GénovaSantos et al. 2015). Such correlationfocused CMBrelated works already provided some useful information regarding the physical properties of the intergalactic medium. Another example of crosscorrelation between optical and radio data sets is given in Amaral et al. (2021), who used this technique to put constraints on the upper bound for largescale magnetic fields in the intergalactic medium.
A direct search for WHIM from optically confirmed filaments, around an individual cluster, can also have complementary scientific value, for instance, to understand the specific properties of galaxies within a filament or to infer mass accretion rate from specific filaments onto clusters. The search of WHIM for individual targets has primarily focused on filaments joining nearby clusters, both in the Xray (Eckert et al. 2015; Reiprich et al. 2021) and in the SZ effect (GénovaSantos et al. 2005; Molnar & Birkinshaw 1998; Planck Collaboration Int. VIII 2013). But these systems cannot probe the nature of WHIM in filaments in an unbiased way, at best they show the transition from the WHIM phase to the denser intracluster medium (ICM) phase. Similarly, important studies for the transition between the unbound and bound gas around galaxy clusters have been made with detailed profile analysis in the cluster outskirts. Examples for this category of work, using the SZ signal, include an analysis of the Virgo cluster outskirts (Planck Collaboration Int. XL 2016) as well as stacking a sample of several hundred clusters (Anbajagane et al. 2022) from South Pole Telescope data.
In this paper, we attempt to measure the SZ signal contributed by WHIM directly within a known and wellstudied system: the filaments of the Virgo cluster. As our nearest massive galaxy clusters, the galaxy distribution along the various filaments connecting to the Virgo cluster has been very well studied (Tully 1982). Our tool is a direct crosscorrelation study in the mapsspace between the Planck SZ effect map and the optical galaxy distribution in the Virgo filaments, as has been described in a more recent study using a larger data set of HyperLeda database (Kim et al. 2016).
The current paper is structured as follows. In Sect. 2, we describe the two data products. In Sect. 3, the details of the crosscorrelation analysis. We simulated the WHIM signal to predict the expected nature of the crosscorrelation and to make null tests for error analysis, which are described in Sect. 4. In Sect. 5, the actual results based on Planck data are presented. A discussion and summary, along with an outline of possibilities for future improvement, are given in Sect. 6.
We assume a standard flat ΛCDM cosmology (Planck 2013 results in Planck Collaboration XVI 2014, Ω_{m} = 0.31, H_{0} = 67.8 km s^{−1} Mpc^{−1}) to convert the redshifts listed in the optical galaxy catalog to distances. We discuss in Sect. 6 whether the resulting choice of the Virgocentric distance could have any influence on our correlation analysis.
2. Data products
2.1. Optical catalog of filament galaxies
The filamentlike structures (i.e., prolate and oblate overdensities of galaxies) around the Virgo cluster, which is the nearest rich and dynamically young cluster of galaxies (Aguerri et al. 2005), were first studied by Tully (1982). The most recent and comprehensive attempt to map the largescale structures around the Virgo cluster exploiting the HyperLeda database (Paturel et al. 2003) had been carried out by Kim et al. (2016), whose filamentgalaxy catalog is used in this work. A visual representation of these various filaments is shown in Fig. 1.
Fig. 1. Galaxy distribution in the seven filaments and one sheet around the Virgo cluster taken from Kim et al. (2016). Different colors represent different structures. Bright (M_{B} < −19) and faint (M_{B} > −19) galaxies are denoted by large open circles and small filled circles, respectively. The large rectangular box is the region of the Virgo cluster enclosed by the EVCC (Kim et al. 2014). The gray dots represent galaxies not associated with the Virgo cluster or a particular structure. 
The catalog of galaxies contains 1100 members, grouped into seven filaments and one sheet, as depicted in Fig. 1. The search volume and radial velocity cutoffs are discussed in Kim et al. (2016), to which we refer for further details. One important criterion is the exclusion of the Virgo cluster members, which were defined by the Extended Virgo Cluster Catalog (EVCC; Kim et al. 2014, see the large rectangular box in Fig. 1) and excluded in the current study to eliminate possible signal contamination from the bound gas to the Virgo cluster itself (see Sect. 3.2 for the detailed description of the masking procedure).
Smoothed galaxy density map template obtained from this catalog is presented on the right side of Fig. 2. With the purpose of suppressing the sparsity effect on the optical data in which only ∼1000 galaxy positions exist within a widefield galaxydensity map is constructed by using the get_interp_weights() function in HEALPY module, which returns the four closest pixel indexes on the two rings above and below the actual index of the galaxy location, together with their corresponding weights (inversely proportional to their distances). In this way, the sparse signal is redistributed over the neighborhood by conserving the total power, which can be considered a preprocessing prior to smoothing. In this way, this procedure provides both the capability for bilinear interpolation along a latitudelongitude grid, as well as a more accurate representation of the galaxy points whose locations are much close to pixel boundaries. Thereafter, smoothing with a Gaussian beam whose FWHM is ≈30 arcmin was applied to the map to break the high discreteness level, for the purpose of allowing for the correlation statistics to be carried out adequately. Using different sizes of beams has a negligible effect on the main results, as mentioned in Sect. 5.3.
Fig. 2. Modeled noisefree ymap of the filamentary structures around the Virgo Cluster, which is also smoothed with the Planck beam (Sect. 4.1), shown on the left. Projection effects for each stated filament can be realized by comparing with Fig. 1. Constructed galaxy density map from the point locations (Sect. 2.1) by using get_interp_weights() function in HEALPY module and smoothing with a Gaussian beam whose FWHM is ≈30 arcmin, shown on the right. Both maps are centered on the Virgo Cluster. 
2.2. Planck allsky ymap
The very large extension of the Virgo cluster and its filaments in the sky means that only a spacebased millimeter and submillimeter telescope can recover the full extended signal of the thermal SunyaevZeldovich (tSZ) effect and, currently, the best data come from the allsky survey of the Planck satellite (Planck Collaboration I 2011). In addition to the sky coverage, the multiple frequencies of the Planck mission enables optimal separation of foregrounds, leading to a cleaner reconstruction of the tSZ signal. We used the allsky Comptonyparameter maps (Planck Collaboration XXII 2016), where the yparameter is the frequencyindependent amplitude of the tSZ signal and is defined as the lineofsight integral of the electron pressure (Sunyaev & Zeldovich 1970):
where σ_{T} is the Thomson crosssection, n_{e} is the electron number density, T_{e} is the electron temperature, k_{B} is the Boltzmann constant, and m_{e}c^{2} is the electron rest mass energy. A visual representation of the ymap, on the half skysphere centered on the Virgo cluster, is shown in the left image of Fig. 3.
Fig. 3. NILC ymap of the halfsphere which is centered on the Virgo Cluster (left). Masked version of the left map, described in Sect. 3.2 (right). The color scale denotes the value of the Comptonization parameter, as inferred from multifrequency Planck data. 
The extraction of the ysignal from multifrequency data of Planck is done via separate enhanced modifications of the Internal Linear Combination (ILC) method, which assumes individual maps as a linear superposition of multiple sky components which are spatially uncorrelated (see Planck Collaboration XII 2013, for an overview). The two main methods are known as MILCA (Hurier et al. 2013) and NILC (Delabrouille et al. 2009). We specifically used the publicly available Planck NILC maps for our analysis, since their lowest multipole window function can be expected to provide a better response to the largestscale emission (see Fig. 1 in Planck Collaboration XXII 2016). However, we also compared our results with MILCA maps and obtained similar upper limits, which is likely not surprising since both these maps have similar noise properties, particularly in high Galactic latitudes and at around ℓ ∼ 40, where we searched for our key crosscorrelation feature.
3. Crosscorrelation analysis
3.1. Spherical harmonics
In the case of a twopoint correlation function, structures can be quantified by the correlator using the expectation value integration:
where the average is taken over all positions y and all orientations of the separation vector, x, on a field function, f. For that definition, the twopoint correlation function would be only dependent on the scale length x by taking advantage of the spherical symmetry feature of the field as a manifestation of homogeneity and isotropy characteristics on large scales. Besides, if the field density preserves the Gaussianity, the whole statistical information on the corresponding probability density function is contained in ξ(x).
Given that the calculation of the twopoint correlation function on real space is computationally burdensome (especially for thousands of maps, see Sect. 4.2), we calculate it through the Fourier transform of the corresponding power spectrum. However, standard Fourier modes are not appropriate for describing functions that are defined over spherical surfaces or over a large sky region, such as the entire Virgo filaments system, which covers roughly onefourth of the celestial sphere. The solution is to use the angular power spectrum, defined via spherical harmonics functions that provide a complete set of orthonormal basis for Fourier transforms on a sphere. That is the reason we used the PYTHON module HEALPY^{1}, based on HEALPix (Gorski et al. 2005), to properly manipulate pixelated data on the sphere and apply spherical harmonics transforms. In addition to this software, we also made use of other cut sky analysis tools such as PyMaster^{2} and PolSpice^{3} which have additional features for generating full sky power spectrum from the pseudo power spectrum. In particular, we used PyMaster (Hivon et al. 2002) to implement statistical null tests, as described in Sect. 4.2.
We compute the angular crosscorrelation function (CCF) from the crosscorrelation power spectrum using the standard definition as follows (see Appendix A for details):
For our case, where we analyze the crosscorrelation signal between two maps, the Fourier amplitudes in Eq. (A.2) can be calculated for both of these maps and used for the estimator of angular crosspower spectrum
where y_{lm} and g_{lm} are the Fourier amplitudes which are calculated from the NILC ymap and galaxydensity map, respectively. This estimator does not handle the impact of pixel masking properly (masked pixels are directly considered zerovalued). If conversion from the power spectrum of the pixelated map to the hypothetical unpixelated one is needed, then we should use effective pixel window function (see Appendix B and Eq. (15) in Gorski et al. 2005; Hivon et al. 2002, respectively).
The effect of masking should be taken into account and this will lower the signal power (e.g., see Fig. D.2 in Planck Collaboration XI 2016). The correction calculation roughly can be done by introducing effective sky fraction, f_{sky}, which is included as a multiplicative factor in Eq. (4) as (Hill & Spergel 2014). However, f_{sky} correction is only true at the first order, and applying a cut sky couples the harmonic modes between them. That is the reason why a modemode coupling kernel is introduced and corrected using PyMaster to obtain an unbiased estimate of the fullsky power spectrum (Eqs. (14) and (15) in Hivon et al. 2002) of a given masked field. The details of the masking procedure are outlined in Sect. 3.2. Ignoring this fact might induce bias in the derived parameters from likelihoods on the estimated .
For checking the statistical significance of the crosscorrelation signal coming from the specified sky patch around the Virgo cluster by implementing statistical nulltests (Sect. 4.2), we only consider a comparison of the observed signal with random correlations arising from the same configuration space: same pixelization, masking, and twopoint statistics. We note that smoothingeffect (instrumental beam and finite pixel size) corrections do not need to be applied here. We then exploited the PyMaster package (Sect. 4.2) for generating random realizations.
3.2. Masking procedure
In our work, we aim to probe lowdensity regions of the cosmic web, namely, diffuse WHIM gas located within the filaments around the Virgo cluster. We thus mask other possible sources of correlation as much as possible. First, to localize our study, we remove the region outside of the filaments’ sky patch, which corresponds to RA ≈[120, 230] and Dec ≈[−30, 55] range of degrees. Then, besides the sky patch cropping, a galactic plane mask that retains 60% of the sky with lowest emission at 353 GHz after CMB subtraction (foreground galactic emission masks)^{4} is applied with a “no apodization” option to avoid residual galactic foreground contribution induced mainly by diffuse thermal dust emission.
As we are only interested in unbound diffuse gas within the filamentary structures located outside of virialized halos, the galaxy clusters detected in the Planck tSZ map (PSZ2; Planck Collaboration XXVII 2016), MCXC catalog (Piffaretti et al. 2011), and SDSS DR9 galaxy clusters optical catalog (Banerjee et al. 2018) are masked out. A galaxy cluster mask was applied with the same strategy as in Tanimura et al. (2019a), by setting the masking radius to three times the cluster radius (3 × R_{500}). Specifically for the Virgo cluster, to eliminate an impact of the dynamical association of the filaments with their parent cluster (e.g., possible reverse shocks that are generated by infalling matter along the Virgo filaments), we can apply a more conservative mask whose radius is three times the angular size of R_{200} of the Virgo cluster (mainly to exclude EVCC region in Fig. 1). For consistency, we adopted the same virial radius R_{vir} as in Urban et al. (2011), corresponding to ∼1.08 Mpc, with a projected radius of ∼3° .9 in the sky. For Planck SZ clusters without assigned radius, a mask with 10 arcmin of radius corresponding to the effective angular resolution (FWHM) of the NILC ymap was exploited.
In addition to the contamination from clusters, there might exist a residual emission due to the point sources. In the reconstructed tSZ maps, radio sources appear as negative peaks while infrared sources show up as positive peaks mimicking the signal from clusters. To avoid contamination from these sources, the union of individual frequency point source masks (HFI and LFI point source foreground masks)^{4} was used as suggested in Planck Collaboration XXII (2016).
The union of each abovementioned mask is applied onto both Compton ymap and galaxy density map before crosscorrelation computation. Mask apodization is redundant for the purpose of spherical harmonics expansion on the cut sky in which case there is no gradient calculation on the field, so an additional apodization was not applied on our mask (in other words, it is a type of tophat mask). Masked pixels are ignored by being interpreted as zeros in HEALPY spherical harmonics transform operations. After all of this filtering, the total unmasked fraction of the sky is f_{sky} ≈ 10%, which corresponds to the half of the proportional area that is enclosed by the filaments around the Virgo cluster (the region with RA ≈[120°, 230°] and Dec ≈[−30°, 55°] covers an area of f_{sky} ≈ 20%, and after carrying out the masking operation, half of this area is filtered out, so then the final reduced region has f_{sky} ≈ 10%). In Fig. 3, the left panel shows the half sphere of NILC ymap which is centered on the Virgo cluster and the right panel shows a masked version of the same map whose mask is graycolored.
4. Simulations
Before analyzing the real data to check for the existence of a crosscorrelation signal between the NILC ymap and the galaxydensity map, it is worth exploring what kind of signal would be expected if there were no noise and the sky signals (in both Comptony and optical) only consisted of the Virgo filaments. In this clear representation of the crosscorrelation signal, there will be no contamination from any other additional source and no random correlation between structures, and it can be used for comparison with actual or boosted (Sect. 5.2) signal and interpreting the effect of noise.
The galaxy density map is “clear” by construction, so pixels outside the smoothed points are set to zero, which is equivalent to masking these regions. However, this is not the case for the ymap, which contains thermal SZ contribution from all the structures in the universe as well as residuals from other contaminating signals, plus noise. Thus, to get a clear representation of the tSZ signal coming from diffuse gas within the Virgo filaments, we constructed a noisefree filament ymap in which the filament topology is directly reflected by the galaxy density map, based on the assumption that galaxies within each filament are distributed around the main spine axis of these structures. A uniform crosssection curved cylinder around each axis is then filled with WHIM for our study. The details of this modeling are explained below.
4.1. Filament Compton ymap reconstruction
4.1.1. Sampling the filament gas particles around spineaxes
For simulating the pressure profiles for the diffuse WHIM within known filaments, the first operation would be to sample the gas particles around the main spine axis of these filaments. The toy model is used in a way that a single filament volume shape is considered as a regular curved cylinder. The algorithm consists of two parts: a 3D fit to extract the main spine axis of the filament and uniform sampling around the main spine axis. Each of these steps is outlined below.
For a 3D fit to positions of galaxies at each filament, a leastsquares method with the option of Huber Loss (Huber 1992) is used, to minimize outliers’ influence in the fitting process. This robust loss function is widely used in the context of machine learning applications, among others. Then, a cubic spline interpolation is applied onto fitted points and a smooth 3D fit for the main spine axis is obtained. For each filament, the spine axis consists of N ordered points, namely, x_{1}, x_{2}, ...x_{N}. To be able to perform suitable rotation and scaling operations, a Cartesian coordinate transformation is applied onto the axis points. Sampling is iteratively done for each point along a spine axis and these points act as the center of constructed discs. In each iteration, particles are uniformly sampled around these centers in the xy plane, and then uniform random perturbation along the zaxis is added onto the sampled points. In this way, a uniform disc whose thickness is equal to x_{n + 1} − x_{n} and whose normal lies along the zaxis is constructed. Then, the constructed disc is rotated by using a 3D rotation matrix and linearly transformed to the exact location x_{n} on the spine axis. At each point interval on the spine axis, 30 000 particles are used for sampling. In the end, we have roughly 40 million particles for all filaments in total.
There is one sheetlike structure besides the filamentary structures in our sample, which is called the WM sheet. It includes galaxies from the W and M clouds of the Virgo cluster (Kim et al. 2016). It has been interpreted as a filament as well, whose spine axis is chosen to be nearly perpendicular to the line of sight to enhance any potential correlation signal close to the EVCC region (Fig. 1). The effect of this approximation for the WM sheetlike structure would be also negligible in our crosscorrelation analysis because most of this structure lies within a radius 3 * R_{200} around the central position of the Virgo cluster and then is mostly masked out.
In the sampling process, the radius of discs, in other words, the filament radius is determined by using scaling relations from Gheller et al. (2015, 2016). Since the filament lengths are already known, we can use the scaling relation between enclosed gas mass and estimated length of filaments, which is shown in Fig. 8 in Gheller et al. (2015) and Fig. 3 in Gheller & Vazza (2019). As can be clearly seen from those figures while the filament length increases, the deviation for the gas mass within the scaling relation reduces considerably. For this reason, the longest filament (NGC 5353/4) whose length is ≈22 Mpc (Kim et al. 2016) is chosen to estimate gas mass range which then can be used to determine reasonable filament radii. For ≈22 Mpc, the relevant mass range is roughly 10^{12.7} − 10^{13.7} M_{⊙}. Then, the relation between the enclosed gas mass and radius of the filaments (Fig. 9 in Gheller et al. 2015) is used to get an estimated radius range for NGC 5353/4, which is roughly 0.4 − 1.2 Mpc. The median value within this range, which equals 0.8 Mpc, is assigned as the filament radius for NGC 5353/4. We note, however, that an accurate assignment of the filament radius is not very important for our study, as we are interested in only the approximate shape and normalization of the crosscorrelation function.
Then, the remaining issue is estimating radius ratios for the filamentary structures and we can directly use the scaling relation between filament gas mass and volume once again. However, this approach may enhance the deviations in relations because of the filaments with short lengths in our sample. Hence, we used another way to get ratios of filament radii following Gheller et al. (2016). In Fig. 18 (Gheller et al. 2016), a tight correlation between the host filament total gas mass and total mass of galaxies, and also between the length of the host filament and total mass of galaxies, can be clearly seen. We used these to recheck the mass range we get for NGC 5353/4. According to the length of 22 Mpc, the related galaxy mass range is roughly 10^{12} − 10^{13} M_{⊙} and this corresponds to filament gas mass range of 10^{12.7} − 10^{13.7} M_{⊙}. It is shown that the mass of filaments tends to grow almost linearly with that of the resident galaxy halos within filaments.
By assuming nearly similar galaxy population properties for all filamentary structures on average (all filaments mostly consist of faint galaxies with M_{B} > −19, which corresponds to ≈88% of the total filamenthosted galaxy sample, and that bright galaxies exhibit uneven distribution with the sheer amount of gaps along their host filament, as underlined in Kim et al. 2016), we can expect a linear correlation between filament gas mass and the number of galaxies hosted by the filament. By taking advantage of this approach, and considering filament shapes as cylindrical, the filament gas mass is related to LR^{2}, where L is backbone length and R is filament radius. This way, we obtain the radius ratios of these filaments whose spine lengths are already known. The rest of the filament radii are determined according to these ratios and also picked reference radius value for NGC 5353/4. In Fig. 4, uniformly sampled gas particles around the fitted spineaxis for LeoIIA filament can be seen.
Fig. 4. LeoIIA filament in threeaxes: right ascension (RA, in degrees), declination (Dec, in degrees), and lineofsight distance (in megaparsecs). Filamenthosted galaxies and uniformly sampled gas particles (to simulate WHIM diffuse gas) are shown with dark and light blue points, respectively. The 3D fitted spine axis is shown with a green curve around which the gas particles are sampled. Distances of filament galaxies are derived from their heliocentric radial velocities (retrieved from HyperLEDA in Kim et al. 2016) based on the assumption of a linear Hubble distance – redshift relationship by applying a standard flat ΛCDM cosmology (Planck 2013 results in Planck Collaboration XVI 2014, Ω_{m} = 0.307, H_{0} = 67.8 km s^{−1} Mpc^{−1}). 
4.1.2. Temperature selection for isothermal filaments
After obtaining uniformly sampled particles, the lineofsight integral of the pressure needs to be calculated to extract Comptony value within each pixel on our HEALPix sky maps. First, an isothermal filament model is assumed, which is also motivated by Fig. 8 in Gheller & Vazza (2019) and Fig. 17 in Gheller et al. (2015) that show transverse profile of gas temperature as a function of the scale height from the spine axis for the representative filaments. It seems that gas temperature remains nearly uniform from the spine to the outer parts for most of the models. Gheller et al. (2015) explains this phenomenon arguing that inner parts of the filaments are thermalized at an almost constant temperature by shock waves propagating from the center outwardly, but the outer part of the filaments can comprise nonshocked cells that lower the average temperature. In addition, it has been observed that filaments follow welldefined scaling relations in the gas temperaturemass plane (Gheller et al. 2015). Within our sample of dynamically active Virgo filaments, the expected mass range is 10^{12.5} − 10^{14} M_{⊙} and the corresponding average massweighted gas temperature range is roughly 10^{6} − 10^{6.8} K, which is shown in Fig. 3 from Gheller & Vazza (2019) and Fig. 10 from Gheller et al. (2015), with a large scatter even for highmass filaments. Because of this broad dispersion, we did not apply scaling relation and used the logarithmic median temperature value within the specified range: 10^{6.4} K is picked as a diffuse gas temperature for each isothermal filament in our sample. The uncertainty in the Compton yparameter calculated over these model parameters is addressed in Sect. 5.3. In the context of normalized crosscorrelation function shape, the temperature selection only affects the scaling – and not the signal shape itself.
4.1.3. Density modeling and projection
The Comptony value in an isothermal gas is simply proportional to the lineofsight integral of the gas density:
Numerically, can be easily obtained from uniform point sampling under a given density function n, and the lineofsight distance can be calculated from the interval between the nearest and furthest particles within each pixel on the map. Regarding the density profile, the betamodel is chosen as the average transverse profile of gas density for selected representative filaments, as shown from simulated models in Fig. 8 from Gheller & Vazza (2019) and Fig. 17 from Gheller et al. (2015). It is clear that for simulations with a radiative scenario, the central mass density becomes roughly 2 to 3 times higher than average, and the profile gets steeper due to the higher compression levels (Gheller et al. 2015). Commonly, all density profiles exhibit a smooth drop by a factor of broadly 2–10 at ∼0.3r_{fil} from the spine axis, then slowly declining outwardly. Uncertainty on the density profile steepness and its effects on the results are discussed in Sect. 5.3. With the purpose of CCF estimation and initial null tests, an almost shallow density profile is picked and sampled with a quasiMonte Carlo integration to obtain an expected density value E(n), which corresponds to . The used isothermal βmodel the diffuse WHIM gas density profile is written as
1: y_map ← initialize() ⊳ null ymap
2: pixel indexes ← ang2pix(spatial positions)
3: unique indexes ← unique(pixel indexes)
4: for each unique indexes do
5: for each structures do
6: windowed particles ← fetch(unique index)
7: n̄ ← mean_density(windowed particles)
8: l ← max_separation(windowed particles)
9: ŷ ← calculate_Compton_y(n̄, l, T) ⊳ Eq. (1)
10: y_map[unique index] += ŷ
11: end for
12: end for
13: y_map smooth(y_map) ⊳ with Planck beam
where r_{c} is the core radius of the distribution and n_{e}(0) is the central electron density value. To be able to construct relatively shallow density profiles with smooth drop at the inner region and nearly flat outwardly for the diffuse gas, we kept r_{c} and β parameters small, namely, r_{c} = 0.06 r_{filament} where r_{filament} is the fixed filament radius, and β = 0.17 retains the relation of ⟨n_{e}⟩≃0.42 n_{e}(0), where ⟨n_{e}⟩ is the transversal (widthwise) electron mean density which corresponds to in Eq. (5).
The final stage to construct a noisefree ymap from the uniformly sampled isothermal gas particles is the projection. The pseudocode for the algorithm is shown in Algorithm 1. First, for each of the sampled particles, host pixel indexes are found by using ang2pix() function in HEALPY module and a unique list of pixel indexes is obtained. Then, for each unique pixel index, particles located inside the boundaries of the pixel with a given index are fetched and mean density (by simply averaging the local density values which are attached to particle positions w.r.t. vertical distance from spine axis) and the lineofsight distance interval is calculated for the obtained set of particles. As a last step, the Compton yparameter is calculated by multiplying these quantities with fixed temperature value and the rest of the constants as given by Eq. (1). This calculation needs to be done separately for each filamentary structure because of the linearly additive behavior of the Compton yparameter in the case of overlap, for instance, the node region of LeoII filaments in Fig. 2. After finishing the main loop to get estimated Compton yparameter values over pixels, the generated map is smoothed with the Planck beam whose FWHM is ∼10 arcmin.
Projection effects can be obviously seen for filamentary structures on the left side of Fig. 2. Leo Minor seems to be the largest filament because of its closeness; indeed, it is one of the smallest filaments in terms of the number of hosted galaxies and also its backbone length. LeoII A and LeoII B are overlapping on the ymap, observed as a multistem structure. NGC 5353/4 is the longest and most distant filament in our sample and extends out from the NGC 5353/4 group, running tangentially past the Virgo cluster rather than pointing toward it which is seen to be very thin (Kim et al. 2016). Additionally, Canes Venatici is one of the filaments running from the vicinity of the NGC 5353/4 group to the Virgo cluster and its maximal elongation component lies along a SGY direction, which is the axis of the supergalactic coordinate system that roughly coincides with the line of sight from our Galaxy.
4.1.4. CCF calculation
Before retrieving the expected correlation function between these two maps shown in Fig. 2, the cumulative mask, extracted in Sect. 3.2, was applied. To calculate the cross(or auto)correlation function, anafast() function is used to obtain cross(or auto)power spectrum and then the bl2beam() function exploits this transfer function in spherical harmonic space to compute a circular beam profile in real space – in other words, the correlation function. In Fig. 5, the normalized crosscorrelation function (CCF) between modeled filament ymap and galaxy density map is shown (red dashed line). There is an obvious bump at large scales around 35° and to test its validity, we also calculated the autocorrelation function of the filament ymap, which can also be seen in Fig. 5 (blue line). The exact coincidence of the two bumps confirms its reality, as the characteristic angular scale of the Virgo filaments, which we explain in Sect. 5.2.
Fig. 5. Expected normalized correlation functions for the masked version (Sect. 3.2) of the noisefree (and also without any contamination) fields that are shown in Fig. 2. The blue line points up the autocorrelation function for the modeled filament ymap (left panel in Fig. 2). Red dashed line shows the crosscorrelation function between filament ymap and galaxy density map (right panel in Fig. 2). The characteristic statistical feature of these filamentary structures in terms of correlation degree is obviously seen as a bump in large scales roughly between 30 and 40 degrees, which is shown with a gray band. This scale corresponds to the average separation between primary elongation occurrences of the filamentary structures, as the bump in correlation function indicates that there should be a specific repetitive interval within the analyzed topology. 
In addition, we performed a null test for this excess correlation signal at a 35° scale. We constructed random maps only consisting of different numbers of straight lines, whose lengths are within a similar range of Virgo filament spine lengths. While drawing these lines on the latitudelongitude grid, position pairs for the endpoints were picked randomly by preserving a filament length scale and pixels residing along this stripe between endpoints were set to one, as the rest of the field is zero initially. The CCF between these random maps and the galaxy density map was calculated, and no bump was observed in these null tests, which confirms our hypothesis as the bump is a characteristic topological feature that we targeted for our CCF signal identification (Sect. 5.2).
4.2. Nulltest for random correlations
With the intent of assessing the significance of the measured CCF signal between NILC ymap (right panel in Fig. 3) and smoothed galaxy density map (right panel in Fig. 2), the range of CCF values need to be calculated when there is no real correlation expected between these two maps. For these null tests, we performed Monte Carlo based simulations by creating random realizations of the actual NILC ymap by using both PyMaster (NaMaster python wrapper) and HEALPY. The pseudocode is shown in Algorithm 2. First, the hypothetical fullsky map is defined and characterized by using NmtField() function (creates an object that contains all the information describing the fields to correlate including their observed maps, masks, and contaminant templates) considering NILC ymap as the main input field as an observed map and specifying the mask that is created in Sect. 3.2 (used in the rest of analysis). Then, the fullsky autopower spectrum of the given masked ymap is estimated by using compute_full_master() function (estimates of the fullsky power spectrum of two masked fields by calculating deprojection bias power spectra and applying decoupling routines). This reference power spectrum C_{l} is used to create 1000 random fullsky maps with synfast() function (applies conversion from power spectrum to map via spherical harmonics) and the abovementioned mask is applied onto these simulated maps afterward. As the last step, CCFs are calculated between these randomly created masked maps and a given smoothed galaxy density map.
1: random_maps ← initialize() ⊳ empty list is assigned
2: NULL_CCFs ← initialize() ⊳ empty list is assigned
3: y_ field ← construct_ field(NILC_y_map, mask)
4: full_sky_C_{l} ← compute_ full_master(y_field, y_field)
5: simulation_number ← 1000
6: for i in range(simulation_number) do
7: random_map ← synfast(full_sky_C_{l})
8: masked_random_map ← mask(random_map, mask)
9: random_x[i] ← masked_random_map
10: NULL_CCFs[i] ← CCF(random_maps[i],
11: galaxy_map)
12: end for
We followed this strategy for the null tests to create random masked maps whose statistics are as much as close to the masked NILC ymap in terms of the power spectrum and also the probability density function (PDF) of unmasked pixels. In Fig. 6, on the left panel, power spectra of the randomly created masked ymaps and the one from masked NILC ymap are compared. On the right, 1D PDF of unmasked pixels of the NILC ymap and one of the randomly created maps via NaMaster are presented. In both of the cases, statistics match well enough and it clearly indicates that our approach is wellsuited for the null test within the context of the CCF significance measure.
Fig. 6. Power spectrum comparison between randomly created maps via Algorithm 2 (blue lines) and actual NILC ymap (red line), including the cumulative mask (Sect. 3.2), shown on the left. 1D probability density functions (PDF) of the unmasked pixel values for the NILC ymap (blue bars), one of the simulated maps (with NaMaster) by using the analyzed NILC ymap field itself (green bars), and one of the simulated maps by using the halfdifference map (red bars), shown on the right. Gaussian fits are sketched for the histograms of randomly simulated maps (black dashed lines). Excess of the signal seen in the NILC ymap (blue bars at the leftmost and rightmost parts of the histogram) compared to the random simulation (green bars) in the right panel indicates some measure of nonGaussianity originated from ∼1/10 000 of the total map pixel number, and therefore limits the statistical representability of the underlying field (NILC ymap, at this case) via random realizations obtained from a crosscorrelation method (Sect. 3.1). Together with the left panel, it’s obvious that usage of the halfdifference map is not suitable for the null tests as it’s not providing the matching statistics. On the other hand, by taking advantage of NaMaster to create simulated maps which are obtained by using the own statistics of the analyzed field (right panel in Fig. 3) appears to be the best option with respect to conducting null tests in our case. This is because it mimics the reference field also statistically by reaching a good agreement with the signal characteristics, which can be directly seen from both 1D PDF and power spectrum comparisons. 
For the purpose of the null test, we also checked the usage of halfdifference maps:
where y_{1} and y_{2} are the ymaps generated out of the NILCpipeline corresponding to the two jackknife halves of the data; hd is the halfdifference of these two ymaps, which can be used in noise characterization as it gives an estimate of map noise (N_{1} and N_{2}, contains all noise components including 1/f residuals, white noise, and other possible systematics), where the astrophysical signal (S) is canceled out in the measured data. Planckymaps exhibit a highly inhomogeneous noise, so there is an expected strong spatial dependency as well.
In Fig. 6, a comparison of the histograms of pixel values between masked NILC ymap and halfdifference map is shown on the right panel, with corresponding Gaussian fits for each of them. As can be seen, 1D PDF of the halfdifference map is well described by the Gaussian fit, with negligible tails in both directions as expected. For the masked NILC ymap, the Gaussian fit is also wellsuited except for a considerable proportion of the endmost pixel values, which indicates that the possible contamination could come from a remaining astrophysical signal from the unmasked structures. The most important outcome from this comparison is that the halfdifference noise map is not really a good choice to be used in the null test analysis in our case because the statistic is significantly unmatched in terms of a 1D PDF. The masked NILC ymap has a higher variance in comparison with the halfdifference map due to either NILCpipeline or residual foreground contamination, and this residual must be included in the CCF robustness estimation. Hence, the measured power spectrum of the masked ymap can be directly used to create statistically identical random realizations providing a characterization of both statistical and systematic errors in a more robust way.
To perform the tests, the Algorithm 2 was applied with the primary input of the halfdifference map as well as the NILC ymap to compare a range of the CCFs arising from random correlations. We generated 1000 random realizations of the masked halfdifference map and NILC ymap via the abovementioned strategy (by using functions provided by PyMaster and HEALPY packages), and corresponding CCFs with the galaxy density map were calculated. In Fig. 7, the cumulative range of nullCCFs is presented for the usage of both NILC ymap and halfdifference map by comparing deviation levels up to 3σ. On all scales, scattering is much broader for the case with direct NILC ymap usage, which can be already expected by just looking at 1D PDF results directly in Fig. 6, as the halfdifference map manifests a smaller amplitude range and also lower spectral power magnitude.
Fig. 7. Measured crosscorrelation function between masked NILC ymap and constructed galaxy density map is shown with a black dashed line, as a potential indicator of the existence of WHIM in the filamentary structures around the Virgo Cluster. Green and red colored regions correspond to the 3σ range around the mean value (the blue line) being extracted from null tests by using masked NILC ymap and halfdifference map, respectively. In other words, these shaded regions correspond to our benchmark to check the significance of the measured signal. It is obvious that the measured signal cannot reject the null hypothesis (onetailed hypothesis test for the expected positive correlation), even when the statistic of the halfdifference map is exploited, as the signal does not exceed the positive side of the 3σ range in the benchmark. 
5. Results
5.1. Nondetection with Planck data
Since there is an expected positive correlation between the modeled filament ymap and constructed galaxy density map as seen in Fig. 5, the significance of the measured CCF between NILC ymap and galaxy density map is checked via a onetailed null hypothesis test. In Fig. 7, the measured signal appears to be within nondetection case, in other words, the null hypothesis cannot be rejected. Even when only the instrumental noise statistic (from halfring maps) is used for the null tests, the signal itself is still too weak, which means that noise is too dominant over the related Compton y signal coming from the diffuse gas within these filaments, if any. As is the common standard, the statistical significance threshold is set to 3σ, and the corresponding range is shown in Fig. 7. The measured CCF is far too small from exceeding this level throughout the entire range of scales (from 0 to 50 degrees). Moreover, it does not exhibit the expected CCF characteristic as a bump around the scales ∼30–45 degrees specified for analyzing the WHIM signal (Fig. 5). In addition, the mean CCF value within the null tests is found to be nearly zero (blue line in Fig. 7), which shows the consistency of the applied method because the expected correlation between random ymap and galaxy density map is zero.
The considerable variation when using the masked ymap rather than a halfdifference map on the correlation amplitudes can be already seen in Fig. 7. Since it preserves both systematic noise sources and the instrumental noise, for further analysis, we go on to consider that the nullCCF range getting extracted by using the masked NILC ymap via Algorithm 2 (a greencolored region in Fig. 7) as a reference threshold for setting the upper limits on the WHIM physical parameters.
5.2. Density upper limits with Planck data
We used the nondetection result to constrain the density upper limits of the hypothesized WHIM phase within the filamentary structures around the Virgo cluster. Based on the nulltest results (Fig. 7), a positive 3σ threshold is picked for the “detection” case, namely, to reject the null hypothesis. Starting with the reference model filament ymap (left panel in Fig. 2), with a relatively low central density value (n_{e}(0)) of the diffuse gas, the strategy is to boost the signal by increasing the axial n_{e}(0) iteratively, until the calculated CCF between boosted ymap and galaxy density map exceeds the 3σ limit which corresponds to the upper boundary of the greencolored region in Fig. 7 (see Sect. 5.1 regarding why only the greencolored region is taken into account as a reference and the redcolored region is thus ignored).
At this point, two options exist regarding the ymap that will be used to calculate the enhanced signal: either the boosted model ymap of the filaments directly or the composite ymap done by adding the boosted model ymap onto the actual measured NILC ymap. The latter approach is used in Brown et al. (2017) for limiting magnetic field values in the cosmic web with diffuse radio emission. Besides, the termination criterion from an iteration is also crucial, namely, the question of which interval should be used as a check when the 3σ limit is exceeded. In Brown et al. (2017), there is no constraint regarding this criteria, the iteration is terminated when the boosted CCF signal exceeds 3σ at any degree. In Tanimura et al. (2019b), the likelihood of the data exceeding a null hypothesis was calculated by comparing the χ^{2} of the measured data to the χ^{2} of the probability distribution function of 1000 null samples.
In our case, the boosted modeled filament ymap can directly be used as a reference. In other words, there is no need to contaminate a reference ymap with the noisedominant ymap (i.e., the masked NILC ymap, Fig. 7). Nevertheless, the derivation of the density threshold numbers is done using both techniques, and we note some important points regarding both. First, the boosted signal from the composite ymap can still have contamination of the correlated noise (instrumental or other systematic), as can be seen in Fig. 8 (turquoise dashdotted line). This can complicate the interpretation of the density limit values. Secondly, a specific characteristic of the expected CCF (the bump between roughly 30 and 45 degrees, Fig. 5) is exploited to eliminate the noise contamination effect, especially when using the composite map (Fig. 8).
Fig. 8. Moment of the null hypothesis rejection, in other words, the “detection” case (Sect. 5.2), as the CCF signal of the boosted modeled filament ymap (red dashed line) exceeds the 3σ level at the positive side of the null hypothesis range (green shaded region) for a single point within the characteristic scale (gray band in Fig. 5), where the actual measured CCF signal (dotted black line) is far below this level. In fact, the measured signal appears to be totally lost as it even could not exceed the 1σ range (1,2, and 3σ intervals for the null hypothesis are shown with different colored bars: orange, blue, and green, respectively). The turquoise dashdotted line corresponds to the CCF signal being extracted from the composite map as a summation of the masked NILC ymap and boosted modeled filament ymap. 
For each constructed bin with a width of 0.2 degrees throughout the angular range (0 to 50 degrees), null samples have a Gaussian distribution centered on zero, so a pvalue calculation according to the chisquare test for nulltest rejection is equal to the direct 3σ thresholding. As a consequence, the optimized boosting algorithm is as follows: the modeled filament ymap is initialized with relatively low central density and then it is iteratively increased until the calculated CCF exceeds the 3σlevel of nullrange in any point between 30° and 45°. The CCF was calculated between the galaxy density map and boosted modeled ymap directly, not combined with NILC ymap at this point; however, adequate results were also obtained with the composite ymap usage for a robustness check. Within the iterative loop, to obtain the density value corresponding to the limit for null hypothesis rejection, an adapted version of the NewtonRaphson method (firstorder in the class of Householder’s methods which are numerical algorithms for rootfinding) is used as an approximated interpolation of the electron density between the two steps of the iterative process with the CCF lying below and above the 3σ level, with a tolerance of Δn_{e}(0) ≤ 3 m^{−3}.
In Fig. 8 (reddashed line), we can see the moment where the CCF between boosted modeled filament ymap and galaxy density map exceeds 3σ threshold at the positive side (detection case) with a tolerance of Δn_{e}(0) ≤ 3 m^{−3}, where the obtained central and mean (transverse) electron densities are n_{e}(0)≈3 × 10^{−4} cm^{−3} and ⟨n_{e}⟩≈1.25 × 10^{−4} cm^{−3}, respectively. Furthermore, this bumpbased iterative boosting algorithm is already providing a similar null hypothesis rejection ability when an overall chisquare test is applied (with a calculation of pvalue according to χ^{2} PDF of null CCF profiles), as the boosted CCF exceeds 3σ in a wide range of scales as well, and at least above 2σ level for the entire range (excluding the very small angles close to the resolution limit of the model map, where the autocorrelation signal rises rapidly).
In the detection case mentioned above, the corresponding CCF profile extracted from the composite ymap (masked NILC ymap + boosted modeled ymap) is also shown in Fig. 8 (turquoise dashdotted line). It is noticeably affected by random correlations within the iterative process, as its shape is not smoothly evolving while boosting (i.e. not conserving its noisedominant profile shape) and this might cause an erroneous interpretation when the profile exceeds 3σ level at any degree. This behavior obviously shows the critical importance of eliminating noise characteristics and instrument dependencies. Nevertheless, the density upper limit value applying the same algorithm is obtained by using this boosting technique (with the composite ymap rather than only the modeled ymap itself). As the intrinsic noise results in confusion (via triggered random damping), to be able to exceed the 3σ level at the bump region the composite map requires more boosting than the default one. However, the difference is not considerable: the resulting upper limit value is roughly 10% larger than the default one (⟨n_{e}⟩≈1.40 × 10^{−4} cm^{−3}). Hence, the main outcome is that the obtained density upper limit does not strongly depend on the chosen boosting technique, which provides a degree of robustness for our obtained limits. Because of this finding, for further analysis (described in Sect. 5.3), a noisefree ymap (boosted modeled filament ymap) is used for null hypothesis rejection.
5.3. Error estimation and robustness check
When computing the power spectrum of a given sky patch, variance in the power value (especially at lower multipoles) must include the effect of cosmicsampling variance. However, in our case, the fullsky C_{l} estimator (NaMaster) is used for obtaining randomlyrealized maps that have similar statistics on the same specific sky patch. Thus, cosmic variance is already included in the error estimates. In fact, as it can be seen from the left panel in Fig. 6, generated power spectra of the simulated maps are already scattering around the measured spectrum of the masked NILC ymap, which thus provides the requisite tolerance of the measured reference power spectrum from cosmicsampling variance.
Therefore, the dominant factors regarding error estimation all concern with the modeling of the WHIM within the filamentary structures (Sect. 4.1): the chosen density profile, fixed filament radius, and assumed temperature isothermality. There is no wellaccepted density profile for the WHIM as it depends on the baryonic physics adopted in the simulations. Hence, a general estimation needs to be done at this point and it ought to be a relatively shallow one which is characteristically seen in many simulations regardless of the used baryonic physics to simulate filaments (Sect. 4.1.3). Then, the βmodel fit is motivated for obtaining a shallow density with a smooth inner region and a constant powerlaw slope outwardly for the diffuse gas. Thus, we just tried to mimic a typical density profile being extracted from hydrodynamical simulations (Sect. 4.1.3), and fitted profile parameters are not altered for the purpose of error estimation for this work. Consequently, the essential drivers for our derivation of density upper limit are the filament radius and the temperature, as these are the main variables affecting the integrated gas pressure of the modeled filaments. Therefore, we constructed a mesh grid for the density values based on these two parameters, which can be seen in Fig. 9. The range of these parameter values was selected according to the noticed scatter in the scaling relations which are mentioned as a part of filament ymap reconstruction (Sect. 4.1).
Fig. 9. Constructed mesh grid being colored by the obtained density upper limits ⟨n_{e}⟩ in m^{−3}. The axes are filament temperature in logarithmic scale (isothermal assumption) versus reference (NGC 5353/4) filament radius in Mpc, which are used while modeling the filaments ymap (left panel in Fig. 2). 
In Fig. 9, the boosting algorithm (Sect. 5.2) is applied for a given parameter combination, and the resulting grid is colored by the obtained mean density upper limits (⟨n_{e}⟩). As expected, the density upper limit increases, while the temperature and radius decrease to be able to produce the same level of correlation signal. The range of obtained ⟨n_{e}⟩ limits is between ≈3 × 10^{−5} cm^{−3} (for T = 10^{6.8} K and r = 1.3 Mpc) and ≈9 × 10^{−4} cm^{−3} (for T = 10^{6} K and r = 0.5 Mpc). Our approach provides a conservative upper limit that enables us to state that it is unlikely to exceed this mean density value for a diffuse WHIM gas within these filaments. To calculate a representative mean, a onetailed test analogy with a 2σ interval was applied by removing the uppermost 5% of the sorted density limit values which represents the highend, leftbottom triangle of the mesh grid in Fig. 9 (greenish and blueish boxes). The calculated mean density limit over the constructed mesh grid after this filtering is ⟨n_{e}⟩≤4 × 10^{−4} cm^{−3}, that is, roughly three times the value given in Sect. 5.2.
The robustness of the density upper limit range is also checked in terms of the Planck Collaboration pipeline that was used for Compton ymap generation and also for the smoothing kernel size (FWHM) in the galaxy density map construction. For the first, MILCA ymap (Planck Collaboration XXII 2016) was used, rather than the NILC ymap, within each step of the overall algorithm. This consisted of masking, simulating randomly realized maps, and boosting until the nullrejection was achieved for each defined mesh grid parameters. Then, a conservative upperdensity limit was obtained by using the abovementioned iterative technique: ⟨n_{e}⟩≤3.9 × 10^{−4} cm^{−3}. For smoothing scales, a galaxy density map is built by using different sizes of kernels spanning from 20′ to 40′ in FWHM and a single density upper limit value is obtained by using a fixed temperature and radius – which are the same as the values noted in Sects. 4.1.1 and 4.1.2 (T = 10^{6.4} K and r = 0.8 Mpc). Then, these upper limits were compared with the reference value given in Sect. 5.2 and the observed deviations are smaller than the scale of ∼0.1 × 10^{−4} cm^{−3}. Accordingly, the resulting differences in the obtained mean density upper limits are minor and, therefore, their effect can be considered negligible because the amount of divergence is smaller than intersquare variations between adjacent cells visualized in Fig. 9. These outcomes allow us to draw the conclusion that the ymap pipeline and smoothing scales do not play a major role in deriving the quantity computed here and our results, presented above, should be robust against such factors.
6. Discussion and conclusions
For the WHIM modeling within the filamentary structures around the Virgo cluster, certain simplifications are made on the characteristic properties whose effects are ignored for a firstorder estimation of the conservative density upper limit values. Below, we give a list of unaccounted factors that might have an imprint on the smallscale correlation strength:

There is a weak correlation between gas mass and temperature in filaments, as described by the MT scaling relation in Gheller et al. (2015), but the temperature is considered independent from the mass of the filaments in our case.

Lengthwise changes in density profile along the spine axis are ignored and homogeneity is assumed throughout our analysis (i.e., no gas clumping).

Fixed density profile was adapted to each filament and the same temperature value was shared to exploit selfsimilar behavior for the filaments, even though Virgo is a dynamically active galaxy cluster. In this way, possible variations on the pressure profile due to internal processes such as feedback from active galactic nuclei (AGNs) and matter accretion shocks (Gheller & Vazza 2019; GalárragaEspinosa et al. 2022) are ignored.

Filaments may have irregular geometries, depending on the environment or their evolutionary stages (Gheller et al. 2015). In our case, they were implemented as smooth cylinders with a fixed diameter in each individual case.

Some amount of the “bound” gas around overdense objects, such as the galaxy halos within the filaments, can also be expected (Gheller & Vazza 2019). This can contribute to an additional correlation signal in the measured CCF. However, this contribution is also well within the noise limit, as already can be seen from Fig. 7. Besides, the derived upper limits in Sect. 5.3 can be considered conservative in the context of this issue, since these values do not account for the extra contribution of bound gases toward the correlation signal.
We assumed a flat ΛCDM cosmology, given by the Planck 2013 results (Planck Collaboration XVI 2014) to convert the redshift listed in the optical galaxy catalog to distances and to construct the 3D filament model. This choice has led to a mean Virgocentric distance of ≈18.3 Mpc. Since there is still some debate around the measured distance value of the Virgo cluster from the Milky Way reference frame (a generally adopted value is ≈16.5 Mpc, see Mei et al. 2007), we ought to consider whether an error in the assumed Virgocentric distance could affect the results from the crosscorrelation analysis. First, we note that even if there is an overall shift in the Virgocentric distance, the sky coordinates remain the same and, hence, our 2D CCF results based on the ymap are unaffected. This also includes the nulltest results for error analysis. However, when converting the results into density limits, there can be a small change in the values due to a resulting change in the volume. We consider this effect to be fairly small, given the roughly 10% difference in the absolute distance scale of Virgo that can be expected and also due to the effect of projection. As such, any changes in the density upperlimit values should be well below the accuracy permitted by the data.
There might be another potential effect from the pipeline used in the identification of galaxies in cosmic filaments around the Virgo cluster. In Kim et al. (2016), 2σ clipping is applied to the galaxies around the fitted spine axis to extract filamenthosted galaxy members. This number of galaxies is taken as a reference parameter in the scaling relations to calculate filament radius rates in Sect. 4.1. Nevertheless, as the galaxy number ratio of filaments is critical at this point rather than the actual member numbers, the clipping procedure does not play a crucial role for the mass proxy via scaling relations in this study. Besides, Gheller et al. (2015) used a criterion based on the mass density threshold in simulation cells for the identification of filaments and their Appendix A shows that the main trends and dispersions obtained for different threshold values are similar and that the scaling relations are consistent as well.
In this study, the main result is a representative range of the density upperlimit values, obtained via varying the crosscorrelation signal amplitude, by means of using models of Compton yparameter, the filament radius, and the temperature. The mean density upper limit at 3σ detection level found in Sect. 5.3 is totally consistent with the WHIM parameter space extracted from simulations such as those shown Fig. 4 in Martizzi et al. (2019) and Fig. 17 in Gheller et al. (2015). The upperlimit range being shown in Fig. 9 is particularly wellmatched with the expected gas density range for filamentary structures (Martizzi et al. 2019; GalárragaEspinosa et al. 2021). This raises the exciting possibility that the next generation of spacebased multiwavelength CMB experiments, such as the proposed LiteBIRD mission (Hazumi et al. 114432F.), would, in fact, be able to measure the Comptony signals from large angular scales and with better precision, as well as to detect the CCF signal that we have discussed from the Virgo filaments (and possibly also from many other nearby cluster systems). That said, this expected improvement is reliant on the final 1/f noise properties of LiteBIRD in intensity. An overall reduction by a factor of 2 − 3 on the combined instrumental noise and systematic noise residuals (astrophysical foregrounds) should be sufficient for that goal, which is within the design criteria of missions such as LiteBIRD. Even if there is a nondetection result as in this paper, the constraining power of the CMB data from these future missions will enable the relaxation of some of the simplifying assumptions mentioned above, to obtain consistent results with the hydrodynamic simulations. Ultimately, this would contribute to the critical evaluation of WHIMformation theories.
Acknowledgments
We thank Franco Vazza for many helpful suggestions regarding baryonic properties of cosmic filaments. We thank Jens Erler, Mathieu Remazeilles, and Frank Bertoldi for their constructive discussions. This work was supported by the National Research Foundation of Korea through grants NRF2022R1A2C1007721 (S.C.R.) and NRF2019R1I1A1A01061237 and NRF2022R1C1C2005539 (S.K.). We specifically thank the anonymous referee for numerous detailed and insightful comments that significantly improved the quality of this paper.
References
 Aguerri, J. A. L., Gerhard, O. E., Arnaboldi, M., et al. 2005, AJ, 129, 2585 [CrossRef] [Google Scholar]
 Alam, S., Albareti, F. D., Allende Prieto, C., et al. 2015, ApJS, 219, 12 [Google Scholar]
 Amaral, A., Vernstrom, T., & Gaensler, B. M. 2021, MNRAS, 503, 2913 [NASA ADS] [CrossRef] [Google Scholar]
 Anbajagane, D., Chang, C., Jain, B., et al. 2022, MNRAS, 514, 1645 [NASA ADS] [CrossRef] [Google Scholar]
 Banerjee, P., Szabo, T., Pierpaoli, E., et al. 2018, New Astron., 58, 61 [NASA ADS] [CrossRef] [Google Scholar]
 Bertone, S., Schaye, J., & Dolag, K. 2008, Space Sci. Rev., 134, 295 [CrossRef] [Google Scholar]
 Brown, S., Vernstrom, T., Carretti, E., et al. 2017, MNRAS, 468, 4246 [NASA ADS] [CrossRef] [Google Scholar]
 Cao, L., Liu, J., & Fang, L.Z. 2007, ApJ, 661, 641 [NASA ADS] [CrossRef] [Google Scholar]
 Cen, R., & Ostriker, J. P. 2006, ApJ, 650, 560 [Google Scholar]
 Danforth, C. W., & Shull, J. M. 2008, ApJ, 679, 194 [NASA ADS] [CrossRef] [Google Scholar]
 de Graaff, A., Cai, Y.C., Heymans, C., & Peacock, J. A. 2019, A&A, 624, A48 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Delabrouille, J., Cardoso, J.F., Le Jeune, M., et al. 2009, A&A, 493, 835 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dunkley, J., Calabrese, E., Sievers, J., et al. 2013, JCAP, 2013, 025 [CrossRef] [Google Scholar]
 Eckert, D., Jauzac, M., Shan, H., et al. 2015, Nature, 528, 105 [Google Scholar]
 Fresco, A. Y., Péroux, C., Merloni, A., Hamanowicz, A., & Szakacs, R. 2020, MNRAS, 499, 5230 [NASA ADS] [CrossRef] [Google Scholar]
 GalárragaEspinosa, D., Aghanim, N., Langer, M., & Tanimura, H. 2021, A&A, 649, A117 [Google Scholar]
 GalárragaEspinosa, D., Langer, M., & Aghanim, N. 2022, A&A, 661, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 GénovaSantos, R., RubiñoMartín, J. A., Rebolo, R., et al. 2005, MNRAS, 363, 79 [CrossRef] [Google Scholar]
 GénovaSantos, R., AtrioBarandela, F., Kitaura, F.S., & Mücket, J. P. 2015, ApJ, 806, 113 [CrossRef] [Google Scholar]
 Gheller, C., & Vazza, F. 2019, MNRAS, 486, 981 [NASA ADS] [CrossRef] [Google Scholar]
 Gheller, C., Vazza, F., Favre, J., & Brüggen, M. 2015, MNRAS, 453, 1164 [Google Scholar]
 Gheller, C., Vazza, F., Brüggen, M., et al. 2016, MNRAS, 462, 448 [Google Scholar]
 Gorski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759 [Google Scholar]
 Hashimoto, D., Nishizawa, A. J., Shirasaki, M., et al. 2019, MNRAS, 484, 5256 [NASA ADS] [CrossRef] [Google Scholar]
 Hazumi, M., Ade, P. A. R., Adler, A., et al. 2020, SPIE Conf. Ser., 11443, 114432F [Google Scholar]
 HernándezMonteagudo, C., GénovaSantos, R., & AtrioBarandela, F. 2004, ApJ, 613, L89 [CrossRef] [Google Scholar]
 Hill, J. C., & Spergel, D. N. 2014, JCAP, 2014, 030 [Google Scholar]
 Hivon, E., Górski, K. M., Netterfield, C. B., et al. 2002, ApJ, 567, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Huber, P. J. 1992, Breakthroughs in Statistics (Springer), 492 [CrossRef] [Google Scholar]
 Hurier, G., MacíasPérez, J., & Hildebrandt, S. 2013, A&A, 558, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kim, S., Rey, S.C., Jerjen, H., et al. 2014, ApJs, 215, 22 [NASA ADS] [CrossRef] [Google Scholar]
 Kim, S., Rey, S.C., Bureau, M., et al. 2016, ApJ, 833, 207 [NASA ADS] [CrossRef] [Google Scholar]
 Martizzi, D., Vogelsberger, M., Artale, M. C., et al. 2019, MNRAS, 486, 3766 [Google Scholar]
 Mei, S., Blakeslee, J. P., Côté, P., et al. 2007, ApJ, 655, 144 [Google Scholar]
 Molnar, S., & Birkinshaw, M. 1998, ApJ, 497, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Mroczkowski, T., Nagai, D., Basu, K., et al. 2019, Space Sci. Rev., 215, 17 [Google Scholar]
 Paturel, G., Petit, C., Prugniel, P., et al. 2003, A&A, 412, 45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Piffaretti, R., Arnaud, M., Pratt, G., Pointecouteau, E., & Melin, J.B. 2011, A&A, 534, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration I. 2011, A&A, 536, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XII. 2013, A&A, 571, A12 [Google Scholar]
 Planck Collaboration XVI. 2014, A&A, 571, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XI. 2016, A&A, 594, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXII. 2016, A&A, 594, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXVII. 2016, A&A, 594, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. VIII. 2013, A&A, 550, A134 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XL. 2016, A&A, 596, A101 [EDP Sciences] [Google Scholar]
 Reiprich, T. H., Veronica, A., Pacaud, F., et al. 2021, A&A, 647, A2 [EDP Sciences] [Google Scholar]
 Shull, J. M., Smith, B. D., & Danforth, C. W. 2012, ApJ, 759, 23 [NASA ADS] [CrossRef] [Google Scholar]
 Sunyaev, R. A., & Zeldovich, Y. B. 1970, Astrophys. Space Sci., 7, 3 [Google Scholar]
 Tanimura, H., Aghanim, N., Douspis, M., Beelen, A., & Bonjean, V. 2019a, A&A, 625, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tanimura, H., Hinshaw, G., McCarthy, I. G., et al. 2019b, MNRAS, 483, 223 [Google Scholar]
 Tanimura, H., Aghanim, N., Bonjean, V., Malavasi, N., & Douspis, M. 2020a, A&A, 637, A41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tanimura, H., Aghanim, N., Kolodzig, A., Douspis, M., & Malavasi, N. 2020b, A&A, 643, L2 [EDP Sciences] [Google Scholar]
 Tanimura, H., Aghanim, N., Douspis, M., & Malavasi, N. 2022, A&A, 667, A161 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tully, R. B. 1982, ApJ, 257, 389 [NASA ADS] [CrossRef] [Google Scholar]
 Urban, O., Werner, N., Simionescu, A., Allen, S., & Böhringer, H. 2011, MNRAS, 414, 2101 [NASA ADS] [CrossRef] [Google Scholar]
 Van Den Busch, J., Hildebrandt, H., Wright, A., et al. 2020, A&A, 642, A200 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Appendix
Correlation function from the angular power spectrum
By following the HEALPix conventions, a bandlimited function f defined over the spherical surface grid can be expanded in spherical harmonics, Y_{lm}, as
where n denotes a unit vector pointing at polar angle θ ∈ [0, π] and azimuth ϕ ∈ [0, 2π). Pixelated f(n) corresponds to sampling f at N_{pix} locations n_{p}, p ∈ [0, N_{pix} − 1]. Then, this sampled function can be used to estimate a_{lm},
where the superscript star denotes a complex conjugation operator and an equal weighting scheme is assumed for each pixel. Then, these coefficients which correspond to Fourier amplitudes can be used to compute estimates of the angular power spectrum as
Using this definition we can rewrite the correlation function in spherical harmonics
where θ is the angle between the two directions and P_{l} is the Legendre polynomial which corresponds to sum over m indexes in Eq. 3 (see Appendix A.4 in Gorski et al. 2005). In HEALPY module, anafast computes the power spectrum (Eqs. A.2 and A.3) and then bl2beam computes the angular profile from its transfer function (Eq. A.4). Besides, we can use synfast to generate a random realization of f(n_{p}) from its input power spectrum C_{l}, which is also used within the null tests in crosscorrelation analysis (Section 4.2).
All Figures
Fig. 1. Galaxy distribution in the seven filaments and one sheet around the Virgo cluster taken from Kim et al. (2016). Different colors represent different structures. Bright (M_{B} < −19) and faint (M_{B} > −19) galaxies are denoted by large open circles and small filled circles, respectively. The large rectangular box is the region of the Virgo cluster enclosed by the EVCC (Kim et al. 2014). The gray dots represent galaxies not associated with the Virgo cluster or a particular structure. 

In the text 
Fig. 2. Modeled noisefree ymap of the filamentary structures around the Virgo Cluster, which is also smoothed with the Planck beam (Sect. 4.1), shown on the left. Projection effects for each stated filament can be realized by comparing with Fig. 1. Constructed galaxy density map from the point locations (Sect. 2.1) by using get_interp_weights() function in HEALPY module and smoothing with a Gaussian beam whose FWHM is ≈30 arcmin, shown on the right. Both maps are centered on the Virgo Cluster. 

In the text 
Fig. 3. NILC ymap of the halfsphere which is centered on the Virgo Cluster (left). Masked version of the left map, described in Sect. 3.2 (right). The color scale denotes the value of the Comptonization parameter, as inferred from multifrequency Planck data. 

In the text 
Fig. 4. LeoIIA filament in threeaxes: right ascension (RA, in degrees), declination (Dec, in degrees), and lineofsight distance (in megaparsecs). Filamenthosted galaxies and uniformly sampled gas particles (to simulate WHIM diffuse gas) are shown with dark and light blue points, respectively. The 3D fitted spine axis is shown with a green curve around which the gas particles are sampled. Distances of filament galaxies are derived from their heliocentric radial velocities (retrieved from HyperLEDA in Kim et al. 2016) based on the assumption of a linear Hubble distance – redshift relationship by applying a standard flat ΛCDM cosmology (Planck 2013 results in Planck Collaboration XVI 2014, Ω_{m} = 0.307, H_{0} = 67.8 km s^{−1} Mpc^{−1}). 

In the text 
Fig. 5. Expected normalized correlation functions for the masked version (Sect. 3.2) of the noisefree (and also without any contamination) fields that are shown in Fig. 2. The blue line points up the autocorrelation function for the modeled filament ymap (left panel in Fig. 2). Red dashed line shows the crosscorrelation function between filament ymap and galaxy density map (right panel in Fig. 2). The characteristic statistical feature of these filamentary structures in terms of correlation degree is obviously seen as a bump in large scales roughly between 30 and 40 degrees, which is shown with a gray band. This scale corresponds to the average separation between primary elongation occurrences of the filamentary structures, as the bump in correlation function indicates that there should be a specific repetitive interval within the analyzed topology. 

In the text 
Fig. 6. Power spectrum comparison between randomly created maps via Algorithm 2 (blue lines) and actual NILC ymap (red line), including the cumulative mask (Sect. 3.2), shown on the left. 1D probability density functions (PDF) of the unmasked pixel values for the NILC ymap (blue bars), one of the simulated maps (with NaMaster) by using the analyzed NILC ymap field itself (green bars), and one of the simulated maps by using the halfdifference map (red bars), shown on the right. Gaussian fits are sketched for the histograms of randomly simulated maps (black dashed lines). Excess of the signal seen in the NILC ymap (blue bars at the leftmost and rightmost parts of the histogram) compared to the random simulation (green bars) in the right panel indicates some measure of nonGaussianity originated from ∼1/10 000 of the total map pixel number, and therefore limits the statistical representability of the underlying field (NILC ymap, at this case) via random realizations obtained from a crosscorrelation method (Sect. 3.1). Together with the left panel, it’s obvious that usage of the halfdifference map is not suitable for the null tests as it’s not providing the matching statistics. On the other hand, by taking advantage of NaMaster to create simulated maps which are obtained by using the own statistics of the analyzed field (right panel in Fig. 3) appears to be the best option with respect to conducting null tests in our case. This is because it mimics the reference field also statistically by reaching a good agreement with the signal characteristics, which can be directly seen from both 1D PDF and power spectrum comparisons. 

In the text 
Fig. 7. Measured crosscorrelation function between masked NILC ymap and constructed galaxy density map is shown with a black dashed line, as a potential indicator of the existence of WHIM in the filamentary structures around the Virgo Cluster. Green and red colored regions correspond to the 3σ range around the mean value (the blue line) being extracted from null tests by using masked NILC ymap and halfdifference map, respectively. In other words, these shaded regions correspond to our benchmark to check the significance of the measured signal. It is obvious that the measured signal cannot reject the null hypothesis (onetailed hypothesis test for the expected positive correlation), even when the statistic of the halfdifference map is exploited, as the signal does not exceed the positive side of the 3σ range in the benchmark. 

In the text 
Fig. 8. Moment of the null hypothesis rejection, in other words, the “detection” case (Sect. 5.2), as the CCF signal of the boosted modeled filament ymap (red dashed line) exceeds the 3σ level at the positive side of the null hypothesis range (green shaded region) for a single point within the characteristic scale (gray band in Fig. 5), where the actual measured CCF signal (dotted black line) is far below this level. In fact, the measured signal appears to be totally lost as it even could not exceed the 1σ range (1,2, and 3σ intervals for the null hypothesis are shown with different colored bars: orange, blue, and green, respectively). The turquoise dashdotted line corresponds to the CCF signal being extracted from the composite map as a summation of the masked NILC ymap and boosted modeled filament ymap. 

In the text 
Fig. 9. Constructed mesh grid being colored by the obtained density upper limits ⟨n_{e}⟩ in m^{−3}. The axes are filament temperature in logarithmic scale (isothermal assumption) versus reference (NGC 5353/4) filament radius in Mpc, which are used while modeling the filaments ymap (left panel in Fig. 2). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.