Equation of state of the intergalactic medium in the early Universe

Spectroscopy of the Ly$\alpha$ forest in quasar spectra proved to be a useful tool for probing the intergalactic gas. We developed the automatic program for Voigt profile fitting of Ly$\alpha$ forest lines. We run this code on 9 high resolution ($\sim 50000$) quasars spectra with a high signal-to-noise ratio ($\sim 50 -100$) from the Keck telescope archive and obtained the sample of single well-fitted Ly$\alpha$ lines. Fitting the joint 2d distribution of column density and Doppler parameter from this sample by physically reasonable model we estimate a power law index $\gamma$ of the equation of state of the intergalactic medium in the redshift range $z\sim 2-3$. We found that our measurement is in an agreement with measurements by other groups obtained with Voigt profile fitting technique.


Introduction
It is well known that the intergalactic medium (IGM) contains a dominant fraction of baryons in the Universe [1]. After the recombination epoch at a redshift z ∼ 1100 the intergalactic medium (IGM) had remained completely neutral until first stars, galaxies and quasars began to produce the ionizing background radiation. The IGM at this epoch mostly consisted from hydrogen and helium and consequently the IGM reionization included two major stages. The first stage is the reionization of H i which was completed at z > 6 [2] and dominated by starlight escaped from the galaxies [3]. At this epoch the helium was already singly ionized. The second stage is the reionization of He ii which occurs at later epochs, z < 3 [4], and is driven by quasars radiation [5]. Despite the considerable progress, our understanding of the reionization process is incomplete. For instance, we still poorly know the properties of the ionizing ultraviolet background, precise redshifts of the reionization events and their impact on the thermal state of IGM. It is assumed [6] that the equation of state -the temperature-density, T − ρ, relationafter the H i reionization epoch is set by the competing processes of photo-heating and adiabatic cooling and obeys a power-law form where T 0 is the temperature at the mean densityρ, ∆ is the overdensity and according to theoretical predictions. Numerical calculations suggest the index γ is about 1.6 well after the reionization [6,7]. The small residual neutral fraction of hydrogen remained after the reionization (at z < 6) leaves traces in the spectra of background quasars observed as series of the absorption features called the Lyα forest ( figure 1). An analysis of these lines allows to probe the thermal state of the IGM and its evolution with redshift. Various methods were suggested based on the statistical properties of the Lyα forest. These include analysis of the flux probability distribution function [8,9], the power-spectrum of the transmitted flux [10,11], the average local curvature [12,13] and wavelet decomposition of the Lyα forest [14,15]. In this paper we employ another widely-used technique which is based on a decomposition of the Lyα forest into individual components and their subsequent Voigt profile analysis [16][17][18][19][20][21][22]. Table 1 collects basic parameters of 9 selected high-resolution quasar spectra with high signal-tonoise ratio (SNR) we used for the analysis. Most of them was taken from KODIAQ database [23] and one spectrum, Q1442+2931, was reduced by hand from the Keck telescope archival data using MAKEE package 1 . An unabsorbed continuum level for the Q1442+2931 spectrum was estimated using B-spline interpolation, spectra from KODIAQ database are already normalized, but in some cases we adjusted continuum manually.

Data and analysis
We developed an automatic procedure for analysis of the H i absorption lines in the Lyα forest. In each quasar spectrum we used the spectral region between Lyα and Lyβ emission lines with offset of 12000 km/s bluewards from Lyα emission lines. This allowed us to avoid confusion with Lyβ forest lines and to minimize proximity effects for quasar. Each individual H i absorber in the Lyα forest was described by the Voigt profile, with the column density, N , Doppler parameter, b, and redshift, z. Instead of fitting each Lyα forest line in specified redshift range with multicomponent profile, we searched for the solitary lines that can be fitted with only one component (see example in the inlay of figure 1). We constructed the dense grid of parameters (z, N , b), see table 2, and used χ 2 criterion to compare model profiles calculated on that grid with observed ones. The Lyα forest lines in the probed redshift range, z ∼ 2 − 3, frequently blend with each other. Therefore, we used a sort of outlier rejection routine to include in analysis also the lines with weak blends in one of the wings. In addition we checked that the higher-order Lyman-series absorption lines are in agreement with the estimated parameters  The choice of steps in the grid was based on the Fisher matrix estimate.
from corresponding Lyα line. Finally, we excluded lines with parameters lying on the grid edges to get rid of the boundary effects on the line sample density estimation in the parameters space.
In this way we obtained the representative sample of "solitary" Lyα forest lines in the specified range of (N , b) parameters. Since our selection procedure is not the line profile fitting, it does not return the line parameter uncertainties σ N and σ b . These were obtained analytically via the Fisher matrix for the χ 2 likelihood function with account for resolution and SNR of each spectra. We checked that this gives reasonably close estimates of the expected uncertainties would the real fit be performed. The resulting sample in (N, b) parameter space is shown in figure 2. The mean redshift of the Lyα lines in our sample isz = 2.35.
The obtained (N, b) distribution shows a prominent lower envelope. Its physical origin is as follows [16]. The parameter b characterizes the width of the velocity distribution of hydrogen atoms, which contains thermal and turbulent motions inside each absorber. Since thermal broadening is unavoidable, for each N there exists a minimal value of b which corresponds to a pure thermal broadening where k B is the Boltzmann constant and m is the hydrogen atomic mass. The position of the low boundary of the N − b distribution can be connected to the temperature-density relation of the IGM, given by eq (1). In assumption of uniform background  radiation, the H i column density N is related to the overdensity ∆ as [24,25] where Γ −12 is the hydrogen photoionization rate in units of 10 −12 s −1 . Eqs (1), (2) and (3) allow to predict the b th − N dependence where the normalization constant b 0 depends on T 0 , Γ −12 and z. Once b 0 and Γ are found, the parameters γ(z) and T 0 (Γ −12 , z) can be inferred. Techniques that employ decomposition of Lyα forest lines are usually based on direct determination of the lower boundary of N − b distribution. This is not a straightforward task. Various iterative techniques are proposed to eliminate the turbulent contribution to width of the lines [17,19,20]. Most popular statistics is based on the so-called 2σ rejection algorithm [16], which in some cases is calibrated with the hydrodynamical simulations.
Here we propose to fit the whole (N, b) distribution. We use microturbulence assumption that thermal and turbulent motions are uncorrelated and can be characterized by the Gaussian distributions with widths b th and b turb , respectively. The convolution of these distributions is also Gaussian with the characteristic width b = b 2 th + b 2 turb . We assume that turbulent motions are independent on the thermal state of the gas and its density. Besides, we suggest that distributions over N and b turb have power-law shapes with indexes β and p, respectively.
Notice that for N distribution this is a standard assumption and index β is well measured, see, e.g., [18]. In this way the 2D probability density distribution has the form where f N and f turb are the probability density distributions over N and b turb , respectively. Thus the function f depends on the parameter vector Θ=(Γ, log b 0 , p, β) in addition to N and b.
To estimate the posterior distributions of the parameters Θ we used the Bayesian framework. Flat priors were chosen for all parameters. Our (N, b) sample probes only a limited part of the full distribution f (N, b), namely the region specified by the searching procedure box (see table 2). The estimated uncertainties σ N i and σ b i in the sample are heteroscedastic, therefore the probability of the data point to fall in this box is different for each line, since this probability depends on the observed (N, b) values and not on the unknown true values ( N , b). As a result, the correctly normalized likelihood function is [26] where f i (N, b| N

Results and discussion
The marginal posterior distributions of the fit parameters are shown in figure 3. According to eq (4), the estimated power-law index in the temperature-density relation (1) of the IGM at the average redshift z = 2.35 is γ − 1 = 0.53 ± 0.07 at the 68 per cent confidence level. Using eq (3) we find the temperature at the mean density log T 0 [K] = (3.98 ± 0.05) − (0.35 ± 0.04) log Γ −12 (7) which, however, depends on the photoionization rate Γ −12 . From the analysis of the mean opacity of the Lyα forest lines it was constrained that (T 0 [K]) −0.72 /Γ −12 = (1.88 ± 0.08) × 10 −3 at z = 2.4 [28]. Their calculations also depend on the equation of state of the IGM, however, their adopted value for γ − 1 = 0.6 is rather close to our 0.53 (see their figure 3), so we do not expect a significant bias. Adopting the result of [28] we can constrain both T 0 = 1.25 +0.05 −0.19 ×10 4 K and Γ −12 = 0.60 +0.08 −0.03 with a strong anticorrelation between these parameters. The sample is much smaller than some samples derived from the multicomponent fit (see e.g. [19,20]), but it allows us to avoid a lack of knowledge of "true" multicomponent models and correlations of the parameters of such compound models. This can be comprehended if one compares our sample with sample, obtained on the same quasar spectra [19]: at each column density region the Rudie et al detect much more lines with low b values than we did. But finally they rejected most of the lines with low b values by some arbitrary routine.
In previous sections we suggest that broadening of absorption lines determined by thermal and turbulent motions. Turbulent broadening is driven by peculiar velocities of the absorbed intergalactic clouds and in our approach it is independent on column density N . But in general real turbulent mechanisms may depend on N . For example we did not take into account expected dependence of b on N due to pressure smoothing scale and differential Hubble flow along to physical extent of intergalactic clouds [29][30][31].