The large-scale environment from cosmological simulations II: The redshift evolution and distributions of baryons

TL;DR

This study examines how large-scale environments and baryon distributions evolve with redshift in cosmological simulations, confirming that baryons minimally impact environment classification and that the gas web remains an unbiased matter tracer.

Contribution

It extends previous analysis to higher redshifts, analyzing baryon effects and the distribution of the warm-hot intergalactic medium across different large-scale environments.

Findings

Baryon models have negligible effect on large-scale environment classification.

The gas web is an unbiased tracer of total matter, especially at high redshifts.

Approximately 40% of gas is in the warm-hot intergalactic medium, decreasing with redshift.

Abstract

Following Cui et al. 2018 (hereafter Paper I) on the classification of large-scale environments (LSE) at z = 0, we push our analysis to higher redshifts and study the evolution of LSE and the baryon distributions in them. Our aim is to investigate how baryons affect the LSE as a function of redshift. In agreement with Paper I, the baryon models have negligible effect on the LSE over all investigated redshifts. We further validate the conclusion obtained in Paper I that the gas web is an unbiased tracer of total matter -- even better at high redshifts. By separating the gas mainly by temperature, we find that about 40 per cent of gas is in the so-called warm-hot intergalactic medium (WHIM). This fraction of gas mass in the WHIM decreases with redshift, especially from z = 1 (29 per cent) to z = 2.1 (10 per cent). By separating the whole WHIM gas mass into the four large-scale environments (i.e. voids, sheets, filaments, and knots), we find that about half of the WHIM gas is located in filaments. Although the total gas mass in WHIM decreases with redshift, the WHIM mass fractions in the different LSE seem unchanged.

Tables

Table 1. Table 1: The cosmological gas and stellar mass fractions in percentiles with respect to the total baryon mass. The fractions from the AGN run are shown in the first followed by the results in brackets from the CSF run. The results from Nicastro et al. ( 2018 ) are at z < 0.5. We further note here that the WHIM fraction from Nicastro et al. ( 2018 ) has a slightly different upper temperature limit (10 6.2 versus 10 7 ). The fractions from Haider et al. ( 2016 ) are coming from the Illustris simulation (Vogelsberger et al., 2014 ) . The same method introduced in Section 4.1.2 is used to estimate these error bars.

Redshift	hot gas	WHIM	cold gas	star
Nicastro et al. (2018)a	9$\pm$4.5	$\gtrsim$24 & $\lesssim$55	29.7$\pm$11	7$\pm$2
Haider et al. (2016)b	6.5	53.9	32.8	-
z = 0	4.6$\pm$0.7 (4.6$\pm$0.1)	41.3$\pm$1.1 (38.3$\pm$1.0)	50.9$\pm$0.1 (50.6$\pm$0.2)	3.2$\pm$0.1 (6.5$\pm$0.2)
z = 0.6	2.4$\pm$0.4 (1.1$\pm$0.2)	34.9$\pm$1.1 (29.7$\pm$1.1)	60.2$\pm$0.1 (65.0$\pm$0.1)	2.5$\pm$0.1 (4.2$\pm$0.2)
z = 1.0	1.3$\pm$0.2 (0.3$\pm$0.0)	28.7$\pm$1.1 (21.9$\pm$1.0)	68.1$\pm$0.1 (74.8$\pm$0.1)	1.9$\pm$0.1 (2.8$\pm$0.1)
z = 2.1	0.2$\pm$0.0 (0.0$\pm$0.0)	10.8$\pm$0.5 (6.2$\pm$0.4)	88.2$\pm$0.1 (92.9$\pm$0.0)	0.8$\pm$0.0 (0.8$\pm$0.0)

Table 2. Table 2: The distribution of different baryon components in the LSE at the four redshifts. The AGN results are shown by the first value with the CSF results follow inside the brackets. These fractions are normalised to the total mass of each baryon component, i.e. the sum of each line equals 1. See Section 4.1.2 for details of how these error bars are calculated.

	Voids	Sheets	Filaments	Knots
		z = 0
$f_{\mathrm{total}\mathrm{gas}}$	0.15$\pm$0.02 (0.16$\pm$0.02)	0.38$\pm$0.01 (0.39$\pm$0.01)	0.37$\pm$0.02 (0.36$\pm$0.02)	0.10$\pm$0.01 (0.09$\pm$0.01)
$f_{\mathrm{star}}$	0.08$\pm$0.01 (0.06$\pm$0.01)	0.33$\pm$0.02 (0.29$\pm$0.02)	0.44$\pm$0.01 (0.48$\pm$0.01)	0.15$\pm$0.02 (0.17$\pm$0.02)
$f_{\mathrm{hot}\mathrm{gas}}$	0.00$\pm$0.00 (0.00$\pm$0.00)	0.04$\pm$0.01 (0.05$\pm$0.01)	0.46$\pm$0.03 (0.50$\pm$0.03)	0.49$\pm$0.03 (0.45$\pm$0.03)
$f_{\mathrm{WHIM}}$	0.05$\pm$0.01 (0.06$\pm$0.01)	0.30$\pm$0.02 (0.31$\pm$0.02)	0.51$\pm$0.01 (0.50$\pm$0.01)	0.14$\pm$0.01 (0.13$\pm$0.01)
$f_{\mathrm{cold}\mathrm{gas}}$	0.25$\pm$0.02 (0.25$\pm$0.02)	0.48$\pm$0.01 (0.48$\pm$0.01)	0.24$\pm$0.02 (0.25$\pm$0.02)	0.03$\pm$0.01 (0.03$\pm$0.01)
		z = 0.6
$f_{\mathrm{total}\mathrm{gas}}$	0.14$\pm$0.02 (0.14$\pm$0.02)	0.39$\pm$0.01 (0.39$\pm$0.01)	0.37$\pm$0.02 (0.36$\pm$0.02)	0.11$\pm$0.01 (0.11$\pm$0.01)
$f_{\mathrm{star}}$	0.05$\pm$0.01 (0.03$\pm$0.00)	0.28$\pm$0.02 (0.22$\pm$0.02)	0.47$\pm$0.01 (0.48$\pm$0.01)	0.21$\pm$0.02 (0.27$\pm$0.02)
$f_{\mathrm{hot}\mathrm{gas}}$	0.00$\pm$0.00 (0.00$\pm$0.00)	0.01$\pm$0.00 (0.00$\pm$0.00)	0.23$\pm$0.01 (0.15$\pm$0.02)	0.76$\pm$0.02 (0.85$\pm$0.02)
$f_{\mathrm{WHIM}}$	0.03$\pm$0.00 (0.03$\pm$0.00)	0.25$\pm$0.02 (0.23$\pm$0.02)	0.53$\pm$0.01 (0.52$\pm$0.01)	0.20$\pm$0.01 (0.22$\pm$0.01)
$f_{\mathrm{cold}\mathrm{gas}}$	0.21$\pm$0.02 (0.19$\pm$0.02)	0.48$\pm$0.01 (0.47$\pm$0.01)	0.28$\pm$0.02 (0.29$\pm$0.02)	0.04$\pm$0.01 (0.04$\pm$0.01)
		z = 1.0
$f_{\mathrm{total}\mathrm{gas}}$	0.14$\pm$0.02 (0.14$\pm$0.02)	0.40$\pm$0.01 (0.40$\pm$0.01)	0.36$\pm$0.02 (0.35$\pm$0.02)	0.10$\pm$0.01 (0.10$\pm$0.01)
$f_{\mathrm{star}}$	0.03$\pm$0.01 (0.02$\pm$0.00)	0.26$\pm$0.02 (0.20$\pm$0.02)	0.48$\pm$0.01 (0.48$\pm$0.01)	0.22$\pm$0.02 (0.29$\pm$0.02)
$f_{\mathrm{hot}\mathrm{gas}}$	0.00$\pm$0.00 (0.00$\pm$0.00)	0.02$\pm$0.01 (0.00$\pm$0.00)	0.22$\pm$0.02 (0.08$\pm$0.02)	0.76$\pm$0.02 (0.91$\pm$0.02)
$f_{\mathrm{WHIM}}$	0.02$\pm$0.00 (0.02$\pm$0.00)	0.24$\pm$0.01 (0.21$\pm$0.01)	0.53$\pm$0.01 (0.52$\pm$0.01)	0.21$\pm$0.01 (0.25$\pm$0.01)
$f_{\mathrm{cold}\mathrm{gas}}$	0.20$\pm$0.02 (0.18$\pm$0.02)	0.48$\pm$0.01 (0.46$\pm$0.01)	0.28$\pm$0.02 (0.31$\pm$0.02)	0.04$\pm$0.01 (0.05$\pm$0.01)
		z = 2.1
$f_{\mathrm{total}\mathrm{gas}}$	0.18$\pm$0.02 (0.18$\pm$0.02)	0.45$\pm$0.01 (0.45$\pm$0.01)	0.31$\pm$0.02 (0.31$\pm$0.02)	0.06$\pm$0.01 (0.06$\pm$0.01)
$f_{\mathrm{star}}$	0.02$\pm$0.00 (0.02$\pm$0.00)	0.23$\pm$0.02 (0.20$\pm$0.01)	0.52$\pm$0.01 (0.50$\pm$0.01)	0.24$\pm$0.02 (0.28$\pm$0.02)
$f_{\mathrm{hot}\mathrm{gas}}$	0.00$\pm$0.00 (0.00$\pm$0.00)	0.07$\pm$0.01 (0.00$\pm$0.00)	0.35$\pm$0.02 (0.04$\pm$0.03)	0.58$\pm$0.03 (0.96$\pm$0.03)
$f_{\mathrm{WHIM}}$	0.02$\pm$0.00 (0.01$\pm$0.00)	0.24$\pm$0.01 (0.18$\pm$0.01)	0.58$\pm$0.01 (0.53$\pm$0.01)	0.23$\pm$0.01 (0.28$\pm$0.02)
$f_{\mathrm{cold}\mathrm{gas}}$	0.20$\pm$0.02 (0.19$\pm$0.02)	0.47$\pm$0.01 (0.47$\pm$0.01)	0.28$\pm$0.02 (0.30$\pm$0.02)	0.04$\pm$0.01 (0.05$\pm$0.01)

Equations

$${\Sigma }_{\alpha \beta }=-\frac{1}{2H_0}\left(\frac{\partial v_{\alpha }}{\partial r_{\beta }}+\frac{\partial v_{\beta }}{\partial r_{\alpha }}\right),$$

Full text

††pubyear: 2019††pagerange: The large-scale environment from cosmological simulations II: The redshift evolution and distributions of baryons–References

1 Introduction

In the early Universe, matter is thought to have been distributed almost homogeneously, with very small fluctuations in the density field. These “initial condtions” are primarily evidenced by the smoothness of the CMB and the power spectrum of its inhomogenities amounting to roughly one part in 105. Due to the nature of gravitational instability, small over-densities in the matter distribution are able to overcome the cosmological expansion and collapse to form potential wells, that may virialize into so-called dark matter halos (see White & Rees, 1978, for example). If the matter budget of the Universe is dominated by Cold Dark Matter (CDM), then smaller dark matter halos attract each other and merge together to form larger halos, growing according to a hierarchical scenario (e.g. Davis et al., 1985). Therefore, the Universe’s largest structures are dynamically young, having either completed their merging recently or are, in fact, still in the process of growing (for example Planelles et al., 2015).

There are a number of ways to qualitatively and quantitatively characterise the large-scale environment (LSE) of galaxies and dark matter haloes (e.g. Stoica et al. 2005; Hahn et al. 2007; Zhang et al. 2009; Sousbie 2011; Hoffman et al. 2012; Cautun et al. 2013; Leclercq et al. 2016 and see Libeskind et al. 2018 for a review of these). One class of methods is based on the construction of tensors which characterise the tidal field, specifically the Hessian of the potential (Pweb) or the velocity shear (Vweb). Such methods segment space into four dynamically distinct environments: knots, filaments, sheets, and voids. Following Hahn et al. (2007) and Hoffman et al. (2012), we applied both methods to classify the LSE in Cui et al. (2018) (hereafter Paper I), which showed consistent results between the two methods.

Pushing our previous analysis (Paper I) into high redshifts, we are able to present a full picture of how the LSE forms and evolves. Further, the next generation of telescopes, such as the Square-Kilometre Array111https://www.skatelescope.org/(SKA), will be able to probe HI galaxies out to redshift $z\approx 2$ (Abdalla et al., 2015). The Hot Universe Baryon Surveyor222http://heat.tsinghua.edu.cn/~hubs/en/index.html (HUBS) is designed to probe the warm-hot intergalactic medium (WHIM) through a soft X-ray telescope and will be able to map the warm-hot gas on large scales. In Paper I, we found that the gas web is an unbiased tracer of the dark matter web. Is this statement still valid at high redshift? And how are the baryons distributed in these LSE?

In the literature, Cautun et al. (2014) studied the evolution of the cosmic web with the NEXUS code. Using a dark-matter-only simulation, they found that the cosmos is dominated by tenuous filaments and sheets at early times. During subsequent evolution, after the first small structures merge together, the present-day web is dominated by fewer, but much more massive, structures. Using hydro-simulations with cooling, but no star forming nor feedback schemes, Zhu & Feng (2017) investigated the evolution of the cosmic web. They found that sheets appear early; filaments supersede sheets as the primary collapsing structures from $z\sim 2-3$. They further found that there is relatively less baryonic matter residing in filaments and knots compared with dark matter. In the recent work of Martizzi et al. (2018), they applied the tidal tensor method to the Illustris-TNG simulation (see Springel et al., 2018, for example) and found the WHIM to be the dominant baryon mass contribution in filaments and knots at redshift $z=0$, but declines faster at higher redshifts ($z>1$). Here we are going to revisit this issue, however, using a velocity based technique (i.e. Vweb) to define the LSE. Similarly, Nevalainen et al. (2015) found that the WHIM is usually associated with the filamentary structures($\sim$70 per cent of the gas in WHIM phase), which is identified with the Bisous model Tempel et al. (2014).

This paper is organised as follows: we first briefly present the suite of hydro-simulations as well as the sub-grid models (§ 2), and then describe the cosmological structure classification method used for this work – the Vweb (§ 3). In Section 4 we present our findings on the evolution of these large-scale structures as well as the baryon distributions. Finally, we summarise our conclusions and comment on the results in § 5.

2 The cosmological simulations

The same series of three cosmological simulations presented in Paper I are used in this work. These simulations, which share the same initial conditions with box-size $410~{}h^{-1}~{}{\rm Mpc}$ and $2\times 1024^{3}$ particles, use a flat $\Lambda$CDM cosmology, with cosmological parameters of $\Omega_{\rm m}=0.24$ for the matter density parameter, $\Omega_{\rm b}=0.0413$ for the baryon contribution, $\sigma_{8}=0.8$ for the power spectrum normalisation, $n_{\rm s}=0.96$ for the primordial spectral index, and $h=0.73$ for the Hubble parameter in units of $100~{}h~{}{\rm km}~{}s^{-1}~{}{\rm Mpc}^{-1}$. The masses of gas and DM particles have a ratio such that to reproduce the cosmic baryon fraction, with $m_{g}=7.36\times 10^{8}~{}h^{-1}\>{\rm M_{\odot}}$ and $m_{DM}=3.54\times 10^{9}~{}h^{-1}\>{\rm M_{\odot}}$, respectively. Gravitational forces have been computed using a Plummer-equivalent softening length which is fixed to $\epsilon_{PL}=7.5h^{-1}$ physical kpc from z = 0 to 2, and fixed in comoving units at higher redshift. The simulations were run with the TreePM-SPH code GADGET-3, an improved version of the public GADGET-2 code (Springel, 2005). We refer to the dark-matter-only simulation as the DM run. The two hydrodynamical simulations both include radiative cooling, star formation and kinetic feedback from supernovae. In one case, we ignore feedback from AGN (which is referred to as the CSF run), while in the other we include it (which is referred to as the AGN run). Radiative cooling for gas is computed for non-zero metallicity using the cooling tables by Sutherland & Dopita (1993) with an evolving ultraviolet background. Star formation and supernova feedback with chemical enrichment are fully described by Tornatore et al. (2007). The model of AGN feedback used in the AGN run is the same as that adopted by Fabjan et al. (2010), which has optimised the parameters for galaxy cluster simulations. Interested readers are referred to Cui et al. (2012); Cui et al. (2014) for details on the baryon models of these simulations, and to Cui et al. (2016, 2017) for details on the statistical samples of galaxy clusters from them. In this work, we only focus on four distinct redshifts, i.e. $z=[2.1,1.0,0.6,0.0]$. We understand that our resolution is much lower than the one employed in the latest hydrodynamical simulations. However, this set of simulations have a much larger volume which suits for our LSE studies.

3 Web classification method – Vweb

Following Paper I, we apply the shear tensor technique Vweb to classify the LSE of the three simulations in this work.

The velocity shear tensor is used to classify the cosmological space into different environments. Following Hoffman et al. (2012), this tensor is defined as

[TABLE]

where, $H_{0}$ is the Hubble constant. The eigenvalues of $\Sigma_{\alpha\beta}$ are denoted as $\lambda_{i}$ ($i$ = 1, 2 and 3).

In practice, the computation of the eigenvalues for the matrix is performed on regular $256^{3}$ grids, corresponding to a cell size of comoving $\sim 1.6~{}h^{-1}~{}{\rm Mpc}$. The velocity of all kinds of particles is assigned to the nearby mesh points using the triangular-shaped cloud (TSC) method. Further, a top-hat smoothing is applied to these cells. For the particular analysis of gas webs, we only use the velocities from gas particles. Then, the velocity shear tensor is calculated for each grid cell by taking derivatives from nearby cells. At last, the eigenvalues of $\Sigma_{\alpha\beta}$ for each cell are given by matrix elements calculated in the previous step. The velocity or density field estimated with different methods (as discussed in Hahn et al., 2015; Cui et al., 2008, respectively) may lead to unacceptable levels of noise and biases. However, the large cell size and additional smoothing should be able to provide us with a robust velocity field in the real space, thereby a reliable shear tensor.

As shown in the Paper I, this cell size does not play an important role. The relative mass and volume fraction differences between the three investigated numbers of cells ($128^{3},256^{3},512^{3}$) are basically within 10 per cent (see Martizzi et al., 2018, for similar results). We also note here that the mesh size should larger than the mean size of the objects found in the simulation to avoid the interiors.

Each individual cell is then classified as either ‘void’, ‘sheet’, ‘filament’, or ‘knot’ according to the eigenvalues $\lambda_{1}>\lambda_{2}>\lambda_{3}$:

void, if $\lambda_{1}<\lambda_{th}$,

• 2.

sheet, if $\lambda_{1}\geq\lambda_{th}>\lambda_{2}$,

• 3.

filament, if $\lambda_{2}\geq\lambda_{th}>\lambda_{3}$,

• 4.

knot, if $\lambda_{3}\geq\lambda_{th}$,

where $\lambda_{th}$ is a free threshold parameter (Hoffman et al., 2012; Libeskind et al., 2012, 2013). Following the discussion of Carlesi et al. (2014), we set $\lambda_{th}=0.1$, which presents better agreement to the visualised density field. As the gas and dark matter generally share the same velocity field, the same threshold ($\lambda_{th}=0.1$) is used for the gas component in the Vweb code. Whether $\lambda_{th}$ depends on redshift or not is still unclear. Zhu & Feng (2017) tried different redshift evolution of the $\lambda_{th}$(z) and argued that the results are in quantitatively agreement with a fixed $\lambda_{th}$. Similarly Libeskind et al. (2014) showed that the directions of the shear tensor eigen-vectors are only weakly dependent on redshift. Therefore, we adopt the same value of $\lambda_{th}$ at all our investigated redshifts for simplicity.

4 Results

4.1 Redshift evolution of the LSE

We first illustrate slices of projected LSE classified by the Vweb method in Fig. 1. Each slice has a thickness of one cell size ($\sim 1.6~{}h^{-1}~{}{\rm Mpc}$). Each cell in this projection plot corresponds to one pixel of the image. The three different runs (DM, CSF, and AGN) show an almost identical classification of the LSE at each redshift. This is not surprising as baryon processes normally leave more influences on scales smaller than the cell size used here to classify the cosmic web. While the result at $z=0$ already are in very good agreement (see Paper I, for details), we expect even fewer differences between the three runs at higher redshift because the mass distribution is less clustered and hence more similar. It is clear from the density evolution (seen in the bottom row of Fig. 1) that the material is condensing into filament and knot environments. However, qualitatively these LSE show less evolution visually (cf. three upper panels). This could be due to the fact that our classification scheme is less sensitive to spatial distributions (see the volume fraction evolution in subsection 4.1.2). Similar results have been found by Martizzi et al. (2018) who used the tidal tensor method. We note here again that the same threshold ($\lambda_{th}=0.1$ for Vweb) is applied over all redshifts.

4.1.1 The evolution of the $\lambda$s

In order to better quantify differences in the LSE and its redshift evolution, we present the evolution of the distributions of the three eigenvalues from the Vweb method, which are directly used for classifying the LSE. Note, the evolution of these three $\lambda$s is directly responsible for any changes in the classification of the LSE. Actually, previous theoretical works (see Doroshkevich, 1989, for example) in linear scales has predicted that the $\lambda_{1}$ corresponds to the formation of the "pancakes" in the high-z, the three direction contraction corresponds to the smallest eigenvalue, $\lambda_{3}$ with the maximum value of $\lambda_{3}$ determines the condensation of the clusters. Furthermore, the divergence of the velocity field is proportional to the sum of the three eigenvalues in linear regimes (e.g. Bernardeau & van de Weygaert, 1996). The evolution of the three eigenvalues directly reflect the changes in the velocity field. The resulting probability distribution functions (PDFs) can be viewed in Fig. 2. The three simulations runs, which are shown by different line styles have almost the same PDFs for the three $\lambda$s at the different redshifts. In agreement with the results shown in Fig. 1, the baryon models have almost no influence on the LSE, i.e. the PDFs for DM, CSF, and AGN are practically indistinguishable (see the residuals in the lower row of Fig. 2). Further, the overall shape and the peak position of these PDFs for all three $\lambda$s show a very weak redshift evolution, apart from the tails, where the maximum values of $\lambda_{1}$ and $\lambda_{2}$ increase as redshift decreases. But there is almost no change for the minimum values. Because the classification of void is solely depending on $\lambda_{1}$ (i.e. $\lambda_{1}<\lambda_{th}$ for void), the left panel clearly indicates a less significant redshift evolution from redshift $z=2.1$ for the volume fractions of the LSE (see also Martizzi et al., 2018) with a negligible contribution from the baryon models. The width of the PDF for $\lambda_{3}$ gets broader as redshift decreases, i.e. the maximum value increases and the minimum value decreases. This indicates that more matter is flowing into the knot environment.

4.1.2 The evolution of the mass and volume fractions

In Fig. 3 we further explore the redshift evolution of the LSE by showing the mass (left panel) and volume fraction (right panel) of each web type (i.e. knot, filament, sheet, void) as a function of redshift. The error bars are estimated as this: the whole simulation box is equally separated into 8 small boxes with size of 205$~{}h^{-1}~{}{\rm Mpc}$; these fractions (e.g. mass fraction, volume fraction) are calculated individually in each sub-box; the error of each fraction is the standard deviation of the 8 sub-box.

For the knots environment, its volume fraction is always around 2 per cent, while its mass fraction increases towards $\sim 10$ per cent at $z=0$. The volume fraction for the filament environment remains unchanged at the level of $\sim 20$ per cent with a slightly lower fraction at $z=0$. On the other hand, the filament’s mass fraction increases from about 30 per cent at $z=2.1$ to $\sim$40 per cent at $z=0.0$. As for the sheet environment, both its corresponding mass and volume fractions decrease with decreasing redshift. As expected from theoretical studies, the void environment occupies most of the space, but looses its mass as matter flows from low-density to high-density regions. The overall evolution of the LSE we obtain here is similar to the findings of Cautun et al. (2014); Zhu & Feng (2017); Martizzi et al. (2018) who used different classification methods and different values of the threshold $\lambda_{th}$. We do not expect an exact match to their results because these fractions precisely depend on the classification method, smoothing scales, etc. (see more discussions in Libeskind et al., 2018).

Even tough it is already clear from the top panels that these fractions are very close to each other when considering the three DM, CSF, and AGN runs, we further quantify the relative differences between the two hydro runs and the DM run in the lower panels by dividing the fractions for each web type and in each hydro model by the value for that web type in the DM run. At $z=0$, both mass and volume fractions in the knot environment show the largest disagreement, around 2 per cent, which is negligible. The difference in other environments at $z=0$ and the fraction differences in all environments at $z>0$ are within 1 per cent. This is mainly because knot environments are the densest and are thus most susceptible to the dynamics of the baryons.

In agreement with previous figures, the baryons leave almost no effect on both fractions over all investigated redshifts. We note here that the error bars in these fractions are generally very small, only a few percents for the largest value. These fraction differences are even smaller at higher redshifts, which can be understood as the accumulation of baryon processes and the matter distribution is identical at the initial redshift for the three simulations. Furthermore, as redshift increases, the fixed co-moving scale examined here ($\sim 1.6~{}h^{-1}~{}{\rm Mpc}$) represents a quasi-linear length scale, especially with respect to the density field. This can be seen by examining the rms of the density field computed on this scale at $z=2$ and at $z=0$ ($\delta_{\rm rms}\approx 0.8$ and $\delta_{\rm rms}\approx 3.0$, respectively, see Fig. 4 of Libeskind et al. 2014).

4.1.3 Consistency between the dark-matter and gas cosmic web

One important finding of Paper I is that gas filaments are an unbiased tracer of the large-scale structures – at least at redshift $z=0$ (see the right column of Fig. 5 in there). Therefore, it is interesting to ask whether that conclusion is still valid at high redshift or not. We now extend our previous analysis to higher redshifts and show in Fig. 4 the relative difference for the mass (left-hand side panel), volume (middle panel) and consistency (right-hand side panel) fractions between the results from the gas-only component and the results from total matter. The consistency fraction is defined as the number of cells that are identified as the same web type in both the gas-only and total matter results divided by the the total number of cells in that web environment found by considering the total matter result. This consistency level best quantifies how the two results agree with each other.

In the left panel, the ratio (the mass fraction classified by the gas component divided by the mass fraction computed from the total matter) for void and sheet environments is always greater than unity (from $\sim$ 1 per cent to 5 per cent, depending on the redshift). This could be because the gas is more smoothly distributed at semi-linear scale (such as the cell size used in this paper) than the collision-less dark matter. Gas follows dark matter to drop into its potential. However, it has additional supports from the pressure to slow down its clumpiness. Therefore, if only the gas component is used to classify the LSE instead of the total, which is dominated by dark matter, more gas tends to be assigned into void and sheet environments, leaving a lower amount allocated in the knot and filament environment. Consistently, the volume fraction is also tiny higher in the void and sheet environments. Nevertheless, it is worth to note that these fractions are very small, which is almost negligible. The agreement (within 1 per cent) for the mass fractions in filaments, sheets and voids becomes better at high redshift, but not for the mass fraction in knots where the relative difference is at about 3.5 - 5 per cent. For all ratios, the error bars are very small which indicates these results are robust.

Unlike the ratio for the mass fraction, the volume ratios in the middle panel show no clear redshift evolution. Besides the volume fraction ratio in the knot environment (from the CSF run at z = 0 and from the AGN run at higher redshift), the difference is much smaller than the mass fraction differences, within about 2 per cent.

The consistency fraction in the right panel is always smaller than 1, which means that the same LSE are identical with both classifications from gas-only component and total matter. It is clear that there is a better agreement at higher redshift for all LSE. Again, the knot environment shows the largest disagreement compared with the other environments. This is because baryon physics plays a more important role in denser environments and on smaller scales. This is especially true for the AGN run, because the gas distribution in knot environments is largely affected by AGN feedback.

Interestingly, for the filament environment, both mass and volume difference fractions are very small (within about $\sim 2$ per cent), which is similar to the $z$=0 results. Over 95 per cent of the cells are cross identified in the gas-only filaments and total matter filaments. These fractions and the consistency ratio show better agreement at higher redshift. Therefore, we can conclude that our previous finding, i.e. that the gas component is an unbiased tracer of the large-scale structure, as expected, is still valid for all redshifts – with an even better agreement at the highest redshift $z$=2.1.

4.2 The evolution and distribution of baryons

The baryons, which occupy less than 5 per cent of the total energy content in the Universe, are the only component of the Universe which can be directly observed. The gas component, which occupies about 90 per cent of total baryons (Nicastro et al., 2018), can be roughly separated into: 1) cold gas (T < 105 K) in diffused or condensed environments; 2) warm-hot gas (105 < T < 107) in the intergalactic medium (WHIM); 3) hot gas (T > 107) mostly in halos (e.g. Davé et al., 2001; Cen & Ostriker, 2006; Tornatore et al., 2010; Planelles et al., 2018). It is believed that the WHIM, which is very hard to be detected, plays a central role to solve the missing baryon problem, namely the mismatch (about one-tenth) in the apparent baryon content between observations at low and high redshift (Persic & Salucci, 1992; Fukugita et al., 1998; Cen & Ostriker, 1999; Bregman, 2007; Shull et al., 2012; Nicastro et al., 2018; Fang, 2018, etc.). Recent works by de Graaff et al. (2017); Tanimura et al. (2019) have applied different methods to search for the missing baryons and found that the WHIM tends to be located in the filaments and between pairs of merging systems. Therefore, in this Section, we detail the evolution and the distribution of baryons and hope to shed some light on this issue from the theoretical point of view.

We first investigate the cosmological baryon fractions in our two hydrodynamical simulations and compare them with the observational results from Nicastro et al. (2018), which are presented in Table 1. Our simulations give a similar WHIM mass fraction as the observations, but the cold gas fractions are a bit higher. The effect of the AGN feedback does not only reduce the total stellar mass by almost half of the CSF run (e.g. Planelles et al., 2013), but it also heats a certain amount of cold gas into the warm/hot phase. The difference of the warm/hot gas mass fraction between the our two simulations is much larger at the two middle redshifts. Therefore, it results in a larger difference in the cold gas fraction as well. This indicates that the AGN feedback plays a more important role (more powerful) at higher redshifts than at z = 0, when most of the AGN feedback is seized due to the galaxies being quenched in the simulation.

It is worth noting that the stellar mass fraction formed in the CSF run is comparable to Nicastro et al. (2018). However, the simulation with AGN feedback, which is believed to be necessary to reproduce the correct stellar mass function, has much lower stellar mass fraction than that seen in observation. This contradiction found in the AGN run, hence such theoretical models (and their numerical implementation), implies that the subgrid baryon models are still challenging.

Compared to the results from Haider et al. (2016), their WHIM fraction is about 12 per cent higher than ours with a much lower mass fraction ($\sim$18 per cent) in cold gas. Note that their gas mass fractions are normalised to the total gas mass instead of total baryonic mass. Further, it is known that the AGN feedback in Illustris is even powerful that gives much more hot gas. Similarly, Martizzi et al. (2018) also report that the mass in the WHIM in Illustris is $\sim$23% larger than in IllustrisTNG.

4.2.1 The evolution of the gas density

We now examine the evolution and distribution of the gas in the simulation as well as in different LSE by showing its density PDF in Fig.5. Instead of directly using the density of each gas particle in the simulation, we recalculate it inside each cell.333We use the same grid for the density calculation as already used for the cosmic web classification with Vweb. As it is well known, the distribution of density in the Universe is expected to be close to log-normal (see Coles & Jones, 1991, and references therein). The PDF of the log of normalized gas density therefore has its peak around 0 with a narrow width at $z=2.1$. As redshift decreases down to zero, the PDF profile gets broader with the peak shifting to smaller values. This can be easily understood as the gas moves from low density regions (voids) to higher dense regions (knots), which will make the void regions even emptier (decreased min density value) and the knot region more compact (increased max density value). This behaviour is well known at least since the seminal work of Coles & Jones (1991) (although see Klypin et al., 2018, for a modern view of the distribution of cosmic density), yet it is worth noting that the two hydrodynamical simulations share the same PDF distributions for low gas density. There is a very small discrepancy (slightly higher density in the AGN run than in the CSF run) at the higher density tail especially at low redshift. That could be caused by the AGN feedback stopping the star formation, leaving more gas inside the halos.

It is not surprising to see the gas with highest density tends to be located in the knot environment. From the filament to the sheet and the void environments, the PDFs of the gas density gradually shift from higher density to lower density with the most under-dense gas located within the void region. Similar to the PDFs of total gas density, the curve of these PDFs in different LSE also have the shape of a log normal distribution with smaller widths at high redshift and broader widths at low redshift. The redshift evolution of these PDFs in different LSE is generally following the total PDFs at the corresponding density range.

Seeing from the lower panel, the AGN feedback seems to leave no influences on the low density gas from the total. The CSF run tends to have less high dense gas compared to the AGN run. This can be understood as the AGN feedback stops the star formation and has more gas left without condensed to stars. It is not surprising to see that the gas PDF in each environment basically follows the total PDFs. It is worth to note that the CSF run tends to have relatively more low-dense gas in the knot and filament environments than the AGN run. This is consistent with the previous explanation as the knot and filament environments can be strongly affected by the AGN feedback than sheet and void environments.

4.2.2 The evolution of the gas temperature

The distribution and evolution of the gas temperature, which is calculated as the mass-weighted mean in the cells, are presented in Fig. 6. We only show the PDFs from the AGN run in upper row with the relative difference in the lower row. Similar to the total gas density PDFs, the total gas temperature distribution evolves along the redshift with a decrease of the peak ($\sim 10^{4}$[k] at $z$=2,1) and an increase of both high and low temperature gas. In knot environment, we find a relatively larger spread of temperature than for the density due to the co-existence of the cold and hot gas in these regions. The gas temperature distribution in the filament environment is also rather wide. However, the low temperature gas starts to contribute more. In both sheet and void regions, the cold gas is the dominant component. There are very few cells with higher temperature $>10^{6}$ K, especially at higher redshift.

Unlike the PDF for the gas density in Fig. 5, there are clear differences between the CSF and the AGN run. The AGN run tends to systemically have higher temperature gas (about a factor of 2 for $T>10^{4}$ with redshift dependence) than the CSF run. While the CSF run exceeds the AGN run by orders of magnitude for the gas at temperature of $\sim 1000$ K. Not surprisingly, the same level of disagreement between the two runs is also showing in each of our LSE. This indicates that the AGN feedback is more efficient in changing the gas temperature than the actual gas distribution.

4.2.3 The evolution of the gas phase diagram

We present the gas phase diagram at different redshifts in the four panels of Fig. 7. The gas density – temperature distribution is separated into 4 different regions by the black vertical and horizontal lines. The vertical lines indicate the limitation of gas density for star formation. Star forming gas particles are assigned to the cold gas phase ignoring their temperature. The larger distribution of the temperature for the star forming gas is because the gas is treated as multi-phase, so as to provide a sub-resolution description of the interstellar medium. For our considered regions below the gas density limitation, the gas is simply separated by its temperature into three parts: hot gas, WHIM and cold gas. It is clear that the cold gas dominates the mass fraction at $z=2.1$. As redshift decreases, more and more gas is heated up into the WHIM and hot region. Detailed fractions can be found in Table 1. There are noticeable differences between the two runs, especially at high redshifts and in the low temperature region. This could be due to the fact that cooling in the CSF runs is more efficient without the feedback from AGN. It is interesting to see that the contours are in good agreement with the colour-map at $z=0$. This means that the AGN feedback is again more powerful at higher redshift and becomes negligible after the galaxies are quenched.

4.2.4 The distribution of baryons in different environments

We now investigate how the gas phase diagram behaves in different LSE (shown by columns) and at different redshifts (shown by rows) in Fig. 8. Note that the colour is normalised to the total gas mass to be consistent with Fig. 7. Again, the colour-map is coming from the AGN run whereas the contours are for the CSF. The gas is separated into different parts by black lines. We refer to Table 2 for detailed fractions in each LSE and redshift.

At $z=0$, voids contain mainly cold gas; sheets are also dominated by cold gas, but there is a similar amount of WHIM; filaments have the largest portion of WHIM with a small fraction of cold gas; knots host the largest amount of hot gas with a small fraction of WHIM. As redshift increases, there is less and less WHIM and hot gas in the void as indicated by the colour map shrinking towards the lower temperature region. A similar conclusion can be drawn for the sheet environment for which the WHIM is also visibly decreasing, especially from $z=1.0$ to $z=2.1$. The WHIM dominates filament environments at $z=0$, but its total amount also decreases with an increase of the cold gas fraction, especially from redshift $z=1$ to $z=2.1$. In the knot environment, the majority of gas remains at the warm/hot phase with a temperature of about 107 K until $z=1.0$, whereas it is taken over by the cold gas at $z=2.1$. This is because both the total WHIM and hot gas decrease very fast from redshift $z=1$ to $z=2.1$ (see Table 1 for details). Besides the difference at the lower temperature region between the two runs, which is already discussed in the previous subsection, the effect of AGN feedback is negligible.

We further detail the fractions of total gas, hot gas, WHIM, cold gas and stars by separating them into different LSE in Table 2. For each gas component and stars, the fraction is normalised to their total mass in that baryon component, i.e. the sum of the fractions in the four LSE is equal to unity. At $z=0$, it is clear that the majority of the gas can be found in sheet and filament environments. The hot gas is mainly divided between filament and knot environments. About half of the WHIM gas is located in filaments, followed by sheets ($\sim$30 per cent). About half of the cold gas is in sheets, the other half is almost equally assigned to filaments and voids. Most stars are in filaments (about 44 per cent), sheets (33 per cent) and knots (15 per cent).

From low to high redshift, more and more (total) gas is found in sheets. This is the essence of structure formation: gas flows from voids to sheets then filaments and finally nodes. Therefore, towards lower redshifts, more and more gas flow into filaments and nodes. The stellar fraction increases mainly in filaments when moving backwards in time, indicative of star formation primarily taking place in filaments at redshifts for which the star formation rate peaks. The hot gas tends to shift from filament towards knot where feedback is prominent. The fractions of WHIM seem stable over redshift with always over half of them in filament. Similar to the WHIM, the fractions for cold gas are also independent of redshift with about half of them in sheets.

At last, there are some pronounced differences between the two hydrodynamical simulations. While they overall give similar fractions, the hot gas fractions at high redshift differ substantially. There is noticeable difference at $z=2.1$, where 96 (58) per cent of hot gas is in knots from the CSF (AGN) run. This is again due to the powerful AGN feedback, which can heat up the gas and push it outside of the halos.

5 Discussion and Conclusions

Using a series of cosmological simulations, which include baryonic processes with different levels of complexity, we study the distribution of matter in the large-scale environments (LSE) at four distinct redshifts. These sets of simulations allow us to explore how baryons affect the classification of the LSE: 1) a dark-matter-only run as a gauge, 2) a simulation that includes gas cooling, star formation and supernova feedback, and 3) a run that additionally models AGN feedback. To analyse all these simulations we employ the velocity shear tensor (Vweb) based cosmic web classification on a regular grid of co-moving spacing $\sim 1.6~{}h^{-1}~{}{\rm Mpc}$. Each cell is assigned a tag ‘knot’, ‘filament’, ‘sheet’, or ‘void’ depending on how many eigenvalues are above a certain threshold $\lambda_{th}=0.1$, which serves as the web classifier.

In this paper, we extend the analysis presented in Paper I to higher redshifts $z>0$ and quantify the detailed baryon distributions. Our main results are summarised as follows

In agreement with the result in Paper I, the baryon models have negligible effect on the LSE classification. Actually, the difference is smaller at redshifts $z>0$. 2. 2.

Using the gas component alone to classify the LSE gives consistent results to the ones obtained accounting for all the matter. Again, the difference tends to be smaller at higher redshift. This (unsurprisingly) extends our initial finding that gas web is an unbiased tracer of total matter out to high redshift, especially for the filament environment. 3. 3.

The gas density PDFs indicate that knots host the densest gas and the mean gas density keeps dropping from filaments to voids. This is independent of the redshift evolution of the total gas density PDFs. 4. 4.

By separating the gas mainly in temperature, we find about 40 per cent of gas is WHIM. The WHIM mass fraction drops with redshift, especially from $z=1$ (29 per cent) to $z=2.1$ (10 per cent). 5. 5.

The detailed gas phase diagram shows the WHIM occupies the largest amount of gas in filaments. By separating the whole WHIM mass into the four LSE, about a half of the WHIM is in filaments and this fraction seems unchanged over redshift. 6. 6.

We do not find significant difference between the two hydrodynamical runs for these fractions, besides the hot gas fraction in the knot environment at high redshift.

We have shown that cold gas and WHIM share a similar fraction of gas mass at $z=0$. However, most cold gas is located in the sheet environment while half of the WHIM mass is in the filament (see also Martizzi et al., 2018). Therefore, we predict that the filament environment is the most promising place to look for WHIM and solving the missing baryon problem. However, because the whole WHIM mass fraction drops with increasing redshift, only $\sim 10$ per cent of total baryonic mass at $z=2.1$, the signal coming from the WHIM (mostly in the soft X-ray band with the emission lines from highly ionised oxygen atoms such as OVII and OVIII, see Fang et al., 2005; Bertone et al., 2008; Bertone et al., 2010; van de Voort & Schaye, 2013, for example), should be much stronger at lower redshift. The largest fraction of gas at high redshift is in cold phase, which should be detected by SKA.

Our findings are generally in agreement with Martizzi et al. (2018)(see also Nevalainen et al. (2015) for a similar result on the WHIM fractions in filaments from hydro-simulation and Eckert et al. (2015) for the results from observations), who applied the tidal tensor method to the IllustrisTNG simulations. For example, the filament volume fractions are around 20 per cent with mass fractions of 40 per cent; the total hot gas mass fractions are 7.5% versus 4.6% at $z=0$, 1.6% versus 1.3% at $z=1.0$ and 0.24% versus 0.2% at $z\sim 2$; the total WIHM mass fractions are 49.9% versus 41.3%, 35% versus 28.7% at $z=1.0$ and 18.5% versus 10.8% at $z\sim 2$; the WHIM mass fractions to the total gas mass in filaments, sheets and knots are 0.27, 0.17, 0.03 in their study, whereas we obtain, respectively, 0.21, 0.12 and 0.06.

Although our simulation resolution is not very high, we believe that our results are almost unaffected by that. This is because our analysis mainly focuses on larger scales. Finer resolution may provide more correct figures at smaller scale. However, the overall baryon distributions, for example baryonic fractions within our mesh cell, are more model dependent than resolution dependent. As van Daalen et al. (2011) (see their figure A2) have shown, the effect of resolution on the power spectrum is within 1 per cent at $k\approx 5~{}h{\rm/Mpc}$, which is the scale corresponding to our cell size.

Acknowledgements

We thank the referee for their thorough and thoughtful review of our paper. The authors would like to thank Giuseppe Murante and Stefano Borgani for preparing these simulations, which have been carried out at the CINECA supercomputing Centre in Bologna, with CPU time assigned through ISCRA proposals and with the support from the PRIN-INAF12 grant ’The Universe in a Box: Multi-scale Simulations of Cosmic Structures’, the PRINMIUR 01278X4FL grant ’Evolution of Cosmic Baryons’, the INDARK INFN grant and ’Consorzio per la Fisica di Trieste’.

This work has made extensive use of the Python packages — Ipython with its Jupyter notebook (Pérez & Granger, 2007), NumPy (van der Walt et al., 2011) and SciPy (Oliphant, 2007; Millman & Aivazis, 2011). All the figures in this paper are plotted using the python matplotlib package (Hunter, 2007). This research has made use of NASA’s Astrophysics Data System and the arXiv preprint server.

WC, AK, GY and RM are supported by the Ministerio de Economía y Competitividad and the Fondo Europeo de Desarrollo Regional (MINECO/FEDER, UE) in Spain through grant AYA2015-63810-P. WC further acknowledges the supported by the European Research Council under grant number 670193. AK further acknowledges funding through the Spanish Red Consolider MultiDark FPA2017-90566-REDC and thanks Blumfeld for l’etat et moi. NIL acknowledges financial support of the Project IDEXLYON at the University of Lyon under the Investments for the Future Program (ANR-16-IDEX-0005). XY is supported by the National Science Foundation of China (NSFC, grant Nos. 11833005, 11621303). WC wishes to acknowledge support by the National Key R&D Program of China (grant No. 2018YFA0404502) and by the NSFC (grant No. 11821303). SP is “Juan de la Cierva” fellow (ref. IJCI-2015-26656) of the Spanish Ministerio de Economía y Competitividad (MINECO) and acknowledges additional support by the MINECO through the grant AYA2016-77237-C3-3-P. HW acknowledges the support from the NSFC under grant No. 11733004 and 11421303.