Simulated metal and HI absorption lines at the conclusion of Reionization

TL;DR

This study uses high-resolution simulations to analyze intergalactic metal absorption lines during and after the Epoch of Reionization, matching several observational data and highlighting discrepancies in ionic relationships at high redshift.

Contribution

It provides a detailed theoretical framework for modeling metal and HI absorption lines at high redshift, incorporating post-processing with CLOUDY and Voigt profile fitting, and compares results with observations.

Findings

Simulations reproduce the evolution of CII and CIV mass densities.

CII shows bimodal distribution with large absorptions near galaxies.

Good agreement with observed HI column density distribution at z=4-6.

Abstract

We present a theoretical study of intergalactic metal absorption lines imprinted in the spectra of distant quasars during and after the Epoch of Reionization (EoR). We use high resolution hydrodynamical simulations at high redshift ($4 <z<8$), assuming a uniform UV background Haardt--Madau 12, post-processing with CLOUDY photoionization models and Voigt profile fitting to accurately calculate column densities of the ions CII, CIV, SiII, SiIV and OI in the intergalactic medium (IGM). In addition, we generate mock observations of neutral Hydrogen (HI) at $z<6$. Our simulations successfully reproduce the evolution of the cosmological mass density ($\Omega$) of CII and CIV, with $\Omega_{CII}$ exceeding $\Omega_{CIV}$ at $z >6$, consistent with the current picture of the tail of the EoR. The simulated CII exhibits a bimodal distribution with large absorptions in and around galaxies, and some traces in the lower density IGM. We find some discrepancies between the observed and simulated column density relationships among different ionic species at $z=6$, probably due to uncertainties in the assumed UV background. Finally, our simulations are in good agreement with observations of the HI column density distribution function at $z = 4$ and the HI cosmological mass density $\Omega_{HI}$ at $4 < z < 6$.

Tables

Table 1. Table 1: Summary of the simulations used in this work. Column 1: run name. Column 2: box size. Column 3: Plummer-equivalent comoving gravitational softening length. Columns 4 and 5: mass of gas and dark matter particles. All the simulations have the same initial number of gas and DM particles ( 2 × 512 3 2 superscript 512 3 2\times 512^{3} ). Column 6: feedback model. Column 7: inclusion of low-temperature metal and molecular cooling (Maio et al., 2007 ; Maio & Tescari, 2015 ) . The first run, Ch 18 512 MDW, is the fiducial model. The second one in the list, Ch 18 512 MDW mol, has exactly the same configuration as the reference run, but includes low-T metal and molecular cooling.

Simulation	Box size	Comoving softening	$M_{\text{gas}}$	$M_{\text{DM}}$	Model for	low-T metal &
	(cMpc/$h$)	(ckpc/$h$)	($\times$ 10${}_{}{}^{5}M_{☉}^{}$/$h$)	($\times$ 10${}_{}{}^{6}M_{☉}^{}$/$h$)	SN–driven winds	molecular cooling
Ch 18 512 MDW	18	1.5	5.86	3.12	Momentum–driven
Ch 18 512 MDW mol	18	1.5	5.86	3.12	Momentum–driven	✓
Ch 18 512 EDW	18	1.5	5.86	3.12	Energy–driven
Ch 18 512 EDW mol	18	1.5	5.86	3.12	Energy–driven	✓
Ch 12 512 MDW mol	12	1.0	1.74	0.925	Momentum–driven	✓
Ch 25 512 MDW mol	25	2.0	15.73	8.48	Momentum–driven	✓

Table 2. Table 2: Contribution from DLA systems to Ω H  i subscript Ω H  i \Omega_{\text{H\,{i}}} .

Simulation	$\alpha =\frac{{\Omega }_{\text{H }\text{i}}}{{\Omega }_{\text{DLA}}}$	%
Ch 18 512 MDW	1.19 $\pm$ 0.03	80.9
Ch 18 512 MDW mol	1.25 $\pm$ 0.02	75.3
Ch 18 512 EDW	1.15 $\pm$ 0.03	84.9
Ch 18 512 EDW mol	1.22 $\pm$ 0.02	77.8
Ch 12 512 MDW mol	1.29 $\pm$ 0.04	70.4
Ch 25 512 MDW mol	1.21 $\pm$ 0.02	78.7

Equations

$$\zeta (m)=\{\begin{array}{ll}\text{0.497}\times m^{-\text{0.2}}& \text{0.1}{M}_{\text{sun}}\le m<\text{0.3}{M}_{\text{sun}},\\ \text{0.241}\times m^{-\text{0.8}}& \text{0.3}{M}_{\text{sun}}\le m<1M_{\text{sun}},\\ \text{0.241}\times m^{-\text{1.3}}& m\ge 1M_{\text{sun}}.\end{array}$$

$${\dot{M}}_w=\eta {\dot{M}}_{\star }.$$

$$v_w=2\sqrt{\frac{G{M}_h}{R_{200}}}=2\times v_{\text{circ}},$$

$$R_{200}=\sqrt[3]{\frac{3}{4\pi }\frac{M_h}{200{\rho }_c{\Omega }_{0m}}}(1+z{)}^{-1},$$

$$\eta =2\times \frac{{\text{600 km s}}^{-1}}{v_w}.$$

$$\eta =2\times {\left(\frac{{\text{600 km s}}^{-1}}{v_w}\right)}^2,$$

$$\frac{{\Gamma }_{\text{phot}}}{{\Gamma }_{\text{UVB}}}=(1-f){\left[1+{\left(\frac{n_{\text{H}}}{n_0}\right)}^{\beta }\right]}^{{\alpha }_1}+f{\left[1+\frac{n_{\text{H}}}{n_0}\right]}^{{\alpha }_2},$$

$${\rho }_X(j)=a^3\sum _{i}X(i)W_{ij},$$

$$(\rho T{)}_X(j)=a^3\sum _{i}X(i)T(i)W_{ij},$$

$$(\rho v{)}_X(j)=a^3\sum _{i}X(i)(a\dot{x}(i)+\dot{a}\left[x(i)-x(j)\right])W_{ij},$$

$$q_{ij}=\frac{a|x(i)-x(j)|}{h_{ij}},$$

$$h_{ij}=\frac{1}{2}\left[h(i)+h(j)\right],$$

$$\tau (u)=\frac{{\sigma }_{0,I}c}{H(z)}{\int }_{-\infty }^{\infty }{n}_I(x)V\left[u-x-v_{\text{pec,||}}^{\text{IGM}}(x),b(x)\right]dx,$$

$$f(\text{N},X)=\frac{n_{\text{sys}}(\text{N},\text{N}+\Delta \text{N})}{n_{\text{lov}}\Delta X},$$

$$\Delta X=\frac{H_0}{H(z)}(1+z{)}^2\Delta z.$$

$$\Delta z=(1+z)\frac{\Delta v}{c},$$

$${\Omega }_{\text{ion}}(z)=\frac{H_0{m}_{\text{ion}}}{c{\rho }_{\text{crit}}}\frac{\sum \text{N}(\text{ion},z)}{n_{\text{lov}}\Delta X},$$

$$f_{\text{Hi}}(\text{N},X)=\frac{n_{\text{sys}}(\text{N},\text{N}+\Delta \text{N})}{\Delta X\cdot n_{\text{cube}}\cdot n_{\text{grid}}^2},$$

$${\Omega }_{\text{Hi}}(z)=\frac{H_0{m}_{\text{Hi}}}{c{\rho }_{\text{crit}}}{\int }_{{\text{N}}_{\text{min}}}^{{\text{N}}_{\text{max}}}{f}_{\text{Hi}}(\text{N},z)\text{N}d\text{N},$$

Full text

1 Introduction

Understanding the epoch of Reionization (EoR) is one of the current challenges of extragalactic astronomy. It will complete the big picture of the thermal history of the Universe (e.g. McQuinn, 2015). The EoR commenced when the first stars switched on, producing UV radiation that transformed the neutral Hydrogen (H i) in the surrounding circumgalactic and intergalactic media (CGM and IGM) into ionized hydrogen (H ii). Observations of Lyman–$\alpha$ photons (Ly$\alpha$ 1216 Å) in absorption towards $z_{\text{em}}\gtrsim$ 6 quasars (QSOs) show fluctuations in flux consistent with Reionization concluding at $z\lesssim$ 6 (e.g. Fan, 2006; Becker et al., 2015b). Measuring Lyman series absorption or UV emissivity at wavelengths blueward of Ly$\alpha$ becomes almost impossible at redshifts greater than 5.5 due to the increasing density of matter and neutral Hydrogen fraction (Becker et al., 2015a).

Metal absorption lines are an alternative proxy to the Lyman series for probing the evolution of the IGM during the EoR and offer many advantages at high redshift. For example, they can be detected even when Hydrogen is completely saturated in the spectra of the background quasars, since some of these transitions occur redward of Ly$\alpha$ and are thus unaffected by Lyman series absorption and Gunn–Peterson troughs (Gunn & Peterson, 1965). Also, ionic transitions in the QSO spectra at high redshift provide important constraints on the model of the ionizing background that drove Reionization (e.g. Furlanetto & Mesinger, 2009; Becker et al., 2011, 2015b; Finlator et al., 2016). Low ionization transitions (O i, C ii, Si ii, Mg ii, Fe ii) trace the location of neutral Hydrogen at high redshift. On the other hand, high ionization states (C iv, Si iv, O vi) have a comparably larger ionization potential than H, therefore the energy required to produce these transitions is reached in regions where Hydrogen is highly ionized at early times.

The direct measurement of absorption features imprinted on the spectra of quasars at high redshift is the best method to infer the ionic ratios, but the observations depend on the wavelength bands where each transition can be found and the observational sensitivity is limited by the decreasing likelihood of an absorption line to be detected. In order to complement the observational techniques, hydrodynamical simulations are used to simulate the physical environment at high redshift where the absorption occurs. The numerical approach can produce a large number of sightlines, improving the statistical estimation of the column densities at redshifts that are not accessible due to current observational limits. Although the field is rapidly expanding, this approach is computationally expensive and the large dynamical range of the underlying physical phenomena makes a true self–consistent simulation impossible. Current simulations are still not able to simultaneously resolve the small scales (turbulence, shocks, fluctuations in the background density, etc.) and large scales (cosmological structures such as clusters and filaments) involved in the progression of the EoR.

An extensive observational effort has been pursued with C iv, the triply ionized state of Carbon. C iv offers many advantages that make its detection easier than other ionic species: large oscillator strength, wavelength redward of Ly$\alpha$ emission and a doublet transition. Observations of C iv in the foreground of quasar spectra at high redshift by Songaila (2001, 2005), Pettini et al. (2003), Ryan-Weber et al. (2006, 2009), Simcoe (2006), Simcoe et al. (2011), D’Odorico et al. (2010); D’Odorico et al. (2013, 2016) and Boksenberg & Sargent (2015) have built the largest sample of metal absorbers as a function of redshift to date. The cosmological mass density $\Omega_{\text{C\,{iv}}}=\frac{\rho_{\text{C\,{iv}}}}{\rho_{\text{crit}}}$ (the density of C iv with respect to the critical density today) drops at high redshifts during the progression of the EoR, due to the changing ionization state of the IGM and a decrease of its metallicity.

Numerical simulations by Oppenheimer & Davé (2006), Oppenheimer et al. (2009), Tescari et al. (2011), Cen & Chisari (2011), Pallottini et al. (2014), Finlator et al. (2015), Rahmati et al. (2016) and a comparison by Keating et al. (2016) have tried to reproduce the evolution of $\Omega_{\text{C\,{iv}}}$, taking into account different feedback prescriptions, photoionization modeling and variations in the UV ionizing background at high redshift. These theoretical efforts provide insight into the physical environments where C iv absorptions occur and show that C iv traces the distribution of the IGM at high temperature up to a few hundred kpc away from galaxies (e.g. Oppenheimer et al., 2009).

On the other hand, low ionization state ions, such as neutral oxygen O i, should trace the distribution of H i at high redshift. In fact, since the ionization potential of O i differs from the neutral Hydrogen one by 0.02 eV, the two ions sit in tight charge-exchange. The detections of O i toward high redshift QSOs in Becker et al. (2006, 2011) are consistent with large variations from one line of view to another, i.e. with an inhomogeneous distribution of the ionizing sources as well as the absorbers. Numerical results from Finlator et al. (2013) showed a tight correlation among H i and O i column densities and studies by Keating et al. (2014) at N${}_{\text{H\,{i}}}>$1017cm*-2* (where the gas is self-shielded) revealed an excellent agreement between O i and H i fractions, regardless of the photoionization model assumed, as well as an increasing incidence rate of O i at higher redshift, consistent with the IGM being more neutral when approaching the EoR (Becker et al., 2011).

Likewise, C ii has been observed by Becker et al. (2006) and modelled by D’Odorico et al. (2013) at $z=$ 5.7 to be best fitted by low density gas, $\delta=\frac{\rho_{\text{gas}}}{\langle\rho\rangle}-1=10$, where ${\langle\rho\rangle}$ is the mean density at the considered redshift. Low ionization metals as C ii and Mg ii are routinely detected at low redshifts in regions close to the centre of galaxies or in damped Ly$\alpha$ systems (DLAs, i.e. absorption systems with N${}_{\text{H\,{i}}}>10^{20.3}$ cm*-2*, Wolfe et al., 1986; Wolfe et al., 2005). These latter are optically thick structures that constitute the main reservoir of neutral Hydrogen after the EoR. DLAs have been studied both observationally (e.g. Péroux et al., 2003; Prochaska et al., 2005; O’Meara et al., 2007; Prochaska & Wolfe, 2009; Crighton et al., 2015) and theoretically/numerically (e.g. Nagamine et al., 2004; Pontzen et al., 2008; Barnes & Haehnelt, 2009; Tescari et al., 2009; Bird et al., 2014; Rahmati et al., 2015; Maio & Tescari, 2015) in order to understand their physical properties and statistical distribution. However, many open questions still remain on the nature of these systems and their connection with the chemical enrichment of the Universe.

The aim of this work is to i) compare high-$z$ observations of metal and H i absorbers with synthetic spectra and ii) use simulations to gain insight into the underlying physical explanation for the observations. This paper is divided as follows: in section 2, we introduce our numerical simulations. In section 3, we explain the method used to include the UV background Haardt & Madau 2012 (hereafter HM12 UVB), the photoionization modeling and the column density calculation procedure. In section 4, we show the evolution with redshift of the physical properties that describe the IGM. Sections 5 and 6 analyze the ionic species in the simulations focussing on the cosmological mass density and the relationships among different ions. In the subsequent section 7, we discuss the evolution of the cosmological mass density of neutral Hydrogen at $z<6$. Finally, in sections 8 and 9 we summarize our results and the implications of the models. Throughout the paper we use the prefix c for comoving and p for physical distances.

2 Cosmological simulations

The numerical simulations used in this work reproduce representative volumes of the Universe at redshift 4 $<z<$ 8 and were run with the smoothed particle hydrodynamics (SPH) code P-GADGET3(XXL) $-$ a customized version of GADGET-3 (Springel, 2005). The model is an extension at high redshift of the AustraliaN GADGET-3 early Universe Simulations (Angus) project (Tescari et al., 2014; Maio & Tescari, 2015; Katsianis et al., 2015, 2016, 2017). Among other technical improvements, the suite of hydrodynamical simulations takes into account: a multiphase star formation criterion from Springel & Hernquist (2003), self-consistent stellar evolution and chemical enrichment modeling (Tornatore et al., 2007), supernova (SN) momentum– and energy–driven galactic winds111Although active galactic nuclei (AGN) feedback is implemented in the code (Springel et al., 2005; Fabjan et al., 2010; Planelles et al., 2013), in this work we do not consider it. We stress that for the analysis presented in the paper the role of AGN feedback is expected to be negligible (see e.g. Keating et al. 2016 and references therein). (Springel & Hernquist, 2003; Puchwein & Springel, 2013), metal–line cooling (Wiersma et al., 2009) and low-temperature cooling by molecules/metals (Maio et al., 2007). In addition, the code is supported by a parallel Friends-of-Friends (FoF) algorithm to identify collapsed structures and SUBFIND to classify substructures within FoF haloes.

Each simulation generates a cosmological box (including periodic boundary conditions) with the same initial number of gas and dark matter (DM) particles for a total of $2\times 512^{3}$. Initial masses of the gas and DM particles are given in Table 1. Whenever the gas density is above a threshold $\rho_{\text{th}}$, there is a probability that a gas particle will turn into a star222This probability is calibrated to reproduce the Kennicut-Schmidt law (Schmidt, 1959).. Stochastically, a new star–type particle is introduced in the simulation. Each star particle represents a simple stellar population with mass 0.1$M_{\sun}\leq m\leq$ 100 $M_{\sun}$. The stars with mass $m\leq$ 40 $M_{\sun}$ explode as SNe before turning into a black hole, while stars that have masses above this threshold collapse into a black hole without passing through the SN stage.

A flat $\Lambda$CDM model with cosmological parameters from the latest release of the Planck Collaboration et al. (2016) is assumed: $\Omega_{0m}=$ 0.307, $\Omega_{0b}=$ 0.049, $\Omega_{\Lambda}=$ 0.693, $n_{s}=$ 0.967, $H_{0}=$ 67.74 km s*-1Mpc-1* (or $h=$ 0.6774) and $\sigma_{8}=$ 0.816. The simulations were calibrated according to the parameters used in Tescari et al. (2014) and Katsianis et al. (2015), and are compatible with observations of the cosmic star formation rate (SFR) density history and the galaxy stellar mass function at $z=$ 6 to 8.

The subgrid scheme takes into account the lifetimes of stars of different mass and follows the evolution of Hydrogen, Helium, molecules and metals (C, Ca, O, N, Ne, Mg, S, Si and Fe) released from SNIa, SNII and low and intermediate mass stars. It is possible to vary the initial mass function (IMF), the lifetime function and stellar yields. The algorithm is ideally suited to modeling IGM enrichment produced by galactic winds blown by “starburst” galaxies at high redshift. Radiative cooling and heating processes are included according to Wiersma et al. (2009) and Maio et al. (2007). The IMF considered is the Chabrier multi–sloped (Chabrier, 2003), as described in equation (1), that produces a large number of intermediate– and high–mass stars, expected to play an important role during Reionization.

[TABLE]

We adopted the following stellar yields:

• •

SNIa: Thielemann et al. (2003). The mass range for the SNIa originating from binary systems is 0.8 $M_{\rm\odot}<m\leq 8$ $M_{\rm\odot}$, with a binary fraction of 7%.

• •

SNII (massive stars): Woosley & Weaver (1995). The mass range for SNII is $m>8$ $M_{\rm\odot}$.

• •

Asymptotic giant branch (low and intermediate mass) stars: van den Hoek & Groenewegen (1997).

2.1 Feedback mechanisms

In order to regulate the star formation and chemically enrich the IGM, the simulations have been set up to account for kinetic supernova–driven winds. According to the nature of the feedback and its effectiveness, galactic winds are classified as momentum– or energy– driven winds.

2.1.1 Momentum–driven galactic winds (MDW)

Momentum–driven galactic winds (Puchwein & Springel, 2013) are the reference feedback model used in this work. We assume that the mass-loss rate associated with the winds $\dot{M}_{w}$ is proportional to the star formation rate $\dot{M}_{\star}$ through $\eta$, the wind mass loading factor that accounts for the efficiency of the wind:

[TABLE]

In this model, the velocity of the winds $v_{w}$ is regulated by the mass of the host halo and $R_{200}$, the radius where the density is 200 times larger than the critical density at redshift $z$:

[TABLE]

where $v_{\text{circ}}$ is the circular velocity and $R_{200}$ is defined as follows:

[TABLE]

with $\Omega_{0m}$ and $\rho_{c}$ the matter density and the critical density today. The conservation of momentum of the winds imposes that $\eta$ is proportional to the inverse of the wind–velocity $v_{w}$:

[TABLE]

The normalization factor was chosen to reproduce the observed cosmic SFR density and galaxy stellar mass function up to $z=8$. According to the model, weak winds ($v_{w}<$ 600 km s*-1*) lead to a large efficiency and more material is expelled from supernovae with respect to strong winds ($v_{w}>$ 600 km s*-1*) that reach further but load less material333Please note that in these simulations the normalization of the wind mass loading factor, $\eta_{0}=600$ km s*-1*, is slightly larger than in previous Angus runs, such as those presented in Tescari et al. (2014) and Katsianis et al. (2015), where $\eta_{0}=450$ km s*-1*. This is due to the fact that these simulations do not include AGN feedback (and in particular the early AGN feedback in low-mass galaxies model), and therefore a recalibration of the wind model was necessary for properly matching the cosmic star formation rate density history and galaxy stellar mass function up to $z=8$.. Stochastically, some gas particles are selected to be part of the wind, and subsequently, decouple from the hydrodynamical scheme for a given amount of time $t_{\text{dec}}$ (Tescari et al., 2014).

2.1.2 Energy–driven galactic winds (EDW)

To test the impact of the adopted feedback model on our results, we also use energy–driven winds (Springel & Hernquist, 2003). The only remarkable difference with momentum–driven winds is the scaling form of $\eta$. As for MDW, the speed of the winds is regulated by $v_{w}=2\times v_{\text{circ}}$, but in this case due to conservation of the energy of the winds, the wind mass–loading factor $\eta$ scales as the square of the inverse of $v_{w}$,

[TABLE]

making EDW more aggressive than MDW, especially in low-mass galaxies. Thermal feedback produced by SNIa and SNII is also considered, in addition to the kinetic feedback just described.

3 Methodology

The set of simulations was run to $z=$ 4, with initial conditions at $z=$ 125, using the Raijin supercluster from the National Computational Infrastructure (NCI) facility444http://nci.org.au. These simulations describe the physical conditions of the gas and the evolution of H, He, C, Ca, O, N, Ne, Mg, S, Si and Fe released from SNIa and SNII. We post-process the simulations introducing a uniform UV ionizing background: a field radiation due to the CMB and the Haardt & Madau (2012) ultraviolet/X-ray background from quasars and galaxies with saw-tooth attenuation that evolves with redshift (hereafter HM12).

Assuming this UVB, we compute the ionization states of each element with Cloudy photoionization code v8.1 (Ferland et al., 2013) for optically thin gas in ionization equilibrium, focusing in particular on the following ions: H i, C ii, C iv, Si ii, Si iv and O i. We consider only the transition with the highest oscillator strength and rest-frame wavelength $\lambda_{\text{rest}}\geq$ 1216 Å, which can be observed at high redshift. The assumed metallicity in the Cloudy tables is solar.

In addition, an effective prescription for H i self-shielding was introduced to accurately describe the regions at the centre of the galaxies, where N${}_{\text{H\,{i}}}$ is significantly higher than in the IGM due to the shielded bubbles that contain pristine Hydrogen. We adopt the parametric function of Rahmati et al. (2013):

[TABLE]

which best–fitting parameter values reproducing the radiative transfer results at $1<z<5$ are: $\alpha_{1}=-2.28\pm 0.31$, $\alpha_{2}=-0.84\pm 0.11$, $\beta=1.64\pm 0.19$, $n_{0}=(1.003\pm 0.005)$ $n_{\text{H,SSh}}$ (where $n_{\text{H,SSh}}$ is the self-shielding density threshold) and $f=0.02\pm 0.0089$. The equilibrium neutral Hydrogen density is obtained using a recombination rate given by Hui & Gnedin (1997) and the photoionization rate as a function of the temperature of Theuns et al. (1998).

In order to mimic real observations and avoid the introduction of bias in the data, we generate random lines of sight inside the simulated box along the three perpendicular directions and extract the relevant physical information taking into account positions, velocities, densities and temperatures of the SPH particles inside each line of sight. Subsequently, we compute a synthetic spectra in density and optical depth/flux as a function of the velocity width (1024 pixels) for each ion in the simulation, according to the procedure introduced by Theuns et al. (1998). For a bin $j$ at a position $x(j)$, the density and density weighted temperature and velocity are calculated from:

[TABLE]

[TABLE]

[TABLE]

where $a$ is the scale factor, $X(i)$ the abundance of the species $X$ of the SPH particle $i$ and $W_{ij}=mW(q_{ij})/h_{ij}^{3}$ the normalized SPH kernel. $W$ is the SPH kernel, $m$ the particle mass and:

[TABLE]

[TABLE]

with $h$ the physical softening scale.

With this information, it is possible to calculate the number density of the ion transition considered $n_{\text{ion}}$. The synthetic flux for any transition at the redshift-space coordinate $u$ is given by $F(u)=\text{exp}\left[-\tau(u)\right]$, with $\tau(u)$:

[TABLE]

where $\sigma_{0,I}$ is the cross-section of the ion transition, $H(z)$ the Hubble parameter at $z$, $x$ the space coordinate in km s*-1*, $b$ the velocity dispersion (in units of $c$) and $V$ the Voigt profile. The spectra can be converted from the velocity space $v$ to the observed wavelength using $\lambda=\lambda_{0}(1+z)(1+v/c)$.

In numerical works, it is a common practice to normalize H i fluxes averaged over a large number of random lines of sight to the observed mean normalized flux of the Ly$\alpha$ forest at a given redshift, $\langle F(z)\rangle=\exp(-\tau_{\text{eff}})$, through a constant rescaling factor $A_{\text{H\,{i}}}$. For consistency, also the fluxes of metal ions are rescaled by the same constant factor. However, at $z>5$ reliable measurements of $\tau_{\text{eff}}$, and therefore $A_{\text{H\,{i}}}$, are not possible, due to the thickening of the Ly$\alpha$ forest. To overcome this problem, we adopted a reversed approach. First, we calculated the C iv column density distribution function at $z=$ 4.8 and 5.6 (see section 5.1) using three different simulations (Ch 18 512 MDW, Ch 18 512 MDW mol and Ch 18 512 EDW) and a range of scaling factors for the C iv optical depths, $\tau_{\text{C\,{iv}}}$. Then, we selected the scaling factor, $A_{\text{C\,{iv}}}$, which provided the best chi-by-eye agreement between the simulated and observed distribution functions. Finally, we rescaled all the ionic fluxes (including H i) in all the simulations and at all the redshifts considered (i.e. also at $z\leq 5$, for consistency) by $A_{\text{C\,{iv}}}$. The best fit value is $A_{\text{C\,{iv}}}=0.85$, very close to unity.

The individual spectra are convolved with Gaussian noise profiles with full width at half maximum $\text{FWHM}=$ 7 km s*-1* to produce final synthetic spectra with a signal-to-noise ratio, $S/N=$ 50, comparable with observations obtained using the UVES spectrograph mounted at the Unit 2 of the Very Large Telescope. Finally, to obtain a fair comparison with real data, we fit all the individual absorption features in the spectra through Voigt profile components in a range of 100 km s*-1* among contiguous systems with the code Vpfit v.10.2 (Carswell & Webb, 2014).

The column density N, equivalent width $EW$ and Doppler parameter $b$ from each spectrum is used to produce a sample for each ion. The estimated errors of the column density $\Delta$N, Doppler parameter $\Delta$b and redshift $\Delta$z of the absorption features are also calculated with Vpfit. The individual fits are selected if the conditions N $>\Delta$N and b $>\Delta$b are fulfilled simultaneously. Otherwise, the component is rejected.

As a final remark, we note that radiative transfer effects are not included in our simulations. The models use the evolving HM12 uniform UVB to quantify the ionization state of the CGM/IGM at a given redshift. The aim of this work is not to follow the progression of Reionization, nor the evolution of the H ii bubbles or their topology. The implicit assumption is that our boxes (that are small compared to the size of the H ii bubbles at the redshifts of interest) represent a region of the Universe already reionized at a level given by the HM12 UVB. At $6<z<8$, chemical enrichment occurs mostly inside and in close proximity of galaxies (interstellar medium, CGM and high density IGM) where, assuming an inside-out progression of Reionization, the gas in which metals lie should be ionized. Although proper RT calculations would be more accurate, they are extremely expensive from the computational point of view. Using a uniform (evolving) UVB is a common approach in the literature (see e.g. Oppenheimer et al., 2009; Keating et al., 2016). Moreover, Finlator et al. (2015) studied the evolution of Carbon absorption from $z=10$ to $z=5$ with cosmological hydrodynamic simulations that include a self-consistent multifrequency, inhomogeneous UVB. They found that the difference between their more realistic UVB and the uniform HM12 is within $\sim 2-4$ times, which according to them is fair agreement given the uncertainties.

We stress that in this work we do not present results on the neutral Hydrogen fraction at $z>6$. In this case, a uniform background slowly ionizing the IGM is a completely wrong picture. For this reason, we only study the H i fraction at $z<6$. As already mentioned before, to accurately describe the regions at the centre of galaxies we also include an effective prescription for H i self-shielding (Rahmati et al., 2013).

4 Evolution of the IGM during Reionization

Figure 1 shows the distribution of the gas particles in the fiducial simulation Ch 18 512 MDW in diagrams of overdensity–temperature at redshifts $z=$ 8, 6 and 4, comparing total metallicity, number of particles and H i fraction (left, middle and right panels, respectively). These redshifts are chosen to study the IGM at i) $z=$ 8 when Reionization is progressing; ii) $z=$ 6, when the overlap of the H ii bubbles is nearly complete at the end of the EoR; and iii) $z=$ 4 when the Reionization of Hydrogen has fully concluded and the only remaining reservoirs of neutral Hydrogen are self-shielded regions of H i and DLA systems.

The diagrams are split in well–defined regions whose sizes evolve with time in terms of the overdensity, $\Delta=\delta+1=\frac{\rho_{\text{gas}}}{\langle\rho\rangle}$, and the temperature $T$:

• •

voids: gas particles at very low densities ($\Delta<1$);

• •

IGM: the gas in this regime follows an adiabatic relation $T=T_{o}\Delta^{\gamma-1}$ and is in ionization equilibrium: cooling is counterbalanced by photoheating. The conditions of the particles in the IGM are $T\leq$ 104 K and $1\leqslant\Delta<10$;

• •

Circumgalactic/interstellar media (CGM/ISM) represent the transition between the IGM and the inside of galaxies, with typical densities of $\Delta\sim$ 101-3. This environment is heated by photoionization fronts coming from galaxies. The sharp feature at $\Delta\sim$ 103 represents the density threshold where gas particles turn stochastically into star–particles in the simulation;

• •

star–forming region: particles in this regime are above $\Delta\sim$ 103, have temperatures higher than 104 K and follow an effective equation of state imposed by the subgrid star formation model (Springel & Hernquist, 2003);

• •

shock–heated gas: a growing region of the IGM that is heated by feedback processes at late times ($T>$ 105 K and low densities).

At $z=$ 6 – 8, most of the gas in the diagrams is located in the so–called diffuse phase $-$ low $T$ ($\leq$ 104 K) and low $\Delta$ ($\leqslant$ 102) $-$ the first stars are being formed and the main processes that raise the temperature of the IGM (and voids) are taking place.

The leftmost panels show the metallicity with respect to solar. The enriched sections of the phase diagram are mainly star-forming regions and some gas at very high temperature, expelled from galaxies through supernova–driven winds. On the other hand, voids exhibit zero metallicity, because the galactic outflows do not reach regions that far from the centre of the galaxies. The CGM and ISM are chemically enriched at high redshift. Instead, the IGM is very metal-poor at $z=$ 8. At late stages (bottom panels), feedback prescriptions spread out material from the high-density gas to regions originally empty at higher redshifts. Chemical enrichment contributes progressively more to the cooling processes.

The right panels show the distribution of Hydrogen in the reference simulation. Yellow indicates Hydrogen completely neutral ($X_{\text{HI}}\geqslant$ 10*-1*), dark green particles containing H ii ($X_{\text{HI}}<$ 10*-7*) and the intermediate colors are used to illustrate the transition between the two regimes. It is clear that most of the reservoirs of H i are located at low temperature and intermediate– to high–density regions, that are presumably where DLAs or isolated self-shielded regions are. On the other hand, the shock–heated gas with temperatures $T>$ 105 K contains a large amount of ionized gas. The fraction of this gas that has been chemically enriched is expected to be found in high ionization states, such as C iv or Si iv, consistent with Cen & Chisari (2011).

Please note that, for the reasons discussed at the end of section 3, the top right panel should not be regarded as an accurate representation of the neutral Hydrogen content in the Universe at $z=8$, since a uniform UVB does not properly describe the progression of Reionization. We included this panel just for illustrative purposes.

5 Evolution of Carbon

In order to study the evolution of the IGM at high redshift, we use synthetic spectra and analyze the evolution of metal absorption lines.

5.1 CIV column density distribution function

C iv, the triply ionized state of Carbon, is the transition most detected and studied in the foreground of quasar spectra at high redshift. We devote this section to its cosmological evolution and compare our theoretical predictions with the available observational archive.

The column density distribution function (CDDF) is an observable that takes into account the statistical distribution of absorption systems with respect to their column densities. Once the absorptions are detected and the redshift path is confirmed to be complete, the CDDF can be used to build the cosmological mass density of a particular ion. By construction, the CDDF, or $f(\text{N},X)$, quantifies the number of absorption systems $n_{\text{sys}}$ in the column density interval $(\text{N},\text{N}+\Delta\text{N})$ in an absorption path $\Delta X$:

[TABLE]

where $n_{\text{lov}}$ is the number of lines of view (lov) considered. The absorption path relates the Hubble parameter at a given redshift $z$ with the correspondent redshift path $\Delta z$ as follows:

[TABLE]

The equivalent redshift path of our cosmological boxes is defined as:

[TABLE]

where $\Delta v$ is the box size in km s*-1* at a given redshift $z$.

We impose a cut–off in column density, $\log$ N${}_{\text{th}}$ (cm*-2*) $=$ 12.5, to the synthetic data to mimic the sensitivity of available observations and to avoid including poorly sampled features. Although a large number of systems below this column density threshold is present ($\sim 80$%), these have not been taken into account in the statistics. Moreover, we limit the calculation of the C iv–CDDF to N${}_{\text{\text{C\,{iv}}}}$ (cm*-2*) $<$ 1015.2, consider bins of 0.4 dex and introduce Poissonian errors for the theoretical sample as $\sqrt{n_{\text{sys}}}$ to fairly compare with the observations of D’Odorico et al. (2013).

Figure 2 shows the predicted column density distribution function of C iv at $z=$ 4.8 and compares with observations by D’Odorico et al. (2013) and two fitting functions proposed by the authors: $f($N$)=B$N*-α* with $B=$ 10.29 $\pm$ 1.72 and $\alpha=$1.75 $\pm$ 0.13 and $f($N$)=f($N${}_{0})($N$/$N${}_{0})^{-\alpha}$ with $f($N${}_{0})=$ 13.56 and $\alpha=$ 1.62 $\pm$ 0.2. Although we used this statistics as a rough guide to calibrate the C iv optical depths in three of our simulations (see the final part of section 3), there is good agreement among all the simulations and the observational data at this redshift, and small deviations from the observations are within the error bars. The error bars in C iv–CDDF shown correspond to the reference run and are representative for all the simulations.

The same observable is computed at $z=$ 5.6 in Figure 3. The values from the numerical simulations are shown in points and they are compared to systems detected by D’Odorico et al. (2013). The observational fitting functions parameters for $f($N$)=B$N*-α* are $B=$ 8.96 $\pm$ 3.31, $\alpha=$ 1.69 $\pm$ 0.24 and for $f($N$)=f($N${}_{0})($N$/$N${}_{0})^{-\alpha}$ are $f($N${}_{0})=$ 14.02 and $\alpha=$ 1.44 $\pm$ 0.3. There is good agreement between the observational and theoretical points except at large column densities, where simulations deviate to values higher than observations. We point out that the function fits at $z=$ 4.8 and 5.6 extend to high column densities, even though observations are only available up to N${}_{\text{\text{C\,{iv}}}}=$ 1014.4 cm*-2*. Therefore, the outcome of the simulations gives additional information in a range of N${}_{\text{C\,{iv}}}$ where the C iv–CDDF has not been measured.

We stress that, at these redshifts, the C iv–CDDF cannot be used to disentangle different physical prescriptions in the simulations, since all the runs are in agreement within the error bars.

Finally, it is important to note that the run Ch 12 512 MDW mol predicts the highest C iv–CDDF at low column density both at $z=$ 4.8 and 5.6. The higher resolution of this run with respect to all the others leads to a more complete and precise estimation of the column density distribution function at low N${}_{\text{\text{C\,{iv}}}}$.

In similar works by Finlator et al. (2015, 2016) and Keating et al. (2016), strong absorbers are highly disfavored for different simulated UV backgrounds and feedback prescriptions. In Finlator et al. (2015, 2016), the feedback mechanisms are quite efficient at enriching the IGM with metals. However, their simulated boxes are relatively small. Therefore, galaxies with large halo masses are suppressed. If these structures host the high column density systems, such simulations will also lack strong absorptions. On the other hand, Keating et al. (2016) explore various feedback prescriptions. Their Sherwood run has a similar hydrodynamics and configuration as our models, whereas HVEL implements a more aggressive version of the energy–driven winds used in this work (including also a minimum wind speed of 600 km s*-1*). Nevertheless, all their models struggle to reproduce the high C iv column density absorbers. This issue is not alleviated by varying the UVB. The presence of rare strong C iv absorbers in our theoretical models results from an appropriate level of resolution on the scale of the absorbers and a post-processing pipeline that closely mimics the method used by observers and accounts for individual features to calculate the column densities from Voigt profile fits to the absorption lines. The latter increases the accuracy in the estimation of the column densities.

5.2 CIV cosmological mass density

The evolution of the total density in C iv ions with respect to the critical density is described by the comoving mass density $\Omega_{\text{C\,{iv}}}$. For any ion, $\Omega_{\text{ion}}$ can be derived by summing the column densities of identified absorbers:

[TABLE]

where $m_{\text{ion}}$ is the mass of the ionic species, $n_{\text{lov}}$ is the number of lines of view (lov), $\rho_{\text{crit}}$ is the critical density today and $\Delta X$ is defined in equation (15).

In the case of C iv, we restricted the column densities in the sum to the range 13.8 $<$ $\log$ N${}_{\text{C\,{iv}}}$(cm*-2*) $<$ 15.0, to compare directly with the results of D’Odorico et al. (2013). The results for all of the simulations are shown in Figure 4. The theoretical results are compared to observations by Pettini et al. (2003) and Ryan-Weber et al. (2009), Songaila (2001, 2005), Simcoe et al. (2011), D’Odorico et al. (2013), Boksenberg & Sargent (2015) and Díaz et al. (in prep). It is worth mentioning that the observations by Pettini, Ryan-Weber and Díaz have been recalibrated to the cosmology adopted in this paper. In the other cases, missing details of the precise pathlength probed do not allow us to convert the observations to the Planck cosmology.

The rise of C iv over the redshift period from 8 to 4 reflects both the increase in the chemical enrichment of the intergalactic medium through galactic feedback mechanisms, and the evolution of the ionization state of the IGM. At $z=$ 8, C iv is largely suppressed and most of the Carbon in the box is in its neutral state. As time passes, more SN events take place and pollute the IGM with metals (as we will show in section 5.4). Simultaneously, the specific intensity of the HM12 UVB around the wavelength where the C iv transition occurs ($\lambda_{\text{C\,{iv}}}=192$ Å) increases by $\sim 4$ orders of magnitude between $z=$ 8 and 4. At all redshifts, each simulation is broadly consistent with current observational data.

Using results from a hybrid wind model, Finlator et al. (2015, 2016) integrated a power law fit to the column density distribution function to calculate the evolution of the comoving mass density of C iv. Instead, we calculate $\Omega_{\text{C\,{iv}}}$ by summing the column density of the absorbers using equation (17). This method is more accurate when there are uncertainties in the slope of the C iv–CDDF, especially at high column densities.

5.3 CII cosmological mass density

We extend the analysis to C ii and compare with the observational lower limits of Becker et al. (2006) in Figure 5, taking into account systems whose column densities are in the range 13.0 $<$ $\log$ N${}_{\text{C\,{ii}}}$(cm*-2*) $<$ 15.0. In the context of Reionization, it is interesting to consider the evolution of C ii (a low ionization state of Carbon with ionization energy of 11.26 eV) with respect to C iv, whose ionization energy is more than 4 times larger (47.89 eV). Since the UVB is softer at earlier times, C ii is expected to be more dominant than higher ionization states of Carbon at high redshifts.

On the left and central panels of Figure 5 we draw a comparison between different physical scenarios among the simulations with boxsize 18 cMpc/$h$: the further left relates simulations with the same feedback prescription (MDW) with and without low temperature metal and molecular cooling. The run with low-T cooling included produces slightly more C iv (and less C ii) than the fiducial run, but the difference is only significant at $z>$ 6.5. In the middle panel, we contrast simulations without low temperature metal and molecular cooling with different feedback prescriptions: MDW and EDW. The different physical scenarios do not have a strong impact on the evolution of C iv, as seen also in Figure 4. The trends for C ii are more dependent on the model, because this low–ionization state is very sensitive to the densities of the absorbers (which are correlated to the strength of the feedback prescription), and the Ch 18 512 EDW run is in good agreement with the observations of Becker et al. (2006).

Although the run with lowest resolution, Ch 25 512 MDW mol, produces slightly degraded results, the evolution of $\Omega_{\text{C\,{iv}}}$ does not change significantly with the spatial resolution or size of the simulations (see the rightmost panel of Figure 5).

On the other hand, C ii trends are more influenced by the resolution. Ch 12 512 MDW mol better reproduces the observed evolution of this ion. The absorption features are affected by the box size (and the comoving softening) for different reasons: the absorption path calculated to create the synthetic spectra grows with the box size, as well as the size of the pixels, while the definition of the individual features decreases for larger configurations. The Ch 12 512 MDW mol run is in better agreement with the observational limits for C ii at high redshift, indicating that higher resolution in the simulations helps to better describe the absorption features, in particular for low ionization states.

The most interesting part of the analysis arises when the C ii and C iv curves are compared for each model: in the fiducial simulation, the amount of C iv exceeds the amount of C ii at $z\sim$ 6, consistent with the fact that the Universe is approaching the tail of the EoR. This result is seen in all of the simulations, with a crossover in the range of $z\sim$ 6–6.5.

Despite the fact that in some runs the synthetic C ii does not reach the lower limits predicted in Becker et al. (2006), the crossover of these ions at $z\sim$ 6–7 is very promising. Other theoretical works that also follow the trend of these ions (e.g. Finlator et al., 2015) obtained a crossover between $\Omega_{\text{C\,{ii}}}$ and $\Omega_{\text{C\,{iv}}}$ at $z\sim$ 8. However, a later crossover of C ii and C iv at $z\sim$ 6 is more consistent with the observations and the current paradigm of the tail of Reionization.

5.4 Total Carbon cosmological mass density

To conclude section 5, in Figure 6 we analyse the evolution of the cosmological mass density of (total) Carbon in the simulated boxes, $\Omega_{\text{C}}$. This is calculated differently than $\Omega_{\text{C\,{iv}}}$ and $\Omega_{\text{C\,{ii}}}$: here we sum the amount of Carbon in each gas particle in the simulation and divide by the comoving volume. Thus, $\Omega_{\text{C}}$ gives an estimation of the total amount of Carbon at high redshift. The trend is consistent with the enrichment history of the Universe, increasing by 1.5 orders of magnitude from $z=$ 8 to 4. At $z=4$, our results are compatible with those of Tescari et al. (2011). Given the different methods used to calculate $\Omega_{\text{C}}$ and $\Omega_{\text{C\,{iv}}}$ & $\Omega_{\text{C\,{ii}}}$, the first quantity should not be directly compared with the other two. Moreover, the evolution of the total Carbon cosmological mass density is a purely theoretical prediction of our simulations, since $\Omega_{\text{C}}$ is not (directly or indirectly) observable.

The overall evolution of C is nearly independent on the nature of the cooling (blue vs red and green vs magenta lines), but depends at low redshift on the feedback mechanism. The EDW models quench the formation of stars in galaxies more effectively than MDW, therefore less Carbon is produced and $\Omega_{\text{C}}$ is lower.

Figure 6 also compares simulations with different spatial resolutions (the comoving softening of the different configurations scales with the size of the box). The Ch 12 512 MDW mol, Ch 18 512 MDW mol and Ch 25 512 MDW mol runs have comoving softening of 1.0, 1.5 and 2.0 ckpc/$h$ (yellow, red and purple lines), respectively. These simulations predict a similar trend for the mass density of Carbon and at low redshift converge to the same $\Omega_{\text{C}}$. However, at high redshift the resolution affects the amount of metals. The run with the highest resolution (Ch 12 512 MDW mol, yellow line) can resolve higher densities at earlier times and therefore better describes star formation in galaxies and produces more Carbon than Ch 18 512 MDW mol (red) and, particularly, Ch 25 512 MDW mol (purple).

6 Column density relationships of metal absorption lines

One of the advantages of the numerical approach is that it reproduces a large number of sightlines, improving the statistical estimation of the column densities at redshifts that are hard to access with current observations. By convolving our synthetic spectra with Gaussian noise and extracting column densities using Voigt profile fitting, our simulated results closely mimic the observed spectra.

As mentioned in section 5.1, this procedure produces a very large number of systems with low column densities, that display non physical linear correlations among the metal absorption features. To avoid this issue, a cut–off in column density, $\log$ N${}_{\text{th}}$ (cm*-2*) $\geq$ 12.5, has been imposed to the synthetic data, and systems below this noise threshold have not been taken into account in the statistics.

We explore relationships in column densities among different ionic species in a uniform UV background HM12 at $z=$ 6 using random lines of sight to emulate the observational method and compare with the observations available to date in Figure 7. In panel A, N${}_{\text{C\,{ii}}}$ vs. N${}_{\text{C\,{iv}}}$ are compared, panel B displays N${}_{\text{Si\,{iv}}}$ vs. N${}_{\text{C\,{iv}}}$ and panel C the ratio of N${}_{\text{C\,{ii}}}$/N${}_{\text{C\,{iv}}}$ vs. N${}_{\text{Si\,{iv}}}$/N${}_{\text{C\,{iv}}}$. All these quantities are compared to systems detected by D’Odorico et al. (2013), while panels D, E and F draw a comparison with low ionization states observed by Becker et al. (2006). The comparisons N${}_{\text{C\,{ii}}}$ vs. N${}_{\text{Si\,{ii}}}$, N${}_{\text{C\,{ii}}}$ vs. N${}_{\text{O\,{i}}}$ and N${}_{\text{Si\,{ii}}}$ vs. N${}_{\text{O\,{i}}}$ are shown in panels D, E and F, respectively. In all cases, we display the simulations with 18 cMpc/$h$ box size, to analyze possible variations among the predicted values with different prescriptions of feedback and cooling. However, the dispersion of the data points do not allow clear distinctions between the models.

The high ionization states are relatively well represented by the simulated data with the UV ionizing background implemented, especially N${}_{\text{C\,{iv}}}$ as seen in panels A and C. On the other hand, the low ionization species are produced far in excess by the simulations and in all the different configurations considered the column densities are overestimated for O i and Si ii.

C ii is a very interesting case in this study. The numerical estimation of C ii column densities do not seem to be well reproduced by the simulations, although in most of the cases the observational data are just a lower limit, and a bimodality in the distribution of column densities appears, which will be analyzed in section 6.1.

The scatter plots in many cases struggle to reproduce the observed column densities, regardless of the feedback model used, possibly due to the UV background implemented. The parameters that characterize HM12 can suppress the low ionization states. Bolton & Viel (2011) claimed that N${}_{\text{Si\,{iv}}}/$N${}_{\text{C\,{iv}}}$ is sensitive to the spectral shape of the UV background, whereas Finlator et al. (2015) showed that N${}_{\text{C\,{ii}}}/$N${}_{\text{C\,{iv}}}$ is sensitive to the overall intensity of the UV background and its normalization. We are currently studying the effect of variations in the normalization and hardness of the HM12 UVB. The results will be presented in a follow-up paper.

O i receives special attention because its ionization energy is similar to the one of H i. Besides the overprediction of N${}_{\text{O\,{i}}}$ in our simulations555There is also a constant trend of N${}_{\text{O\,{i}}}$ at low column densities, possibly indicating that in low-density environments Oxygen is mostly found in higher ionization states due to the effect of the UVB., the large column densities in this ionic species may indicate that there is a huge amount of O i self-shielded that it is not being detected due to the low probability of tracing a high–density region with the current observational methods. It is worth noting that, to date, there is not a self-consistent scheme in the simulation literature that accounts for the self-shielding of the ions (corrections for self-shielding at $z<$ 5 can be found in Bird et al. 2014 and Bolton et al. 2017). This could make the estimated column densities of the low ionization states less accurate.

6.1 Bimodality in CII column densities

The location of the low-ionization states of metals at $z\sim$ 3, for instance C ii and O i, traces the regions close to the centre of the galaxies or DLA systems. This correlation of gas ionization state and proximity to the centre of the nearest galaxy may no longer be valid at high redshift, as C ii may be tracing either DLAs or the IGM at low temperature.

Observations at $z\sim$ 6 exhibit N${}_{\text{C\,{ii}}}\leq$ 1014 cm*-2*, while at higher redshift only lower limits in C ii are available. In all the simulations, we found a remarkable bimodal distribution in the N${}_{\text{C\,{ii}}}$ calculated using random sightlines through each box (see e.g. panel A in Figure 7). In order to test the nature of this bimodality, in Figure 8 we plot the column densities of C ii and C iv, using lines of sight at different impact parameters from the centre of the galaxies in the fiducial run Ch 18 512 MDW ($d=$ 20, 100 and 500 ckpc/$h$ in yellow, red and green points, respectively). At each $d$, we extracted a thousand lov around halos with masses in the range 10${}^{9-10}M_{\sun}$/$h$. The blue points represent the column densities along 1000 random lines of sight through the box shown in panel A of Figure 7. The black stars are observational data from D’Odorico et al. (2013). The top panel in Figure 8 displays the distribution of the C ii absorption features and the right–most window the corresponding distribution of C iv column densities above the noise threshold.

When the impact parameter $d$ approaches the galaxies, the peak at high column densities becomes sharper, indicating that there are more C ii absorbers in the CGM where the temperature and density are larger than in the IGM. For $d=$ 20 ckpc/$h$, $\sim$ 3 pkpc at this redshift, the absorbers are located inside or in the outskirts of the galaxies (spectroscopic confirmation of galaxies at $z=$ 6–8 with the Hubble Space Telescope shows that the radii of these galaxies are in the range of 0.6–1.1 pkpc, Jiang et al. 2013). Absorbers are gravitationally bounded to the halos in the simulated sample, which have virial radii of 3–30 pkpc.

The distribution of the C ii column densities at $d=$ 20 ckpc/$h$ indicates that the few observations towards high-redshift QSOs currently available have not yet detected rare systems with large column densities in C ii, due to the very low likelihood of intersecting a galaxy at a low impact parameter (that would produce a deep absorption feature). Furthermore, D’Odorico et al. (2013) have only identified C ii absorptions that were first detected as C iv doublets. Thus, any low ionization gas that may have produced a strong C ii absorber and no associated C iv has not been reported. However, future observations of QSOs will raise the number of metal absorption lines detected and cover a larger range in column densities.

The results for N${}_{\text{C\,{ii}}}$ with random lines of sight (blue line) agree with numerical results by Oppenheimer et al. (2009), Keating et al. (2016) and a photoionization comparison by D’Odorico et al. (2013), predicting that C ii lies mostly in overdense regions with $\delta=\frac{\rho_{\text{gas}}}{\langle\rho\rangle}$ - 1 $>$ 10, and there are some traces in low temperature regions (mostly the IGM probed by current observations) where its column densities are of the same order of N${}_{\text{C\,{iv}}}$.

Finally, numerical results in the CGM (yellow and red cases in Figure 8) lead to a bimodality in C iv and a clear excess of high column density systems at low impact parameters. Oppenheimer et al. (2009) showed that the properties of the C iv absorbers depend mostly on the evolution of the environment of their hosting galaxy (distance, mass and metallicity). By tracing the surroundings of the galaxies in halos with masses 10${}^{9-10}M_{\sun}$/$h$, we confirm that there is a correlation between C iv column densities and the proximity of the galaxies that host the absorbers.

7 Cosmological mass density of HI at redshift lower than 6

The estimation of the amount of neutral Hydrogen can be derived using DLAs. It is worth noting that at high redshift ($z>$ 5) the nature of DLAs is ambiguous (e.g. Simcoe et al., 2012), since the definition of a DLA is by column density only (N${}_{\text{H\,{i}}}>$ 1020.3 cm*-2*) and does not discriminate between dense pockets of neutral Hydrogen in the IGM, where Reionization is still progressing, and DLAs associated with self-shielded gas in galaxies (like at $z<$ 5).

However, it is possible to simulate DLA systems to calculate their corresponding comoving mass density $\Omega_{\text{DLA}}$ and from it, estimate $\Omega_{\text{H\,{i}}}$. We identify the halos with masses above 10${}^{9}M_{\sun}/h$ (using a halo finder routine, the halos are sorted by mass) and extract cubes of size 200 ckpc/h around the center of mass of one thousands of them (randomly distributed in mass). In this way, we guarantee that the simulated DLAs are not affected by a selection bias in density. Then, the neutral Hydrogen density in the cube is integrated along the line of sight and projected into a 2D grid of N${}_{\text{H\,{i}}}$. With these column densities, we finally calculate the H i–CDDF, $f_{\text{\text{H\,{i}}}}(\text{N},X)$, as:

[TABLE]

where $n_{\text{sys}}$ is the number of systems in the interval $(\text{N},\text{N}+\Delta\text{N})$, $\Delta X$ the absorption path defined in equation (15), $n_{\text{cube}}$ the number of cubes generated and $n_{\text{grid}}$ the number of bins in the grid. The parameters selected to achieve a good resolution of the systems that contain neutral Hydrogen are: $n_{\text{cube}}=$ 1000 and $n_{\text{grid}}=$ 64 (i.e. the grid size is 3.1 ckpc/h).

Figure 9 shows the H i–CDDF, with column densities in the range 12 $<$ $\log$ N${}_{\text{H\,{i}}}$(cm*-2*) $<$ 22, for all the simulations at $z=$ 4 and compares with the observations of Prochaska et al. (2005), Crighton et al. (2015), Bird et al. (2017) (all in the column densities corresponding to DLA systems at $z\geq$ 4) and O’Meara et al. (2007) in the super Lyman limit systems (LLS) range at $z\sim$ 3.4666Although these systems are at $z<4$, they are a reliable benchmark for our simulated data, since there is evidence of low to no evolution in this redshift range.. Note that there is a very large number of systems with low column densities in all the models and DLAs are extremely rare. The agreement with the observational data is good, indicating that the method implemented to find neutral systems works well at these redshifts.

Simulations with smaller comoving softening (or smaller box size), better resolve the Ly$\alpha$ forest. In Figure 9 the yellow line (simulation Ch 12 512 MDW mol) is above the red line, corresponding to the run Ch 18 512 MDW mol, and the latter is above the purple line (Ch 25 512 MDW mol, the run with the lowest resolution in the suite of simulations) at $\log$ N${}_{\text{H\,{i}}}$(cm*-2*) $<$ 19.

At DLA column densities, the distributions follow a double–power law, as discussed by Zwaan & Prochaska (2006) and Prochaska & Wolfe (2009), with a knee at $\log$ N${}_{\text{H\,{i}}}$ (cm*-2*) $\sim$ 21 in all of the cases considered. Also, at N${}_{\text{H\,{i}}}>$ 1021 cm*-2*, where molecular cooling plays an important role in the star forming regions, the H i–CDDF in models with molecular cooling implemented drops off more rapidly than in the other models.

In Figure 10 a zoom of the DLA region is shown ($\log$ N${}_{\text{H\,{i}}}$(cm*-2*) in the range $[$20.3 , 22.0$]$). Overplotted are the most recent observational data of H i absorbers at high redshift by Prochaska et al. (2005), Crighton et al. (2015), Bird et al. (2017) and the fitting function proposed by Prochaska & Wolfe (2009) for DLA systems at redshift 4.0–5.5. The simulations are compatible with the observational data, especially runs with no molecular cooling.

The main distinction between the models with and without molecular cooling is driven by the conversion of neutral Hydrogen to H2, that becomes important in high density regions where new stars are being formed. Molecules boost the cooling of the surrounding gas and the formation of stars. As a result, the amount of atomic Hydrogen decreases.

At this point, we calculate the comoving mass density of H i, using the following equation:

[TABLE]

with $m_{\text{H\,{i}}}$ the mass of the Hydrogen atom and $f_{\text{H\,{i}}}(\text{N},z)$ defined in equation (18). In theory, the integral should be computed with an upper limit $\text{N}_{\text{max}}=\infty$, but in reality this value is set by the maximum column density detected. Thus, the limits considered for $\Omega_{\text{H\,{i}}}$ are 12 $<$ $\log$ N${}_{\text{H\,{i}}}$ (cm*-2*) $<$ 22.

We compare the predictions from the theoretical models to observations at high redshift by Prochaska et al. (2005), Prochaska & Wolfe (2009), Zafar et al. (2013) and Crighton et al. (2015) in the left panel of Figure 11.

The amount of H i at $z=$ 4 in the Ch 18 512 MDW mol (red line) and Ch 25 512 MDW mol (purple line) runs is in excellent agreement with the observational data, probably because the inclusion of molecular cooling prevents the production of too much neutral Hydrogen (by converting it into molecular Hydrogen). On the other hand, simulations without molecular cooling overproduce $\Omega_{\text{H\,{i}}}$, especially in the EDW case (green line). All the models predict an increasing amount of H i at $z>$ 5.5. However, we stress again that our assumption of a uniform UVB is no longer valid to describe the diffuse H i in the IGM when approaching $z=$ 6. At these redshifts, our simulated $\Omega_{\text{H\,{i}}}$ only describes the amount of neutral Hydrogen in collapsed/self-shielded systems. Please note that for this calculation we also extended the H i self-shielding prescription of Rahmati et al. (2013), which was originally proposed only up to redshift $z=$ 5.

In order to estimate the contribution of DLAs to the overall H i budget, the right panel of Figure 11 compares the prediction of $\Omega_{\text{H\,{i}}}$ for all the simulations (solid lines) with the DLA comoving mass density $\Omega_{\text{DLA}}$ (dashed lines), obtained by integrating the H i–CDDF in the range 20.3 $<$ $\log$ N${}_{\text{H\,{i}}}$(cm*-2*) $<$ 22. The theoretical estimations are compared with recent data from SDSS (Bird et al., 2017) in black stars. In this case, simulations without molecular cooling included are in better agreement with observations (as for the H i–CDDF shown in Figure 10).

We investigate the assumption $\Omega_{\text{H\,{i}}}=\alpha\times\Omega_{\text{DLA}}$, where $\alpha$ is a factor that accounts for the numbers of systems that contain neutral Hydrogen with N${}_{\text{H\,{i}}}\geq$ 1020.3 cm*-2*. The values for $\alpha$ in each simulation (averaged in the redshift range 4–6) are in Table 2. In the observational work of Crighton et al. (2015), the authors estimate a 20 per cent contribution from systems with column density below the DLA threshold (or 80% contribution from DLAs) at $z<5$. This assumption seems to be validated with the numerical results here reported.

8 Discussion and summary

The purpose of this work is to investigate the physical environment of the IGM in proximity of galaxies at the end of the Epoch of Reionization by producing mock spectra and measuring ionic column densities of metal transitions with a method that closely resembles current observational techniques. Taking advantage of a suite of simulations with a very robust model for chemical enrichment, SN–driven feedback that reproduces the cosmic SFR density and the galaxy stellar mass function at high redshift and a suitable comoving softening at the IGM scale (1–2 ckpc/$h$), we reproduce well some observables using metal and H i absorption lines at high redshift. In particular, the comoving mass density $\Omega$ of C iv and C ii in the redshift range of 4 $\leq z\leq$ 8, the H i column density distribution function at $z=$ 4 and $\Omega_{\text{H\,{i}}}$ & $\Omega_{\text{DLA}}$ at $z<6$. In all cases, we compare the theoretical results with the latest release of observational data.

8.1 Effect of wind feedback

The simulations we used in this work are based on the analysis of Tescari et al. (2014), where there is an extensive discussion of the feedback mechanisms implemented and how they successfully reproduce the observations of the SFR at high redshift (e.g. Madau & Dickinson, 2014). When these models are used to study the evolution of different ionic species, we find that it is hard to disentangle the effect of different SN-driven galactic wind prescriptions for some of the statistics considered, in particular the C iv and H i column density distribution functions. The reasons are mostly numerical: probably higher resolution is needed to better resolve patches of enriched gas at relatively low-density in the CGM/IGM (as suggested by our resolution tests). Moreover, the SPH algorithm is not very effective in the mixing of metals. Therefore, only extreme feedback models would produce significant differences in the observables at high redshift. As a result, when calculating the column densities of metal absorption lines with Voigt profile fitting at $z>$ 5, it is not possible to distinguish between momentum–driven and energy–driven winds (MDW and EDW). The dispersion in the relationships among the metal absorption lines (see Figure 7) is too high to provide definitive hints on the evolution of the ionization states due to one particular wind model.

Nonetheless, a clear difference between the numerical estimation of the total Carbon content, $\Omega_{\text{C}}$, in EDW and MDW models arises at $z\leqslant$ 6 (Figure 6). This is due to the fact that EDW quenches the overall star formation in galaxies more effectively than MDW, and consequently the production of Carbon in the simulation is 0.2–0.3 dex lower at $z=$ 4.

8.2 Influence of the low temperature metal and molecular cooling

As expected, the implementation of low temperature metal and molecular cooling affects mostly neutral Hydrogen statistics. In particular, molecular cooling suppresses H i in regions where the majority of it is self-shielded, mainly because of the conversion of H i to H2 at high densities. This effect is important when computing the H i column density distribution function and cosmological mass density. Interestingly, the introduction of molecular cooling improves the estimation of $\Omega_{\text{H\,{i}}}$ at $z\sim$ 4–5 (left panel of Figure 11), whereas, when the H i–CDDF is computed in the DLA regime, 20.3 $<$ $\log$ N${}_{\text{H\,{i}}}$ (cm*-2*) $<$ 22, runs without molecular cooling provide a better match to new SDSS observations (Bird et al., 2017) of $\Omega_{\text{DLA}}$ (right panel of Figure 11) and the distribution function at $z=$ 4 (Figure 10).

Once the feeback prescription is fixed, simulations with low-T metal and molecular cooling included produce slightly more C iv than the other runs (Figure 4 and left panel of Figure 5). The effect is less important for C ii and not visible for the total Carbon content and the dispersion in the relationships among the metal absorption lines (Figures 6 and 7).

8.3 Limitations

Even though our models are in good agreement with observations of ionic species at high redshift, there are some caveats to consider.

First of all, radiative transfer effects are not included. The simulations are run in a relatively small volume (the maximum size considered is 25 cMpc/$h$), whereas a typical H ii bubble at the tail of EoR is larger than 100 cMpc/$h$. Therefore, we assume that our boxes represent a region of the Universe already reionized by a uniform HM12 UV ionizing background at $z\leq 8$. Although full RT calculations would be more accurate, this is a fair approximation to describe the ionization state of metals (Finlator et al., 2015), since chemical enrichment at these redshifts occurs mostly inside and in close proximity of galaxies. On the other hand, even in small volumes the distribution of H i optical depths cannot be properly described with a uniform UVB or a simple model that assumes galaxies and quasars as ionizing sources and uses a fixed mean free path for the ionizing photons (Becker et al., 2015b). For this reason, we have studied neutral Hydrogen statistics only at $z<6$.

A necessary step to complement the analysis of this work will be to vary the spectral hardness and the normalization of the UV background, in order to refine the calculation of the column densities of some ionic species, that are not well represented with the assumed uniform HM12 UVB. We plan to investigate this in a subsequent paper.

In recent years, numerical modeling of absorption features have improved enormously, but the current generation of simulations still struggles to reproduce low ionization states of metals at high redshift. This is mainly due to insufficient resolution, but it may also be influenced by the lack of a consistent implementation of self-shielding prescriptions for the ions. This work (like many others) just considers the effect of H i self-shielding (Rahmati et al., 2013), but does not introduce any self-shielding of low ionization absorbers (which lay in clumpy structures), that has been proposed for DLA systems at lower redshift (Bird et al., 2015).

Finally, there is a fundamental problem with the lack of reliable measurements of the mean normalized flux of the Lyman–$\alpha$ forest at $z\geq$ 5, when the QSOs spectra display very large Gunn–Peterson troughs, which leads to the absence of an effective parametrization of the optical depth of H i to re-scale or calibrate the simulated Hydrogen and ion spectra. Different approaches have been used in the literature to reproduce the evolution of the metal ions independently from H i. In this work, we calibrate our simulations using the C iv column density distribution function at $z=$ 4.8 and 5.6.

9 Conclusions

We have introduced high-resolution hydrodynamical simulations to study the physical and chemical state of the IGM at the end of the Epoch of Reionization. We have tested the effect of a low temperature metal and molecular cooling model and two prescriptions for galactic wind feedback on the evolution of metal absorption lines in the redshift range 4 $<z<$ 8. We have also studied the H i cosmological mass density at $z<6$. Our theoretical predictions are consistent with the available observations at high redshift.

The evolution of the high and low ionization states of Carbon is the main focus of our work. The drop in the cosmological mass density of C iv, $\Omega_{\text{C\,{iv}}}$, from $z=$ 4 to 8 is due to the combined effect of a change in the ionization state of the gas and the decreased metallicity of the IGM (Figure 4). In fact, the total Carbon comoving mass density at $z=$ 8 is more than a factor of 10 lower than at $z=$ 4 (Figure 6), and the simulations show a consistent transition to more neutral states of metals at high redshift (as shown in Figure 5). Most notably, the crossover between $\Omega_{\text{C\,{iv}}}$ and $\Omega_{\text{C\,{ii}}}$ happens at $z\sim$ 6–6.5, in agreement with available observations and the current paradigm of the tail of Reionization.

In comparison to other numerical works in the literature, we are able to produce C iv absorbers in the IGM with large column densities. This is mostly due to a combination of efficient galactic winds and a post-processing pipeline that mimics observational methods and, in particular, accounts for individual features to calculate the column densities from Voigt profile fits to the absorption lines.

The simulated C ii exhibits a bimodal distribution with large absorptions in and around galaxies, and some traces in the lower density IGM. These latter correspond to systems detected with current observations (D’Odorico et al., 2013; Becker et al., 2006). We predict that, at high redshift, most of the high column density C ii, with $\log$ N${}_{\text{C\,{ii}}}$ (cm*-2*) $>$ 15, should lie in proximity of galaxies at impact parameters of order 20 ckpc/$h$ (Figure 8) and has not been detected yet due to the low likelihood of reaching these overdense regions with lines of sight towards distant QSOs.

We have studied the column density relationships among different ionic species (Figure 7). High ionization states (like N${}_{\text{C\,{iv}}}$) are reasonably well described by our simulations with the adopted HM12 UV background. On the other hand, a comparison with observations in the literature shows that the low ionization states are not well represented, regardless of the feeback model implemented.

Finally, our simulations are in good agreement with observations of the H i column density distribution function at $z=$ 4 (Figures 9 and 10) and the H i cosmological mass density, $\Omega_{\text{H\,{i}}}$, at 4 $<z<$ 6 (Figure 11). We validate the estimate made in the observational work of Crighton et al. (2015) that DLA systems contribute to $\sim 80$% of $\Omega_{\text{H\,{i}}}$ (see Table 2).

In a companion paper (García et al., 2017), we have explored the likelihood of reproducing the observed Ly$\alpha$ emitter galaxy – C iv absorption pair detected by Díaz et al. (2015) and studied the physical processes that produced the metal enrichment in the IGM at $z\geqslant$ 5.6. Future work will be focused on analyzing more metals and ions to have a broader perspective on the evolution of the IGM at the end of the EoR. Some of the column density relationships are extremely sensitive to the ionizing background, therefore varying the normalization and hardness of the UVB is the next step to obtain a more realistic description of this cosmic era.

Acknowledgements

Parts of this research were conducted by the Australian Research Council Centre of Excellence for All-sky Astrophysics (CAASTRO), through project number CE110001020. The numerical simulations were run using the Raijin distributed-memory cluster from the National Computational Infrastructure (NCI) facility, post-processed with the Edward HPC machine from the University of Melbourne and analyzed with the gSTAR supercomputer at Swinburne University of Technology. This work was supported by the Flagship Allocation Scheme of the NCI National Facility at the ANU. The authors acknowledge CAASTRO for funding and allocating time for the project Diagnosing Hydrogen Reionization with metal absorption line ratios (fy6) during 2015 and 2016. L.A. García thanks Valentina D’Odorico and Neil Crighton for providing their observational data in private communications and Jeff Cooke for the insightful discussions. E. Ryan–Weber acknowledges ARC DP 1095600. We also thank the anonimous referee for providing valuable comments that considerably improved our manuscript and Klaus Dolag for his help with the technical details of the SPH code.