© ESO 2021

1. Introduction

Since the pioneering work of Cen & Ostriker (1999), large-scale structure simulations have systematically predicted that a significant fraction of the cosmic hydrogen at z < 1 resides in a warm-hot (log T(K) ≈ 5 − 7), low density (δ_b ≈ 0–100) phase (e.g. Davé et al. 2001; Dolag et al. 2006; Branchini et al. 2009; Tepper-García et al. 2012; Cui et al. 2012, 2018; Martizzi et al. 2019). As this phase resides in the intergalactic space, it is usually called the warm-hot intergalactic medium (WHIM). Specifically, this WHIM is expected to reside within the filaments of the cosmic web (Bond et al. 1996). This cosmic web consists of a complex and intricate configuration of galaxies, gas, and dark matter, forming a large-scale structure of clusters, filaments, and sheets surrounding large, empty voids. These structures originated as small initial density perturbations that grew driven by gravity to form the non-linear web-like structure we observe today (Zel’Dovich 1970; Shandarin & Zeldovich 1989). For a recent overview on the cosmic web and its study using a variety of methods to detect the large-scale structure, see Libeskind et al. (2018) and the references therein.

According to the above cosmological simulations, the primordial WHIM has been polluted by metals expelled from the galaxies via supernova-driven galactic winds and active galactic nuclei (AGN) feedback. In fact, the low temperature phase (log T(K) ≈ 5 − 5.5) of the predicted WHIM (‘warm WHIM’ hereafter) has been robustly detected via metal ion absorption in the far ultraviolet (FUV) band. For example, a recent analysis of Hubble Space Telescope/Cosmic Origins Spectrograph (HST/COS) data (Danforth et al. 2016) of 82 UV-bright AGN measured hundreds of intergalactic C IV, N V, and O VI absorbers at z < 1, probing a redshift path of Δ_z ∼ 30, tracing ∼11% of the cosmic baryon budget.

The situation is different for the hot (log T(K) ≈ 5.5 − 7) phase of the WHIM (‘hot WHIM’ hereafter). The light elements in the hot WHIM are fully ionised and thus not observable via absorption. The expected X-ray emission of such low-density plasma, and the emission and absorption lines due to the embedded metal ions are too faint to be significantly detected for a sufficientlys large sample of individual filaments with current instruments (see Nevalainen et al. 2019 for a review of the few statistically significant X-ray absorption measurements possibly originating from the hot WHIM and Ahoranta et al. 2020 for a new possible case). Thus, it is generally assumed that the hot WHIM, when robustly detected, will aid in completing the puzzle of the cosmic baryon budget.

Indeed, there is a need for such a component. When accounting for all the detected cosmic baryon components, the total amount falls short of that predicted by the concordance cosmology by ∼30–40% according to Shull et al. (2012). Thus, the hot WHIM is a promising candidate for the missing baryons, which is our basic assumption in this work.

Our aim in this work is to improve the search for the missing baryons by utilising the state-of-the-art EAGLE simulation (Schaye et al. 2015). We explore the spatial distribution of the basic thermodynamic quantities (the temperature and the density) of the intergalactic medium (IGM) with the goal of localising the hottest IGM in observations.

We base our search of the missing baryons on a scenario, whereby the hot WHIM has been heated up via the accretion shocks resulting from the mass flow of baryons towards the gravitational potential maxima of the cosmic web structures containing most of the mass, namely the dark matter filaments (Ryu et al. 2003; Kang et al. 2005; Ganeshaiah Veena et al. 2019). The enhanced density also renders the galaxy formation more effective within the filaments. Thus, the scenario suggests that the missing baryons and galaxies are co-located within the filaments, as indicated by the cosmological simulations mentioned above. We utilise this suggestion in the current work by detecting the filaments based on the galaxies and then examining the diffuse baryon content captured within the filament volumes.

Any cosmic filament detection algorithm faces the fundamental question of what is a cosmic filament. In observations, the cosmic web manifests itself as a network of galaxy filaments connecting dense nodes (galaxy clusters) between empty void regions (e.g. Joeveer et al. 1978; Geller & Huchra 1989). According to theory, the cosmic web is a result of the non-linear evolution of the primordial dark matter density field via gravitational instability. Once the linear growth phase reached a density contrast of the order of one, it triggered the formation of low density (sheet-like) structures. Next, the non-linear accumulation of matter formed caustics which provided the skeleton for the subsequent anisotropic accretion of mass, leading to the formation of filaments and dense virialised nodes. While the nodes are structures defined clearly by the spherical collapse theory, the filaments are not virialised objects and therefore the spherical collapse theory does not help to define them quantitatively. Thus, there is a wide variety of (heuristic) filament finders, employing different methods and assumptions and acting on different mass components (dark matter, gas, or galaxies). Libeskind et al. (2018) have shown consequent significant differences in the mass and volume fractions of dark matter in the filaments detected by a large set of different methods.

In the current work we employ the Bisous and NEXUS+ formalisms as our primary tools for detecting the filamentary network (Stoica et al. 2007, 2010; Tempel et al. 2016; Cautun et al. 2013). Along with the DisPerSE filament finder (Sousbie 2011), the Bisous method has been established as a useful tool for filament detection when working with the spectroscopic galaxy survey data (e.g. Tempel et al. 2014; Libeskind et al. 2015; Poudel et al. 2017; Kuutma et al. 2017; Guo et al. 2015). Also, importantly for our aim of deriving observationally relevant results from the in-depth analyses of the cosmological simulations, the Bisous formalism can be applied to both observational and simulated data (for applications of the Bisous method on the simulations, see e.g. Nevalainen et al. 2015; Libeskind et al. 2018; Ganeshaiah Veena et al. 2019). However, we had to make sure that our results on the properties of the WHIM derived from the Bisous method are not artefacts of the structure finding algorithm. To this end we repeated the relevant parts of our analysis using the filaments obtained by the MMF/NEXUS+ method (Aragón-Calvo et al. 2007; Cautun et al. 2013). For the purpose of this paper, we use the NEXUS+ cosmic web environments detected using the dark matter density field.

Recently, the Bisous and NEXUS+ methods were extensively compared for the case of spin alignments in Ganeshaiah Veena et al. (2019). They showed that filaments detected by both of these formalisms occupy only ≈5–6% of the total volume of the simulation at z = 0, while containing 40–50% of the diffuse baryons. We extend the above work and focus on finding the fraction of the missing hot phase of the cosmic baryons within these filaments.

A study of a single filament in cosmological simulations (Gheller & Vazza 2019) indicated that the density of the diffuse baryons in the full temperature range peaks within the central ∼1 Mpc from the filament axis, by a factor of a few above the cosmic mean. The temperature of the diffuse baryons within the central ∼2 Mpc radii is higher than at larger radii. Our work differs from the above in the sense that we investigate the spatial behaviour of the hot phase of the diffuse baryons for a sample of filaments. This is necessary for estimating the filament-to-filament variation of the radial distribution of the missing baryon properties. This enables a proper error propagation to the estimates of the missing baryon density and temperature at a given distance from a filament detected in the observational data. This in turn will provide a realistic tool for planning the missing baryon observations via X-ray emission and Sunyaev-Zeldovich (SZ) effect.

We use Ω_m = 0.31, Ω_Λ = 0.69, Ω_b = 0.048, H₀ = 67.8 km s⁻¹ Mpc⁻¹, σ₈ = 0.83 and mean baryon density 〈ρ_b〉 = 6.2 × 10⁹ M_⊙ Mpc⁻³ = 4.2 × 10⁻³¹ g cm⁻³. We define the baryon overdensity as δ_b = (ρ_b − ⟨ρ_b⟩)/⟨ρ_b⟩ and the overdensity of the luminosity density field as δ_LD = (LD − ⟨LD⟩)/⟨LD⟩. For plotting logarithmic scales we use the baryon and luminosity density contrasts Δ_b = 1 + δ_b and Δ_LD = 1 + δ_LD, respectively.

2. The EAGLE simulations

2.1. Diffuse baryons

In this work we investigate the diffuse baryons (i.e. gas) in the Evolution and Assembly of GaLaxies and their Environments (EAGLE) simulations (Schaye et al. 2015; Crain et al. 2015; McAlpine et al. 2016). They consist of a suite of hydrodynamical simulations within periodic cubic volumes with a Lambda cold dark matter (ΛCDM) cosmology, assuming the parameters described above (Planck Collaboration I 2014). The simulations used a modified version of the N-body GADGET3 code (Springel 2005), with adjustments to the smoothed particle hydrodynamics (SPH), time steps, and subgrid physics. In particular, the subgrid stellar feedback was calibrated to the observed galaxy stellar mass function and galaxy sizes at redshift z ∼ 0. This, in turn, yielded a good agreement with many other observable properties not used in the calibration.

As our objective is to investigate the baryons within the large-scale structure of the Universe, we selected the gas particles in the largest simulation, REFL0100N1504. Within a co-moving volume of 100³ Mpc³, it followed 1504³ dark matter and baryon particles, with initial masses of m_DM = 9.7 × 10⁶ M_⊙ and m_b = 1.81 × 10⁶ M_⊙, respectively. While the number of dark matter particles and their masses remained the same throughout the simulation, the mass of baryon particles were allowed to vary. Some baryon gas particles were converted into stellar particles or swallowed by black holes, while other gas particles received mass from stellar particles via feedback effects. In addition to mass, each of the particles carried a variety of information throughout the simulation. In this work we utilise the coordinates, masses, temperatures, and densities of each SPH particle. Furthermore, each particle holds information on the bound structures they belong to, if any, such as galaxies and groups of galaxies.

In order to obtain a general overview of the large-scale distribution of the EAGLE baryons, we first computed the baryon densities within a representative 5 Mpc thick slice in a grid with a cell size of 0.1 Mpc. We subsequently projected the densities into a plane to obtain a 2D density distribution (see Fig. 1). Consistently with earlier simulation works, the distribution indicates the large underdense voids and the filamentary cosmic web at overdensities of δ ∼ 10–100.

Fig. 1. Projected diffuse baryon density within a representative 5 Mpc slice (colour map) above the mean baryon density and galaxies brighter than Mr = −18.4 (blue dots) in the EAGLE simulation. The light blue colour denotes the regions below the mean baryon density. The density was computed by co-adding the masses of individual gas particles and dividing by the projected volume in each pixel.

Throughout this work all the 2D projections were computed using the aforementioned grid. Unless otherwise stated, all values, distributions, and profiles have been computed directly from the particle data.

2.2. Galaxies

We utilised the EAGLE galaxy data (McAlpine et al. 2016) for experimenting with the filament detection (see Sect. 1). For each galaxy we selected the coordinates, the dust corrected r band magnitudes, and the virial radius R₂₀₀ of the FoF group the galaxy belongs to (the radius where the mean internal density is 200 times over the critical density of the Universe). Within the virial radius the galaxies contain a stellar mass which constitutes ≈3.4% of the total baryon budget in the EAGLE simulation.

Since our aim is to provide observationally relevant results, we must consider that EAGLE is volume-limited while the observational surveys are typically flux-limited. Our particular aim in the near future is to apply the results of this work to the Sloan Digital Sky Survey (SDSS, Ahn et al. 2014) up to redshift z = 0.05, since our adopted Bisous filament detection algorithm works well up to this distance (Tempel et al. 2014). At z = 0.05, the limiting SDSS magnitude corresponds to M_r ≈ −19. Unfortunately the luminosity distribution in EAGLE is somewhat different from SDSS (see Fig. 2). Namely, the EAGLE distribution decreases faster than that of SDSS with higher luminosities. The behaviour is similar as when comparing the EAGLE simulations with the GAMA data at z = 0.1 (Trayford et al. 2015). Thus, the EAGLE galaxy density is smaller than that of SDSS, if using the same magnitude limit. However, the galaxy density is the primary parameter determining the performance of our filament detection algorithm (M. Muru, priv. comm.). Thus, to approximate the SDSS survey at z  ≤ 0.05, we have to include fainter EAGLE galaxies by using M_r = −18.4 to match the SDSS M_r = −19 galaxy density. With this fainter cut we obtained a sample of filaments for the EAGLE simulation with similar statistical properties as those of the SDSS filaments. We assume that the slight mismatch of the EAGLE and SDSS magnitude limits does not produce complications in our analysis, since the spatial distributions of the non-overlapping populations are expected to be similar: the galaxies are mostly in groups. A 2 Mpc smoothing length when computing the luminosity density field spreads the entire group luminosity into the IGM regions that are of interest in this work (see Sect. 5.3).

Fig. 2. Luminosity distributions of EAGLE (yellow points and error bars) and SDSS (red points and error bars) galaxies. The vertical dashed lines indicate the selected magnitude cuts of −18.4 and −19 for EAGLE and SDSS galaxies, respectively.

Visual inspection indicates that the galaxy distribution follows the filamentary diffuse baryon distribution well (see Fig. 1). This is consistent with Nevalainen et al. (2015), who reported a correlation between the density of the WHIM and the galaxy luminosity density in another simulation. This behaviour is consistent with the role of the dark matter in our adopted scenario, whereby the dominating dark matter gravity at the filaments enhances the baryon and galaxy density.

3. The missing baryons

3.1. The problem

Shull et al. (2012) presented a popular census of the observational status of the cosmic baryon budget at low redshift, which indicated that ∼30 ± 10% of the baryons predicted by the concordance cosmology have not been detected. Later developments in the FUV measurements of AGN with HST/COS (Danforth et al. 2016) have refined the above account. The H I-traced phase (Lyman-Alpha forest and broad Lyman-Alpha absorbers (BLAs)) constitutes 25 ± 5% of the baryon budget at z = 0.0 − 0.5 (see Table 1), which is smaller than estimated in Shull et al. (2012) (42 ± 13%). Assuming that the observational status of other significant baryon components (i.e. stellar mass and CGM in galaxies and the diffuse intergalactic baryons within galaxy groups and clusters) is as estimated in Shull et al. (2012), the Danforth et al. (2016) estimate of the Lyman-Alpha (Lyα) abundance yields a directly observational baryon budget of ≈40% of the expected cosmic density.

In addition to the above directly observational cosmic hydrogen, there is also observational evidence of warm-hot log T(K) ≈ 5.0 − 5.5 cosmic baryons via extragalactic O VI absorbers (e.g. Danforth et al. 2016). In the above work the cosmic contribution of the baryons traced by O VI is reported as 11 ± 1%. There are several uncertainties with the usage of this result in the cosmic baryon budget.

First, the above conclusion requires information on the temperature and metal abundance of the gas so that the O VI measurements can be converted into the underlying dominating baryon component, the hydrogen mass. This information is not attainable via O VI line measurements only. The solution used in Danforth et al. (2016) (and in similar earlier works by Tilton et al. 2012 and Shull et al. 2012) is to utilise the correlation between the ion fraction, metal abundance, and the ion column density from simulations. This complicates the interpretation of the derived hydrogen population as observed cosmic baryons.

Second, the reported ≈10% statistical precision indicates that the uncertainties of the temperatures and metallicities are relatively small. Observational estimates do not exist and the scatter of the correlation due to cosmic variation, available in the simulations, is not propagated to the uncertainties. Thus, the reported uncertainties are underestimated.

Third, O VI traces hydrogen within the same temperature range as the BLA measurements (i.e. log T(K) = 5.0 − 5.5). If the O VI-traced baryons and those directly detected via BLA are treated independently, we double-count the overlapping, uncertain population.

Thus, the O VI-traced contribution lies between 0% (full overlap with the Lyα observations) and ≈10% (completely independent component), yielding a total missing baryon fraction within the range of ≈50 − 60%. This would be a significant overestimate if the true O VI – traced cosmic baryon density is much higher than the reported 11%, which is possible given the underestimated uncertainties.

3.2. Our approach

In this work we utilise the hypothesis that a large fraction of the missing baryons reside in the poorly observed hot gas phase. Assuming collisional ionisation equilibrium, the ionisation fraction of the most important FUV ion O VI peaks at log T(K) ≈ 5.5 and decreases rapidly with higher temperatures where the most important X-ray ion O VII dominates. The transition is not clear-cut, as shown by Nelson et al. (2017) in the overlap of temperature-dependent O VI and O VII ion mass distributions around log T(K) = 5.5 at typical filament baryon overdensitisses of 10–100 in IllustrisTNG simulations. In order to minimise the robustly detected FUV phase (Danforth et al. 2016), in our work we adopt log T(K) ≥ 5.5 as one criterion for the missing baryons. Approximately 42% of the EAGLE baryon content is in the form of gas above this temperature limit (see Table 2).

3.2.1. The bound phase

Previous simulations (e.g. Davé et al. 2001; Dolag et al. 2006; Branchini et al. 2009; Tepper-García et al. 2012; Cui et al. 2012, 2018; Martizzi et al. 2019) have indicated that a significant hot baryon reservoir is located within the filaments of the cosmic web. We utilise the results of the Friends-of-Friends (FoF) procedure to separate the more easily detectable hot and dense phase bound to galaxies and galaxy groups¹ from the purely intergalactic phase as follows.

The FoF procedure has been applied by the EAGLE team to the dark matter particle distribution to define through constant density boundaries the dark matter haloes. A gas particle whose neighbouring dark matter particle belongs to a FoF halo has been defined as a halo particle. We assume that such particles form the bound phase of baryons. We separate the bound phase into two environments: (1) the virial phase, that is gas particles located within the virial radius R₂₀₀ of a FoF halo given in the public EAGLE database (McAlpine et al. 2016) and (2) the outskirts phase, gas particles occupying the FoF halo outside R₂₀₀.

In order to understand the baryon properties of the bound phase, we examined the mass distribution of the gas within FoF haloes as a function of the temperature. We divided the temperature range of log T(K) = 3 − 8 into 50 bins of log T(K) = 0.1 width. We then added the mass of gas particles whose temperatures were within the boundaries of a given bin. Subsequently we divided the resulting masses with the temperature bin size, thus obtaining the distribution dM/dT(T). We normalised the distribution to the total gas mass (see Fig. 3). We repeated the process separately for the gas within the virial radius and the gas in the outskirts of galaxies and groups of galaxies to compare the different environments.

Fig. 3. Left panel: distribution of the full EAGLE baryon gas mass content as a function of the temperature, normalised to unity at log T(K) = 3. Right panel: corresponding distributions of different components, relative to the full distribution: the bound phase (orange dashed line), the virial phase within R200 (blue dotted line), and the gas on the outskirts of galaxies and groups of galaxies (outside R200 but within FoF halo, green dash-dotted line). The black vertical line indicates the approximate lower limit for the hot WHIM (i.e. log T(K) = 5.5).

The temperature distribution of the bound phase indicates that its contribution to the total baryon budget at log T(K)  ≤  6 is smaller than 10% (see Fig. 3, right panel). Towards higher temperatures the relative importance of the bound phase increases: at log T(K) > 6.5 the bound phase contains ≈75% of the EAGLE baryon mass. We assume that this is due to the gravitational collapse shocks involved in heating the bound baryons being more energetic than the accretion shocks heating the intergalactic baryons. However, since the absolute amount of baryons decreases rapidly towards the highest temperatures (see Fig. 3, left panel), the contribution of the bound phase to the baryons in the missing baryon temperature range (log T(K) > 5.5) is only ≈44%. Thus, there is ample room for a large fraction of the missing baryons to reside in other environments (particularly in filaments, see Sect. 5).

At the low end of the missing baryon temperature range (log T(K) = 5.5 − 6.0) the outskirts contribute ≈50% to the bound phase (see Fig. 3). At higher temperatures the gas within R₂₀₀ dominates the bound phase. Thus, according to the EAGLE simulation, the contributions of the outskirts and R₂₀₀ to the baryon budget in the missing baryon temperature range (log T(K) > 5.5) are ≈12% and ≈32%, respectively. Their respective contributions to the full EAGLE baryon budget are ≈5% and ≈13% (see Table 2). Stars within R₂₀₀ contribute an additional ≈3% to the total baryon mass. In our filament analysis below, we include the hot baryons of the outskirt regions in our missing baryon sample and exclude only the baryons residing within R₂₀₀.

To understand the bound phase better, we produced temperature-density phase diagrams as follows. We binned the density contrast and temperature into 70 bins each (see Fig. 4). The range for the overdensities was Δ_b = 0.01–10⁵ (bin width of log Δ_b = 0.1) and for the temperature log T(K) = 2 − 9 (bin width of log T(K) = 0.1). Then, we summed the masses of the particles that fell within the ranges of the bins. Finally, the mass in each bin was divided by the total mass of all the gas particles to calculate the corresponding mass fractions. In addition to the full EAGLE gas content, we performed the above analysis for particles within R₂₀₀ and outskirts separately.

Fig. 4. Baryon density contrast – temperature phase diagrams (for the convenience of plotting we use the density contrast Δb = 1 + δb rather than the overdensity δb). All the gas in the full EAGLE volume is shown in the top left panel, the gas within R200 in the top right panel and the gas within the Bisous filaments excluding gas within R200 in the bottom left panel. For each temperature and density bin we co-added the mass of particles and divided by the total mass of the full box. We cut the colour scale at 1% level of the maximum value for clarity. The bottom right panel indicates the isomass curves containing 68% of the baryon mass in Bisous filaments excluding R200 (blue line), within the outskirts (black line), and within R200 (red line). The low-density, low-temperature track seen in the upper left panel corresponds to the photo-heated IGM.

The temperature-density phase diagram of the bound phase indicates the expected behaviour that the denser gas is hotter (see Fig. 4). The concentration of virial gas above log T(K) = 6 (discussed above) corresponds to a density concentration above Δ_b ∼ 100. This phase is relatively compact and dominates the hottest and densest domain. The outskirt gas has similar behaviour, but at a bit lower temperatures and densities (see Fig. 4). The outskirts gas phase matches the dense, high-temperature part of the Bisous filament phase (see Sect. 6.2 for more discussion). A small fraction of the diffuse baryons (1.8%) is locked within the central regions inside and around galaxies. This low temperature – high density phase is seen as a low-level deviation from the main component of the virialised gas (see upper right panel of Fig. 4).

3.2.2. The missing baryon sample in EAGLE

The virial phase behaves smoothly in the temperature-density plane (see Fig. 4) and dominates the hottest and densest domain. Thus, the exclusion of the virial phase from the full EAGLE diffuse baryon population would provide a simple way of separating the more easily detectable high density and high temperature gas from the missing baryons. After the additional low temperature (log T(K) < 5.5) cut, only ≈7% of the total baryon mass is outside the canonical WHIM baryon density interval of 1–100 times the mean. Similarly, only ≈2% of the total baryon mass exceeds the canonical WHIM upper temperature limit of log T(K) = 7.

While being simple, the exclusion of all the hot virial gas in EAGLE from the missing baryon category would yield a somewhat underestimated missing baryon fraction. Namely, assuming that the extrapolation of the observed hot gas mass in the central regions of galaxies to the virial radii (e.g. Bregman et al. 2018) is accurate, most of the virial gas within galaxies is currently unobserved (i.e. missing). Assuming that (1) the observed CGM contribution of ≈5% to the cosmic baryon budget (Shull et al. 2012) accurately corresponds to the observed hot gas in galaxies discussed in Bregman et al. (2018), and that (2) EAGLE accurately describes the contribution of the hot gas within R₂₀₀ to the cosmic baryon density (≈13%, see Sect. 3.2.1), by excluding the gas within R₂₀₀ we would ignore a missing baryon component amounting to ≈8% of the cosmic baryon budget.

However, our focus in this work is the intergalactic WHIM and thus we do not analyse the possible missing baryons inside R₂₀₀ further (for a thorough analysis of the CGM in EAGLE see Wijers et al. 2020). For the purposes of this work, we adopt the following definition of a missing baryon: a) its temperature is higher than 10^5.5 K and b) it is outside the virial radius (R₂₀₀) of collapsed structures (i.e. we consider that the outskirts phase belongs to the missing baryons category).

The baryons satisfying these criteria (P_true) correspond to ≈29% of the total baryon mass in the EAGLE simulation. This substantial reservoir of baryons that are currently difficult to observe is a significant share of the missing baryon fraction of ≈50–60% estimated with the baryon census of Shull et al. (2012) and Danforth et al. (2016) (see Sect. 3.1). This supports our assumption that the hot WHIM is a major component of the missing baryon population. The currently non-observable baryon population within R₂₀₀ of galaxies, if amounting to ≈8% of the cosmic baryon density (see above), would increase the missing baryon fraction in the form of hot gas to ≈37% according to EAGLE. While this level is lower than the formal missing baryon fraction census of ≈50–60% (Shull et al. 2012; Danforth et al. 2016), the missing baryon census might be significantly overestimated due to the possibly underestimated O VI-traced cosmic baryon density (see Sect. 3.1).

4. Overall spatial distribution of thermodynamic properties

In this section we utilise the information on the spatial distribution of the baryons to improve the search for the missing ones.

4.1. Temperature map

Using the same representative 5 Mpc thick slice as in in Fig. 1, we computed the mass-weighted mean temperatures of each particle in a projected 2D grid with a resolution of 0.1 Mpc. In order to visualise the temperatures of only the intergalactic medium, we excluded all particles within the virial radius R₂₀₀ of haloes. The resulting temperature map (Fig. 5) reveals a significant structure. Most of the volume appears to have low temperatures (log T(K) ≤ 5), presumably corresponding to the voids and sheets.

Fig. 5. Mass-weighted mean temperature distribution (in shades of blue) of the intergalactic gas (i.e. gas outside R200) within the same slice as Fig. 1. Red dots denote the location of galaxies brighter than Mr = −18.4.

Importantly for the current work, at the WHIM temperatures (log T(K) = 5 − 7) the cosmic web – like structure appears. This is expected from previous simulations (e.g. Davé et al. 2001; Dolag et al. 2006; Branchini et al. 2009; Tepper-García et al. 2012; Cui et al. 2012, 2018; Martizzi et al. 2019) and supports our approach of searching for the hottest WHIM within the cosmic filaments. The temperature structures have very sharp edges, revealing the accretion shock heating processes which are essential in our adopted theoretical scenario.

At the highest temperatures (log T(K) ≥ 7), the mass is concentrated around the crossing points of the apparent filaments traced by the log T(K) = 5 − 7 gas. They are the sites of the most energetic gravitational collapse induced shock heating.

4.2. Densities of different phases

We repeated the procedure of using the representative slice of the simulation box from Sect. 2.1, but this time constraining the baryon temperatures of gas outside R₂₀₀ into different ranges. Visual inspection (see Fig. 6) yields the following suggestions. We investigate the suggested behaviour of the temperature distribution in different cosmic web environments quantitatively with NEXUS+ as a function of temperature in Sect. 6.3. In order to roughly approximate the volume filling fractions, we computed the mass-weighted mean temperatures of the particles in the representative slice (see Fig. 1) in a 3D grid, with a separation of 0.1 Mpc between grid points. The volume filling fraction for each environment was then calculated by dividing the number of grid points associated with the respective temperature range by the total number of grid points in the representative slice.

Fig. 6. Projected diffuse baryon mass density contrast within the same 5 Mpc slice as Fig. 1, but here in different temperature ranges. Top left: log T(K) < 3. Top right: log T(K) = 3 − 4. Middle left: log T(K) = 4 − 5. Middle right: log T(K) = 5 − 5.5 (the warm WHIM phase). Bottom left: log T(K) > 5.5 (the missing baryon phase). Bottom right: Gas within the virial radii R200 of haloes, excluded from the other panels.

The baryons with log T(K) < 3 trace the lowest-density parts of the voids. The absence of such cold gas between the voids suggests the cosmic web structure. This phase occupies the largest volume at ∼55% of the total volume. At log T(K) = 3 − 4 the volume filling factor is smaller, ∼23%, consistent with the assumption that this temperature range presents the hottest (i.e. rare) void gas together with sheets which are less volume filling but warmer than the voids. The baryons with log T(K) = 4.0 − 5.0 correspond to the bulk of the intergalactic phase detected via the Lyman alpha forest. The volume filling fraction (∼12%) is smaller than at the lower temperatures since the neutral hydrogen, accumulated towards the more compact filaments and sheets, becomes important.

In the warm FUV-detected range of log T(K) = 5.0 − 5.5 the mass is much more concentrated, as demonstrated by the smaller volume filling fraction of ∼5%. While the sheets still contribute to this temperature range, the mass is now more concentrated around the apparent galaxy filaments which are warmer than the sheets due to the low energy part of the filament accretion shock-heated gas.

Finally, in the missing baryon range (i.e. log T(K) > 5.5 excluding gas within R₂₀₀) the diffuse gas is compact with a volume filling fraction of ∼5% and very strongly concentrated to the apparent filaments. This is consistent with our adopted scenario whereby this phase is produced by the strongest accretion shocks due to mass flow towards filaments. The virial structures as defined in Sect. 3.2.1 (dominating the mass budget at the highest temperatures of log T(K) > 6) apparently trace the same general large-scale web structure as the WHIM, but have a very narrow and clumpy spatial distribution within the apparent filaments. This phase is produced by the gravitational collapse and mergers of dark matter haloes.

5. The filaments

5.1. The detection method (Bisous)

Our goal is to investigate the properties of the missing baryons in the environment of the cosmic filaments. To enable this, we applied the Bisous formalism to EAGLE galaxies to detect the filaments (see Stoica et al. 2007, 2010; Tempel et al. 2016 for a detailed description of the method). It is designed for detecting the large-scale filamentary networks based on spectroscopic galaxy catalogues (both observational and simulated). We used the software to model the distribution of the EAGLE galaxies brighter than M_r = −18.4 (see Sect. 2.2) by fitting connected and aligned cylinders to it. We repeated this 1000 times to account for the stochastic nature of the Bisous formalism. We thus obtained a map for each run, indicating the spatial volumes covered by the fitted cylinders (hereafter we refer to this map of filament coverage as a ‘visit map’). The regions which get covered more often in the individual fits are more likely to belong to a filament. Our experimental definition of a spatial location belonging to a filament is that it is covered by more than 5% of the individual visit maps, since this choice adequately yielded the statistical properties of the SDSS filaments (Tempel et al. 2014). We define the filament spines as the local maxima of the combined visit map, in other words where most of the separate visit maps overlap.

In this work, we employed the same version of the Bisous as used for detecting the filamentary galaxy network in SDSS (Tempel et al. 2014) (see the above references for details of the method). Here, however, the cylinder radius is not fixed to 0.7 Mpc, as in Tempel et al. (2014), but is selected to be in the range 0.5–1.0 Mpc, as in Kuutma et al. (2020). Thus, the Bisous approach adopted in this paper is not optimised for detecting gas or the missing baryons, but to construct the filamentary network from the distribution of galaxies. In fact, a recent study by Kooistra et al. (2019) found that the aforementioned radius for the Bisous filaments was appropriate to obtain a reliable signal-to-noise ratio for neutral H I observations. Our aim in this work is to evaluate the efficiency of the current version in tracing the hot WHIM and, if promising, we will develop it for observations in future work.

5.2. Results

The visit map values (see Sect. 5.1) were stored in a 3D grid with a separation between grid points of 0.4 Mpc². The filament volume covered by the visit map ≥0.05 occupies ≈5% of the EAGLE volume (see Fig. 7), hosting ≈53% of the baryon gas (including gas within collapsed structures). Excluding baryons within R₂₀₀, the Bisous filaments contain ≈40% of the total baryon budget. The filament spines were stored as a catalogue of spine points with their respective xyz-coordinates. In order to avoid possible edge effects, we excluded all the filament spines for which any spine point is closer than 3 Mpc to the edges of the simulated volume. The resulting sample contains a total of 779 filament spines (see Fig. 7).

Fig. 7. Same slice as in Fig. 1, but for a visualisation of the Bisous formalism. Top left: the galaxies (blue dots) and the visit map ≥0.05 (grey). Top right: the spatial distribution of the missing baryons (colour map) and the filament volumes (grey areas). Bottom left: the spatial distribution of the missing baryons captured with the Bisous filament finder. Bottom right: the spatial distribution of the missing baryons (colour map), the filament spines of our high luminosity overdensity δLD sample (blue), and filaments excluded from the high δLD sample (magenta).

The current Bisous set-up is tuned for detecting filaments at an approximate scale of radius r ≈ 1 Mpc (Tempel et al. 2014). We used the data in hand to check the accuracy of the filament scale. Visual inspection of the Bisous volumes (see Fig. 7) indicates quite a complex geometry. Filaments cross and curve three dimensionally with both narrow and wide regions. Thus, to investigate the filament scale issue, we computed the stacked ratio of the grid points determined to be within a filament to the total number of grid points as a function of distance from the filament spines (see Fig. 8). We found that at r = 1 Mpc distance from the filament spine on average ≈70% of the volume probed by the grid points is within filaments. Subsequently the distribution experiences a steep drop, and at r = 1.5 Mpc only ≈10% of the grid points remain within the Bisous filaments. This quantitative result verifies that the current Bisous set-up is correctly fine-tuned for detecting filaments at scales of r ∼ 1 Mpc.

Fig. 8. Fraction of grid points within Bisous volumes (visit map value > 0.05) as a function of distance from the filament spines. All spines were stacked together.

5.3. Classification

The filament length is one of the few properties which can be derived relatively easily from the spectroscopic galaxy surveys. However, the different methods for detecting the filaments yield systematic differences in the basic filament properties (see e.g. Libeskind et al. 2018; Malavasi et al. 2020). Also, by varying the parameters of a given filament detection method, one may obtain systematic differences in the extracted populations. Thus, the length of a filament spine is generally not a robust way of classifying the filaments.

Particularly, the Bisous formalism we are employing in this work with the parameter set optimised for SDSS may break a given structure, defined as a single filament by the DisPerSE method (Sousbie 2011), into several shorter segments. Also, the EAGLE simulation box size of 100 Mpc will suppress structure formation modes larger than this. Consequently, the maximum Bisous filament spine length in this work is ∼35 Mpc, shorter than in Malavasi et al. (2020).

The gas density is a better filament measure, especially for studying the WHIM, as it is assumed to follow the local DM density. While the gas density is not directly available from the basic filament detection analysis, we can use the easiest and most direct measurable property of the cosmic filaments, the galaxy luminosity, as a proxy. Namely, Nevalainen et al. (2015) proposed that the luminosity density (LD) field of the galaxy population related to a given filament can be used as a low-scatter proxy for the WHIM density. The LD field method consists of three-dimensional smoothing of the spatial distribution of the luminosities of a given galaxy sample and evaluating the field at desired locations (see Liivamägi et al. 2012 for details of the method).

Indeed, to first order, we expect galaxies to trace the gas in filaments because they form in haloes collapsed into the filaments due to the same gravitational force as the more diffuse IGM. At the same time, the galaxy luminosities trace the masses of their parent haloes, and haloes of different masses form in regions with different large-scale over- and underdensities. Therefore, we utilised this scenario by employing the LD value as a filament classification criterion. We applied the LD field method to the r-band luminosities of the EAGLE galaxies brighter than M_r = −18.4 (see Sect. 2.2), using a smoothing kernel width of 2 Mpc (see Fig. 9). We then evaluated the LD field at the locations along the spines of the 779 filaments detected in the EAGLE box. We divided the LD values by the average LD of the full EAGLE box (〈LD〉 = 7.1 × 10⁷ L_⊙ Mpc⁻³) and consequently discuss the luminosity overdensity δ_LD = LD/〈LD〉 − 1. For each filament we then calculated the δ_LD value by averaging the individual δ_LD values of the spine points.

Fig. 9. Galaxy luminosity density contrast (ΔLD) in the same slice as in Fig. 1, using a 2 Mpc smoothing kernel. Filament spines of the high δLD sample are shown in blue, the excluded filaments are shown in red.

The resulting δ_LD distribution (see Fig. 10) exhibits a well-behaved decreasing trend towards high values up to δ_LD ∼ 50. Above this, there is a low level tail extending to δ_LD ∼ 200. Most of the highest values are due to filament spine fragments in galaxy groups. Our Bisous formalism is not optimal for working in such environments: some spines are spuriously created due to the local overdensity of galaxies without any filamentary structures being present. Thus, we experimented by limiting the further analysis to filaments with δ_LD ≤ 50.

Fig. 10. Distribution of filaments as a function of the luminosity overdensity in SDSS (green line) and in EAGLE (purple line, normalised to SDSS volume). The vertical lines indicate the limits for the low (δLD = 0 − 10), medium (δLD = 10 − 20) and high (δLD = 20 − 50) luminosity overdensity groups (see Sect. 5.3).

We then experimented by combining a modest length cut for the filament spines together with the luminosity overdensity δ_LD. We found that a 2 Mpc cut is optimal: a deeper length cut reduces the number of selected filaments, while a more relaxed cut increases scatter due to spurious filament spines.

After excluding all the filament spines shorter than 2 Mpc, we divided the remaining spines into three groups based on their δ_LD: low (δ_LD < 10), medium (δ_LD = 10 − 20), and high (δ_LD = 20 − 50). The low, medium, and high overdensity groups contain 432, 134, and 77 filament spines, respectively. The relative overdensity distribution shifts towards higher values with higher luminosity density (see Fig. 11), consistent with the relation between filament gas density and luminosity density, as described by Nevalainen et al. (2015). Since the observational signal improves with higher density, the above indicates that it is beneficial to concentrate the observational efforts on the high δ_LD filaments.

Fig. 11. Top: absolute (left panel, normalised to unity at the maximum value of the full baryon sample, left panel in Fig. 3) and relative distributions (right panel) of the baryon gas mass distribution as a function of mass-weighted temperature in filaments with low (orange dashed line), medium (blue dotted line), and high (green dash-dotted line) luminosity overdensity. For comparison, mass within all Bisous filaments is shown with a black line. Bottom: same as top panel, but for the gas mass distribution as a function of gas density contrast. Gas within R200 has been excluded in all panels.

Bearing in mind that we want to develop methods to be applicable to observational data, we experimented whether the above LD classification scheme would work for the filaments detected in SDSS (Tempel et al. 2014). We selected the same sample of filaments as published in Kuutma et al. (2020), except for the detail that we excluded filaments shorter than 2 Mpc, not 5 Mpc (see the above work for the details of the analysis). Adjusting linearly with the difference of the volumes covered by SDSS in the adopted redshift range of 0.02–0.05 and EAGLE, we found that the δ_LD distributions for SDSS and EAGLE are very similar (see Fig. 10). At the lowest luminosity densities the SDSS filament numbers agree with the scaled EAGLE values within 15%. While both samples display a decreasing trend with higher luminosity density, SDSS drops a bit faster, resulting in a ∼20% deficit of the filaments in the high LD group. The larger difference in the luminosity distributions between EAGLE and SDSS galaxies (see Fig. 2) is not seen in the filament distribution, since most of the galaxies have magnitudes in the range where the two luminosity distributions agree. Thus, we assume that the deficit in SDSS filaments is small enough that the conclusions derived for the high δ_LD group in this paper are relevant for SDSS and that our filament classification scheme is applicable to SDSS at low redshifts. From a total of 5992 filament spines, the number of the SDSS filaments in the three δ_LD groups are 4247 for the low, 1190 for the medium, and 498 for the high overdensity. We will perform a stacking analysis of the Planck data related to these filament groups in a future work.

6. Capturing the missing baryons with the filament finder

6.1. Capture fraction

We extracted the total EAGLE baryon population within the full filament sample (see Sect. 5) for closer inspection. Applying our criteria for the missing baryons (log T(K) > 5.5 and outside R₂₀₀, see Sect. 3.2.2) to the above sample we then selected the missing baryon population captured by the Bisous filament finder. We thus define the ratio of the total mass of this sample to that of the total missing baryons in the full EAGLE volume (P_true, see Sect. 3.2.2) as the capture fraction, that is the content of missing baryons that reside in the intergalactic space within filaments. According to EAGLE and the Bisous analysis, this capture fraction is ≈87%.

Since the filaments occupy only ≈5% of the full simulation volume (see Sect. 5), the above results on the EAGLE data demonstrate that the missing baryons are strongly concentrated towards the filaments. This is consistent with the accretion shock heating scenario (Ryu et al. 2003; Kang et al. 2005). Since the Bisous formalism is applicable to observational galaxy surveys, our tests indicate that by focusing on the galaxy filaments detected in optical data one can localise most of the missing baryons.

The fraction of the captured missing baryons is comfortably high, especially considering that our current implementation of the Bisous formalism is optimised for the SDSS galactic filament detection and not for gas detection (see Sect. 5). In particular, we are missing some of the more extended hot gas in the outskirts of the thickest filaments which is less strongly correlated with the galaxies (see Fig. 7).

6.2. Thermodynamic properties

To investigate the thermodynamic properties of the baryon population within Bisous filaments, we analysed the corresponding mass distribution dM/dT and phase diagram (similarly to the collapsed phase, see Sect. 3.2.1). While the collapsed phase dominates the higher temperatures, the gas captured by the Bisous filaments is most abundant at temperatures log T(K) = 5 − 6.5 (Fig. 11). Within this temperature range, ≈80% of the gas in filaments resides outside R₂₀₀. The distribution peaks within the missing baryon phase (described in Sect. 3.2.2), at a temperature of log T(K) ∼ 5.9. Above the missing baryon temperature limit (log T(K) = 5.5) and outside R₂₀₀, the Bisous filaments contain ≈25% of the total baryon budget (see Table 2). Since this diffuse gas is hard to observe, it fits well in the role of the missing baryons. Its detection would constitute a crucial step towards the solution of the cosmological missing baryons problem.

The temperature-density phase diagram of the Bisous filaments (Fig. 4) indicates that while the density contrast varies greatly from the lower end (log Δ_b ∼ 0) to the higher end (log Δ_b ∼ 4), most of the gas is concentrated within Δ_b ∼ 10 − 150. This density concentration corresponds to the temperature concentration at log T(K) = 5 − 6.5 reported above, rendering this phase quite compact and smooth in the temperature-density plane.

The hottest and densest phase within Bisous filaments resides within the outskirts of galaxies and groups of galaxies (see Sect. 3.2.1). The lower density and temperature phase of the Bisous baryons forms a progressive transition from the high density and temperature virial and outskirt phases (see Fig. 4), albeit with some overlap between the different phases.

The baryon sample extracted within the filaments also contains already robustly detected lower temperature phases (see Fig. 11). The ratio of the missing baryon mass to the total baryon mass within the filaments, for our purposes the purity, in the full EAGLE filament sample is ≈63%. Our experiments with the EAGLE data and Bisous filaments indicated that by excluding the shorter filaments we do not improve the purity significantly, while the capture fraction decreases. This is due to the lack of correlation between the gas temperature and the Bisous filament spine length. This indicates that the filament length (which can be derived relatively easily from the galaxy spectroscopic measurements via the Bisous filament finder, see more discussion in Sect. 5.3) alone is not useful for the selection of the most missing-baryon rich filaments.

Thus, we analysed the thermodynamic properties of the δ_LD groups described in Sect. 5.3. We found that the purity increases with δ_LD, reaching a level of ≈82% for the high density class. This is due to the increasing relative contribution from the hotter gas with higher δ_LD (see Fig. 11). This re-enforces our previous statement about the need to focus on the filaments with a high luminosity density when searching for the missing baryons (Sect. 5.3).

6.3. Extending the view with NEXUS+

When deriving observational implications for the diffuse hot cosmic baryons with our approach, we should estimate the systematic uncertainties related to our essential assumptions that (1) the Bisous filament finder method we employed does accurately capture the hot baryon population related to the cosmic filaments and (2) the data in the EAGLE simulation we utilised does form an accurate representation of the hot baryons in the local Universe. In order to isolate the two effects one should perform filament detection with several methods (including Bisous) to several simulations (including EAGLE) and investigate the hot baryon content extracted from the different filament samples. Such a substantial amount of work is beyond the scope of this paper.

To get an idea of the level of the systematics in our results due to the filament detection methodology, we performed one test on the same EAGLE data we have been using throughout this paper, but applying a very different method, NEXUS+ (Cautun et al. 2013) instead of Bisous, to dark matter density fields instead of the galaxies. NEXUS+ uses the morphology of the density field deformation tensor (Hessian) as a basis for assigning each location of the large-scale environment to a void, wall, filament, or node. On the other hand, the Bisous method uses the spatial distribution of galaxies to detect the filamentary network. An important difference between the two methods is that NEXUS+ is a scale-free algorithm, whereas Bisous requires the specification of scales on which to search for the filaments.

The definition of the boundary of a filament is essential when capturing the hot diffuse baryons embedded in the filaments, and it is very different in NEXUS+ and Bisous. The main setup of NEXUS+ includes the following steps. First, a log-Gaussian filter is applied to a given density field at several smoothing scales. Second, the Hessian matrix of this field is calculated at each smoothing scale and the eigenvalues λ₁ ≤ λ₂ ≤ λ₃ of this Hessian are obtained. In the third step each point is assigned a morphology based on the inequalities of the eigenvalues. Nodes are defined as regions with λ₁ ≈ λ₂ ≈ λ₃ < 0, filaments as regions with similar densities along two axes, λ₁ ≈ λ₂ < 0 with λ₂ ≪ λ₃, and walls are regions with λ₁ ≪ λ₂ and λ₁ < 0. Any volume element that cannot be characterised as a filament, wall, or node is classified as a void. This process is repeated for every smoothing scale and in the final step, the information is combined to obtain a scale-independent signature. Morphological features with a signature above this threshold are considered as valid structures. This scale-space analysis ensures that the morphological features prevalent on several spatial scales are identified. As a result, NEXUS+ gives a volume-filling field with each point labelled by the local morphology: node, filament, wall, and void.

On the other hand, in our current implementation of the Bisous method the location of the filament boundary essentially depends on our choice of a spatial location being assigned to a filament if it is covered by more than 5% of the visit maps, indicating the locations covered by the fitted cylinders (see the full description in Sect. 5). If we defined the filament volume by requiring a higher number of Bisous cylinders to cover it, we would shrink the filament boundaries closer to the filament axis and thus capture smaller fraction of missing baryons. We revisit this issue when deriving conclusions below.

While the two particular methods we are comparing here work with distinct mass components (NEXUS+ traces dark matter and gas; Bisous works on galaxies), they are expected to be spatially well correlated. Namely, the dark matter dominates the formation of the large-scale structure and drives the formation of galaxies embedded within the cosmic web (e.g. Cautun et al. 2013). Also, the diffuse baryons are expected to be concentrated around the dark matter density maxima of the large-scale structure, that is the axes of the cosmic web filaments. In fact, the volumes covered and the amount of dark matter, gas, and stellar mass contained by the three different filament samples were consistent within 10%. These agreements suggest that (1) the NEXUS+ and Bisous methods quite accurately detect the main filamentary network and (2) dark matter, gas and galaxies quite accurately trace the same filaments.

Our specific focus in the current paper is on the hot phase of the intergalactic gas in the filaments. Therefore, the above comparison of the full gas population is not adequate for quantitative assessment in our case. We thus utilised the above published characterisation of the large-scale environments obtained by applying the NEXUS+ method on the dark matter density fields in the EAGLE simulations (Ganeshaiah Veena et al. 2019). We repeated our analysis of the thermodynamic properties of the diffuse EAGLE baryons (excluding the gas within R₂₀₀, see Figs. 3 and 11), but this time using NEXUS+ voids, walls, filaments, and nodes. The resulting distributions of the mass as a function of temperature and overdensity indicate that through the voids-walls-filaments-nodes sequence the temperature and density increase systematically (see Fig. 12). This is expected in the scenario of gravitational collapse and shock heating, where sheets, filaments, and nodes are collapsed in 1, 2, and 3 dimensions, respectively. Each collapsed dimension increases the density and causes more shocks, which in turn induces more heating and higher temperatures. Voids, on the other hand, experience no collapse and thus remain at the low end of the distributions.

Fig. 12. Top: absolute (left panel, normalised to unity at the maximum value of the full baryon gas sample, same as Figs. 3 and 11) and relative distributions (right panel) of the gas mass as a function of the mass-weighted temperature for Bisous filaments (black continuous line), NEXUS+ filaments (orange, dashed line), NEXUS+ nodes (blue, dotted line), NEXUS+ walls (green, dash-dotted line), and NEXUS+ voids (yellow, dot-dashed-dotted line). Bottom: same as the top panel, but for the density contrast. The gas within R200 has been excluded in all panels.

We found that the filaments detected by applying NEXUS+ to the dark matter density fields of EAGLE captured ≈79% of the missing baryons (i.e. diffuse baryons with log T(K) > 5.5, outside R₂₀₀). This is somewhat smaller than the corresponding value obtained with Bisous (≈87%, see Sect. 6.1). At the lower temperature end of the missing baryon phase (log T(K) ∼ 5.5 − 6, see Fig. 12), the difference can be attributed to walls. As discussed in Ganeshaiah Veena et al. (2019), ≈30% of the Bisous filament galaxies are located in regions classified as walls by NEXUS+ in the EAGLE simulation. The temperature distribution of the EAGLE baryons (see Fig. 12) shows that the NEXUS+ walls contain a significant amount of the hot baryons. The downside is that most of the gas in walls is cooler than log T(K) = 5.5 and thus its inclusion unavoidably reduces the purity (i.e. the mass fraction of the hottest WHIM to all the IGM captured within the filaments).

At the high temperature end (log T(K) > 6.8) the relative mass distribution within the Bisous filaments has a bump, due to the outskirts of galaxies and groups of galaxies (see Fig. 12). NEXUS+ filaments do not have this feature, probably because the regions we defined as outskirts (using the FoF and R₂₀₀ information, see Sect. 3.2.1) are situated within NEXUS+ nodes, and thus the baryons in those regions are not captured within the NEXUS+ filaments. This explains partly the lower missing baryon capture fraction compared to that obtained with the Bisous filaments. The gas in the outskirts has very high purity which very closely compensates for the lower purity in the walls.

The mass distribution as a function of temperature in different NEXUS+ environments corroborates the visual inspection carried out in Sect. 4.2 (see Fig. 6). The lowest temperature range (log T(K) < 3) is dominated by voids. Between log T(K) = 3 − 4 voids are still the largest component, but with increasing contribution from walls. At log T(K) = 4 − 5 the walls dominate, while at log T(K) > 5 filaments emerge as the most prominent structures.

The density distributions of the intergalactic gas (outside R₂₀₀) in NEXUS+ and Bisous filaments are very similar (see Fig. 12). As with the temperature distribution, on the high end of the density distribution the gas density within NEXUS+ filaments drops slightly faster than the density within Bisous filaments. This is explained by the gas contained in NEXUS+ nodes becoming more prominent towards higher densities.

In summary, assuming that the Bisous and NEXUS+ methods bracket the systematic effect related to the filament detection methodology, the total mass of the hot intergalactic baryons captured in the filaments may vary by ≈10%. We consider this to be a good match, considering the very different approaches these two methods use to obtain filaments. While the purity remains almost the same, the Bisous method captures a higher amount of missing baryons resulting from the different galaxy populations belonging to the filaments traced by the two methods. Nonetheless, both methods indicate that from all the hot gas in EAGLE (≈42% at log T(K) > 5.5, see Table 2), ≈54%–59% reside in the IGM (outside R₂₀₀) within filaments.

7. Profiles

We next quantitatively analysed the visual regularities described in the previous section, whereby hotter and denser gas is closer to the filament spines, by deriving radial temperature and density profiles. For this end we used a simple, cylindrical geometry, where the gas properties were projected into a one dimensional distribution as a function of distance from the filament spine. For this aim we binned the filament gas within concentric hollow cylinders centred on each spine (for the radial binning, see Sect. 7.1). We did not apply the filament volume definition (the visit map in Sect. 5) here, as it would truncate the profiles at ∼1 Mpc (the number of particles per bin drops dramatically beyond this radius, see Sect. 5.2). Instead, we were interested in seeing how the thermodynamic properties of the gas change when moving further away from the filaments.

We co-added the baryon masses and calculated the median baryon temperatures using the gas particles located within each hollow cylinder. We subsequently divided the masses by the hollow cylinder volumes to obtain densities in each distance bin. We repeated the above procedure for all filaments, thus obtaining a sample of 779 radial density and temperature profiles. We investigated the filament-to-filament variation of the density and temperature at each radius and include the scatter in our profiles (see Sects. 7.1 and 7.2).

We also performed a null test by repeating the profile extraction but this time changing the locations and orientations of the filament spines. We used the same 779 filament spines as before, but rotated each spine in a random direction and displaced it a random distance within the simulation volume. This way we eliminated any possible co-spatial connection between the filament spines and the gas distribution.

7.1. Inner temperature profiles

The resulting temperature distributions of one value per filament in each radial bin have a well-defined maximum and a symmetric distribution around it in logarithmic space. We thus fitted the distributions of the full sample with a log-normal model and used the best-fit parameters to construct the radial profiles of the most likely values and the 68% confidence intervals of the temperatures (see Fig. 13). The resulting profile peaks at r = 0 and decreases smoothly with increasing radius.

Fig. 13. Temperature profile as a function of distance from filament spines. The pink, blue, and brown bars indicate the 1σ interval of the log-normal distribution of the median values of individual filaments within each distance bin for the high δLD, full, and random filament samples, respectively. The red and purple squares mark the peak of each respective distribution. The solid yellow line is the best-fit temperature profile model fitted to the high δLD sample. The brown dashed line indicates the background temperature, while the black dashed lines mark the temperature limits of the WHIM (log T(K) = 5) and the hot WHIM (i.e. missing baryons, log T(K) = 5.5).

The random profile sample (see above) indicated no significant radial temperature variation (see Fig. 13) demonstrating that the central temperature increase in the filament sample is indeed associated with the filament spines. The small increasing trend may be explained by the larger volumes at the outer radii more frequently capturing structures including hot gas. This also serves as an indication that the Bisous formalism we adopted for the filament detection is robust. Since the median background value (T ≈ 2 × 10³ K) is below 10% of the profile values at all radii, we ignored it in the profile fitting below.

The temperature profile derived from the full filament sample has a maximum of log T(K) ≈ 5 at r = 0. Thus, it is very unlikely to randomly encounter the hot WHIM temperatures. We therefore repeated the above procedure using only the sub-sample of the filaments with the highest luminosity overdensity values (δ_LD = 20–50, defined in Sect. 5.3).

In order to ensure proper counting statistics, we required that each hollow cylinder must contain more than 1000 gas particles. Additionally, in order to ensure proper sampling, we required that the cylinder width be at least 5 times the mean interparticle distance d given by the formula d = (100/1504)/(1+δ_b)^1/3 (in Mpc), where δ_b is the baryon overdensity computed for each particle. Considering both requirements, we increased the hollow cylinder width with increasing radius, from 0.1 Mpc in the centre to 0.4 Mpc at the largest radii.

As above, the maximum temperature occurs at r = 0, but this time at a temperature of log T(K) ≈ 6 (see Fig. 13). This shows the importance of the high δ_LD sample: we probe the hot WHIM (i.e. the missing baryons) more efficiently, when focusing the observational efforts on the filaments with the highest luminosity overdensities.

Visual inspection of the temperature profiles (see Fig. 13) indicated a quite regular radial behaviour, justifying the employment of a radial profile model representing the high δ_LD sample. After experimenting with analytical functions, we settled on the commonly used single-β model with a modification. We noticed that the model fits the data better if we change the exponent of the $\frac{r}{r_{c,T}}$ term from 2 to 3, our fitting function thus being

$$\begin{array}{c}\hfill T(r)=T_0\times {[1+{\left(\frac{r}{r_{c,T}}\right)}^3]}^{(-\frac{3}{2}{\beta }_T)}.\end{array}$$(1)

Since our profiles only extend to few times the core radius r_c, T (see below), we could not meaningfully constrain the slope parameter β_T. There is a high degree of degeneracy between r_c, T and β_T and thus, after trial and error, we fixed β_T to 1.0. We then fitted T₀ and r_c, T to the data in the radial range r = 0–4 Mpc following this procedure: (1) We used the best-fit log-normal models of the temperature distributions at each radius (obtained above) to derive random variables and thus constructed 1000 randomised temperature profiles³. Such a sample statistically reflects the filament-to-filament variation. (2) We fitted each randomised temperature profile with the above model using a simulated annealing procedure (Press et al. 2007) to minimise the residuals between the model prediction and the simulated temperature value. (3) We fitted the distributions of the best-fit T₀ and r_c, T values with a log-normal model (see Fig. 13), thus yielding the most likely values and the 68% confidence intervals.

The procedure yielded $T_0=1.2_{-0.2}^{+0.3}\times {10}^6$ K, $r_{c,T}=1.4_{-0.2}^{+0.3}$ Mpc (see Fig. 13). Our temperature profile is qualitatively similar to that of Gheller & Vazza (2019), in that it has a rather flat central part and a steeper drop at larger radii. Quantitatively the match is not good: our values are lower at all radii and our core radius is smaller.

The temperature of the best-fit model decreases (Fig. 13) below the approximate warm and hot WHIM border of 10^5.5 K at r ≈ 1.5 Mpc. This demonstrates the importance of focusing the missing baryon search close to the filament spines.

We then used the observational results obtained by Tanimura et al. (2020) on the radial distribution of the filament gas temperatures for comparison with our results derived from the EAGLE simulation. Using Planck (Planck Collaboration XXII 2016) thermal Sunyaev-Zeldovich (tSZ) and cosmic microwave background (CMB) lensing data, Tanimura et al. (2020) derived a constant temperature profile of $T=1.4_{-0.4}^{+0.4}\times {10}^6$ K up to r = 5 Mpc. While Tanimura et al. (2020) used a different filament detection method (DisPerSE) for a different sample of SDSS filaments, their results are consistent with ours in the central r ≈ 1 Mpc region (see Fig. 14). We further discuss the similarities and differences in Sect. 8.

Fig. 14. Temperature as a function of distance from the filament spines. The solid yellow line is the best-fit profile for the high δLD sample with its respective distribution indicated by the pink bars (same as in Fig. 13). The black dash-dotted line is the temperature derived by Tanimura et al. (2020) from Planck tSZ and CMB-lensing observations, with the grey area showing the error. The dashed line at log T(K) = 5.5 indicates the lower limit for the hot WHIM.

7.2. Inner density profiles

Similarly to the temperatures, the densities in each radial bin follow a log-normal distribution well. Following the above procedure, we determined the profiles of the most likely density and the 68% confidence intervals for the full sample (see Fig. 15). The resulting profile has a maximum at r = 0 and decreases smoothly with the increasing radius.

Fig. 15. Gas density profile as a function of distance from filament spines. The 1σ interval of the log-normal distribution of individual filament densities at each distance bin for the high δLD, full, and random filament samples is represented by pink, blue, and brown bars, respectively, with the squares indicating the peak of the distribution. The yellow line is the best-fit density profile model of the high δLD sample. The descending brown dotted line is the best-fit density profile model without the effect of the background. The black dashed line indicates the mean cosmic baryon density, while the brown dotted line shows the background level.

As in the case of the temperatures above, the null test (i.e. the random profile sample) indicates no significant radial density variation (see Fig. 15), reinforcing the conclusions about the association of the physical WHIM properties with the filaments and the Bisous robustness. We adopt the mean value ρ ≈ 6.5 × 10⁸ M_⊙ Mpc⁻³ as the background level for the profile fitting below.

The maximum baryon overdensity derived from the full filament sample is very low, δ_b≈ 5 at r = 0. Thus, considering the full filament sample, it is very unlikely to randomly encounter the highest expected baryon overdensities in the cosmic filaments (∼100). We therefore repeated the above procedure using only the high δ_LD, defined in Sect. 5.3. As above, the baryon overdensity is largest at r = 0, but this time at a much higher value, δ_b ≈ 40 (see Fig. 15). Thus, by focusing observational efforts on the filaments with the highest luminosity densities one probes the densest phase of the hot WHIM, namely the one with the highest X-ray and SZ signal.

Visual inspection of the density profiles indicated a quite regular radial behaviour, justifying the employment of a radial density profile model representing the high δ_LD sample. We fitted the density profile using a single-β model, usually used for the intracluster gas density:

$$\begin{array}{c}\hfill \rho (r)={\rho }_0\times {[1+{\left(\frac{r}{r_{c,\rho }}\right)}^2]}^{(-\frac{3}{2}{\beta }_{\rho })}+\mathrm{bkg}\end{array}$$(2)

using the same randomisation scheme as above in the case of the temperature profiles. Following the arguments given for the case of the temperature profiles, we fixed the slope parameter β_ρ, this time to 2.0. Since the density of the filament sample approaches the background level at the largest radii, we included the constant background level (bkg) derived above in our density profile model.

The procedure yielded ${\rho }_0=2.3_{-0.5}^{+0.6}\times {10}^{11}$ M_⊙ Mpc⁻³ and $r_{c,\rho }=1.2_{-0.2}^{+0.3}$ Mpc. Our gas mass density profile is different from the gas particle density profile in the Gheller & Vazza (2019) simulations. Our profile peaks at higher values at r = 0 and decreases to underdense values at r > 2 Mpc, while the Gheller & Vazza (2019) profile approaches a constant overdensity level of approximately 3 at r = 2 − 3 Mpc. In any case, this further demonstrates the importance of focusing the missing baryon search close to the filament spines.

As in the case of the temperatures above, we compared our results from the EAGLE simulation with the observational results from Tanimura et al. (2020), but this time for the radial distribution of the density of the filament gas. Again, within the central r ≈ 1 Mpc both the simulations and the observations agree (see Fig. 16). This lends confidence on the results and strengthens the simulation-based insight that it is beneficial to focus the missing baryon search to the central r ≈ 1 Mpc regions around the filament axes. For further discussion see Sect. 8.

Fig. 16. Comparison between the baryon density profile of our high δLD sample (yellow line and pink bars) and the baryon density profile derived by Tanimura et al. (2020) (black dash-dotted line for the best-fit model and grey area for the error). In order to derive their density profile, Tanimura et al. (2020) fitted a beta-model profile to the tSZ and CMB-lensing data observed by Planck. The mean cosmic baryon density is represented by the black dashed line.

7.3. Outside the filament core regions

The diffuse baryon mass density profile of the high δ_LD filament sample is significantly above the background level in the inner r = 1 Mpc region. The background density level, obtained with the randomised profiles, should be dominated by the voids due to their large volume filling fraction of gas (76%) in EAGLE, based on a NEXUS+ analysis (Ganeshaiah Veena et al. 2019). Indeed, the background level is ∼10% of the cosmic mean, consistent with the density distribution in the NEXUS+ voids (see Fig. 12). This indicates that on average at r ≈ 3 Mpc distance from a filament spine we enter the void domain.

Similarly, the background level of the temperatures, log T(K) ∼3.0–3.5, is consistent with the temperature distribution in the NEXUS+ voids (see Fig. 12). However, the filament temperature profile indicates higher values of log T(K) ≈ 4.0–4.5 at r ≥ 3 Mpc (see Fig. 13). Analysis of the temperature distribution of the particles in the full filament sample at these radii indicated that in addition to the maximum at the expected void value of log T(K) ≈ 3, ≈20% of the particles have extremely high temperatures (in excess of log T(K) = 5) which drive the medians of the temperature distributions of individual filaments to high values. As the full filament distribution is calculated using a log-normal fit of the individual medians, the whole distribution is thus biased towards high temperatures. The imaging analysis indicated that this population of gas resides in extended hot regions around massive haloes, probably due to the most energetic shocks propagating to larger distances. Their importance is negligible in the randomised profile sample due to the very small volume filling fraction of the haloes (Ganeshaiah Veena et al. 2019). The consequent ≈20% excess in the particle density at r ≈ 3 Mpc above the approximate void level (≈10% of the cosmic mean) is negligible compared to the scatter in the density profile. In summary, combined Bisous and NEXUS+ analysis indicates that at r ≈ 3 Mpc distance from a filament we enter the void domain, unless we are in a neighbourhood of a massive halo.

8. Discussion

8.1. Limitations of EAGLE

We showed that hot and diffuse intergalactic gas within cosmic filaments makes up a significant fraction of the baryons in the EAGLE simulation at redshift z ∼ 0. However, several EAGLE related limitations could potentially have significant consequences on our results.

The volume covered by the simulation is relatively small (100³ Mpc³), suppressing the possible larger structure formation modes. Thus we were not able to study the potential variation of the cosmic baryon heating related to the longest filaments (> 50 Mpc) missing in EAGLE. However, a recent filament study on large volume simulations (Illustris TNG300, 300³ Mpc³, and Magneticum, 500³ Mpc³, Galárraga-Espinosa et al. 2020) has shown that the number of filaments decreases rapidly as a function of filament length, with much less than 1% of filaments being larger than 50 Mpc. Ignoring such a small sample will not significantly affect our results.

On the other hand, the resolution of the simulation limits the minimum scale of structures. The mass resolution in EAGLE is enough to resolve the Jeans scales for the warm interstellar medium (of the order of kiloparsecs, Schaye et al. 2015), whereas the filamentary structure and the shocks heating the IGM are part of the large-scale structure (of the order of megaparsecs, see e.g. Ryu et al. 2003). Thus, we do not expect the resolution limits to affect our results.

In addition, simulations may fail to reproduce all the relevant thermodynamic interactions that take place during structure formation. While the IGM in filaments is presumably heated mainly via accretion shocks driven by gravity (e.g. Ryu et al. 2003), the winds driven by feedback processes from star formation and AGN that are included in EAGLE are also thought to affect the WHIM (e.g. Tepper-García et al. 2012). However, these galactic winds are poorly constrained observationally. Cosmic rays, which are not modelled in EAGLE, may or may not have a significant impact on the heating process (Hopkins et al. 2020). Nonetheless, while cosmic rays could potentially change the thermal structure of the CGM around the galaxies they are expelled from, the extent of these rays is limited. Indeed, cross-correlation of the extragalactic gamma-ray background with SZ maps imply that cosmic rays cannot affect for a significant fraction of the total pressure in clusters of galaxies (Shirasaki et al. 2020). Thus, cosmic rays are not expected to alter the intergalactic gas studied here.

Another factor that might have an effect on the results presented here is the cooling of the IGM. While the subgrid radiative cooling is implemented in EAGLE in a sophisticated, element-by-element manner, uncertainties in the intergalactic metallicities may lead to an under- or overestimate of the cooling rates. If, in reality, the metal abundances in filaments are much higher than simulated in EAGLE, the cooling rate of the IGM would be higher. This could lead to an overestimate of the missing baryon population in filaments based on EAGLE. However, the intergalactic metallicity in filaments is not well known. Adopting the common assumption of 0.1 times the Solar metallicity, filamentary gas in the hot WHIM temperature (log T(K) > 5.5) and density range (δ ≈ 1 − 50, see Fig. 15) has a cooling time longer than the Hubble time (Wiersma et al. 2009). Assuming the extreme case of Solar metallicities in the intergalactic gas, the cooling time scale is still above or comparable to the Hubble time. Thus, any filamentary IGM that reaches these high temperatures is expected to remain hot.

Further reassurance is given by the good agreement of our temperature and density profiles with recent observations of the tSZ effect in filaments (Tanimura et al. 2020, see Figs. 14 and 16). This agreement would be highly unlikely if the EAGLE filaments were unrealistic. Thus, albeit future simulations may implement improvements to account for all the limitations mentioned above, we do not expect significant changes in our results.

8.2. Bisous as a missing baryon finder

We found that with our current version of the Bisous formalism, we can potentially locate a large fraction of the missing baryons by tracing them with the cosmic filaments detected using the observational spectroscopic galaxy catalogues. Based on the analysis of the EAGLE data it is beneficial to concentrate the observational missing baryon search (via e.g. the SZ effect, X-ray absorption and emission) within ∼1 Mpc distance from the axes of the cosmic filaments with the highest (observable) galaxy luminosity densities, since in these regions the diffuse baryon density exceeds ten times the cosmic mean and the baryon temperature is significantly above the approximate hot WHIM threshold of log T(K) = 5.5 ( see Figs. 13 and 15). On the other hand, the average Bisous filament radius (based on galaxy distributions only, i.e. independent of the gas, see Sect. 5.2) in our current implementation of the Bisous method is also ∼1 Mpc. Thus, there is not much room for relaxing the Bisous filament volume boundary condition (i.e. lowering the visit map value from 0.05, see Sect. 6.3) since it would add very little hot WHIM while increasing the amount of gas with low densities and temperatures. On the other hand, a stricter visit map criterion would result in using smaller extraction radii and consequently in less hot WHIM being covered. Thus, the current Bisous setting is close to optimal for the missing baryon search.

8.3. Simulations versus observations

We modelled and reported the spatial distribution of the basic thermodynamic quantities (density and temperature) of the hot diffuse gas in the filamentary environments as indicated by the EAGLE simulations. This may serve as a starting point for producing predictions for observational signals to be used for planning future observations and interpreting the measured signals. We will carry out such work in the near future for the optimal SDSS filament sample identified in this work (see Sect. 5.3).

8.3.1. The SZ effect

Perhaps the most promising observational avenue for detecting the missing baryons is to utilise the SZ effect imposed on the cosmic microwave background by the hot gas in the filaments (e.g. de Graaff et al. 2019; Tanimura et al. 2019, 2020). The data from the Planck mission (Planck Collaboration XXII 2016) currently has the best combination of the sky coverage and sensitivity for this work. Stacking the Planck data around the axes of the cosmic filaments for a large enough sample may enable a statistically significant signal to be measured. Among the current spectroscopic galaxy surveys, the Sloan Digital Sky Survey has the optimal combination of sky coverage and sensitivity for producing a large data base for detecting the cosmic filaments.

The stacking of Planck tSZ and CMB-lensing data on a different sample of SDSS filaments than ours was carried out by Tanimura et al. (2020). While Tanimura et al. (2020) used filaments constructed with the DisPerSE method and a length criterion of l = 30 − 100 Mpc, we used the Bisous method which yielded a maximum filament length of 35 Mpc. This may lead to a comparison between two distinct filament populations, as Galárraga-Espinosa et al. (2020) showed that in simulations longer (l > 20 Mpc) DisPerSE filaments are thinner and have lower galaxy densities. Thus, the Bisous high-luminosity density sample might correspond to high galaxy density, in other words to shorter and wider filaments. We keep this caveat in mind when comparing the two works below.

Within the central r ≈ 1 Mpc, the EAGLE temperatures and baryon densities are consistent with those derived from the Planck data. This raises the question whether this is compatible with the fact that the two filament samples are chosen differently. Perhaps the central galaxy density difference between the short and long filaments (Galárraga-Espinosa et al. 2020) is relatively small, and consequently the difference in gas properties is within the larger error bars coming from the Planck and EAGLE analysis. In that case, the consistency between EAGLE and Planck is not a co-incidence but rather lends support to the results.

However, at r  ≳  1 Mpc there are clear differences: the EAGLE density and temperature profiles are steeper. In particular, EAGLE does not support the isothermality assumption beyond r ≈ 1 Mpc, while the isothermal profile derived by Tanimura et al. (2020) extends up to r = 5 Mpc. If one was to use a radially decreasing temperature profile in the Planck analysis, the density profile should become even flatter to yield the same SZ signal. Possibly, the difference in the samples of filaments may explain the discrepancy at larger radii. At larger radii the observed SZ signal is dominated by the point spread function scatter of the signal originating from the central r ≈ 1 Mpc in the filaments. Also, the total signal at larger radii is background-dominated. Thus, the shape of the SZ profile may not be very accurately determined observationally. The possible bias in the temperature and density profiles at larger radii may remain unnoticed due to their very small statistical weight on the fitting procedure.

8.3.2. X-rays

With ATHENA it may become possible to measure a significantly large sample of X-ray absorption signals from high ions (e.g. O VII) embedded in the WHIM within filaments in the spectra of blazars behind them. Analyses of the EAGLE and IllustrisTNG simulations (Wijers et al. 2019, 2020; Nelson et al. 2017) have shown that the hot gas tracers O VII and O VIII appear densely packed around the galaxies, while the lower density phase follows the large-scale structure. A recent work on SIMBA simulations (Borrow et al. 2019) found that AGN jets may carry 10% of the baryons up to 4 Mpc distances from the host galaxies and thus enrich the intergalactic medium with metal ions. We will extend the above works and pave the way for ATHENA by examining and modelling the spatial distributions of the high ion densities in the EAGLE simulation within the cosmic filaments in a future work.

Once we have mapped the metals within the filaments in the EAGLE simulation, we can produce predictions for the soft X-ray emission. The currently operational eROSITA satellite may have the potential to detect the relatively faint X-ray emission form a large number of stacked 4MOST filaments.

9. Conclusions

In this work we presented a novel method for localising the missing baryons around cosmic filaments and estimating their spatio-thermal properties. For developing our method, we chose to use the state-of-the-art hydrodynamical EAGLE simulation, as it provides information on both galaxies and baryon gas. Here we summarise our results:

• We identified the filamentary structure within the EAGLE simulation using the previously published NEXUS+ classification of the large-scale structure environment based on the dark matter density field and by applying the Bisous model to an SDSS-like magnitude limited galaxy sample.

• The diffuse hot intergalactic medium (log T(K) > 5.5, outside R₂₀₀) captured by filaments amounts to ≈23%−25% of the total baryon budget, or ≈79%−87% of all the hot WHIM. This accounts for a large fraction of the missing baryons.

• Filaments fill only ≈5% of the total volume, indicating that the hot WHIM is very tightly concentrated along the filament spines. Thus our tests indicate that by focusing on the galaxy filaments detected in optical data, one can localise most of the missing baryons.

• The hot WHIM mass in NEXUS+ and Bisous filaments differs by only ∼10%. This indicates that the missing baryons can be reliably traced with these filament finding methods.

• Based on the galaxy luminosity density we selected the optimal ∼10% of the filament spines for more detailed study. The purity (the mass fraction of the missing baryons) of this high δ_LD sample is ≈82%.

• For our high δ_LD sample we modelled the density and temperature profiles as a function of distance from the filament spine with analytic functions. The profiles show that in order to obtain the strongest observational signal, it is beneficial to focus on the central ∼1 Mpc from the spine.

• Our temperature and density profiles of the filament gas are consistent with those derived from Planck tSZ and CMB lensing observations (Tanimura et al. 2020) within the central ∼1 Mpc regions of the filaments.

These results allow for some optimism in locating the hot WHIM. The filament detection methods together with the baryon density and temperature profiles derived here are applicable to observations and may thus be useful for planning future research on the missing baryons.

Acknowledgments

We acknowledge the support by the Estonian Research Council grants PUT246, IUT40-2, PRG1006, and by the European Regional Development Fund (TK133). We thank Stuart McAlpine and Till Sawala for help with the EAGLE data. Thanks to Davide Martizzi for his help on IllustrisTNG. We acknowledge the Virgo Consortium for making their simulation data available. The eagle simulations were performed using the DiRAC-2 facility at Durham, managed by the ICC, and the PRACE facility Curie based in France at TGCC, CEA, Bruyèresle-Châtel. We thank Rien van de Weygaert, Bernard J. T. Jones and Marius Cautun for their permission to use data of the NEXUS+ analysis. We also thank Nabila Aghanim and Hideki Tanimura for providing the density and temperature profiles derived from observations.

References

All Tables

All Figures

	Fig. 1. Projected diffuse baryon density within a representative 5 Mpc slice (colour map) above the mean baryon density and galaxies brighter than Mr = −18.4 (blue dots) in the EAGLE simulation. The light blue colour denotes the regions below the mean baryon density. The density was computed by co-adding the masses of individual gas particles and dividing by the projected volume in each pixel.
In the text

	Fig. 2. Luminosity distributions of EAGLE (yellow points and error bars) and SDSS (red points and error bars) galaxies. The vertical dashed lines indicate the selected magnitude cuts of −18.4 and −19 for EAGLE and SDSS galaxies, respectively.
In the text

Fig. 3. Left panel: distribution of the full EAGLE baryon gas mass content as a function of the temperature, normalised to unity at log T(K) = 3. Right panel: corresponding distributions of different components, relative to the full distribution: the bound phase (orange dashed line), the virial phase within R200 (blue dotted line), and the gas on the outskirts of galaxies and groups of galaxies (outside R200 but within FoF halo, green dash-dotted line). The black vertical line indicates the approximate lower limit for the hot WHIM (i.e. log T(K) = 5.5).

In the text

Fig. 4. Baryon density contrast – temperature phase diagrams (for the convenience of plotting we use the density contrast Δb = 1 + δb rather than the overdensity δb). All the gas in the full EAGLE volume is shown in the top left panel, the gas within R200 in the top right panel and the gas within the Bisous filaments excluding gas within R200 in the bottom left panel. For each temperature and density bin we co-added the mass of particles and divided by the total mass of the full box. We cut the colour scale at 1% level of the maximum value for clarity. The bottom right panel indicates the isomass curves containing 68% of the baryon mass in Bisous filaments excluding R200 (blue line), within the outskirts (black line), and within R200 (red line). The low-density, low-temperature track seen in the upper left panel corresponds to the photo-heated IGM.

In the text

	Fig. 5. Mass-weighted mean temperature distribution (in shades of blue) of the intergalactic gas (i.e. gas outside R200) within the same slice as Fig. 1. Red dots denote the location of galaxies brighter than Mr = −18.4.
In the text

Fig. 6. Projected diffuse baryon mass density contrast within the same 5 Mpc slice as Fig. 1, but here in different temperature ranges. Top left: log T(K) < 3. Top right: log T(K) = 3 − 4. Middle left: log T(K) = 4 − 5. Middle right: log T(K) = 5 − 5.5 (the warm WHIM phase). Bottom left: log T(K) > 5.5 (the missing baryon phase). Bottom right: Gas within the virial radii R200 of haloes, excluded from the other panels.

In the text

Fig. 7. Same slice as in Fig. 1, but for a visualisation of the Bisous formalism. Top left: the galaxies (blue dots) and the visit map ≥0.05 (grey). Top right: the spatial distribution of the missing baryons (colour map) and the filament volumes (grey areas). Bottom left: the spatial distribution of the missing baryons captured with the Bisous filament finder. Bottom right: the spatial distribution of the missing baryons (colour map), the filament spines of our high luminosity overdensity δLD sample (blue), and filaments excluded from the high δLD sample (magenta).

In the text

	Fig. 8. Fraction of grid points within Bisous volumes (visit map value > 0.05) as a function of distance from the filament spines. All spines were stacked together.
In the text

	Fig. 9. Galaxy luminosity density contrast (ΔLD) in the same slice as in Fig. 1, using a 2 Mpc smoothing kernel. Filament spines of the high δLD sample are shown in blue, the excluded filaments are shown in red.
In the text

	Fig. 10. Distribution of filaments as a function of the luminosity overdensity in SDSS (green line) and in EAGLE (purple line, normalised to SDSS volume). The vertical lines indicate the limits for the low (δLD = 0 − 10), medium (δLD = 10 − 20) and high (δLD = 20 − 50) luminosity overdensity groups (see Sect. 5.3).
In the text

Fig. 11. Top: absolute (left panel, normalised to unity at the maximum value of the full baryon sample, left panel in Fig. 3) and relative distributions (right panel) of the baryon gas mass distribution as a function of mass-weighted temperature in filaments with low (orange dashed line), medium (blue dotted line), and high (green dash-dotted line) luminosity overdensity. For comparison, mass within all Bisous filaments is shown with a black line. Bottom: same as top panel, but for the gas mass distribution as a function of gas density contrast. Gas within R200 has been excluded in all panels.

In the text

Fig. 12. Top: absolute (left panel, normalised to unity at the maximum value of the full baryon gas sample, same as Figs. 3 and 11) and relative distributions (right panel) of the gas mass as a function of the mass-weighted temperature for Bisous filaments (black continuous line), NEXUS+ filaments (orange, dashed line), NEXUS+ nodes (blue, dotted line), NEXUS+ walls (green, dash-dotted line), and NEXUS+ voids (yellow, dot-dashed-dotted line). Bottom: same as the top panel, but for the density contrast. The gas within R200 has been excluded in all panels.

In the text

Fig. 13. Temperature profile as a function of distance from filament spines. The pink, blue, and brown bars indicate the 1σ interval of the log-normal distribution of the median values of individual filaments within each distance bin for the high δLD, full, and random filament samples, respectively. The red and purple squares mark the peak of each respective distribution. The solid yellow line is the best-fit temperature profile model fitted to the high δLD sample. The brown dashed line indicates the background temperature, while the black dashed lines mark the temperature limits of the WHIM (log T(K) = 5) and the hot WHIM (i.e. missing baryons, log T(K) = 5.5).

In the text

Fig. 14. Temperature as a function of distance from the filament spines. The solid yellow line is the best-fit profile for the high δLD sample with its respective distribution indicated by the pink bars (same as in Fig. 13). The black dash-dotted line is the temperature derived by Tanimura et al. (2020) from Planck tSZ and CMB-lensing observations, with the grey area showing the error. The dashed line at log T(K) = 5.5 indicates the lower limit for the hot WHIM.

In the text

Fig. 15. Gas density profile as a function of distance from filament spines. The 1σ interval of the log-normal distribution of individual filament densities at each distance bin for the high δLD, full, and random filament samples is represented by pink, blue, and brown bars, respectively, with the squares indicating the peak of the distribution. The yellow line is the best-fit density profile model of the high δLD sample. The descending brown dotted line is the best-fit density profile model without the effect of the background. The black dashed line indicates the mean cosmic baryon density, while the brown dotted line shows the background level.

In the text

Fig. 16. Comparison between the baryon density profile of our high δLD sample (yellow line and pink bars) and the baryon density profile derived by Tanimura et al. (2020) (black dash-dotted line for the best-fit model and grey area for the error). In order to derive their density profile, Tanimura et al. (2020) fitted a beta-model profile to the tSZ and CMB-lensing data observed by Planck. The mean cosmic baryon density is represented by the black dashed line.

In the text