© The Authors 2024

1 Introduction

Clusters of galaxies lie at the intersection of the filamentary structure that forms the cosmic web (Bond et al. 1996). While cluster properties have also been intensively studied in a broader cosmological context, the filamentary structure, that is, both the large-scale cosmic filaments and bridges connecting pairs of clusters, is still poorly understood. The complex evolution of the cosmic web can be studied by describing the properties of the large-scale filaments and shorter bridges between galaxy clusters.

In the past few decades, the search for large-scale filaments has taken a big leap forwards with a compilation of a number of filament catalogues (usually based on the galaxy distribution), for example, within the Sloan Digital Sky Survey (SDSS; Carrón Duque et al. 2022), the two-degree Field Galaxy Redshift Survey (2dF GRS; Pimbblet et al. 2004), the Galaxy and Mass Assembly survey (GAMA; Alpaslan et al. 2013), in the Cosmic Evolution Survey (COSMOS) field (Luber et al. 2019), the Sydney-AAO Multi-Object IFS (SAMI) Galaxy Survey (Welker et al. 2020), and the VIMOS Public Extragalactic Redshift Survey (VIPERS; Malavasi et al. 2017). As for the gas component tracing the cosmic filaments, it is currently difficult to detect due to the low density and low emissivity that diffuse baryons have within filaments. The gas component in the filamentary structures, such as intercluster filaments and larger cosmic filaments, is studied by combining different data and using different techniques. For example, ROSAT X-ray data (Tanimura et al. 2020b) and stacked Compton-y maps from the Planck satellite (Tanimura et al. 2020a) were used to obtain gas temperatures and over-densities for a large number of cosmic filaments at intermediate redshifts (detected previously using SDSS data) having lengths above 30 Mpc. A study of the network of filaments around the Coma cluster was performed by Malavasi et al. (2020), and it combined SDSS data with the cosmic web detection algorithm DisPerSe (Sousbie et al. 2011) and obtained information about cluster connectivity as well a tentative detection of a signal of Sunyaev-Zeldovich (SZ) within filaments using Planck data. To analyse the signal of hot gas in extended structures, Lokken et al. (2022) combined Compton-y maps from the Atacama Cosmology Telescope (ACT) stacked on redMaPPer cluster positions from the optical Dark Energy Survey (DES).

Gas in bridges between clusters is commonly studied in X-rays (e.g. Sugawara et al. 2017; Alvarez et al. 2018, among others) combining data from Chandra, XMM-Newton, Suzaku and recently from eRosita. Some bridges (between Abell 399 and Abell 401) were also studied using thermal SZ effect and performing a multiwavelength analysis (e.g. Bonjean et al. 2018; Hincks et al. 2022, for example).

According to numerical simulations focusing on the inter-galactic gas at present epoch, ≃40–50% of the baryons hosted within filaments are in the form of a warm-hot intergalactic medium (WHIM) (Cen & Ostriker 2006; Davé et al. 2001). In the last decade, a number of observations have focused on the detection and characterisation of the diffuse WHIM in the vicinity of clusters (e.g. Takei et al. 2007; Akamatsu et al. 2011; Eckert et al. 2015; Bulbul et al. 2016) and between pairs of clusters or cluster systems, for example, A222/223 (Werner et al. 2008); A3391/A3395 (Reiprich et al. 2021 ; Veronica et al. 2024); A2029/A2033 (Mirakhor et al. 2022); Abell 98 triple merging system (Alvarez et al. 2022); and A399/401 (Akamatsu et al. 2017; Bonjean et al. 2018). Recently, Reiprich et al. (2021) analysed the complex galaxy cluster system Abell 3391/95 using SRG/eROSITA data and found that it contains (among other structures) a warm-hot emission filament 15 Mpc long. Zhang et al. (2024) used data from the eRASS All-Sky X-ray survey combined with the SDSS optical filament catalogue, detecting WHIM by stacking around 8000 X-ray filaments. Mild positive detections of the WHIM in cosmic filaments were reported using the X-ray stacking analysis of cosmic filaments (Tanimura et al. 2020b, 2022) and the stacked SZ effect signal from the intergalactic medium (de Graaff et al. 2019; Tanimura et al. 2020a).

In simulations, the detection of filaments can be done with different methods. Usually, four main components are detected as parts of the cosmic web: knots, filaments, sheets, and voids. In the past years, a number of different methods have been employed to trace the large-scale structure and its components. Some of the methods are able to detect all of the objects of the cosmic web, for example, NEXUS+ (Cautun et al. 2013), T-web (Forero-Romero et al. 2009), and V-web (Hoffman et al. 2012). Similarly, DisPerSE (Sousbie et al. 2011) and Spineweb (Aragón-Calvo et al. 2010) detect all of the cosmic web elements, but instead of knots, they map topological nodes. Some of the methods are able to obtain all the elements except knots (e.g. MMF-2 (Aragon-Calvo & Yang 2014)), while others focus only on filaments (e.g. SCMS (Chen et al. 2015a), MST (Alpaslan et al. 2013), and Bisous (Tempel et al. 2014)). Despite the different approaches, all of the methods give results that are comparable overall and have very similar structures detected. (For a detailed comparison and a more in-depth analysis of the main methods, see Libeskind et al. 2018.)

With the large number of available methods to detect the cosmic web, cosmological simulations are currently used to characterise the properties of filaments (e.g. Cautun et al. 2014; Galárraga-Espinosa et al. 2020, 2021), their evolution (e.g. Zhu et al. 2021; Galárraga-Espinosa et al. 2024), the properties of galaxies within filaments (e.g. Lee et al. 2021; Zakharova et al. 2023) and the connection between the cluster outskirts and filaments (e.g. Rost et al. 2021; Kuchner et al. 2021). Angelinelli et al. (2021) studied clumps in and around galaxy clusters in non-radiative simulations. Since high-density clumps are easier to detect with X-ray telescopes, the authors expected them to be tracers of the filaments that compose the cosmic web. They found the filament temperature to (mildly) correlate with the mass of the main cluster. Besides studies of the properties of the cosmic web, simulations are also used to investigate the origin of the gas in the observed filamentary regions (e.g. Biffi et al. 2022) and to predict future observational strategies (e.g. Churazov et al. 2023, for the detection of WHIM in the soft X-ray band).

In this work, we used a sample of re-simulated volumes extracted from a larger cosmological simulation to investigate the intrinsic properties of the filaments within those regions. By selecting regions with different accretion histories and comparing simulations where different feedback effects are at work, we aim to study the mechanisms that impact the thermal and chemical properties of the intergalactic gas within the detected filaments, focusing on the WHIM phase.

In Sect. 2, we describe the simulations used for this work. Sec. 3 is devoted to the detailed description of the methods used to detect filaments with the subspace constrained mean shift (SCMS) algorithm and how to determine their properties. In Sec. 4, we explore the physical properties of the filaments and their evolution in time. The last Sec. 5 provides a summary of our conclusions.

2 Data

2.1 The simulation code

In this work, we focus on a subset of Lagrangian regions that are centred around massive galaxy clusters and are part of the Dianoga set of simulations. These zoomed-in regions were extracted from a parent dark matter-only cosmological simulation (of h⁻³ 1 Gpc³) originally described in Bonafede et al. (2011) and later re-simulated with an improved resolution and with the addition of the baryonic component. The set of simulated regions adopted for this work was performed with GADGET-3, which is a modern version of the Tree-PM Gadget code (Springel et al. 2005) and accounts for a new description of the smooth particle hydrodynamics that includes higher order interpolation kernels (Wendland C⁴ kernel with 200 neighbours) and advanced formulations for artificial viscosity and thermal diffusion (Beck et al. 2016). The full set of simulations was described initially in Rasia et al. (2015) and extensively studied in a number of later works (see, e.g. Planelles et al. 2017; Biffi et al. 2017, 2018; Truong et al. 2018).

Simulations presented here are based on a ΛCDM model, with cosmological parameters consistent with 7-year WMAP measurements (see Komatsu et al. 2011): Ω_m = 0.24 and Ω_b = 0.04 for the density parameters of matter and baryons, H₀ = h₀ 100 km s⁻¹ Mpc⁻¹ and h₀ = 0.72 for the present day Hubble parameter, n_s = 0.96 for the primordial spectral index and σ₈ = 0.8 for the amplitude of the power spectrum of the density fluctuations. The selected Lagrangian regions are extracted and re-simulated with the zoomed-initial technique (ZIC) described in Tormen et al. (1997) by increasing the resolution in mass. In the high-resolution region, gravity is calculated with a Plummer-equivalent softening length of ϵ = 2 h⁻¹ kpc for stars and black holes and ϵ = 3.75 h⁻¹ kpc for DM and gas particles. The softening is fixed to comoving coordinates for all except DM particles, to which below z = 2 is given in physical units. The highest mass resolution of DM particles in this set of simulations is m_DM = 8.47 × 10⁸ h⁻¹ M_⊙, while the initial mass of a gas particle is m_gas = 1.53 × 10⁸ h⁻¹ M_⊙.

The set of Lagrangian regions was re-simulated with two different baryon physics models. The main difference between the two is the presence (or absence) of feedback from Active Galactic Nuclei (AGN). In particular:

CSF (Cooling and Star Formation) model. This model accounts for radiative cooling and subsequent formation of star particles. Star particles that describe the evolution of a stellar population include in the subgrid model also the feedback from supernovae (SN) and the metal enrichment from different stars (SN Ia, SN II and asymptotic giant branch AGB stars) (see, e.g. Biffi et al. 2017, 2018, for a detailed explanation). Stellar evolution and metal enrichment models are described in Tornatore et al. (2007) and allow one to follow the distribution of 15 chemical species (H, He, C, Ca, O, N, Ne, Mg, S, Si, Fe, Na, Al, Ar and Ni). These elements also contribute to the cooling. Models by Wiersma et al. (2009) and Haardt & Madau (2001) account for rates of metal-dependent radiative cooling and for the effects of the UV/X-ray background emission. Prescriptions for star formation are based on the work by Springel & Hernquist (2003), where the velocity of galactic winds v_w originated by SN driven outflows, which are fixed at 350 km/s.

AGN (Active Galactic Nuclei). These simulations rely on the same prescriptions described in the CSF case, and, in addition, they include the feedback effect from AGNs. The AGN feedback is based on a subgrid model of accretion on supermassive black holes (SMBH), where the mechanical and radiative outflows are accounted for as thermal feedback in both cases. Kinetic feedback from jets would be required to resolve subkiloparsec scales and is not modelled explicitly (see Steinborn et al. 2015). The gas accretion rate (Eddington-limited) and the SMBH mass are both included in the efficiencies of the two outflows and provide a smooth transition between the radio and quasar mode. The model can account for hot and cold accretion, but in Dianoga simulations (see Rasia et al. 2015), only cold gas accretion is considered. We refer to Steinborn et al. (2015) for a detailed description of the model and its performance in general.

2.2 The cluster sample

Table 1 contains the main characteristics of the high-density regions selected for the analysis. Three of the regions (namely D1, D6 and D22) contain a central cluster with mass M₂₀₀ > 1.7 · 10¹⁵ M_⊙ h⁻¹, while the main clusters in the other two regions (D5 and D9) have central clusters with approximately ten times smaller masses. The two regions D5 and D9 were selected as ‘isolated regions’, and their present-day volume is three times smaller than that of the largest and denser regions. We also checked the merging histories of the central clusters. In all but one case, the last major merger¹ happened at a lookback time of about 6.5 Gyr or greater (corresponding to z > 0.7). The only exception is region D6, where the last major merger happened 3.77 Gyrs ago (at z ≃ 0.3).

Inside each Lagrangian region, the code identified the main halos using a Friend-of-Friend (FoF) algorithm, while substructures within haloes were detected with the SUBFIND algorithm (Springel et al. 2001). The identification of self-bound substructures inside FoF halos was defined in its original form for DM-only simulations and extended to simulations with baryon physics by Dolag et al. (2009). The sample of main haloes was obtained with a FoF algorithm with a linking length equal to 0.16 times the mean dark matter particle separation. The centre of each halo was defined by the most bound particle (the particle with the lowest gravitational potential). Radii and masses at different overdensities were then calculated around each cluster centre, while substructures were identified within the cluster virial radius. Filaments were detected within each region using the centres of substructures.

3 Determining the filaments and their properties

3.1 The subspace constrained mean shift algorithm

The SCMS algorithm, described by Ozertem & Erdogmus (2011), was modified and used as a method for filament detection by Chen et al. (2015a). This method models filaments as ridges of the galaxy probability density function. The same algorithm has been successfully used for addressing different problems: reconstructing real cosmic filaments using SDSS (Carrón Duque et al. 2022), studying galaxy-filament alignment (Chen et al. 2019), filaments in the hydrodynamic simulation MassiveBlack-II² (Chen et al. 2015b), weak lensing maps (Moews et al. 2021) and studying velocity around stellar filaments (Chen et al. 2020).

SCMS is a multiple-step algorithm that determines filaments based on the density function p(x), where x is the spatial coordinate. The density function can be calculated with the standard kernel density estimator (KDE): $$p(x)=\frac{1}{n{h}^d}\sum _{i=1}^{n}K\left(\frac{\parallel x-X_i\parallel }{h}\right),$$(1)

where n is the number of tracers with coordinates X_i and ||x − X_i|| is the Euclidean distance between the i tracer and the location where KDE is evaluated; d is number of spatial dimensions (e.g. d = 2 or d = 3), h is the smoothing bandwidth, K is the smoothing (e.g. Gaussian) kernel. The details of the algorithm can be found in Appendix A and in Chen et al. (2015a), while in this section we specify the input values that we chose for this work:

• Tracer coordinates: As previously said, we aim at identifying the filaments from the galaxy probability function to follow an observational-like approach; thus, our tracers are the substructures identified with SUBFIND. Since we selected them by applying a cut in total mass, they can be associated with galaxies (there will be a subtle difference in the selection function). Specifically, we considered all substructures with total masses >10¹¹ h⁻¹ M_⊙, which is a good compromise between the resolution of our simulations and the need for large statistics of tracers.

• Smoothing bandwidth h: This parameter controls the smoothing of the distribution of tracers and thus depends on the local property of the density field. We utilised the expression provided by Chen et al. (2015a): $$h=A_0{\left(\frac{1}{d+2}\frac{1}{n}\right)}^{\frac{1}{d+4}}{\sigma }_{min},$$(2)

  where n and d have the same meaning as before, σ_min is the minimal value of the three standard deviations associated with the spatial coordinate, and A₀ is the only parameter that does not depend on the local environment and needs to be calibrated according to the analysed simulation. In our case, we chose A₀ = 0.5 after visually inspecting the overlap of the resulting skeleton and the tracer density field (see Appendix B). The usual value of the smoothing bandwidth h is in the range 1–1.6 Mpc.

• Threshold parameter τ: When defining a filament in a selected region, we can a priori avoid the computation in lowest-dense volumes. This translates to discarding all x locations whose density function is below a certain threshold, p(x) < τ. Applying this selection avoids false filament detections, reduces overall noise, and speeds up the process. The value of the threshold parameter again depends on the local properties of the environment. Namely, it is defined as the difference between the local, $\hat{p}$, and the mean density, $\overline{p}$: $$\tau =\sigma (\hat{p})\equiv {({\int }_{\mathbb{K}}(\hat{p}(x)-\overline{p}(\mathbb{K}){)}^2\mathrm{d}x)}^{1/2}\sim \hat{p}-\overline{p},$$(3)

  where $\mathbb{K}$ is the region to which we are applying the algorithm, and the local density is computed in 1 cubic cell with a size equal to the size of $\mathbb{K}$ divided by 100 (see definition of ℳ in Appendix A).

The SCMS algorithm works on a mesh of points that are shifted towards the density ridges defined by the tracers (see Appendix A for a description). The final output of the algorithm is a collection of points. In Fig. 1, we show one example of the output of the SCMS algorithm in 2D space. The skeleton is extracted from a thin slice (~1.2 Mpc thick) in region D6 at redshift z = 0. The black dots represent the skeleton of the cosmic web. They are overimposed on the colour-coded KDE and the tracing substructures, which are marked with green dots. A thin cyan line separates the regions with p(x) > τ, where τ = 9.78 · 10⁻⁶. The blue circle represents the virial radius of the massive clusters in the region. We note that some skeleton points lie inside their virial radius. A large part of the skeleton follows the densest regions and the tracer distribution. We noticed, though, that there are some filaments that, while in the dense region, do not appear to be close to tracers. Most of them are very short (some are even single points), and since they are the result of the noise, at later stages, they will be removed.

Fig. 1 Example of the final result of the SCMS algorithm in 2D space. Black dots mark the points of the skeleton. Along with the skeleton, the colour-coded kernel density estimator is plotted. Tracers are marked with green points. The contour line is plotted with cyan colour at the threshold τ value. The blue circles represent the virial radius of the two massive clusters (Mvir > 1013 M⊙ in the 2D space. This skeleton was extracted from a thin slice (~1.2 Mpc thick) in region D6 at redshift z = 0.

3.2 Sequential chain algorithm for resolving filaments

As shown and presented above, the outcome of the SCMS algorithm is a skeleton consisting of an ensemble of points that need to be automatically sorted and organised into filaments. For this purpose, we developed a sequential chain algorithm for resolving filaments (SCARF).

The SCARF algorithm is a two-step process: the initial step involves combining neighbouring points into a singular chain (see Sect. 3.2.1), while the second step breaks down the single chain into separate filaments (see Sect. 3.2.2). A chain is a sequential arrangement of points where each point is next to its closest neighbour, forming a connected series like a linked chain. This two-step approach allows SCARF to organize and sort points of the skeleton into filaments in an effective and automated manner.

3.2.1 Constructing a single chain

To facilitate the explanation, we start by presenting the main concepts used in the procedure, which are also displayed in Fig. 2. Considering the entire collection of points forming the skeleton of the cosmic web or, in other words, the outcome of SCMS, the first step is to identify the closest neighbour to the point of reference. This neighbour point is always referred to as g. To identify the chain, we use the letter c. The final result of SCARF is one unique chain for all points in the skeleton. However, at the intermediate steps (one of them is shown in Fig. 2), more chains can be built, and these are distinguished with a superscript indicating their numberings, for example c⁰ and c¹ would be the first and the second chain that is created, and so on. In each step, the initial and final points of the i-th chain are respectively called $c_0^i$ and $c_{-1}^i$. These points change through time as the chain becomes longer. In some instances, the chain could be constituted of only one point and then $c_0^i\equiv c_{-1}^i$.

The procedure starts by selecting a random point, which by definition is associated with the chain c⁰. Subsequently, the closest companion, the point g, is identified, and its distance with respect to c⁰ is computed. The closest companion to g is searched, and if it is farther away than the distance between g and c⁰, then g is associated with the first chain c⁰, which now constitutes two points, and the search for another neighbour continues. Otherwise, the point g is not connected to c⁰, but it is considered the first point of a second chain, c¹.

Once two separate chains are identified, the procedure focuses on the closest neighbour of c¹, called again g. The program then compares the distance between the new g and c¹ with the distance between c¹ and the end point of the chain c⁰. If the latter is smaller than the former, then c¹ is connected to c⁰, and the two separate chains become one, named c⁰. In the other case, g is attached to c¹ only if it does not have any other closest neighbours. Indeed, in this circumstance it is the starting point of a third chain c².

Generalising the concept and the notation, whenever a new neighbour g is found, the program computes its distance with its closest chain cⁱ as the minimum distance between g and the chain final points: $$d_{g,c}=min\{\parallel \mathbf{g}-{\mathbf{c}}_0^i\parallel ,\parallel \mathbf{g}-{\mathbf{c}}_{-1}^i\parallel \}.$$(4)

This distance is compared with the distance between two previously identified subsequent chains, defined as the minimum distance between their final points: $$d_{c,c}=min\{\parallel {\mathbf{c}}_0^i-{\mathbf{c}}_0^{i-1}\parallel ,\parallel {\mathbf{c}}_0^i-{\mathbf{c}}_{-1}^{i-1}\parallel \parallel ,\parallel {\mathbf{c}}_{-1}^i-{\mathbf{c}}_0^{i-1}\parallel ,\parallel {\mathbf{c}}_{-1}^i-{\mathbf{c}}_{-1}^{i-1}\parallel \}.$$(5)

If d_c,c is smaller, then the two chains are fused into one. Otherwise, the program checks whether g has any other companion at a distance smaller than d_g,c. If it does not, then g is associated with cⁱ; otherwise, g will be the starting point of the new chain cⁱ⁺¹.

The procedure continues until all points from G are connected in one single chain.

Fig. 2 Visual presentation of an intermediate step of the first part of SCARF showing the key concepts on which the algorithm is based: the skeleton G (all points); chains c0, c1, and c2; chain ends ${\mathbf{c}}_0^0$ and ${\mathbf{c}}_{-1}^0$; and nearest neighbour g (red dot). The chain c1 only has one point, which is fairly common when sorting points. Grey points are part of the skeleton G but are yet to be sorted.

Fig. 3 Final result of the SCARF algorithm applied to the same 2D region of Fig. 1. Different colours mark different filaments, while in a thin grey line, we show the chain identified in the first step of SCARF.

3.2.2 Resolving the filaments

Once the entire skeleton is regrouped into a single chain, the second step of SCARF fragments the chain into individual filaments. The procedure starts by computing the vectors between two subsequent points: ξ_i; = c_i+1 − c_i. Each vector has two main properties, length ||ξ|| and direction, from which we can evaluate the change of direction of two neighbouring vectors as the angle: $${\theta }_i=\arccos ({\mathit{\xi }}_i\cdot {\mathit{\xi }}_{i-1}).$$(6)

The criteria used to break the chain into filaments are that either ||ξ|| > 1 Mpc or θ_i < π/6. We select the first threshold to be of the order of the cluster radii, while the second limit is the minimum angle that still allows for smooth transitions of directions whenever the angle is below 30 degrees.

In Fig. 3, we show the final results of the SCARF algorithm in 2D space where each filament is shown with a different colour. The points that lie inside the virial radius of the clusters with mass M > 10¹³ M_⊙ (in this case, two) are removed. A low-opacity line represents a chain that connects all the skeleton points detected by the SCSM algorithm. Filaments arising from the noise are marked with small, low-opacity points. Because of their shortness, they are easily removed. In this particular 2D case, SCARF detected 21 filaments. The choices limiting the vector length and direction change angle are applied to all 3D skeletons for all the regions analysed.

3.3 Mapping the filament

A filament consists of N points in successive order f₀, f₁,..., f_N. We now presume that the filament is shaped like a bent cylinder. We denote a particle’s position as p. The longitudinal distance l for a particle in p (e.g. gas or dark matter particle) is defined as the distance between the filament’s starting point f₀ and the particle’s projected point on the filament’s spine and is calculated as $$l=\sum _{m=0}^n\parallel f_{m+1}-f_m\parallel +\frac{(\mathbf{p}-f_n)\cdot (f_{n+1}-f_n)}{\parallel f_{n+1}-f_n\parallel },$$(7)

where n is the sequence number of the point in the filament. Points f_n and f_n+1 are the closest filaments’ points to the projected point p. The radial distance r is the shortest distance between the particle and the spine of the filament and is calculated as $$r=\frac{\parallel (\mathbf{p}-f_n)\times (f_{n+1}-f_n)\parallel }{\parallel f_{n+1}-f_n\parallel }.$$(8)

3.4 Length and radius of the filament

We defined the length of the filament as the sum of distances between neighbouring points: $$L=\sum _{n=0}^N\parallel f_{n+1}-f_n\parallel ,$$(9)

where N is the total number of points in the filament’s spine. This differs from the usual definition of length, where filaments are connecting the nodes. SCMS does not define nodes, and this is the reason why we also removed points that lie inside the virial radii of groups and clusters, as described in Sect. 3.2.2. Consequently, the filaments may be slightly shorter than expected compared to other studies. Nonetheless, this definition still offers a reliable estimate of the filaments’ geometrical properties.

To estimate the filaments’ radius, we first mapped the gas, dark matter and star particles surrounding the filament, as discussed in Sect. 3.3, up to 5 Mpc from the filament spine. We obtained the particle distribution around the filament in l and r coordinate system. Our goal is to determine the filament’s radius based on its overdensity, defined as $$\delta =\frac{\rho }{{\rho }_{\text{crit }}}-1,$$(10)

where ρ_crit is the critical density³ at the given redshift. This allowed us to normalise the overdensity over all different regions, making them comparable.

We create an overdensity map that is described using the longitudinal and radial coordinates, l and r. The map is divided into cells with a height and width of approximately 0.1 Mpc. We mapped the particles at the same radial distance r from the spine in the cell at the correspondent longitudinal distance l. We applied the Gaussian filter with σ = 0.2 Mpc, equivalent to the size of two neighbouring cells, to smooth the map and damp sharp irregularities. We then categorised the cells as overdense if their overdensity is δ > 0 and underdense otherwise.

Ideally, the filament would include most of the overdense cells close to the spine relative to the entire map. We can take two approaches: either determine a constant radius using the cells of the entire map or calculate a variable radius R(l) using cells associated only with that l. In practice, since R(l) can change rapidly, we applied a Gaussian filter with σ = 1 Mpc, allowing for a more gradual change in the filament’s radius.

The computation of both radii is shown in the following example. The main panel of Fig. 4 illustrates an example of the colour-coded overdensity map around a single filament, with length L ≈ 6.2 Mpc, in the D6 region (AGN simulation). The purple colour represents the underdense cells. The solid blue line shows the variable radius R(l). We also evaluated a constant radius, indicated by the dashed blue line, in the following way. In the right panel, for each radius r we computed the frequency of overdense cells for the upper area (red, above r) and lower area (green, below r). The constant radius is shown as the radius where the difference in the frequency of overdense cells is maximum (yellow line). In the lower panel, the calculated variable radius is represented by a solid line, while the dashed line illustrates the radius after applying the Gaussian filter. This example shows that the constant radius cuts off part of the overdense region, while the filament’s variable radius traces the filament’s radial boundaries better.

However, issues can arise when no overdense cells are present at l, thus getting R(l) = 0, or when objects, such as galaxies or small groups, are located outside the filament (red clumps in the overdensity map in Fig. 4). In that case, we can encounter a lone overdense cell (or a few) located far away from the filament’s spine; because of that, the radius will extend beyond the actual filament. This can be mitigated by imposing a hard limit, where if the first few cells from the filament’s spine are under-dense, then the radius at that point is R(l) = 0. This overcomes, to some extent, the issue. In the cases where the average over-density of the filament is below the threshold δ < 0, we chose to remove the filaments from further analysis. On average, we had to remove two additional filaments per snapshot.

Fig. 4 Map of overdensity. The purple regions are below the selected overdensity threshold. The map also contains the computed constant radius, dashed blue line, and the variable radius R(l), indicated with a solid blue line. In the right panel, at each r, we have plotted the frequencies of overdense cells (red line for cells above that r and green line for cells below r) and their difference (yellow line). The constant radius is shown as the radius where the difference in overdense cells is maximum. The variable radius R(l) is plotted at the bottom of the figure. The grey solid line is the calculated radius, while the black dashed line is the radius with the applied Gaussian filter. The dark and red clumps at larger r correspond to a larger overdensity of particles related to substructures. We obtained this map for one of the filaments in the D6 region (AGN simulation).

Fig. 5 Phase diagram of hydrogen number density nH (x-axis) and temperature T (y-axis) for all gas particles in all regions at redshift z = 0 in AGN (left) and CSF (right) simulations. The dashed lines mark the limits for each gas phase defined in the text (Hot/WHIM/other). The colour-code (right) refers to the probability density.

3.5 Physical properties of the filament

To see if the filaments are the fair containers of cosmic baryons, we computed the gas, stellar and baryon depletion factors. They are calculated as in Planelles et al. (2013): $$Y_{\mathrm{x}}=\frac{M_{\mathrm{x}}/M_{\mathrm{t}\mathrm{o}\mathrm{t}}}{{\mathrm{\Omega }}_{\mathrm{b}}/{\mathrm{\Omega }}_{\mathrm{m}}},$$(11)

where x can be gas, stars or baryons (gas and stars combined), M_x is the total mass of the x component, M_tot total mass of the filament, and Ω_b and Ω_m the density parameters used in the simulations.

We analysed the properties of the gas included in a larger region around the spine of the detected filaments. For the purpose of this work, we studied the gas components, dividing them into three phases based on temperature and hydrogen number density n_H, as done, for example, in Martizzi et al. (2019). The hot gas phase corresponds to gas having temperatures above T > 10⁷ K and any number density. This phase is common in the potential wells of massive clusters, but as we later show, we can find this phase in filaments as well, especially near the filament ends (if they are near a cluster). For the WHIM gas phase, we used the range 10⁵ K < T < 10⁷ K for temperature and n_H > 10⁻⁴(1 + z) cm⁻³ for the hydrogen number density. We expect most of the gas to be found in this phase. The remaining gas was treated separately and includes the colder gas associated with the ISM of galaxies and cold IGM, as well as the warm CGM created by shock heating and feedback processes near galaxies. In Fig. 5 we plotted the phase diagram of hydrogen number density n_H and temperature T for all gas particles in all regions combined at redshift z = 0 in AGN and CSF simulations. Dashed lines separate the different gas phases that we have defined previously: hot, WHIM and other. Here, we can also note the difference between AGN and CSF simulations. The biggest difference is seen at the tail-end of the diagram, where high-density and low-temperature gas is located. This is the star-forming gas, which is more abundant in CSF simulations, while it is efficiently removed from the star-forming phase by AGN feedback.

To show the different behaviour of each gas component inside the filamentary structures, we calculated the mass fraction for each gas phase as M_x/M_gas, where x is the mass fraction of the desired gas phase and M_gas the total gas mass inside the filament.

To explore the metallicity content of filaments, we obtained the iron metallicity Z_Fe for each gas phase in the filament as $$Z_{\mathrm{F}\mathrm{e}}=\frac{\sum _i{m}_{\mathrm{F}\mathrm{e},\mathrm{i}}}{\sum _i{m}_{\mathrm{H},\mathrm{i}}},$$(12)

where the sum is done over all the particles contained in the region of the filament and where the solar abundance is Z_Fe,⊙ = 1.77 · 10⁻³ as in Grevesse & Sauval (1998). To be consistent with the observed temperatures, we also used the definition of spectroscopic-like temperature T_SL as described in Mazzotta et al. (2004) (see also Rasia et al. 2005): $$T_{\mathrm{S}\mathrm{L}}=\frac{\sum _i{\rho }_{\mathrm{i}}{T}_{\mathrm{i}}^{0.25}{m}_{\mathrm{i}}}{\sum _i{\rho }_{\mathrm{i}}{T}_{\mathrm{i}}^{-0.75}{m}_{\mathrm{i}}},$$(13)

where T_i is the temperature of the gas particle, ρ_i its density and m_i its mass.

4 Results

First, we applied the SCMS algorithm on the selected set of the Dianoga regions, obtained with two different physics (AGN and CSF), as described in Sec. 2.1. To study the evolution of gas properties, we selected redshifts z = 0, 0.25, 0.51, 0.76, 1, 1.26 and 1 .48. We removed all the points of the skeleton that fall inside the virial radii of the groups and clusters with masses above 10¹³ M_⊙. After that, we ran the SCARF algorithm on the remaining skeleton. The final output is one catalogue of filaments per region, redshift and physical model.

We calculated the physical length L of all the filaments as specified in Sec. 3.4. We only kept filaments with length L > 2 Mpc for further analysis. Shorter filaments could be, at best, considered galaxy bridges and were entirely disregarded. With this, we also removed the false filaments, as discussed in Sec. 3.1. Almost all of the detected filaments are inside the high-resolution regions entirely, as they are either connecting two massive halos within this denser region or are situated in between filaments. A few of the filaments do not have a cluster or group at one of the ends. We still consider these filaments, as our main focus is the study of diffuse gas.

4.1 Geometrical properties

Before analysing the geometrical properties of filaments, we provide in Table 2 the number of filaments with L > 2 Mpc found per region, redshift and simulation type. We do expect that the number of filaments is related to the volume of the region. However, we see a steady increase in time of the filament absolute number in all regions, the numbers growing by ~60% for the smallest (“isolated”) regions, D5 and D9, and ~40% for the other three largest regions. This is expected since the density contrast increases with time, especially in isolated regions (see, e.g. Galárraga-Espinosa et al. 2024). We checked that, in general, filaments become longer with lower redshift in proper coordinates, though not considerably. There is no significant difference in the evolution of the number of filaments between AGN and CSF simulations. We do observe slight deviations between AGN and CSF simulations, but the difference is, in most cases, below 5%. Since the formation of the cosmic web is mostly determined by gravity, a similar number of filaments is expected as physical models used in AGN and CSF simulations mostly affect the gas component.

Along with the filaments’ length, we can also observe their shape (whether the filament is straight or curved). Therefore, we computed the distance D between the filament’s first and last point and compared it with the filament length L.

Since, in this case, we are not interested in environmental, model or redshift dependency, we combined all filaments with L ≥ 2 Mpc from all regions and redshifts. In Fig. 6, the 2D histogram between D and L is plotted, with dashed lines representing the values of D/L and solid black line representing the median values of D at given L. On the side are the histograms for L and D. The ratio D/L represents the shape of the filament; the lower the value, the more curved the filament. As expected, since we are working with a low-volume simulation box, most of the filaments have lengths below 9 Mpc, with number counts falling with length. Similarly, this happens with distance D since this quantity is correlated with the length of the filament. According to multiple works (Galárraga-Espinosa et al. 2021, 2024; Wang et al. 2024), filament populations can be divided into short and long filaments. Though there is no universally set limit, in all of the works, long filaments have a length of at least 20 Mpc. This means that in this work, we exclusively study short filaments. The median line (solid black line in Fig. 6) indicates that the longer the filament, the more likely it is to be curved. The median value is close to D/L = 0.9 for short filaments. With larger lengths, the median value is moving towards the lines with lower values for D/L, indicating that the longer filaments are more likely to be curved since they do interact tidally with clusters located close to the filaments (Colberg et al. 2005; González & Padilla 2010; Cautun et al. 2014). Even so, there are a number of longer filaments that are straight and smaller filaments that are curved.

We explored how the shape of the filaments changes with redshift. Shown in Fig. 7 are median values plotted against redshift z for AGN and CSF simulations. The median values were calculated at each redshift combined for all the filaments in all regions. One can observe the filaments’ shape is consistently rising with lower redshift, which means that the filaments become straighter with time. A similar conclusion was reached in Cautun et al. (2014), where the filament’s shape was similarly defined. As expected, the evolution to straighter filaments is similar for both AGN and CSF simulations.

Fig. 6 A 2D histogram between L and D, with N as the colour-coded number of counts. Dashed lines represent the fraction D/L with coefficients of 1, 0.9, 0.8, 0.7, 0.6 and 0.5 from top to bottom, respectively. The black line represents the median value of D at a given length L. On the sides are histograms of L and D. This is for all detected filaments across all regions, redshifts and simulations with L ≥ 2 Mpc. This is a distribution for AGN simulations, but the distribution for CSF is almost identical.

Fig. 7 Median value of shape D/L of the detected filaments in all regions with respect to redshift z for AGN (solid blue) and CSF (dashed red) simulations. The shaded regions are values between the 16th and 84th percentile.

4.2 Global properties

We calculated the radial profile for each filament and then averaged all of them into a single radial profile at a given redshift, the same way as in Galárraga-Espinosa et al. (2020, 2024). The results are shown in Fig. 8, where we plotted the dark matter radial profiles for all filaments from regions D1, D6 and D22 (left panel) and all filaments from regions D5 and D9 (right panel), separately. These radial profiles are in agreement with the results from other works (e.g. Tanimura et al. 2020a; Galárraga-Espinosa et al. 2020, 2024; Wang et al. 2024). For all of the regions, as the redshift decreases, we observe a decline in physical density. It is also clear that the filaments around the more isolated objects have a higher density close to their spine. Still, they fall quicker to values similar to the other three regions at larger radii. In the same figure, we also observe the widening of radial profiles with decreasing redshift. This indicates that, on average, the radius of filaments increases in time. This aligns with the findings in Wang et al. (2024), where they examined the radii of short, medium, and long filaments in the MilleniumTNG simulations. For short filaments (median value 7.56 Mpc), which most of the filaments in our study are, the radius decreases until redshift z = 2 and then begins to increase. We observed similar trends for medium and long filaments, with the turning point occurring at redshift z = 1. However, for long filaments, the radius remains relatively unchanged after this redshift.

The following computations were performed for all identified filaments across five regions at seven distinct redshifts. We begin by calculating for each filament the mean radius, denoted as ⟨R⟩, by taking the mean of the variable radii. We verified that, similarly to what is shown in Fig. 8, ⟨R⟩ increases with redshift. Subsequently, we gathered all particles located within the filaments’ variable radius and collected their masses to obtain the total mass of the filament. Given that we know the filament’s mass and volume, we computed its average density ρ and its average overdensity δ.

In Fig. 9, the top-left panel displays the distribution of total filament mass against its length for filaments at redshift z = 0 in AGN simulations. The results for CSF simulations are similar. The graph includes median mass values within specific length ranges, revealing a consistent increase in mass. A power-law fit indicates a relation of M ∞ L^1.7. Cautun et al. (2014) determined a similar power-law relation, ∞ L^2.2, based on filaments in Millennium simulations detected with a different algorithm, NEXUS+. Notably, their simulations featured longer and, consequently, more massive filaments, a larger sample size, and exclusively dark matter-filled filaments. The existence of this correlation spanning a different range of lengths and masses implies that the longest filaments are not a fortunate connection of short filaments but rather clearly defined and possibly persistent structures. The bottom-left panel illustrates median lines for seven redshifts. Filaments consistently adhere to the same mass-length relation but exhibit decreasing mass with lower redshifts.

Our results might seem to go against the expectation that filaments become thinner as time passes, becoming more massive and diminishing radii with redshift. However, we should keep in mind that we are not studying large-scale filaments from cosmological boxes but relatively shorter filaments in denser regions populated with clusters. This aspect is likely influencing the evolution of filaments. Massive clusters will grow thanks to the mass flow from filaments, but we expect the filaments will not accrete mass from surroundings at the same rate. This would require a deeper investigation that is out of the scope of this paper. Moreover, as noted also in Cautun et al. (2014), the correlation between the width and density of filaments depends on evolutionary processes. They found that thin filaments are usually found in underdense regions (areas where large-scale filaments form and persist), while thick ones are in overdense areas. The area we are studying is rather dense and populated with a large number of clusters that act as gravitational attractors.

A similar plot is presented for mass versus mean radius in the middle panels of Fig. 9. Most of the mean radii are below 1 Mpc, which is a standard radius for the filaments, most recently determined by Wang et al. (2024). Still, the trend indicates that larger mean radii correspond to more massive filaments, which also correspond to longer filaments. Additionally, the radius of a filament appears to increase with redshift at the same mass, which is consistent with Fig. 8. The panels on the right depict mass versus overdensity. Initially, the mass appears constant at lower overdensity but rises steadily with overdensity until it becomes nearly constant again. We note that, as seen in the upper right panel, most of the filaments have overdensity larger than δ > 2.

Fig. 8 Mean radial profiles of dark matter for all filaments at different redshifts. The left graph has filaments from regions D1, D6, and D22, while the right has filaments from regions D5 and D9. These radial profiles were extracted from AGN simulations but are the same in CSF simulations. The radius and density are both in proper units.

Fig. 9 Filament mass with respect to its length, mean radius, and overdensity at redshift z = 0 in AGN simulations. Top panels: distribution of mass of the filament M versus its length L (left), mean radius ⟨R⟩ (middle) and overdensity δ + 1 (right) at redshift z = 0. A dashed curve represents a median mass at a given filament range in x-axes. A fitted power law curve is plotted on a mass versus length graph, with a solid line representing the relation M ∞ L1.7. The red dots represent the discarded filaments with δ < 0. Bottom panels: median mass values of filaments with respect to their length L (left), mean radius ⟨R⟩ (middle), and overdensity δ + 1 (right) at different redshifts. For redshifts, the different colours and lines are reported in the legend (middle plot). All the filaments on the plot are from the AGN simulation.

4.3 WHIM and hot gas phase properties

To understand the broader evolution of the mass in filaments, we computed the fractions of the main matter components: dark matter, gas and stars within the boundary of the variable radius R(l). We noticed that the filaments located at the borders of the high-resolution regions have a very low number of gas particles. Therefore, we removed those with $\frac{M_{\text{gas }}}{M_{\text{tot }}}<0.05$ from further analysis. We then calculated the depletion factor for gas, stars and baryons using Eq. (11) for both AGN and CSF simulations. The results at redshift z = 0 are shown in Table 3, where one can notice that filaments are not fair containers of the baryons in both simulations since Y_baryons is around ~0.79 in AGN and ~0.83 in CSF simulations lower than unity. This is broadly consistent with the result by Galárraga-Espinosa et al. (2022), where they observed short filaments in IllustrisTNG simulations having Y_baryons ~0.9. They concluded that there is a baryon deficiency at r < 0.7 Mpc, which coincides with the typical radius of the filaments we study. The gas in the filaments within AGN simulations has a comparable fraction with that of the CSF simulations. In turn, Y_stars in CSF filaments is almost twice the value of that in AGN simulations, both having a large error. This is expected since CSF simulations are characterised by overcooling and missing the AGN feedback needed to regulate the star formation (e.g. AGN feedback), and overdense regions have higher production of stars. This increases the abundance of stars in the filaments in CSF simulations while simultaneously reducing their gas abundance. On the other hand, the reduced star formation in the AGN simulations makes more baryons available to be part of the diffuse medium, but such extra baryons are, in fact, displaced by the same AGN feedback. As a result, while Y_stars in AGN simulations decreases, there is almost a compensation of the two effects in determining the value of Y_gas. One can notice that at the end, the overall baryon content of filaments is slightly reduced in AGN simulations.

We checked that the mean mass fractions of baryons and dark matter do not change much in time (around 5% in AGN and 10% in CSF simulations with respect to the present day value). The most notable rise of the mass fraction is for the stars in the CSF simulations, which is expected due to the high star-formation rate. Because of the relatively small mass fraction in stars, this does not visibly impact the gas and dark matter mass fraction. We then calculated the mass fractions of the three gas phases, defined in Sec. 3.5 and calculated the median mass fraction for each gas phase and the 16th and 84th percentiles for AGN and CSF simulations at each redshift. The results are shown in Fig. 10. We immediately notice that the gas phases follow slightly different evolutions depending on the simulation. This is broadly consistent with what was found by Tornatore et al. (2010), where slightly different simulations were analysed in a broader context of the WHIM content inside smaller cosmological boxes. In the AGN simulation, the WHIM gas phase consistently stands out as the most abundant across different redshifts, while the hot gas phase remains at very low levels. We recall that the filaments in this work are outside the virial radii of groups and clusters. Notably, the WHIM gas phase becomes even more abundant as redshift decreases. At redshift z = 0, its median mass fraction reaches almost 0.6–0.7 in the two physical models. In similar works, this value is between 0.5–0.7 (Martizzi et al. 2019; Galárraga-Espinosa et al. 2021). Meanwhile, the presence of other gas phases gradually decreases. Though the abundance of the hot gas phase is barely noticeable, it does increase slightly at lower redshift. This indicates that the gas temperature in the filaments is slowly rising with lower red-shift, likely because of the shock heating (Martizzi et al. 2019). Due to that, the filaments would possibly be easier to detect with X-ray telescopes at lower redshifts. The two models show very similar trends with redshift. However, the fractions of the WHIM and the colder gas reach equality earlier for the AGN simulation.

One of the most striking differences between the two models can be seen in the filament’s metal content. Furthermore, metallicity will be important in the observational detection of filaments. Therefore, we calculated the metallicity Z_Fe for each gas phase in every filament as well as the median values at each redshift. The results are shown in Fig. 11 for AGN and CSF simulations, respectively.

In AGN simulations, the hot gas phase exhibits the highest metallicity among the three considered phases, with the other gas phases ranking second and the WHIM displaying the lowest metallicity. Over time, the metallicity of the hot gas gradually decreases, while both the WHIM and other gas phases experience a slow but steady rise in metallicity. At redshift z = 0, the metallicity of the hot gas phase falls below the value of other gas phases. These results align with those of Martizzi et al. (2019), who found that the decrease in metallicity at lower redshifts is attributed to increasing temperatures caused by shock heating. This results in the transition of the metal-poor WHIM phase (as well as the warm CGM phase) to the hot phase. Since hot gas is not abundant in the filaments, as shown in Fig. 10, the addition of metal-poor gas decreases the median metallicity.

Conversely, the hot gas phase in CSF simulations shows much lower metallicity than in AGN simulations. It has barely larger values than WHIM, which also has lower metallicity than in AGN simulations. Other gas phases, though, exhibit similar values at all redshifts in the two physical models. Notably, across all three gas phases in CSF simulations, there is a collective increase in metallicity over time, marking a distinctive difference from the trend observed in AGN simulations, likely due to a higher star-formation rate, resulting in larger metal production. The overall metallicity values are higher in AGN simulations compared to CSF simulations, highlighting the effect of the AGN feedback model in simulations. This is consistent with the previous results shown in Biffi et al. (2018) where the AGN feedback, especially active at z = 2, ejected metal-enriched gas from the small potential well of high-redshift galaxies. It is, therefore, expected that in AGN simulations, filaments retain gas with high metallicity, funnelling it towards clusters.

Focusing on the AGN model, although the hot gas inside filaments is not a relevant part of the gas mass, the filaments seem to be enriched with high metallicity hot gas. Gas in the WHIM phase is instead reaching the highest fraction of gas mass at redshift z = 0 while having an almost constant metallicity with time. At the redshift of ~1.5, WHIM is already enriched. Its gas mass then grows by encompassing gas that gradually reaches the temperatures typical of this phase while being enriched as well. The same is not true for CSF simulations, where both the mass fraction of WHIM and its metallicity are growing in time. In this case, the WHIM phase is gradually enriched by the accumulation of gas with higher metallicities.

To better explore the relation between the metal content and the temperature of the gas, we show in Fig. 12 the metallicity Z_Fe versus spectroscopic-like temperature T_SL for the WHIM and the hot gas phase separately. The spectroscopic-like temperature was calculated according to Eq. (13), where only appropriate gas particles that fall into the designated temperature and density range are included, i.e. only WHIM particles are used for the calculation of the T_SL for WHIM phase.

On the upper left panel of Fig. 12 (for AGN simulations), we see that the WHIM undergoes a swift rise in metallicity from low temperatures up to approximately T_SL ≈ 3 · 10⁵ K. The metallicity stabilises for larger temperatures, resulting in a relatively constant metallicity. Both trends can be attributed to shock heating. According to Martizzi et al. (2019), the only gas phase with metallicity lower than the WHIM is the diffuse IGM, which has a similar density to the WHIM but with temperatures below 10⁵ K. As filaments undergo shock heating, their temperature increases, causing the metal-poor diffuse IGM to transition into the WHIM phase. This transition is most evident at lower temperatures but diminishes rapidly at slightly higher temperatures. Since shock heating does not add metals to the gas, the metallicity is expected to remain constant as the temperature rises. The metallicity of the WHIM is constant across all redshifts, consistent with results in Fig. 11. The hot gas phase in AGN simulations (upper right panel) displays variable median values at different redshifts. At higher redshifts, we observe a slight rise in metallicity with temperature and a gradual decline with respect to temperature at lower redshifts. We note that the majority of the filaments have a spectroscopic-like temperature below 2 · 10⁷ K. Similar to the discussion above, the metallicity of the hot gas is decreasing due to the shock heating. This is especially apparent at z = 0, where we notice that the median metallicity is only slightly higher than that of WHIM and is constant at all temperatures.

On the lower panel of Fig. 12 (CSF simulations), the metallicity of the WHIM gas phase shows a similar trend as in AGN simulations. The metallicity sharply rises until the temperature is T_SL ≈ 2 · 10⁵ K. At higher temperatures, the metallicity rises very slowly with temperature. In time, the metallicity steadily rises at all temperatures. For the hot gas phase, the metallicity rises with higher temperatures but seemingly shows less dependence on redshift.

For the WHIM phase, we find a significant diversity between the AGN and CSF simulations. In AGN simulations, the iron abundance within filaments remains unaffected by redshift, while CSF simulations show a slow rise of the values in time. This difference can be attributed to the role of AGN feedback, which was explored in the same simulations by Biffi et al. (2018), accordingly to what was discussed for Figs. 11 and 12. The AGN feedback is able to pre-enrich gas at z > 2, and consequently, the metallicity level is already high in the AGN model for the WHIM phase.

Fig. 10 Median values (solid lines) of gas mass fractions where x is WHIM, hot, and other gas phases. The left panel shows results for AGN simulations and the right panel results for CSF simulations.

Fig. 11 Median metallicity ZFe for WHIM, hot, and other gas phases in redshift dependence. The left panel shows results for AGN simulations, and the right panel shows results for CSF simulations.

Fig. 12 Median values of metallicity ZFe as a function of spectroscopic-like temperature TSL for both WHIM and the hot gas phase at seven different redshifts. The top panel is for the AGN simulations, and the bottom is for the CSF simulations.

5 Conclusions

In this work, we have studied the filaments in a set of highresolution regions belonging to the Dianoga simulations that were centred around massive clusters. The set is composed of five different regions extracted from a parent cosmological simulation. Two of the regions were selected as ‘isolated’, with a less massive cluster in their centre, while three of them are centred on massive clusters and densely populated with satellite clusters.

We used the SCMS algorithm Chen et al. (2015a) and a SCARF algorithm to extract the single filaments. From there, we were able to quantify several geometrical and physical properties: length, shape, mass content for dark matter, gas and stars, mean radius, and mean density. We focused on the diffuse components (WHIM and hot gas) in the filaments, extracting their mass fractions and determining the mass-weighted metallicity and spectroscopic-like temperature.

We summarize our main findings as follows:

1. On average, the filaments inside the selected regions become longer in time. We found that longer filaments are more likely to be curved. With a lower redshift, the filaments become straighter on average.

2. The dark matter radial profiles in filaments have different central densities for regions that evolved as isolated with respect to the densest populated ones. Moreover, we noticed that the mean radius of filaments increases during evolution, which is consistent with recent results (Wang et al. 2024). This behaviour is related to the fact that our regions are regions around massive galaxy clusters instead of cosmological boxes covering a fair sample of environments.

3. We find that the mass of filaments is correlated with their length and their size. Filaments’ mass increases with its length according to the power law M ∞ L^1.7, similar to what was found by Cautun et al. (2014). This indicates that the longest filaments are clearly defined and not a fortunate connection of shorter filaments. The mass also increases with a larger mean radius ⟨R⟩. Both correlations show a redshift dependence: at a fixed length or radius, filament masses decrease in time.

4. We explored the evolution of the gas found in hot, WHIM, and colder gas phases inside the detected filaments. The predominant gas phase in filaments is the WHIM phase. Looking at the evolution of the WHIM mass fraction, it steadily rises towards lower redshifts. We find that the physical models do not impact the amount of baryon mass in the WHIM phase.

5. In general, we find that the metallicities of the WHIM and hot gas phases are higher in the AGN model than in the CSF one.

6. The evolution of the iron abundance in filaments depends on the physical model used in simulations. In the case of AGN simulations, the WHIM metallicity remains constant with redshift. This is connected with the AGN feedback model, which is able to displace metals at very high redshifts (e.g. Biffi et al. 2018). The already enriched gas populates the WHIM gas phase, resulting in an almost constant metallicity. The opposite is true for CSF simulations, where enriched gas is not circulated by powerful outflows, increasing the iron abundance over time.

7. In the case of the WHIM, metallicity is positively correlated to gas temperature, up to T_sl ~ 3 · 10⁵ K. At higher temperatures, the metallicity remains constant.

In this work, we have shown that in denser regions around massive galaxy clusters, a number of filaments can be traced. The different baryon processes included in the Dianoga simulation can be used to understand the impact of AGN feedback. The trends reported in this work, especially about the gas phases and their evolution, should be investigated in more detail in future studies. To better understand the evolution of filament properties, we plan to select single filaments that persist with redshift. These filaments, possibly connected on both ends to large virialized structures, will give us a deeper insight into the regions that connect cluster outskirts and the large-scale structure. Moreover, detecting filaments near galaxy clusters in simulations will give us the opportunity to compare them with different observational data (e.g. X-ray and/or SZ observations of cluster bridges) that are currently becoming available.

Acknowledgements

We thank the anonymous referee for the constructive comments that helped improving the paper overall, especially the presentation of the method and of the results. S.I. acknowledges the support from the Slovenian national research agency ARRS through grant MR-53649. S.I. and D.F. acknowledge financial support from the Slovenian Research Agency (research core funding no. P1-0188). S.B. is supported by: the Italian Research Center on High Performance Computing Big Data and Quantum Computing (ICSC), project funded by European Union - NextGenerationEU - and National Recovery and Resilience Plan (NRRP) - Mission 4 Component 2, within the activities of Spoke 3, Astrophysics and Cosmos Observations; by the PRIN 2022 PNRR project (202259YAF) “Space-based cosmology with Euclid: the role of High-Performance Computing”. S.B. acknowledges partial financial support from the INFN Indark Grant. KD acknowledges support by the COMPLEX project from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program grant agreement ERC-2019-AdG 882679 as well as support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy - EXC-2094 - 390783311.

Appendix A The SCMS algorithm steps

Based on the density function p(x) defined by Eq. 1, the SCMS algorithm shifts points to the closest density ridge, converging towards them, forming a skeleton. Here, we summarize the following steps of the SCMS algorithm (see Chen et al. 2015a, for a detailed description):

1. Compute the kernel density estimator $\hat{p}(x)$ via Eq. 1.

2. Select a mesh ℳ of points that the SCMS will shift towards density ridges. To fully cover the box volume, we used the uniformly distributed grid, with a total of 10⁶ points.

3. Remove x ∈ ℳ if p(x) < τ, where τ is a thresholding parameter. This is called thresholding and is used to remove points in low-density regions. This significantly reduces the clutter noise and number of falsely detected filaments.

4. For each remaining x, perform the SCMS:

   • (a)

     Compute the Hessian matrix H(x) with

     $$H(x)=\frac{1}{n}\sum _{i=1}^n{c}_i({\mu }_i{\mu }_i^T-\frac{1}{h^2}\mathbb{I}),$$(A.1)

     where

     $${\mu }_i=\frac{x-X_i}{h^2},c_i=K\left(\frac{\parallel x-X_i\parallel }{h}\right).$$(A.2)

   • (b)

     Perform spectral decomposition on H(x).

   • (c)

     Obtain eigenvectors corresponding to the smallest d − 1 eigenvalues, to form V(x) = (v₂(x),v_d(x)).

   • (d)

     Update x ← V(x)V(x)^Tm(x) + x until convergence, where

     $$m(x)=\frac{\sum _{i=1}^n{c}_i{X}_i}{\sum _{i=1}^n{c}_i}-x$$(A.3)

     is the mean shift vector.

5. The final output is a collection of points that form the skeleton of the cosmic web.

Appendix B Selecting A₀

In this appendix, we show the calibration of the value of the free parameter A₀ = 0.5 and, subsequently, the value of smoothing bandwidth h. Smoothing bandwidth h controls how much we smooth the distribution of tracers when calculating p(x). We are following the guidelines as written in the Appendix of Chen et al. (2015a).

We illustrate this by using a thin slice (around 1.2 Mpc thick) of the D6 region in the AGN simulation at redshift z = 0 (same slice as used in Fig. 1), centred around the largest halo in the region. The number of tracers (substructures) in this slice is n = 351. In Fig. B.1, we show the resulting skeletons of the SCMS for different values of free parameter A₀. In the top left, top right and bottom left panels are the skeleton for A₀ = 1, 0.5 and 0.25, respectively. These values correspond to a smoothing bandwidth of about 2.38 h⁻¹ Mpc, 1.19 h⁻¹ Mpc and 0.59 h⁻¹ Mpc, respectively. The contour lines encircle the areas where p(x) > τ. In the bottom right panel is the direct comparison between all three skeletons. The resulting skeletons show several distinctive differences.

Using a lower value for A₀ better follows the denser regions and more reliably detects filaments in lower-density areas. Still, using a value that is too low will detect more false filaments and possibly oversaturate points at the isolated tracers. We find that a good compromise is reached by adopting the value of A₀ = 0.5. Using values close to the A₀ = 0.5 does not change the skeleton significantly.

Fig. B.1 Skeletons, produced by SCMS, for different values of A0 and therefore smoothing bandwidth h. The top left panel is for A0 = 1 (corresponding to smoothing bandwidth of 2.38 h−1 Mpc), top right panel for A0 = 0.5 (smoothing bandwidth of 1.19 h−1 Mpc) and bottom left panel for A0 = 0.25 (smoothing bandwidth of 0.59 h−1Mpc). Tracer locations used in SCMS are marked with black dots. Contour lines in grey outline regions where p(x) > τ. The bottom right panel provides a direct comparison among the three A0 values, with each value distinguished by a different colour.

References

All Tables

All Figures

Fig. 1 Example of the final result of the SCMS algorithm in 2D space. Black dots mark the points of the skeleton. Along with the skeleton, the colour-coded kernel density estimator is plotted. Tracers are marked with green points. The contour line is plotted with cyan colour at the threshold τ value. The blue circles represent the virial radius of the two massive clusters (Mvir > 1013 M⊙ in the 2D space. This skeleton was extracted from a thin slice (~1.2 Mpc thick) in region D6 at redshift z = 0.

In the text

Fig. 2 Visual presentation of an intermediate step of the first part of SCARF showing the key concepts on which the algorithm is based: the skeleton G (all points); chains c0, c1, and c2; chain ends ${\mathbf{c}}_0^0$ and ${\mathbf{c}}_{-1}^0$; and nearest neighbour g (red dot). The chain c1 only has one point, which is fairly common when sorting points. Grey points are part of the skeleton G but are yet to be sorted.

In the text

	Fig. 3 Final result of the SCARF algorithm applied to the same 2D region of Fig. 1. Different colours mark different filaments, while in a thin grey line, we show the chain identified in the first step of SCARF.
In the text

Fig. 4 Map of overdensity. The purple regions are below the selected overdensity threshold. The map also contains the computed constant radius, dashed blue line, and the variable radius R(l), indicated with a solid blue line. In the right panel, at each r, we have plotted the frequencies of overdense cells (red line for cells above that r and green line for cells below r) and their difference (yellow line). The constant radius is shown as the radius where the difference in overdense cells is maximum. The variable radius R(l) is plotted at the bottom of the figure. The grey solid line is the calculated radius, while the black dashed line is the radius with the applied Gaussian filter. The dark and red clumps at larger r correspond to a larger overdensity of particles related to substructures. We obtained this map for one of the filaments in the D6 region (AGN simulation).

In the text

	Fig. 5 Phase diagram of hydrogen number density nH (x-axis) and temperature T (y-axis) for all gas particles in all regions at redshift z = 0 in AGN (left) and CSF (right) simulations. The dashed lines mark the limits for each gas phase defined in the text (Hot/WHIM/other). The colour-code (right) refers to the probability density.
In the text

Fig. 6 A 2D histogram between L and D, with N as the colour-coded number of counts. Dashed lines represent the fraction D/L with coefficients of 1, 0.9, 0.8, 0.7, 0.6 and 0.5 from top to bottom, respectively. The black line represents the median value of D at a given length L. On the sides are histograms of L and D. This is for all detected filaments across all regions, redshifts and simulations with L ≥ 2 Mpc. This is a distribution for AGN simulations, but the distribution for CSF is almost identical.

In the text

	Fig. 7 Median value of shape D/L of the detected filaments in all regions with respect to redshift z for AGN (solid blue) and CSF (dashed red) simulations. The shaded regions are values between the 16th and 84th percentile.
In the text

	Fig. 8 Mean radial profiles of dark matter for all filaments at different redshifts. The left graph has filaments from regions D1, D6, and D22, while the right has filaments from regions D5 and D9. These radial profiles were extracted from AGN simulations but are the same in CSF simulations. The radius and density are both in proper units.
In the text

Fig. 9 Filament mass with respect to its length, mean radius, and overdensity at redshift z = 0 in AGN simulations. Top panels: distribution of mass of the filament M versus its length L (left), mean radius ⟨R⟩ (middle) and overdensity δ + 1 (right) at redshift z = 0. A dashed curve represents a median mass at a given filament range in x-axes. A fitted power law curve is plotted on a mass versus length graph, with a solid line representing the relation M ∞ L1.7. The red dots represent the discarded filaments with δ < 0. Bottom panels: median mass values of filaments with respect to their length L (left), mean radius ⟨R⟩ (middle), and overdensity δ + 1 (right) at different redshifts. For redshifts, the different colours and lines are reported in the legend (middle plot). All the filaments on the plot are from the AGN simulation.

In the text

	Fig. 10 Median values (solid lines) of gas mass fractions where x is WHIM, hot, and other gas phases. The left panel shows results for AGN simulations and the right panel results for CSF simulations.
In the text

	Fig. 11 Median metallicity ZFe for WHIM, hot, and other gas phases in redshift dependence. The left panel shows results for AGN simulations, and the right panel shows results for CSF simulations.
In the text

	Fig. 12 Median values of metallicity ZFe as a function of spectroscopic-like temperature TSL for both WHIM and the hot gas phase at seven different redshifts. The top panel is for the AGN simulations, and the bottom is for the CSF simulations.
In the text

Fig. B.1 Skeletons, produced by SCMS, for different values of A0 and therefore smoothing bandwidth h. The top left panel is for A0 = 1 (corresponding to smoothing bandwidth of 2.38 h−1 Mpc), top right panel for A0 = 0.5 (smoothing bandwidth of 1.19 h−1 Mpc) and bottom left panel for A0 = 0.25 (smoothing bandwidth of 0.59 h−1Mpc). Tracer locations used in SCMS are marked with black dots. Contour lines in grey outline regions where p(x) > τ. The bottom right panel provides a direct comparison among the three A0 values, with each value distinguished by a different colour.

In the text