A survey of the thermal and non-thermal properties of cosmic filaments

ABSTRACT

In this paper, we exploit a large suite of enzo cosmological magneto-hydrodynamical simulations adopting uniform mesh resolution, to investigate the properties of cosmic filaments under different baryonic physics and magnetogenesis scenarios. We exploit an isovolume-based algorithm to identify filaments and determine their attributes from the continuous distribution of gas mass density in the simulated volumes. The global (e.g. mass, size, mean temperature and magnetic field strength, and enclosed baryon fraction) and internal (e.g. density, temperature, velocity, and magnetic field profiles) properties of filaments in our volume are calculated across almost four orders of magnitude in mass. The inclusion of variations in non-gravitational physical processes (radiative cooling, star formation, feedback from star forming regions, and active galactic nuclei) as well as in the seeding scenarios for magnetic fields (early magnetization by primordial process versus later seeding by galaxies) allows us to study both the large-scale thermodynamics and the magnetic properties of the Warm–Hot Intergalactic Medium (WHIM) with an unprecedented detail. We show how the impact of non-gravitational physics on the global thermodynamical properties of filaments is modest, with the exception of the densest gas environment surrounding galaxies in filaments. Conversely, the magnetic properties of the WHIM in filament are found to dramatically vary as different seeding scenarios are considered. We study the correlation between the properties of galaxy-sized haloes and their host filaments, as well as between the haloes and the local WHIM in which they lie. Significant general statistical trends are reported.

1 INTRODUCTION

Galaxy surveys and numerical simulations show that the large-scale structure of the Universe is organized in the form of filaments, haloes, and voids. Numerical simulations predict also that a large fraction of the baryonic matter (around 50 per cent) resides in the form of plasma in such cosmic web, at densities ∼10–100 times the average cosmic value and 10⁵–10⁷ K temperatures, forming the ‘Warm–Hot Intergalactic Medium’ (WHIM, see, e.g. Bond, Kofman & Pogosyan 1996; Davé et al. 2001). Direct observations of the cosmic web are challenging, due to its extremely low mass density. The cooler phases of the intergalactic gas (⁠|${\sim } 10^4\, \rm K$|⁠) contain a small fraction of neutral hydrogen, producing the characteristic Ly-α absorption. Indications of possible detection of filamentary structures emerged from the analysis of soft X-ray (e.g. Finoguenov, Briel & Henry 2003; Werner et al. 2008; Nicastro et al. 2010; Nicastro 2016), or of the Sunyaev–Zeldovich effect (e.g. Planck Collaboration et al. 2013; Tanimura et al. 2017; de Graaff et al. 2017).

The first X-ray imaging of the terminal part of four filaments connected to the virial radius of cluster A2744 has been reported by Eckert et al. (2015) using XMM–Newton. By using the Sardinia Radio Telescope (SRT), Vacca et al. (2018) reported the possible detection of large-scale diffuse emission around a giant filament connecting clusters. However, the emission could be associated to the outer regions of clusters rather than to the WHIM in filaments. Finally, the possible detection of a Faraday Rotation excess produced by filaments overlapping the polarized emission of a giant radio galaxy at z = 0.34 has been proposed by O’Sullivan et al. (2018), based on LOFAR observations. Statistical techniques via cross-correlation analysis have also attempted to detect the signature from cosmic filaments emitting in radio, but reported only upper limits (e.g. Brown et al. 2017; Vernstrom et al. 2017). Complementary to this, statistical studies of the Faraday Rotation measurement from background sources were used to limit the magnetization of the cosmic web intervening on to high-redshift polarized sources (Blasi, Burles & Olinto 1999; Xu & Han 2015; Oppermann et al. 2015; Pshirkov, Tinyakov & Urban 2016).

In the radio domain, the interest in detecting filaments is rapidly growing, triggered by the fast development of precursors and pathfinders of the Square Kilometer Array (SKA), which will have the potential of scratching the surface of the magnetic cosmic web (e.g. Brown 2011; Vazza et al. 2015a). Current radio observatories like LOFAR, ASKAP, MWA, and MeerKAT may already lead to early detection in a non-negligible number of targets (Vazza et al. 2015b), and sophisticated techniques for the automated detection of very faint and diffuse emission on |${\sim } \rm ~degree$| scales in such surveys are already available (e.g. Gheller, Vazza & Bonafede 2018).

Filaments represent an ideal environment to investigate the past epochs that led to present cosmic structures. In fact, their dynamics is less violent and complex than that of galaxy clusters or of galaxies, with adiabatic physics (besides gravity) dominating the gas dynamics, but with other physical processes influencing the behaviour and the properties of the gas component. Among these processes, energy injection from supernova explosion or AGN jet heating have already proved to be crucial to explain observed features of dense environments (e.g. preventing runaway cooling in dense galaxy clusters cores, or regulating the star formation efficiency (see Borgani et al. 2008; Martizzi et al. 2018) and references therein). In filaments, thanks to their more gentle evolution, the interplay of the different physical processes can be better understood and their effects on observable signatures effectively disentangled. Furthermore, filamentary structures are expected to preserve many traces of the original environment in which the process of gravitational clustering started. This is particularly relevant when magnetic fields come into play. Although all the cosmological structures are expected to be substantially magnetized, the origin and the evolution of such magnetic fields is currently unsettled. Filaments should retain memory of the initial magnetic seed fields since they should not host strong dynamo amplification (e.g. Ryu et al. 2008; Donnert et al. 2009; Vazza et al. 2014a). Furthermore, although astrophysical processes can contribute to the evolution of such primordial seeds, in close connection with the dynamical history of the cosmic structures, their signatures should remain confined to the highest density regions, influencing only partially the overall filaments’ volume.

When first accreted on to filaments, the cosmic gas is shock-heated by strong shocks, |$\mathcal {M} \sim 10-100$| (e.g. Ryu et al. 2003; Pfrommer et al. 2006). Downstream of these shocks, supersonic turbulence is injected (e.g. Kang et al. 2007; Ryu et al. 2008; Wittor, Vazza & Brüggen 2017b). When such turbulence gets dissipated, a tiny fraction of the resulting energy is expected to feed the amplification of magnetic fields (Ryu et al. 2008; Federrath et al. 2011; Jones et al. 2011; Schleicher et al. 2013; Porter, Jones & Ryu 2015). Most of cosmological MHD simulations targeting filaments have found little presence of volume-filling dynamo amplification of seed fields in filaments (Brüggen et al. 2005; Vazza et al. 2014a; Marinacci et al. 2015; Vazza et al. 2017), which is well explained by the little time available for amplification in such environment, magnetic eddies being rapidly advected on to the surrounding cluster environment before being significantly amplified (e.g. Vazza et al. 2014a).

In this scenario, most of filaments’ volume must be filled with magnetic fields whose evolution mostly come from compression/rarefaction of primordial magnetic field lines, following |$B \approx B_{0} \cdot ({n_e} /{\langle n \rangle)}^{\alpha _B}$|⁠, where B₀ is the seed field, 〈n〉 is the cosmic mean (gas) density, and α_B ≈ 2/3 for isotropic gas compression. If this picture is qualitatively correct, then from the observation of the present-day magnetic field level of filaments, it would be possible to infer the amplitude of primordial magnetic fields to a good approximation.

In a different model, the predominant mechanism for the origin of cosmic magnetic fields is the ‘magnetic pollution’ by active galactic nuclei and galactic activities (e.g. Donnert et al. 2009). Relatively low magnetization is expected in the vast majority of filaments, due to a combination of the adiabatic decrease of the field, and of the scarcity of magnetic field sources there (e.g. Xu et al. 2009; Marinacci et al. 2015; Vazza et al. 2017). Furthermore, magnetic fields would be originated following the distribution of galaxies and AGN, filling in an intermittent and discontinuous way the filament’s volume. This would make their characterization even more challenging.

Assessing the magnetization level of filaments is of great interest also for other studies, from the origin and composition of ultra high-energy cosmic rays (UHECRs) propagating across the Universe (e.g. Sigl, Miniati & Ensslin 2003; Dolag et al. 2003; Hackstein et al. 2016, 2018) to the study of rotation measure for cosmological sources, like quasars and Fast Radio Bursts (e.g. Akahori, Ryu & Gaensler 2016; Vazza et al. 2018a). For example, it has been recently argued that the excess of |${\sim } 6 \cdot 10^{19} \, \rm eV$| events reported by the Telescope Array may be due to the trapping of UHECRs by filaments local to the Virgo cluster (Kim et al. 2019).

In this work, we address several of the above topics by studying in detail the thermal and non-thermal (magnetic) properties of the gas within cosmic filaments by means of cosmological magneto-hydrodynamics (MHD) simulations. Different physical models have been considered, selecting a number of simulations from the recent ‘Chronos++ suite’ Vazza et al. 2017), a collection of runs aimed at exploring a number of scenarios for the emergence of large-scale magnetic fields, including different physical and astrophysical processes, encompassing primordial, dynamo, and astrophysical magnetogenesis scenarios.

Filaments are identified within the simulated data sets using the methodology already adopted in Gheller et al. (2015) and Gheller et al. (2016). Such methodology is based on the calculation of baryonic matter density isovolumes to identify filamentary structures. In this, it is different from other approaches used for N-body simulations, that can be broadly grouped in geometric and tessellation methods, based on the topological analysis of the density field by means of sophisticated mathematical approaches (e.g. Stoica et al. (2005); Sousbie et al. (2008); Aragon-Calvo, Shandarin & Szalay (2010); González & Padilla (2010); Chen et al. (2015)) and morphological methods, that classify the 3D distribution of matter based on the density Hessian and/or the tidal or velocity shear fields (e.g. Hahn et al. (2007); Aragón-Calvo et al. (2007); Cautun et al. (2014)). Exploiting the same methodology, it was also possible to extract both the galaxy haloes, in order to correlate their properties with those of the hosting filaments, and the local environment where haloes reside, in order to investigate their characteristics compared to those of the resident haloes. In particular, the cosmological missing baryons problem is addressed, by analysing the baryon fraction of the haloes, of the local environments, and of the overall filament.

The details of the simulations, the code used for the runs, the models selected from the Chronos++ suite are presented in Section 2, together with a summary of the adopted filament extraction methodology. Section 3 is dedicated to the discussion of the statistical and internal properties of the filaments extracted from the different runs, while Section 4 focuses on the statistical properties of the haloes and their correlations with the hosting filament and local environment. In Section 5, we discuss the various results, while Section 6 draws the conclusions.

2 THE NUMERICAL APPROACH

The simulations presented in this paper have been performed with the cosmological Eulerian code enzo (Bryan et al. 2014), with a fixed-mesh resolution. The code has been customized by our group mainly with the purpose of including different mechanisms for the seeding of magnetic fields in cosmology, as explained in detail in Vazza et al. (2017).

The magneto-hydrodynamics (MHD) solver adopted in our simulations implements the conservative Dedner formulation (Dedner et al. 2002), which uses hyperbolic divergence cleaning to enforce the |$\nabla \cdot \vec{B} = 0$| condition. The Dedner cleaning method has already been tested by several works in the literature, showing that despite the relatively large rate of dissipation introduced by its ‘cleaning waves’, it always converges to the correct solution as resolution is increased, at variance with other possible ‘divergence cleaning’ methods (Stasyszyn, Dolag & Beck 2013; Hopkins & Raives 2016; Tricco, Price & Federrath 2016). We refer the reader to more recent reviews for a broader discussion of the resolution and accuracy of different MHD schemes in properly resolving the dynamo in cosmological simulations (Donnert et al. 2018).

The MHD solver adopts the Piecewise Linear Method reconstruction technique; fluxes at cell interfaces are calculated using the Harten–Lax–Van Leer (HLL) approximate Riemann solver. Time integration is performed using the total variation diminishing (TVD) second-order Runge–Kutta (RK) scheme (Shu & Osher 1988). We used the GPU-accelerated MHD version of enzo by Wang, Abel & Kaehler (2010), which overall ensures a ∼×4 speedup compared to the more strandard CPU version of the code when applied to the large, 1024³ uniform grids used here. Constant spatial resolution has the advantage of providing the best-resolved description of magnetic fields even in low-density regions, which would typically not be refined by AMR schemes.

A known limitation in our numerical model is the diffusivity of the HLL scheme combined with the Dedner cleaning, which reduce the actual dynamical range. Compared to other MHD methods such as as the Constrained Transport, the Dedner scheme is more affected by small-scale dissipation of magnetic fields, due to the |$\nabla \cdot \vec{B}$| cleaning waves it generates to keep the numerical divergence under control. Nevertheless, several works in the literature have shown that the method is robust and accurate for most tests in MHD, also including magnetic turbulence in idealized setups (e.g. Wang & Abel 2009; Wang et al. 2010; Kritsuk et al. 2011; Bryan et al. 2014). Various astrophysical applications showed that this method converges to the right solution in control tests, unlike other common cleaning or |$\nabla \cdot \vec{B}$| preserving techniques (Stasyszyn et al. 2013; Hopkins & Raives 2016; Tricco et al. 2016). Moreover, resolution tests presented by our group in Vazza et al. (2014a) and Vazza et al. (2018b) have convincingly shown how at sufficiently high resolution, the scheme can properly model small-scale dynamo in a large-Reynolds-number regime. Given the fairly limited spatial resolution of our simulations here (83.3 kpc/cell), the negligible impact of the small-scale dynamo here may be exacerbated by the limited resolution. Based on a dedicated suite of higher resolution simulations, in Vazza et al. (2014a), we concluded that the development of the small-scale dynamo in filaments is hampered by little dynamical time available for amplification in this fast moving flow, hence the lack of dynamo amplification here looks physically motivated. However, the simplistic hydro-MHD view cannot exclude that in real plasmas ‘microscopic’ small-scale instabilities arising from kinetic plasma effects (e.g. Mogavero & Schekochihin 2014, for a review) can further promote a small-scale dynamo starting from sub-kpc scales. For this reason, a subset of our models (‘subgrid dynamo models’, see below) has been explicitly designed to bracket the maximal amount of magnetic field amplification that would be possible in our filaments, which is basically set by the available kinetic energy.

Our study uses the suite of cosmological MHD simulations at a fixed-mesh resolution, known as the ‘Chronos++ suite’.¹ It encompasses a large number of models which explore various plausible scenarios for the evolution of large-scale magnetic fields, with the goal of bracketing the distribution of magnetic fields outside of virialized haloes. In this work, we focus on one ‘baseline’ non-radiative simulation with a simple primordial seeding of magnetic fields, and we contrast it with 10 additional resimulations out of the full suite of 24 models presented in Vazza et al. (2017), which overall give a good match to the observed cosmic average star formation rate as well as to observed galaxy cluster scaling relations, as discussed there.

The Chronos++ suite includes three basic classes of models: primordial, dynamo, and astrophysical magnetogenesis scenarios (or a mix of the them). In detail:

• primordial models: we assumed the existence of weak and volume-filling magnetic fields at the beginning of the simulation. These are the only seed of magnetization in our simulated universe. We set either a spatially uniform seed field value for every component of the initial magnetic field B₀ (Baseline model in Table 1), or a primordial field with orientation orthogonal to the initial 3D velocity field from the ‘Zeldovich approximation’, used by the standard cosmological initial condition generators (Z2 model in Table 1) which ensures a |$\nabla \cdot \vec{B}=0$| field, being the Zeldovich approximation irrotational by construction. We ensured by construction that (〈B²〉)^1/2 = B₀, which allows us to closely compare the above two models. We consider here |$B_0=1 ~\rm nG$| (comoving) at z_in = 38, consistent with the most recent upper limits derived from the analysis of the Cosmic Microwave Background by Ade et al. (2015), of order |$B_0 \le 5 ~\rm nG$| for fields with a coherence scale of |${\sim } ~\rm Mpc$|⁠, although even lower limits have been suggested (e.g. Trivedi, Subramanian & Seshadri 2014).

• dynamo models: we estimated at run-time via sub-grid modelling the small-scale dynamo amplification of seed field of primordial origin (⁠|$B_0 = 10^{-9} \, \rm nG$| in runs DYN5² and DYN7 or |$B_0=0.01 \rm ~nG$| in run DYN8). The dissipation of solenoidal turbulence is estimated following Federrath et al. (2014), by extrapolating the information resolved at our fixed |$83.3 ~\rm kpc$| cell resolution on smaller scales. These models attempt to bracket the possible residual amount of dynamo amplification which can be attained, on energy grounds, in our cosmic structures, but may be lost because of a lack of resolution and/or for the onset of small-scale plasma instabilities that cannot be captured from first principles in our simple hydro-MHD model. In this work, we consider models which give a reasonable match to the magnetic field strength in galaxy clusters. More details on this procedure can be found in Vazza et al. (2017), as well as in Appendix  B.

• astrophysical models: magnetic fields are seeded by a) stellar activity and/or b) feedback by supermassive black holes (SMBH), simulated at run time using prescriptions available in enzo (e.g. Kravtsov 2003; Kim et al. 2011; Bryan et al. 2014), with tuneable star formation efficiency, timescales, and feedback parameters. Our assumed accretion model for SMBH follows from the spherical Bondi–Hoyle formula, without taking into account the Eddington accretion limit (run CSFBH3 and CSFBH5) or assuming a fixed |$0.01 ~\rm M_{\odot }\,yr^{-1}$| accretion rate (CSFBH2) which are all options available in the enzo public version, that we have tuned to produce realistic accretion and feedback results at the spatial and mass resolution covered here. In all cases, we consider a ‘boost’ factor to the mass growth rate of SMBH (α_Bondi = 100 in CSFBH5 or =1000 otherwise), again to account for the effect of coarse resolution in properly resolving the mass accretion rate on to our simulated SMBH particles. We have extended enzo coupling thermal feedback to the direct injection of additional magnetic energy via bipolar jets, with an efficiency with respect to the feedback energy computed at run-time, ϵ_{SF, b} and ϵ_{BH, b} for the stellar and supermassive black hole feedback, respectively. For the subset of runs analysed in this work, we used |$\epsilon _{\rm SF,b}=10{{\ \rm per\ cent}}$| and |$\epsilon _{\rm BH,b}=1\div 10{{\ \rm per\ cent}}$| for the magnetic feedback, while the different values for the feedback efficiency (referred to the |$\epsilon _{\rm SF, BH} = \dot{M} c^2$| energy accreted by star forming or black hole particles) are listed in Table 1. Further details can be found in Vazza et al. (2017). The initial seed field in all the runs inset to |$B_0=10^{-9} \rm ~nG$|⁠, with the exception of run CSFB11 which started from |$B_0 = 0.01 ~\rm nG$|⁠. A more detailed overview of the above methods is given in Appendix  A. Clearly, the fixed spatial resolution of our runs is, in principle, too coarse to properly model galaxy formation processes with sufficient detail (even massive galaxies are resolved within a few cells at most). However, in Vazza et al. (2017), we showed how our star-formation recipes (based on Kravtsov 2003) are specifically calibrated to properly reproduce the cosmic star formation history in the volume (see Appendix of Vazza et al. 2017). Moreover, the combination of cooling and stellar/AGN feedback has been tuned so that the scaling relations of galaxy clusters and groups are well reproduced (see Appendix of Vazza et al. 2017). In Vazza, Brüggen & Gheller (2013a) and Vazza et al. (2016) we showed that these models can also fairly describe the innermost density, temperature, and entropy profiles of galaxy clusters when going to a higher resolution. In summary, our previous works suggest that the feedback recipes implemented in our simulations can effectively mimic the large-scale effect and energetic of galaxy formation recipes at a higher resolution; hence, they can be effectively used to model the large-scale effects on thermal and magnetic feedback on the surrounding distribution of filaments.

All the runs adopted the ΛCDM cosmology, with density parameters Ω_BM = 0.0478 (BM representing the Baryonic Matter), Ω_DM = 0.2602 (DM being the Dark Matter), Ω_Λ = 0.692 (Λ, being the cosmological constant) and a Hubble constant H₀ = 67.8 km s⁻¹ Mpc⁻¹ (Planck Collaboration et al. 2016). The initial redshift is z = 38, the spatial resolution is |$83.3 ~\rm kpc/cell$| (comoving), and the constant mass resolution of |$m_{\rm DM}=6.19 \cdot 10^{7}\, \mathrm{M}_{\odot }$| for dark matter particles.

2.1 Filaments identification

Filaments are identified in the data produced by the different simulations according to the procedure presented in Gheller et al. (2015), based on the VisIt data visualization and analysis framework (Childs et al. 2011). Here, we summarize the main steps of such a procedure.

The mass density of the baryonic component is used in order to separate over- from under-dense regions. The Isovolume-based approach, described in Meredith (2004), has been adopted in order to accurately trace the boundaries between the warm/hot phase which characterizes collapsed cosmological structures and the cold under-dense phase typical of voids. Isovolumes are then segmented into distinct objects applying a Connected Components filter. Identified overdense isovolumes undergo a further cleaning procedure in order to eliminate objects that either cannot be classified as filaments or that can be affected by large numerical errors, due, in particular, to the limited spatial resolution. The main steps of the identification procedure are:

• Identification and removal of large clumps. Galaxy clusters are identified using the Isovolume algorithm out of the highest peaks of the matter distribution. They are discarded by removing from the data all cells identified as part of the cluster.

• Filament identification. The Isovolume algorithm is used once more on the residual cells to identify the volumes with:

  \begin{eqnarray*} \delta _{\rm BM} = {\varrho _{\rm BM}\over \varrho _0} \ge a_{\rm fil}, \end{eqnarray*}

  (1)

  where δ_BM is the baryonic matter (BM) over-density, ϱ_BM is the cell’s baryonic mass density, a_fil is a proper threshold, whose value has been set to 1 after extensive experimentation (Gheller et al. 2015), and ϱ₀ is the critical density. The Connected Components algorithm is then used to combine cells belonging to distinct (i.e. non-intersecting) filaments, assigning to each filament an Id (an integer number) and tagging each cell belonging to a filament with that Id. The same identification criterion, but using the total (baryonic matter  + dark matter) density, is used to extract the galactic haloes (see Section 4).
• Shape selection. We use the following two-stage cleaning procedure to remove possible round-shaped, isolated structures that could be still present in the identified objects. First, we retain only elongated objects, with:

  \begin{eqnarray*} {\rm MAX}(r_{xy}, r_{xz}, r_{yz}) \gt \alpha , \end{eqnarray*}

  (2)

  where r_ab is the ratio between axes a and b of the bounding box and α was set to 2. Second, among the remaining objects, we accept as filaments all those whose volume is smaller than a given fraction of the corresponding bounding box volume.

The resulting methodology depends on six parameters whose influence and optimal setting are extensively discussed in Gheller et al. (2015).

The whole procedure is implemented within the VisIt framework, which provides all the required numerical methods and supports big data processing through a combination of optimized algorithms, the exploitation of high-performance computing architectures, in particular through parallel processing, and the support of client-server capabilities.

3 PROPERTIES OF THE FILAMENTS

In Fig. 1, we give the visual impression of the 3D distribution of gas density, temperature, and magnetic field magnitude in a |$10^3 ~\rm Mpc^3$| sub-volume for three reference runs in our sample: the baseline non-radiative primordial run, the sub-grid dynamo model (DYN5), and the cooling + feedback run with astrophysical seeding of magnetic fields (CSFBH3). The overall morphology of thermal gas properties on |$\ge \rm Mpc$| scales is very similar in all runs, even if the gas distribution in the radiative run is significantly clumpier. As generally found also in the outskirts of galaxy clusters (e.g. Nagai & Lau 2011; Vazza et al. 2013b), also in this case, the combination of radiative cooling and feedback heating enhances the level of gas clumpiness ((〈ρ²〉)^0.5/〈ρ〉), at least sufficiently away from the active sources of feedback. The gas density in galaxy groups in the field is generally smoother in the radiative run while their core is hotter, as an effect of ongoing thermal feedback from the SMBH active in it. On the other hand, the distribution of magnetic fields in the volume is dramatically different at all scales, becoming progressively less diffuse and more concentrated on to haloes going from the primordial scenario to the dynamo to the astrophysical one. While in the first case filaments have a floor magnetization of |${\sim } 1{-} 10 ~\rm nG$| resulting from the compression of the primordial magnetic field lines, in the dynamo case they get significantly magnetized only where the flow is turbulent enough to boost a dynamo, which happens only in the proximity of substructures within filaments. Finally, in the astrophysical scenario, the activity of star-formation and SMBH feedback could only magnetize ‘bubbles’ around haloes, leaving most of the filaments volume un-magnetized (down to the |${\sim } 10^{-18} ~\rm G$| floor value assumed here). The variety of possible distributions of magnetic fields across eight different simulations is shown in Fig. 2. In galaxy clusters, most of the models generate magnetic field whose magnitude is of the order of observed values (⁠|${\sim } 0.1{-}1 ~\rm \mu G$|⁠). In the vast majority of the cosmic web, instead, the differences in magnitude and volume filling factors of magnetic fields can be large. We will analyse such differences in detail in the following sections.

Table 2 shows the number of filaments identified in the different models using the procedure described in Section 2.1. Considering all masses, the number of detected filaments is highly variable, due to the influence of small objects strongly affected by the local physics and the resolution effects. Selecting objects with |$M_{\rm BM} \ge 5 \cdot 10^{12} \, \mathrm{M}_{\odot }$|⁠, the number counts become homogeneous, with the exception of the models with strong feedback. The feedback, contrasting the infall of matter into the forming structures, leads to a decrease of the total number of filaments.

3.1 Statistical properties

Similar to our previous works (Gheller et al. 2015, 2016), we studied the scaling relations of global quantities of the filaments in our sample as a function of physical variations. For each filament, we have calculated the total gas mass, the filament length, approximated as the diagonal of the bounding box enclosing the filament as determined by our filament finder, the average baryon fraction within each filament (normalized to the cosmic baryon fraction) and the mass-weighted mean gas temperature and magnetic field for all our identified objects at z = 0. All quantities are calculated considering the cells selected by our finder as part of the filament. Cells at the boundaries are described as variable number of faces polyhedra, accurately adapting to the complex geometry of the filament.

Fig. 3 gives the scaling relations between the total gas mass in filaments, the filament length, the baryon fraction within each filament (normalized to the cosmic baryon fraction) and the mass-weighted mean gas temperature for all our identified objects at z = 0. We also give in Tables 3–5 the best-fit relations between the baryonic mass³ of filaments and length, temperature, and baryon fraction, using a simple log₁₀(Y) = α + β log₁₀(X) relation, where X = M_BM. The same fitting relation, for different quantities X and Y, will be adopted throughout the paper, fitting parameters being always referred as α and β.

The same relations have also been investigated in Gheller et al. (2016). We now extend the analysis investigating the role played by different thermal and non-thermal physical processes on the population of cosmic filaments. The lowest mass filaments in our distributions (⁠|$\le 5 \cdot 10^{12} \, \mathrm{M}_{\odot }$|⁠) typically present a large scatter in all scaling relations, which makes our fitting procedure unstable. The origin of this scatter is twofold. On one hand, the smallest filaments can be affected by numerical effects, i.e. the limited numerical resolution makes it difficult to properly resolve the properties of the smallest filaments, their size getting comparable to the numerical resolution. For these small objects, the gravitational potential is under-resolved, numerical pressure affects the force balance, and the finite DM mass resolution can introduce spurious graininess in the mass distribution (and hence in the gravitational potential). At the same time, small-size filaments are more easily affected by the random effect of cooling and feedback from haloes around and within them, which introduce an additional strong scatter in their thermodynamic properties. For example, run CSFBH5, which presents overall the strongest effect of feedback from star formation, due to the low threshold density used to form stars, is found to produce the largest scatter in the temperature, baryon fraction, and magnetic field strength for |$\le 5 \cdot 10^{12} \, \mathrm{M}_{\odot }$| filaments. This second mechanism has a physical origin and changes with the model. The related dependence of the level of scatter from the baryonic mass, different for each model, may in principle offer an (observationally challenging) way of testing magnetic field scenarios in the radio band, as we discuss in Section  5. Here, however, we want to extract robust estimates of the scaling relations of our objects, in a regime in which the random impact of single stellar and AGN activity is reduced, and therefore in the following we discuss our best-fit estimates restricted to |$M_{\rm BM} \ge 5 \cdot 10^{12} \, \mathrm{M}_{\odot }$| filaments.

Our most important findings, shown in Figs 5 and 6 and listed in Tables 3–7, can be summarized as follows.

First, the impact of model variations on the geometric and baryonic properties of filaments is overall modest, and the filaments population looks remarkably similar across 4–5 orders of magnitude in mass, despite the variety of investigated physical setups. The relation between the length of the filaments and their baryonic mass is always close to |$L \propto M_{\rm BM}^{0.5}$|⁠, consistent with Gheller et al. (2016), while only modest differences in normalization are found between models. As expected, the overall evolution of the large scale structure of the universe is guided by gravity and pressure forces, other local physical processes playing only a minor role. The different physical models analysed in these paper cannot be disentangled through the study of global geometric or thermodynamics properties.

Second, the mean gas temperature of all our objects falls within the expected range for the WHIM: |$10^5 ~\rm K\le T \le 10^{7} ~\rm K$|⁠. In models with cooling and feedback, the temperature of filaments with mass |$M_{\rm BM} \approx 10^{14}\, \mathrm{M}_{\odot }$| is a factor ∼2 higher than in non-radiative runs. At the low mass end of the distribution (⁠|$10^{10}\!-\!10^{11} \, \mathrm{M}_{\odot }$|⁠), cooling and feedback models are a factor ∼2 colder. This leads to a relation between gas temperature and baryonic mass always steeper for cooling and feedback models compared to non-radiative runs (e.g. |$T \propto M_{\rm BM}^{0.35}$| versus |$T \propto M_{\rm BM}^{0.3}$|⁠). Moreover, we measure a significant scatter in the temperature values towards the lowest masses as feedback events from neighbouring AGNs or star-forming regions (either within filaments or next to them) affect the thermodynamics of the WHIM, contributing a significant amount of energy per particle. The long-term effect of cooling is to trigger thermal and magnetic feedback, that can indeed partially affect the WHIM thermodynamics on scales of several |${\sim } \, \rm Mpc$|⁠. The relative trend of temperature in filaments simulated with radiative and non radiative setups is qualitatively similar to our earlier analysis in Gheller et al. (2015). A quantitative comparison is challenging due to the different resolution and feedback scheme (which is here more sophisticated). Comparing our non-radiative runs to those with the same mass and spatial resolution presented in in Gheller et al. (2016), we find that the temperature–mass relation for our baseline model is slightly flatter (⁠|$T \propto M_{\rm BM}^{0.3}$| versus |$T \propto M_{\rm BM}^{0.431}$|⁠). This is connected to the more restrictive mass range we considered here, |$M_{\rm BM} \ge 5 \cdot 10^{12} \, \mathrm{M}_{\odot }$|⁠, which minimizes the large scatter produced by low-mass objects, which are potentially the most affected by numerical effects due to the fixed spatial and force resolution. We refer to Gheller et al. (2016) for further discussion on the influence of numerical resolution on the scaling relations.

Third, the baryon fraction enclosed in filaments displays in all models a flat relation with gas mass, i.e. |$f_{\rm b}/f_{b,\rm cosm}\propto M_{\rm BM}^{(-0.004 \div 0.03)}$| in the |$\ge 5 \cdot 10^{12} \, \mathrm{M}_{\odot }$| mass range. The normalization significantly varies across models, with the tendency of filaments in non-radiative to retain slightly more baryons in filaments than the cosmic average. In the strongest feedback models, on the other hand, the average baryon fraction drops to |${\sim } 80\!-\!90{{\ \rm per\ cent}}$| of the cosmic mean. This can be ascribed to the combined effect of cooling and feedback, that on one hand condenses more gas from filaments into haloes and away from the WHIM phase, and, on the other hand, can push baryons even outside filaments, following strong feedback events (Cen & Ostriker 2006). Once more, at the low-mass end of our distribution the scatter is large, as already found in Gheller et al. (2016), again due to a combination of the more extreme effect of AGN feedback on to small-mass systems (e.g. Oppenheimer & Davé 2006), as well as to the limited mass resolution for dark matter. The number distribution of baryon fraction across our entire sample is shown in Fig. 4. In general, the scatter at high baryon fraction, typically found in runs with stronger feedback effects, suggest that large variations due to the AGN-feedback cycle are an important physical mechanism to vary the baryonic (and possibly chemical) content of cosmic filaments across the cosmic volume (Cen & Ostriker 2006; Haider et al. 2016; Martizzi et al. 2018).

Finally, the magnetic properties of the WHIM in filaments across our sample and for the different models as a function of the baryonic mass show the most striking differences. Figs 5 and 6 present the scaling relations between the filament’s baryonic mass and the corresponding maximum and mean (mass-weighted) magnetic field strength, as well as the number distribution of magnetic field strength in our runs. The mean magnetic field in filaments is roughly constant with mass in the high-mass range (i.e. |$\langle |B| \rangle \propto M_{\rm BM}^{-0.04 \div 0.2}$|⁠) with a normalization dependent on the seeding model. This is because, without efficient small-scale dynamo amplification, the average distribution of magnetic fields in the volume of filaments mainly depends on the input seed field, modulo the compression factor (Donnert et al. 2009; Vazza et al. 2017). The normalization in the primordial scenarios investigated here is ∼10² − 10³ higher than in the dynamo or astrophysical scenarios. However, the magnetization of the largest filaments is similar, i.e. |$\langle |B|\rangle \sim 50\, \rm nG$| , as an effect of the |${\propto} M_{\rm BM}^{0.2}$| scaling relation. On the other hand, for |$M_{\rm BM}\le 10^{12} \, \mathrm{M}_{\odot }$| objects, the differences are large.

The differences in magnetic fields grow in radiative runs. Such differences follow from the trend of mean magnetic fields and matter overdensity (Vazza et al. 2017): the astrophysical seeding scales with the number of feedback sources in the environment (and the number of galaxies strongly depends on the size of filaments, e.g. Gheller et al. 2016). Only for the highest efficiency scenario investigated here (CSF5, with an assumed |$\epsilon _{\rm SF,b}\sim 10{{\ \rm per\ cent}}$| magnetization efficiency from star formation feedback) the magnetization in |$M_{\rm BM} \ge 5\cdot 10^{12} \, \mathrm{M}_{\odot }$| filaments reaches the |$\langle |B|\rangle \sim 10\, \rm nG$| level, while in the same mass range this is |${\sim } 0.1 ~\rm nG$| for the lowest efficiency models. We observe steeper relations and (obviously) larger differences in normalization for the relation between the maximum magnetic field and the host baryonic mass of filaments, which is helpful to highlight the presence of rare sub-regions in filaments where small-scale dynamo amplification can be present at some degree. In this case, we observe |$B_{\rm max} \sim 0.5 ~\rm \mu G$| in the non-radiative scenarios, and up to |${\sim } 1 ~\rm \mu G$| in the strong feedback scenario CSF5.

While the above trends apply to the statistical properties of the filaments population in the cosmic volume, in the following section, we will focus on their internal properties, in a range of scales that may, in principle, be probed by future X-ray or radio observations.

3.2 Resolved properties of filaments

Characterizing the trends of thermodynamic and magnetic properties along the major and minor axis of filaments can yield interesting information on their internal structure, similar to what is routinely done for galaxy clusters. However, filaments have complex shapes, often tilted, folded, and with secondary ‘branches’, that pose a number of practical problems while attempting to identify the ‘spine’ of a filament, making it difficult to compute such profiles for a large number of objects (e.g. Gheller et al. 2015).

We have initially focused on a |$M_{BM} \sim 7.04 \times 10^{13} \, \mathrm{M}_{\odot }$| filament that could be identified in all the different simulations, with a simple-enough geometry that allows us to remove secondary branches. The resulting object, whose isodensity contours are shown in Fig. 7, is suitable for standard profile analysis assuming approximate cylindrical geometry. Fields of interest have been averaged within cylindrical shells at increasing distance from the filament’s spine. The location of the spine is identified computing the centre of total mass of horizontal parallel planes crossing the filament, once this is rotated along the z-axis: X_s = [∑x_i · (m_gas + m_DM)_i]/[(m_gas + m_DM)_i] and Y_s = [∑y_i · (m_gas + m_DM)_i]/[(m_gas + m_DM)_i] for the X_s and Y_s coordinates, respectively.

The calculated profiles are presented in Fig. 8 and are referred to as the height, h, defined as the distance of the cylindrical shell from the filament’s spine on the cutting planes. The cells associated to haloes within the same volume have been removed (haloes appear as white objects embedded in the filament in Fig. 7). Small-scale differences in the profiles are likely an effect of the different time-scales of feedback processes in the different runs, and they are not considered significant. However, we observe the tendency of the central gas density and temperature to be higher in cooling and feedback runs. This is a combined effect of the slightly enhanced compression of gas on to the filament’s spine promoted by cooling with the non-gravitational heating effect by AGN and stellar winds from surrounding haloes, which mostly affect the temperature profile at |$h \ge 1.5 ~\rm Mpc$|⁠. The velocity profiles are remarkably similar across the different models, with a rather constant value comprised between |$v_f\sim 400\!-\!600 ~\rm km\, s^{-1}$| at heights greater than 1 Mpc, and a significant drop towards the filament’s spine. The temperature is uniform from the spine to about 2 Mpc at a mean temperature around 10⁶ K for most of the models. At this temperature, the sound speed in the WHIM is |$c_s \sim 250\!-\!300 ~\rm km\, s^{-1}$|⁠, which means that gas motions internal to the filament are transonic, v_f/c_s ∼ 1.5–2. The magnetic field profile is extremely flat in all models, and ranges from |${\sim } 50 \,\rm nG$| in primordial and dynamo scenarios, to |${\sim } 0.01 ~\rm nG$| in the CSF2 scenario.

Given the monotonic dependency of the gas density from the height h, following Gheller et al. (2016), we compute the median value of the quantities of interest as a function of this variable. In this way, we can bypass the limitations related to the complex geometry of the filaments and extend the analysis to a larger number of objects. Fig. 9 gives the 33–66 per cent range around the median relation of gas temperature, magnetic field strength, velocity module, and baryon fraction for 10 filaments identified in all runs and in the |$M_{\rm BM}=2-5 \cdot 10^{13} \, \mathrm{M}_{\odot }$| mass range. In the case of gas temperature, we used the (〈T〉, M_BM) scaling relations measured in the previous section which are basically identical in this mass range to normalize temperatures for different masses, while we leave the normalization free to change for all other quantities, considering that the average magnetic field and the baryon fraction do not show a clear trend with mass, while the gas overdensity at which filaments form is always the same and shall not be further renormalized.

We find that the differences between models are overall mild, and get more pronounced at high-enough gas particle number density (n/〈n〉 > 50), which mark region close to the filament spine but also close to overdense substructures in the filaments’volume. These substructures are located where the sources of feedback are, hence the differences mostly reflects the different impact of feedback processes to the WHIM surrounding galaxies in filaments (see next section). All runs consistently show that in this mass range the baryon fraction is ∼1–3 times the cosmic average for (2 ≤ n/〈n〉 ≤ 50), regardless of the details of gas physics, while at higher densities the impact of different physical prescriptions is large. The (B, n/〈n〉) relation shows instead larger differences at all overdensities, with a 10² difference at high overdensities, and a ∼10⁴ difference in the low-density range. This trend is consistent with the global (B, n/〈n〉) reported in Vazza et al. (2017), although there, they were referring to the entire distribution of baryons in the cosmic volume, and not specific to filaments.

4 GALACTIC HALOS IN FILAMENTS

Adopting the selection procedure described in Section 2.1, it is possible to extract over-dense objects tracing the population of haloes within filaments. The adopted overdensity selection criterion uses the total (dark + gas matter) density. The overdensity threshold to excise haloes is set to the value of a_gal = 100, which defines haloes undergoing virialization, in line with our previous results in Gheller et al. (2016). The simulations presented in this paper cannot follow the internal properties of forming galaxies in a robust way, in particular at the lowest masses of our halo distribution. Galaxy formation physics, in fact, is missing in the non-radiative runs and, in general, the spatial resolution in all the simulations is not sufficient to trustworthily resolve the innermost structure of galaxies. However, the statistical and physical properties of the hosting haloes are well resolved and described by our approach, as shown in Gheller et al. (2016), where the same procedure has been adopted to compare the resulting halo catalogues with the GAMA data (e.g. Alpaslan et al. 2014), reporting a significant degree of similarity between the properties of the haloes identified in our simulations and real galaxies in the survey. We have also calculated the 3D two-point correlation function of the haloes, presented in Fig. 10, which resulted in following for all models the expected power law with an exponent averaged over the different runs β = −1.839 ± 0.002 in the range 0.1 to 10 Mpc and slightly different normalizations (in arbitrary units), further showing that the properties of our halo catalogues are consistent with those expected for galaxies in the standard ΛCDM model. The similarity of the correlation curves for the different models, indicates also that the clustering properties of the haloes are only mildly influenced by the different physical setup characterizing the various simulations.

4.1 Statistical properties of the haloes in the filaments

The number of haloes identified in the different runs is presented in Table 8. Severely under-sampled objects, with volume smaller than 2³ computational cells, have been removed. The resulting halo masses are typically |$M_{\rm h}\ge 10^{11} \, \mathrm{M}_{\odot }$|⁠, although few objects at smaller masses are found. The number of haloes tends to be larger in models with radiative physics, showing that these processes can favour the process of galaxy formation.

The distributions of the number of haloes with given mass and average temperature (top row), average magnetic field strength and mean velocity (bottom row) are given in Fig. 11. For the mass function above |${\sim } 2\cdot 10^{11} \, \mathrm{M}_{\odot }$|⁠, the number of objects decreases with increasing mass. The distributions are rather similar for the different models. However, models CSF5 and CSFBH5 present a tail with masses around |$10^{14} \, \mathrm{M}_{\odot }$|⁠. In these runs, the density threshold for star formation is five times smaller than in the other runs. This leads to a more efficient star formation, removing gas pressure support from within the haloes, enhancing the cooling of further gas on to the halo volume and, finally, increasing the haloes gas content. The sharp drop in the number counts at masses below |$10^{11} \, \mathrm{M}_{\odot }$| is instead due to the adopted selection criterion, which removes most of the least massive objects. Also, for the temperature, the distributions follow similar patterns in all the models, with average temperatures between 10^4.5 and 10^7.5 K. Even when cooling and feedback are included in the simulation, they tend to self-regulate, leading to average temperature whose values are close to those of runs without any cooling or heating active. Models CSF5 and CSFBH5 show their peculiar behaviour also in the temperature distribution, in particular, at the highest temperatures. Once more, the highly efficient feedback processes influence the haloes properties, leading to an increase in the energy input in the gas, hence to objects with temperature higher than in all the other models. Differences among the runs can be found also at the lowest temperatures. The primordial and dynamo models have haloes in a smaller temperature range compared to the astrophysical models, their thermodynamics being controlled exclusively by adiabatic compression and shock waves, that heat all the gas at temperatures between 10⁵ and 10⁷ K. The models CSF2 and CSFB11 are among the least efficient star feedback models and they do not include feedback from supermassive black holes, hence the associated thermal feedback is smaller than for the other astrophysical models. At the same time, cooling is active, decreasing the temperature of the shocked gas. This results in overall colder objects compared to all the other runs.

The mean magnetic field distributions show remarkable differences between the simulations. The two primordial models have similar trends, with the highest values of the magnetic field peaked between 10⁻² and |$10^{-1} \, \mu$|G and with a narrow distribution. The three dynamo models have a similar shape, with peak between 10⁻³ and |$10^{-2}\, \mu$|G and a broader distribution, encompassing magnetic fields from 10⁻⁴ to slightly more than |$10^{-1}\, \mu$|G. When cooling and feedback are active, the distributions result to be even broader, with a tail of few objects with magnetic field B|$\lt 10^{-7}\, \mu$|G. The number of haloes grows with B until reaching a maximum which is different for the various models and then it declines, sharply for the models CSF5 and CSFBH5 and more gently in the other models. The CSFB11 model is an exception: due to the high primordial magnetic seed and the low star formation efficiency, its behaviour is close to that of the dynamo models at a slightly lower range of |B|. We remark that, while these trends are perfectly in line with those presented for the filaments, for haloes the lack of resolution and physical details in the description of such overdense environment hampers any significant small-scale dynamo amplification. Our flows within haloes in fact, have a magnetic Reynolds number R_m ≪ 100 (defined as the ratio between magnetic induction terms and magnetic diffusivity, i.e. U · L/η_m, where U is the characteristic flow velocity on the scale L and η_m is the magnetic diffusivity), which severely prevents the formation of the small-scale dynamo, at variance with what is commonly accepted to happen in real galaxies (e.g. Beck et al. 2013; Schleicher et al. 2013; Schober, Schleicher & Klessen 2013; Rieder & Teyssier 2016; Marinacci et al. 2017).

Finally, the velocity distribution is found to behave similarly in all models for most of the haloes. As expected, the different physics in the various models does not affect the overall kinematic properties of the haloes. At low velocity values large differences are present, but the low number statistics does not allow reaching firm conclusions.

Halo properties, in particular masses, can be correlated to those of the hosting filaments. The dependence of the total mass, the length, the average temperature and magnetic field of the filaments from the total mass of embedded haloes is shown in Fig. 12. All distributions have been fitted with the same power law introduced above. The mass and length distributions become steeper at masses above |$M^*_{TOT,h} = 10^{13}\, \mathrm{M}_{\odot }$|⁠, hence fitting curves have been calculated considering both all the masses and only objects with mass above |$M_{TOT,h}^*$|⁠. Since the trends are very similar in all models, average fitting parameters could be calculated. The results are presented in Table 9.

The mass of the filaments tends to grow almost linearly with that of the resident haloes. The trend is even closer to linearity when only high total halo masses are considered. The length of the filaments scales approximately as |$L_{\rm fil} \propto M_{\rm TOT,h}^{1/3}$| including all objects, while it tends to scale as |$\propto M_{\rm TOT,h}^{1/2}$| restricting to higher masses, due to the tendency of bigger filaments to be thinner and more elongated than the smaller ones. For the mean temperature, no changes in β with the halo mass are found, with a trend given by the power law |$T_{\rm fil} \propto M_{\rm TOT,h}^{1/3}$| at all scales. For the magnetic field results change significantly with the class of models, with the same trends we already observed in Fig. 5.

The above scaling are consistent with what was previously found with non-radiative simulations in Gheller et al. (2016), confirming that the total halo mass is trustworthy and effective proxy of the host filament’s properties.

4.2 Halos and their local environment

The distribution of galaxies naturally provides a clue of how baryonic matter is distributed in the cosmic web (e.g. Alpaslan et al. 2014). The overdense and relatively hot environment surrounding galaxies and the associated haloes in filaments is a natural target for observations. The characterization of the properties of such environment is therefore meaningful in order to study the correlation between galaxies, the local gas distribution, and the overall properties of the hosting filaments.

In order to identify the galaxies’ local environment, we have extracted spherical regions centred on each halo, with radius ten times that of the halo itself. Fig. 13 shows the relation of the properties of the resulting environment volumes and the haloes lying within. In general, the properties within and outside of identified haloes are correlated, albeit with a large scatter, exhibiting trends that, once more, have been fit as |$\rm log_{\rm 10}(X_e) = \alpha + \beta ~ \rm log_{\rm 10}(X_h)$|⁠, where X_e and X_h are various average quantities calculated for the environment and for the corresponding halo, respectively.

In all models, the average temperature of the external WHIM (top left panel) scales with the mean halo temperature, with almost the same exponent, despite the large scatter. Its value, averaged among the different runs, is β = 0.894 ± 0.029. The average normalization is α = 0.547 ± 0.175. Halos result to be at slightly higher temperature that the surrounding medium, this difference slowly growing with the halo mass. Similar considerations can be made for the velocity (bottom left panel), with a comparable spread of the points distribution around the average best-fit curve, with β = 0.742 ± 0.057 and α = 0.620 ± 0.146. The only exception is the DYN8 model, which has a fitting curve remarkably different from all the others, showing that, if sufficiently strong, the magnetic field can affect the dynamical properties of the gas (e.g. Marinacci & Vogelsberger 2016). In general, these trends show that haloes move consistently with the embedding environment. Once more, the magnetic field presents peculiar signatures associated with the different physical properties of the various models, which affects both the normalization and the scaling. The fitting parameters α and β are presented in Table 10. In primordial and dynamo models the magnetic field intensity of both haloes and their surrounding environment is comprised between 0.0001 and 0.2μG. Halos with the highest magnetic field are embedded in environments with approximately the same value of |B|. For smaller values of the haloes magnetic field, the environment tends to have slightly bigger values of |B|, as confirmed by the slope of the fitting curves, always less than one. This trend is even more evident in astrophysical models, that typically encompass a broader range of |B|, with both slopes and normalization significantly lower than that of primordial and dynamo models. However, haloes in astrophysical models at the lowest values of the magnetic field lies in environments with an even lower magnetisation.

Finally, the bottom right panel of Fig. 13 compares the baryonic fraction of the haloes to that of their neighbourhood. The enclosed baryonic fraction in the haloes is referred to the hot/warm gas component, while we do not consider here the mass locked into stars and SMBH particles. In all the models, the baryon fraction within haloes is, as expected, lower than the cosmic average, with values around |${\sim } 10-50 {{\ \rm per\ cent}}$| of such average. These values are typically higher than those for galaxy groups (e.g. Giodini et al. 2009; Andreon 2010) and galaxies (e.g. Hoekstra et al. 2005; Papastergis et al. 2012; Zaritsky et al. 2014) reported in literature. This is the effect of the limited spatial resolution combined to our halo identification procedure, which lead to extract volume significantly larger than those typical for those kind of objects, including matter at higher baryon fraction from the local environment, which, in fact, is found to be at baryon fraction f_b/f_{b, cosm} ≈ 1, significantly higher than that of the halo but also of the overall filament. The gas matter expelled from the halo enriches the surrounding environment for all models and at all masses (e.g. Oppenheimer & Davé 2006), with a mild trend of the baryonic fraction to increase with the mass for non-feedback models, since no cooling mechanisms are present capable of decreasing the gas pressure in the halo, allowing it to fall in the forming object. The scatter around the fitting curves is large. However, the environmental baryon fraction almost never gets smaller than f_b ≈ 0.5 f_{b, cosm}.

In Fig. 14, we analyse the correlation between the halo masses and the average density (top left panel), temperature (top right panel), magnetic field (bottom left panel) and baryon fraction (bottom right panel) of the embedding environment. For all the models the distributions are highly scattered, especially at low masses. The dispersion tends to decrease for masses above |$10^{13}\, \mathrm{M}_{\odot }$|⁠, where the fitting curves are calculated. Once more, for all the models the temperature follows similar power laws. The average fitting parameters (calculated as above) result to be β_T = 0.454 ± 0.005, α_T = 0.431 ± 0.077. At the highest masses this relation provides a precise estimate of the temperature dependence from the halo mass. At lower masses, due to large scatter, it gives a lower limit of the environmental temperature for a given halo mass. The density does not show a well defined trend, and in this case the lower limit is imposed by our filament identification procedure. A slight trend of the environmental average density to increase with the halo mass is identifiable, but, due to the large dispersion a clear behaviour cannot be assessed. The magnetic field shows a clear increasing trend with the halo mass, which is sharper in astrophysical seeding models compared to primordial models, which also present larger normalization. The baryon fraction presents a similar behaviour, its value slightly growing with the hosted halo mass. Once more, model CSFBH5 has a much larger scatter of values, and in general features the lowest baryon fraction values in the external medium, at all given halo masses, due to the extreme feedback strength here (which also explains the large values of magnetic fields measured external to haloes). In this case, the feedback is capable of expelling a significant fraction of baryons even outside of the volume of filaments, as already noticed in Fig. 3.

5 DISCUSSION

While we defer to future work the detailed analysis of the observational signatures in the X-ray and radio bands of filaments in different physical scenarios, some general trends can be anticipated based on these new results. Due to the wide spectrum of models available in the Chronos++ suite, in this work we could extensively compare the effect of different physical processes on the thermodynamics and the magnetic properties of the gas in cosmological filaments. Based on our extraction procedure, we could separate three different types of target objects: filaments as a whole, galactic haloes within filaments, and the local environment where haloes lies.

For all the presented correlations, well defined trends emerge for the most massive objects (i.e. |$M^*_{BM} \ge 5\cdot 10^{12} \, \mathrm{M}_{\odot }$| for filaments and |$M^*_{TOT,h} \ge 10^{13} \, \mathrm{M}_{\odot }$| for the haloes).

For the mass–temperature relation, we found only mild differences among different models, indicating that the final overall thermodynamic status of the gas is comparable in all cases, even if the local effect of heating and cooling can be considerably different, as shown by the broad dispersion in the environmental temperature correlated with the mass of the corresponding haloes. This is true considering, both, the dependency of the filament temperature from the total baryonic mass of the filament and from the total mass of hosted haloes. The slope of the various models tends to decrease for increasing total filament’s gas mass. It remains instead almost uniform at all scales when halo masses are considered. Small mass filaments show indeed a large scatter in their thermodynamic properties, due to the feedback from galaxies and AGN around them. These objects will be challenging to observe both in the soft X-ray band, having |$T \le 10^{6} ~\rm K$|⁠, and in the radio band, having |$\langle B \rangle \sim 10^{-3}\div 1 ~\rm nG$| (depending on the model). Due to the large scatter of scaling relations in these objects, robustly constraining the thermal and magnetic properties of the WHIM in such filaments may be possible only with statistics of hundreds of objects.

The thermal properties of the largest filament (⁠|$M_{\rm BM} \ge 5 \cdot 10^{12} \, \mathrm{M}_{\odot }$|⁠) do not allow an easy discrimination between the different scenarios considered in this work. Such discrimination would be even more challenging for observational data, due to their inherent large uncertainties. The only relevant difference between model are found in the high density range (n/〈n〉 ∼ 10²), i.e. close to the central spine of filaments and in the proximity of haloes, where feedback effects are more effective. On the other hand, no systematic differences among models are found for the majority of the volume filled by WHIM gas.

Significant differences among models are found when the magnetic field is considered. In primordial models, where the only mechanism of amplification of the magnetic field is compression, the normalization of the initial seed field has to be set high enough to get the observed level of magnetisation in galaxy clusters. At current time, this results in objects with magnetic field of the order of |$10^{-1}\, \mu$|G at all scales. Halos have a slightly higher values of |B|, while for filaments such value tends to mildly grow with the baryonic mass. Filaments from dynamo models show a stronger dependency from the mass, with values of |$|B| \approx 10^{-2}\div 10^{-3}\, \mu$|G at low masses, reaching values around |$10^{-1}\, \mu$|G at gas masses above |$10^{13} \, \mathrm{M}_{\odot }$|⁠. This follows from the efficiency of the dynamo mechanisms, which increases with the overdensity. A similar dependency from the mass is shown from the magnetic field in astrophysical models. At masses below |$M^*_{BM}$|⁠, |B| is always smaller than |$10^{-2}\, \mu$|G, with values below |$10^{-4}\, \mu$|G reached by various models. Above |$M^*_{BM}$|⁠, no more dependence from the mass is shown, and only the CSF5 models reaches values for |B| close to |$10^{-1}\, \mu$|G. The distinction between dynamo and astrophysical runs is larger when the local environment is considered, only the CSFBH5 model showing some overlap with DYN7. For embedded haloes with masses above |$M^*_{TOT,h}$|⁠, in all models the dispersion in the values of |B| is small. Therefore, the local volume hosting massive haloes appears to be a promising target environment to characterise the mechanisms driving the evolution of magnetic fields in the cosmic web.

Our results suggest that the magnetic properties of the WHIM can effectively shed light on the mechanisms behind the evolution of cosmic filaments. For instance, the combination of the trend of magnetic fields and volume with filaments’ mass found in our sample, indicates that the difference in Faraday Rotation (⁠|$\rm RM \propto L_{\rm LOS} ~n_e ~B$|⁠, where L_LOS is the length of the line of sight crossing the filament volume) across the filament population varies by orders of magnitude, if extragalactic magnetic fields are mostly contributed by galaxy-formation related processes (i.e. in the astrophysical scenario). Conversely, the scatter across the population is expected to the of order ∼10, if magnetic fields have a primordial origin. In this case, in fact, the effect of compression is rather similar across filaments of very different masses. This addresses a viable test of magnetogenesis with filaments, provided that a large number of objects can be detected in the future using the SKA (Vazza et al. 2018a). Several tens of detected RM are in fact required to remove the contamination by RM intrinsic to sources (e.g. Locatelli, Vazza & Domínguez-Fernández 2018). While these techniques rely on the signature by magnetic fields to detect the WHIM, we expect that with the SKA it should be possible to attempt also a detection of the coldest gas (⁠|${\sim } 10^4 ~\rm K$|⁠) in the cosmic web using the HI hyperfine transition, even if the required sensitivity and statistics represent an even bigger challenge than that for synchrotron and Faraday Rotation studies (e.g. Popping et al. 2015; Horii et al. 2017).

Depending on their true magnetisation level, a significant role of filaments in the production of anisotropies in the propagation of UHECRs in the local Universe can be expected (Sigl et al. 2003; Dolag et al. 2003; Hackstein et al. 2016, 2018) even if a level of |${\ge} 10 \rm ~ nG$| appears necessary (e.g. Kim et al. 2019).

The baryon fraction results to be higher than in galaxies for most of the filaments, but in general less than the cosmic value in all models at almost all scales, except for the largest masses, where f_b/f_{b, cosm} is of the order of the unity or even bigger. In low mass objects, we measure a |${\sim } 5{-}10 {{\ \rm per\ cent}}$| significantly larger baryon fraction in low-mass filaments (⁠|$M_{\rm BM}\le 10^{12} \, \mathrm{M}_{\odot }$|⁠) simulated with cooling and feedback compared to non-radiative simulations. This highlights the potential role of AGN feedback and star winds from surrounding galaxies in expelling gas from their haloes on to filaments (e.g. Cen & Ostriker 2006; Haider et al. 2016). In this case, also metal-rich gas is expected to be deposited in filaments, enriching their composition and boosting their line emission and absorption signatures (e.g. Oppenheimer & Davé 2006; Martizzi et al. 2018). Critical to further explore these features with numerical simulation is the access to even higher resolution and more complex chemical networks (e.g. Hummels et al. 2018; Popping et al. 2019). The analysis of the local volumes hosting haloes shows, despite a large scatter, the tendency of the gas surrounding galaxies to have a large baryon fraction, often significantly above the cosmic value. On the contrary, the vast majority of the haloes have 0.1 ≤ f_b/f_{b, cosm} ≤ 0.5. Differently from the whole filament, no meaningful differences are found between the different models. Irrespective to the details of the galaxy formation physics, the gas in the surrounding environments is fed by the matter ejected from the corresponding haloes, confirming the expectation of metal-enrichment discussed above.

In agreement with Gheller et al. (2016), we have found several scaling laws that relates the mass, the length and the temperature of a filament with the properties of the embedded haloes. These scaling laws have small dispersion with slightly different trends at low and high (⁠|$\gt M^*_{TOT,h}$|⁠) masses, the latter having a higher slope. These relations proved to be substantially the same in all models, providing a general tool to infer the overall properties of a filaments starting from those of galaxies, useful aid in their detection and identification. In addition, a sharp lower limit in the temperature of the WHIM surrounding haloes of a given mass has been identified. Also this information can be exploited to constrain the properties and help in the detection of the filament’s gas pointing to the neighbourhood of the resident haloes.

6 CONCLUSIONS

Cosmic filaments are a fundamental building block of the cosmic web (e.g. Bond et al. 1996; Davé et al. 2001; Cautun et al. 2014; Martizzi et al. 2018) and potentially outstanding instruments to access the past evolution of cosmic gas, preserving traces of the fundamental processes which took place several |$\rm Gyr$| in the past, from the enrichment of complex chemical elements and metals (e.g. Oppenheimer & Davé 2006; Smith et al. 2011; Nicastro 2016), to the injection and evolution of non-thermal particles in large-scale structures(e.g. Pfrommer et al. 2007; Vazza, Gheller & Brüggen 2014b), to the seeding and dynamo amplification of extragalactic magnetic fields (e.g. Ryu et al. 2008; Vazza et al. 2014a; Marinacci et al. 2015).

In this new work, due to a large statistics of cosmological simulation, we accomplished an extensive survey the thermal and non-thermal properties of filaments in the cosmic web.

In the following summary, we list those characteristics that can be considered robust predictions of our simulations. We start from those properties that show only a weak dependency on the investigated variations of baryonic physics or magnetic field seeding scenario:

• the global structural properties of filaments (in particular their length, mass and average temperature);

• the global baryon fraction enclosed within the filaments volume;

• the average profiles of gas and dark matter, temperature and baryon fraction with respect to the central axis of filaments in the |$\ge 10^{13} \, \mathrm{M}_{\odot }$| mass range;

• the mass, temperature and velocity distribution of galaxy haloes in filaments, as well as the clustering properties, as shown by the 3D 2-points correlation function of haloes within filaments;

• the temperature and velocity of the haloes, which are consistent with the environment in which they live. In particular, the velocity of the haloes traces the velocity of the fluid;

• the baryon fraction, for which different values, significantly smaller than the unity in haloes, of the order of the cosmic value or even larger for the local halo’s environment and slightly smaller than one in the rest of the filament, can be identified.

Our study highlighted also several properties which appear to significantly depend on the assumed baryonic physics or magnetic physical models:

• the |$T \propto M_{fil,BM}^{\rm \alpha }$| relation is steeper in radiative runs with feedback, compared to non-radiative runs (i.e. α ≈ 0.35 versus α ≈ 0.3), and the most massive filaments are on average a factor ∼2 hotter if feedback from galaxies and AGN is active. Differences between models are smaller when the total halo mass (instead of the total baryonic mass) is considered. In this case the relation appears to be bi-modal with the slope becoming higher above |$10^{13} \, \mathrm{M}_{\odot }$|⁠;

• the relation between the mass and the mean magnetic fields of filaments hugely varies in slope and normalization if different seeding scenarios are assumed: |$B_M \propto M_{fil,BM}^{0.15}$| in primordial models, |$B_M \propto M^{0.4}_{fil,BM}$| in dynamo models and |$B_M \propto M^{1.8}_{fil,BM}$| in astrophysical seeding models. Similar differences are found considering the halo mass. Galactic haloes present a larger scatter in their magnetic fields in astrophysical seeding models, as well as larger variations with respect to their environment;

• the velocity profile and the magnetic field profile with respect to the central axis of filaments in the |$\ge 10^{13} \, \mathrm{M}_{\odot }$| mass range is different in the different scenarios. Feedback produces a flatter velocity profile out to |${\sim } 4 ~\rm Mpc$| from the axis of filaments due to outflows, while the magnetic field profile in astrophysical seeding models significantly drops for ≥2 Mpc from the filaments’ central axis, due to the scarcity of magnetisation sources.

While this study attempted an exhaustive study of the large-scale properties of filaments in the cosmic web, based on the most complete to date suite of physical variations of WHIM, a natural (yet challenging) follow-up is represented by the access to a higher spatial and mass resolution, besides the adoption of more sophisticated and accurate description of the physics underlying the evolution of cosmic filaments and galactic haloes. This will allow us to include and probe galaxy formation processes, providing a direct link between galaxies and the large scale distribution of the gas.

Furthermore, our current results will be extended in order to investigate in details how the properties of filaments maps to observables, in particular in radio, potentially detectable by the coming generations of radiotelescopes, focusing on synchrotron and HI emission and Faraday Rotation measurements. In addition other interesting phenomena, like the potential impact of filaments on UHECRs emission, can be subject of additional in-depth analysis.

ACKNOWLEDGEMENTS

The cosmological simulations were performed with the enzo code (http://enzo-project.org), which is the product of a collaborative effort of scientists at many universities and national laboratories. We gratefully acknowledge the enzo development group for providing extremely helpful and well-maintained on-line documentation and tutorials. F.V. acknowledges financial support from the ERC Starting Grant ‘MAGCOW’, no. 714196. The simulations on which this work is based have been produced on Piz Daint supercomputer at CSCS-ETHZ (Lugano, Switzerland) under projects s701 and s805 and on the Jülich Supercomputing Centre (JFZ) under project HHH42, in numerical projects with F.V. as PI. We also acknowledge the usage of online storage tools kindly provided by the Inaf Astronomica Archive (IA2) initiave (http://www.ia2.inaf.it). We thank Marcus Brüggen for his support in the first production of the Chronos++ suite of simulations employed in this work, and to Jean Favre for his support in the first design of the Visit filament finder employed in this work.

Footnotes

REFERENCES

APPENDIX A: SUB-GRID MODELS FOR GALAXY FORMATION PHYSICS

Our runs adopt the Kravtsov (2003) star formation model implemented in enzo, with a few ad-hoc modifications to better perform in large unigrid simulations. Star-forming particles (actually stellar populations) in enzo are formed on the fly whenever specific gas conditions are met in cells (Bryan et al. 2014). Following Kravtsov (2003), four criteria must be met to form a star particle: a) the gas in the cell must exceed a threshold density, n ≥ n_*; b) the cell is part of a converging flow (⁠|$\nabla \cdot \vec{v} \lt 0$|⁠); c) the local cooling time is smaller than the dynamical timescale (t_cool ≤ t_*); d) the baryonic mass is larger than the minimum star mass (m_b = ρΔx³ ≥ m_*). When all these conditions are met in one cell, enzo forms at run-time a stellar particle with a mass m_* = m_bΔt/t_*, where Δt is the time-step.

Stars can also return thermal energy, gas mass and metals back into the surrounding gas. The feedback from each star particle (coming from supernovae explosions) depends on the assumed fractions of energy/momentum/mass ejected per each formed star particles, E_SN = ϵ_SFm_*c². In our runs and considering the resolution at which we work, the feedback energy is entirely released as thermal (i.e. hot supernovae-driven winds), which generally promotes pressure-driven winds around stellar particles (v_wind ∼ 10–10² km s⁻¹). Given the relatively large time steps in our simulations (⁠|$\Delta t \approx 5{-}10 ~ \rm Myr$| at low redshift) we decided to model the feedback effect in the same time-step in which each stellar particle form, in an instantaneous recycling assumption, which bypasses unnecessary particle-to grid loops that kills the performances for large unigrid simulations. The Kravtsov (2003) model is suitable for uniform grid simulations and has been designed to reproduce the observed Kennicutt’s law (Kennicutt 1998). In Vazza et al. (2017) we have further shown how in our runs it can well reproduce the observed cosmic star formation history, compiled from infrared and ultraviolet observations by Madau & Dickinson (2014), if an appropriate set of model parameters is set (see Table 1).

In the CHRONOS++ suite, we also modified enzo to support the injection of the magnetic fields by stars, introducing magnetic dipoles at each feedback episode. The dipoles are randomly oriented along one of the simulation axis. The total injected energy scales with the feedback thermal energy, i.e. E_{b, SN} = ϵ_{SF, b} · E_SN. For the subset of simulations presented in this work, we used |$\epsilon _{\rm SF,b}~\sim 1{-}10 {{\ \rm per\ cent}}$| (see Table 1). We notice that the combination of spatial and mass resolution (Δx = 83.3 kpc and |$m_{\rm DM}=6.19 \cdot 10^{7}\, \mathrm{M}_{\odot }$|⁠, respectively) should quench the formation of galaxies at the low mass end |$\le 10^8 \, \mathrm{M}_{\odot }$|⁠, and therefore our runs still lacks a contributor to the cosmic star formation rate (e.g. Genel et al. 2014). However, the contribution to both the chemical, thermal and magnetic enrichment of the intergalactic medium by the implemented feedback from |${\sim } 10^{10}{-}10^{11} \, \mathrm{M}_{\odot }$| galaxies is well resolved and is expected to dominate in the full cosmic volume (e.g. Samui, Subramanian & Srianand 2017).

The injection, growth, and feedback from supermassive black holes (SMBH) on the surrounding cosmic gas has been included as well, implemented with few modifications on top of existing enzo routines, based on Kim et al. (2011). Unlike star-forming particles, SMBH particles in enzo must be seeded at a given redshift (we chose z = 4 in all runs) with an initial mass (⁠|$M_{\rm BH}=10^4 \, \mathrm{M}_{\odot }$| in our case), and from that moment on, they are advected at run time (similar to dark matter particles as well as star particles, they are assumed to be collisionless and only subject to the local gravity acceleration) grow in mass and produce feedback events. We assumed for each SMBH a Bondi-Hoyle accretion rate, with the additional model assumptions (motivated by the low resolution), as in Table 1: a) the gas accretes on to SMBH assuming a fixed temperature of |$3 \cdot 10^5 \rm K$| ; b) the Bondi accretion rate gets boosted by a factor α_Bondi = 10²–10³ if the gas density which can be resolved around each SMBH is too low.

SMBH particles release thermal energy to the surrounding gas with an efficiency ϵ_BH = ΔMc²Δt/E_BH, ΔM being the accreted gas mass and E_BH being the feedback energy. We have implemented the release of magnetic feedback energy from SMBH in enzo in the form of magnetic dipoles, assuming E_{b, AGN} = ϵ_{b, AGN}E_BH. In this paper we explored values in the range |$\epsilon _{\rm b,AGN}=1{-}10 {{\ \rm per\ cent}}$|⁠. Our tests in Vazza et al. (2017) show that such treatment of SMBH and associated feedback (combined with star formation) is reasonably capable of modifying the mass-temperature scaling relations of galaxy groups and clusters, in the |$10^{13} \, \mathrm{M}_{\odot } \le M \le 10^{15} \, \mathrm{M}_{\odot }$| mass range, towards X-ray observations (Reichert et al. 2011; Eckmiller, Hudson & Reiprich 2011), by increasing the typical temperature of |$M_{\rm 500} \le 10^{14}\, \mathrm{M}_{\odot }$| systems in a realistic way.

APPENDIX B: SUB-GRID MODELS FOR MAGNETIC DYNAMO AMPLIFICATION

Our run-time implementation of a sub-grid (SG) model for small-scale magnetic dynamo in enzo attempts to compute the maximal energy that can be channelled into magnetic fields, based on the dissipation of solenoidal turbulent motions that we can measure at run-time in the simulation. This is done in order to bracket, at least in as set of ad-hoc models, the maximum effect of turbulent dynamo amplification on scales that our simulation cannot resolve, and that will probably remain challenging to resolve for a long time (e.g. Donnert et al. 2018). Our extended resolution studies on magnetic amplification in filaments (Vazza et al. 2014a) tend to exclude the possibility of a significant development of a volume filling dynamo in filaments, induced by the decay of fast and supersonic turbulent motions measured there. However, our simple MHD view cannot exclude that in real plasmas ‘microscopic’ small-scale instabilities arising from kinetic plasma effects (e.g. Mogavero & Schekochihin 2014, for a review) can further amplify the field, starting from sub-kpc scales, even if there are presently no observational reason to think this is the case. The maximum energy that can be channelled into magnetic fields during a root-grid timestep, Δt (⁠|${\sim } 10 ~\rm Myr$|⁠) is of the order of |$E_{\rm B,dyn} = \epsilon _{\rm dyn}(\mathcal {M})F_{\rm turb} \Delta t$|⁠, where |$F_{\rm turb} \simeq \eta _t \rho \epsilon _{\omega }^3/L$| is the turbulent kinetic energy flux, with |$\epsilon _{\omega }=(\nabla \times \vec{v})^2$| is the flow enstrophy which we measure at run-time on a L = 3 cells stencil, while η_t is the estimated dissipation rate on to magnetic fields, which must be guessed or calibrated based on other simulations (Jones et al. 2011; Porter et al. 2015; Beresnyak & Miniati 2016; Vazza et al. 2017; Wittor et al. 2017a).

In our runs, we adopted a general approach in which η_t can scale with the local plasma parameters, as its value can critically change going from the sub-sonic to the supersonic regime (Ryu et al. 2008), based on Federrath et al. (2014), who simulated small-scale dynamo in a variety of conditions. We could estimate the saturation level and the typical growth time of magnetic fields as a function of the local turbulent Mach number of the flow, |$\mathcal {M_{\rm turb}}$|⁠. In particular, |$\epsilon _{\rm dyn}=\epsilon _{\rm d}(\mathcal {M_{\rm turb}})$| can be approximated as |$\epsilon _d(\mathcal {M}_{\rm turb}) \approx (E_B/E_k) \Gamma \Delta t$|⁠, where E_B/E_k is the estimated ratio between magnetic and kinetic energy at saturation, and Γ is the typical growth rate, which we both take from Federrath et al. (2014) as a function of the flow Mach number. Once we have computed the amplified magnetic energy, E_{B, dyn}, we can generate at run-time a magnetic field vector, |$\vec{\delta B}$|⁠, to be added to the already existing field, imposing it to be parallel to the local direction of the gas vorticity, so that the new generated field is also solenoidal by construction. A corresponding budget in kinetic energy is removed from cells. Momentum is subtracted assuming an isotropic dissipation of the small-scale velocity vectors. Our procedure is manifestly simpler than more sophisticated SG models on the market (Grete et al. 2016) and the resulting topology of magnetic fields is not consistently attached to the underlying turbulent flow. However, the overall energy budget in amplified magnetic fields well matches the prediction derived in other works (e.g. Ryu et al. 2008) and is overall suitable to bracket the (possible) amplification of magnetic fields in the WHIM under extreme conditions. In our CHRONOS++ suite (Vazza et al. 2017) and in the simulation presented here (see Table 1) we explored a few variations in the amplitude of seed fields before the dynamo acts, as well as in the normalization of the |$\epsilon _d(\mathcal {M}_{\rm turb})$| dynamo efficiency, finding only mild differences.