# Thermal state of the intergalactic medium near to the optical limit for the
Ly-$\alpha$ forest

Tomáš Ondro1,2 and Rudolf Gális2
1Department of Agricultural, Food and Environmental Engineering, Faculty of
AgriSciences, Mendel University in Brno,
Zemědělská 1, 613 00 Brno, Czech Republic
2Institute of Physics, Faculty of Science, Pavol Jozef Šafárik University in
Košice, Park Angelinum 9, 040 01 Košice, Slovakia
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

In this article, the temperature-density relation of the intergalactic medium
was studied in the region $1.6\leq z<2.0$ divided into two bins. For this
purpose, the Ly-$\alpha$ forest decomposition into individual absorption
profiles was used for the study of 35 publicly available quasar spectra
obtained by the Ultraviolet and Visual Echelle Spectrograph (UVES) on the Very
Large Telescope (ESO) and by the High Resolution Echelle Spectrometer (HIRES)
on the Keck Telescope. For the determination of the thermal state sensitive
cut-off position in the $b-N_{\text{H\,{i}}}$ distribution, the iterative
fitting procedure was adopted. The measurements were calibrated using mock
Ly-$\alpha$ forest data generated by 23 hydrodynamical simulations with
different thermal histories. The value of the temperature at mean density
corresponds to the decreasing trend predicted by various models at the lower
redshifts. In the case of power law index, determined values are close to 1.6,
which is expected after all reionization events in various models assuming the
balance of photoheating with adiabatic cooling.

###### keywords:

cosmology: early Universe – quasars: absorption lines – intergalactic medium

††pubyear: 2020††pagerange: Thermal state of the intergalactic medium near to
the optical limit for the Ly-$\alpha$ forest–12

##  1 Introduction

The physical conditions of the intergalactic medium (IGM) could be studied by
analysis of the so-called Ly-$\alpha$ forest in quasar (QSO) spectra. Hui &
Gnedin (1997) showed that for gas overdensities $\delta\lesssim 10$, the
temperature is related to the density throughout a power-law in the form

$T=T_{0}(1+\delta)^{\gamma-1},$ (1)

where $T_{0}$ is the temperature of the IGM at the mean density and $\gamma$
is the power-law index. If we determine the parameters $T_{0}$ and $\gamma$ as
a function of redshift, we can describe the thermal history of the IGM.

Various approaches were used for the characterization of the $T-\rho$ relation
of the IGM from Lyman-$\alpha$ absorption forest. As an example, the approach
that treats the Ly-$\alpha$ forest as a superposition of discrete absorption
profiles (Schaye et al., 1999; Ricotti et al., 2000; Rorai et al., 2018; Hiss
et al., 2018), analysis of the flux probability distribution function by
Bolton et al. (2008) and Viel et al. (2009), the power spectrum of the
transmitted flux (Zaldarriaga et al., 2001; Theuns et al., 2002; Walther et
al., 2019), the average local curvature (Becker et al., 2011) and Boera et al.
(2014) and wavelet decomposition (Lidz et al., 2010).

Most studies deal with the characterization of the $T-\rho$ relation of the
IGM around $z\sim 3$, which is interesting due to the He II reionization
phase. Unlike of this redshift, for which is possible to find $\sim 300$
Ly-$\alpha$ absorption lines in QSO spectra, for $z\sim 2$, there are
typically less than 100 absorption lines with $\log{N_{\text{H\,{\sc
i}}}}\in(12.5-14.5)$ between the Ly-$\alpha$ and Ly-$\beta$ emission lines
(Schaye et al., 1999). Generally, the number density of the Ly-$\alpha$ forest
decreases with the decreasing redshift, and $z\simeq 1.5$ is considered as the
optical limit for the Ly-$\alpha$ forest (Boera et al., 2014).

The temperature in this lower redshift region was firstly studied by Boera et
al. (2014) using the curvature measurement. For the mean redshift $z=1.63$,
the authors determined the temperature $(20.66\pm 2.03)\times 10^{3}$ K and
$(13.00\pm 1.27)\times 10^{3}$ K at the mean density assuming the values of
$\gamma\sim$ 1.3 and 1.5, respectively. For the higher value of mean redshift
$(z=1.82)$, the authors determined the values of teplerature $T_{0}(\gamma\sim
1.3)=(17.61\pm 0.76)\times 10^{3}$ K and $T_{0}(\gamma\sim 1.5)=(11.79\pm
0.51)\times 10^{3}$ K. Based on the results from the redshift range
$1.5\lesssim z\lesssim 2.8$, the authors of this study found a decrease in
temperature with decreasing redshift, which could be interpreted as a trace of
completion of the reheating process related to the He ii reionization.

Walther et al. (2019) presented the fiducial evolution of thermal parameters
in the region $1.8<z<5.4$ using the Ly-$\alpha$ flux power spectrum. Assuming
the strong prior on mean transmitted flux $\overline{F}$, the corresponding
values of the temperature $T_{0}$ and power-law index $\gamma$ are
$(7.68^{+3.69}_{-2.18})\times 10^{3}$ K and $1.63^{+0.16}_{-0.25}$ at the
redshift $z=1.8$, respectively.

It is worth noting that Schaye et al. (2000) also studied the region $1.85\leq
z\leq 2.09$ using the spectrum of only one quasar (Q1100-264). For the median
value of the redshift $z=1.96$ , the authors determined value of $T_{0}\sim
11\,000$ K and $\gamma=1.4$. In this case, the analysis of IGM was based on
the Voigt profile decomposition of the Ly-$\alpha$ forest into the set of
individual absorption lines.

The aim of this study is to determine the coefficients in the $T-\rho$
relation of the IGM in the redshift range of $1.6\leq z<2.0$. The article is
organized as follows: In Section 2 we present the basic information about the
used data and the analysis of the spectra based on the Voigt profile fitting.
This section also contains the description of the metal line and narrow line
rejection, as well as the cut-off fitting procedure. The information about the
hydrodynamical simulations is given in Section 3. The THERMAL suite (Thermal
History and Evolution in Reionization Models of Absorption Lines) is also
introduce in this section. A description of the determination of coefficients
in the $T-\rho$ relation and the calibration procedure are presented in
Section 4. In Section 5 we present our results and their comparison with the
previously published ones and with simulations. Finally, our conclusions are
given in Section 6.

Table 1: List of QSOs whose spectra were used in this study. The S/N ratio was calculated according to Stoehr et al. (2008) for the spectral regions where the absorbers were parameterized. Object | $\alpha_{2000}$ | $\delta_{2000}$ | $z_{\rm{em}}$ | S/N | Source
---|---|---|---|---|---
| [h m s] | [d m s] |  |  |
J000344-232355 | 00 03 44.91 | –23 23 55.3 | 2.28 | 49 | UVES
J000443-555044 | 00 04 43.28 | –55 50 44.6 | 2.1 | 10 | UVES
J001602-001225 | 00 16 02.40 | –00 12 25.1 | 2.085 | 18 | UVES
J012417-374423 | 01 24 17.36 | –37 44 23.0 | 2.2 | 37 | UVES
J014333-391700 | 01 43 33.63 | –39 17 00.1 | 1.807 | 30 | UVES
J014944+150106 | 01 49 44.43 | 15 01 06.70 | 2.06 | 12 | KODIAQ
J022620-285750 | 02 26 20.49 | –28 57 50.7 | 2.171 | 10 | UVES
J022853-033737 | 02 28 53.21 | –03 37 37.1 | 2.066 | 11 | KODIAQ
J024008-230915 | 02 40 08.16 | –23 09 15.6 | 2.223 | 83 | UVES
J025634-401300 | 02 56 34.00 | –40 13 00.3 | 2.29 | 10 | UVES
J042707-130253 | 04 27 07.29 | –13 02 53.6 | 2.159 | 10 | UVES
J092129-261843 | 09 21 29.34 | –26 18 43.2 | 2.3 | 14 | UVES
J101939+524627 | 10 19 39.15 | 52 46 27.80 | 2.17 | 10 | KODIAQ
J103921-271916 | 10 39 21.84 | –27 19 16.4 | 2.23 | 25 | UVES
J110325-264515 | 11 03 25.29 | –26 45 15.8 | 2.145 | 61 | UVES
J110610+640009 | 11 06 10.74 | 64 00 09.60 | 2.203 | 29 | KODIAQ
J121140+103002 | 12 11 40.59 | 10 30 02.0 | 2.191 | 11 | UVES
J122824+312837 | 12 28 24.96 | 31 28 37.60 | 2.2 | 46 | KODIAQ
J124913-055919 | 12 49 13.86 | –05 59 19.1 | 2.247 | 13 | UVES
J124924-023339 | 12 49 24.86 | –02 33 39.7 | 2.12 | 11 | UVES
J131011+460124 | 13 10 11.61 | 46 01 24.50 | 2.133 | 12 | KODIAQ
J133335+164903 | 13 33 35.78 | 16 49 03.9 | 2.089 | 52 | UVES
J134427-103541 | 13 44 27.06 | –10 35 41.7 | 2.134 | 47 | UVES
J141719+413237 | 14 17 19.23 | 41 32 37.00 | 2.024 | 12 | KODIAQ
J141906+592312 | 14 19 06.31 | 59 23 12.20 | 2.321 | 18 | KODIAQ
J144653+011355 | 14 46 53.05 | 01 13 55.9 | 2.216 | 17 | UVES
J145102-232930 | 14 51 02.49 | –23 29 30.9 | 2.215 | 53 | UVES
J152156+520238 | 15 21 56.48 | 52 02 38.40 | 2.208 | 36 | KODIAQ
J162645+642655 | 16 26 45.69 | 64 26 55.20 | 2.32 | 20 | KODIAQ
J212329-005052 | 21 23 29.46 | –00 50 52.9 | 2.262 | 28 | KODIAQ
J221531-174408 | 22 15 31.65 | –17 44 08.2 | 2.217 | 19 | UVES
J222756-224302 | 22 27 56.92 | –22 43 02.5 | 1.891 | 32 | UVES
J225719-100104 | 22 57 19.04 | –10 01 04.7 | 2.08 | 12 | UVES
J231324+003444 | 23 13 24.45 | 00 34 44.50 | 2.083 | 11 | KODIAQ
J234023-005327 | 23 40 23.66 | –00 53 27.0 | 2.085 | 18 | KODIAQ

##  2 Data and analysis

In this study, the sample of 35 publicly available quasar spectra were used
for analysis. The first part of the sample consists of 13 QSO spectra from
KODIAQ (Keck Observatory Database of Ionized Absorption toward Quasars) survey
(Lehner et al., 2014; O’Meara et al., 2015, 2017), which is a repository of
reduced, continuum normalized spectra obtained by the High-Resolution Echelle
Spectrometer (HIRES) on the Keck Telescope. These spectra were observed
between years 1995–2012 and subsequently uniformly reduced and continuum
fitted by eye by KODIAQ team (O’Meara et al., 2015). The second part of the
sample contains of 22 QSO spectra obtained by the Ultraviolet and Visual
Echelle Spectrograph (UVES) on the VLT/ESO (Murphy et al., 2019). The UVES
Spectral Quasar Absorption Database comprises fully reduced, continuum-fitted
high-resolution spectra of quasars in the redshift range $0<z<5$. From the
whole dataset we selected only spectra that meet the following criteria:

  1. 1.

the QSO spectrum covers at least $\Delta z$ = 0.1 of at least one considered
redshift bin,

  2. 2.

the spectrum does not contain large gaps that would not allow the correct
identification of metal lines,

  3. 3.

the signal-to-noise (S/N) ratio of the spectrum is higher than 10 in the
studied spectral region.

The list of 35 QSOs whose spectra were used in this study with their basic
characteristics is presented in Tab. 1. The redshift coverage of the analysed
QSO spectra is shown in Fig. 1. We studied the effects of continuum
misplacement on the achieved results and showed its effect is negligible in
the case of our data. A detailed description of the analysis can be found in
Appendix A.

Figure 1: The redshift coverage of the dataset used in this study.

###  2.1 Voigt profile fitting

In this study, the VPFIT code written by R. F. Carswell and J. K. Webb
(Carswell & Webb, 2014) was used. The QSO spectra were analyzed within the
following parameter space: $b=1-300\,\rm km\,s^{-1}$ and
$\log{N_{\text{H\,{i}}}}=11.5-16.0$. The rest-frame wavelengths 1050 – 1180 Å
inside the Ly-$\alpha$ forest were chosen for the Voigt profile fitting to
avoid proximity effects, i.e. a region affected by the local QSO radiation
(exclude the Ly-$\beta$ and O VI $\lambda 1035$ emission lines) rather than
the metagalactic ultraviolet background (UVB; Palanque-Delabrouille et al.,
2013; Walther et al., 2018; Hiss et al., 2018).

The Ly-$\alpha$ absorbers corresponding to the damped Ly-$\alpha$ (DLA)
systems with $\log{N_{\text{H\,{i}}}}\gtrsim 20$ and also sub-damped
Ly-$\alpha$ (sub-DLA) systems were identified by eye and excluded from the
analysis. Each spectrum was also visually observed for bad points and gaps,
which were subsequently masked and cubically interpolated. Relative errors of
1% were given to such spectral regions, so the Voigt profile fitting procedure
was not influenced (Hiss et al., 2018). An example of the decomposition of the
Ly-$\alpha$ forest into individual absorption lines in the spectrum of QSO
J122824+312837 is shown in Fig. 2.

Figure 2: Resulting fit of the Ly-$\alpha$ forest in the spectrum of QSO
J122824+312837. The original spectrum (grey line) is well described by the
superposition of Voigt profiles fitted using the VPFIT code (red line).

###  2.2 Metal lines rejection

Quasar spectra also contain various metal absorption lines, which can affect
the obtained results. They are usually related to the strong H i absorption,
and can be identified via associations with other ionic metal lines. From this
reason, we identified DLA and sub-DLA systems and determined their redshifts
with help of the associated C iv, Si iv, Mg ii metal lines redward of the
Ly-$\alpha$ emission peak. It is worth noting that in the case of metal
absorption lines not related with DLA systems (mostly within the data
coverage), we used the doublets of Si iv, C iv and Mg ii to determine the
redshift of metal absorption systems. If the redshifts were known, the other
metal lines (see Tab. 2) were determined based on their characteristic
$\Delta\lambda$.

After the identification, the procedure for the metal line rejection in both
cases can be, in general, described as follows: firstly, we fitted the studied
intervals unmasked of identified lines to minimise the impact of possible
adverse effects of masking on the Voigt profile fitting. After this step, we
excluded the previously identified metal lines from the VPFIT output so that
they do not affect the results of our analysis. Example of the metal lines
rejection procedure based on the sub-DLA system at $z\approx 1.839$ in the
spectrum of QSO J110325-264515 is shown in Fig. 3.

As was described by Hiss et al. (2018), we also tested if the remaining
absorbers below the lower envelope of the $b-N_{\text{H\,{\sc I}}}$
distribution belong to the metal lines listed in Tab. 2. We would like to note
that if the absorber candidate was identified only as a single or doubtful
doublet feature, we did not consider it as a metal line.

Figure 3:  Example of the metal lines rejection procedure based on the sub-DLA
system at $z\approx 1.839$ in the spectrum of QSO J110325-264515 (upper
panel). The associated C iv 1548, 1550 Å spectral lines redward of the
Ly-$\alpha$ emission peak (bottom left panel) were used for more precise
determination of the sub-DLA redshift. The Si ii 1190, 1193 Å spectral lines
(bottom right panel) are related to the strong H i absorption depicted in
upper panel. The red dashed lines represent the continuum level. Table 2: List
of identified metal lines with their oscillator strength $f$.

Absorber |  $\lambda_{\text{rest}}$ [Å] | $f$ | Reference
---|---|---|---
O vi | 1031.9261 | 0.13290 | 1
C ii | 1036.3367 | 0.12310 | 1
O vi | 1037.6167 | 0.06609 | 1
N ii | 1083.9900 | 0.10310 | 1
Fe iii | 1122.5260 | 0.16200 | 2
Fe ii | 1144.9379 | 0.10600 | 3
Si ii | 1190.4158 | 0.25020 | 1
Si ii | 1193.2897 | 0.49910 | 1
N i | 1200.2233 | 0.08849 | 1
Si iii | 1206.5000 | 1.66000 | 1
N v | 1238.8210 | 0.15700 | 1
N v | 1242.8040 | 0.07823 | 1
Si ii | 1260.4221 | 1.00700 | 1
O i | 1302.1685 | 0.04887 | 1
Si ii | 1304.3702 | 0.09400 | 4
C ii | 1334.5323 | 0.12780 | 1
C ii* | 1335.7077 | 0.11490 | 1
Si iv | 1393.7550 | 0.52800 | 1
Si iv | 1402.7700 | 0.26200 | 1
Si ii | 1526.7066 | 0.12700 | 5
C iv | 1548.1950 | 0.19080 | 1
C iv | 1550.7700 | 0.09522 | 1
Al ii | 1670.7874 | 1.88000 | 1
Al iii | 1854.7164 | 0.53900 | 1
Al iii | 1862.7895 | 0.26800 | 1
Mg ii | 2796.3520 | 0.61230 | 6
Mg ii | 2803.5310 | 0.30540 | 6

References: (1) Morton (1991), (2) Prochaska et al. (2001), (3) Howk et al.
(2000), (4) Tripp et al. (1996), (5) Schectman et al. (1998), (6) Verner et
al. (1996).

###  2.3 The $b-N_{\text{H\,{i}}}$ distributions

The VPFIT outputs were used to generate the $b-N_{\text{H\,{\sc I}}}$
distributions, from which parameters in the $T-\rho$ relation could be
determined. The $b-N_{\text{H\,{\sc I}}}$ distribution of the parameterized
absorbers in the redshift range of $1.6\leq z<2.0$ is plotted in Fig. 4.
Firstly, we divided the results into two bins with the redshift range of
$\langle 1.6-1.8)$, and $\langle 1.8-2.0)$. Then, the lower limit of $b$ was
chosen to $10$ km s-1 mainly due to the lines with lower values could be
attributed to the ionic metal-line transitions (Hiss et al., 2018; Rauch et
al., 1997). The lines with $b>100$ km s-1 were also excluded, which is a
convention used for example by Rudie et al. (2012) and Hiss et al. (2018). The
reason for such a choice is that the turbulent broadening dominates over the
thermal one in these lines. In the case of the column density, the chosen
range of $12.5\leq\log{N_{\text{H\,{i}}}}\leq 14.5$ corresponds to the gas
density range for which the equation of state is well fitted by a power law
(Schaye et al., 1999).

Figure 4: The $b-N_{\text{H\,{\sc I}}}$ distribution of the parameterized
absorbers in the redshift range of $1.6\leq z<2.0$ with relative errors in the
Doppler widths and column densities less than 50%.

###  2.4 Narrow lines rejection

Although the aforementioned rejection procedure was applied, many narrow
absorption lines in blends and unidentified metal lines with $b\geq 10\,\rm
km\,s^{-1}$ still remain in the results. Due to this fact, we used the similar
iterative rejection algorithm as proposed by Hiss et al. (2018), which can be
described as follows: firstly, we divided the absorbers with $b<40\,\rm
kms^{-1}$ into eight $\log{N_{\text{H\,{i}}}}$ bins of the same size. Then,
the mean and standard deviation of the absorbers in every bin were calculated.
After this step, we excluded all points below $2\sigma$ of the mean. This
procedure is iterated until no points are excluded. After the last iteration,
we fitted the line to the $\log{b_{2\sigma}}$ values of each
$\log{N_{\text{H\,{i}}}}$ bin. Finally, if the position of this line is
determined, all absorbers below it are excluded.

Figure 5: $b-N_{\text{H\,{i}}}$ distributions are shown in left panels, where
the red solid and blue dashed lines correspond to the best cut-off fits and
the results of narrow-line rejection procedure, respectively. The overplotted
red circles correspond to the $\log{N_{\text{H\,{i}},0}}$. The right panels
show the corresponding PDFs $p(b_{0},\Gamma)$ , where the 68% and 95%
confidence levels are plotted by dark and light blue colors, respectively. The
black points correspond to the medians of the marginal distributions of
$b_{0}$ and $\Gamma$.

###  2.5 Fitting the cut-off in the $b-N_{\text{H\,{i}}}$ distribution

In this study, the iterative fitting procedure proposed by Schaye et al.
(1999) was used to determine the position of the thermal state sensitive cut-
off. It is based on the following equation

$\log{b_{\rm
th}}=\log{b_{0}}+(\Gamma-1)\log{N_{\text{H\,{i}}}/N_{\text{H\,{i}},0}}.$ (2)

where $b_{\rm th}$ corresponds to the thermal Doppler broadening, $b_{0}$ is
the minimal broadening value at column density $N_{\text{H\,{i}},0}$ and
$\Gamma$ is the index of this power-law relation. Firstly, we fit the Eq. (2)
to the points in the $b-N_{\text{H\,{i}}}$ distribution using the least
absolute deviations method. After this step, the mean absolute deviation is
calculated as

$|\delta\log{b}|=\frac{1}{n}\sum_{i=1}^{n}|\log{b_{i}}-\log{b_{\rm
th}(N_{\text{H\,{i}},\rm i})|}$ (3)

where $n$ is the sample size. Subsequently, all points for which the value of
the Doppler parameter $b$ fulfills the condition

$\log{b}>\log{b_{\rm th}}+|\delta\log{b}|$ (4)

are excluded in the next iteration. All these steps are repeated until no
points are more than one absolute mean deviation above the fit. The last step
is the exclusion of the points that are below the one absolute mean deviation
of the last fit. The remaining points are then used for the final fit, from
which $b_{0}$ and $(\Gamma-1)$ are determined.

###  2.6 Estimation of the $N_{\text{H\,{i}},0}$

As we mentioned before, the value of $N_{\text{H\,{i}},0}$ is just a
normalization. However, it is useful to choose its value that corresponds to
the column density of a typical absorber at the mean density of the IGM (Hiss
et al., 2018). Under the assumption of the local hydrostatic equilibrium in
the low-density cloud, the relationship between $N_{\text{H\,{i}}}$ and local
overdensity $\Delta=\rho/\overline{\rho}$ can be written as (Schaye, 2001;
Rudie et al., 2012)

$N_{\text{H\,{i}},0}\simeq
10^{13.23}\Delta^{3/2}\frac{T_{4}^{-0.22}}{\Gamma_{\text{ion},\text{H\,{i}}}}\left(\frac{1+z}{3.4}\right)^{9/2}[\text{cm}^{-2}],$
(5)

where $\Gamma_{\text{ion},\text{H\,{i}}}$ is the photoionization rate of H i
in units of $10^{-12}$ s-1 and $T_{4}$ is the temperature of the absorbing gas
in units of $10^{4}$ K.

In the case of simulations, the normalization factor
$N_{\text{H\,{i}},0}=N_{\text{H\,{i}}}(\Delta=1)$ was calculated using the Eq.
(5) for every single thermal model, which was used for calibration (see
below). We would like to note that we used the effective UVB
$\Gamma_{\text{ion},\text{H\,{i}}}$ =
$\Gamma_{\text{ion},\text{H\,{i}},\text{sim}}/A_{r}$ as the used simulations
were rescaled using the scaling factor $A_{r}$ (see Section 3.1). By including
this step, we reduced the variation of $N_{\text{H\,{i}},0}$ with thermal
parameters. In the case of simulation, the average value of
$\log{N_{\text{H\,{i}},0}}=(12.8247\pm 0.0872)\rm\,cm^{-2}$ over all of the
thermal models was used for the cut-off fitting procedure.

In the case of analysed redshift bins, we calculate the normalization factor
according to the equation (Hiss et al., 2018)

$\log{N_{\text{H\,{i}},0}(z)}=0.6225(1+z)+11.1068.$ (6)

The reason for this choice is that we used the THERMAL suite (see below),
which was also used in the aforementioned study to obtain the calibration
between $b_{0}-T_{0}$ and $\Gamma-\gamma$.

###  2.7 Cut-off fitting results

For the determination of the final cut-off parameters, we used the bootstrap
resampling method with replacement. We generate 2 000 datasets by
bootstrapping of the cut-off fitting procedure, from which the probability
distribution functions (PDFs) $p(b_{0},\Gamma)$ were obtained. The medians of
these distributions was used as the best estimates of the $b_{0}$ and
$(\Gamma-1)$ parameters. We would like to note that the uncertainties
correspond to the 16th and 84th percentiles of the probability distribution
functions. The resulting $b-N_{\text{H\,{i}}}$ distributions together with the
kernel density estimation of $p(b_{0},\Gamma)$ for both studied redshift bins
($1.6\leq z<1.8$ and $1.8\leq z<2.0$) are plotted in Fig. 5.

##  3 Simulations

To obtained more accurate values of the temperature $T_{0}$ and power-law
index $\gamma$, Bolton et al. (2014) suggested that calibration of the cut-off
fitting measurements based on the simulations is needed. To reach this goal,
we used the set of simulated skewers at redshift $z=1.8$ from the THERMAL
(Thermal History and Evolution in Reionization Models of Absorption Lines)
111thermal.joseonorbe.com suite (for more details see Oñorbe et al., 2017).
This dataset includes skewers from 88 Nyx hydrodynamical simulations with
different combinations of underlying thermal parameters $T_{0}$, $\gamma$ and
pressure smoothing scale $\lambda_{P}$ (see details in Oñorbe et al., 2017;
Hiss et al., 2018) on a box size $L_{\rm box}=20$ Mpc/$h$ and 10243 cells. The
cosmological parameters used in the simulations were based on the results of
the Planck mission (Planck Collaboration et al., 2014):
$\Omega_{\Lambda}=0.6808$, $\Omega_{m}=0.3192$, $\sigma_{8}=0.826$,
$\Omega_{b}=0.04965$, $n_{s}=0.9655$, and $h=0.6704$. The parameters $T_{0}$
and $\gamma$ were extracted from the simulations by fitting the $T-\rho$
relation power-law to the distribution of gas cells as described in Lukić et
al. (2015). In order to determine $\lambda_{\rm P}$, the cut-off in the power
spectrum of the real-space Lyman-$\alpha$ flux was fitted.

###  3.1 Skewer generation

Firstly, random Ly-$\alpha$ optical depth $(\tau)$ skewers were transformed
into the corresponding flux skewer $F$, i.e. a transmission spectrum along the
line of sight, according to the equation $F=F_{c}\exp{(-A_{r}\tau(\lambda))}$,
where continuum flux $F_{c}$ was set up to unity and the scaling factor
$A_{r}$ allowed the match the lines of sight to observed mean flux values.

Its value was determined by comparing the mean flux of the simulations with
observational one, in this case with the value that match the mean flux
evolution shown in Oñorbe et al. (2017) based on precise measurements by Fan
et al. (2006); Becker et al. (2007); Kirkman et al. (2007); Faucher-Giguère et
al. (2008); Becker & Bolton (2013). We would like to note that we did not
consider the impact of uncertainties on the scaling factor value. However, the
discussion of mean flux rescaling in the models is presented in Appendix B.

###  3.2 Thermal parameter grid

In this study, we used 23 skewers with a different combination of thermal
parameters and pressure smoothing scale. The combinations of parameters
$T_{0},\gamma$ and $\lambda_{\rm P}$ are depicted in the Fig. 6. The grid of
thermal parameters covers intervals: $6\,200<T_{0}<21\,200$ [K] and
$1.06<\gamma<1.91$. In the case of the pressure smoothing scale, we chose
range $46<\lambda_{\rm P}<121$ [kpc].

Figure 6: Combinations of parameters $T_{0}$, $\gamma$ and $\lambda_{\rm P}$
used for the calibration procedure.

###  3.3 Modeling Noise, Resolution and Voigt Profile Fitting

The mock spectra for calibration were prepared by adding the effects of
resolution and noise to the simulated skewers. The quantities of both effects
were selected to match the used spectra. We mimicked the instrumental
resolution by convolving the skewers with a Gaussian with FWHM = 6 km s-1 and
rebinning to 3 km s-1 pixels afterward. In the case of the S/N ratio, we
selected the value of 20. We would like to note that the metal contamination
was not considered in this case.

When the mock spectra were prepared, we used the same Voigt profile fitting
procedure as was described in Section 2.1. An example of the decomposition of
the Ly-$\alpha$ forest from the simulation into individual absorption lines
using VPFIT code is shown in Fig. 7. In this case, we also used the narrow
line rejection procedure even if the effect of the rejected lines on the
results is neglected.

Figure 7: Resulting Voigt profile fit of the Ly-$\alpha$ forest from the
simulation with the thermal parameters $\log(T_{0})=3.97$ [K], $\gamma=1.61$,
and smoothing scale $\lambda_{\rm P}=74$ kpc. The red line is the spectrum
fitted by the VPFIT code and the underlying grey solid line is the simulated
spectrum.

##  4 Determination of the coefficients in the $T-\rho$ relation

In this section, we deal with the relationship between the cut-off parameters
$b_{0}$, $(\Gamma-1)$ and the thermal parameters $T_{0}$, $\gamma$. Firstly,
we rewrote the Eq. (1) as

$\log{T}=\log{T_{0}}+(\gamma-1)\log{\left(\frac{\rho}{\overline{\rho}}\right)},$
(7)

where $\rho$ is the baryon density and $\overline{\rho}$ is its mean value
($\rho=\overline{\rho}(1+\delta)$). Schaye et al. (1999) showed that the
relations between the density vs temperature and column density vs Doppler
parameter of the absorbers near the cut-off can be also fitted by the power
laws as

$\log{\left(\frac{\rho}{\overline{\rho}}\right)}=A+B\log{\left(\frac{N_{\text{H\,{i}}}}{N_{\text{H\,{i}},0}}\right)},$
(8) $\log{T}=C+D\log{b_{\rm th}},$ (9)

where the coefficients are given by following equations

$\log{b_{0}}=\frac{1}{D}[\log{T_{0}}-C+(\gamma-1)A],$ (10)
$(\Gamma-1)=\frac{B}{D}(\gamma-1).$ (11)

It is worth noting that the Eq. (2) characterizes the minimal broadening at a
given $N_{\text{H\,{i}}}$, where the absorbers are thermally broadened
(therefore $b_{\rm th}$ in the relation). The chosen value of
$N_{\text{H\,{i}},0}$ represents the column density of a cloud with mean
density, therefore $A=0$, and the dependency on $\gamma$ disappears from the
Eq. (10) (Hiss et al., 2018). Based on these assumptions, we can rewrite the
Eqs. (10) and (11) as

$\log{T_{0}}=C+D\log{b_{0}},$ (12) $(\gamma-1)=\kappa(\Gamma-1),$ (13)

where $\kappa=D/B$.

###  4.1 Calibration Using Simulations

To generate the calibrations between $b_{0}-T_{0}$ and $\Gamma-\gamma$, we
used 23 selected thermal models at the redshift $z=1.8$. In this case, the
cut-off fitting algorithm was used on simulated $b-N_{\text{H\,{i}}}$
distributions, each of which was constructed using 120 mock spectra from all
used thermal models. The results are shown in Fig. 8, where depicted points
correspond to the median values of $b_{0}$ and $\Gamma-1$ from 2 000 datasets
generated by bootstrapping. This is the same approach as we used for the
observational data.

To include the additional effects, such as various values of $\lambda_{\rm P}$
and using of the same value of the $\log{N_{\text{H\,{i}}}}$ for all models,
we fitted Eqs. (12) and (13) to the 2 000 bootstrap realizations of the points
in the $\log{(T_{0})}-\log{(b_{o})}$ and $(\gamma-1)-(\Gamma-1)$ diagrams (see
Fig. 8) with replacements. In the $\log{(b_{0})}$ vs $\log{(T_{0})}$ diagram,
we also included the case when the value of $b_{0}$ is caused purely by
thermal broadening, i.e. $b_{0}=\sqrt{2k_{\rm B}T_{0}/m_{p}}$, where $k_{\rm
B}$ is the Boltzmann constant and $m_{p}$ is the proton mass.

The best fit values of the coefficients $C$, $D$ and $\kappa$ correspond to
the medians of bootstrap distributions. Using this approach, we obtained the
values: $C=1.39^{+0.20}_{-0.24}$, $D=2.33^{+0.20}_{-0.17}$, and
$\kappa=3.32^{+0.17}_{-0.17}$, where the uncertainties correspond to the 16th
and 84th percentiles of the $p(D,C)$ and $p(\kappa)$.

Figure 8: Calibrations of the $\log{b_{0}}$ vs $\log{T_{0}}$ (upper panel) and
$(\Gamma-1)$ vs $(\gamma-1)$ (bottom panel) relations. The black lines
correspond to the best linear fits to the points. The case, when the value of
$b_{0}$ corresponds to pure thermal broadening is depicted by the green dash
line.

##  5 Results and discussion

The $b-N_{\text{H\,{i}}}$ distributions with the corresponding cut-off fits
and PDFs $p(b_{0},\Gamma)$ for both studied redshift bins ($1.6\leq z<1.8$ and
$1.8\leq z<2.0$) are depicted in Fig. 5. The PDFs indicate the anticorrelation
between $b_{0}$ and $(\Gamma-1)$.

To include uncertainties which arose during the individual steps implemented
in the analysis, we used the approach described by Hiss et al. (2018). We
combined the 2 000 bootstrapped $b_{0}$ and $(\Gamma-1)$ pairs in
$p(b_{0},\Gamma)$ with every point in the boostrapped calibration PDFs
$p(D,C)$ and $p(\kappa)$ from simulations using the Eqs. (12) and (13). The
resulting PDFs $p(T_{0},\gamma)$ for both studied redshift intervals are
depicted in Fig. 9. The medians of bootstrap distributions were used as the
best estimates of the $T_{0}$ and $\gamma$. The derived parameter values for
the redshift intervals $1.6\leq z<1.8$ and $1.8\leq z<2.0$ are listed in Tab
3. The uncertainties correspond to the 16th and 84th percentiles of the
probability distribution function $p(T_{0},\gamma)$. The comparison of the
results obtained in this study with previously published ones is shown in Fig.
10.

Figure 9: Resulting PDFs $p(T_{0},\gamma)$ for the redshift range $1.6\leq
z<1.8$ (left panel) and $1.8\leq z<2.0$ (right panel). The 68% and 95%
confidence levels are plotted by dark and light blue color, respectively. The
black points correspond to the median of the marginal distributions of $T_{0}$
and $\gamma$.

Schaye et al. (1999) showed that the effect of the S/N ratio on the results is
negligible if its value is higher than $\sim 20$. However, the median of this
parameter for the spectra used in this study is equal to 18. We tested the
possibility that the lower value of $T_{0}$ and higher value of $\gamma$ for
mean value of redshift $z=1.8$ could be affected by the lower S/N ratio value,
but our analysis showed that this effect is possible to neglect (see Appendix
C).

We also considered the possibility that the obtained values of the calibration
parameters could be influenced by using the calibration for the redshift
$z=1.8$, whereas the centers of the redshift bins analysed in this study are
1.7 and 1.9. Based on the dependencies of the calibration parameters $C$, $D$,
and $\kappa$ on the redshift presented by Hiss et al. (2018) (Fig. 14
therein), we can conclude that changes of their values for the redshift
differences $\Delta z=0.1$ are significantly lower than their uncertainties
determined in this study, and this effect can be considered negligible.

The extended uncertainties of the derived parameters for the redshift range
$1.6\leq z<1.8$ originate from the multimodal distribution $p(b_{0},\Gamma)$
(see Figs. 5). It is worth noting that the simulated $b-N_{\text{H\,{i}}}$
distributions led to the unimodal solutions. Therefore, this behavior can be
attributed to some systematics related with the used spectra, for which the
cut-off fitting led to the artifacts. We assume that the effect is caused by
the number of absorbers ($\sim 550$) used for this redshift bin.

Table 3: Summary of the derived parameters for the studied redshift bins. Redshift | $b_{0}$ | $T_{0}$ | $\Gamma$ | $\gamma$
---|---|---|---|---
range |  $\rm[km\,s^{-1}$] |  $[10^{3}\rm\,K$] |  |
$\langle 1.6,1.8)$ | $12.26^{+5.67}_{-1.49}$ | $8.69^{+12.05}_{-2.48}$ | $1.17^{+0.05}_{-0.06}$ | $1.55^{+0.20}_{-0.18}$
$\langle 1.8,2.0)$ | $11.78^{+0.29}_{-0.83}$ | $7.51^{+0.81}_{-1.09}$ | $1.22^{+0.03}_{-0.02}$ | $1.72^{+0.11}_{-0.08}$

###  5.1 Comparison with previous studies

By comparing the values of $T_{0}$ and $\gamma$ obtained in this study with
previously published ones (Fig. 10), it can be concluded that even there are
some discrepancies, the results correspond to each other. An exception is the
study by Boera et al. (2014), in which the temperature at the mean density
assuming the value of $\gamma\sim 1.3$ is significantly higher for both mean
redshifts ($z=1.63$ and $1.82$) even compared to simulations. The differences
can be attributed to the different methods used to estimate the values of
$T_{0}$ and $\gamma$.

Only Schaye et al. (2000) studied the IGM in the redshift interval $1.85\leq
z\leq 2.09$ using the Voigt profile decomposition of the Ly-$\alpha$ forest
into the set of individual absorption lines, as was used in this work. For the
median value of the redshift $z=1.96$, the authors determined value of
$T_{0}\sim 11\,000$ K and $\gamma=1.4$, which differ from the ones obtained in
our study for the redshift range $\langle 1.8,2.0)$. It is worth noting that
Schaye et al. (2000) used the spectrum of only one quasar (Q1100-264), the
different value of $\log{N_{\text{H\,{i}},0}}=13$, and their median redshift
value was slightly different from the one in this study ($z=1.9$).

###  5.2 Comparison with models of the IGM thermal evolution

We compared the values of $T_{0}$ and $\gamma$ obtained in this study with
results of various models of the IGM thermal evolution (see Fig. 10). Namely,
the non-equilibrium simulation by Puchwein et al. (2015), the semi-analytical
models by Upton Sanderbeck et al. (2016) where the He II reionization has and
has not occurred, and hydrodynamical simulation by Oñorbe et al. (2017). The
semi-analytical models by Upton Sanderbeck et al. (2016) were used only for
qualitative comparison as these models depend on a number of parameters such
as the ionizing background spectral index, the quasar spectral index, duration
of the reionization, clumping etc. In addition, these models are not
independent of observational data because the authors used the temperature
measurements of Becker et al. (2011) and Boera et al. (2014) for constraining
the model curves.

Figure 10: Comparison of the results obtained in this study with previously
published ones and various models. The black and red curves correspond to the
semi-analytical models by Upton Sanderbeck et al. (2016) where the He II
reionization has and has not occurred, respectively. The blue and violet
curves correspond to the non-equilibrium simulation by Puchwein et al. (2015)
and the hydrodynamical simulation by Oñorbe et al. (2017), respectively.

The comparison shows that the obtained temperatures at mean density for both
redshift bins $\langle 1.6,1.8)$ and $\langle 1.8,2.0)$ are slightly lower
than predicted by various models. In the case of redshift bin $\langle
1.6,1.8)$, the value of $T_{0}$ corresponds to the simulations within its
error, although the uncertainty is considerably large.

The same conclusion can also be stated in the case of the power-law index, for
which its value for the redshift bin $\langle 1.6,1.8)$ matches the non-
equilibrium simulation of Puchwein et al. (2015), and within the error
corresponds to all presented models. On the other hand, for the redshift bin
$\langle 1.8,2.0)$, we obtained the value of the power-law index which is
higher than the expected one ($\gamma\sim 1.6$) resulting from the balance of
photoheating and adiabatic cooling after all reionization events (Hui &
Gnedin, 1997; Upton Sanderbeck et al., 2016).

##  6 Conclusions

In this study, the sample of 35 quasar spectra was used for inferring
parameters of the temperature-density relation in the IGM for the redshift
interval $\langle 1.6,2.0)$. For this purpose, we used the approach that
treats the Ly-$\alpha$ forest as a superposition of discrete absorption
profiles. Based on the fitting, the $b-N_{\text{H\,{i}}}$ distributions were
constructed, from which the position of thermal state sensitive cut-off was
determined. We also calibrate our measurements using 23 combinations of
thermal parameters and pressure smoothing scale from the THERMAL suite of
hydrodynamical simulations. The main results can be summarized as follows:

  * •

The obtained values correspond to the decreasing trend of the temperature at
mean density predicted by the various models toward the lower redshifts.

  * •

For the redshift interval $\langle 1.6,1.8)$, we obtained the power law index
$\gamma=1.55$. The value is close to $1.6$, as expected by various models, and
results from the balance of photoheating with adiabatic cooling after all
reionization events. However, the determined value $\gamma=1.72$ is higher for
the redshift interval $\langle 1.8,2.0)$.

  * •

Our determinations of $T_{0}$ and $\gamma$ at the lower redshift interval
$\langle 1.6,2.0)$ indicate the completion of the reheating process associated
to He ii reionization.

## Acknowledgements

This research is based on the use of the Keck Observatory Archive (KOA), which
is operated by the W. M. Keck Observatory and the NASA Exoplanet Science
Institute (NExScI), under contract with the National Aeronautics and Space
Administration, and also on data products created from observations collected
at the European Organisation for Astronomical Research in the Southern
Hemisphere. The authors would also like to thank Jose Oñorbe for fruitful
discussion.

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Becker & Bolton (2013) Becker G. D., Bolton J. S., 2013, MNRAS, 436, 1023
  * Becker et al. (2007) Becker G. D., Rauch M., Sargent W. L. W., 2007, ApJ, 662, 72
  * Becker et al. (2011) Becker G. D., Bolton J. S., Haehnelt M. G., Sargent W. L. W., 2011, MNRAS, 410, 1096
  * Boera et al. (2014) Boera E., Murphy M. T., Becker G. D., Bolton J. S., 2014, MNRAS, 441, 1916
  * Bolton et al. (2008) Bolton J. S., Viel M., Kim T. S., Haehnelt M. G., Carswell R. F., 2008, MNRAS, 386, 1131
  * Bolton et al. (2014) Bolton J. S., Becker G. D., Haehnelt M. G., Viel M., 2014, MNRAS, 438, 2499
  * Carswell & Webb (2014) Carswell R. F., Webb J. K., 2014, VPFIT: Voigt profile fitting program (ascl:1408.015)
  * Fan et al. (2006) Fan X., et al., 2006, AJ, 132, 117
  * Faucher-Giguère et al. (2008) Faucher-Giguère C.-A., Prochaska J. X., Lidz A., Hernquist L., Zaldarriaga M., 2008, ApJ, 681, 831
  * Hiss et al. (2018) Hiss H., Walther M., Hennawi J. F., Oñorbe J., O’Meara J. M., Rorai A., Lukić Z., 2018, The Astrophysical Journal, 865, 42
  * Howk et al. (2000) Howk J. C., Sembach K. R., Roth K. C., Kruk J. W., 2000, ApJ, 544, 867
  * Hui & Gnedin (1997) Hui L., Gnedin N. Y., 1997, MNRAS, 292, 27
  * Kim et al. (2007) Kim T. S., Bolton J. S., Viel M., Haehnelt M. G., Carswell R. F., 2007, MNRAS, 382, 1657
  * Kirkman et al. (2005) Kirkman D., et al., 2005, MNRAS, 360, 1373
  * Kirkman et al. (2007) Kirkman D., Tytler D., Lubin D., Charlton J., 2007, MNRAS, 376, 1227
  * Lehner et al. (2014) Lehner N., O’Meara J. M., Fox A. J., Howk J. C., Prochaska J. X., Burns V., Armstrong A. A., 2014, The Astrophysical Journal, 788, 119
  * Lidz et al. (2010) Lidz A., Faucher-Giguère C.-A., Dall’Aglio A., McQuinn M., Fechner C., Zaldarriaga M., Hernquist L., Dutta S., 2010, ApJ, 718, 199
  * Lukić et al. (2015) Lukić Z., Stark C. W., Nugent P., White M., Meiksin A. A., Almgren A., 2015, MNRAS, 446, 3697
  * McDonald et al. (2000) McDonald P., Miralda-Escudé J., Rauch M., Sargent W. L. W., Barlow T. A., Cen R., Ostriker J. P., 2000, ApJ, 543, 1
  * Morton (1991) Morton D. C., 1991, ApJS, 77, 119
  * Murphy et al. (2019) Murphy M. T., Kacprzak G. G., Savorgnan G. A. D., Carswell R. F., 2019, MNRAS, 482, 3458
  * Oñorbe et al. (2017) Oñorbe J., Hennawi J. F., Lukić Z., 2017, ApJ, 837, 106
  * O’Meara et al. (2015) O’Meara J. M., et al., 2015, The Astronomical Journal, 150, 111
  * O’Meara et al. (2017) O’Meara J. M., Lehner N., Howk J. C., Prochaska J. X., Fox A. J., Peeples M. S., Tumlinson J., O’Shea B. W., 2017, The Astronomical Journal, 154, 114
  * Padmanabhan (2002) Padmanabhan T., 2002, Theoretical Astrophysics, Volume III: Galaxies and Cosmology. Cambridge University Press
  * Palanque-Delabrouille et al. (2013) Palanque-Delabrouille N., et al., 2013, Astronomy & Astrophysics, 559, A85
  * Planck Collaboration et al. (2014) Planck Collaboration et al., 2014, A&A, 571
  * Prochaska et al. (2001) Prochaska J. X., et al., 2001, ApJS, 137, 21
  * Puchwein et al. (2015) Puchwein E., Bolton J. S., Haehnelt M. G., Madau P., Becker G. D., Haardt F., 2015, MNRAS, 450, 4081
  * Rauch et al. (1997) Rauch M., et al., 1997, Astrophysical Journal, 489, 7
  * Ricotti et al. (2000) Ricotti M., Gnedin N. Y., Shull J. M., 2000, ApJ, 534, 41
  * Rorai et al. (2018) Rorai A., Carswell R. F., Haehnelt M. G., Becker G. D., Bolton J. S., Murphy M. T., 2018, MNRAS, 474, 2871
  * Rudie et al. (2012) Rudie G. C., Steidel C. C., Pettini M., 2012, The Astrophysical Journal Letters,, 757, L30
  * Schaye (2001) Schaye J., 2001, The Astrophysical Journal, 559, 507
  * Schaye et al. (1999) Schaye J., Theuns T., Leonard A., Efstathiou G., 1999, MNRAS, 310, 57
  * Schaye et al. (2000) Schaye J., Theuns T., Rauch M., Efstathiou G., Sargent W. L. W., 2000, MNRAS, 318, 817
  * Schectman et al. (1998) Schectman R. M., Povolny H. S., Curtis L. J., 1998, ApJ, 504, 921
  * Stoehr et al. (2008) Stoehr F., et al., 2008, in Argyle R. W., Bunclark P. S., Lewis J. R., eds, Astronomical Society of the Pacific Conference Series Vol. 394, Astronomical Data Analysis Software and Systems XVII. p. 505
  * Theuns et al. (2002) Theuns T., Schaye J., Zaroubi S., Kim T.-S., Tzanavaris P., Carswell B., 2002, ApJ, 567, L103
  * Tripp et al. (1996) Tripp T. M., Lu L., Savage B. D., 1996, ApJS, 102, 239
  * Upton Sanderbeck et al. (2016) Upton Sanderbeck P. R., D’Aloisio A., McQuinn M. J., 2016, MNRAS, 460, 1885
  * Verner et al. (1996) Verner D. A., Verner E. M., Ferland G. J., 1996, Atomic Data and Nuclear Data Tables, 64, 1
  * Viel et al. (2009) Viel M., Bolton J. S., Haehnelt M. G., 2009, MNRAS, 399, L39
  * Walther et al. (2018) Walther M., Hennawi J. F., Hiss H., Oñorbe J., Lee K.-G., Rorai A., O’Meara J., 2018, The Astrophysical Journal,, 852, 22
  * Walther et al. (2019) Walther M., Oñorbe J., Hennawi J. F., Lukić Z., 2019, ApJ, 872, 13
  * Zaldarriaga et al. (2001) Zaldarriaga M., Hui L., Tegmark M., 2001, ApJ, 557, 519

##  Appendix A Effects of continuum misplacement

In this section, we discuss the effect of continuum misplacement in used data.
In the case of high-resolution, high S/N spectra, the continuum is fitted
locally by connecting apparent absorption-free spectral intervals (McDonald et
al., 2000; Kim et al., 2007). As was noted in Kim et al. (2007), this method
is applicable for the redshift $1.5<z_{\text{em}}<3.5$, where the continuum
placement statistical uncertainty is of the order of a few percent.

As the continuum placement affects the corresponding optical depth of the
spectral line, and consequently the line profile parameters, we address the
effect by calculating the optical depth $\tau$ in the line center as
(Padmanabhan, 2002)

$\tau_{\text{lc}}({\text{H\,{i}}})=1.497\times
10^{-15}\frac{Nf\lambda_{0}}{b},$ (14)

where $\lambda_{0}=1215.67$ Å, $N\,[\rm cm^{-2}]$ and $b\,[\rm km\,s^{-1}]$.
For the calculation $\tau_{\text{lc}}({\text{H\,{i}}})$ using Eq. (14), we
used the typical value of the Doppler parameter in our sample $b=13\,\rm
km\,s^{-1}$ for spectral lines with various column densities. Then, we
calculated the flux at line center according to equation
$F_{\text{lc}}=F_{\text{c}}\exp{(-\tau_{\text{lc}}({\text{H\,{i}}}))}$, where
$F_{\text{c}}$ was set up to unity, and the value of $F_{\text{lc}}$ was
shifted to mimic the effect of misplacement of the continuum. We assumed that
our typical continuum uncertainty is $\sim$5%. The final step was to calculate
the corresponding column density by the reverse calculation with the fixed
value of Doppler parameter.

The results showed that for the continuum shift of 5%, the corresponding shift
of $\log{N_{\text{H\,{i}}}}$ is comparable to the VPFIT uncertainty at
$\log{N_{\text{H\,{i}}}}\approx 12.5$ and exceeds it at lower values of column
densities. From this result we can concluded that the effect of the continuum
misplacement is possible to neglect within the cut-off fitting range used in
this study.

##  Appendix B Mean flux rescaling in the models

As noted in Section 3.1, we rescaled the fluxes obtained from the simulations
in order to match the observed mean flux $\overline{F}$. In this work, we
chose the value corresponding the mean flux evolution in Oñorbe et al. (2017),
which is based on precise measurements by Fan et al. (2006); Becker et al.
(2007); Kirkman et al. (2007); Faucher-Giguère et al. (2008); Becker & Bolton
(2013). The $\overline{F}$ value used in the case of simulations corresponds
to the mean transmitted flux determined based on the KODIAQ/UVES spectra. The
comparison is plotted in Fig. 11 together with the values determined from the
studies of Kirkman et al. (2005) and Kim et al. (2007). We direct the reader
to the study of Hiss et al. (2018) for the discussion of the effects of mean
flux rescaling in the models.

Figure 11: The comparison of the mean flux evolutions from various studies
with the mean flux of the data used in this work.

##  Appendix C Effect of the S/N ratio on the results of cut-off procedure

To test the effects of the S/N ratio for values $\lesssim 20$ (Schaye et al.,
1999), we simulated two $b-N_{\text{H\,{i}}}$ distributions at $z=1.8$ with
different S/N ratios applied to lines of sight (see Fig. 12). The first one
corresponds to the highest S/N value of used spectra (83), and the second one
corresponds to the value, which was used in the case of simulations (20).
After application of the same cut-off fitting procedure as was used in the
case of spectra and simulations, we obtained the values
$b_{0}=14.18^{+0.30}_{-0.38}$, $(\Gamma-1)=0.18^{+0.01}_{-0.02}$ and
$b_{0}=14.31^{+0.32}_{-0.67}$, $(\Gamma-1)=0.17^{+0.01}_{-0.01}$ for the S/N =
20 and S/N = 83, respectively. Since the values of the parameters are same
within their errors, we can conclude that the effect of the S/N ratio is
negligible.

Figure 12: Simulated $b-N_{\text{H\,{i}}}$ distributions at $z=1.8$ with
various S/N ratio values applied to lines of sight. The distributions were
generated by Voigt profile fitting the same skewers as were used in the case
of calibration process. The red points depict the distribution with S/N = 83
which corresponds to the highest S/N ratio value of used spectra. The black
points correspond to the distribution with S/N = 20 characterizing the S/N
ratio used in the case of the simulations. The thermal parameters used in
these mock distributions are $\log(T_{0})=3.97$ [K], $\gamma=1.61$, and
smoothing scale $\lambda_{\rm P}=74$ kpc.


# Non-Blocking Batch A* (Technical Report)

Rishi Veerapaneni, Maxim Likhachev

###### Abstract

Heuristic search has traditionally relied on hand-crafted or programmatically
derived heuristics. Neural networks (NNs) are newer powerful tools which can
be used to learn complex mappings from states to cost-to-go heuristics.
However, their slow single inference time is a large overhead that can
substantially slow down planning time in optimized heuristic search
implementations. Several recent works have described ways to take advantage of
NN’s batch computations to decrease overhead in planning, while retaining
bounds on (sub)optimality. However, all these methods have used the NN
heuristic in a “blocking” manner while building up their batches, and have
ignored possible fast-to-compute admissible heuristics (e.g. existing
classically derived heuristics) that are usually available to use. We
introduce Non-Blocking Batch A* (NBBA*), a bounded suboptimal method which
lazily computes the NN heuristic in batches while allowing expansions informed
by a non-NN heuristic. We show how this subtle but important change can lead
to substantial reductions in expansions compared to the current blocking
alternative, and see that the performance is related to the information
difference between the batch computed NN and fast non-NN heuristic.

## Introduction

Heuristic search is a common deterministic motion planning class of algorithms
that plans paths from start and goal positions by searching over different
paths. Heuristic search algorithms use heuristics, along with other
mechanisms, to speed up the path planning process.

Modern machine learning techniques (e.g. reinforcement, imitation and
supervised learning) learn complex mappings from states to values that can be
viewed as cost-to-go heuristics. Recent work have started to use these
techniques to learn neural network heuristics on combinatorial search and
pathfinding problems which are better able to guide search compared to
existing heuristics (Agostinelli et al. 2019, 2021; Kim and An 2020; Jabbari
Arfaee, Zilles, and Holte 2011; Chrestien et al. 2021).

A major bottleneck of using neural networks (NNs) during search is that the
computational overhead to compute a single heuristic value is significantly
higher than non-neural network approaches. Neural network inferences are known
to be faster when computed on batches of inputs allowing them to take
advantage of their single instruction multiple data (SIMD) nature and
specialized SIMD hardware (i.e. GPUs) rather than computed many times
individually.

Recently, there have been several works that describe different ways for Batch
A*. Agostinelli et al. (2019) describes “Batch Weighted A*”, BWAS, which
expands K states from their OPEN queue. This unfortunately does not maintain
any suboptimality bounds with the inadmissible learnt heuristics. Greco et al.
(2022) introduces “K-Focal Search” which improves upon BWAS by utilizing a
fast admissible heuristic ($h_{fast}$) and a $w_{so}$ bounded suboptimal Focal
Search and expanding K nodes at once from the FOCAL queue that have f-values
$\leq\min_{n}g(n)+h_{fast}(n)$. The admissible $h_{fast}$ (which is usually a
classically derived heuristic, or the zero heuristic if none are available)
maintains bounded suboptimality in the Focal queue and enables batch
computation of the inadmissible $h_{NN}$ during each batch expansion. Their
method shows large runtime benefits by increasing the batch size B to decrease
the runtime overhead of $h_{NN}$. Another recent work Li et al. (2022)
describes a batch A* mechanism that maintains optimality, but contains several
restrictions on its use case as it requires an admissible learnt heuristic and
a problem where many states may share the same heuristic values.

Our method, Non-Blocking Batch A*, is closest to Greco et al. (2022) as we
utilize a $w_{so}$ bounded suboptimal search, an admissible fast heuristic
$h_{fast}$, and an inadmissible neural network heuristic $h_{NN}$. Our two key
differences are that instead of expanding K states at once, we lazily compute
the NN heuristic in batches and that we eliminate waiting by allowing the
search to expand states without requiring the NN heuristic to be computed. We
show that these are required to sustain reduction in nodes generated and
obtain speed-ups.

## Non-Blocking Batch A*

“K-Focal Search” is constructed using a NN function $h_{succ}()$ that takes in
a state $n$ and returns the $h$ values of all its successors $n^{\prime}$.
Thus, expanding K nodes at once allowed them to compute the heuristics of all
their successors at once and insert them into queues accordingly. This
requirement on the NN function is restrictive as the same NN cannot work with
different action spaces as the successor states and corresponding heuristic
values would change when using different action spaces. Additionally this
requirement generally increases the complexity of the learning problem as we
must map one parent state to many successor heuristic values. We instead
formulate our work with a NN heuristic function $h_{state}$ which returns the
heuristic value of one state which simplifies the learning problem and is
independent of the action space. Note we can easily construct $h_{succ}()$ by
using this single-state NN $h_{state}$ internally.

The K-Focal Search analogue using $h_{state}$ is to expand K nodes, collect
their successor states into a single large batch B, and pass the batch into
$h_{state}$ to generate the heuristic values. However, this does not maintain
a fixed batch size B as expanded states can have different numbers of valid
successors. An iterative approach can mitigate this issue by expanding states
sequentially, adding successor states onto a waitlist until a target batchsize
of B is approximately reached, and then calling $h_{state}$ and inserting the
nodes into the queues.

We show that this approach has two fundamental problem; 1) the waitlist blocks
“deeper” expansions, and 2) the search only uses the non-neural network fast
admissible heuristic ($h_{fast}$) to maintain a lower bound and never utilizes
it to guide search. Specifically, once a node $n$ is expanded, its successors
are placed in the waitlist until the waitlist is full, preventing any
successor from being expanded in the meanwhile. This prevents rapid expansions
of subsequent states, essentially forcing us away from a greedy heuristic
search and into a K-breadth greedy search. This construction also ignores the
fact that we may have a fairly informative $h_{fast}$ that we can use along
with $h_{state}$, instead of solely using $h_{state}$ for guiding search
efforts. Additionally, finding a solution path of length $L$ (not cost) will
now require a minimum of $L$ batches which can be a large problem if B and L
are large.

Our alternative removes this blocking behaviour by utilizing $h_{fast}$ along
with $h_{state}$, allowing our search algorithm to expand consecutive nodes
guided by $h_{fast}$ without needing to wait for the expensive batch
computation.

Algorithm 1 shows the full pseudo-code of our method. The main changes are
commented in blue. Specifically, line 18 inserts nodes into FOCAL based on
$h_{fast}$ instead of solely residing in the waitlist as in “K-Focal Search”.
Lines 19-22 call $h_{NN}$ on a batch of states once the batchsize has
surpassed $B$, and updates $h_{fast}$ to use the these values (as they can be
cached into $h_{fast}$ and accessed quickly). Lines 7-9 check if $h_{fast}$
has a new updated value (i.e. from an earlier $h_{NN}$ batch call) and updates
$n$’s priority and reinserts accordingly.

Algorithm 1 Non-Blocking Batch A*

Hyper-parameters: Suboptimality bound $w_{so}$, Heuristic weight $w_{h}$
Input: $n_{start}$, atGoal(), Batch NN heuristic $h_{NN}()$, Fast admissible
heuristic $h_{fast}()$
Output: Path from $n_{start}$ with suboptimality $\leq w_{so}$

1:OPEN = FOCAL = $\\{n_{start}\\}$, Waitlist = $\\{\\}$

2:while FOCAL $\neq\emptyset$ do

3: $n\leftarrow$ FOCAL.pop()

4: OPEN.remove($n$)

5: if atGoal($n$) then

6: return Solution backtracking from $n$

7:$\triangleright$ Lazily update $h_{focal}$ and reinsert if needed

8: if $n.h_{focal}\neq h_{fast}(n)$ then

9: $n.h_{focal}\leftarrow h_{fast}(n)$

10: AddToQueues($n$)

11: else$\triangleright$ Generate successors

12: for $n^{\prime}$ $\in$ succ($n$) do

13: $n^{\prime}.g\leftarrow n.g+cost(n,n^{\prime})$

14: $n^{\prime}.h_{open},n.h_{focal}\leftarrow h_{fast}(n^{\prime})$

15: AddToQueues($n$)

16:$\triangleright$ Update FOCAL

17: for all $n^{\prime}\in OPEN,\notin FOCAL$ do

18: if $n^{\prime}.F_{open}\leq w_{so}*n.F_{open}$ then

19: Waitlist.insert($n^{\prime}$)

20: FOCAL.insert($n^{\prime}$) $\triangleright$ Non-blocking

21: if Waitlist.size() $\geq B$ then

22:$\triangleright$ Compute batch heuristic and update $h_{fast}$

23: $h_{1:K}\leftarrow h_{NN}(n_{i}\in\text{Waitlist})$

24: $\forall{n_{i}\in\text{Waitlist}},h_{fast}(n_{i})\leftarrow h_{i}$

25: Waitlist $=\\{\\}$

26:return No solution

27:

28:procedure AddToQueues($n$):

29: $n.F_{open}\leftarrow n.g+n.h_{open}$

30: $n.F_{focal}\leftarrow n.g+w_{h}*n.h_{focal}$

31: OPEN.insert($n$)

32: if $n.F_{open}\leq w_{so}*\min_{n^{\prime}}n^{\prime}.F_{open}$ then

33: FOCAL.insert($n$)

Figure 1: Node expansions — We see the numbers of nodes expanded during search
by our method with various fidelity $h_{fast}$ (blue, orange, red) and compare
it against the blocking version (red) across different batch sizes. Note that
the blocking version does not depend on $h_{fast}$ in FOCAL so varying
$h_{fast}$ did not change performance much. The last column shows the results
of just running the noisy $h_{fast}$; we see from the orange and green
baselines struggle. We see that increasing batch size increases the number of
expansions as predicted, and that our non-blocking version generally
outperforms the blocking version.

## NBBA* Experimental Results

### Set-up

We show results on path-finding in a 2D Map using a Neural Network. Our
planning space is $(x,y,\theta)$ with the following three actions: 1) forward
in respect to the $\theta$ heading, 2) stationary turn left 90$\degree$, 3)
turn right stationary turn left 90$\degree$. Our goals are $(x_{g},y_{g})$
tuples and our point mass agent starts at positions $(x_{s},y_{s},0)$. We test
on 30 instances on a large 8192x8192 binary “sand-trap” obstacle map (5%
obstacles randomly generated) where all transitions in the map are valid and
unit cost except for transitions exiting sand-traps which cost 100. Figures
show the mean and standard deviation across the 30 instances. We run methods
with $w_{so}=w_{h}=2.5$ until they find the goal or expand a maximum of 1
million nodes. This set-up was arbitrarily used from an existing project, but
from spot-checking briefly on different parameters the trends from these
experiments seem to hold across different 2D grid world scenarios (and likely
to other domains as well).

Our method requires two heuristics, the base fast heuristic $h_{fast}$ and the
batch NN heuristic $h_{NN}$. To relate the difference between the quality of
the heuristics provided by $h_{fast}$ and $h_{NN}$, we make both $h_{NN}$
noisy versions of the Manhattan distance of $(x,y)$. Concretely, let
$h_{m}(s)$ by the true Manhattan distance heuristic. We define a “noisy”
version with $k\in[0,1]$ as of $h_{m_{k}}=h_{m}*N$ where $N\sim U[1-k,1]$. A
noisy $h_{m_{k}}$ with $k=0$ is a “strong” heuristic, perfectly noiseless and
equal to $h_{m}$, while increase $k$ makes $h_{m_{k}}$ increasingly weak and
noisy. Note that although $h_{m_{0}}$ has a clear meaning, $h_{m_{1}}$ does
not and is not the same as a perfectly uninformative heuristic.

Since our goal is to demonstrate the effectiveness of our non-blocking nature,
we do not train a neural network to output $h_{m}$. Instead we use a dummy
$h_{NN}$ which creates a random tensor input of size $(B,242)$ and passes it
into a small 3-layer neural network for a timing purposes, but actually
returns a noisy $h_{m_{0.01}}$. We compare our method against K-Focal which
uses $h_{fast}$ for OPEN queue and $h_{NN}$ in FOCAL, and implement the
iterative batch building variant for a clear comparison against our method. We
label this K-Focal variant as “blocking” in our experiments.

### Results

Figure 1 shows the number of nodes expanded during search by our method as
well as the blocking version with varying $h_{fast}$ and batch sizes. The
right most column shows the results of running search with just $h_{fast}$
where we see that a noisy $n$ level of $0.005$ runs fast (blue), $n=0.05$
struggles but still manages to find paths with an average of 500k expansions
(orange), while $n=0.5$ is not informative and results in a majority of
failures. We first observe that in general, increasing batch sizes increases
the number of nodes expanded, which makes sense as increasing the batch size
requires expanding more nodes until $h_{NN}$ is called again and the updated
heuristic values can be used. However we observe that our non-blocking version
outperforms the blocking version across most instances. If $h_{fast}$ gets
stronger, the performance gap increases as our non-blocking algorithm allows
successive expansions guided by $h_{fast}$.

Figure 2: Neural network overhead — We plot the total inference time for
computing the heuristic compared to the total runtime. We first observe that
increasing batchsize reduces runtime as expected. The exact ratio changes
based on the strength of $h_{fast}$. As $h_{fast}$ becomes stronger, each
$h_{NN}$ call becomes less useful and more of a relative overhead. Figure 3:
Runtime — We compare the runtime of blocking versus non-blocking methods. We
see that for larger batchsizes (e.g. 25, 125), our non-blocking variant offers
runtime benefits at the cost of increased variance due to the use of
$h_{fast}$. Note that both the green and red methods reach the maximum node
expansion and terminate early on the rightmost side (batchsize of 625), but
that there our orange non-blocking method performs much faster even though the
orange heuristic itself is pretty weak (see Figure 1).

Figure 2 highlights the ratio of inference time that $h_{NN}$ takes compared
to the overall search time. Increasing batch size decreases the ratio as
expected. However, we see that the absolute value depends on the strength of
$h_{fast}$. We are unclear why this exactly occurs but can offer some initial
thoughts. As $h_{fast}$ becomes stronger, the search advances deeper between
$h_{NN}$ calls, generating more unique children and populating the waitlist
faster (which we confirmed experimentally, although it does not explain all
the difference). Each $h_{NN}$ call occurs relatively more frequent, but at
the same time adds marginally less information than $h_{fast}$. We feel this
“marginally less” information has a side-effect effecting the runtime ratio
but are unable to find concrete cause and effect with our initial effort.

Finally, we show the aggregate runtime of NBBA* and the blocking K-Focal. At a
low batch size and higher $h_{fast}$ noise, the blocking behaviour can be
beneficial as it eliminates expanding extra nodes and reduces overall runtime.
However with larger batchsizes this blocking behaviour becomes a bottleneck
resulting in our method being faster on average.

## Conclusion

We show a general method for incorporating a learned heuristic $h_{NN}$ in
bounded suboptimal manner that leverages batch computation as well as the fast
but not as informative base heuristic $h_{fast}$. Non-Blocking Batch A*
(NBBA*) lazily creates batches for $h_{NN}$ while permitting expansions based
on $h_{fast}$. We show this can provide speed ups over the current state of
the art K-Focal Search depending on the strength of $h_{fast}$ and batch size.
This technical report has several limitations; namely that we do not conduct
robust experimentation across different domains and neural networks sizes.
Further analysis could also be done into determining how the strength of
$h_{fast}$ changes the ratio of time spent on $h_{NN}$’s inference time as
mentioned in the previous section. However we believe and hope that NBBA* will
be useful to future practitioners incorporating learned heuristics with
search. Acknowledgements:  Rishi V. would like to thank Shivam S. for their
helpful discussion in proofreading this report.

## References

  * Agostinelli et al. (2019) Agostinelli, F.; McAleer, S.; Shmakov, A.; and Baldi, P. 2019.  Solving the Rubik’s cube with deep reinforcement learning and search.  _Nature Machine Intelligence_ , 1–8.
  * Agostinelli et al. (2021) Agostinelli, F.; Shmakov, A.; McAleer, S.; Fox, R.; and Baldi, P. 2021.  A* Search Without Expansions: Learning Heuristic Functions with Deep Q-Networks.  _CoRR_ , abs/2102.04518.
  * Chrestien et al. (2021) Chrestien, L.; Pevný, T.; Komenda, A.; and Edelkamp, S. 2021.  Heuristic Search Planning with Deep Neural Networks using Imitation, Attention and Curriculum Learning.  _CoRR_ , abs/2112.01918.
  * Greco et al. (2022) Greco, M.; Toro, J.; Ulloa, C. H.; and Baier, J. A. 2022.  K-Focal Search for Slow Learned Heuristics (Extended Abstract).  In Chrpa, L.; and Saetti, A., eds., _Proceedings of the Fifteenth International Symposium on Combinatorial Search, SOCS 2022, Vienna, Austria, July 21-23, 2022_ , 279–281. AAAI Press.
  * Jabbari Arfaee, Zilles, and Holte (2011) Jabbari Arfaee, S.; Zilles, S.; and Holte, R. C. 2011.  Learning heuristic functions for large state spaces.  _Artificial Intelligence_ , 175(16): 2075–2098.
  * Kim and An (2020) Kim, S.; and An, B. 2020.  Learning Heuristic A: Efficient Graph Search using Neural Network.  In _2020 IEEE International Conference on Robotics and Automation (ICRA)_ , 9542–9547.
  * Li et al. (2022) Li, T.; Chen, R.; Mavrin, B.; Sturtevant, N. R.; Nadav, D.; and Felner, A. 2022\.  Optimal Search with Neural Networks: Challenges and Approaches.  In Chrpa, L.; and Saetti, A., eds., _Proceedings of the Fifteenth International Symposium on Combinatorial Search, SOCS 2022, Vienna, Austria, July 21-23, 2022_ , 109–117. AAAI Press.