© The Authors 2023

1. Introduction

The cosmic web (de Lapparent et al. 1986; Bond et al. 1996) is a network of structures present throughout the visible Universe. It is formed by the gravitational collapse of matter, which starts in the primordial Universe from tiny density perturbations (Zel’dovich 1970a; Zel’dovich 1970b). In this process, matter departs from underdense regions, which ultimately become the voids in the final matter distribution. Voids are surrounded by walls (two-dimensional structures). At the intersection of walls, filaments are found (one-dimensional, elongated structures). Matter flows inside walls at first and then in filaments (this is especially true at high redshifts) before finally reaching clusters. Galaxy clusters are located at the nodes of the cosmic web and are its most-studied and best-understood components.

The advent of large-scale numerical simulations, both N-body (Springel 2005) and hydrodynamical, such as Illustris (Vogelsberger et al. 2014), IllustrisTNG (Naiman et al. 2018; Marinacci et al. 2018; Springel et al. 2018; Nelson et al. 2018; Pillepich et al. 2018), MAGNETICUM (Hirschmann et al. 2014; Dolag 2015), EAGLE (Schaye et al. 2015), and HorizonAGN (Dubois et al. 2014), has allowed us to start investigating the properties of filaments from the theoretical side and to match observational results. Indeed, the filaments of the cosmic web and their impact on galaxies have been extensively investigated in numerical simulations (see e.g., Aragón-Calvo et al. 2010, 2019; Cautun et al. 2014; Laigle et al. 2015, 2018; Ganeshaiah Veena et al. 2018, 2019, 2021; Kraljic et al. 2019, 2020a; Galárraga-Espinosa et al. 2020, 2021, 2022; Gouin et al. 2021; Malavasi et al. 2022). Results range from the study of the shape, length, and dynamics of filaments, to the matter and volume repartition among the various cosmic web structures, to the connection of the filaments to the clusters and their impact on galaxy properties such as mass, star formation activity, and the direction of their angular momentum vector. These studies have found counterparts in similar ones coming from observations (Brouwer et al. 2016; Kuutma et al. 2017; Malavasi et al. 2017, 2020a; Laigle et al. 2018; Kraljic et al. 2018; Vulcani et al. 2019; Chen et al. 2019; Krolewski et al. 2019; Tanimura et al. 2020a,b; Bonjean et al. 2020; Rost et al. 2020; Welker et al. 2020).

Among numerical simulations, constrained ones (such as the CLUES project, Gottloeber et al. 2010, and its most recent development: CLONES, Sorce 2015, 2018; Sorce et al. 2016a, 2021, used e.g., in HESTIA, Libeskind et al. 2020) are gathering increasing importance. Although limited to the local Universe for the moment (at least when peculiar velocities are used as constraints; other techniques that involve the use of densities and whose description is outside the scope of this introduction allow much higher redshifts to be reached, see e.g., Ata et al. 2021), these kinds of simulations are based on real observations of galaxies and allow existing structures in a cosmic volume to be reproduced in a realistic way. They have proven increasingly important to interpret the observed properties of clusters in the cosmic web environment, including their evolution, contributing to unveil information normally unavailable through the analysis of clusters found in nonconstrained simulations (Sorce et al. 2020).

An important hot topic that has developed following recent advancements in the study of the cosmic web is the connection of the filaments to the galaxy clusters. Filaments have indeed been detected around clusters, both from the gas phase and the galaxy distribution. This ranges from the detection of the tips of the filaments connecting to the cluster in the X-rays (Eckert et al. 2015), to the detection of bridges of matter connecting pairs of close clusters (Akamatsu et al. 2017; Bonjean et al. 2018; Govoni et al. 2019; Reiprich et al. 2021; Biffi et al. 2022), to a systematic detection of filaments connected to clusters in large-scale simulations (Kuchner et al. 2020, 2021; Gouin et al. 2021), galaxy surveys (Sarron et al. 2019; Darragh Ford et al. 2019; Salerno et al. 2019, 2020), and individual cluster observations (Malavasi et al. 2020b; Einasto et al. 2020; Castignani et al. 2022).

It has been demonstrated that the number of filaments connected to a cluster has an impact on the cluster mass accretion process. More connected clusters are more massive (Codis et al. 2018; Sarron et al. 2019; Darragh Ford et al. 2019; Gouin et al. 2021) and in a more disturbed dynamical state (Gouin et al. 2021, 2022). Moreover, galaxy evolution in connected structures proceeded faster, as more connected galaxies (which we can consider to be the central galaxies of more connected clusters) are redder, less star-forming, and/or more massive (Darragh Ford et al. 2019; Kraljic et al. 2020b). Moreover, filaments can accrete matter across the virial shock of halos, bringing cold gas that fuels star formation directly on the forming halo (this is particularly true at high redshift, Dekel et al. 2009; Bennett & Sijacki 2020).

Among the structures for which the connection to the cosmic web has been investigated in the most detail is the Coma cluster of galaxies. The Coma cluster has been studied for several decades (Hubble & Humason 1931; Biviano 1998), with extensive data sets in all wavelengths, from the X-rays (Briel et al. 1992; Neumann et al. 2001, 2003) to the Sunyaev-Zel’dovich signal (Planck Collaboration Int. X 2013), radio (Brown & Rudnick 2011), and optical (den Brok et al. 2011; Adami et al. 2005). The connections of the Coma cluster to other structures of the surrounding cosmic web are already evident from the galaxy distribution and the position of the surrounding clusters (e.g., A1367, West 1998; Mahajan et al. 2018). In Malavasi et al. (2020b) we explored the connections of the Coma cluster to the filaments of the cosmic web. Filaments (Malavasi et al. 2020a) were identified in the Sloan Digital Sky Survey (SDSS, Abazajian et al. 2009; Strauss et al. 2002) with the Discrete Persistent Structure Extractor (DisPerSE, Sousbie 2011; Sousbie et al. 2011). We identified three filaments connected to the cluster, two of which are consistent with the position of the known features of the cluster itself. In particular, we identified a filament connecting to Coma on the west side of the cluster (as projected on the plane of the sky) whose position is broadly consistent with both the location of a shock feature visible in the X-ray map of the cluster (Neumann et al. 2003; Lyskova et al. 2019) as well as with a low-turbulence region visible in the entropy map of the cluster (Mirakhor & Walker 2020). Another secure filament was identified instead on the northeast side of the cluster (as projected on the plane of the sky), whose position is broadly consistent with the supposed direction of infall of the galaxy NGC 4839 and its associated group (Lyskova et al. 2019) and with a region of high gas velocity and low temperature and metallicity detected in XMM-Newton X-ray maps (Sanders et al. 2020).

The analysis performed in Malavasi et al. (2020b) and the supporting evidence exposed in the cited works all point to the possibility that accretion of matter on the Coma cluster has happened through the filaments and that it is ongoing. As observations can only offer a snapshot of the Coma cluster dynamics fixed in time, we decided to resort to numerical simulations to investigate the accretion of matter on the cluster coming from the filaments. The goal of this paper series is to provide further insight into the evolution and accretion of matter from the filaments onto the Coma cluster. In this first installment, we matched the simulated cosmic web to the observed one and we explore the spatial and velocity distributions of matter around the filaments and the simulated Coma cluster at z = 0¹.

This paper is structured as follows: we summarize the observational results of Malavasi et al. (2020b) in Sect. 2.1 and we introduce the constrained cosmological simulation that we used to reproduce the Coma system in Sect. 2.2. We describe the algorithm we used to detect the filaments in Sect. 4 and we extensively compare the observed and simulated cosmic web around Coma in Sect. 5. We then discuss the connectivity of the simulated Coma cluster in Sect. 6.1 and we analyze the dynamics of matter around the cluster and its connected filaments in Sect. 6.2. We discuss our results and draw our conclusions in Sects. 7 and 8, respectively.

Throughout this paper we use a Planck Collaboration XVI (2014) cosmology with H₀ = 100 ⋅ h = 67.8 km Mpc⁻¹ s⁻¹, Ω_m = 0.307, and Ω_Λ = 0.693 to be consistent with Sorce (2018) and Sorce et al. (2023). We note that this is different from the cosmology used in Malavasi et al. (2020b). We specify in the text whenever a different cosmology from Planck Collaboration XVI (2014) was used and we switch between the two whenever needed. Lengths are provided in units of Mpc h⁻¹. The only exceptions are for lengths that explicitly refer to values provided in Malavasi et al. (2020b), for example the ±75 Mpc radius used to explore the cosmic web around Coma in that work. In this case, we subsequently provide the corresponding value in Mpc h⁻¹ in parentheses for reference.

2. Data and simulations

In this section we describe our constrained cosmological simulation, which we use to reproduce the Coma cluster and analyze the physical properties of its surrounding LSS, and the observational data with which we compare our simulation. As this work heavily relies on what previously done in Malavasi et al. (2020a, b), our observational data essentially amount to what has been used in those works. We recapitulate and summarize the essential information in Sect. 2.1, where we also introduce further data sets we used for specific purposes in this analysis.

2.1. Observational data

The analysis by Malavasi et al. (2020b) identified the filaments connected to the Coma cluster among those present in the skeleton reconstructed in the Sloan Digital Sky Survey (Abazajian et al. 2009) with the Discrete Persistent Structure Extractor (DisPerSE, Sousbie 2011; Sousbie et al. 2011) and presented in Malavasi et al. (2020a). This work used a selection of 566 452 galaxies from the SDSS DR7 Main Galaxy Sample (MGS, Strauss et al. 2002) which have Petrosian r-band magnitude r_𝒫 ≤ 17.77, r-band half-light surface brightness μ₅₀ ≤ 24.5 mag arcsec⁻², secure spectroscopic redshift (ZWARNING = 0, ZCONFFINAL > 0.35, and ZFINAL > 0), and located in a contiguous region in the northern hemisphere. We refer the reader to Malavasi et al. (2020a, b) and Strauss et al. (2002) for further details. For the rest of this paper we refer to this selection of objects as the Legacy MGS. In this work, we have enhanced this data set with a measurement of the galaxy masses and K-band magnitudes. In the following, we mainly make use of K-band apparent magnitudes, which we compare with the same quantity available for the halos of our numerical simulation. We use this comparison to inform our selection of a galaxy sample from our halo population. As for galaxy masses, they are derived as shown in Appendix A and used only as a further confirmation of our simulated galaxy selection from our halo population and not to derive physical conclusions.

K-band apparent magnitudes for SDSS galaxies were obtained from the New York University Value-Added Galaxy Catalog (NYU-VAGC, Blanton et al. 2005; Adelman-McCarthy et al. 2008; Padmanabhan et al. 2008). This catalog provides K-band magnitudes for SDSS galaxies from the 2-Micron All Sky Survey (2MASS, Skrutskie et al. 2006). Starting from an initial sample of 2 506 754 galaxies, we select those belonging to the MGS (PRIMTARGET & 64, VAGC_SELECT & 4). The selection of galaxies in the NYU-VAGC adopted the same criteria as the MGS of Strauss et al. (2002), but used less stringent thresholds in terms of Petrosian magnitude cut, star-galaxy separation, fiber magnitude cut, and rejection of bright objects. We therefore reimplement the stricter thresholds of Strauss et al. (2002) and we select, as before, sources with good redshift measurements (ZWARNING = 0, z > 0). This leaves us with a sample of 685 813 galaxies which we match in position to the Legacy MGS with a tolerance of 0.5″. We check that the redshift of the matched objects is consistent and we identify 252 965 galaxies with measured K-band magnitude.

2.2. The numerical simulation

Constrained simulations are numerical simulations whose initial conditions are constrained from a sample of observed galaxies. These initial conditions, evolved with a N-body code, result in a simulation that reproduces the existing structures of the local Universe. The simulation for this work has been developed following the technique detailed in Sorce et al. (2016b) and Sorce (2018). The initial conditions are created from catalogs of distances to galaxies and groups (Tully et al. 2013; Sorce & Tempel 2017) which are converted to peculiar velocities following Sorce & Tempel (2018) and Sorce et al. (2016a) and bias-minimized (Sorce 2015). We stress that only the information on the peculiar velocities is used to constrain the initial conditions, without any use of the observed density field information.

The dark matter-only simulation used in this work was built on initial conditions containing 2048³ particles in a box of side 500 Mpc h⁻¹. The dark matter particle mass resolution offered by the simulation is of 10⁹ M_⊙ h⁻¹, which allows us to resolve dark matter halos of 10¹¹ M_⊙ h⁻¹ with 100 particles at z = 0. The simulation was run from z = 120 to z = 0 with the adaptive mesh refinement code RAMSES (Teyssier 2002). The best-achieved spatial resolution is of ∼1.89 kpc h⁻¹ thanks to subsequent refinements of the mesh (a level is refined if the total density in a cell is larger than that of a cell containing 8 dark matter particles).

We identify dark matter halos and subhalos with the HALOMAKER software (Aubert et al. 2004; Tweed et al. 2009), modified to operate with 2048³ (i.e., more than 2³¹) particles. Dark matter halos are identified in real space with the local maxima of the dark matter particle density field. Their boundary is defined as the point where the dark matter mass over-density is lower than 80 times the background density. The result is a sample of 3 666 018 dark matter halos in the simulation box.

2.3. Identifying the galaxy population in the numerical simulation

We identify a suitable subset of halos to mimic our galaxy population in two ways: with a mass cut and with an apparent K-band magnitude cut. In the rest of this paper we refer to these halos as “galaxies”. We do not make a distinction between whether a halo is a main halo or the subhalo of a main one. We use these halos to run DisPerSE to detect the cosmic web. Our mass-selected galaxy sample is composed of all 749 960 halos with mass in the range 10¹² ≤ M_vir(M_⊙)≤10¹³. To create our magnitude-selected galaxy sample we first computed the K-band absolute magnitude for every halo in the simulation by setting it proportional to the circular velocity of each halo (i.e., by inverting the Tully–Fisher relation, as obtained with fits to observed or simulated samples, see e.g., Eq. (7) of Mathis et al. 2002 and Eq. (3) of van den Bosch 2000):

(1)

The numerical values in this relation have been provided by Dolag (priv. comm.). The absolute K-band magnitude is then converted to apparent magnitude by means of the distance modulus: m_K = M_K + 25 + 5log(d). Based on our analysis of the K-band magnitude of SDSS galaxies we set a threshold at m_K = 14 and select all halos brighter than this limit (263 272) as galaxies.

Figure 1 shows the distribution of the apparent K-band magnitude for all halos in the simulation (as computed with Eq. (1)), those selected as galaxies and for galaxies in the Legacy MGS sample (by matching with the NYU-VAGC, as explained in Sect. 2.1). The distribution of Legacy MGS galaxies sharply drops around the value we selected as magnitude limit for our simulated galaxy sample, therefore providing supporting evidence for our choice. The distribution for all halos, however, peaks at much fainter magnitudes (due to the simulation being complete to the resolution limit for the halos), resulting in us selecting a small fraction of the available halos.

Fig. 1. Apparent K-band magnitude distributions for the halos in our simulation and for SDSS galaxies. The orange line refers to the Legacy MGS galaxies with apparent K-band magnitude from the NYU-VAGC sample. The light blue line refers to the mK distribution for all halos in our simulation (with mK computed with Eq. (1)). The green line shows the distribution for our magnitude-selected simulated galaxy sample. The vertical black dashed line shows our magnitude threshold of mK = 14.

We note that our magnitude selected sample of simulated galaxies includes halos with large values of M_vir (up to M_vir ∼ 10¹⁴ M_⊙). This is visible, for example, in the mass distribution of magnitude selected galaxies shown in Fig. A.1. As our galaxy samples include both main and sub-halos, we identify this population of galaxies residing in high-mass halos to be galaxies at the center of groups and clusters (e.g., BCGs). Indeed, deriving the mass distribution of only those magnitude selected simulated galaxies that are not main halos, it shifts to lower masses, more consistent with observations (not shown here). Given that BCGs are not eliminated from the Legacy MGS either, that the K-band apparent magnitude distributions of real and simulated galaxies are in agreement, and that DisPerSE is run only considering galaxy positions (without information about their mass or K-band magnitude required), we do not consider this to be a problem in our detection of the cosmic web.

Conversely, we note that the presence of an upper mass limit of 10¹³ M_⊙ in our mass-selected simulated galaxy sample largely excludes these BCGs in massive halos from the sample. To understand whether this may be a problem, we have performed our analysis also including these BCGs in massive halos in our mass-selected galaxy sample (in the number of 85 953 halos of mass M_vir > 10¹³ M_⊙ added to those already present in the mass-selected galaxy sample), that is by running DisPerSE on a sample of simulated galaxies selected to have M_vir ≥ 10¹² M_⊙. We do not find substantial differences with the conclusions reported in this paper with our current definition of a mass-selected galaxy sample (10¹² ≤ M_vir(M_⊙)≤10¹³) and we do not show the results of this additional test in the paper. In the following, we use a mass-selected galaxy sample with both a lower and upper mass limit to highlight how the cosmic web in the vicinity of Coma is relatively insensitive of our definition for a galaxy population (as long as it is reasonably similar to the Legacy MGS). We postpone to future papers the task of identifying an even more realistic simulated galaxy sample in our simulation through semi-analytical models.

2.4. Identifying the Coma cluster in the numerical simulation

We identify the halo in our simulation corresponding to the Coma cluster following Sorce (2018). First we select all halos more massive than 10^14.5 M_⊙ and whose distance on the Line-of-Sight (LoS) is within ±30% of the distance to the Coma cluster, that is |d_h − d_Coma|/d_Coma ≤ 0.3 where d_h is the distance to the considered halo and d_Coma is the distance to the Coma cluster. We derive d_Coma = 68.6 Mpc h⁻¹ by computing the comoving distance corresponding to z = 0.023 (which is the redshift of the Coma cluster adopted for the analysis of Malavasi et al. 2020b, see Table 1). Of the three candidates thus selected, we identify Coma with the most massive, which is also the one whose coordinates are the closest to the coordinates of the real Coma. In the rest of the paper we refer to this halo as “simulated Coma” and to the real cluster as “real Coma”. We note that there is still a residual difference in the supergalactic coordinates of the real and simulated Coma of ∼7 ÷ 14 Mpc h⁻¹. In the following, we work under the assumption that this residual difference in the coordinates of the simulated and real Coma is due to the simulation and not to the astrometric uncertainty in the coordinates of the real Coma, which we assume to be correct. We therefore apply a corresponding shift to the reference frame of the simulated box (i.e., to the coordinates of all halos including the one corresponding to Coma and of simulated galaxies) so that the position of the simulated and real Coma are identical. As this shift does not distort the box nor change its dimensions, it does not affect our conclusions. The physical parameters of the simulated Coma are therefore M_vir = 1.259 × 10¹⁵ M_⊙ h⁻¹, M₂₀₀ = 9.441 × 10¹⁴ M_⊙ h⁻¹, and r_vir = 1.75 Mpc h⁻¹. We report a value for the mass which is higher for the simulated Coma than for the real one used in Malavasi et al. (2020b, real Coma: M₂₀₀ = 3.57 × 10¹⁴ M_⊙, from Gavazzi et al. 2009, M_vir = 9.8 × 10¹⁴ M_⊙ h⁻¹, from Łokas & Mamon 2003). This difference could be due to the performance of the halo and sub-halo detection software, which may include in the mass determination for the halo corresponding to the simulated Coma unresolved substructure or particles which are actually unbound (see e.g., Knebe et al. 2011 for a comparison of halo detection algorithms), or to uncertainties in the Coma mass determination from observations. The mass of our simulated Coma is in better agreement with (albeit slightly lower than) the mass for the real Coma cluster reported by Tully (2015, M_vir = 1.7 × 10¹⁵ M_⊙ h⁻¹). In our case, the difference between the mass of the real and simulated Coma is relevant when addressing the connectivity-mass relation in Sect. 6.1. Indeed this difference in the mass prevents us from decisively concluding whether the connectivity measurement for the simulated Coma is in agreement with the observed value. However, including the uncertainty on the observed mass measurement, we can conclude that the connectivity of simulated Coma is broadly in line with the expected trends from observations and numerical simulations in the literature (see Sect. 6.1).

3. Inclusion of the redshift space distortions

In order to perform a more realistic comparison between the cosmic web we detect in simulations and the one we detect in observations, we have introduced redshift space distortions (Finger of God effect, FoG) in the position of our simulated galaxies. The FoG effect is a distortion of galaxy positions along the LoS in the vicinity of massive structures due to the fact that galaxy peculiar velocities introduce an error in the measurement of their redshift (Kaiser 1987). Due to this effect, clusters appear as elongated structures (and not spherical ones) in the LoS direction, and they can be easily mistaken for a filament by algorithms such as DisPerSE. As the filaments detected around Coma in observations are recovered with the FoG effect being present in the Legacy MGS (see the discussion of this point in Malavasi et al. 2020a,b) while this effect is absent in our constrained simulation, we introduced it in order to make the comparison more realistic.

We introduced the FoG effect in two ways, both of which follow a general procedure: first of all we computed peculiar velocities for all the galaxies in our simulation, then we converted peculiar velocities to redshift, and finally we added the peculiar redshift to the cosmological one for each galaxy. The peculiar velocity (v) is provided in cartesian supergalactic coordinates for each galaxy as a direct output of our simulation (see also Sect. 6.2). We convert it into supergalactic spherical coordinates, which provides us with a proper motion for each galaxy on the plane of the sky (which we ignore) and a radial velocity (to which we refer as v_pec in the following). This velocity is along the LoS (negative or positive depending whether a galaxy is approaching or receding). While this radial velocity is directly provided by the simulation, we also obtained v_pec in a different way for completeness and to test a different method: for each galaxy which is a main halo (i.e., only for halos which are not substructures of other main halos) we derive the velocity dispersion from the virial mass, using the formula (see also Malavasi et al. 2020a). For each sub-halo (i.e., only for halos which are substructures of other main halos), we then randomly extracted a value of v_pec from a Gaussian distribution centered on zero and with a dispersion of σ_v. As this approach of determining v_pec provided results consistent with using the peculiar velocities directly provided by the simulation, we do not discuss it further and we rather focus on the first method.

Once a value of v_pec was determined, we converted it to redshift using the formula z_pec = v_pec/c, with c being the speed of light. We then added this peculiar redshift to the true redshift of each galaxy (as derived from its distance, obtained from its supergalactic coordinates) as described in Tully et al. (2013) and Sorce et al. (2023): z_obs = (1 + z_true)(1 + z_pec)−1. This observed redshift (z_obs) was then used to compute a new distance for each galaxy, which was then used to compute new supergalactic cartesian coordinates (assuming the angular position of each galaxy on the plane of the sky to remain unchanged). These new positions were then used to run DisPerSE in the same way we did for the galaxy sample where the FoG effect was not introduced (see Sect. 4). We visually checked that the prominence of the FoG for the simulated Coma cluster constructed in this way is consistent with the FoG for the observed Coma present in both the Legacy MGS and the sample of galaxies of Tully et al. (2013).

4. Detecting the cosmic web

We detect the cosmic web in the simulation by applying the same cosmic web detection algorithm used in Malavasi et al. (2020a,b) to our population of simulated galaxies.

The Discrete Persistent Structure Extractor (DisPerSE, Sousbie 2011; Sousbie et al. 2011) extracts the filaments of the cosmic web by means of their topology. It works on discrete sets of points: galaxy distributions from surveys, and halo distributions in simulations. For the application of DisPerSE to the Legacy MGS we refer the reader to Malavasi et al. (2020a). In this work, we apply DisPerSE to our galaxy populations selected through the mass and K-band magnitude cuts. The starting point is a measurement of the galaxy density field, obtained through the Delaunay Tessellation Field Estimator algorithm (Schaap & van de Weygaert 2000; van de Weygaert & Schaap 2009). While not necessary for the filament extraction, the galaxy density field can be smoothed by recursively averaging the density measured at the position of each galaxy with the density of all galaxies connected to it by edges of the tetrahedrons of the Delaunay tessellation. For our analysis we work with either an un-smoothed density field (referred to as “SD0”) or a density field smoothed only once (“SD1”). DisPerSE then computes the discrete gradient of the density field and identifies points (called critical points) where the gradient is null (this is done by applying the discrete Morse theory). These critical points can be maxima, minima, and two type of saddles (1-saddles, local density minima bound to walls, and 2-saddles, local density minima bound to filaments). Filaments consistently connect maxima and 2-saddles together through lines of constant gradient (ridges of the density field). The DisPerSE algorithm implements a selection based on the persistence of the filaments (tied to the density contrast of the critical points at their extrema and measured based on the persistent homology theory). The distribution of persistence values for the pairs of critical points detected in the data is compared to the same distribution for pairs of critical points detected in a Gaussian random field (which models the set of filaments and critical points which would be obtained if only noise and no real data were present). All filaments connected to critical points closer (in terms of their persistence value) to the noise distribution than a certain number of standard deviations are removed via topological simplification. In this work we consider persistence thresholds for the elimination of spurious filaments of 2σ and 3σ. The skeleton (i.e., the set of filaments and critical points) thus constructed is then post-processed: artificial critical points (called bifurcations) are inserted where two or more filaments intersect without a maximum being present. The filaments are composed of small segments (the size of the edges of the tetrahedrons of the Delaunay tessellation). The filament shape is smoothed by recursively averaging the positions of the extremes of each of these segments with the positions of the extremes of the segments directly attached to it (without modifying the positions of the critical points at the extremes of the filaments). In this work, this smoothing is done only once.

5. Comparing the cosmic web in simulated and observed data

The approach to detect the cosmic web used in this work is the same as the one used in Malavasi et al. (2020b). As a reference, from observations we select the skeleton extracted in the Legacy MGS with the combination of parameters SD1 (1 smoothing cycle) and a persistence threshold of 3σ (see the top panel of Fig. 2 of Malavasi et al. 2020b). The first step in our analysis is represented by finding among the parameter combinations explored in this work (DisPerSE parameters and galaxy mass vs. magnitude selection) those which yield the best match between the skeleton detected in the constrained simulation and the one detected in observations. As a perfect match between the simulation and observations is very unlikely due to the residual cosmic variance of the former², we aim at excluding some of the parameter combinations by means of this comparison. In the rest of this section we proceed by first comparing visually the simulated and observed skeleton, then by comparing some global properties of the filaments and critical points and finally by comparing some simulated and observed properties of the Coma cluster and its filaments. We begin by comparing the observed skeleton with the one derived in the constrained simulation without the inclusion of the FoG effect. We provide a comparison with the inclusion of the FoG effect in Sect. 5.3.

5.1. Visual comparison of the filament configuration

We start by identifying the critical points associated with Coma, in a similar way to what we did in Malavasi et al. (2020b). For each combination of DisPerSE parameters, we identify all the critical points within 1 r_vir of the simulated Coma. Table 1 reports our findings, broken down by critical point type, as well as our results from the parameter combination SD1−3σ from Table 2 of Malavasi et al. (2020b), for reference.

Table 1 shows how the simulated Coma cluster is consistently identified with a maximum regardless of the combination of DisPerSE parameters or galaxy selection method. In the case of the mass selection for the simulated galaxy population, an additional bifurcation point is found inside r_vir, similar to what happens in the case of Malavasi et al. (2020b) in some of the DisPerSE parameter combinations. In the case of magnitude-selected galaxies, the amount of critical points detected within the r_vir of Coma is larger. When no smoothing and a low persistence threshold are chosen (SD0−2σ) both 1-saddles and 2-saddles (local minima bound to filaments and sheets) are also found within the r_vir of Coma.

Starting from this selection of critical points we recursively find all the filaments attached to them and those attached to the critical points at their ends, iterating the search for filaments up to a radius of 75 comoving Mpc (50.8 Mpc h⁻¹), as it was done in Malavasi et al. (2020b). In agreement with that work, we call the filaments connected to the critical points within the virial radius of Coma “first generation” and those connected to the critical points at the other end of first generation filaments (and not within the virial radius of Coma) “second generation”. “Third generation” filaments are then connected to the critical points at the other end of second generation filaments (and which are not connected to first generation filaments), “fourth generation” filaments are connected to third generation, etc. We show first and second generation real and simulated filaments in Fig. 2. This figure shows a supergalactic coordinate projection (SGZ-SGY) in which first and second generation filaments connected to the real and simulated Coma are compared. In this figure, the redshift direction is aligned with the SGY axis, while the SGZ axis is located in the plane of the sky. The northeast and west filaments from the observations (in blue) are highlighted in this plot, clearly visible intersecting at the cluster position in an almost straight axis (in the proximity of the cluster). Perpendicular to this axis is an elongated filament which is located roughly along the Finger of God (FoG) due to the Coma cluster. While such a filament is most probably unphysical and the result of the FoG distortion, in Malavasi et al. (2020b) we performed tests to check whether the fact of not having compressed the FoGs had any bearing on our results related to the connection of the filaments to the cluster, while in Malavasi et al. (2020a) we demonstrated that the filament reconstruction in the SDSS, without compressing the FoG, did not statistically affect the analysis we intended to perform. In this work we aim to identify the match between observations and simulations on the basis of the positions of the northeast and west filaments alone, whose detection is more secure as outlined in Malavasi et al. (2020b). Given the observed filaments’ orientation (i.e., the fact that the northeast and west filaments are along an almost-straight axis along the SGZ direction), we believe this is the best projection to identify possible matches between observed and simulated filaments. We provide an extra projection in Appendix B along with a short discussion.

Fig. 2. Filaments in the constrained simulation compared to filaments in observations. The left column shows the case of a mass selected galaxy sample in the constrained simulation, the right column shows the case of a magnitude selected galaxy sample (see text). In each panel, simulated galaxies are shown as gray points in a slice of thickness 50 Mpc h−1 centered on Coma. Blue lines are filaments from Malavasi et al. (2020b) in supergalactic coordinates, the red square marks the position of the Coma cluster. Dark red lines are filaments obtained in the simulation with a variety of DisPerSE parameter combinations, namely: SD0−2σ (top row), SD0−3σ (second row), SD1−2σ (third row), SD1−3σ (bottom row). The observed northeast and west filaments are marked for reference. Circle and triangles mark the positions of critical points (in observations and simulations, respectively) and are color-coded according to their type (red: maxima, green: type 2 saddles, brown: bifurcations). This figure shows the SGZ-SGY projection.

Several of the DisPerSE parameter combinations offer filament configurations around Coma which are similar to what observed in Malavasi et al. (2020b). In particular, SD0−2σ in the case of magnitude selected galaxies and SD1−2σ and SD1−3σ in the case of mass selected galaxies all show prominent northeast and west filaments aligned in an axis at around the same location as the observed ones. In the case of magnitude selected galaxies, the west filament is present in all DisPerSE parameter configurations, while the northeast filament is present also in the case of SD1−2σ, although shorter in length. Based on these considerations, it seems that in the case of magnitude selected galaxies, the most important DisPerSE parameter is the persistence threshold (which has to be on the order of 2σ), while in the case of the mass selected galaxies smoothing seems the most important parameter.

5.2. Statistical comparison of the skeleton properties

We then check how the selected skeletons compare in terms of their general properties. We stress that this comparison is less constraining in terms of determining which simulated skeleton is similar to the observed one as several statistical quantities of the skeleton (e.g., the density of the critical points and the length of the filaments) strongly depend on the density of tracers (higher in the simulations and dependent on redshift for the observations) and on the region in which the skeleton is extracted. However, such comparisons (in particular for the connectivity) may still be interesting to understand better the differences between a simulated and observed skeleton and to have a general view of the properties of filaments explored.

In Fig. 3 we begin by showing the distributions of filament lengths for the observed and simulated skeletons. To derive this and the following distributions, both for observations and simulations, we have focused only on filaments and critical points within a slice of thickness ±75 Mpc (50.8 Mpc h⁻¹) from the location of the (simulated or observed) Coma cluster. We do so in order to limit in redshift the filaments from the observations for which we derive properties (and ensure that they are extracted in a region where the density of tracers is more or less constant) and to limit in density the region we explore (the environment of a rich cluster is different from the average density field and the filaments there may not be representative of the average filament in the field we explore).

Fig. 3. Filament length distribution. The black bold distribution refers to observed filaments from Malavasi et al. (2020b) with SD1−3σ. The colored lines refer to simulated filaments with parameter combinations SD0−2σ and SD0−3σ (for magnitude selected galaxies, in purple and red, respectively), and SD1−2σ and SD1−3σ (for mass selected galaxies, in orange and green, respectively). Dashed gray lines represent the rest of the parameter combinations and provide an idea of the range of lengths explored. Only filaments within a slice of thickness ±75 Mpc (50.8 Mpc h−1) from the location of the (simulated or observed) Coma cluster were used to derive the distributions.

The observed and simulated skeletons have largely similar length distributions, showing that although the tracer density may differ, the filaments have comparable lengths. In this and the following figures, we highlight with colored lines the distributions corresponding to the four simulated skeletons selected as those that best match the observations via visual comparison in the previous step (i.e., the combinations SD0−2σ and SD0−3σ for magnitude selected galaxies and SD1−2σ and SD1−3σ for mass selected galaxies). For completeness, we report in gray the distributions corresponding to the simulated skeletons obtained with other parameter combinations not selected via visual comparison in Fig. 2. We do so to show the full ranges of skeleton properties explored by our simulated skeletons while highlighting the DisPerSE parameter combinations for the simulated skeletons that we use in our analysis. Figure 3 shows how these distributions are actually those with the largest systematic difference in terms of filament length distributions from the observed skeleton. Simulated skeletons with lower values of smoothing (SD0) are characterized by systematically shorter filaments than observations, while skeletons with larger values of smoothing (SD1) have filaments which are on average longer³. Our selection therefore allows to explore two extrema in the relation between observed and simulated filaments (at least with respect to their length distribution) and to obtain a better idea of the general behavior.

We then focus on the connectivity of maxima and bifurcations. In this case, we define the critical point connectivity as the number of filaments connected to a given maximum or bifurcation. We explore the connectivity distribution for filaments in the same distance range of ±75 Mpc (50.8 Mpc h⁻¹) around Coma as for the length distribution. The connectivity distribution for maxima and bifurcations combined is shown in Fig. 4. In the following of this paper we refer to maxima and bifurcations together as nodes. We stress that this term is not used in a topological sense (unlike maxima and bifurcations), but in an astrophysical one to refer to dense structures of the cosmic web.

Fig. 4. Node connectivity distribution. The black bold distribution refers to observed nodes from Malavasi et al. (2020b) with SD1−3σ. The colored lines refer to simulated nodes with parameter combinations SD0−2σ and SD0−3σ (for magnitude selected galaxies, in purple and red, respectively), and SD1−2σ and SD1−3σ (for mass selected galaxies, in orange and green, respectively). Dashed gray lines represent the rest of the parameter combinations and provide an idea of the range of connectivities explored. Only nodes within a slice of thickness ±75 Mpc (50.8 Mpc h−1) from the location of the (simulated or observed) Coma cluster were used to derive the distributions.

This figure shows that the majority of points in the vicinity of Coma have a connectivity of three and this is consistent between observations and simulations. For reference, for the Coma cluster, Malavasi et al. (2020b) reported a connectivity of 2.5. All the DisPerSE parameter combinations we explored are in agreement with the observed connectivity distribution, having a peak at a value of κ = 3. This figure shows the same situation as the length distribution: the parameter combinations we selected based on visual identification (highlighted in color in the figure) are at the outskirts of the range covered by the various connectivity distributions from the simulations. Parameter combinations with a larger smoothing (SD1) have a larger fraction of nodes with a connectivity of three (the combination SD1−3σ for mass selected galaxies is the closest to the observed distribution), while combinations with a lower smoothing have a lower fraction of nodes with a connectivity of three and a larger fraction of nodes with a connectivity of two.

We then further refine our selection of the nodes for the comparison of statistical properties of the skeleton. First, we select only nodes in a slice of thickness ±75 Mpc (50.8 Mpc h⁻¹) from the location of the (simulated or observed) Coma cluster and with an overdensity of log(1 + δ) = log(1 + δ_Coma)±0.5, where log(1 + δ_Coma) = 1.64 is the average overdensity of the critical points found within the virial radius of the cluster in observations. We show the connectivity distribution for these nodes in Fig. 5. In this case, the peak at κ = 3 is less evident and the low connectivity values are more populated. We find that in this case, low persistence (2σ) parameter combinations are closer to observations.

Fig. 5. Node connectivity distribution. The black bold distribution refers to observed nodes from Malavasi et al. (2020b) with SD1−3σ. The colored lines refer to simulated nodes with parameter combinations SD0−2σ and SD0−3σ (for magnitude selected galaxies, in purple and red, respectively), and SD1−2σ and SD1−3σ (for mass selected galaxies, in orange and green, respectively). Dashed gray lines represent the rest of the parameter combinations and provide an idea of the range of connectivities explored. Only nodes within a slice of thickness ±75 Mpc (50.8 Mpc h−1) from the location of the (simulated or observed) Coma cluster and with overdensity log(1 + δ) = log(1 + δComa)±0.5 (see text) were used to derive the distributions.

As a last test, we select nodes in a slice of thickness ±75 Mpc (50.8 Mpc h⁻¹) from the location of the (simulated or observed) Coma cluster and with a connectivity value of κ = 2 ÷ 3, that is in line with the connectivity of the Coma cluster from observations. We measure the 1 + δ distribution for maxima and bifurcations in this connectivity range, which we show in Fig. 6. This figure shows that nodes in observations are distributed at higher densities than nodes in simulations. While no density distribution from the simulated skeletons is fully in the range of densities covered by observations, skeletons with large smoothing levels (SD1) are closer to the observed skeleton. In this case too, the highlighted parameter combinations are found at the extrema of the range of behaviors shown by the simulated skeletons.

Fig. 6. Node density distribution. The black bold distribution refers to observed nodes from Malavasi et al. (2020b) with SD1−3σ. The colored lines refer to simulated nodes with parameter combinations SD0−2σ and SD0−3σ (for magnitude selected galaxies, in purple and red, respectively), and SD1−2σ and SD1−3σ (for mass selected galaxies, in orange and green, respectively). Dashed gray lines represent the rest of the parameter combinations and provide an idea of the range of connectivities explored. Only nodes within a slice of thickness ±75 Mpc (50.8 Mpc h−1) from the location of the (simulated or observed) Coma cluster and with κ = 2 ÷ 3 are used to derive the distributions.

5.3. The cosmic web around Coma with the FoG effect

In Table 2 we report the version of Table 1 we obtain when we introduce FoG distortions in our galaxy population. We stress that the FoG distortions are introduced using all galaxies in the box, therefore for all clusters, not just Coma. This distorts the entire skeleton we recover in the box and not just the cosmic web around Coma. The fact that we still detect a realistic cosmic web connected to the Coma cluster is thus a success of our method. One of the differences between the cosmic web detected around Coma with and without FoG distortion is that we had to increase the radius of the sphere centered on Coma inside which we match critical points to the cluster. While this radius was set to r_vir for the case without FoG effect for both the mass selected and the magnitude selected sample, we resolve to set this search radius to 1.5 × r_vir (for the combination SD0−2σ), 2.5 × r_vir (for the combination SD0−3σ), and 3.5 × r_vir (for the combinations SD1−2σ and 3σ) for mass selected galaxies. This is due to the fact that no critical points can be found within one virial radius of our simulated Coma in the case of mass selected galaxies. The value of the search radius is still set to r_vir for magnitude selected galaxies⁴.

Figure 7 shows the new skeleton detected in the vicinity of the Coma cluster when the FoG effect is included in the analysis of our simulation.

Fig. 7. Filaments in the constrained simulation compared to filaments in observations when the FoG effect is introduced in our simulation. The left column shows the case of a mass selected galaxy sample in the constrained simulation, the right column shows the case of a magnitude selected galaxy sample (see text). In each panel, simulated galaxies are shown as gray points in a slice of thickness 50 Mpc h−1 centered on Coma. Blue lines are filaments from Malavasi et al. (2020b) in supergalactic coordinates, the red square marks the position of the Coma cluster. Dark red lines are filaments obtained in the simulation with a variety of DisPerSE parameter combinations, namely: SD0−2σ (top row), SD0−3σ (second row), SD1−2σ (third row), SD1−3σ (bottom row). The observed northeast and west filaments are marked for reference. Circle and triangles mark the positions of critical points (in observations and simulations, respectively) and are color-coded according to their type (red: maxima, green: type 2 saddles, brown: bifurcations). This figure shows the SGZ-SGY projection.

This figure shows that the cosmic web we detect introducing FoGs is indeed realistic and close to the observed one. In the case of magnitude selected galaxies, for all smoothing levels and persistence threshold we recover a filament to the west and one to the northeast, more prominent for higher smoothing levels. The presence of these filaments is consistent with the case in which no FoGs are introduced. Moreover, filaments are now detected along the FoG direction (which were not present in the case without FoGs), confirming our assumption that this filament is due to the FoG distortion effect.

In the case of mass selected galaxies, the situation is more complex. The large amount of virial radii that we need to reach in order to match critical points to the simulated Coma cluster indicates that the reconstruction of the cosmic web around the simulated Coma in this case is more uncertain. In the case of mass selected galaxies we recover a filament along the FoG direction for Coma only in the case SD0−2σ and only hints of this in the other smoothing level and persistence combinations. In all combinations we do, however, recover both the west and the northeast filaments. Similar conclusions (albeit less evident) are reached also when FoGs are introduced using random extracted velocities for galaxies.

6. Results

In the previous sections we have proved that we can reconstruct a realistic simulated skeleton in the vicinity of the Coma cluster using our constrained simulation and we have tested that the reconstruction is robust against the fact that redshift space distortions are present in our observed data. In the following we use the simulated skeleton derived without the introduction of the FoG effect to perform our analysis. We do so as we consider this skeleton to be more reliable, all the while we make sure that our results do not change if we use the skeleton with the inclusion of the FoG effect.

Based on our previous analysis we identify in the combination of DisPerSE parameters SD0−2σ, galaxy selection based on apparent K-band magnitude and skeleton computed without introducing redshift space distortions our reference simulated filaments around Coma. In the following we use this reference skeleton for visual representation, while we include all DisPerSE parameter combinations and galaxy selection criteria (namely SD0−2σ and SD0−3σ for magnitude selected galaxies and SD1−2σ and SD1−3σ for mass selected galaxies) in our statistical measurements.

In the remainder of the paper we analyze the properties of the cosmic web connected to the simulated Coma cluster in comparison with observations. We start by studying the connectivity of the cluster (i.e., the number of connected filaments) and then move on to the velocity distribution of the halos around the cluster and the filaments.

6.1. The connectivity of the Coma cluster

In Malavasi et al. (2020b) we measured the connectivity of the real Coma cluster by counting the number of filaments crossing a sphere of radius 1.5 × r_vir centered on the cluster. We detected a connectivity κ = 2 ÷ 3 (median connectivity κ = 2.5). This is in line with the average connectivity of the nodes in both our simulated box and in observations (see e.g., Figs. 4 and 5).

We have measured the connectivity for the simulated Coma in the same way. Starting from the list of critical points associated with the cluster (Table 1), we have counted all the filaments with one of those critical points as an extreme and the other one located outside a sphere of 1.5 × r_vir. We report a median connectivity for simulated Coma of κ = 3.0 (with the connectivity in the range κ = 3 ÷ 4 depending on the DisPerSE parameter combination and galaxy selection criterion). This is in good agreement with our observations, with the slightly larger number of connected filaments possibly due to the larger amount of tracers and better definition of the filaments in the constrained simulation. In Fig. 8 we show a rendition of Fig. 5 of Malavasi et al. (2020b), which shows the connectivity of the real Coma, simulated Coma, and other cluster samples from observations and simulations as a function of mass. The other theoretical and observed samples of clusters for which the connectivity has been measured are the AMASCFI clusters (Sarron et al. 2018, with the connectivity measured by Sarron et al. 2019), the groups detected in COSMOS by Darragh Ford et al. (2019), the connectivity-mass relation identified in N-body simulations by Aragón-Calvo et al. (2010) and the theoretical one derived by Codis et al. (2018) renormalized to pass through the point corresponding to the real Coma as done in Malavasi et al. (2020b). To these we have added the connectivity-mass trends derived by Gouin et al. (2021) in the IllustrisTNG simulation (Naiman et al. 2018; Marinacci et al. 2018; Springel et al. 2018; Nelson et al. 2018; Pillepich et al. 2018), for relaxed old clusters (RO), and unrelaxed old and young ones (UO and UY, respectively, see Gouin et al. 2021 for the distinction between relaxed and unrelaxed, old and young clusters). As for the measurements related to the real and simulated Coma, they are placed on this plot based on the measured median connectivity (with error bars corresponding to the measured connectivity range) and according to their M₂₀₀ value (measured by Gavazzi et al. 2009 in the case of the real Coma). We have converted the M₂₀₀ value for the simulated Coma to the same cosmology we adopted in Malavasi et al. (2020b). We report the error bars on the measurement for M₂₀₀ in the case of the real Coma to show that, although the measurement of the mass is lower than what we find in the simulation, there could still be room for agreement thanks to the uncertainty on the observed value (see also the discussion in Sect. 2.4).

Fig. 8. Connectivity of the real and simulated Coma clusters as a function of mass. The cyan square refers to the real Coma cluster as analyzed in Malavasi et al. (2020b), while the red square corresponds to the simulated Coma cluster as detected in this work. Black crosses are the observed connectivity values from the AMASCFI Clusters (Sarron et al. 2018, 2019), green triangles are the observed connectivity values for the groups in COSMOS (Darragh Ford et al. 2019). Gray lines and points are connectivity measurements derived in the IllustrisTNG simulation (Gouin et al. 2021) for clusters and groups which are unrelaxed and old (UO, solid line), relaxed and old (RO, dotted line), and unrelaxed and young (UY, dashed line). The purple line and shaded region are the relation from the numerical N-body simulations of Aragón-Calvo et al. (2010), and its corresponding 1σ uncertainty, while the solid orange line is the theoretical relation of Codis et al. (2018) renormalized to pass through our measurement for the real Coma cluster, so as to provide a better comparison of trends. We note that the mass on the x-axis of the plot is M200 for the real and simulated Coma, the measurements by Darragh Ford et al. (2019), and those of Sarron et al. (2019), while it is a value close to Mvir for the Aragón-Calvo et al. (2010) relation.

Even if the caveats on the connectivity measurements as a function of mass for observed and simulated clusters from the literature still apply to this situation⁵, we find agreement between the connectivity value at the mass of the simulated Coma and the values from the literature reported in this plot. The connectivity we measured for the simulated Coma is at the lower edge of the general trend identified by the literature values, indicating moderate agreement or a slightly lower value. As this is the same conclusion we found for the real Coma, investigating why these massive clusters have lower connectivity values than what expected given their mass, whether this is a real effect or due to the way in which connectivity is measured, and how this picture evolves with redshift strongly calls for a deeper investigation of the connectivity-relation from low to high redshift and from the group to the cluster (and potentially super-cluster) regime, which will be the subject of a future work.

6.2. Accretion on the Coma cluster through the cosmic web

We explore the dynamics of matter in the vicinity of the simulated Coma cluster and its connected filaments. To do this we make use of the velocities of the dark matter halos in the simulation box. Each halo is characterized by a velocity vector with components v = (v_x, v_y, v_z). For each halo (including those selected as galaxies), we subtract the velocity vector of the halo identified as the simulated Coma: v_h = v − v_Coma. We then compute the dot product of the velocity for each halo and the difference in position between each halo and the position of the simulated Coma (r_h = r − r_Coma). The dot product

(2)

where r_h is the norm of the vector r_h, provides the module of the velocity of each halo, projected along the line connecting the simulated Coma cluster with that halo, in a reference frame where Coma is at the origin and at rest. We then construct the 3D components of such an array in the form:

(3)

We show maps of these quantities in relation to the filament positions in the surroundings of Coma in Fig. 9. This figure shows three projections of the 3D system formed by Coma and its surrounding filaments. In each projection, the color map shows a 2D histogram in which in each bin the average v_rad, C is computed. Black arrows, instead, show the average direction computed only from the components relevant to the considered slice (e.g., in the x − y projection each black arrow is constructed from the average in the 2D histogram bin of v_{rad, C, x} and v_{rad, C, y}). This gives an idea of the module and direction of the velocity of the halos in the simulation box in proximity of the simulated Coma and its filaments. Only the case of magnitude selected galaxies and SD0−2σ is shown for reference in this figure, while all four parameter combinations described above as the closest to observations are considered when the infall on Coma from the filaments is quantified in the following.

Fig. 9. Maps of halo velocities around the simulated Coma cluster. The gray square at the center marks the position of the Coma cluster. Orange lines are first, second, and third generation filaments (the case SD0−2σ for magnitude selected galaxies is reported here as an example, see brown lines in Fig. 2 and text). The top left panel shows the x − y projection, the top right panel shows the x − z projection, and the bottom left panel shows the z − y projection (shown in the first panel on the right column of Fig. 2). In each slice, the color map shows a 2D histogram in which in each bin the average vrad, C is measured (red stands for halos moving away from Coma, blue for halos moving toward Coma, white pixels stand for empty bins). Black arrows are constructed from the average in each bin of the vrad, C components relative to the projection shown (i.e., vrad, C, x and vrad, C, y in the top left panel, vrad, C, x and vrad, C, z in the top right panel, vrad, C, z and vrad, C, y in the bottom left panel). The thickness of the considered slice is of 20 Mpc h−1, all halos in the slice are considered, regardless of mass or whether they are selected as galaxies. The labels mark the position of the simulated northeast and west filaments.

This figure shows how there is a significant number of halos moving toward Coma from the large-scale structure. The velocity of the halos seems to correlate fairly well with the filament position. This is particularly evident for example in the x − y and z − y projections, especially for the filaments to the northeast and the west of the cluster in the latter.

We try to quantify the relation between halo velocity (whether infalling or outflowing) and filament direction. In Fig. 10 we report the value of v_rad, C as a function of halo distance from Coma. We consider all the halos which are outside of the simulated Coma virial radius (so as to reduce the influence of the cluster on the galaxy motions) and within 10 Mpc (6.774 Mpc h⁻¹) from the cluster center. We chose this limit as we see from Fig. 9 that even third generation filaments do not generally extend beyond such a radius. The general distribution of halos in this region shows that infall velocity tends to increase with decreasing distance from the simulated Coma center, except in the close vicinity of the cluster, where galaxies tend to have a larger scatter around zero, with some having large outflow velocities. This could be due to galaxies that have already gone through the cluster once and are coming out the other side of the cluster with a positive velocity or to galaxies that are close to being virialized and that are within the zero-velocity surface of the cluster. At distances from Coma larger than 2.5 Mpc h⁻¹, velocities are generally negative, implying infall on the cluster. In the same figure, we also show the velocity that a body would reach as a function of distance from Coma by starting at 10 Mpc (6.774 Mpc h⁻¹), with a velocity equal to the value of the last bin of the running mean as derived for the total halo population, and being subject only to the gravitational attraction of the Coma cluster. The velocity distribution for the total halo population closely follows this relation (except close to the virial radius of the cluster), implying that no additional effect besides the overall gravitational attraction from the cluster is at play here. We derive the relation between infall velocity and distance from the simulated cluster also only for those halos that are closer to the axes of the filaments than a certain threshold. For each halo, we measure the distance between it and the axis of the closest filament (d_fil) by using the common formula for the distance between a point and a line (where the line in this situation is the closest segment of the closest filament to the halo, using the DisPerSE formalism)⁶. As a threshold to consider a halo close enough from the axis of a filament to be reported on Fig. 10, we chose the 25th percentile of the d_fil distribution.

Fig. 10. Halo velocity as a function of distance from Coma. Each panel shows the halo velocity vrad, C as a function of the distance between each halo and the center of the simulated Coma cluster dh − Coma = |rh|. Only halos with rvir ≤ dh − Coma ≤ 10 Mpc (6.774 Mpc h−1) have been considered. Blue squares refer to all halos within this distance range, orange triangles are halos with dfil smaller than the 25th percentile of the dfil distribution. Solid lines refer to the running means of the distributions (cyan: all halos, red: halos close to filaments), error bars encompass the region between the 16th and the 84th percentile of the distributions. Different panels refer to different combinations of DisPerSE parameters and galaxy selections. In each panel, the black line shows the velocity that a body would reach as a function of distance from Coma by starting at 10 Mpc (6.774 Mpc h−1), with a velocity equal to the value of the last bin of the running mean as derived for the total halo population, and being subject only to the gravitational attraction of the Coma cluster.

There is no obvious difference in the velocity distribution as a function of halo distance from Coma for halos in filaments compared to the general halo population. In both the case of filaments extracted with magnitude-selected and mass-selected galaxies, even the averages of the populations are largely overlapping, showing no difference. Only in the case of mass selected galaxies, a few of the bins located at distances from Coma of ∼4 ÷ 6 Mpc show more negative values, suggesting larger inflow velocities. However, from this figure we are bound to conclude that the amplitude of the velocity of halos is somewhat uncorrelated from their position with respect to filaments, and it is mostly sensitive to the large-scale inflow of matter on the Coma cluster.

We therefore explored the relation between the distance of the halos from the axes of the filaments d_fil and the angle between the direction of the halo projected velocity (v_rad, C = (v_{rad, C, x}, v_{rad, C, y}, v_{rad, C, z}) and the local direction of the closest filament (a_fil). The local direction of the closest filament to a given halo is consistently chosen to be positive when pointing away from the Coma cluster, so as to be consistent with our definition of the projected velocity. In this reference frame, halos with velocities infalling toward Coma have angles closer to 180 deg, while halos outflowing from the cluster have an angle closer to 0 deg.

Figure 11 shows the relation between a_fil and d_fil. The majority of the halos have angles close to 180 deg, as it is expected given the significant infall of material on the cluster. However, there is also a trend visible of a_fil becoming closer to 180 deg with decreasing d_fil. This provides an indication of the fact that the flux of matter infalling on Coma becomes more collimated around the location of the filaments. This is the first time that such a result is reported. Moreover, while this effect is small at low redshift, we expect it to become stronger at higher redshift when the accretion on the cluster from the filaments was higher. We will investigate this in a future work in this series.

Fig. 11. Angle afil between the projected velocity of a halo (vrad, C) and the local direction of the closest filament as a function of the halo distance from the axis of the filaments dfil. Points represent the actual measurements, solid lines represent the running means of afil as a function of dfil, error bars encompass the region between the 16th and the 84th percentile of the distributions. Different panels refer to different combinations of DisPerSE parameters and galaxy selections.

We did perform the same analysis (in particular the one shown in Figs. 10 and 11) also on the cosmic web around our simulated Coma extracted when the FoG effect is included (see Sect. 3). In this case too, we do not see a difference in the velocity distribution as a function of halo distance from Coma for halos in filaments compared to the general halo population, which allows us to reach the conclusion that halos in filaments are not faster in terms of accretion velocity than the general halo population. When we analyze the a_fil − d_fil relation, we do not recover trends of a_fil to become closer to 180 deg with decreasing distance from the axes of the filaments, except for the case of magnitude selected galaxies and lower persistence threshold. However the value of a_fil remains stable at a large value close to 120 ÷ 140 deg also at large values of d_fil. This indicates that the flux of matter is collimated close to the filaments also when the FoG effect is included in the extraction of the skeleton.

7. Discussion

The first challenge we face in this paper, is determining how similar the filaments close to the simulated Coma are to the filaments around the real Coma. Our approach to this problem is different from what usually done in the literature due to the peculiar nature of our constrained simulations. For example, one approach commonly adopted is to apply different algorithms to detect the filaments to the same observed or simulated data set. Libeskind et al. (2018), Rost et al. (2020), and Bonnaire et al. (2020) performed these kinds of comparisons by means of distributions of global quantities for the filaments (e.g., length, radius, luminosity, mass, volume, redshift) or, in the case of simulations, by also comparing which components are associated to the same structure by different algorithms. We adopted a similar approach in Malavasi et al. (2020a), where we compared filament catalogs obtained from the same galaxy survey with the same algorithm, only run with varying parameters. In the case of the present work, this approach is less constraining. In fact, several of the filament properties (e.g., their length) could be influenced by the volume density of the tracers (i.e., galaxies), which is different in simulations and observations.

A complementary approach is to run the same algorithm on simulated samples of galaxies extracted from simulations light-cones tuned to reproduce a given galaxy survey, where the mock samples are slightly modified each time to test a different source of uncertainty in the filament reconstruction. In this case, a close correspondence between the various filament catalogs can be expected. More sophisticated approaches can then be undertaken to compare the catalogs, such as it was done in Malavasi et al. (2017), Laigle et al. (2018), and Kraljic et al. (2018) where the distances of all filaments from one skeleton realization to the filaments of another one can be derived and compared. Also this approach is not viable in our case as, although our simulations are constrained to reproduce existing structures, the residual cosmic variance does not ensure a perfect correspondence.

We have therefore to resort mainly to visual comparison of the samples of filaments and to set for a match between observations and simulations in terms of number and spatial distribution of the filaments around the simulated cluster. We do conclude that our filament reconstruction is overall satisfactory. In the majority of cases we identify prominent northeast and west filaments, connecting to the Coma cluster. This is important, as those are the filaments that observational evidence indicates as being the most linked to features of the Coma cluster connected to matter accretion.

When investigating the dynamics of matter around the simulated Coma cluster and all the filaments in its vicinity, we found that there is limited evidence for the halo velocity to be higher in the proximity of the filaments, but stronger trends for the angle between the direction of the matter flow and the direction of the filament to be smaller. This points us toward the conclusion that the matter flow close to filaments is more collimated. This is the first work to report such a result in the literature.

This conclusion is in agreement with what found in several works in the literature. In particular, Dekel et al. (2009), for example, highlighted the presence of high-flux, low-entropy streams accreting cold gas on a forming halo (M_vir > 10¹² M_⊙) at z = 2.5. Although this is a different mass and redshift range than what we explored in our paper, Fig. 1 of Dekel et al. (2009) shows that while the velocity field magnitude does not seem to be different inside and outside the filaments, the flux is higher and the entropy lower inside the streams. This is in agreement with Danovich et al. (2012) who analyzed 350 halos of 10¹² M_⊙ at z = 2.5, finding that the streams of accreting gas are narrow and covering a small portion of the virial shell of the accreting halos. This is also in agreement with what found by Bennett & Sijacki (2020) for a 10¹² M_⊙ halo at z = 6, where they identify filaments of cold gas piercing through the shock at the virial radius of the halo. Although the radial velocity magnitude of the gas in the filaments is not much lower than the velocity of other accreting gas from the large scale, its turbulent velocity component is low, meaning that these stream of gas are collimated. This is supported also by the velocity dispersion of the gas in the filaments being lower than in the surrounding large-scale structure (see their Fig. 8). Cold gas filaments that accrete matter on the forming halo inside the virial radius are also found by Valentini et al. (2021), for a 10¹² M_⊙ halo at z = 6. Coherent streams that accrete matter onto forming halos at high redshift are also highlighted in the simulations by Aragón-Calvo et al. (2019, see their Fig. 1), who also introduces the idea that when halos are detached from such flows galaxies can quench their star formation (see also Moutard et al. 2020 where this effect, called “cosmic web detachment” is introduced, when induced by halo-halo mergers, as an explanation for the larger fraction of heavily-obscured, X-ray emitting AGNs in quenching galaxies and the change in spin direction of galaxies flowing along the cosmic web). Such accretion flows have also been found in observations. For example Martin et al. (2016) reports a filament with a collimated flow having a very low velocity dispersion (∼50 km s⁻¹) feeding a proto-disk in a 4 × 10¹² M_⊙ halo at z = 2.843.

We stress the fact that the accretion flows of gas and dark matter are expected to be different, a fact that our N-body simulations are unable to probe. In fact, our analysis of the dark matter dynamics around Coma is based on the velocity of halos. Moreover, accretion of cold gas onto forming halos through filaments (a process explored in the literature cited above which uses hydrodynamical simulations) is a phenomenon limited to the high redshift Universe. Indeed, cold mode accretion is not the main mode of accretion onto halos today (clusters at z = 0 are expected to accrete mostly hot gas, inefficiently, through isotropic accretion more than from filaments of the cosmic web, see e.g., Fig. 11 of Overzier 2016; Dekel & Birnboim 2006). What we detect around Coma may be only minor accretion from filaments still happening at z = 0. However, we also conclude that our findings provide a compelling science case for the investigation of the evolution of the cosmic web around this simulated Coma cluster and how the accretion of matter from the filaments changed as a function of redshift and of the different configurations of the filaments themselves (e.g., because of the possible merging of the filaments, Cadiou et al. 2020). This will be the subject of upcoming papers in this series.

8. Conclusions

In this work, we explore the cosmic web around a reproduction of the Coma cluster of galaxies, identified in a box of 500 Mpc h⁻¹ on a side simulated with a dark matter only run of the code RAMSES. The initial conditions of this simulation have been constrained using the velocity of galaxies so that they reproduce existing structures in the local Universe.

We identified the halo corresponding to the Coma cluster in our simulation and we detected filaments connected to it with the DisPerSE algorithm. We then analyzed the filaments connected to the simulated cluster, comparing their number and spatial distribution with the observations of Malavasi et al. (2020b), and the velocity distribution of the halos around these filaments to study the accretion of matter onto the cluster. We summarize our findings in the following way:

1. We report a good agreement between the cosmic web detected around the simulated Coma cluster and the large-scale structure configuration in observations. We identified several combinations of the DisPerSE parameter choices and galaxy selection criteria where prominent northeast and west filaments were detected at the correct observed position.

2. The connectivity around the simulated Coma is in line with the value for the real Coma. However, given the difference between the mass of the real and simulated Coma (higher for the latter), the placement of the simulated Coma on the mass-connectivity relation is more uncertain. The value of κ = 3 is only in partial agreement with other values from the literature, both from observational and simulated samples.

3. The magnitude of the radial infall velocity components of the halos in a 10 Mpc (6.774 Mpc h⁻¹) radius around the cluster shows no significant difference between the general halo population and a selection of only those halos that are close to the axis of the surrounding filaments. However, the angle between the radial infall velocity component vector and the local direction of the closest filament to a given halo become increasingly aligned with decreasing distance of the halos from the filaments’ axis. We conclude that the filaments convey a collimated flow of matter onto the simulated Coma cluster. This is the first time that such a result is reported for dark matter at low redshift.

This paper is the first of a series in which we aim to establish a correspondence between our constrained simulation and the cosmic web around the real Coma cluster. The additional power of this simulation is that it is capable of offering us an evolutionary perspective on the Coma system. In future papers we will explore the evolution of the accretion of matter onto the Coma cluster from the filaments with redshift and how it impacted the cluster formation.

Acknowledgments

The authors would like to praise and thank the anonymous referee, whose on-point and in-depth comments helped greatly improve the quality of the paper. Part of this research has been supported by the funding for the ByoPiC project from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme grant agreement ERC-2015-AdG 695561. K.D. 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. J.S., K.D., and N.A. acknowledge support by the grant agreements ANR-21-CE31-0019/490702358 from the French Agence Nationale de la Recherche/DFG for the LOCALIZATION project. The authors acknowledge the Gauss Centre for Supercomputing e.V. (http://www.gauss-centre.eu) and GENCI (https://www.genci.fr/) for funding this project by providing computing time on the GCS Supercomputer SuperMUC-NG at Leibniz Supercomputing Centre (http://www.lrz.de) and Joliot-Curie at TGCC (http://www-hpc.cea.fr). Funding for the SDSS and SDSS-II has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, the US Department of Energy, the National Aeronautics and Space Administration, the Japanese Monbukagakusho, the Max Planck Society, and the Higher Education Funding Council for England. The SDSS Web Site is http://www.sdss.org/. The SDSS is managed by the Astrophysical Research Consortium for the Participating Institutions. The Participating Institutions are the American Museum of Natural History, Astrophysical Institute Potsdam, University of Basel, University of Cambridge, Case Western Reserve University, University of Chicago, Drexel University, Fermilab, the Institute for Advanced Study, the Japan Participation Group, Johns Hopkins University, the Joint Institute for Nuclear Astrophysics, the Kavli Institute for Particle Astrophysics and Cosmology, the Korean Scientist Group, the Chinese Academy of Sciences (LAMOST), Los Alamos National Laboratory, the Max-Planck-Institute for Astronomy (MPIA), the Max-Planck-Institute for Astrophysics (MPA), New Mexico State University, Ohio State University, University of Pittsburgh, University of Portsmouth, Princeton University, the United States Naval Observatory, and the University of Washington. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. This work has made use of CosmoHub. CosmoHub has been developed by the Port d’Informació Científica (PIC), maintained through a collaboration of the Institut de Física d’Altes Energies (IFAE) and the Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT) and the Institute of Space Sciences (CSIC & IEEC), and was partially funded by the “Plan Estatal de Investigación Científica y Técnica y de Innovación” program of the Spanish government. The CosmoSim database used in this paper is a service by the Leibniz-Institute for Astrophysics Potsdam (AIP). The MultiDark database was developed in cooperation with the Spanish MultiDark Consolider Project CSD2009-00064. The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (http://www.gauss-centre.eu) and the Partnership for Advanced Supercomputing in Europe (PRACE, http://www.prace-ri.eu) for funding the MultiDark simulation project by providing computing time on the GCS Supercomputer SuperMUC at Leibniz Supercomputing Centre (LRZ, http://www.lrz.de). The Bolshoi simulations have been performed within the Bolshoi project of the University of California High-Performance AstroComputing Center (UC-HiPACC) and were run at the NASA Ames Research Center.

References

Appendix A: Masses for SDSS galaxies

Galaxy masses were obtained from the Max Planck Institute for Astrophysics-Johns Hopkins University Value Added Catalog (MPA-JHU VAC⁷). Spectral quantities in this catalog were obtained with the GALSPEC software, following Brinchmann et al. (2004), Kauffmann et al. (2003), and Tremonti et al. (2004). Masses were derived with a fit to the ugriz photometry corrected with spectral information, following Kauffmann et al. (2003) and Salim et al. (2007). We make use of the median value of the mass distribution output from the Bayesian fit to the photometry, available in the catalog. We start with an initial sample of 927 552 galaxies in the MPA-JHU VAC, from which we remove duplicate observations and galaxies not in the MGS or in special releases (PRIMTARGET & 64, RELEASE ≠ “extra”, “special”, and “extraspecial”). We select only objects classified as galaxies spectroscopically (SPECTROTYPE = “galaxy”) and with good redshift measurements (ZWARNING = 0 and z > 0). The result of these selections is a sample of 654 324 galaxies which we match in position to the Legacy MGS with a tolerance of 0.5″. We check that the redshift of the matched objects is consistent between the MPA-JHU VAC and the Legacy MGS and we identify 520 682 galaxies with measured stellar mass. We stress that in this work we make use of galaxy masses only to calibrate our assumptions for the selection of halos corresponding to galaxies in the numerical simulation (see Section 2.3) and we do not derive physical conclusions on the Coma system or its filaments based on this data set.

Figure A.1 shows the mass distribution for the halos and galaxies in our simulation and for the Legacy MGS galaxies (with mass derived by matching the Legacy MGS sample with the MPA-JHU VAC sample). The mass distribution of our mass-selected galaxy sample is largely consistent with the expected mass distribution for the halos hosting the galaxies of the Legacy MGS (this distribution was derived by shifting the mass distribution of Legacy MGS galaxies by one order of magnitude to higher masses, for reference). As for the mass distribution of our magnitude-selected simulated galaxies, it is somewhat overlapping with the expected mass distribution of the halos hosting the Legacy MGS galaxies, although marginally. However we stress that the matching of the Legacy MGS with the MPA-JHU VAC is itself somewhat uncertain, while the expected distribution of the halo mass of the Legacy MGS galaxies is only for reference. We therefore use this figure only to check our parameter selection rather than to drive it.

Fig. A.1. Mass distributions for the halos in our simulation and for SDSS galaxies. The solid light blue line shows the mass distribution for Legacy MGS galaxies, while the dashed light blue line is the same distribution shifted by one order of magnitude to provide a reference for the expected position of the mass distribution of the halos hosting the Legacy MGS galaxies. The red line is the mass distribution for simulated galaxies selected through a cut in mass and the green line is the mass distribution of simulated galaxies selected through a cut in magnitude. We note that the mass reported on the x-axis is the stellar mass for SDSS galaxies and the virial mass for the halos and galaxies in our simulation.

The expected mass distribution for the halos hosting the galaxies of the Legacy MGS was derived in this case assuming a Stellar-to-Halo Mass Relation (SHMR) expressed in terms of a fixed ratio of halo to stellar mass of 1 dex. We do note that this approach, while sufficient for us to check our parameter selection for the extraction of a simulated galaxy population from our parent halo population, may not reflect the true distribution of mass for halos hosting the Legacy MGS galaxies that could be obtained using variable SHMRs found in the literature. We have therefore compared the halo mass distribution of our mass and magnitude selected galaxies with the expected mass distributions for halos hosting Legacy MGS galaxies obtained using SHMRs from the literature.

In particular, we used Equation 6 (with best-fit values from the first row of Table 1) of Girelli et al. (2020) to derive the stellar mass distribution for halos in our mass-selected and magnitude-selected simulated galaxy samples. We find these to be in very good agreement with the stellar mass distribution for Legacy MGS galaxies obtained through the matching with the MPA-JHU VAC.

We also used the Data Release 1 (DR1) of the UNIVERSEMACHINE code applied to the Bolshoi-Planck simulation (Klypin et al. 2016; Rodríguez-Puebla et al. 2016) by Behroozi et al. (2019)⁸. We make use of the light-cones created to reproduce the COSMOS field. While we are aware that this means a limited area and therefore a restricted number of halos, we can use this light-cone to select objects at z ≤ 0.3 (to match the redshift limit of the Legacy MGS) and in two mass ranges: 10¹² ≤ M_vir/M_⊙ ≤ 10¹³ (to simulate our mass-selected simulated galaxy sample) and M_vir/M_⊙ ≥ 10¹² (to simulate the absence of an upper mass limit as is the case for our magnitude-selected simulated galaxy sample, although the actual mass distribution may differ). We then compare the stellar mass distribution for these halo samples, as derived by Behroozi et al. (2019). We find good agreement with the stellar mass distribution for Legacy MGS galaxies obtained through the matching with the MPA-JHU VAC for both cases.

Additionally, we use the simulated galaxies generated by Cora et al. (2018) by applying the semi-analytic model of galaxy formation SAG (Cora 2006; Lagos et al. 2008; Orsi et al. 2014; Gargiulo et al. 2015; Muñoz Arancibia et al. 2015; Cora et al. 2018) to the Multidark MDPL2 simulation (Klypin et al. 2016)⁹. We select halos with M_200c ≥ 10¹² M_⊙ and z ≤ 0.3. We then derived their total stellar mass by summing the stellar mass of the disk and spheroid components and we compared its distribution to the stellar mass distribution for Legacy MGS galaxies, finding good agreement.

Finally, we use the MICECATv1.0¹⁰ catalog of Carretero et al. (2015), created from the Marenostrum Institut de Ciències de l’Espai (MICE) Grand Challenge (MICE-GC) simulation (Crocce et al. 2015; Fosalba et al. 2015a,b). From this halo catalog with photometric quantities computed by Carretero et al. (2015), we selected halos with z ≤ 0.3 and apparent r-band magnitude r ≤ 17.77 to mimic the SDSS selection. We then compared the halo mass distribution to the halo mass distribution for mass-selected and magnitude-selected galaxies. We find close agreement.

Our conclusion is that although choosing a fixed SHMR may not be the most accurate way to convert between halo and stellar mass, it provides a good enough approximation for our goal of checking whether our parameter choice for the selection of a simulated galaxy sample is correct. We postpone to future works the use of a semi-analytical model to populate halos with galaxies in a realistic way.

Appendix B: Extra projection

In this section we provide an additional projection to the one showed in Figure 2. This projection is located on the plane of the sky and shows the northeast and west filaments system in the vicinity of Coma. This figure shows filaments leaving the simulated Coma cluster more or less at the correct directions.

Fig. B.1. Filaments in the constrained simulation compared to filaments in observations. The left column shows the case of a mass selected galaxy sample in the constrained simulation, the right column shows the case of a magnitude selected galaxy sample. In each panel, simulated galaxies are shown as gray points. Blue lines are filaments from Malavasi et al. (2020b) in supergalactic coordinates, the red square marks the position of the Coma cluster. Dark red lines are filaments obtained in the simulation with a variety of DisPerSE parameter combinations, namely: SD0−2σ (top row), SD0−3σ (second row), SD1−2σ (third row), SD1−3σ (bottom row). Circle and triangles mark the positions of critical points (in observations and simulations, respectively) and are color-coded according to their type (red: maxima, green: type 2 saddles, brown: bifurcations). This figure shows the SGX-SGZ projection (in this projection, the increasing redshift direction is along the y axis). The observed northeast and west filaments are marked for reference.

All Tables

All Figures

Fig. 1. Apparent K-band magnitude distributions for the halos in our simulation and for SDSS galaxies. The orange line refers to the Legacy MGS galaxies with apparent K-band magnitude from the NYU-VAGC sample. The light blue line refers to the mK distribution for all halos in our simulation (with mK computed with Eq. (1)). The green line shows the distribution for our magnitude-selected simulated galaxy sample. The vertical black dashed line shows our magnitude threshold of mK = 14.

In the text

Fig. 2. Filaments in the constrained simulation compared to filaments in observations. The left column shows the case of a mass selected galaxy sample in the constrained simulation, the right column shows the case of a magnitude selected galaxy sample (see text). In each panel, simulated galaxies are shown as gray points in a slice of thickness 50 Mpc h−1 centered on Coma. Blue lines are filaments from Malavasi et al. (2020b) in supergalactic coordinates, the red square marks the position of the Coma cluster. Dark red lines are filaments obtained in the simulation with a variety of DisPerSE parameter combinations, namely: SD0−2σ (top row), SD0−3σ (second row), SD1−2σ (third row), SD1−3σ (bottom row). The observed northeast and west filaments are marked for reference. Circle and triangles mark the positions of critical points (in observations and simulations, respectively) and are color-coded according to their type (red: maxima, green: type 2 saddles, brown: bifurcations). This figure shows the SGZ-SGY projection.

In the text

Fig. 3. Filament length distribution. The black bold distribution refers to observed filaments from Malavasi et al. (2020b) with SD1−3σ. The colored lines refer to simulated filaments with parameter combinations SD0−2σ and SD0−3σ (for magnitude selected galaxies, in purple and red, respectively), and SD1−2σ and SD1−3σ (for mass selected galaxies, in orange and green, respectively). Dashed gray lines represent the rest of the parameter combinations and provide an idea of the range of lengths explored. Only filaments within a slice of thickness ±75 Mpc (50.8 Mpc h−1) from the location of the (simulated or observed) Coma cluster were used to derive the distributions.

In the text

Fig. 4. Node connectivity distribution. The black bold distribution refers to observed nodes from Malavasi et al. (2020b) with SD1−3σ. The colored lines refer to simulated nodes with parameter combinations SD0−2σ and SD0−3σ (for magnitude selected galaxies, in purple and red, respectively), and SD1−2σ and SD1−3σ (for mass selected galaxies, in orange and green, respectively). Dashed gray lines represent the rest of the parameter combinations and provide an idea of the range of connectivities explored. Only nodes within a slice of thickness ±75 Mpc (50.8 Mpc h−1) from the location of the (simulated or observed) Coma cluster were used to derive the distributions.

In the text

Fig. 5. Node connectivity distribution. The black bold distribution refers to observed nodes from Malavasi et al. (2020b) with SD1−3σ. The colored lines refer to simulated nodes with parameter combinations SD0−2σ and SD0−3σ (for magnitude selected galaxies, in purple and red, respectively), and SD1−2σ and SD1−3σ (for mass selected galaxies, in orange and green, respectively). Dashed gray lines represent the rest of the parameter combinations and provide an idea of the range of connectivities explored. Only nodes within a slice of thickness ±75 Mpc (50.8 Mpc h−1) from the location of the (simulated or observed) Coma cluster and with overdensity log(1 + δ) = log(1 + δComa)±0.5 (see text) were used to derive the distributions.

In the text

Fig. 6. Node density distribution. The black bold distribution refers to observed nodes from Malavasi et al. (2020b) with SD1−3σ. The colored lines refer to simulated nodes with parameter combinations SD0−2σ and SD0−3σ (for magnitude selected galaxies, in purple and red, respectively), and SD1−2σ and SD1−3σ (for mass selected galaxies, in orange and green, respectively). Dashed gray lines represent the rest of the parameter combinations and provide an idea of the range of connectivities explored. Only nodes within a slice of thickness ±75 Mpc (50.8 Mpc h−1) from the location of the (simulated or observed) Coma cluster and with κ = 2 ÷ 3 are used to derive the distributions.

In the text

Fig. 7. Filaments in the constrained simulation compared to filaments in observations when the FoG effect is introduced in our simulation. The left column shows the case of a mass selected galaxy sample in the constrained simulation, the right column shows the case of a magnitude selected galaxy sample (see text). In each panel, simulated galaxies are shown as gray points in a slice of thickness 50 Mpc h−1 centered on Coma. Blue lines are filaments from Malavasi et al. (2020b) in supergalactic coordinates, the red square marks the position of the Coma cluster. Dark red lines are filaments obtained in the simulation with a variety of DisPerSE parameter combinations, namely: SD0−2σ (top row), SD0−3σ (second row), SD1−2σ (third row), SD1−3σ (bottom row). The observed northeast and west filaments are marked for reference. Circle and triangles mark the positions of critical points (in observations and simulations, respectively) and are color-coded according to their type (red: maxima, green: type 2 saddles, brown: bifurcations). This figure shows the SGZ-SGY projection.

In the text

Fig. 8. Connectivity of the real and simulated Coma clusters as a function of mass. The cyan square refers to the real Coma cluster as analyzed in Malavasi et al. (2020b), while the red square corresponds to the simulated Coma cluster as detected in this work. Black crosses are the observed connectivity values from the AMASCFI Clusters (Sarron et al. 2018, 2019), green triangles are the observed connectivity values for the groups in COSMOS (Darragh Ford et al. 2019). Gray lines and points are connectivity measurements derived in the IllustrisTNG simulation (Gouin et al. 2021) for clusters and groups which are unrelaxed and old (UO, solid line), relaxed and old (RO, dotted line), and unrelaxed and young (UY, dashed line). The purple line and shaded region are the relation from the numerical N-body simulations of Aragón-Calvo et al. (2010), and its corresponding 1σ uncertainty, while the solid orange line is the theoretical relation of Codis et al. (2018) renormalized to pass through our measurement for the real Coma cluster, so as to provide a better comparison of trends. We note that the mass on the x-axis of the plot is M200 for the real and simulated Coma, the measurements by Darragh Ford et al. (2019), and those of Sarron et al. (2019), while it is a value close to Mvir for the Aragón-Calvo et al. (2010) relation.

In the text

Fig. 9. Maps of halo velocities around the simulated Coma cluster. The gray square at the center marks the position of the Coma cluster. Orange lines are first, second, and third generation filaments (the case SD0−2σ for magnitude selected galaxies is reported here as an example, see brown lines in Fig. 2 and text). The top left panel shows the x − y projection, the top right panel shows the x − z projection, and the bottom left panel shows the z − y projection (shown in the first panel on the right column of Fig. 2). In each slice, the color map shows a 2D histogram in which in each bin the average vrad, C is measured (red stands for halos moving away from Coma, blue for halos moving toward Coma, white pixels stand for empty bins). Black arrows are constructed from the average in each bin of the vrad, C components relative to the projection shown (i.e., vrad, C, x and vrad, C, y in the top left panel, vrad, C, x and vrad, C, z in the top right panel, vrad, C, z and vrad, C, y in the bottom left panel). The thickness of the considered slice is of 20 Mpc h−1, all halos in the slice are considered, regardless of mass or whether they are selected as galaxies. The labels mark the position of the simulated northeast and west filaments.

In the text

Fig. 10. Halo velocity as a function of distance from Coma. Each panel shows the halo velocity vrad, C as a function of the distance between each halo and the center of the simulated Coma cluster dh − Coma = |rh|. Only halos with rvir ≤ dh − Coma ≤ 10 Mpc (6.774 Mpc h−1) have been considered. Blue squares refer to all halos within this distance range, orange triangles are halos with dfil smaller than the 25th percentile of the dfil distribution. Solid lines refer to the running means of the distributions (cyan: all halos, red: halos close to filaments), error bars encompass the region between the 16th and the 84th percentile of the distributions. Different panels refer to different combinations of DisPerSE parameters and galaxy selections. In each panel, the black line shows the velocity that a body would reach as a function of distance from Coma by starting at 10 Mpc (6.774 Mpc h−1), with a velocity equal to the value of the last bin of the running mean as derived for the total halo population, and being subject only to the gravitational attraction of the Coma cluster.

In the text

Fig. 11. Angle afil between the projected velocity of a halo (vrad, C) and the local direction of the closest filament as a function of the halo distance from the axis of the filaments dfil. Points represent the actual measurements, solid lines represent the running means of afil as a function of dfil, error bars encompass the region between the 16th and the 84th percentile of the distributions. Different panels refer to different combinations of DisPerSE parameters and galaxy selections.

In the text

Fig. A.1. Mass distributions for the halos in our simulation and for SDSS galaxies. The solid light blue line shows the mass distribution for Legacy MGS galaxies, while the dashed light blue line is the same distribution shifted by one order of magnitude to provide a reference for the expected position of the mass distribution of the halos hosting the Legacy MGS galaxies. The red line is the mass distribution for simulated galaxies selected through a cut in mass and the green line is the mass distribution of simulated galaxies selected through a cut in magnitude. We note that the mass reported on the x-axis is the stellar mass for SDSS galaxies and the virial mass for the halos and galaxies in our simulation.

In the text

Fig. B.1. Filaments in the constrained simulation compared to filaments in observations. The left column shows the case of a mass selected galaxy sample in the constrained simulation, the right column shows the case of a magnitude selected galaxy sample. In each panel, simulated galaxies are shown as gray points. Blue lines are filaments from Malavasi et al. (2020b) in supergalactic coordinates, the red square marks the position of the Coma cluster. Dark red lines are filaments obtained in the simulation with a variety of DisPerSE parameter combinations, namely: SD0−2σ (top row), SD0−3σ (second row), SD1−2σ (third row), SD1−3σ (bottom row). Circle and triangles mark the positions of critical points (in observations and simulations, respectively) and are color-coded according to their type (red: maxima, green: type 2 saddles, brown: bifurcations). This figure shows the SGX-SGZ projection (in this projection, the increasing redshift direction is along the y axis). The observed northeast and west filaments are marked for reference.

In the text