1.INTRODUCTION

Unveiling the end of the epoch of reionization (EoR) and the sources that complete the budget of ionizing photons is currently a key topic in Astronomy. The phase transition of neutral hydrogen (HI) into its ionized state (HII) is caused by the radiation released by the first stars , the second generation of stars (POP III, ^{Abel et al. 2002}; ^{Bromm et al. 2002}; ^{Yoshida et al. 2003}), and quasars . Other candidates are proposed, such as miniquasars, with masses around 10^3-6M_ʘ^{Dijkstra et al. 2004}; ^{Hassan et al. 2018}). Other candidates are proposed, such as miniquasars, with masses around 10^3-6M_ʘ (^{Mortlock et al. 2011}; ^{Bolton et al. 2011}; ^{Smith et al. 2017}),decaying or self-annihilating dark matter particles or decaying cosmic strings. Nonetheless, the latter objects seem to be unlikely to ionize the Universe by themselves.

Understanding the EoR is intimately tied to the evolution of the ultraviolet background (UVB): the grand sum of all photons that have escaped from quasars and galaxies. Its spectral energy distribution is reasonably well determined at z < 5 (^{Bolton et al. 2005}) and modelled (^{Haardt et al. 2001}) used a cosmological 1D radiative transfer model that follows the propagation of H and He Lyman continuum radiation in a clumpy ionized intergalactic medium (IGM), and uses mean free path and hydrogen photoionization rate decreasing with redshift ¹. However, as z ≥ 6 is approached, the population of UV sources is not well determined (^{Haardt et al. 1999}). The uncertainty in estimating the UV photon emissivity from each type of object is caused by the lack of knowledge on the star formation rate, clumping factor and UV escape fraction (^{Cooke et al. 2014}) at the redshift of interest, which are strongly model-dependent. Measuring Lyman series absorption or UV emissivity in a spectrum blueward of Lyα at 1216 Å is rendered almost impossible by the increasing density of matter and neutral hydrogen fraction at redshifts greater than 5.5.

On the other hand, the assumption of a uniform UV radiation field breaks down close to and during the EoR, when the interaction of the ionizing sources with the IGM requires a very accurate description (^{Lidz et al. 2006}). A real-time reionization simulation should first ionize high density regions and fill some regions before others, leading to a multiphase IGM with spatial fluctuations (^{Lidz et al. 2016}).

Alternatively, intervening metal absorption lines detected in the spectra of high redshift quasars offer a completely different method to calibrate the UVB at high redshift. A growing number of absorption systems has been detected (^{Bosman et al. 2017}; ^{Codoreanu et al. 2018}; ^{Meyer et al. 2019}; ^{Becker et al. 2019}) and with the advent of bigger telescopes (e.g. GMT), the expectation is that the sample of absorption lines detected will significantly increase.

An increase in the sample of absorption lines (of at least an order of magnitude compared to the current observational sample) is possible with numerical simulations. A number of works have directed their efforts to describing the evolution of metal absorption lines in the intergalactic and circumgalactic medium (CGM). These simulations take into account different feedback prescriptions, photoionization modelling and variations in the UV ionizing background in the high redshift Universe (^{Oppenheimer et al. 2006}, ²⁰⁰⁹; ^{Tescari et al. 2011}; ^{Cen et al. 2011}; ^{Finlator et al. 2013}; ^{Pallottini et al. 2014}; ^{Keating et al. 2014}; ^{Finlator et al. 2015}; ^{Rahmati et al. 2016}; ^{Keating et al. 2016}; ^{García et al. 2017}b; ^{Doughty et al. 2018}, ²⁰¹⁹). The methods employed in each of these works, as well as the set-up of the hydrodynamical simulations, show advantages for the description of the IGM.

However, the result of (^{García et al. 2017b}) with a uniform UV background (^{Haardt et al. 2012}) show that the calculated column densities of the low ionization states (CII, SiII, OI) in the simulations have difficulty matching the values observed by ^{Becker et al (2011)}. The lack of spatial resolution on the scale of the low ionization absorbers evidences that further work needs to be done to reach a better description of the environment of these absorbers. Nonetheless, the uncertainties on the assumed UVB at high z suggest that varying its normalization could alleviate the discrepancy between simulated results and the current observations. Works from ^{Finlator et al. (2015}, ²⁰¹⁶⁾; ^{Doughty et al. (2018)} also show that alternative models to the UVB can reduce the gap between the observables associated to metal lines calculated with their simulations and the observations. Their models account for simulated UVB with contribution from galaxies + quasars and quasar-only.

The triply ionized state of silicon (SiIV) offers a unique avenue of investigation. Although it is classified as a high ionization state, its ionization potential energy is significantly lower than CIV; it does not necessarily lie in the same environment and it exhibits the same physical conditions as CIV. Detections at high redshift of this ion have been made by ^{Songaila et al. (2001}, 2005); ^{D'Odorico et al. (2013}), ^{Boksenberg et al. (2015}), and more recently by ^{Codoreanu et al. (2018}). The latter authors identify 7 systems across a redshift path of 16.4 over the redshift range 4.92 < z < 6.13. They are ≈ 50% complete down to a column density of log N_sys (cm^-2) of 12.50. This limiting column density allows them to study the identified SiIV population across the column density range [12.5, 14.0] over the full redshift path of their survey. In addition, Codoreanu et al. (2018) show that the fiducial configuration in ^{García et al. (2017}b) at z = 5.6 is compatible with the SiIV observations.

On the other hand, the self-shielding (SSh) of HI gas in very high density regions (above 10¹⁷ cm^-2) is also a component that needs to be refined in the description, specifically when HI statistics are made. Current studies implement the self-shielding prescription proposed by ^{Chardin et al. (2018)}, but unfortunately, this is only valid up to z = 5. At redshifts when reionization is concluding, this treatment is no longer valid. A new scheme has recently proposed by ^{Rahmati et al. (2013}a), using radiative transfer calculations to find the best fitting parameters from the functional form of the photoionization rate Γphot described in ^{Rahmati et al. (2013a)}. Different SSh prescriptions could have a different outcome in the HI statistics. HI column densities are commonly classified in three regimes: Lyman-α forest (12 < log N^HIcm-2 < 17.2), Lyman limit systems (or LLS with 17.2 < log N_HI cm^-2 < 20.3) and damped Lyman-α absorbers (DLAs) with N_HI > 10^20.33 cm^-2. Works carried out by on DLAs showed that DLAs are the main contributors to the cosmological mass density of HI. The HI self-shielding and molecular cooling prescriptions are important factors as these absorbers reside in low temperature and high density environments.

This paper, in particular, builds on previous findings of ^{García et al. (2017}b). The aim of this work is to explore and discuss to a first approximation variations in the assumed UVB and the HI self-shielding prescriptions.

The paper is presented as follows: § 2 describes the simulations and the method used to post-process them. In § 3 we explore two variations to the uniform assumed in ^{García et al. (2017}b). Additionally, we propose an alternative method to compare the current observations with the synthetic sample of metal absorbers, in contrast with previous works that compare two or more ions at once. § 4 shows results for HI statistics when two different HI self-shielding prescriptions are implemented in post-process. Finally, § 5 summarizes the findings of this paper and explores the limitations encountered.

2. The Numerical Simulations and Post-Process

The results presented in this work are a follow-up to ^{García et al. (2017}b), and are based on the simulations and the methodology presented in that paper. The suite of numerical simulations uses a customized version of GADGET-3 (^{Springel 2005}): P-GADGET3(XXL). The model was first tested in the context of the Angus project. In ^{Tescari et al. (2014)} and ^{Katsianis et al. (2015)}; ^{Katsianis et al. (2016}, ²⁰¹⁷⁾, the authors showed that their simulations were compatible with observations of the cosmic star formation rate density and the galaxy stellar mass function at 1 < z < 7. The model takes into account the following physical processes: a multiphase star formation criterion; self-consistent stellar evolution and chemical enrichment modeling; supernova (SN) momentum- and energy-driven galactic winds; AGN feedback, metal-line cooling; low-temperature cooling by molecules and metals (^{Maio et al. 2007}, ²⁰¹⁵). Moreover, the model is supported by: a parallel Friends-of-Friends FoF) algorithm to identify collapsed structures, and a parallel SUBFIND algorithm to identify substructures within FoF halos.

The numerical model self-consistently follows the evolution of hydrogen, helium and 9 metal elements (C, Ca, O, N, Ne, Mg, S, Si and Fe) released from supernovae (SNIa and SNII) and low and intermediate mass stars. The chemical evolution scheme is based on the stochastic star formation model implemented in the simulations (^{Tornatore et al. 2007}). It accounts for the age of stars of different mass; hence, the amount of metals released over time varies with the mass of the stars.

The lifetime function from for stars with mass m is adopted. The stellar yields quantify the amount of different metals which is released during the evolution of the stellar population, as follows: (i) SNIa: ^{Thielemann et al. (2003)}; (ii) SNII (massive stars): ^{Woosley & Weaver (1995)}; (iii) low and intermediate mass (AGB) stars: van den Hoek & Groenewegen. As one of the main contributors to the reionization of the Universe is the POP III stars (a very massive and short-life population), they are best described by a ^{Chabrier (2003)} initial mass function (IMF).

Galactic scale winds were introduced in simulations by Springel & Hernquist (2003)to regulate the star formation, spread metals from the galaxies and high-density regions to the IGM and shock-heated gas, and prevent the overcooling of gas. The phenomenological model for energy-driven wind feedback is presented in ^{Springel & Hernquist (2003)}. It assumes that the mass-loss rate associated with the winds, M_w, is proportional to the star formation rate M*̇, such that Mẇ=ηM*̇, with η the wind mass loading factor that accounts for the efficiency of the wind. The kinetic energy of the wind is related to the energy input of the supernova. The velocity of the wind is given by the expression vw=2GMhR200=2×vcirc. Due to the conservation of the wind energy, the velocity of the wind goes as the square of the inverse of the loading factor η=2600km/svw2.

However, ^{Puchwein et al. (2013)} suggest that the the mass carried by the wind is not necessarily proportional to the SFR of the galaxy. In such case, it would be more natural to assume that there is a relation between the momentum flux (instead of the energy flux) of the wind and the SFR of the galaxy, thus η is proportional to the inverse of the wind velocity v _w , such that η=2×600km/svw.

Chemical pollution caused by star formation contributes to the cooling of gas. Some metal line cooling efficiencies peak at T ≈ 10⁴ K (mostly low ionization transitions, ^{Gnat & Ferland 2012}). These transitions are privileged in metal-poor high-density environments, as DLAs, where H is mostly neutral or in its molecular form. As discussed in ^{Maio et al. (2007)}, molecular and low temperature metal cooling is particularly important when collapsed structures reach temperatures T < 10⁴K due to the formation of molecules. At this temperature, atomic cooling is not efficient and is highly suppressed, yet, molecular H continues to cool the gas.

The assumed cosmology is a flat ΛCDM model with cosmological parameters from ^{Planck Collaboration (2015)}:Ω_om = 0.307, Ω_ob = 0.049, ΩΛ=0.693 and H₀= 67.74 km s^-1Mpc^-1 (or h = 0.6774). The simulations considered in the paper are described in Table 1, with comoving box size and softening of 18 Mpc/h and 1.5 kpc/h, respectively. All the simulations have the same initial number of gas and DM particles (2 x 512³), with masses of gas and dark matter particles of Mgas=5.86×105M⊙/h and MDM=3.12×106M⊙/h. Moreover, we include molecular and low-temperature metal cooling in our simulations. The fiducial run is labelled Ch 18 512 MDW.

Table 1 Summary of the simulations used in this work 

Simulation	Box size (cMpc/h)	Comoving softening (ckpc/h)	Model for SN-driven winds	Low-T metal % molecular cooling
Ch 18 512 MDW	18	1.5	Momentum-driven
Ch 18 512 MDW mol	18	1.5	Momentum-drive
Ch 18 512 EDW	18	1.5	Energy-driven

Column 1: Run name. Column 2: Box size. Column 3: Plummer-equivalent comoving gravitational softening length. Column 4: Feedback model. Column 5: Inclusion of low-temperature metal and molecular cooling (^{Maio et al. 2007}, ²⁰¹⁵). The first run, Ch 18 512 MDW, is the fiducial model. The second one, Ch 18 512 MDW mol, has exactly the same configuration as the reference run, but includes low-T metal and molecular cooling.

The numerical simulations are post-processed to recreate the observations of high redshift quasars and recover synthetic spectra for each ion along each sightline. The pipeline derived for this purpose relies on the physical conditions of the gas (reproduced with the hydrodynamical simulations), on top of a uniform field radiation that accounts for the cosmic microwave background (CMB) and the ultraviolet/X-ray background from quasars and galaxies HM12 (^{Haard et al. 2012}). The photoionization modeling for the metal transitions is computed with v8.1 (^{Ferland et al. 2013}) for optically thin gas. In addition, a HI self-shielding prescription is imposed to the simulations to account for neutral hydrogen inside high-density regions where gas is optically thick. We choose a thousand random lines of sight along to the three perpendicular directions inside the box, and in each one of these sightlines we calculate a simulated spectrum, containing relevant physical information: HI flux/optical depth, density and temperature of the gas, and the number density of all the ions considered in the analysis, among other quantities. The box size Δv at a given redshift is translated to the equivalent redshift path through the relationship ∆z=(1+z)∆vc. Once the synthetic spectra are computed, they are convolved with Gaussian noise profiles with full width at half maximum FWHM=7kms^-1. Finally, the individual absorption line features are fitted with Voigt profiles with the code v.10.2 (^{Carswell et al. 2014}). We focus our attention on the ionic transitions shown in Table 2.

3.Variations of the UVB in Post-Process

Investigations carried out with a uniform UV background (^{Haard et al. 2012}) in showed that the calculated column densities of the low ionization states (CII, SiII, OI) and their corresponding observables (comoving mass density, column density distribution function, etc.) leave room for improvement when compared with observations. There is general agreement that current simulations do not have enough resolution on the scale of the absorbers with low ionization energies. However, uncertainties in the high-z UVB suggest that varying its normalization is a first step towards a better agreement with the observations at high redshift.

Here we explore the sensitivity of the results presented in to the presence of different ultraviolet/X-ray ionizing backgrounds by modifying the normalization factor at 1 Ryd of the uniform UVB in post-process. Preliminary work allows us to conclude that decreasing the UVB intensity at z = 6 is equivalent to imposing an ionizing background at times earlier than redshift 6. Therefore, low ionization ions would prevail in the early stages of the epoch of reionization. Hence, reductions of more than one order of magnitude in the UVB are very aggressive for high ionization species and, simultaneously, lead to an overproduction of the low ionization ones. Consequently, a variation of 1 dex below the fiducial emissivity in HM12 is conservative but offers a non-negligible imprint on the calculated ions.

Variations of the UVB spectrum (quasars+galaxies model) from the original (see Figure 1) require modified input files to run new tables.

Fig. 1 UV emissivity for the uniform background at three different redshifts: 𝑧= 8, 6 and 4 (blue, dark red and green, respectively) compared with the wavelength of the radiation, in the wavelength range where the ion transitions occur. The color figure can be viewed online. 

The procedure followed here is explained in § 2 and it has been tuned and applied in , with some small adjustments: the normalization parameter at 1 Ryd is reduced by an order of magnitude below the fiducial value in at all redshifts. This leads to a softer UVB. Hereafter, the test is refered as log J_v -1 (see § 3.1.).

3.1.Change in the Normalization of HM12

It is a reasonable expectation that the presence of a softer UVB input than the uniform in the photoionization model would favor low ionization states and more neutral states would show large incidence rates. In order to test this hypothesis, the UV emissivity J_v at 1 Ryd is reduced by one dex compared with the value defined by the uniform at all redshifts. This leads to a softer UVB.

In order to avoid introducing more variables to this test, the box size has been fixed to 18 Mpc/h. The simulations used to recover the observables are Ch 18 512 MDW, Ch 18 512 MDW mol and Ch 18 512 EDW. The convergence and resolution tests are not included in this document, but they can be found in ^{García et al. (2017b)}.

The first comparison with observations explored here is the CIV column density distribution function (CDDF), defined in equation 1 as follows:

f(N,X)=nsys(N,N+ΔN)nlovΔX. (1)

Here, nlov is the number of lines of view considered. The absorption path ΔX=H0H(z)(1+z)2Δz relates the Hubble parameter at a given redshift 𝑧 with the correspondent redshift path Δz=(1+z)Δvc. The term Δv is the box size in km s^-1.

At z = 4.8 and 5.6, the CCDFs are compared with observations from and in Figures 2 and 3, respectively. At z = 6.4, the simulated values of the CDDF are compared with upper limits from in Figure 4.

Fig. 2 CIV column density distribution function at z= 4.8 and comparison with observational data by ^{D'Odorico et al. (2013)}, orange diamonds and ^{Codoreanu et al. (2018)}, yellow squares. The black dashed line represents the fitting function f(N)=BN^-a with B= 10.29 ± 10.29 ± 1.72 and α = 1.75 ± 0.13 and the dotted-dashed line f(N) = f(N₀)(N/N₀)^-a with f(N₀) = 13.56 and α =1.62 ± 1.62 ±0.2, from ^{D'Odorico et al. (2013)}. The error bars are the Poissonian errors for the reference run and are a good representation of the errors in the other models. The plot illustrates results with the test , for the simulations Ch 18 512 MDW, Ch 18 512 MDW mol and Ch 18 512 EDW. In simulations without molecular cooling implemented, the number of CIV absorbers in the log J_v- 1 case is under-represented in the range of column densities considered. The color figure can be viewed online. 

Fig. 3 CIV column density distribution function at z = 5.6 and comparison with observational data by ^{D'Odorico et al. (2013)}, orange diamonds and , yellow squares. The black dashed line represents the fitting function f(N) = BN^-α with B =8.96 ± 3.31 and α =1.69 ± 0.24 and the dotted-dashed line f(N) = f(N₀)(N/N₀) ^-α with f(N₀) = 14.02 and α =1.44 ± 0.3, from work. The blue error bars are the Poissonian errors for the reference run and are a good representation of the errors in the other models. The diagram shows results with the test log J_v-1, for the simulations Ch 18 512 MDW, Ch 18 512 MDW mol and Ch 18 512 EDW. The color figure can be viewed online. 

Fig. 4 CIV column density distribution function at z = 6.4 and comparison with observational data by ^{Bosman et al. (2017)}, orange diamonds. The error bars are the Poissonian errors for the reference run and are a good representation of the errors in the other models. The left panel shows results with the uniform HM12 and the right panel the results with the test log J_v - 1, both for the simulations Ch 18 512 MDW, Ch 18 512 MDW mol and Ch 18 512 EDW. The color figure can be viewed online. 

The key feature of Figures 2, 3, 4 is that all of them show a notable underproduction of the CIV absorbers at all redshifts when the emissivity is decreased. The calculated CDDFs are always below the observed values, and for high column densities, there is a clear departure from the fitting functions provided by ^{D´Odorico et al. (2013)}. Nonetheless, the synthetic CDDFs show a closer match with the observational values from ^{Codoreanu et al. (2018)}, in particular at high column densities. It is worth noting that ^{Bosman et al. (2017)} data are just upper limits for the CIV-CDDF at 6.2 < z < 7.0. Nevertheless, the computed values in this test are significantly underrepresented.

There is a change in the number of CIV absorbers when the UVB emissivity is lower than the original HM12. This result gives a hint in regard to the number of the CIV absorbers: the test log J_v - 1 is too aggressive, with high ionization states of C lying in the IGM.

Next, we compute the cosmological mass density Ω_CIV. The cosmological mass density of an ion is defined as:

Ωion(z)=H0mioncρcrit∑N(ion,z)nlovΔX, (2)

With m_ion the mass of the ionic species, n_lov the number of lines of view, p_crit the critical density today, and ΔX is expressed above.

Figure 5 shows a comparison of the calculated CIV cosmological mass density at 4 < z < 8 for synthetic absorbers in the range 13.8 < (log NCIV/cm-2) < 15.0 and observations by ^{Pettini et al. (2003)} and ^{Ryan-Weber et al. (2009)} as orange circles, ^{Codoreanu et al. (2018)} as yellow circles, ^{Songaila et al. (2001}, ²⁰⁰⁵⁾ as cyan triangles, ^{Meyer et al. (2019)} as magenta triangles, ^{Simcoe et al. (2011)} as dark green inverted triangles, ^{D'Odorico et al. (2013)} as pink squares, ^{Boksenberg et al. (2015)} as a grey diamond, upper limits from ^{Bosman et al. (2017)} as a purple star and ^{Díaz et al. (2016)} as black pentagons

Fig. 5 CIV cosmological mass density at 4 < z < 8 for 13.8 < log NCIV(cm^-2) < 15.0. Comparison between the simulated data and observations by ^{Pettini et al. (2003)} and ^{Ryan-Weber et al. (2009)}, orange circles; Odoreanu et al. (2018), yellow circles; ^{Songaila et al. (2001}, ²⁰⁰⁵⁾, cyan triangles; ^{Meyer et al. (2019)}, magenta triangles; ^{Simcoe et al. (2011)}, dark green inverted triangles; ^{D'Odorico et al. (2013)}, pink squares; ^{Boksenberg et al. (2015)}, grey diamond; upper limits from ^{Bosman et al. (2017)}, purple star and ^{Díaz et al. (2016)}, black pentagons. Pettini, Ryan-Weber, Codoreanu and Díaz measurements are converted to the Planck cosmology, while for the others this recalibration was not possible due to missing details of the precise pathlength probed. On the left panel the results with the uniform HM12 are presented. The right panel shows CIV in the framework of the test log J^-1. In both cases, the simulations used are Ch 18 512 MDW, Ch 18 512 MDW mol and Ch 18 512 EDW. As a consequence of the low number of absorbers in this column density range, when the UVB normalization is varied, CIV is at least an order of magnitude lower than the case with the original UVB. The color figure can be viewed online. 

As a consequence of the discrepancy presented in Figures 2, 3 and 4, the CIV cosmological mass density in the right hand panel of Figure 5 is at least an order of magnitude below the reference case on the left panel, because the number of absorbers in the range 13.8 < log N_CIV(cm^-2) < 15.0 is underproduced by the simulations in the framework of the modified UVB. The most remarkable difference is visible in the reference run Ch 18 512 MDW, with an order of magnitude shift between the blue curve (on the left) and the light blue one (on the right).

Therefore, strong variations of the UVB around 1 Ryd (specifically, normalization changes in the wavelength range where the transition occurs) seem to have a large impact on the number of CIV absorbers at all redshifts. Figures 6 and 7 draw a comparison of the evolution of CII and CIV in the redshift range 4 < z < 8, with the original (on the left panel) and the test log J_v -1 (on the right). As discussed above, the cosmological mass density of CIV is underrepresented. However, the resulting mass density of CII significantly improves with a softer UVB, indicating that the hypothesis made to perform this test is well-motivated. Effectively, the number of large column density CII absorbers increases and the right panels of Figures 6 and 7 are now compatible with the limits measured by and the most recent estimates made by ^{Cooper et al. (2019)}.

Fig. 6 Evolution of the CII and CIV cosmological mass density when the normalization of the UVB is varied at 1 Ryd (comparison of molecular cooling content). On the left panel the results with the uniform HM12 are presented. The right panel shows results of the test log J_v ^-1. In both cases, the runs used are Ch 18 512 MDW and Ch 18 512 MDW mol. The solid lines mark the evolution of CIV for 13.8 < log NCIV(cm^-2) < 15.0, and the dashed lines CII in the range 13.0 < log NCII(cm^-2) < 15.0. The orange points with errors represent the observational lower limits for CII from ^{Becker et al. (2006)} and the grey arrows the corresponding estimates made by ^{Cooper et al. (2019)}. The latter estimates have been done for z > 5.7 and z < 5.7, shown as right and left arrows, respectively. Although in the case with softer UVB there is no crossover of CII and CIV (due to the low number of CIV absorbers), the mass density of CII matches the limits from ^{Becker et al. (2006)} and ^{Cooper et al. (2019)} in both simulations, and CIV matches the observational detection at z = 5.7 by ^{Díaz et al. (2016)} in the molecular cooling run. The color figure can be viewed online. 

Fig. 7 Evolution of the CII and CIV cosmological mass density when the normalization of the UVB is varied at 1 Ryd (comparison of MDW and EDW feedback prescriptions). On the left panel the results with the uniform HM12 are presented. The right panel shows results of the test log J_v ^-1. In both cases, the runs used are Ch 18 512 MDW and Ch 18 512 EDW. The solid lines mark the evolution of CIV for 13.8 < log NCIV(cm^-2) < 15.0, and the dashed lines CII in the range 13.0 < log NCII(cm^-2) < 15.0. The orange points with errors represent the observational lower limits for CII from ^{Becker et al. (2006)}, and the grey arrows, the estimates from ^{Cooper et al. (2019)}. There is no crossover of CII and CIV at any redshift, because of the different orders of magnitude of the mass densities of these ions. In addition, different feedback prescriptions do not seem to give rise to a remarkable distinction in the evolution of CII. Yet, the plot reveals that a softer UVB effectively favours low ionization states as CII, and brings down the gap between the observations from ^{Becker et al. (2006)} and ^{Cooper et al. (2019)} and the simulated column densities. The color figure can be viewed online. 

The feedback prescription does not play a major role in the evolution of CII in this redshift range, while the molecular cooling run Ch 18 512 MDW mol shows a relatively good agreement with the observational data.

Due to the different orders of magnitude between the calculated mass densities of CII and CIV, the right panels of Figures 6 and 7 show no crossover of the low and high ionization states of carbon. A natural conclusion from this could be that decreasing the intensity of the UV background leads to an improvement in the low ionization states at the expense of a poor calculation of CIV absorbers, which are traditionally well constrained by observations.

As pointed out before, CIV is not well represented by this variation of the UVB, but there is a good improvement in the column densities of CII. However, it is difficult to draw definitive conclusions from the column density relationships, because the number of absorbers depends strongly on the ion, with less CIV synthetic absorbers in the case of log J_v -1 . We suggest variations in the UVB smaller than the original HM12 emissivity, but not below an order of magnitude, so as to preserve the improvements in the low-ionization states while keeping positive results in CIV.

4. Variation of the assumed HI self-shielding prescription

A final test that can be done with our simulations is a variation of the HI self-shielding prescription. ^{García et al. (2017}b) briefly discuss the need for low ionization states self-shielding (SSh) treatment. However, here we try to quantify the impact of a HI self-shielding prescription different from the extensively used one by ^{Rahmati et al. (2013}a) with the HM12 model. For this purpose, we use the new HI SSh model described by ^{Chardin et al. (2018)}. Their best fitting parameters as a function of redshift were obtained with radiative transfer simulations (calibrated with Lyα forest data after the EoR). One of the merits of these models is that they focused on redshifts corresponding to reionization. Instead, the ^{Rahmati et al. (2013a)} HI SSh prescription is valid up to z = 5,which constitutes a limitation when applying this formulation in our models, as commented in ^{García et al. (2017b)}. Here we compare results for the metal ions with the HI SSh treatments by ^{Rahmati et al. (2013a)} and ^{Chardin et al. (2018)} -hereafter R13 and C17, respectively-.

The evolution of the photoinization rate Γphot is calculated using RT codes, such that:

where ΓUVB is the photoionization rate as a function of redshift and it is assumed from the UVB field. Here, f α ₁ , α ₂ , β, and n₀ are free parameters of the model, calculated with RT; the number density of hydrogen n_H and the temperature T are taken directly from the numerical simulation used.

The best fitting parameters of equation 3 as a function of redshift found by R13 and C17 are shown in Table 3.

In contrast to the outcome for the metal ions studied in the previous section, one would expect that HI would be more sensitive to a variation of the self-shielding prescription adopted for this transition. In fact, works by and self-consistently calculate the distribution of neutral hydrogen with RT codes, based on the number density of hydrogen. We compute the following observables for HI: the column density distribution function (CDDF) f_HI, the HI cosmological mass density Ω_HI, and the mass density associated to DLA systems Ω_DLA.

In Figure 8 is shown the HI-CDDF at z = 4, comparing the two self-shielding prescriptions (R13 and C17) and simulations with different molecular cooling contents, Ch 18 512 MDW and Ch 18 512 MDW mol. In addition, we compare our theoretical predictions with observational detections of HI-CDDF at z around 4 in two regimes: the range of column densities 12 < log N_HI (cm^-2) < 22 in the left panel, and a zoom around the DLAs region, 20.3 < log N_HI (cm^-2) < 22 , in the right panel. In the first case, observations by are plotted in grey, in black, in orange, and in purple. Instead, in the DLA zoom, a comparison is drawn with the fitting function proposed by for DLA systems at redshift 4.0-5.5.

Fig. 8 HI column density distribution function at z = 4 comparing two self-shielding prescriptions (^{Rahmati et al. 2013}a; ^{Chardin et al. 2018}) in simulations with different molecular cooling content (Ch 18 512 MDW and Ch 18 512 MDW mol) and the HM12 model. On the left side is shown the simulated HI{CDDF in the range 12 < log NHI (cm^-2º < 22 and a comparison with observations by ^{Prochaska et al. (2005}) in grey, ^{O'Meara et al. (2007}) in black, ^{Crighton et al. (2015}) in orange and ^{Bird et al. (2017}) in purple. In the right panel, the CDDFs are limited to the DLA regime (20:3 < log NHI (cm^-2º < 22) and compared with the fitting function by ^{Prochaska et al. (2009)} for DLA systems at redshift 4.0-5.5 (black dashed line). The color figure can be viewed online. 

One can see a difference between the two self-shielding prescription implemented at z = 4. The ^{Chardin et al. (2018)} model predicts a larger amount of neutral hydrogen hosted in systems in the Lyman-α forest and fewer systems in the sub-DLA and DLA regimes, indicating that the HI mass density should be larger with the ^{Rahmati et al. (2013}a) SSh prescription. In ^{García et al. (2017}b), we find that the largest contribution to Ω_HI comes from systems with large column densities, in particular, DLAs. We predict higher values of Ω_HI with the R13 HI SSh formulation than C17, regardless of the molecular cooling model considered. It is important to remember that at this redshift (z = 4), both formulations are valid, and their best fitting parameters are calibrated with observations. Thus, our conclusions are not limited by different constraints of the HI SSh modelling.

Interestingly, the number of systems in a column density bin at a given absorption path is barely affected by the SSh prescription, but it depends strongly on the chemistry of the molecules; the gap between the number of systems is approximately fixed when comparing two simulations with and without molecules in the Lyman α forest regime (blue and dark red, and light blue and magenta cases, respectively).

Figure 9 (left panel) displays the cosmic mass density of neutral hydrogen with the self-shielding prescriptions by (blue and dark red for the runs Ch 18 512 MDW and Ch 18 512 MDW mol, respectively) and ^{Chardin et al. (2018)} (light blue and magenta corresponding to Ch 18 512 MDW and Ch 18 512 MDW mol, respectively) and compares them with observations by ^{Prochaska et al. (2005)} and ^{Prochaska et al. (2009)}, grey inverted triangles; ^{Zafar et al. (2013)}, pink square; and ^{Crighton et al. (2015}) red stars. As predicted above, the SSh prescription by ^{Rahmati et al. (2013}a) gives rise to a larger amount of neutral hydrogen when compared with the results for Ω_HI of the ^{Chardin et al. (2018)} model. Interestingly, the introduction of this new self-shielding treatment reduces the tension between our models without molecular cooling and the observations at z = 4, and accurately predicts the HI mass density when molecular cooling is taken into account.

Fig. 9 Cosmological mass density of HI with different HI self-shielding treatments: ^{Rahmati et al. (2013}a) and ^{Chardin et al. (2018)}. The left diagram shows the prediction of the HI mass density from simulations with specific molecular cooling content (Ch 18 512 MDW and Ch 18 512 MDW mol) and a comparison with ^{Prochaska et al. (2005}, 2009), grey inverted triangles; ^{Zafar et al. (2013)}, pink square; and ^{Crighton et al. (2015)} red stars. The right panel displays the cosmological mass density associated to neutral hydrogen and DLAs, HI and DLA (solid and dashed lines, respectively) and the observational data for DLA by ^{Bird et al. (2017)}, black stars. The color figure can be viewed online. 

This is by far the most important effect of the implementation of a different self-shielding prescription in our simulations: the amount of HI mass density at redshifts between 4-6 is in better agreement with observational detections.

Finally, we point out that models including molecular cooling give a better prediction of Ω_HI compared with observations, because they take into account the conversion of atomic to molecular hydrogen at very high densities, where the self-shielding of HI is occurring.

On the right hand side of Figure 9; we show a comparison of Ω_HI and Ω_DLA (solid and dashed lines, respectively) with observations for Ω_DLA by ^{Bird et al. (2017)}. The amount of HI mass density hosted in DLAs (dashed lines) converges at high redshift (z ≈ 6) for models with a different SSh prescription when the molecular cooling content is fixed. This is indeed quite an interesting result because the R13 SSh prescription has not been calibrated at this redshift, although it is extensively used in the literature at redshifts higher than 5. One can say that the use of the model is justified at z = 6.

Additionally, there is agreement between the predicted trends for Ω_DLA (with different self-shielding prescriptions and molecular cooling content models) and observations of this quantity by ^{Bird et al. (2017)}.

It is worth noting that the column density distribution functions in Figure 8 show that 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. The molecules cool the surrounding gas and atomic hydrogen is less abundant. This effect is particularly important in the regime of DLAs (^{García et al. 2017}b). Interestingly, molecular chemistry plays a more relevant role in the amount of neutral hydrogen than the HI self-shielding model itself, although the evolution of the ^{Chardin et al. (2018}) model is more dynamic at high redshift than the ^{Rahmati et al. (2013}a) one.

Figure 9 leads to the same conclusion in a different way: simulations without molecular cooling (Ch 18 512 MDW) boost the amount of HI compared with Ch 18 512 MDW mol, regardless of the SSh prescription considered. Then, the cosmological mass density increases by a factor of 2 compared with simulations containing molecular cooling. The effect is more pronounced contrasting HI and DLA cosmic mass density: there are two processes involved, the molecular cooling that converts HI into H2 and reduces the amount of HI in the DLA column density range, but also the evolution of the HI-SSh modeling, that becomes more critical at high redshift ^{Rahmati et al. (2013}a). SSh stops being valid at z = 6, whilst the ^{Chardin et al. (2018)} model is valid up to z = 10.

5. Discussion and Conclusions

This work contains a compilation of variations to the uniform ionizing field assumed in the photoionization models to obtain the results presented in ^{García et al. (2017}b) and ^{García et al. (2017b)}. We find that the large uncertainties in the UVB, especially at high redshift, require better constraints for several probes, and of course, improved UVB models.

In a similar vein, observational detections of metal ionic species could be used to constrain the UVB spectral normalization in the wavelength range where these states occur (100-1000 Å). A careful fine-tuning is required, but most of the steps are clear: a reduction in the normalization of the UV emissivity improves the number of absorbers with low ionization potential energy.

When the emissivity of the assumed UVB is reduced by one dex, we find that the total number of CIV absorbers significantly decreases compared with the original HM12, and that absorbers in the range of the detections by ^{D’Odorico et al. (2013)} are definitely not well represented (e.g. there is a discrepancy of the CIV comoving mass density and CIV CDDF with the observations). Instead, there is a moderate improvement in the calculated number of absorbers for SiIV with a softer UVB, in agreement with findings by ^{Bolton et al. (2011)}.

It is particularly encouraging that the number of absorbers of CII, SiII and OI in the range of the observations also rises when the UVB implemented is softer, independently of the simulation used; therefore, the estimate of the column densities of these low ionization states improves, justifying the hypothesis that led us to do this test. Additional work can be done in this direction, with a more moderate reduction in the emissivity of the UVB that would give better predictions for the column density of low ionization states, and be more compatible with the CIV absorbers incidence rate inferred from the ^{D´Odorico et al. (2013)} observations and our mock spectra for that ion.

We also computed some hydrogen statistics up to z = 6. We limited our calculations to this redshift for two reasons: (1) at high redshift a uniform UVB does not provide a good representation of the rapid evolution of HI during reionization; (2) the self-shielding prescription from ^{Rahmati et al (2013}a) is calibrated with the photoionization rates of HM12. To alleviate the absence of a self-consistent SSh at higher redshifts we compared our previous results with a recent SSh prescription by ^{Chardin et al. (2018)}, whose best-fit parameters were calibrated up to z = 10.

We find that at z = 4, the HI CDDF predicted by the ^{Chardin et al. (2018)} model produces more neutral hydrogen in systems with low column densities (the Lyman-α forest) and fewer systems in the sub-DLAs and DLAs range. This result leads to higher mass densities of HI, Ω_HI, with the ^{Rahmati et al. (2013}a) SSh prescription than with Chardin et al. (2018), independently of the molecular cooling model considered. At z = 4, both formulations (Rahmati et al. 2013a; Chardin et al. 2018) were calibrated with observations. Therefore, our results for the HI CDDF and the cosmological mass density of neutral hydrogen are valid, regardless of the assumptions of each HI self-shielding model.

Regardless of the intensity of the UVB field considered, we find a non-negligible number of ion absorbers with column densities above log N ≥ 15 cm^-2. This prediction from our models can be tested with future observations that, in principle, will be able to detect rarer high column density systems. In addition, an appreciable difference in the number of high column density absorbers would be detected if we traced high gas overdensity, but our pipeline is currently not tracing gas in this region for statistical reasons, since we produce random lines of sight inside the box.

As discussed in , we confirm that different feedback prescriptions (EDW, MDW) do not make a significant difference in the calculation of the number of absorbers of the ionic transitions. The same conclusion is also true for different models with/without molecular cooling for observables associated to the metal ions. Instead, the low temperature and molecular cooling module has proven to be quite important to match the calculated HI CDDF, Ω_HI and Ω_DLA with observations available at z = 4 - 6..

The numerical runs and the pipeline of this paper were calibrated to resemble the physical conditions of the gas residing in the CGM and IGM at z ≈ 6. Most of the observables at that time of the Universe are well represented by our theoretical models. Nonetheless, there are a few caveats in this work: first, more observations of metal absorption lines in the spectra of high redshift quasars are required. Although we have included the largest catalogue of HI and ionic detections at high redshift, more observational data are needed to better constrain the emissivity of the UVB. Second, the numerical resolution is still not enough to fully trace the environments where the low ionization states lie, and simultaneously describe the IGM. A new generation of numerical simulations will improve the description of these species, without sacrificing the physics occurring in the inter- and circumgalactic medium. Finally, we stress that radiative transfer effects are not included in our models and, therefore, we do not follow the progression of reionization, nor the evolution of the HII bubbles or their topology. The implicit assumption is that our boxes (that are small compared to the size of the HII 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, but could be done in the future.

A fundamental issue persists in the field from the numerical point of view: it is extremely challenging to model the low ionization states present in the gas and to provide a good description of the environment where they lie, mainly due to insufficient resolution and to the lack of a proper self-shielding treatment for the ions. We considered only the effect of HI self-shielding (^{Rahmati et al. 2013}a; ^{Chardin et al. 2018}), but did not introduce any self-shielding of low ionization absorbers (which lie in clumpy structures). A first attempt was proposed for DLA systems at low redshift by ^{Bird et al. (2017)}. However, current works miss this component at high redshift because it is still not well understood how high density regions self-shield the gas during the progression of reionization.

In summary, we have compared HI and ion observables available at high redshift and found that most of the results discussed are compatible in the redshift range 4 < z < 8. When discrepancies between observations and synthetic calculations arised, we provided a physical explanation of their nature. It is worth noting that all results from mock spectra will be improved in the future with more observational detections of ion absorbers in the high redshift quasar spectra.