The effect of stellar and AGN feedback on the low redshift Lyman-$\alpha$ forest in the Sherwood simulation suite

TL;DR

This study investigates how different feedback mechanisms in hydrodynamical simulations affect the low redshift Lyman-alpha forest, revealing persistent discrepancies with observations that suggest the need for increased gas temperatures or additional turbulence.

Contribution

It systematically assesses the impact of stellar and AGN feedback on the low redshift Lyα forest within the Sherwood simulation suite, highlighting the limitations of current feedback models.

Findings

Feedback has limited effect on the CDDF and velocity widths.

Discrepancies with COS data suggest the need for higher gas temperatures.

Additional turbulence or heating mechanisms may resolve the observed differences.

Abstract

We study the effect of different feedback prescriptions on the properties of the low redshift ($z\leq1.6$) Ly$\alpha$ forest using a selection of hydrodynamical simulations drawn from the Sherwood simulation suite. The simulations incorporate stellar feedback, AGN feedback and a simplified scheme for efficiently modelling the low column density Ly$\alpha$ forest. We confirm a discrepancy remains between Cosmic Origins Spectrograph (COS) observations of the Ly$\alpha$ forest column density distribution function (CDDF) at $z \simeq 0.1$ for high column density systems ($N_{\rm HI}>10^{14}\rm\,cm^{-2}$), as well as Ly$\alpha$ velocity widths that are too narrow compared to the COS data. Stellar or AGN feedback -- as currently implemented in our simulations -- have only a small effect on the CDDF and velocity width distribution. We conclude that resolving the discrepancy between the COS data and simulations requires an increase in the temperature of overdense gas with $\Delta=4$--$40$, either through additional He$\,\rm \scriptstyle II\ $ photo-heating at $z>2$ or fine-tuned feedback that ejects overdense gas into the IGM at just the right temperature for it to still contribute significantly to the Ly$\alpha$ forest. Alternatively a larger, currently unresolved turbulent component to the line width could resolve the discrepancy.

Tables

Table 1. Table 1: The hydrodynamical simulations used in this work. The first column uses the naming convention L-N-param . Here L 𝐿 L is the box size in units of h − 1 ​ cMpc superscript ℎ 1 cMpc h^{-1}{\rm cMpc}\, , N 𝑁 N is the cube root of the total number of gas particles and -param is associated with the subgrid treatment of star formation and feedback. Here − p ​ s ​ 13 𝑝 𝑠 13 -ps13 refers to the star formation and energy driven outflow model of Puchwein & Springel ( 2013 ) , while -ps13+agn also includes AGN feedback based on Sijacki et al. ( 2007 ) . All other models instead convert all cold ( T < 10 5 ​ K 𝑇 superscript 10 5 K T<10^{5}\rm\,K ), dense ( Δ > 1000 Δ 1000 \Delta>1000 ) gas into collisionless particles. The remaining columns give the mass of each dark matter and gas particle in h − 1 ​ M ☉ superscript ℎ 1 subscript 𝑀 ☉ h^{-1}M_{\sun} and the gravitational softening length in h − 1 ​ ckpc superscript ℎ 1 ckpc h^{-1}\rm{ckpc} . All simulations are performed to z = 0 𝑧 0 z=0 .

Name	$M_{\mathrm{dm}}$	$M_{\mathrm{gas}}$	$l_{\mathrm{soft}}$
	$\left[h^{-1}{M}_{☉}\right]$	$\left[h^{-1}{M}_{☉}\right]$	$\left[h^{-1}\mathrm{ckpc}\right]$
80-512-ps13	$2.75\times {10}^8$	$5.10\times {10}^7$	$6.25$
80-512-ps13+agn	$2.75\times {10}^8$	$5.10\times {10}^7$	$6.25$
80-512	$2.75\times {10}^8$	$5.10\times {10}^7$	$6.25$
80-1024	$3.44\times {10}^7$	$6.38\times {10}^6$	$3.13$
80-256	$2.20\times {10}^9$	$4.08\times {10}^8$	$12.50$
40-512	$3.44\times {10}^7$	$6.38\times {10}^6$	$3.13$
20-256	$3.44\times {10}^7$	$6.38\times {10}^6$	$3.13$

Table 2. Table 2: A comparison of baryon phases in the simulations with different treatments for cold, dense gas and feedback at z = 0 , 0.1 𝑧 0 0.1 z=0\,,0.1 and 1.6 1.6 1.6 where Δ th = 97.2 ,  65.9 subscript Δ th 97.2 65.9 \Delta_{\rm th}=97.2,\,65.9 and 62 62 62 , respectively. The percentage of the total gas mass is divided into four categories following Davé et al. ( 2010 ) .

Name	Model	Diffuse	WHIM	Hot Halo	Condensed	Stars
		K		$\Delta >{\Delta }_{\mathrm{th}},T>{10}^5$K	K
$z=0.1$	80-512	36.5	21.3	17.1	0.9	24.2
	80-512-ps13	36.6	24.3	22.4	9.1	7.7
	80-512-ps13+agn	36.5	36.9	20.0	3.3	3.3
$z=1.0$	80-512	51.7	15.4	12.6	1.9	18.4
	80-512-ps13	51.9	17.0	16.3	10.6	4.3
	80-512-ps13+agn	51.6	19.2	17.7	8.5	3.0
$z=1.6$	80-512	61.8	11.3	8.9	2.8	15.1
	80-512-ps13	62.0	12.3	11.3	11.5	2.8
	80-512-ps13+agn	62.0	12.7	11.7	11.1	2.4

Table 3. Table 3: As for Table 2 , but now including all remaining models using the QUICKLYA prescription for the removal of dense gas from the simulations. We repeat the entry for the 80-512 model here for ease of comparison.

Name	Model	Diffuse	WHIM	Hot Halo	Condensed	Stars
		K		$\Delta >{\Delta }_{\mathrm{th}},T>{10}^5$K	K
$z=0.1$	80-1024	35.0	20.1	15.3	0.7	28.8
	80-512	36.5	21.3	17.1	0.9	24.2
	80-256	39.7	22.6	20.0	1.2	16.6
	40-512	34.6	20.5	15.3	0.8	28.8
	20-256	35.3	19.4	15.2	0.7	29.3
$z=1.0$	80-1024	49.6	14.5	11.1	1.6	23.2
	80-512	51.7	15.4	12.6	1.9	18.4
	80-256	55.4	16.4	14.9	2.1	11.2
	40-512	49.2	15.0	11.1	1.5	23.2
	20-256	49.7	13.5	11.5	1.6	23.7
$z=1.6$	80-1024	59.3	10.7	7.7	2.3	20.0
	80-512	61.8	11.3	8.9	2.8	15.1
	80-256	65.7	12.3	10.6	2.9	8.4
	40-512	58.9	11.1	7.8	2.2	19.9
	20-256	59.2	10.2	7.9	2.3	20.4

Equations

$$N_{\mathrm{HI}}=\frac{2.44\times 10^{13}\mathrm{c}{\mathrm{m}}^{-2}}{{\Gamma }_{-13}}{\left(\frac{\Delta }{10}\right)}^{3/2}{\left(\frac{T_4}{2}\right)}^{-0.26}{\left(\frac{1+z}{1.1}\right)}^{9/2}.$$

$$b_{\mathrm{HI}}=18.2\mathrm{km}{\mathrm{s}}^{-1}{\left(\frac{{\mathrm{T}}_4}{2}\right)}^{1/2}.$$

Full text

††pagerange: The effect of stellar and AGN feedback on the low redshift Lyman-$\alpha$ forest in the Sherwood simulation suite–LABEL:lastpage

1 Introduction

The Ly$\rm{\alpha}$ forest is an indispensable tool for probing the evolution and distribution of gas in the intergalactic medium (see Meiksin, 2009; McQuinn, 2016, for reviews). At intermediate redshifts ($2\la z\la 5$), the majority of the intergalactic medium (IGM) by volume is in the form of warm ($T\sim 10^{4}\rm\,K$) diffuse gas that is kept photo-ionised by the metagalactic ultraviolet background (UVB). The latest hydrodynamical simulations of the IGM that model the Ly$\rm{\alpha}$ forest of absorption produced by this gas are in very good agreement with a wide range of spectroscopic data from optical, ground-based telescopes (e.g. Bolton et al., 2017). However, at lower redshifts ($z\la 2$), due to the interplay between cosmic expansion, the growth of structure and the declining intensity of the UVB, the Ly$\rm{\alpha}$ forest thins out and becomes more transmissive (Danforth et al., 2016). By $z\sim 0$ a significant fraction of the baryons in the IGM are in a collisionally ionised, shock heated phase with temperatures $T\sim 10^{5}$–$10^{7}\rm\,K$ (Cen & Ostriker, 1999). In the present day Universe, the Ly$\rm{\alpha}$ forest therefore traces less than one half of the remaining diffuse gas in the IGM (Lehner et al., 2007; Shull et al., 2012).

Over the past two decades, the properties of Ly$\rm{\alpha}$ absorption arising from the low redshift IGM have been studied with observations made by a series of space-based UV spectrographs (Bahcall et al., 1993; Weymann et al., 1998; Janknecht et al., 2006; Kirkman et al., 2007; Williger et al., 2010; Tilton et al., 2012). Most recently, the Cosmic Origins Spectrograph (COS, Green et al., 2012), mounted on the Hubble Space Telescope (HST) has collected many more active galactic nuclei (AGN) spectra with improved signal-to-noise, providing new opportunities for absorption line studies of the low redshift Ly$\rm{\alpha}$ forest (Shull et al., 2014; Danforth et al., 2016; Pachat et al., 2016). Some of the key questions addressed by recent investigations are related to the nature of these Ly$\rm{\alpha}$ absorbers, and their relationship to the evolution of the UVB and galactic feedback (see e.g. Davé et al., 2010; Tepper-García et al., 2012; Kollmeier et al., 2014; Shull et al., 2015; Khaire et al., 2016; Viel et al., 2017; Gurvich et al., 2017).

The mostly widely used measurements of Ly$\rm{\alpha}$ absorbers at $z<2$ are the column density distribution function (CDDF), the velocity (Doppler) widths of the absorption lines and the evolution of the absorption line number density with redshift. The CDDF shape and normalisation is determined by the density distribution of absorbers and the intensity of the UVB. The velocity width distribution provides constraints on the gas temperature. Furthermore, both observables are in principle sensitive to galactic feedback. In this context, hydrodynamical simulations have significantly aided our understanding and interpretation of these observational data (Rauch et al., 1997; Theuns et al., 1998a; Davé et al., 1999; Richter et al., 2006; Paschos et al., 2009; Oppenheimer & Davé, 2009; Davé et al., 2010; Rahmati et al., 2013).

Of particular note is the recent comparison of the CDDF for Ly$\rm{\alpha}$ absorbers at $z\simeq 0.1$ observed with COS (Danforth et al., 2016) to hydrodynamical simulations. Kollmeier et al. (2014) reported that the UVB intensity needed to match the COS measurement of the CDDF is a factor of five larger than predicted by empirically calibrated UVB models (Haardt & Madau, 2012, HM12). The claimed “photon underproduction crisis” prompted further investigation into this apparent discrepancy. Large uncertainties in the production rate of ionising photons by quasars and star forming galaxies, and possibly also the role of AGN feedback, have been suggested as the likely cause (Khaire & Srianand, 2015; Shull et al., 2015; Wakker et al., 2015; Gurvich et al., 2017; Gaikwad et al., 2017a).

Viel et al. (2017) (hereafter V17) also recently considered this problem using hydrodynamical simulations including Illustris (Vogelsberger et al., 2014) and a sub-set of the Sherwood simulations (Bolton et al., 2017). They concluded the HM12 photo-ionisation rate needed to be increased by a factor of two in the hydrodynamical simulations to match the observed CDDF, smaller than the factor of five found by Kollmeier et al. (2014). However, in contrast to other recent studies, V17 also focussed on the distribution of Ly$\rm{\alpha}$ line velocity widths (see also Gaikwad et al., 2017b). V17 found that the simulations struggled to reproduce the Ly$\rm{\alpha}$ velocity width distribution at $z\simeq 0.1$, producing line widths that were too narrow compared to the COS data. This was attributed to gas in the simulations that was either too cold to produce the observed line widths, or too highly (collisionally) ionised to produce Ly$\rm{\alpha}$ absorption due to vigorous AGN feedback heating the gas to $T>10^{5}\rm\,K$ (particularly in the Illustris simulation). V17 concluded that the Ly$\rm{\alpha}$ line width distribution at $z\simeq 0.1$ is a valuable diagnostic for models of feedback in the low redshift Universe.

The goal of the present work is to the extend the V17 analysis, and present in full the low redshift ($0<z<2$) IGM models that V17 used from the Sherwood simulation suite (Bolton et al., 2017). Following V17, we perform the same Voigt profile analysis on the simulations and COS data to minimise potential biases in the comparison due to line fitting. This furthermore enables a detailed assessment of the effect that our stellar and AGN feedback implementation has on the Ly$\rm{\alpha}$ forest.

The paper is organised as follows. In Section 2, we briefly describe the simulations used in this work. In Section 3, we compare the simulations to current observational data at $z<2$ and assess the impact of different galactic feedback prescriptions on Ly$\rm{\alpha}$ forest observables, including the the Ly$\rm{\alpha}$ velocity width distribution. The physical origin of the absorbers is examined in Section 4. Finally, we summarise our results in Section 5. A series of convergence tests are given in the Appendix. Throughout this paper we refer to comoving distances with the prefix “c”. The cosmological parameters used throughout this work are $\Omega_{\rm m}=0.308$, $\Omega_{\rm\Lambda}=0.692$, $h=0.678$, $\Omega_{\rm b}=0.0482$, $\sigma_{\rm 8}=0.829$ and $n=0.961$ consistent with the best fit $\Lambda{\rm CDM}+Planck+WP+highL+BAO$ cosmological parameters (Planck Collaboration et al., 2014).

2 Methodology

2.1 Hydrodynamical simulations

The simulations used in this work are summarised in Table 1. These were performed using P-Gadget-3, a modified version of the publicly available code Gadget-2 (Springel, 2005). The initial conditions were generated at $z=99$ on a Cartesian grid using the N-GenIC code (Springel et al., 2005) using transfer functions generated by CAMB (Lewis et al., 2000). A detailed description of these models may be found in the Sherwood overview paper (Bolton et al., 2017).

The photo-ionisation and photo-heating rates are calculated using a spatially uniform UVB model for emission from star forming galaxies and quasars (Haardt & Madau, 2012) applied in the optically thin limit. A small increase to the He\,\rm\scriptstyle II\photo-heating rate, $\epsilon_{\rm HeII}=1.7\epsilon_{\rm HeII}^{\rm HM12}$ at $2.2<z<3.4$ was applied to better match observational measurements of the IGM temperature at $z>2$ (Becker et al., 2011). The ionised fractions for hydrogen and helium are obtained assuming ionisation equilibrium. Metal line cooling is not included, although this is not expected to impact significantly on the Ly$\rm{\alpha}$ forest at the redshifts considered here (Tepper-García et al., 2012). In five of the simulations, gas particles with temperatures $T<10^{5}$K and densities $\Delta=\rho/\langle\rho\rangle>1000$ are converted into collisionless star particles (the QUICKLYA method, see Viel et al., 2004) to speed up the calculation. We will show this choice makes very little difference to Ly$\rm{\alpha}$ forest systems at low column densities, $N_{\rm HI}<10^{14.5}\rm\,cm^{-2}$, in our analysis. However, galactic winds and AGN feedback may have an impact on high column density lines that probe the (circumgalactic) gas in the vicinity of dark matter haloes. Therefore, two models that follow star formation and feedback are included. The first model (80-512-ps13) includes energy driven galactic outflows following Puchwein & Springel (2013). The star formation model is based on Springel & Hernquist (2003), but assumes a Chabrier rather than Salpeter initial mass function (IMF) and a galactic wind velocity that is directly proportional to the escape velocity of the galaxy. These choices increase the available supernovae feedback energy by a factor of two, as well as raising the mass-loading of winds in low-mass galaxies. The second run (80-512-ps13+agn) also incorporates AGN feedback, again following Puchwein & Springel (2013). This AGN feedback model is based on Sijacki et al. (2007). In the “quasar” mode, when accretion rates are above $0.01$ of the Eddington rate, only $0.5$ per cent of the accreted rest mass energy is thermally coupled to the surrounding gas. For lower accretion rates a more efficient coupling of AGN jets in the “radio” mode is assumed and $2$ per cent of the rest mass energy is used for recurrently injecting hot AGN bubbles.

The evolution of the low-density IGM thermal history at $z<2$ in a sub-set of the simulations is displayed in Figure 1. This is parameterised as a power law temperature-density ($T-\Delta$) relation for the low density ($\Delta<1$–$10$) IGM, $T=T_{0}\Delta^{\gamma-1}$. Here $\Delta=\rho/\langle\rho\rangle$ is the gas density normalised by the background density, $T_{0}$ is the gas temperature at mean density and $\gamma-1$ is the slope of the relation (Gnedin & Hui, 1998; Upton Sanderbeck et al., 2016). Note that the 80-512-ps13, 80-512-ps13+agn and 80-512 models have very similar thermal histories; this suggests gas at the mean background density is not strongly impacted by feedback. Additionally, by comparing these models to the 80-1024 simulation (solid black line), it is clear that the lower resolution simulations are not converged, with volume weighted temperatures at mean density $\sim 800\rm\,K$ hotter in the 80-512 runs (further details may be found in the Appendix).

There are currently no detailed observational constraints on the temperature of the low density IGM at $z<1.6$, but there are a few measurements from optical data at $1.6\leq z\leq 2$. In the upper left panel of Figure 1 the red triangles show the IGM temperature at mean density from Boera et al. (2014). This study uses the curvature of the Ly$\rm{\alpha}$ forest transmitted flux to measure the gas temperature at the characteristic density, $T(\bar{\Delta})$, probed by the Ly$\rm{\alpha}$ forest absorption. We derive values for $T_{0}$ from these data assuming the $\gamma-1$ for the 80-1024 model shown in the lower left panel, such that $T_{0}=T(\bar{\Delta})/\bar{\Delta}^{\gamma-1}$. Additionally, the cyan squares display the Boera et al. (2014) measurements obtained assuming the effective optical depth, $\tau_{\rm eff}\,$ (where \tau_{\rm eff}\,$$=-\ln\langle F\rangle and $\langle F\rangle$ is the mean transmission in the Ly$\rm{\alpha}$ forest), evolution used in an earlier curvature analysis performed by Becker et al. (2011) (E. Boera, private communication). The effective optical depth sets the characteristic gas density probed by the Ly$\rm{\alpha}$ forest, and this comparison gives an indication of the likely systematic uncertainty in $\tau_{\rm eff}$ arising from line of sight variance. For completeness, we also directly compare the simulations to the $T(\bar{\Delta})$ measurements in the right panel of Figure 1. The measurements bracket all the simulations at $z>1.6$. The models are also formally consistent with the constraints on $\gamma-1$ from Schaye et al. (2000) and Ricotti et al. (2000) shown in the lower left panel of Figure 1. These were obtained from an analysis of the lower envelope of the column density – velocity width plane obtained from Voigt profile fits to the Ly$\rm{\alpha}$ forest, although the error bars on these measurements are very large. Note also the simulations do not follow non-equilibrium ionisation effects, and therefore produce a temperature-density relation slope that is steeper by $\Delta(\gamma-1)\sim 0.1$ at $z=1$–$2$ compared to the recent study by Puchwein et al. (2015).

In all cases the simulated $T_{0}$ monotonically decreases toward lower redshift as the IGM cools due to adiabatic expansion, reaching $\sim 4000\rm\,K$ by $z=0$. The agreement with measurements of the IGM thermal state at $1.6<z<2$ suggests the photo-heating rates used in the UVB model are broadly correct at these (and higher) redshifts. Significantly higher temperatures at mean density than those predicted by the models at $z<1.6$ would require an additional, unidentified source of heating in the low density IGM.

2.2 Mock Ly$\rm{\alpha}$ forest spectra

Mock Ly$\rm{\alpha}$ forest spectra were obtained by extracting $5000$ lines of sight parallel to each of the box axes on regularly spaced grids of $50^{2}$, $40^{2}$ and $30^{2}$, each with $2048$ pixels. The H\,\rm\scriptstyle I\Ly$\rm{\alpha}$ optical depths, $\tau^{\alpha}_{\rm HI}$, are obtained using the method described by Theuns et al. (1998b) combined with the Voigt profile approximation described by Tepper-García (2006). The transmitted flux, $F$, in each pixel is then $F=e^{-\tau^{\rm\alpha}_{\rm HI}}$, and the effective optical depth evolution with redshift is averaged over all sight lines. In our analysis we also make use of the fact that the pixel optical depths are inversely proportional to the H\,\rm\scriptstyle I\photo-ionisation rate, $\tau^{\alpha}_{\rm HI}\propto\Gamma_{\rm HI}^{-1}$. This enables the optical depths (and hence the $\tau_{\rm eff}\,$ evolution) to be rescaled to match different UVB models in post-processing (e.g. Rauch et al., 1997).

In Figure 2 (left panel), the $\tau_{\rm eff}$ evolution for a variety of sub-grid physics and UVB models are compared with Ly$\rm{\alpha}$ forest observations. The results shown by the black curves use the HM12 UVB model, whereas the red triple-dot dashed and blue dashed curves are obtained by scaling the simulated optical depths to match the H\,\rm\scriptstyle I\photo-ionisation rates in the earlier Haardt & Madau (2001, HM01) and Faucher-Giguère et al. (2009) UVB models for emission from star-forming galaxies and quasars. The evolution of $\tau_{\rm eff}\,$ in the 80-1024 model when using the HM01 UVB is in reasonable agreement with observational measurements from Kirkman et al. (2007) (except the dip at $z\simeq 1.4$) and Kim et al. (2007) at $1.7<z<2.0$, but underpredicts $\tau_{\rm eff}\,$ at $z<0.4$. In contrast, the fiducial HM12 UVB predicts $\tau_{\rm eff}\,$ values that are too high, particularly $z<1.5$; a closely related result was noted by Kollmeier et al. (2014) when analysing the H\,\rm\scriptstyle I\CDDF. None of the UVB models is able to reproduce a $\tau_{\rm eff}\,$ evolution that is in good agreement with the more recent COS data from Danforth et al. (2016) at $z\leq 0.4$. There is also a difference between these data and the measurements from Kirkman et al. (2007) at $z\sim 0.4$, possibly due to systematic differences in the measurements and cosmic variance. The Faucher-Giguère et al. (2009) UVB is closer to these data than HM12 at $z<0.5$, but overpredicts $\tau_{\rm eff}\,$ at all redshifts including $z=2$. Interestingly, the differences between models with different feedback implementations are small, indicating that it is the uncertain amplitude of the UVB rather than galactic feedback that is largely responsible for this discrepancy (see also figure 1 in Davé et al., 2010).

This is further illustrated in the right panel of Figure 2 where observational estimates of the H\,\rm\scriptstyle I\photoionisation rate, $\Gamma_{\rm-13}=\Gamma_{\rm HI}/10^{\rm-13}{\rm s}^{\rm-1}$, are compared to the three UVB models. The measurements by Shull et al. (2015); Gaikwad et al. (2017a) and V17, based on an analysis of the COS Ly$\rm{\alpha}$ CDDF, are in mutual agreement and predict photo-ionisation rates at $z<0.4$ intermediate between HM01 and HM12. The value proposed by Kollmeier et al. (2014) (cyan triangle) lies very close to the HM01 UVB at $z=0.1$. The red inverted triangles with large error bars display earlier constraints on $\Gamma_{\rm HI}$ obtained from the proximity effect by Scott et al. (2002). The blue rectangle corresponds to the recent, independent constraint from Fumagalli et al. (2017) using observations of H$\alpha$ fluorescence from H\,\rm\scriptstyle I\in a nearby galactic disk (see also Adams et al., 2011).

Differing assumptions regarding the relative contributions of star forming galaxies and quasars and their spectral shape are largely responsible for differences between existing UVB models. Given these considerable differences, efforts to calibrate the H\,\rm\scriptstyle I\photo-ionisation rate to match $\tau_{\rm eff}\,$ and CDDF measurements at $z<0.4$ (e.g. Khaire & Srianand, 2015; Madau & Haardt, 2015), as well as renewed attempts to obtain measurements of $\Gamma_{\rm HI}$ in the redshift interval $0.4<z<2$, are desirable.

3 Voigt profile analysis

We now proceed to present the results of the Voigt profile analysis of our mock Ly$\rm{\alpha}$ forest spectra. We fit lines at three redshifts, $z=[0.1,\,1,\,1.6]$. When comparing simulations to each other or to data compilations that have a range of resolution and signal-to-noise properties, the mock spectra are post-processed by convolving them with an instrument profile with a Full Width Half Maximum (FWHM) of $7$ ${\rm km\,s}^{-1}$and rebinned to a pixel size of 3 ${\rm km\,s}^{-1}$. A uniform Gaussian distributed noise is added with a signal-to-noise of S/N$=50$ per pixel. The mean transmission of the spectra are rescaled to match \tau_{\rm eff}\,$$=[0.012,\,0.043,\,0.086] at $z=[0.1,\,1,\,1.6]$. These values correspond to the 80-1024 model $\tau_{\rm eff}\,$ evolution for the HM01 UVB for emission from star forming galaxies and quasars, which is in good agreement with the Kirkman et al. (2007) and Kim et al. (2007) $\tau_{\rm eff}$ measurements in the left panel of Figure 2.

The exception to this is when we directly compare the simulations to the V17 COS measurements of the CDDF and velocity width distribution at $z=0.1$. In this case, we follow V17 and convolve the spectra with the LF1 G130M COS line spread function including scattering at wavelength $1350\,$Å111http://www.stsci.edu/hst/cos/performance/spectral_resolution/. The spectra are rebinned to pixels of $\sim 7.2\rm\,km\,s^{-1}$ and a Gaussian distributed signal-to-noise of S/N$=30$ per resolution element (or equivalently $\sim 20$ per pixel) is then added. The effective optical depth is scaled to \tau_{\rm eff}\,$$=0.021 ($\langle F\rangle=0.979$). This corresponds to the value V17 required to match the observed CDDF in the range $10^{13}\leq N_{\rm HI}/{\rm cm}^{-2}\leq 10^{14}$.

The COS observations of the CDDF and velocity width distribution at $z=0.1$ are as described by V17; further details regarding the data reduction and AGN spectra can also be found in Wakker et al. (2015) and Kim et al. (in prep). The data comprises of 704 H\,\rm\scriptstyle I\lines within column densities $10^{12.5}\leq N_{\rm HI}/\rm cm^{-2}\leq 10^{14.5}$ with a total redshift path length of $\Delta z=4.991$, covering the Ly$\rm{\alpha}$ forest at $0<z<0.2$. The Voigt profile fits are obtained using Ly$\rm{\alpha}$ absorption only.

The fitting of Voigt profiles to the mock Ly$\rm{\alpha}$ spectra is performed with VPFIT (Carswell & Webb, 2014), which deconvolves the (already convolved) mock spectra with the instrument profile to obtain intrinsic line widths. Importantly, this approach matches that used to fit the COS observational data at $z=0.1$. A total path length of 10^{5}$$h^{-1}{\rm cMpc}\, is fitted for each model in segments of 20$$h^{-1}{\rm cMpc}\, to obtain our simulated line lists. All lines within 50\,$${\rm km\,s}^{-1} of the start and end of each segment are ignored to avoid spurious line fits due to edge effects. Finally, to estimate the sample variance in the COS observational data, we adopt the approach of Rollinde et al. (2013) and bootstrap resample the mock spectra with replacement over the observed path length $\Delta z=4.991$. The sampling is performed 1000 times over the total simulated path length of $10^{5}h^{-1}\rm\,cMpc$.

3.1 The CDDF

Figure 3 displays the comparison between the CDDF, defined as the number of absorbers per unit log column density per unit redshift, $f(N_{\rm HI},z)=\partial^{2}N/\partial\log N_{\rm HI}\partial z$, from the Sherwood models and the COS observations described by V17 at $z=0.1$ (see also figure 1 in V17). Note that only the 80-512-ps13+agn model was examined in detail by V17; in this work we present a comparison of the additional Sherwood models in which star formation, feedback and numerical resolution have been changed while other parameters remain fixed. We also display the independent measurements of the CDDF from Danforth et al. (2016) (red triangles) and Gaikwad et al. (2017a) (orange circles), although we caution that differences in the line fitting methodologies used in these studies will bias any direct comparison with our simulations.

The simulations fail to capture the correct number of absorbers at $N_{\rm HI}<10^{13.2}\rm\,cm^{-2}$, by as much as a factor $\sim 2$ at $10^{13}\rm\,cm^{-2}$. Note that the simulations are well converged with mass resolution – this may be assessed by comparing the 80-512 (blue dashed curve) and 80-1024 (orange dot-dashed curve) models. However, the CDDF is very sensitive to the signal-to-noise and spectral resolution at low column densities (cf. Figure 4 where S/N$=50$ and FWHM=$7\rm\,km\,s^{-1}$), and the observational data are incomplete at $N_{\rm HI}<10^{13}\rm\,cm^{-2}$. This discrepancy therefore reflects differences between the true noise on the data and the idealised, uniform noise added to the mock spectra and is probably not significant. In contrast, the agreement between the V17 data and simulations at $10^{13.2}<N_{\rm HI}/\rm cm^{-2}<10^{14}$ is within the $2\sigma$ uncertainty from sample variance, shown by the grey shading in the lower panel. At $N_{\rm HI}>10^{14}\rm\,cm^{-2}$, however, the number of absorbers are underpredicted by a factor $\sim 2$ and the slope of the observed CDDF is shallower than the simulated prediction. A similar discrepancy was present in the recent analysis by Shull et al. (2015) using simulations performed with the Eulerian hydrodynamical code enzo, as well as in the study by Gurvich et al. (2017) with the moving mesh code Arepo, suggesting this discrepancy is not unique to P-Gadget-3 simulations. In contrast Gaikwad et al. (2017b) recently found better agreement with the CDDF at $N_{\rm HI}>10^{14}\rm\,cm^{-2}$ using Gadget-2 adiabatic hydrodynamical simulations post-processed to incorporate radiative cooling, photo-heating and ionisation. Note, however, the effect of gas pressure on the H\,\rm\scriptstyle I\distribution will not be correctly captured in these models, complicating interpretation of this result.

Lastly, it is clear that the various feedback scenarios explored by the simulations impact very little overall on the CDDF at low column densities, particularly at $N_{\rm HI}<10^{14}\rm\,cm^{-2}$. Outflows have a similarly small impact on the low column Ly$\rm{\alpha}$ forest at $z>2$; these systems are associated with the diffuse IGM rather than circumgalactic gas (Theuns et al., 2002; Viel et al., 2013, and see Section 4). Any differences between the feedback models are within the expected 1$\sigma$ sample variance.

This is further emphasised in Figure 4 where we examine the Ly$\rm{\alpha}$ absorption systems at column densities $10^{12.2}\leq N_{\rm HI}/\rm{cm}^{-2}\leq 10^{17}$ at $z=[0.1,\,1.0,\,1.6]$. Note that we adopt S/N$=50$ per pixel and a instrument resolution $\rm FWHM=7$ ${\rm km\,s}^{-1}$ here. The lower panels display the residual with respect to the 80-512-ps13 model, and grey shading shows the Poisson error. The QUICKLYA method (80-512, blue dashed-line) agrees well with the more sophisticated sub-grid physics models for weak absorption systems with $N_{\rm HI}<10^{13.5}\rm\,cm^{-2}$ at $z=0.1$ and $N_{\rm HI}<10^{14.5}\rm\,cm^{-2}$ at $z=1.6$. However, higher column density lines are either underpredicted or absent relative to the other models. Multi-phase star formation and stellar feedback (80-512-ps13, black curve) results in a factor of 2-5 more absorption systems at $N_{\rm HI}\ga 10^{14}\rm\,cm^{-2}$. A similar result was noted by Bolton et al. (2017) for the CDDF at $2<z<3$, although the H\,\rm\scriptstyle I\column density threshold where QUICKLYA underpredicted the CDDF was higher, $N_{\rm HI}\simeq 10^{14.5-15}\rm\,cm^{-2}$. The lack of $N_{\rm HI}\ga 10^{14}\rm\,cm^{-2}$ absorbers is in part because of missing cold dense gas in the QUICKLYA model, but also because winds redistribute gas into the diffuse IGM and introduce peculiar velocities that can increase the equivalent widths of the lines (Davé et al., 2010; Meiksin et al., 2015).

Finally, including AGN feedback in the Sherwood models yields similar results for the CDDF relative to the stellar feedback only model, except at $z=0.1$ where the incidence of high column density systems is slightly lower (although the total number of absorbers here is small). A stronger effect was noted by Gurvich et al. (2017), who found that the more vigorous AGN feedback in the Illustris simulation suppresses the CDDF at $N_{\rm HI}=10^{14}$–$10^{15}\rm\,cm^{-2}$ by a factor of two due to changes in the density and temperature of the gas associated with these absorbers. The higher temperatures reduce the Ly$\rm{\alpha}$ opacity by lowering the radiative recombination rate and (for $T\ga 10^{5}\rm\,K$) collisionally ionising the gas. In contrast, Tepper-García et al. (2012, see their figure C1) find stellar and AGN feedback have a negligible effect on the CDDF at $z=0.25$ using simulations drawn from the OWLS project (Schaye et al., 2010). These differences likely reflect variations in the sub-grid physics, as the precise impact of feedback will depend on how efficiently gas is ejected from dark matter haloes.

In contrast to the present study, Tepper-García et al. (2012) also found reasonably good agreement with the H\,\rm\scriptstyle I\CDDF measurement at $10^{14}\leq N_{\rm HI}/\rm cm^{-2}\leq 10^{14.5}$ from Lehner et al. (2007) based on an earlier compilation of FUSE (Far Ultraviolet Spectroscopic Explorer, FWHM$\sim 20\rm\,km\,s^{-1}$) and HST/STIS (Space Telescope Imaging Spectrograph, FWHM$\sim 7\rm\,km\,s^{-1}$) observations at $z\la 0.4$. However, the models in Tepper-García et al. (2012) underpredicted the incidence of saturated lines at $N_{\rm HI}\sim 10^{14.5}\rm\,cm^{-2}$; their Ly$\rm{\alpha}$ only line fits underestimated the column density of these absorbers since the Lehner et al. (2007) HST/STIS data are obtained from Ly$\rm{\alpha}$ and higher order Lyman lines. The Lehner et al. (2007) H\,\rm\scriptstyle I\CDDF furthermore only includes absorbers with velocity widths b_{\rm HI}<40$${\rm km\,s}^{-1} and relative errors less than 40 per cent on $N_{\rm HI}$ and $b_{\rm HI}$ (i.e. figure 14 in Lehner et al., 2008). We have confirmed we obtain similar agreement with the Lehner et al. (2007) CDDF measurements and our simulations (for $\rm S/N=50$ and FWHM=$7\rm\,km\,s^{-1}$) to that obtained by Tepper-García et al. (2012) when selecting only absorbers with $b_{\rm HI}<40\rm\,km\,s^{-1}$.

3.2 The velocity width distribution

We next turn to consider the distribution of Ly$\rm{\alpha}$ line velocity widths, $b_{\rm HI}$. Figure 5 displays a comparison between the V17 COS measurements and simulations at $z=0.1$. We again show the independent measurements from Danforth et al. (2016) (red triangles) and Gaikwad et al. (2017a) (orange circles) for comparison. Overall, the models slightly underpredict the incidence of broad lines at $b_{\rm HI}=40$–$60$ ${\rm km\,s}^{-1}$ and overpredict at b_{\rm HI}<30$${\rm km\,s}^{-1}, although generally the data lie within the expected 95 per cent scatter arising from sample variance in the reference model. It is important to note, however, that unlike the CDDF (where the effect of feedback is comparable to or smaller than the effect of mass resolution) the line widths are not well converged with mass resolution (see also the Appendix later). The distribution from the 80-1024 simulation (orange dot-dashed curve) differs significantly from the 80-512 model (blue dashed curve), particularly at b_{\rm HI}\la 30\,$${\rm km\,s}^{-1}, with a median value of 30.8$${\rm km\,s}^{-1} (cf. 33.5$${\rm km\,s}^{-1} for the 80-512 model and 36.2$${\rm km\,s}^{-1} for the COS data). Higher resolution simulations will therefore predict line widths that are narrower than observed. As pointed out by V17, this suggests there may be additional physics that broadens the absorption lines that is not incorporated into the simulations at present.

We may nevertheless still examine the relative difference between the different feedback models. The incorporation of multi-phase star formation and supernovae-driven winds (black curve) makes very little difference compared to the QUICKLYA model. The addition of AGN feedback (red dotted curve) also has a small effect, increasing the number of broad lines at $b_{\rm HI}>50\rm\,km\,s^{-1}$ with a median222This is consistent with V17, who reported an increase by roughly $2$ ${\rm km\,s}^{-1}$ in the median b-parameter when including AGN feedback. Note that V17 incorrectly reported this quantity as the peak of the velocity width distribution rather than the median. b-parameter 35.3$${\rm km\,s}^{-1}.

In comparison to earlier work, Davé et al. (2010) also reported a small difference of \sim 1$${\rm km\,s}^{-1} in the median $b_{\rm HI}$ values predicted by their constant wind and momentum-driven wind models at $z=0$. However, their velocity widths are generally higher then we find here, with a median value of \sim 43$${\rm km\,s}^{-1}. The mass resolution of their simulations is around factor of two higher than this work ($M_{\rm gas}=2.49\times 10^{7}h^{-1}M_{\odot}$) although performed in smaller volumes $\,(48h^{-1}\rm\,cMpc)^{3}$. This difference is therefore most likely because Davé et al. (2010) do not deconvolve their Voigt profile fits with the instrument profile (see the first paragraph in section 4 of that paper). Tepper-García et al. (2012) similarly find that AGN feedback produces a small increase in the typical Ly$\rm{\alpha}$ line widths at $z=0.25$ using the OWLS simulations, and report a median b-parameter of 29.4\,$${\rm km\,s}^{-1} from mock spectra with FWHM=7$${\rm km\,s}^{-1} and a S/N$=50$ per pixel. This is similar to the value of 27.7$${\rm km\,s}^{-1} we obtain from our AGN feedback (80-512-ps13+agn) model using the same signal-to-noise and instrument resolution. As Tepper-García et al. (2012) also use SPH simulations performed at slightly lower mass resolution compared to this work ($M_{\rm gas}=8.7\times 10^{7}h^{-1}M{\odot}$), their line widths will be similarly under-resolved. More recently, Gaikwad et al. (2017b) also found their simulated velocity width distribution predicts lines that are narrower than the COS data. It is possible that the lack of any coupling between radiative cooling, photo-ionisation and the hydrodynamical response of the gas in their simulations may exaggerate this difference further, however.

The impact of feedback on the velocity distribution at $z=[0.1,\,1.0,\,1.6]$ is displayed in Figure 6, for FWHM$=7\rm\,km\,s^{-1}$ and $\rm S/N=50$ per pixel. Again it is clear that the velocity width distribution is not converged with mass resolution. By comparing the 80-1024 (orange dot dashed line) and 80-512 (blue dashed line) the line widths differ significantly at $b_{\rm HI}<30\rm\,km\,s^{-1}$. This implies a comparison between simulations and observations requires a mass resolution of at least $m_{\rm gas}\simeq 6\times 10^{6}h^{-1}M_{\odot}$. However, as for the CDDF, the simulations with varying sub-grid physics performed at the same mass resolution are very similar. The only exception is at $z=0.1$ when hot gas from AGN feedback produces slightly broader lines.

3.3 The evolution of number density of absorbers

Lastly we examine the number of H\,\rm\scriptstyle I\absorbers per unit redshift, $dN/dz$. Figure 7 displays $dN/dz$ for low ($10^{13.1}\leq N_{\rm HI}/\rm{cm}^{-2}\leq 10^{14}$, upper panel) and high ($10^{14}<N_{\rm HI}/\rm{cm}^{-2}\leq 10^{17}$, lower panel) column density absorbers, overlaid with a compilation333The individual measurements are from Hu et al. (1995); Lu et al. (1996); Kim et al. (1997); Kirkman & Tytler (1997); Savaglio et al. (1999); Kim et al. (2001); Sembach et al. (2004); Williger et al. (2006); Aracil et al. (2006); Janknecht et al. (2006); Lehner et al. (2007); Williger et al. (2010). of observational data from Kim et al. (2013). In this instance we do not attempt to match the resolution and signal-to-noise of the simulations to the data compilation, which was obtained with a variety of instruments. We instead use FWHM=7\,$${\rm km\,s}^{-1} and S/N$=50$ per pixel and assume the HM01 UVB for emission from star forming galaxies and quasars. We also show the more recent COS measurements444Note the Danforth et al. (2016) measurements in the upper panel of Fig. 7 are for the slightly different column density range $10^{13}\leq N_{\rm HI}/\rm{cm}^{-2}\leq 10^{14}$, rather than $10^{13.1}\leq N_{\rm HI}/\rm{cm}^{-2}\leq 10^{14}$ used for the Kim et al. (2013) data compilation. from Danforth et al. (2016) at $0\leq z\leq 0.47$. These are formally consistent with the earlier data compilation at $z\leq 0.3$, but with an increased incidence of lines at $0.3<z<0.5$. Note also we use Ly$\rm{\alpha}$ only fits for comparison to the observations, which may complicate the comparison to the stronger, saturated lines with column densities inferred using higher order Lyman lines. We have therefore tested this by also directly integrating the H\,\rm\scriptstyle I\number densities along the simulated sight-lines to obtain $dN/dz$ (see also Gurvich et al., 2017). This procedure yields similar results to the Voigt profile fits.

The number of low column density Ly$\rm{\alpha}$ forest absorbers (top panel) is not sensitive to the differences in star formation and feedback prescriptions in the models. This is consistent with the CDDF, where galactic winds and AGN feedback do not significantly impact on the gas distribution at the low column densities that probe the diffuse IGM (see also Davé et al., 2010). The simulations are in good agreement with the observational data at $z>1$, but underpredict the number of weak lines at $z=0.1$. To investigate this further, we have rescaled the mean transmission in the 80-512-ps13 model to $\langle F\rangle=0.979$ at $z=0.1$ (blue triangle), which V17 find provides a better match to CDDF at $10^{13}<N_{\rm HI}/\rm{cm}^{-2}<10^{14}$. The improved agreement again demonstrates this discrepancy is due to the uncertain UVB amplitude; the HM01 model overpredicts the H\,\rm\scriptstyle I\photo-ionisation rate at $z\la 0.5$.

In contrast to the low column density systems, the simulated predictions for higher column density Ly$\rm{\alpha}$ forest absorbers exhibit some differences with the choice of star formation and feedback model at $z<1.5$. The 80-512-ps13 model produces a factor of $\sim 2$ higher $dN/dz$ relative to 80-512 at $z\simeq 0$, with the AGN feedback intermediate between the two. However, all models fail to reproduce the observational data at $z\leq 1.5$. Rescaling the mean transmission at $z=0.1$ slightly improves agreement but does not remove the discrepancy. In combination with the CDDF, this strongly suggests the simulations are not correctly capturing the saturated Ly$\rm{\alpha}$ absorption systems that are still optically thin to Lyman continuum radiation.

4 Identifying the gas responsible for low redshift Ly$\rm{\alpha}$ forest absorbers

In order to more closely examine the physical origin of the absorption in our models and assess the effect of feedback, we consider the typical temperature and density of gas in the QUICKLYA model (80-512) and the Puchwein & Springel (2013) energy driven winds model with (80-512-ps13+agn) and without (80-512-ps13) AGN feedback. We have divided the gas into four categories based on temperature, $T_{\rm th}=10^{5}\rm\,K$, and normalised density thresholds, $\Delta_{\rm th}$, following Davé et al. (2010). These are i) diffuse gas ($\Delta<\Delta_{\rm th},T<10^{5}$K) ii) the warm-hot IGM (WHIM, $\Delta<\Delta_{\rm th},T>10^{5}$K) iii) hot halo gas ($\Delta>\Delta_{\rm th},T>10^{5}$K, and iv) condensed gas ($\Delta>\Delta_{\rm th},T<10^{5}$K). In Table 2, the percentage of the total baryon mass in each category is shown at $z=0.1,\,1$ and $1.6$ where $\Delta_{\rm th}=97.2,\,65.9$ and $62$, respectively (here $\Delta_{\rm th}$ is computed using equations 1 and 2 in Davé et al., 2010). We have also included the percentage of gas that has been converted into collsionless star particles.

The redshift evolution of gas in the diffuse phase is very similar for all models, with no noticeable difference due to feedback. This diffuse gas follows a power law temperature-density relationship, $T=T_{0}\Delta^{\gamma-1}$ (see Figure 1). However, there is a significant increase in the mass fraction in the WHIM for the model with AGN feedback at $z=0.1$, with a corresponding decrease in the condensed phase fraction and stellar mass. Including multi-phase star formation and stellar feedback produces more gas in the hot halo and condensed phases compared to the QUICKLYA model at $z=0.1$. In contrast, the QUICKLYA model produces significantly more stars at $z\leq 1$ and – while capturing the diffuse IGM reasonably well – will produce unreliable results for gravitationally bound gas at $\Delta>\Delta_{\rm th}$.

The mass fractions for the 80-512-ps13 model are broadly similar to the “vzw” momentum driven winds model described in Davé et al. (2010), with slightly more ($\sim 5$ per cent) gas in the dense phases with $\Delta>\Delta_{\rm th}$. There is also reasonable agreement between the 80-512-ps13 model and the post-processed adiabatic simulations described by Gaikwad et al. (2017a) at $z=0$, with around 4 per cent more (less) gas in the WHIM (hot halo) phase. This suggests the bulk of the hot gas with $T>10^{5}\rm\,K$ in the models without AGN feedback is produced by gravitational infall. The largest difference is the condensed phase which contains 18.8 per cent of the mass at $z=0$ in Gaikwad et al. (2017a); this is because dense gas is not converted into stars in their adiabatic simulations.

The connection between the density and temperature of the gas and the column density and velocity widths of the Ly$\rm{\alpha}$ absorbers at $z=0.1$ is summarised in Figure 8. Along each simulated line of sight we have computed the optical depth weighted temperature, $T_{\tau}$, and density, $\Delta_{\tau}$, where $X_{\tau,\rm j}=\sum_{\rm i}\tau_{\rm i,j}X_{\rm i}/\sum_{\rm i}\tau_{\rm i,j}$ and $\tau_{\rm i,j}$ is optical depth at the $i^{\rm th}$ pixel in real space contributing to the $j^{\rm th}$ pixel in velocity space (Schaye et al., 1999). We then associate $T_{\tau}$ and $\Delta_{\tau}$ with $N_{\rm HI}\,$ and $b_{\rm HI}\,$ at the location of the line centres identified by VPFIT. This approach takes redshift space distortions due to the peculiar motion of the gas into account. As before, following V17 we only consider systems with relative errors less then 50 per cent on $b_{\rm HI}$ and $N_{\rm HI}$.

The top panel of Figure 8 displays the relationship between $N_{\rm HI}$ and $\Delta_{\rm\tau}$ in the simulations. The data are coloured according to their velocity widths. The red dashed line shows the power-law relationship between $N_{\rm HI}\,$and $\Delta$ predicted by Schaye (2001), assuming the typical size of an H\,\rm\scriptstyle I\absorber is the local Jeans scale:

[TABLE]

Here $\Gamma_{-13}=\Gamma_{\rm HI}/10^{-13}\rm\,s^{-1}$ is the H\,\rm\scriptstyle I\photo-ionisation rate (where $\Gamma_{-13}=0.78$ at $z=0.1$ in the 80-512-ps13 simulation for $\langle F\rangle=0.979$) and $T_{4}=T/10^{4}\rm\,K$ is the gas temperature.

This provides a good description of the scaling between gas density and H\,\rm\scriptstyle I\column density (see also Davé et al., 1999, 2010; Tepper-García et al., 2012). The majority of the absorption lines at $N_{\rm HI}=10^{13}$–$10^{14}\rm\,cm^{-2}$ correspond to gas densities $\Delta=4$–$40$ in the simulations (V17), with typical velocity widths of $b_{\rm HI}=20$–$40\rm\,km\,s^{-1}$. In all models there are a small number of broad lines, $b_{\rm HI}>40\rm\,km\,s^{-1}$, with column densities $N_{\rm HI}=10^{12.5}$–$10^{13.5}\rm\,cm^{-2}$ that are associated with hot $T\sim 10^{5}\rm\,K$ gas in the WHIM with $\Delta\simeq 1$–$100$ (see also Tepper-García et al., 2012). The effect of feedback on the diffuse, low column density gas is rather subtle; there are slightly more broad lines lying above the $N_{\rm HI}-\Delta_{\rm\tau}$ correlation for the AGN feedback model, consistent with additional heating. Although there is significantly more hot gas in the WHIM by mass at $z=0.1$ when including AGN feedback, much of this is also collisionally ionised. In all instances, however, the higher column density systems, $N_{\rm HI}>10^{14}\rm\,cm^{-2}$ correspond to gas densities at $\Delta\ga 20$. The simulations fail to correctly capture the number density of these systems, some of which are at densities associated with circumgalactic gas.

The lower panel in Figure 8 shows the corresponding values of $T_{\rm\tau}$ at the line centres against the velocity widths. The red dashed line displays the expected velocity width under the assumption of pure thermal broadening:

[TABLE]

In general, although the minimum line widths trace the expectation for thermal broadening, most of the absorbers are significantly broader, particularly for low column density lines originating from low density gas (see also Garzilli et al., 2015, for a recent discussion at $z\sim 3$). Many of the absorbers are therefore broadened through the additional effects of gas pressure (i.e. Jeans smoothing). The absorbers with $N_{\rm HI}=10^{13}$–$10^{14}\rm\,cm^{-2}$ are associated with gas with a broad range of temperatures, but most are in the range $T=10^{4}$–$10^{4.5}\rm\,K$. The velocity widths of these lines remain too narrow with respect to the COS data. Again, the effect of stellar and AGN feedback on the temperature of the absorbers is small.

5 Conclusions

We have investigated the impact of different feedback prescriptions on the properties of the low redshift Ly$\rm{\alpha}$ forest, using a selection of hydrodynamical simulations drawn from the Sherwood simulation suite. The simulations adopt stellar feedback, AGN feedback and a simplified scheme for modelling the low column density Ly$\rm{\alpha}$ forest that converts all gas with $\Delta>1000$ and $T<10^{5}\rm\,K$ into collisionless star particles. We have examined the CDDF, velocity width distribution, the redshift evolution in the line number density, $dN/dz$, and the relationship between the column density, velocity widths, and the underlying gas density and temperature.

Our analysis confirms the finding by V17 that the velocity widths of simulated Ly$\rm{\alpha}$ absorption lines with gas densities $\Delta=4$–$40$ are too narrow when compared with COS data. We also stress that this difference will be exacerbated further for simulations performed at higher mass resolution than considered here; the line widths are under-resolved at our fiducial resolution of $M_{\rm gas}=5.1\times 10^{7}h^{-1}M_{\odot}$ (see also the Appendix). We also find that absorbers with column densities $N_{\rm HI}>10^{14}\rm\,cm^{-2}$, which correspond to gas with $\Delta\ga 20$ in our simulations, are underestimated, in agreement with earlier work (Shull et al., 2015; Gurvich et al., 2017).

On the other hand, using simulations that vary the feedback prescription while keeping other parameters fixed, we find that the impact of galactic winds and AGN feedback is generally very modest. AGN feedback only impacts on the Ly$\rm{\alpha}$ absorbers at very low redshift, $z=0.1$, by producing some broader lines and slightly reducing the incidence of systems with $N_{\rm HI}>10^{14}\rm cm^{-2}$ with respect to a model with only stellar feedback. There is no noticeable impact on the CDDF due to stellar or AGN feedback at $N_{\rm HI}<10^{14}\rm\,cm^{-2}$. Furthermore, we find a simplified scheme that ignores feedback and star formation altogether adequately capture the properties of absorption lines with $N_{\rm HI}<10^{14}\rm\,cm^{-2}$, although the incidence of stronger systems is underestimated compared to models that include feedback. Similarly, the Ly$\rm{\alpha}$ line number density and its evolution with redshift, $dN/dz$, is in good agreement with observational data for low column density systems, $10^{13.1}\leq N_{\rm HI}/\rm cm^{-2}\leq 10^{14}$, once the amplitude of the UVB is rescaled to match the CDDF. In contrast, for high column density systems, $10^{14}\leq N_{\rm HI}/\rm cm^{-2}\leq 10^{17}$, the simulations underpredict the incidence of absorbers at $z=0.1$ and $z=1$ by a factor of $\simeq 2.2$ and $3$, respectively.

This demonstrates that stellar and AGN feedback – as currently implemented in our simulations – have a limited impact on the Ly$\rm{\alpha}$ forest. We conclude that resolving the discrepancy between COS data and the simulations requires an increase in the temperature of gas with $\Delta=4$–$40$, either through (i) additional He\,\rm\scriptstyle II\photo-heating at $z>2$, (ii) additional fine-tuned feedback that ejects overdense gas into the IGM at just the right temperature for it to still contribute significantly to the Ly$\rm{\alpha}$ forest or alternatively (iii) a larger, currently unresolved turbulent component to the line widths (Oppenheimer & Davé, 2009; Iapichino et al., 2013; Gaikwad et al., 2017b). Note, however, the first possibility would require a UVB with a harder spectrum at $z>2$, potentially producing gas temperatures that are too high with respect to existing observational constraints from the Ly$\rm{\alpha}$ forest at $1.6\leq z\leq 3$ (Becker et al., 2011; Boera et al., 2014; Bolton et al., 2014). Our results provide further motivation for studying the impact of state-of-the-art feedback models on the properties of the low redshift Ly$\rm{\alpha}$ forest.

Acknowledgments

The hydrodynamical simulations used in this work were performed with supercomputer time awarded by the Partnership for Advanced Computing in Europe (PRACE) 8th Call. We acknowledge PRACE for awarding us access to the Curie supercomputer, based in France at the Tre Grand Centre de Calcul (TGCC). This work also made use of the DiRAC High Performance Computing System (HPCS) and the COSMOS shared memory service at the University of Cambridge. These are operated on behalf of the STFC DiRAC HPC facility. This equipment is funded by BIS National E-infrastructure capital grant ST/J005673/1 and STFC grants ST/H008586/1, ST/K00333X/1. We thank Volker Springel for making P-GADGET-3 available. FN is supported by a Vice-Chancellor’s Scholarship for Research Excellence. JSB acknowledges the support of a Royal Society University Research Fellowship. MV and TSK are supported by the FP7 ERC grant “cosmoIGM” and the INFN/PD51 grant. TSK also acknowledges the NSF-AST118913. MGH and EP acknowledge support from the FP7 ERC Grant Emergence-320596, and EP gratefully acknowledges support by the Kavli Foundation. DS acknowledges support by the STFC and the ERC starting grant 638707 “Black holes and their host galaxies: co-evolution across cosmic time”.

Appendix A numerical convergence

We have performed convergence tests with mass resolution and box size for the Ly$\rm{\alpha}$ forest CDDF and velocity width distribution at the redshifts considered in this paper. These are displayed in Figures 9-12. The CDDF is well converged with mass resolution and box size within the range $10^{12.7}\leq\ N_{\rm HI}/\rm{cm}^{-2}\leq 10^{14}$. However, at $N_{\rm HI}>10^{14}\rm\,cm^{-2}$ ($N_{\rm HI}<10^{12.7}\rm\,cm^{-2}$) the number of lines are over (under) estimated at lower mass resolution. This suggests that numerical resolution is unlikely to be the cause of the discrepancy between the COS observations and our simulations at $N_{\rm HI}>10^{14}\rm\,cm^{-2}$.

In Figure 11, however, it is clear that mass resolution significantly impacts on the velocity width distributions, with line widths that are systematically overestimated at our fiducial resolution of $M_{\rm gas}=5.1\times 10^{7}h^{-1}M_{\odot}$ (80-512). A mass resolution of at least $M_{\rm gas}=6.4\times 10^{6}h^{-1}M_{\odot}$ is required for an acceptable level of convergence. This means the discrepancy between the COS observations of the velocity width distribution will be even larger for higher resolution simulations. The Ly$\rm{\alpha}$ line width distribution is generally well converged with box size, as shown Figure 12, with the exception of lines at b_{\rm HI}<20$${\rm km\,s}^{-1}in the $20h^{-1}\rm\,cMpc$ box.

Finally, Figures 11-12 and Table 3 display convergence tests for the volume weighted IGM thermal state with box size and mass resolution. The results are well converged with box size, but demonstrate the temperature of the IGM at mean density is overestimated by $10$–$20$ per cent at our fiducial resolution (80-512). We have checked, however, that the mass weighted temperature obtained from the gas particles is well converged with mass resolution, suggesting this difference originates from the larger smoothing lengths used in the lower resolution models. Note also that $T(\bar{\Delta})$, shown in the middle panels of Figures 11-12, corresponds to $T(\bar{\Delta})=T_{0}\bar{\Delta}^{\gamma-1}$, and therefore reflects the convergence of both $T_{0}$ and $\gamma-1$. Additionally, we observe from Table 3 that as mass resolution is increased using the QUICKLYA model, a greater proportion by mass of the baryons are converted to collisionless star particles. Note, however, this is already a significant overestimate and the results of the QUICKLYA model will be unreliable for gas at $\Delta>\Delta_{\rm th}$.