# Hydrodynamical Simulations of the Lyman Alpha Forest: Model Comparisons

## 1 Introduction

Several numerical simulations of the Ly$`\alpha `$forest in cold dark matter (CDM) dominated cosmologies have been performed in recent years and compared with observations (Cen et al. 1994; Zhang, Anninos, & Norman 1995; Hernquist et al. 1996; Miralda–Escudé et al. 1996; Zhang et al. 1996; Davé et al. 1997; Bond & Wadsley 1997; Zhang et al. 1998). Remarkably all the simulations have been able to reproduce the measured neutral hydrogen column density distribution, the size of the absorbers (Charlton et al. 1997), and the line number evolution reasonably well, despite the differences in the cosmological models used: Cen et al. adopt a $`\mathrm{\Lambda }`$CDM model; Zhang et al. investigate sCDM models with both an unbiased and a cluster scale normalization; Hernquist et al. (1996) evolve an sCDM model with a cluster scale normalization. The distribution of Doppler parameters has fared somewhat less well: the predicted distribution peaks toward lower values than observed when the simulations are performed with adequate resolution (Bryan et al. 1998; Theuns et al. 1998). Nonetheless, the generally good agreement with observations of the Ly$`\alpha `$forest suggests that the models are capturing the essential physical properties of the absorbers. This has prompted recent work by Croft et al. (1998) aimed at using flux statistics of the observational data to extract the fluctuation spectrum of the underlying cosmology. We are thus encouraged to investigate the possibility that differences in the statistical properties of the Ly$`\alpha `$forest predicted by different cosmological models may provide a means of testing the models.

The objective of this paper is to compare the Ly$`\alpha `$forest statistics derived from simulations in different cosmological models and to investigate what key properties of the cosmological models control a given statistic. We present results from nine numerical simulations using five different background cosmological models, three of which are flat with no cosmological constant, one is open, and one is flat with a nonzero cosmological constant. For five of the simulations, which we will refer to as the model comparison study, the parameters of the cosmological models have been selected by their ability to match the local or low redshift observations, although all of these models except the standard cold dark matter (sCDM) model are also consistent with COBE measurements of the cosmic microwave background. A tilted CDM model is further designed to match COBE constraints on the normalization of the power spectrum on large scales. In the remaining four simulations we keep the underlying cosmology fixed (sCDM) while varying the normalization of the fluctuation power spectrum in order to clarify the dependence of the Ly$`\alpha `$statistics on this parameter. The radiation field is normalized to the absorption properties of the Ly$`\alpha `$forest as measured at high redshifts. While the emphasis in this paper is on comparing cosmological models, we also test how well the models are doing using the tabulated statistics of the Ly$`\alpha `$forest as determined primarily by Kim et al. (1997) for several QSO lines-of-sight. A more complete comparison with existing data from several observational groups will be presented in Meiksin et al. (1999).

The paper is organized in the following way. In §2 we describe the cosmological models and simulation technique. In §3 we investigate the model differences, power dependence, and redshift evolution in the raw opacity data as characterized by nonparametric statistics of the flux and optical depth. In §4 we present a line analysis of the spectra generated by the various simulations focusing on the column density distribution and line number evolution statistics. In §5 we discuss the Doppler b parameter distributions and related nonparametric statistics, and in §6 we present model predictions for He $`\mathrm{II}`$absorption. We summarize our results in §7.

## 2 The Models and Simulations

All the model background spacetimes we consider are in the context of Cold Dark Matter (CDM) dominated cosmologies. We examine the following five models: a standard critical density flat CDM model (sCDM), a flat CDM model with a nonvanishing cosmological constant ($`\mathrm{\Lambda }`$CDM ), a topologically open CDM model (OCDM), the standard CDM model but with the power spectrum of the density perturbations tilted (tCDM) to match the normalization on large scales as determined from the COBE measurements of the Cosmic Microwave Background (Bunn & White 1997), and a flat critical density mixed dark matter model with a hot component added to the CDM (CHDM). There are several important and well–established astrophysical measurements which constrain the various combinations of cosmological parameters. The parameters for each model, which we list in Table 1, have been determined to provide good consistency with these observations. For example, the combination $`\mathrm{\Omega }_Bh^2`$ is restricted by Big Bang nucleosynthesis constraints and the measured abundance of primordial deuterium to lie in the range $`0.015`$$`0.025`$ (Copi, Schramm, & Turner 1995; Burles & Tytler 1998). In addition, because the $`\mathrm{I}`$column density scales approximately as $`(\mathrm{\Omega }_bh^2)^2`$ for a fixed UV radiation intensity, we choose $`\mathrm{\Omega }_b`$ and $`h`$ so that $`\mathrm{\Omega }_bh^2`$ is the same for three of the models sCDM, $`\mathrm{\Lambda }`$CDM and OCDM. The fluctuation normalization in a sphere of $`8h^1`$ Mpc is defined to match observations of the number density of galaxy clusters (White, Efstathiou, & Frenk 1993; Bond & Myers 1996) in all the models. In addition, a tilt has been applied to the CDM power spectrum in the tCDM model in order to approximately match the amplitude of the CMB quadrupole as measured by COBE (Bunn & White 1997). The cosmological constant in the $`\mathrm{\Lambda }`$CDM case is consistent with the upper limit ($`\mathrm{\Omega }_\mathrm{\Lambda }<0.7`$) of Maoz and Rix (1993) and the best fit parameters of Ostriker and Steinhardt (1995). One of the major problems with the sCDM model is its difficulty in matching observations of the large scale structures in the universe. Since the standard CDM model is historically one of the most studied models, however, we use it as our canonical model to which the perhaps more viable additional models considered here may be compared and through which we investigate the dependence of the Ly$`\alpha `$statistics on the fluctuation power spectrum. We refer the reader to Zhang et al. (1995, 1997, 1998) for further details and results from our previous sCDM simulations.

The initial data were generated using COSMICS (Bertschinger 1995) with the BBKS transfer function (Bardeen et al. 1986) to compute the starting redshifts and the initial particle positions and velocity perturbations appropriate for all models except CHDM. We used CMBFAST (Seljak & Zaldarriaga 1996) to solve the linearized Boltzmann equations to set the initial conditions for CHDM. For the comoving box size adopted (9.6 Mpc) and the corresponding comoving grid cell size (37.5 kpc in our high resolution runs), the relevant wavenumber domain of the simulations at $`z=0`$ is $`168>k>0.65`$ Mpc<sup>-1</sup>, where $`k=2\pi /\mathrm{}`$ and $`\mathrm{}`$ is the length scale. Over this domain, the sCDM, $`\mathrm{\Lambda }`$CDM and OCDM models all have a similar power distribution. The tCDM and CHDM models, on the other hand, have an overall lower normalization (see Table 1) in addition to a steeper slope that drops slightly more sharply than the other models over the smaller scales. In Figure 1 we show the linear power spectra for these models evolved to $`z=3`$, the redshift at which we present many of our results. Since previous work (Zhang et al. 1998) indicates that the sizes of the low column density absorbers at $`z3`$ are $`100`$ kpc, it is useful to characterize the models in terms of their power at these small fluctuation scales. A useful measure of this power, introduced by Gnedin (1998), is

$$\sigma _{34}^2=_0^{\mathrm{}}P(k)e^{2k^2/k_{34}^2}\frac{k^2dk}{2\pi ^2}$$
(1)   
(where $`k=2\pi /\mathrm{}`$, $`P(k)`$ is the linear power spectrum at $`z=3`$ and $`k_{34}=34\mathrm{\Omega }_0^{1/2}h`$ Mpc<sup>-1</sup>). This is also listed for each model in Table 1.

In addition to specifying a cosmological model, it is also necessary to include a background UV radiation field to ionize the IGM. We have implemented the spectrum computed by Haardt & Madau (1996) for a flat universe on the basis of radiation transfer in a clumpy universe and the measured luminosity function of QSOs, accounting for QSO sources, absorption by the Ly$`\alpha `$forest and Lyman limit systems, and re-emission of the recombination radiation from the absorbing clouds. This spectrum reionizes the universe between redshifts $`7`$ and $`6`$ and peaks at about $`z=2`$. Because the OCDM model corresponds to a cosmology not considered by Haardt & Madau (1996), we use the field for a flat universe for this model also, noting that the neutral fractions are in any case rescaled, as described below, to match observations. We also note that only clouds which are optically thin at the Lyman edge are considered in this paper. Hence the optically thin limit is a good approximation and it is not necessary to account for self–shielding and radiative transfer of the external ionizing radiation field.

The numerical computations were performed using two different numerical codes, Kronos and Hercules, each in a simulation box of length $`9.6\mathrm{Mpc}`$ comoving with the universal expansion. Kronos (Bryan et al. 1995) is a single grid Eulerian code that uses a particle-mesh algorithm to follow the dark matter and the piecewise parabolic method (PPM) to simulate the gas dynamics. Since non–equilibrium chemistry and cooling processes can be important, six particle species ($`\mathrm{I}`$, $`\mathrm{II}`$, He $`\mathrm{I}`$, He $`\mathrm{II}`$, He $`\mathrm{III}`$and the electron density) are followed with a sub-stepped backward finite-difference technique (Abel et al. 1997; Anninos et al. 1997). This is the same non–equilibrium chemistry and cooling model used in our previous studies of the Ly$`\alpha `$Forest (Zhang et al. 1995, 1997, 1998; Charlton et al. 1997; Bryan et al. 1998). For sCDM, $`\mathrm{\Lambda }`$CDM , OCDM, and tCDM we use $`256^3`$ grid cells in the simulation box to follow the evolution of $`128^3`$ dark matter particles. This results in a spatial resolution of $`\mathrm{\Delta }x=37.5\mathrm{kpc}`$. For CHDM $`128^3`$ grid cells are used with $`64^3(128^3)`$ cold (hot) dark matter particles, respectively, resulting in a lower spatial resolution, $`\mathrm{\Delta }x=75\mathrm{kpc}`$, for this model. For the sCDM simulations with varying $`\sigma _{8h^1}`$ we present simulations with both $`128^3`$ and $`256^3`$ grid cells resulting in both low and high spatial resolutions of $`75\mathrm{kpc}`$ and $`37.5\mathrm{kpc}`$, respectively.

We have also simulated three of our models (sCDM, $`\mathrm{\Lambda }`$CDM , OCDM) with a different numerical code, Hercules (Anninos, Norman & Clarke 1994; Anninos et al. 1997). Hercules is a nested grid code that utilizes a multiscale PM method for the dark matter, artificial viscosity methods for the baryonic fluid, and the same non–equilibrium chemistry and cooling model as above. The simulations produced from this code use $`128^3`$ particles and $`128^3`$ cells for both the nested and parent grids. However, in order to derive a more representative sample for statistics, the results discussed in this paper are extracted from the parent grid only. Thus these simulations are of lower spatial resolution than most of the Kronos simulations, although the dark matter mass resolution is the same. For statistics that are insensitive to spatial resolution, a comparison of the results of the two codes is useful to insure that simulation results are robust against changes in numerical technique.

Synthetic spectra are generated along $`300(900)`$ random lines of sight through the Kronos (Hercules) simulated volume using the method of Zhang et al. (1997) including the effects of peculiar velocity and thermal broadening of the gas. (We have verified that decreasing the sample size from 900 to 300 for the Kronos data does not affect the results except for a slight increase in the scatter of the line properties for the highest, optically thick column density systems, a regime where our results become unreliable anyway because of the absence of radiative transfer in the codes.) Since we are primarily concerned in this paper with a comparison of model predictions, we have not included noise or continuum fitting in the analysis. Furthermore the resolution of the spectra, $`1.2\mathrm{km}/\mathrm{s}`$, is the same for all the simulations, a value that is smaller than current observations. However, we have shown elsewhere (Bryan et al. 1998) that, as long as we restrict ourselves to high quality observational data, the impact of not including these observational difficulties is small. In addition to analyzing the raw optical depth and flux distributions, line lists are extracted from the data using a Voigt profile fitting procedure. This is described in more detail elsewhere (Zhang et al. 1997), but we outline it briefly here. First, maxima in the optical depth distribution are identified as line centers. Then Voigt profiles are fit, using a non-linear minimization, to the part of the spectrum which is above $`\tau _{HI}=0.05`$ and between neighbouring minima. This results in the same spectral threshold $`F_te^{\tau _t}=0.95`$ as the high resolution Keck HIRES spectrometer. Each line of sight chosen produces a sample spectrum with on the order of $`10`$$`100`$ lines per redshift interval $`\delta z=0.1`$ depending on the redshift and cosmological model. The statistics of these linelists are discussed in §4 and §5.

The amplitudes of the distributions found in the models cannot be used as a basis for comparing the models since they may be arbitrarily re-scaled for any individual model using the ionization bias factor $`b_{\mathrm{ion}}=\mathrm{\Omega }_B^2/\mathrm{\Gamma }`$, where $`\mathrm{\Omega }_B`$ is the fraction of the critical density carried in baryons and $`\mathrm{\Gamma }`$ is the Haardt-Madau (1996) parameterization of the metagalactic UV ionizing background extracted from the observed distribution of quasars. It is important to normalize all the models consistently before comparing the shapes of any of the distributions. This may be done in a variety of ways. We do so by matching the mean $`\mathrm{I}`$opacity in each simulation to the measured intergalactic $`\mathrm{I}`$opacity at $`z=3`$. In Zhang et al. (1997), we found the opacity measurements of Steidel & Sargent (1987) and Zuo & Lin (1993) gave a mean $`\mathrm{I}`$opacity at $`z=3`$ of $`\overline{\tau }_\alpha =0.270.35`$, although values as much as 30–60% larger have been claimed (Press et al. 1993; Rauch et al. 1997). Because of the uncertainty in this measurement, we also require consistency with the number density of lines observed above a threshold of $`\mathrm{log}N_{\mathrm{HI}}=13.5`$. Using the three quasars in Hu et al. (1995) for which lines in the full redshift range $`3<z<3.1`$ are listed, we find a total of 61 lines for the three lines–of–sight in this redshift interval with $`\mathrm{log}N_{\mathrm{HI}}>13.5`$, for which the line lists should be complete. (An estimate based on using the available lines for all four QSO line lists in Hu et al. in the redshift interval $`2.9<z<3.1`$ gave essentially the same line density.) Normalized to $`\overline{\tau }_\alpha =0.30`$, the CHDM, sCDM, $`\mathrm{\Lambda }`$CDM , OCDM, and tCDM models predict, respectively, 60.8, 62.1, 62.7, 63.7, and 59.5 lines, in close agreement with the observed number. Normalizing to $`\overline{\tau }_\alpha =0.35`$, the respective numbers of predicted lines are 73.7, 74.3, 73.9, 75.2, and 72.8. While these are not badly inconsistent with the observed number, they are all fairly high. We normalize the spectra according to $`\overline{\tau }_\alpha =0.30`$ throughout this paper, noting that this value is still not well agreed upon. In Figure 2 we plot a related statistic, $`\tau _{eff}`$ (Zhang et al. 1997) for the normalized spectra of our models and compare to recent data by Kirkman & Tytler (1997). After normalization all of our models are consistent with the data over the redshift range $`2z4`$ considered by this paper.

## 3 Direct Optical Depth and Flux Measurements

Historically Ly$`\alpha `$absorption spectra have been analyzed in terms of the statistics of spectral line features and as such have been plagued with difficulties of the line fitting procedure such as line identification and blending. Many of these difficulties become increasingly severe at higher redshifts making the results of the analysis uncertain. It is thus natural to ask whether statistics dependent directly on the observed flux and optical depth without recourse to line fitting might be of use in describing the forest and discriminating among competing models. Statistics of this kind have recently been proposed by several authors (Miralda–Escudé et al. 1997; Rauch et al. 1997; Cen 1997). Since these nonparametric measures are also easier to relate theoretically to the physical state of the absorbing gas, we begin our discussion with them.

### 3.1 Optical Depth Probability Distribution Function

The optical depth $`\tau `$ is related to the transmitted flux $`F`$ by $`F=exp(\tau )`$. We define the optical depth probability distribution $`dP/d\tau `$ as the probability that a pixel will have optical depth between $`\tau `$ and $`\tau +d\tau `$. In Figure 3 we use spectra generated from the sCDM high resolution simulation to show $`\tau dP/d\tau `$ versus $`\tau `$ for redshifts $`z=2,3`$ and $`4`$ (top panel). Although the peak of the distribution decreases and the distribution broadens slightly with decreasing redshift, the principal contributor to the redshift evolution seen in Figure 3 is the evolution of the optical depth $`\tau `$. Hui, Gnedin & Zhang (1997) discuss in detail the dependence of $`\tau `$ on the distribution and properties of neutral hydrogen along the line of sight in an expanding universe. Since we would like to understand the redshift evolution of the optical depth in terms of simple scaling laws, we repeat some of their discussion here in order to the isolate the key factors controlling this redshift evolution and clarify the scaling law assumptions. The optical depth is defined as

$$\tau (\nu _o)=_{x_a}^{x_b}n_{HI}\sigma _\alpha \frac{dx}{1+z}$$
(2)   
where $`\nu _o`$ is the observed frequency, $`n_{HI}`$ is the number density of neutral hydrogen, $`z`$ is the redshift of the absorbing gas, $`\sigma _\alpha `$ is the absorption cross section for Ly$`\alpha `$, and the integral is over the line of sight between the quasar ($`x_a`$) and the observer ($`x_b`$) in comoving coordinates. In practice the form of the Ly$`\alpha `$absorption cross section limits the integration range per absorber to a small portion of the line of sight. It is thus useful to make a change of variable to velocity coordinates $`u`$ about some characteristic average redshift $`\overline{z}`$ in the problem. For example, for simulated data the redshift $`\overline{z}`$ might be a given output redshift for the simulation. The observed frequency $`\nu _o`$ and the frequency $`\nu `$ of the radiation in the absorber rest frame are then related by

$$\nu =\nu _o(1+\overline{z})(1+u/c)$$
(3)   
where

$$u\frac{H(\overline{z})(x\overline{x})}{1+\overline{z}}+v_{pec}(x).$$
(4)   
$`\overline{x}`$ is the comoving position along the line of sight whose redshift is exactly $`\overline{z}`$, $`v_{pec}`$ is the physical velocity of the gas, and $`H(\overline{z})`$ is the Hubble parameter defined by

$$H(\overline{z})=H_0\sqrt{\mathrm{\Omega }_m(1+\overline{z})^3+(1\mathrm{\Omega }_m\mathrm{\Omega }_\mathrm{\Lambda })(1+\overline{z})^2+\mathrm{\Omega }_\mathrm{\Lambda }}.$$
(5)   
The first term in Equation 4 represents the contribution of the residual Hubble flow about the mean while the second term is due to the physical bulk flow of the gas. We assume $`u/c<<1`$ and neglect contributions from turbulent flows since they would be unlikely in the low column density regions we are considering. Under this change of variable the Ly$`\alpha `$cross section becomes

$$\sigma _\alpha =\frac{\sigma _{\alpha 0}c}{b\sqrt{\pi }}e^{(uu_0)^2/b^2}$$
(6)   
where $`\sigma _{\alpha 0}=4.5\times 10^{18}`$ cm<sup>2</sup> sets the scale of the absorption cross section in terms of fundamental constants, $`u_0`$ is the velocity $`u`$ for which the frequency $`\nu `$ in the rest frame of the absorbing gas is equal to the Ly$`\alpha `$frequency $`\nu _\alpha `$ and $`b=\sqrt{2k_BT/m_p}`$ is the thermal width. For absorption lines of neutral hydrogen with column densities $`N_{\mathrm{HI}}<10^{17}`$ cm<sup>-2</sup> the thermal profile dominates the cross section so we neglect the contribution of the natural line width to $`\sigma _\alpha `$. The optical depth $`\tau `$ can now be written as

$$\tau =\frac{\sigma _{\alpha 0}c}{\sqrt{\pi }}\underset{streams}{}\frac{n_{HI}}{b(1+\overline{z})}\left|\frac{du}{dx}\right|^1e^{(uu_0)^2/b^2}𝑑u$$
(7)   
The sum over streams represents the possibility that a given velocity $`u`$ corresponds to more than one position $`x`$. Although the integration formally runs over the full line of sight from quasar to observer, the Gaussian form for the cross section effectively limits the $`u`$ integration to a narrow range around $`u_0`$ (thus justifying our replacement of $`z`$ everywhere by $`\overline{z}`$.) To simplify notation we drop the bar letting $`z`$ represent $`\overline{z}`$ in what follows.

We assume that the number density of hydrogen traces the baryon gas density well. (There has been little metal production in these low density regions and there is no interaction that would cause the helium and hydrogen to separate.) Thus the number density of neutral hydrogen is $`n_{HI}=\rho _bX_{HI}`$ where $`X_{HI}`$ is the neutral fraction and $`\rho _b`$ the gas density. In ionization equilibrium (which is well satisfied except for the period of initial reionization) the neutral fraction of hydrogen is $`X_{HI}\rho _bT^{0.7}`$ such that the number density of neutral hydrogen (relevant to Ly$`\alpha `$absorption) scales as

$$n_{HI}(\mathrm{\Omega }_bh^2)^2\mathrm{\Gamma }^1(z)(1+z)^6(1+\delta _b)^2T^{0.7}$$
(8)   
where $`\delta _b`$ is the baryon overdensity. Studies (Hui & Gnedin, 1996; Weinberg et al. , 1996) of the equation of state for the gas find that for unshocked gas at low to moderate baryon overdensities ($`\delta _b5`$) the equation of state is well fit by a power law:

$$T(1+z)^{1.7}(1+\delta _b)^{\gamma 1}.$$
(9)   
Thus

$$n_{HI}(\mathrm{\Omega }_bh^2)^2\mathrm{\Gamma }^1(z)(1+z)^{4.8}(1+\delta _b)^{2.70.7\gamma }.$$
(10)   
For a uniform radiation field and reionization that occurs before $`z=5`$, as is the case in our simulations, $`\gamma 1.4`$. This is in agreement with the value $`\gamma 1.5`$ found by Zhang et al. (1998) for clouds with column densities in the range $`12.5<\mathrm{log}N_{\mathrm{HI}}<14.5`$. Furthermore the assumption that most of the optical depth arises from low column density absorbers, large structures whose overdensities and peculiar velocities are slowly varying compared to the thermal profiles, means that multiple streaming is rare, the sum over streams in Equation 7 can be dropped, and $`\left|\frac{du}{dx}\right|\frac{H}{1+z}`$. We then integrate over the thermal profile to obtain (Croft et al. 1997)

$$\tau \frac{c\sigma _{\alpha 0}(\mathrm{\Omega }_Bh^2)^2}{\mathrm{\Gamma }(z)H}(1+z)^{4.8}(1+\delta _b)^{1.7}.$$
(11)   
Note that in this limit $`\tau `$ need no longer have a thermal profile about its maximum (Hui, Gnedin & Zhang, 1997). If $`\delta _b`$ is evolving slowly over this redshift range, $`\tau `$ should scale as

$$\tau \frac{(1+z)^{4.8}}{\mathrm{\Gamma }(z)H}.$$
(12)   
In the middle panel of Figure 3 we use this simple scaling law to rescale the $`z=2`$ and $`z=4`$ sCDM distributions from the top panel to $`z=3`$, the redshift at which all the models are normalized. We do this in order to test how well the simulations obey this simple scaling relation: if they followed it exactly then all three curves would overlap. Most, but not all, of the redshift evolution of this distribution is accounted for by the scaling of $`\tau `$ given in Equation 12. Since the evolution of the metagalactic UV radiation field $`\mathrm{\Gamma }`$ is relatively slight over this redshift range, we are left with the remarkable conclusion that most of the evolution of the Ly$`\alpha `$forest is a direct consequence of the universal expansion. The direct numerical results of Zhang et al. (1998) support this conclusion. If we include the evolution of the baryon overdensity, as shown in Figure 4 for the sCDM simulation, and shift the overdensity distribution until the peaks overlap, the $`(1+\delta _b)^{1.7}`$ dependence in Equation 11 for the optical depth distributions predicts additional scaling factors of $`1.64`$ ($`0.77`$) at $`z=2`$ ($`z=4`$), respectively, for $`\tau `$ that bring the distributions (shown in the bottom panel of Figure 3) into close agreement. The remaining small differences, the slight broadening of the distribution and a reduction in its peak amplitude with decreasing redshift, most probably reflect the fact that the shape of the baryon overdensity distribution is also evolving slowly with $`z`$.

The top panel of Figure 5, shows $`\tau dP/d\tau `$ for the five sCDM simulations with varying power and spatial resolution. From this we can see that the optical depth PDF at a given redshift is insensitive to the spatial resolution of the simulation. The shape of the distribution, however, is strongly dependent on the amount of small scale power present. Models with less power at these scales produce narrow, sharply peaked distributions. As the power increases, the distribution flattens and broadens. In the lower panel of Figure 5 we show $`\tau dP/d\tau `$ versus $`\tau `$ for the simulations in the model comparison study. These distributions again display a clear dependence on the power spectrum of the model with OCDM (our model with the most small-scale power) producing the broadest distribution, and CHDM and tCDM (models with the least small-scale power) producing the most sharply peaked distributions. Thus this statistic is particularly promising as a model discriminator in that these differences between models are significant in the range $`0.02<\tau <4`$ that should be accessible to observers.

We quantify this relation between the shape of the $`\tau `$ distribution and the amplitude of the power spectrum by fitting a log-normal to the curves:

$$\tau \frac{dP}{d\tau }e^{(\mathrm{ln}\tau \mathrm{ln}\tau _0)^2/2\sigma _\tau ^2}.$$
(13)   
Although this does not fit the profiles in Figure 5 in detail, it does provide an adequate description as long as we restrict the range of optical depths used in the fitting. Here we adopt $`0.02<\tau <4`$, corresponding roughly to the observable range. A different range or a different fitting function changes the details, but not the nature of our result. In Figure 6, we show the correlation between $`\sigma _\tau `$, a measure of the width of the distribution, and $`\sigma _{34}`$, the amplitude of the linear power spectrum on small scales as defined in Equation 1. The strength of the correlation is striking. The low scatter around the power law relation shown in this figure bolsters our claim that the shape of the $`\tau `$ distribution function is insensitive to other cosmological parameters. To give an idea of the uncertainty in each point, we fit both the high and low resolution simulations for the sCDM $`\sigma _8=0.3`$ and $`\sigma _8=0.7`$ models. In both of these cases $`\sigma _\tau `$ differs by less than 10 %.

### 3.2 Flux Probability Distribution

Although the optical depth PDF is easier to model theoretically, the flux PDF (where $`dP/dF`$ is the probability that a pixel will have transmitted flux between $`F`$ and $`F+dF`$) is closer to actual observation. The top panel of Figure 7 shows the flux probability distribution functions for the high spatial resolution sCDM model with $`\sigma _8=0.7`$ at $`z=2,3`$ and $`4`$. The bottom panel of Figure 7 shows the prediction of the simple $`\tau `$ scaling given in Equation 12 applied to the flux and these same flux probability distributions. Again we attempt to rescale the $`z=4`$ and $`z=2`$ distributions to $`z=3`$ in order to test the scaling. This results in a highly nonlinear mapping of the flux and the flux PDF from $`z=j`$ to $`z=3`$ given by $`dP_j/dF\eta F^{11/\eta }dP_j/dF`$ and $`F_jF_j^{1/\eta }`$ for $`\tau _j\tau _j/\eta `$, where $`\eta =0.356(6.511)`$ for $`j=2(4)`$, respectively. While the shapes of the distributions in the top panel appear quite different, much of the $`z`$ evolution of the flux probability distributions is explained by this simple scaling, the remainder representing mostly the effect of the evolution of the baryon density in the cosmological model. We do not plot the scaled distribution for $`z=4`$ below the scaled flux of $`0.5`$ because this already corresponds to an unscaled flux of $`0.015`$, close to saturation and most likely noise dominated in the observations.

The flux PDF depends only weakly on simulation grid resolution (Bryan et al. 1998). Its shape is strongly dependent on the power spectrum of the underlying cosmology. In the top panel of Figure 8 we show the flux probability distributions in the sCDM model (spatial resolution $`\mathrm{\Delta }x=37.5`$ kpc) for cluster scale normalizations $`\sigma _8=0.3`$ and $`0.7`$. The dependence on the normalization of the power spectrum is clear. The number of pixels found with flux in the central flux range $`0.3<F<0.9`$ is greater for models with less power ($`\sigma _8=0.3`$); while the number of pixels with flux in the low ($`F<0.3`$) and high ($`F>0.9`$) flux ends of the distribution are less than for models with greater power ($`\sigma _8=0.7`$). This is in qualitative agreement with Croft et al. (1997a). We note, however, that our result (using the Kronos code) does not require any smoothing of the simulations as was the case for their TreeSPH simulations. In the lower panel of Figure 8 we present the $`z=3`$ flux PDFs for the five models of the model comparison study. Models with lower power at small scales (tCDM, CHDM) have a larger flux PDF for $`0.3<F<0.9`$ than sCDM and $`\mathrm{\Lambda }`$CDM , while the low density model (OCDM) with the highest spectral power at these scales has the smallest flux PDF in this range, as expected. Furthermore, the differences between models can be substantial. For example, at $`F=0.6`$ the OCDM results lie $`10`$ % below the $`\mathrm{\Lambda }`$CDM model result while the CHDM result lies above the $`\mathrm{\Lambda }`$CDM result by about a factor of $`1.4`$. We remind the reader that the mean of the distribution has been fixed to match observations. Thus this statistic should be useful to constrain competing models.

### 3.3 Fraction of high Ly$`\alpha `$opacity

Another possibly useful statistic for discriminating models is the fraction of a quasar spectrum with Ly$`\alpha `$optical depth greater than a specified value $`\tau _0`$, i.e. the cumulative distribution in optical depth. Small differences in the amplitude normalization of the primordial power spectrum may be enhanced in the cumulative opacity data (Cen 1997).

Figure 9 shows the linear correlation between the opacity at line center and the column density of absorption features in the sCDM model, ranging from the optically thin to thick at the Lyman edge. The nearly unbroken relation $`\tau _cN_{HI}`$, which exists down to the incompleteness density of $`10^{12}`$ cm<sup>-2</sup> is attributed to the weak correlation between the Doppler parameter and column density since, in general, $`N_{HI}b\tau _c`$. The lower bound on the opacity ($`\tau _c>0.05`$) is set by the transmission or spectral threshold $`F_t=e^{\tau _t}=0.95`$ used in the line identification procedure. Using Figure 9 as a guide, we investigate the cumulative opacity distribution with the following minimum opacity thresholds: $`\tau >`$ 0.1, 1, and 7 which, if associated with the line centers, would correspond roughly to column densities of $`\mathrm{log}N_{HI}=`$ 12.5, 13.5 and 14.5 respectively. The distributions $`P(\tau >\tau _0)`$ for the above minimum opacity thresholds are plotted in Figure 10 at redshifts $`z=2`$, $`3`$, and $`4`$ for the models in the model comparison study. In comparing groups with the same $`\tau _0`$, the smaller threshold curves are more highly clustered and less sensitive to the background cosmological model parameters. This is especially evident in Figure 11 where we show the cumulative distributions of the optical depth at redshift $`z=3`$ for these models.

## 4 Line Parameter Statistics

In this section we present a line analysis of the spectra generated from the various model simulations. We compare and contrast the cosmological models based on the column density distribution and the evolution of line number.

### 4.1 $`\mathrm{I}`$Column Density Distribution

One of the most robust line statistics used in the analysis of the Ly$`\alpha `$forest is the $`\mathrm{I}`$column density distribution, which is well converged by simulation box sizes of $`9.6`$ Mpc and is insensitive to changes in the simulation grid resolution or treatment of gas hydrodynamics (Bryan et al. 1998, Zhang et al. 1997). The $`\mathrm{I}`$column density distribution, defined to be $`N_{HI}=_{x_A}^{x_B}\frac{n_{HI}}{1+z}𝑑x`$, is closely related to the optical depth $`\tau `$ through the dependence of each on the number density of neutral hydrogen. Thus using Equation 10 and the approximations that led to Equation 11 we expect the $`\mathrm{I}`$column density to scale as

$$N_{HI}\frac{(\mathrm{\Omega }_Bh^2)^2}{\mathrm{\Gamma }(z)H}(1+z)^{4.8}(1+\delta _b)^{1.7}𝑑u$$
(14)   
In the top panel of Figure 12 we show the raw (uncut) $`\mathrm{I}`$column density distribution for the high resolution sCDM model for redshifts $`z=2,3`$ and $`4`$. In the bottom panel of Figure 12 we see that the column density scaling

$$N_{HI}\frac{(\mathrm{\Omega }_Bh^2)^2}{\mathrm{\Gamma }(z)H}(1+z)^{4.8},$$
(15)   
the same relation as the naive scaling relation given in Equation 12 for the optical depth $`\tau `$, accounts for the redshift evolution of the column density distribution amazingly well. This demonstrates that the column density, an integrated quantity, is much less sensitive than the optical depth distributions to the redshift evolution of the gas overdensity within an absorbing structure. The differences seen in the low column density end of the distributions, particularly for $`z=4`$, may be a result of the simulation spatial resolution (Bryan et al. 1998), while the differences observed in the high column density end may partially be due to shot noise in the high $`z`$ data. In Figure 13 (top) we explore the dependence of this distribution at a given redshift ($`z=3`$) on the power spectrum of the underlying cosmology. We find qualitative agreement with semi-analytic arguments (Gnedin 1998; Hui, Gnedin, & Zhang 1997) in that models with less power on small scales (such as sCDM with $`\sigma _8=0.3`$ and $`\sigma _{34}=0.812`$) have $`\mathrm{I}`$column distributions with significantly steeper slopes than models (such as sCDM with $`\sigma _8=0.7`$ and $`\sigma _{34}=1.89`$) with more power at these scales. However, as we discuss in more detail below and in Table 2, quantitative agreement between the simulations and the predictions of these semi-analytic arguments seems more difficult to achieve.

In Figure 13 (bottom), we show the $`\mathrm{I}`$column density distribution at redshift $`z=3`$ for (Kronos) simulated spectra in the model comparison study and compare the simulated data with data from Kirkman & Tytler (1997) and the fits provided by Kim et al. (1997). The distributions are conventionally quantified by fitting them to power laws, $`dN/dN_{\mathrm{HI}}N_{\mathrm{HI}}^\beta `$. We use the same sets of column density cuts on the simulated data in the model comparison study as Kim et al. in order to expedite comparison with the data and use a direct unweighted least squares fit (all quoted errors are $`2\sigma `$) to extract the slope $`\beta `$ from the simulated data. Our results are summarized in Table 2. We find again the expected dependence on the fluctuation power spectrum. For the column density range $`13.7<\mathrm{log}N_{\mathrm{HI}}<14.3`$ (given by the column labeled $`\beta _\mathrm{h}`$) the most shallow slope is for OCDM, the low matter density model with $`\sigma _{34}=2.50`$ while CHDM and tCDM with $`\sigma _{34}=1.14`$ and $`1.09`$, respectively, give the steepest distributions (see Table 1). The predicted column density distributions generally also steepen with time (decreasing redshift). Kim et al. find $`\beta =1.46\pm 0.07`$ ($`2\sigma `$) for this column density range at $`z=2.85`$. This is formally inconsistent with all the models at the $`3\sigma `$ level except for OCDM, although it is marginally consistent with $`\mathrm{\Lambda }`$CDM .

Results for the column density range $`12.8<\mathrm{log}N_{\mathrm{HI}}<14.3`$ in Table 2 are shown in the column labeled $`\beta _{\mathrm{}}`$. The average distributions are generally shallower when extended to lower column densities, showing that the distributions are curved. Kim et al. similarly find a shallower distribution over this column density range with their results for lines at $`z=2.31`$, $`z=2.85`$, and $`z=3.35`$ shown in the last row of Table 2. The result at $`z=2.31`$ is inconsistent with all of the simulation results at $`z=2`$, but note that the quoted uncertainty in the observation is eight times smaller than at $`z=2.85`$, despite comparable numbers of absorbers. By contrast, the result at $`z=2.85`$ is formally consistent at the $`3\sigma `$ level with all the models. At $`z=3.35`$, Kim et al. find $`\beta =1.59\pm 0.13`$ ($`2\sigma `$). This result is consistent with the simulation results for the OCDM model and marginally consistent (at the $`3\sigma `$ level) for the sCDM and $`\mathrm{\Lambda }`$CDM models. The observational data, however, also suggest a weak steepening of the distribution with increasing redshift, contrary to our findings. These discrepancies might indicate that the redshift evolution of the ionizing radiation field may be somewhat different from that of the Haardt & Madau spectrum assumed in the simulations.

We may compare the simulation results with the semi-analytic predictions of Hui, Gnedin, & Zhang (1997) to understand the trend of changing steepness with power spectrum. We provide the predicted values of $`\beta `$ according to the prescription of Hui et al. in Table 2. We assume $`T\rho _B^{0.5}`$, as found by Zhang et al. (1998) for this column density range in an sCDM simulation. The uncertainty in the $`T\rho _B`$ relation introduces only a 10% uncertainty in the prediction for $`\beta 1`$, so it seems reasonable to retain it for the other models as well for this purpose. The predicted values of $`\beta `$ for the sCDM, tCDM, and CHDM models at $`z=3`$ match the simulation values to within $`1\sigma `$, in agreement with the comparison in Hui et al. with one of our earlier sCDM models. However, the predictions for OCDM and $`\mathrm{\Lambda }`$CDM at $`z=3`$ are in disagreement with the semi-analytic arguments giving too steep a slope. The predictions do less well for all models at $`z=2`$. In particular, the simulation results show a steepening of the column density distribution toward decreasing redshifts for all the models, opposite to the predicted trend.

Over the wider column density range $`10^{12.8}<N_{\mathrm{HI}}<10^{16}`$ (summarized in the column labeled $`\beta _\mathrm{f}`$ in Table 2), we see that the average distributions continue to steepen toward higher column densities. Kim et al. obtain $`\beta =1.46`$ for this column density range at $`z=2.85`$, with a steepening to $`\beta =1.55`$ at $`z=3.7`$. The results for the tCDM and CHDM models ($`1.95\pm 0.06`$ and $`1.92\pm 0.06`$ at $`z=3`$, respectively) are substantially steeper than these values. Because the distribution deviates from a pure power law at the low column density end, it is useful to split the simulation samples into two halves, fitting each to a power law. These results are given in the last two columns of Table 2 where $`\beta _1`$ is the slope of the column density distribution for lines at the low column density end ($`10^{12.8}<N_{\mathrm{HI}}<10^{14}`$) and $`\beta _2`$ is the slope of the column density distribution for lines at the high end ($`10^{14}<N_{\mathrm{HI}}<10^{16}`$). Giallongo et al. (1996) obtain $`\beta =1.8`$ for systems with $`N_{\mathrm{HI}}>10^{14}`$ and $`2.8<z<4.1`$.

Finally we note that the analogous column density distributions derived from the lower resolution Hercules runs give similar results and slopes as the Kronos data. For example over the full column density range at $`z=3`$, Hercules data give slopes for the column density distribution of $`1.71`$, $`1.66`$ and $`1.62`$ for the sCDM, $`\mathrm{\Lambda }`$CDM and OCDM models, respectively, consistent within errors with the Kronos results. This suggests that the distribution function is a robust diagnostic, being relatively insensitive to grid resolution and numerical method. A preliminary comparison with the data favors models with more power at these scales than in our CHDM or tCDM cosmologies. However, there appears to be some discordance in the observations, so a more definitive comparison will require more work.

### 4.2 Line Number Evolution

The number of Ly$`\alpha `$lines at a particular redshift reveals how many intergalactic absorbers exist at that time between the quasar and observer and, given certain assumptions on their geometry, the size and volume filling factor of the absorbers can also be deduced. Since the column density of Ly$`\alpha `$lines corresponds to the mean overdensity and size of the clouds fairly well (Charlton et al. 1997; Zhang et al. 1997), it is useful to see how the number of lines evolves with different column density cutoffs, as this will track the evolution of morphologically distinct small scale structures in the universe.

Figure 14 shows the evolution of the number of lines with $`\mathrm{I}`$column densities greater than $`10^{13}`$ cm<sup>-2</sup> , $`10^{13.5}`$ cm<sup>-2</sup> , and $`10^{14}`$ cm<sup>-2</sup> , respectively, comparing results for the models in the model comparison study with the observed data from Kulkarni et al. (1996) at $`z2`$, Hu et al. (1995) at $`z3`$ and Lu et al. (1997) at $`z4`$. For a fixed transmission cutoff (here $`F_t=0.95`$) and column density threshold, the total number of lines per unit redshift decreases with time because the opacity of the universe decreases from both the increasing flux of radiation and the expansion of the universe. With the exceptions of the tCDM and CHDM models for the lowest column density threshold $`N_{\mathrm{HI}}>10^{13}`$ where incompleteness due to line blending becomes significant at higher $`z`$, the deviation from a fixed power law behavior tracks predominantly the behavior of the radiation flux. Fitting the evolutions to the form $`dN/dz(1+z)^\gamma `$ over the range $`2<z<4`$, we find the exponents are fairly similar in the different models. We summarize these results in Table 3 (all errors are $`2\sigma `$). To compare these simulated results with observational data we fit the combined line lists from Kulkarni et al. , Hu et al. , Kirkman & Tytler, and Lu et al. to the same power law behavior and display those results in the row labeled “combined” in Table 3 again with $`2\sigma `$ errors. (Lines near the QSO emission redshift were avoided because of the proximity effect, as were lines associated with metal systems.) Kim et al. obtain $`\gamma =2.78\pm 1.42`$ ($`2\sigma `$), fit over $`2<z<3.5`$ for systems with column densities $`10^{13.77}<N_{\mathrm{HI}}<10^{16}`$ cm<sup>-2</sup> . Using the same column density cuts and simulation data for $`z=2`$ and $`3`$ only, we find power law exponents (labeled $`\gamma _\mathrm{h}`$ in Table 3) for the sCDM, $`\mathrm{\Lambda }`$CDM , OCDM, tCDM and CHDM models of the comparison study in good agreement with the observational result. For the lower column density range $`10^{13.1}<N_{\mathrm{HI}}<10^{14}`$ cm<sup>-2</sup> , the power law exponents are labeled $`\gamma _{\mathrm{}}`$ in Table 3. Kim et al. (1997) obtain $`\gamma =1.29\pm 0.90`$ ($`2\sigma `$) fit over $`2<z<4`$. Our sCDM, $`\mathrm{\Lambda }`$CDM and OCDM model predictions are in good agreement with this observational result, although tCDM and CHDM show a somewhat stronger evolution.

All the models yield comparable levels of evolution at each of the column density cutoffs and, in fact, the evolution slopes in all the models agree for the most part within errors with the observed values. Two trends in the model predictions are apparent. The first is that at a given column density threshold the slope $`\gamma `$ of the line number evolution is correlated to the slope $`\beta `$ of the column density distribution with $`\gamma `$ increasing for models with larger $`\beta `$ (i.e. models with less power in the fluctuation spectrum at small scales). The second trend is that stronger column density lines for a given model exhibit a greater rate of redshift evolution. This is not unexpected since the evolution for a fixed transmission threshold is essentially determined by the radiation field (which is the same in all the models) and, to a lesser degree, the expansion of the universe. In our previous studies (Zhang et al. 1997) we have found little intrinsic cloud evolution over these redshift intervals. If we assume that the evolution of the gas overdensity does not contribute significantly to the evolution of line number or column density at these redshifts we can use Equation 15 and the power law dependence of the column density distribution $`dN/(dN_{\mathrm{HI}}dz)N_{\mathrm{HI}}^\beta `$ to predict the number of lines above a fixed column density threshold as a function of redshift $`z`$. We find

$$\frac{dn}{dz}(>N_{HI})\left(\frac{(1+z)^{4.8}}{\mathrm{\Gamma }(z)H(z)}\right)^{(\beta 1)}$$
(16)   
In Figure 15 we compare the scaling predictions of Equation 16 to the simulation results using the high resolution sCDM $`\sigma _8=0.7`$ model for column density thresholds $`N_{\mathrm{HI}}>10^{13}`$, $`10^{13.5}`$, and $`10^{14}`$ cm<sup>-2</sup> , respectively. Since $`\beta `$ also evolves weakly with $`z`$, we use $`\beta `$ at $`z=3`$ as representative of the average for $`2<z<4`$ in Equation 16. For the lowest two column density thresholds we use the single power law fit to the column density distribution over the range $`10^{12.8}<N_{\mathrm{HI}}<10^{16}`$ cm<sup>-2</sup> , while we use $`\beta _2`$ from the two power law fit for the high ($`N_{\mathrm{HI}}>10^{14}`$ cm<sup>-2</sup> ) column density threshold. We normalize the scaling predictions to the simulated number of lines at $`z=3`$ because that is where all the models in our study were normalized to the observational data. For the lower two column density thresholds, the scaling prediction tends to overestimate the number of lines at $`z=4`$. For the lowest column density threshold ($`N_{\mathrm{HI}}>10^{13}`$) this again is partly due to incompleteness in the simulation line lists caused by line blending, an effect that becomes more severe for low column densities at high $`z`$. Furthermore the slope of the low column density end softens for $`N_{\mathrm{HI}}<10^{14}`$ with the break at $`N_{\mathrm{HI}}10^{14}`$ probably reflecting a change in the absorbers from low density structures evolving primarily with the universal expansion to structures undergoing gravitational collapse (Bryan et al. 1998). This deviation of the column density distribution from the pure power law assumed in the scaling relation would also cause the scaling law to overproduce the low column density lines. For the column density threshold $`N_{\mathrm{HI}}>10^{13.5}`$ the discrepancy between the scaling prediction and the simulation results at $`z=4`$ is reduced. This is to be expected since the high spatial resolution sCDM linelists should be complete at this column density threshold. For lines with $`N_{\mathrm{HI}}>10^{14}`$, where the absorbers share a common morphological type and the column density distribution is well fit by a single power law, agreement between the scaling prediction and the simulations is good. In the lower panel of Figure 15 we compare the scaling predictions with the simulation results for the models in the model comparison study for the high column density threshold case. The scaling predictions for all models agree reasonably well with the simulations.

Encouraged by these results we solve Equation 16 for the shape of the UV ionizing background $`\mathrm{\Gamma }(z)`$ in terms of the Hubble parameter $`H(z)`$ modeling the universal expansion and the (in principle) measurable quantities $`dn/dz(N_{\mathrm{HI}}>10^{14})`$ and $`\beta _2`$, the slope of the column density distribution over this column density range.

$$\mathrm{\Gamma }(z)\frac{(1+z)^{4.8}}{H(z)}\left(\frac{dn}{dz}\right)^{1/(1\beta _2)}$$
(17)   
We use Equation 17 to compute $`\mathrm{\Gamma }(z)`$ with simulation data from the model comparison study and compare these predictions to the Haardt–Madau spectrum actually used in Figure 16. Although the prediction is highly sensitive to the slope of the column density distribution used (whose errors are still quite large), it is gratifying, given the simplicity of the scaling relations, that all of the models reproduce the assumed Haardt-Madau evolutionary trend.

## 5 Doppler $`b`$ Parameter

Recent papers (Bryan et al. 1998, Theuns et al. 1998) have shown that both the Doppler $`b`$ parameter and a related nonparametric statistic, the mean flux difference as a function of velocity, require high simulation spatial resolution to model properly. In this section we investigate the dependencies of these statistics on the properties of the cosmological model. Because these statistics are highly sensitive to the spatial resolution of the simulation, we present results only for those models (sCDM, $`\mathrm{\Lambda }`$CDM , OCDM, and tCDM) simulated with our highest spatial resolution $`\mathrm{\Delta }x=37.5`$ kpc.

The Doppler $`b`$ parameter measures the amount of line broadening due to thermal broadening, physical velocities, Hubble expansion broadening and the shape of the absorber density profile (Bryan et al. 1998). Both Hubble and thermal broadening are significant for the lower column density lines that arise from structures found in voids that are still expanding in absolute coordinates. The thermal contribution only becomes dominant for the higher column density lines that have turned around and are gravitationally collapsing. Furthermore, the $`b`$ parameter is highly sensitive to the simulation spatial resolution. Lower resolution simulations numerically thicken the lines causing the width of the lines, and thus $`b`$ to be overestimated (Bryan et al. 1998; Theuns et al. 1998). In our previous work (Bryan et al. 1998) we argued that the shape of the $`b`$–distribution was in rough agreement with observation and particularly that the high $`b`$ power law tail of the distributions arises naturally in hierarchical models when quasar lines of sight pass obliquely through the filamentary absorbing structures (Rutledge 1998). However the median of the simulated $`b`$ parameter distribution for the sCDM model, calculated from simulations with high spatial resolution, was now substantially below the $`30`$ km/s median seen in the observations. Thus the sCDM model, which previously had appeared to be in agreement with observations, is now discrepant. In Figure 17 we show the Doppler $`b`$ parameter distributions extracted from the high grid resolution ($`\mathrm{\Delta }x=37.5`$ kpc) Kronos simulations at redshift $`z=3`$ for the sCDM, $`\mathrm{\Lambda }`$CDM , tCDM, and OCDM models for lines with column densities between $`10^{13.1}`$ cm<sup>-2</sup>$`<N_{\mathrm{HI}}<10^{14}`$ cm<sup>-2</sup> . We present for comparison data from Kim et al. (1997) for $`z=3.35`$. The $`\mathrm{\Lambda }`$CDM and OCDM models, like sCDM, have their $`b`$ distributions shifted too much to the left (to low $`b`$ values) to agree with observation. Only tCDM, the model with the least fluctuation power at small scales and thus broader density structures at this redshift, has a median $`b`$ approaching the observational values. We explore this dependence on the fluctuation power spectrum with the two highest resolution sCDM models (with $`\sigma _8=0.7`$ and $`0.3`$, respectively) in the lower panel of Figure 17 and see that indeed the model with lower spectral power produces a $`b`$ parameter distribution shifted towards higher $`b`$ (as predicted by Hui & Rutledge 1998). The increase in $`b`$ for models with less fluctuation power at small scales may be partly due to line blending effects at these low column densities. However as shown below, the shift to higher $`b`$ values for these models persists for lines with higher column densities as well where line blending should not be as significant and thus can not be explained by line blending alone.

To facilitate a better comparison of the models with observations we plot the median Doppler parameters as a function of redshift in Figure 18 where we have imposed the same column density cuts on the lines as those used by Kim et al. (1997). The median $`b`$ for lines with column densities $`10^{13.8}`$ cm<sup>-2</sup>$`<N_{\mathrm{HI}}<10^{16}`$ cm<sup>-2</sup> and $`10^{13.1}`$ cm<sup>-2</sup>$`<N_{\mathrm{HI}}<10^{14}`$ cm<sup>-2</sup> are shown in the top and bottom panels, respectively. While the $`\mathrm{\Lambda }`$CDM , sCDM, and OCDM models predict roughly the observed evolutionary trend for both sets of column density cuts, the median $`b`$ values lie systematically more than $`6`$ km/s below the observational data. OCDM, the model with the most power at these scales, is the most discrepant. Although tCDM predicts median $`b`$ parameters more consistent with observation, the redshift evolution predicted by this model appears to be in disagreement with the data. We can compare these results with other recent data sets. Confining to lines with $`N_{\mathrm{HI}}>10^{13}`$ cm<sup>-2</sup> , we obtain from the published line lists, for $`1.9<z<2.0`$ ($`1\sigma `$ errors), $`(b_{\mathrm{mean}},b_{\mathrm{median}})=(32.1\pm 2.6,29.7\pm 3.3)`$ km s<sup>-1</sup> ($`1\sigma `$) (Kulkarni et al. ), $`3<z<3.1`$ $`(b_{\mathrm{mean}},b_{\mathrm{median}})=(38.0\pm 1.6,33.6\pm 2.0)`$ km s<sup>-1</sup> (Hu et al. ), $`(b_{\mathrm{mean}},b_{\mathrm{median}})=(27.3\pm 1.9,25.9\pm 1.2)`$ km s<sup>-1</sup> (Kirkman & Tytler), and at $`4<z<4.1`$ $`(b_{\mathrm{mean}},b_{\mathrm{median}})=(32.6\pm 2.4,25.9\pm 3.1)`$ km s<sup>-1</sup> (Lu et al. ), again clearly discrepant with the model predictions. Thus none of the models considered here can restore agreement with the observational data.

We also argued in Bryan et al. (1998) that this discrepancy is not a result of the particular choice of line fitting algorithm, but appears for sCDM in the nonparametric moments of the two point flux distribution functions as well. The two-point function $`P_2(F_1,F_2,\mathrm{\Delta }v)`$ gives the probability that two pixels with separation $`\mathrm{\Delta }v`$ will have flux $`F_1`$ and $`F_2`$. We plot the normalized moments of this function averaged over the flux range $`F_a`$ to $`F_b`$ given by

$$\frac{_{F_a}^{F_b}𝑑F_1_0^1𝑑F_2P_2(F_1,F_2,\mathrm{\Delta }v)(F_1F_2)}{_{F_a}^{F_b}𝑑F_1_0^1𝑑F_2P_2(F_1,F_2,\mathrm{\Delta }v)}$$
(18)   
which represents the average flux diffence as a function of velocity for pixels in the range $`F_a`$ to $`F_b`$. In Figure 19 we plot the above statistic at $`z=3`$ as a function of velocity for several flux ranges for the high resolution models (sCDM, $`\mathrm{\Lambda }`$CDM , OCDM, tCDM) of the model comparison study (lower panel) and study its dependence on the small scale fluctuation power (top panel) using sCDM with $`\sigma _8=0.7`$ and $`0.3`$, respectively. There is little difference for low velocities, independent of flux level, due to the high coherence of the lines. At very large velocity differences there is no coherence and the value is just the difference between the mean value of the transmitted flux and the mean flux in a given flux interval (Bryan et al. 1998). It is at intermediate velocity separations where the statistic is heavily influenced by the structure of the lines. There $`\mathrm{\Lambda }`$CDM , OCDM, and sCDM with $`\sigma _8=0.7`$, whose power spectra are very similar, produce very similar distributions; while the tCDM model is quite distinct. We may quantify these model differences by determining at what $`\mathrm{\Delta }v`$ the model prediction passes through a given average flux difference. For the flux interval $`0<F<0.1`$ the simulation predictions pass through the mean flux difference of $`0.3`$ at $`\mathrm{\Delta }v35`$ km/s for $`\mathrm{\Lambda }`$CDM , OCDM, and sCDM $`\sigma _8=0.7`$ and at $`\mathrm{\Delta }v45`$ km/s for tCDM. Although observational data is limited, these are all lower than the $`\mathrm{\Delta }v55`$ km/s from Figure 3 of Miralde-Escudé et al. (1997).

It is important to ask what is needed to restore agreement between the simulations and observations. Although we can not completely rule out the possibility that the line fitting algorithm contributes to differences in the simulated and observed $`b`$-parameter distributions, we argue that its effect should not be significant because the discrepancy is seen at a comparable level (Bryan et al. 1998) in the fit–independent two–point distribution of the flux as well. The mean optical depth of our models was scaled to agree with observations, but this normalization is in some dispute. However, changing this normalization has little effect on the median of the $`b`$-distribution. For example, using sCDM ($`\sigma _8=0.7`$), an increase in $`\overline{\tau }_\alpha `$ from $`0.225`$ to $`0.35`$ at $`z=3`$ causes the median $`b`$ value to decrease from $`20.8`$ km/s to $`20.1`$ km/s, a change of only $`1`$ km/s. One possibility might be to change the ionization history of the universe such that pressure broadening would widen the absorbing structures. Another might be to change the density structure through the power spectrum of the cosmology itself. However, with the suite of models considered here it seems difficult for a single model to give good agreement with both the column density and $`b`$-parameter data.

## 6 Flux Statistics for Helium II

Previous work (Zhang et al. 1997,1998; Croft et al. 1997b) indicates that He $`\mathrm{II}`$Ly$`\alpha `$(304 Å) absorption may be significant in regions where $`\mathrm{I}`$Ly$`\alpha `$is not. Thus the study of He $`\mathrm{II}`$Ly$`\alpha `$absorption in quasar spectra provides a unique probe of structure in the lowest density regions of the universe. Comparison of both $`\mathrm{I}`$and He $`\mathrm{II}`$absorption within the context of a given cosmological model may also yield important information about the spectral shape of the metagalactic UV radiation field and its redshift evolution. While current observations still struggle to obtain sufficient resolution to detect any but the broadest individual He $`\mathrm{II}`$Ly$`\alpha `$lines, it is still possible to determine mean statistics of the He $`\mathrm{II}`$flux and optical depth which are not so sensitive to instrumental resolution. We define the mean optical depth $`\overline{\tau }_{HeII}=\mathrm{ln}F`$ where $`F`$ is the mean transmitted flux ($`F=1`$ signifying complete transmission). In Figure 20 we present $`\overline{\tau }_{HeII}`$ as a function of redshift for the sCDM models with varying power normalizations (top) and for the models in the model comparison study (bottom). Several trends are apparent. First all models produce a rapid rise in mean optical depth with increasing redshift (roughly a factor two between $`z=2`$ and $`z=3`$), with tCDM and CHDM rising slightly more steeply. This is consistent with previous work on a smaller number of hierarchical cosmologies (Zhang et al. 1997; Croft et al. 1997b) and with the interpretation that the observed optical depth is due primarily to absorption by gas in underdense regions. The redshift evolution of the optical depth is thus dominated by the change in the gas density due to universal expansion and (to a lesser degree) by the shape of the UV metagalactic ionizing background (here assumed to be that of Haardt-Madau (1996) with frequency dependence $`\nu ^{1.8}`$). Second, for a given redshift $`z`$, models with less power on small fluctuation scales have progressively larger optical depths. This is again consistent with the interpretation that the absorption is due to gas in predominantly underdense regions since less gas in these low power models will have turned around and collapsed.

The first observation of a flux decrement at the wavelength where the He $`\mathrm{II}`$Ly$`\alpha `$absorption should occur was made by Jakobson et al. (1994) using the HST Faint Object Camera to observe quasar Q0302-003. They obtained a 90% confidence lower limit of $`\overline{\tau }_{HeII}>1.7`$ at $`z=3.286`$. Subsequently improved measurements using spectra from this same quasar were made by Hogan et al. (1996) with the Goddard High Resolution Spectrograph on the HST and by Heap et al. (1998) using the Space Telescope Imaging Specrometer (STIS). STIS provides better sensitivity and background determinations than previous measurements. Davidsen, Kriss & Zheng (1996) used the Hopkins Ultraviolet Telescope to study the average He $`\mathrm{II}`$opacity in the spectrum of quasar HS1700+64 over the redshift interval $`2.2<z<2.6`$ (lower than that available with HST). They find $`\overline{\tau }_{HeII}=1.0\pm 0.07`$, although as shown in Figure 3 of Croft et al. (1997b) there is considerable scatter when the wavelength range is divided into $`10\AA `$ bins. Measurements of He $`\mathrm{II}`$absorption have also been made by Anderson et al. (1998) with STIS using the spectrum of quasar PKS 1935-692. Although the number of lines of sight studied so far are limited and thus a detailed comparison of observations with our model simulations (that average over 300 lines of sight) is premature, these data are also presented in Figure 20. On face value these data favor higher optical depths and thus models with lower fluctuation power. However, none of the simulation models presented here can reproduce the apparent break at $`z=3`$ in the optical depth observed by Heap et al. (1998). If this break persists it would most likely signal a departure from the Haardt-Madau quasar reionization spectrum assumed here.

For completeness and comparison with previous work (Croft et al. 1997b) we show in Figure 21 (from top to bottom) the He $`\mathrm{II}`$Ly$`\alpha `$flux probability distribution functions at $`z=4,3,2`$, respectively, for the models in the model comparison study. The distribution functions are calculated from the flux smoothed with window functions corresponding to the same Full Width at Half Maximum (FWHM) as STIS with high, 50 km s<sup>-1</sup> resolution (left column) and low, 500 km s<sup>-1</sup> resolution (right column). As Figure 21 shows, the shape of the flux PDF is highly dependent on the smoothing. We see however, that models with less fluctuation power on small scales have far fewer truly transparent regions (pixels with $`F`$ near $`1`$). For the high $`\mathrm{\Delta }v=50`$ km s<sup>-1</sup> resolution case, all models converge in the fully saturated regime ($`F<0.05`$) as expected.

## 7 Summary

We have performed several simulations of the Ly$`\alpha `$forest using different background cosmological models, numerical codes and grid resolutions. Five different cosmological models were considered here: the standard flat critical density cold dark matter model (sCDM), a flat CDM model with a nonzero cosmological constant ($`\mathrm{\Lambda }`$CDM ), an open CDM model (OCDM), a flat critical density CDM model with a tilted power spectrum matching both the COBE amplitude and small scale clustering contraints (tCDM), and a flat critical density mixed dark matter model (CHDM). The high resolution shock capturing code Kronos was used with grid resolution $`\mathrm{\Delta }x=37.5`$ kpc ($`\mathrm{\Delta }x=75`$ kpc for CHDM) for the benchmark calculations presented in this paper. Three of the models (sCDM, $`\mathrm{\Lambda }`$CDM , OCDM with identical parameters) were also evolved with the artificial viscosity based code Hercules at the lower grid resolution $`\mathrm{\Delta }x=75`$ kpc. Both simulation techniques give similar results for statistics, such as the slope of the column density distribution, that are insensitive to grid resolution.

We have presented results from several statistical analyses of absorption features present in the Ly$`\alpha `$spectra, both from the unprocessed optical depth data and from the reduced line lists. Explicitly we have considered the optical depth and transmitted flux probability distribution functions, the cumulative optical depth distributions, the $`\mathrm{I}`$column density distributions, line number evolution, Doppler $`b`$ parameter distribution, the average flux difference as a function of velocity (first moment of the 2-point flux distribution function), and the mean optical depth and flux probability distribution functions for He $`\mathrm{II}`$absorption. We find:

1. Simple scaling laws describe the redshift evolution of the optical depth, flux PDF, the $`\mathrm{I}`$column density distribution and, in conjunction with the slope of the column density distribution, the line number evolution remarkably well. This demonstrates that most of the evolution of the Ly$`\alpha `$forest is a direct consequence of universal expansion.
2. The shape of the optical depth PDF is strongly correlated to the amplitude of the density fluctuation spectrum. Differences between models my be significant in the observationally accessible region $`0.02<\tau <4`$. Thus this statistic may be a useful discriminator among models. Similar conclusions hold for the related flux PDF.
3. Cumulative opacity distributions for the models are strongly clustered at low optical depth thresholds and high $`z`$. Significant differences do occur for optical depth thresholds $`\tau _0>1`$, but these may be more difficult to observe.
4. The column density distribution function is a robust statistic relatively insensitive to grid resolution and numerical method. Its redshift evolution is described well by the same naive scaling law that describes the evolution of the optical depth. The slope of the column density distribution is sensitive to the amplitude of the power spectrum on scales roughly the size of the absorbers ($`100`$ kpc). Models with less power at these scales produce steeper distributions in qualitative agreement with semi- analytic arguments (Hui, Gnedin & Zhang 1997). A preliminary comparison with data favors models with more power (sCDM, $`\mathrm{\Lambda }`$CDM , OCDM) over those with less power (tCDM, CHDM).
5. All models show comparable evolution for the number of lines above a given $`\mathrm{I}`$column density threshold in reasonable agreement with the data. Thus this statistic is not a sensitive discriminator among models.
6. Although the shape of the Doppler $`b`$ parameter distribution is well reproduced by all the models, the median of the distribution for sCDM, $`\mathrm{\Lambda }`$CDM , and OCDM models is well below observed values. The median of the $`b`$ parameter increases for models with less power on small scales. Thus the observations favor low power models, such as tCDM, making it difficult for any model considered in this study to simultaneously give good agreement with both the $`\mathrm{I}`$column density and $`b`$ parameter data. This discrepancy is confirmed as well in the nonparametric first moment of the two-point flux distribution function and so is not solely the result of the line-fitting algorithm employed. The solution to this problem may require modification of the reionization history of the universe to produce more pressure broadening of the absorbing structures or a modification of the power spectrum of the underlying cosmology itself.
7. All of the models simulated in this study produce a rapid rise in He $`\mathrm{II}`$mean optical depth with increasing redshift consistent with the interpretation by previous work (Zhang et al. 1997, Croft et al. 1997b) that the observed optical depth is due to absorption by gas in underdense regions where universal expansion dominates the evolution of the gas density. Models with less power on small scales (tCDM, CHDM) produce larger mean He $`\mathrm{II}`$optical depths. Preliminary comparison with the data tends to favor these low power models. However, none of the models can reproduce the break seen by Heap et al. near $`z=3`$. If this break persists in the data, it would most likely reflect that the Haardt- Madau (1996) form for the metagalactic UV ionizing background, based on homogeneous reionization by quasars alone in a clumpy medium, must be modified.

###### Acknowledgements.

This work is supported in part by NSF grant AST-9803137 under the auspices of the Grand Challenge Cosmology Consortium (GC<sup>3</sup>). NASA also supported this work through Hubble Fellowship grant HF-0110401-98A from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc. under NASA contract NAS5-26555. M.M. acknowledges support from the Research and Scholarship Development Fund of Northeastern University. The computations were performed on the Convex C3880, the SGI Power Challenge, and the Thinking Machines CM5 at the National Center for Supercomputing Applications, and the Cray C90 at the Pittsburgh Supercomputing Center under grant AST950004P.