© The Authors 2023

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 warm-hot intergalactic medium (WHIM) following its expected temperature and density ranges: T_e  ∈  [10⁵ − 10⁷] K and n_e  ∈  [10⁻⁵ − 10⁻⁸] cm⁻³ (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 X-ray bands. Hence, the UV and X-ray 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 so-called “missing baryon problem”.

The onset of wide-area millimeter and submillimeter surveys for detecting the Sunyaev-Zeldovich (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 high-energy 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 X-ray 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 full-sky 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 signal-to-noise 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 Sunyaev-Zeldovich (SZ) map (Tanimura et al. 2020a). Also, detections have been made in the X-ray bands, with wide-area X-ray 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. Cross-correlation techniques can be used in a multi-purposed 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 cross-correlation technique for searching thermal SZ signature induced by hot gas (bound or unbound to clusters) by using CMB data (Hernández-Monteagudo et al. 2004; Génova-Santos et al. 2015). Such correlation-focused CMB-related works already provided some useful information regarding the physical properties of the intergalactic medium. Another example of cross-correlation 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 large-scale 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 X-ray (Eckert et al. 2015; Reiprich et al. 2021) and in the SZ effect (Génova-Santos 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 well-studied 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 cross-correlation study in the maps-space 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 cross-correlation analysis. We simulated the WHIM signal to predict the expected nature of the cross-correlation 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₀ = 67.8 km s⁻¹ Mpc⁻¹) to convert the redshifts listed in the optical galaxy catalog to distances. We discuss in Sect. 6 whether the resulting choice of the Virgo-centric distance could have any influence on our correlation analysis.

2. Data products

2.1. Optical catalog of filament galaxies

The filament-like 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 large-scale structures around the Virgo cluster exploiting the HyperLeda database (Paturel et al. 2003) had been carried out by Kim et al. (2016), whose filament-galaxy 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 (MB < −19) and faint (MB > −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 cut-offs 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 wide-field galaxy-density 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 latitude-longitude 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 noise-free y-map 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 all-sky y-map

The very large extension of the Virgo cluster and its filaments in the sky means that only a space-based millimeter and submillimeter telescope can recover the full extended signal of the thermal Sunyaev-Zeldovich (tSZ) effect and, currently, the best data come from the all-sky 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 all-sky Compton-y-parameter maps (Planck Collaboration XXII 2016), where the y-parameter is the frequency-independent amplitude of the tSZ signal and is defined as the line-of-sight integral of the electron pressure (Sunyaev & Zeldovich 1970):

$$\begin{array}{c}\hfill y={\displaystyle \int }{n}_{\mathrm{e}}\frac{k_{\mathrm{B}}{T}_{\mathrm{e}}}{m_{\mathrm{e}}{c}^2}{\sigma }_T\mathrm{d}l,\end{array}$$(1)

where σ_T is the Thomson cross-section, n_e is the electron number density, T_e is the electron temperature, k_B is the Boltzmann constant, and m_ec² is the electron rest mass energy. A visual representation of the y-map, on the half sky-sphere centered on the Virgo cluster, is shown in the left image of Fig. 3.

Fig. 3. NILC y-map of the half-sphere 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 y-signal 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 largest-scale 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 cross-correlation feature.

3. Cross-correlation analysis

3.1. Spherical harmonics

In the case of a two-point correlation function, structures can be quantified by the correlator using the expectation value integration:

$$\begin{array}{c}\hfill \xi (x)\equiv ⟨f\left(\mathit{y}\right)f(\mathit{x}+\mathit{y})⟩,\end{array}$$(2)

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 two-point 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 two-point 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 one-fourth 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¹, 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² and PolSpice³ 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 cross-correlation function (CCF) from the cross-correlation power spectrum using the standard definition as follows (see Appendix A for details):

$$\begin{array}{c}\hfill \xi \left(\mathit{x}\right)\equiv ⟨f\left(n\right)f\left(n^{\prime }\right)⟩={\displaystyle \sum _l\frac{2l+1}{4\pi }{C}_l{P}_l(cos(\theta )).}\end{array}$$(3)

For our case, where we analyze the cross-correlation 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 cross-power spectrum

$$\begin{array}{c}\hfill {\widehat{C}}_l^{\mathit{yg}}=\frac{1}{2l+1}{\displaystyle \sum _m{\widehat{y}}_{\mathit{lm}}{\widehat{g}}_{\mathit{lm}}^{\ast },}\end{array}$$(4)

where y_lm and g_lm are the Fourier amplitudes which are calculated from the NILC y-map and galaxy-density map, respectively. This estimator does not handle the impact of pixel masking properly (masked pixels are directly considered zero-valued). If conversion from the power spectrum of the pixelated map $C_l^{\text{pix}}$ to the hypothetical unpixelated one $C_l^{\text{unpix}}$ 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 $f_{\mathrm{sky}}^{-1}$ (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 mode-mode coupling kernel is introduced and corrected using PyMaster to obtain an unbiased estimate of the full-sky 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 $\widehat{C_l}$.

For checking the statistical significance of the cross-correlation signal coming from the specified sky patch around the Virgo cluster by implementing statistical null-tests (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 two-point statistics. We note that smoothing-effect (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 low-density 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)⁴ 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₅₀₀). 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₂₀₀ 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 y-map 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)⁴ was used as suggested in Planck Collaboration XXII (2016).

The union of each above-mentioned mask is applied onto both Compton y-map and galaxy density map before cross-correlation 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 top-hat 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 y-map which is centered on the Virgo cluster and the right panel shows a masked version of the same map whose mask is gray-colored.

4. Simulations

Before analyzing the real data to check for the existence of a cross-correlation signal between the NILC y-map and the galaxy-density map, it is worth exploring what kind of signal would be expected if there were no noise and the sky signals (in both Compton-y and optical) only consisted of the Virgo filaments. In this clear representation of the cross-correlation 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 y-map, 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 noise-free filament y-map 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 cross-section 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 y-map reconstruction

4.1.1. Sampling the filament gas particles around spine-axes

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 least-squares 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₁, x₂, ...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 z-axis 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 z-axis 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 sheet-like structure besides the filamentary structures in our sample, which is called the W-M 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 W-M sheet-like structure would be also negligible in our cross-correlation analysis because most of this structure lies within a radius 3 * R₂₀₀ 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 cross-correlation 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 re-check the mass range we get for NGC 5353/4. According to the length of 22 Mpc, the related galaxy mass range is roughly 10¹² − 10¹³ 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 filament-hosted 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², 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 spine-axis for LeoII-A filament can be seen.

Fig. 4. LeoII-A filament in three-axes: right ascension (RA, in degrees), declination (Dec, in degrees), and line-of-sight distance (in megaparsecs). Filament-hosted 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, H0 = 67.8 km s−1 Mpc−1).

4.1.2. Temperature selection for isothermal filaments

After obtaining uniformly sampled particles, the line-of-sight integral of the pressure needs to be calculated to extract Compton-y 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 non-shocked cells that lower the average temperature. In addition, it has been observed that filaments follow well-defined scaling relations in the gas temperature-mass plane (Gheller et al. 2015). Within our sample of dynamically active Virgo filaments, the expected mass range is 10^12.5 − 10¹⁴ M_⊙ and the corresponding average mass-weighted gas temperature range is roughly 10⁶ − 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 high-mass 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 y-parameter calculated over these model parameters is addressed in Sect. 5.3. In the context of normalized cross-correlation function shape, the temperature selection only affects the scaling – and not the signal shape itself.

4.1.3. Density modeling and projection

The Compton-y value in an isothermal gas is simply proportional to the line-of-sight integral of the gas density:

$$\begin{array}{c}\hfill \overline{n}=\frac{\int n\mathrm{d}l}{\int \mathrm{d}l}.\end{array}$$(5)

Numerically, $\overline{n}$ can be easily obtained from uniform point sampling under a given density function n, and the line-of-sight distance can be calculated from the interval between the nearest and furthest particles within each pixel on the map. Regarding the density profile, the beta-model 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 quasi-Monte Carlo integration to obtain an expected density value E(n), which corresponds to $\overline{n}$. The used isothermal β-model the diffuse WHIM gas density profile is written as

$$\begin{array}{c}\hfill n_{\mathrm{e}}(r)=n_{\mathrm{e}}(0){[1+{\left(\frac{r}{r_c}\right)}^2]}^{-3\beta /2},\end{array}$$(6)

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 $\overline{n}$ in Eq. (5).

The final stage to construct a noise-free y-map 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 line-of-sight distance interval is calculated for the obtained set of particles. As a last step, the Compton y-parameter 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 y-parameter 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 y-parameter 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 y-map, observed as a multi-stem 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 cross-correlation function (CCF) between modeled filament y-map 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 auto-correlation function of the filament y-map, 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 noise-free (and also without any contamination) fields that are shown in Fig. 2. The blue line points up the auto-correlation function for the modeled filament y-map (left panel in Fig. 2). Red dashed line shows the cross-correlation function between filament y-map 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 latitude-longitude grid, position pairs for the endpoints were picked randomly by preserving a filament length scale and pixels residing along this stripe between end-points 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. Null-test for random correlations

With the intent of assessing the significance of the measured CCF signal between NILC y-map (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 y-map by using both PyMaster (NaMaster python wrapper) and HEALPY. The pseudocode is shown in Algorithm 2. First, the hypothetical full-sky 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 y-map 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 full-sky auto-power spectrum of the given masked y-map is estimated by using compute_full_master() function (estimates of the full-sky 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 full-sky maps with synfast() function (applies conversion from power spectrum to map via spherical harmonics) and the above-mentioned 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.

We followed this strategy for the null tests to create random masked maps whose statistics are as much as close to the masked NILC y-map 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 y-maps and the one from masked NILC y-map are compared. On the right, 1D PDF of unmasked pixels of the NILC y-map 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 well-suited 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 y-map (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 y-map (blue bars), one of the simulated maps (with NaMaster) by using the analyzed NILC y-map field itself (green bars), and one of the simulated maps by using the half-difference 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 y-map (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 non-Gaussianity originated from ∼1/10 000 of the total map pixel number, and therefore limits the statistical representability of the underlying field (NILC y-map, at this case) via random realizations obtained from a cross-correlation method (Sect. 3.1). Together with the left panel, it’s obvious that usage of the half-difference 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 half-difference maps:

$$\begin{array}{c}\hfill hd=\frac{y_1-y_2}{2}=\frac{(S+N_1)-(S+N_2)}{2},\end{array}$$(7)

where y₁ and y₂ are the y-maps generated out of the NILC-pipeline corresponding to the two jack-knife halves of the data; hd is the half-difference of these two y-maps, which can be used in noise characterization as it gives an estimate of map noise (N₁ and N₂, 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. Plancky-maps 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 y-map and half-difference map is shown on the right panel, with corresponding Gaussian fits for each of them. As can be seen, 1D PDF of the half-difference map is well described by the Gaussian fit, with negligible tails in both directions as expected. For the masked NILC y-map, the Gaussian fit is also well-suited 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 half-difference 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 y-map has a higher variance in comparison with the half-difference map due to either NILC-pipeline or residual foreground contamination, and this residual must be included in the CCF robustness estimation. Hence, the measured power spectrum of the masked y-map 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 half-difference map as well as the NILC y-map to compare a range of the CCFs arising from random correlations. We generated 1000 random realizations of the masked half-difference map and NILC y-map via the above-mentioned 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 null-CCFs is presented for the usage of both NILC y-map and half-difference map by comparing deviation levels up to 3σ. On all scales, scattering is much broader for the case with direct NILC y-map usage, which can be already expected by just looking at 1D PDF results directly in Fig. 6, as the half-difference map manifests a smaller amplitude range and also lower spectral power magnitude.

Fig. 7. Measured cross-correlation function between masked NILC y-map 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 y-map and half-difference 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 (one-tailed hypothesis test for the expected positive correlation), even when the statistic of the half-difference map is exploited, as the signal does not exceed the positive side of the 3σ range in the benchmark.

5. Results

5.1. Non-detection with Planck data

Since there is an expected positive correlation between the modeled filament y-map and constructed galaxy density map as seen in Fig. 5, the significance of the measured CCF between NILC y-map and galaxy density map is checked via a one-tailed null hypothesis test. In Fig. 7, the measured signal appears to be within non-detection case, in other words, the null hypothesis cannot be rejected. Even when only the instrumental noise statistic (from half-ring 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 y-map and galaxy density map is zero.

The considerable variation when using the masked y-map rather than a half-difference 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 null-CCF range getting extracted by using the masked NILC y-map via Algorithm 2 (a green-colored 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 non-detection result to constrain the density upper limits of the hypothesized WHIM phase within the filamentary structures around the Virgo cluster. Based on the null-test 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 y-map (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 y-map and galaxy density map exceeds the 3σ limit which corresponds to the upper boundary of the green-colored region in Fig. 7 (see Sect. 5.1 regarding why only the green-colored region is taken into account as a reference and the red-colored region is thus ignored).

At this point, two options exist regarding the y-map that will be used to calculate the enhanced signal: either the boosted model y-map of the filaments directly or the composite y-map done by adding the boosted model y-map onto the actual measured NILC y-map. 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 χ² of the measured data to the χ² of the probability distribution function of 1000 null samples.

In our case, the boosted modeled filament y-map can directly be used as a reference. In other words, there is no need to contaminate a reference y-map with the noise-dominant y-map (i.e., the masked NILC y-map, 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 y-map can still have contamination of the correlated noise (instrumental or other systematic), as can be seen in Fig. 8 (turquoise dash-dotted 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 y-map (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 dash-dotted line corresponds to the CCF signal being extracted from the composite map as a summation of the masked NILC y-map and boosted modeled filament y-map.

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 p-value calculation according to the chi-square test for null-test rejection is equal to the direct 3σ thresholding. As a consequence, the optimized boosting algorithm is as follows: the modeled filament y-map is initialized with relatively low central density and then it is iteratively increased until the calculated CCF exceeds the 3σ-level of null-range in any point between 30° and 45°. The CCF was calculated between the galaxy density map and boosted modeled y-map directly, not combined with NILC y-map at this point; however, adequate results were also obtained with the composite y-map 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 Newton-Raphson method (first-order in the class of Householder’s methods which are numerical algorithms for root-finding) 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⁻³.

In Fig. 8 (red-dashed line), we can see the moment where the CCF between boosted modeled filament y-map and galaxy density map exceeds 3σ threshold at the positive side (detection case) with a tolerance of Δn_e(0) ≤ 3 m⁻³, where the obtained central and mean (transverse) electron densities are n_e(0)≈3 × 10⁻⁴ cm⁻³ and ⟨n_e⟩≈1.25 × 10⁻⁴ cm⁻³, respectively. Furthermore, this bump-based iterative boosting algorithm is already providing a similar null hypothesis rejection ability when an overall chi-square test is applied (with a calculation of p-value according to χ² 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 y-map (masked NILC y-map + boosted modeled y-map) is also shown in Fig. 8 (turquoise dash-dotted 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 noise-dominant 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 y-map rather than only the modeled y-map 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⁻⁴ cm⁻³). 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 noise-free y-map (boosted modeled filament y-map) 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 cosmic-sampling variance. However, in our case, the full-sky C_l estimator (NaMaster) is used for obtaining randomly-realized 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 y-map, which thus provides the requisite tolerance of the measured reference power spectrum from cosmic-sampling 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 well-accepted 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 power-law 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 y-map reconstruction (Sect. 4.1).

Fig. 9. Constructed mesh grid being colored by the obtained density upper limits ⟨ne⟩ 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 y-map (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⁻⁵ cm⁻³ (for T = 10^6.8 K and r = 1.3 Mpc) and ≈9 × 10⁻⁴ cm⁻³ (for T = 10⁶ 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 one-tailed test analogy with a 2σ interval was applied by removing the uppermost 5% of the sorted density limit values which represents the high-end, left-bottom 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⁻⁴ cm⁻³, 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 y-map generation and also for the smoothing kernel size (FWHM) in the galaxy density map construction. For the first, MILCA y-map (Planck Collaboration XXII 2016) was used, rather than the NILC y-map, within each step of the overall algorithm. This consisted of masking, simulating randomly realized maps, and boosting until the null-rejection was achieved for each defined mesh grid parameters. Then, a conservative upper-density limit was obtained by using the above-mentioned iterative technique: ⟨n_e⟩≤3.9 × 10⁻⁴ cm⁻³. 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⁻⁴ cm⁻³. 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 inter-square variations between adjacent cells visualized in Fig. 9. These outcomes allow us to draw the conclusion that the y-map 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 first-order estimation of the conservative density upper limit values. Below, we give a list of unaccounted factors that might have an imprint on the small-scale correlation strength:

• There is a weak correlation between gas mass and temperature in filaments, as described by the M-T 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 self-similar 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árraga-Espinosa 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 over-dense 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 Virgo-centric 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 Virgo-centric distance could affect the results from the cross-correlation analysis. First, we note that even if there is an overall shift in the Virgo-centric distance, the sky coordinates remain the same and, hence, our 2D CCF results based on the y-map are unaffected. This also includes the null-test 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 upper-limit 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 filament-hosted 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 upper-limit values, obtained via varying the cross-correlation signal amplitude, by means of using models of Compton y-parameter, 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 upper-limit range being shown in Fig. 9 is particularly well-matched with the expected gas density range for filamentary structures (Martizzi et al. 2019; Galárraga-Espinosa et al. 2021). This raises the exciting possibility that the next generation of space-based multiwavelength CMB experiments, such as the proposed LiteBIRD mission (Hazumi et al. 114432F.), would, in fact, be able to measure the Compton-y 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 non-detection 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 WHIM-formation 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 NRF-2022R1A2C1007721 (S.C.R.) and NRF-2019R1I1A1A01061237 and NRF-2022R1C1C2005539 (S.K.). We specifically thank the anonymous referee for numerous detailed and insightful comments that significantly improved the quality of this paper.

References

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

$$\begin{array}{c}\hfill f(n)={\displaystyle \sum _{l=0}^{l_{\mathit{max}}}\sum _m{a}_{\mathit{lm}}{Y}_{\mathit{lm}}(n)}\end{array}$$(A.1)

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,

$$\begin{array}{c}\hfill {\widehat{a}}_{\mathit{lm}}=\frac{4\pi }{N_{\mathit{pix}}}{\displaystyle \sum _{p=0}^{N_{\mathit{pix}}-1}{Y}_{\mathit{lm}}^{\ast }(n_p)f(n_p)}\end{array}$$(A.2)

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 ${\widehat{C}}_l$ as

$$\begin{array}{c}\hfill {\widehat{C}}_l=⟨a_{\mathit{lm}}{a}_{\mathit{lm}}^{\ast }⟩=\frac{1}{2l+1}{\displaystyle \sum _m{|{\widehat{a}}_{\mathit{lm}}|}^2}\end{array}$$(A.3)

Using this definition we can rewrite the correlation function in spherical harmonics

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \xi \left(x\right)\equiv ⟨f\left(n\right)f\left(n^{\prime }\right)⟩& ={\displaystyle \sum _l{C}_l\sum _m{Y}_{\mathit{lm}}(n)Y_{\mathit{lm}}^{\ast }(n^{\prime })}\hfill \\ \hfill & ={\displaystyle \sum _l\frac{2l+1}{4\pi }{C}_l{P}_l(cos(\theta ))}\hfill \end{array}\end{array}$$(A.4)

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 ${\widehat{C}}_l$ (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 cross-correlation 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 (MB < −19) and faint (MB > −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 noise-free y-map 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 y-map of the half-sphere 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. LeoII-A filament in three-axes: right ascension (RA, in degrees), declination (Dec, in degrees), and line-of-sight distance (in megaparsecs). Filament-hosted 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, H0 = 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 noise-free (and also without any contamination) fields that are shown in Fig. 2. The blue line points up the auto-correlation function for the modeled filament y-map (left panel in Fig. 2). Red dashed line shows the cross-correlation function between filament y-map 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 y-map (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 y-map (blue bars), one of the simulated maps (with NaMaster) by using the analyzed NILC y-map field itself (green bars), and one of the simulated maps by using the half-difference 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 y-map (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 non-Gaussianity originated from ∼1/10 000 of the total map pixel number, and therefore limits the statistical representability of the underlying field (NILC y-map, at this case) via random realizations obtained from a cross-correlation method (Sect. 3.1). Together with the left panel, it’s obvious that usage of the half-difference 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 cross-correlation function between masked NILC y-map 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 y-map and half-difference 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 (one-tailed hypothesis test for the expected positive correlation), even when the statistic of the half-difference 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 y-map (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 dash-dotted line corresponds to the CCF signal being extracted from the composite map as a summation of the masked NILC y-map and boosted modeled filament y-map.

In the text

	Fig. 9. Constructed mesh grid being colored by the obtained density upper limits ⟨ne⟩ 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 y-map (left panel in Fig. 2).
In the text