# A consistent and robust measurement of the thermal state of the IGM at
$2\leq z\leq 4$ from a large sample of Ly$\alpha$ forest spectra: Evidence for
late and rapid HeII reionization

Prakash Gaikwad1,2, Raghunathan Srianand3, Martin G. Haehnelt1,2 and
Tirthankar Roy Choudhury4
1Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge,
CB3 0HA, UK
2Kavli Institute for Cosmology, University of Cambridge, Madingley Road,
Cambridge, CB3 0HA, UK
3Inter-University Centre for Astronomy and Astrophysics (IUCAA), Post Bag 4,
Pune 411007, India
4National Centre for Radio Astrophysics, Tata Institute of Fundamental
Research, Pune 411007, India E-mail<EMAIL_ADDRESS>
###### Abstract

We characterise the thermal state of the intergalactic medium (IGM) in ten
redshift bins in the redshift range $2\leq z\leq 4$ with a sample of 103 high
resolution, high S/N Ly$\alpha$ forest spectra using four different flux
distribution statistics. Our measurements are calibrated with mock spectra
from a large suite of hydrodynamical simulations post-processed with our
thermal IGM evolution code cite, finely sampling amplitude and slope of the
expected temperature-density relation. The thermal parameters inferred from
our measurements of the flux power spectrum, Doppler parameter distribution,
as well as wavelet and curvature statistics agree well within their respective
errors and all clearly show the peak in temperature and minimum in slope of
the temperature density relation expected from He ii reionization. Combining
our measurements from the different flux statistics gives $T_{0}=(14750\pm
1322)$K for the peak temperature at mean density and a corresponding minimum
slope $\gamma=1.225\pm 0.120$. The peak in the temperature evolution occurs at
$z\approx 3$, in agreement with previous measurements that had suggested the
presence of such a peak, albeit with a large scatter. Using cite, we also
calculate the thermal state of the IGM predicted by five widely used
(spatially homogeneous) UV-background models. The rather rapid thermal
evolution inferred by our measurements is well reproduced by two of the
models, if we assume (physically well motivated) non-equilibrium evolution
with photo-heating rates that are reduced by a moderate factor of $\sim
0.7-0.8$. The other three models predict He ii reionization to be more
extended with a somewhat earlier as well as higher temperature peak than our
measurements suggest.

###### keywords:

cosmology: large-scale structure of Universe - methods: numerical - galaxies:
intergalactic medium - QSOs: absorption lines

††pubyear: 2019††pagerange: A consistent and robust measurement of the thermal
state of the IGM at $2\leq z\leq 4$ from a large sample of Ly$\alpha$ forest
spectra: Evidence for late and rapid HeII reionization–I

##  1 Introduction

Ly$\alpha$ absorption seen in the spectra of distant bright QSOs (Quasi
Stellar Objects) allow one to probe the thermal and ionization history of the
Intergalactic Medium (IGM) in addition to constraining cosmological
parameters. The thermal state of the IGM is often characterized by
normalization $(T_{0})$ and slope $(\gamma)$ of the temperature-density
relation (TDR, $T=T_{0}\Delta^{\gamma-1}$; Hui & Gnedin, 1997), while the
ionization state of the IGM is characterized by H i and He ii photo-ionization
rates (e.g. Haardt & Madau, 1996). These parameters have been measured using
high-resolution, high signal-to-noise spectroscopic observations in
conjunction with high-resolution hydrodynamical simulations.

Accurate measurements of thermal parameters and photo-ionization rates of the
IGM have been used to place constraints on (i) the epoch and extent of H i and
He ii reionization (Worseck et al., 2011; Puchwein et al., 2015; Worseck et
al., 2016; Upton Sanderbeck et al., 2016; Gaikwad et al., 2019; Worseck et
al., 2019; Upton Sanderbeck & Bird, 2020), (ii) ionizing ultra-violet
background (UVB) models that are important inputs for cosmological
hydrodynamical simulations (Haardt & Madau, 2012; Khaire & Srianand, 2015;
Oñorbe et al., 2017; Khaire & Srianand, 2019; Puchwein et al., 2019; Faucher-
Giguère, 2020, hereafter ; ; ; ; ;  respectively), (iii) the escape fraction
of H i ionizing photons from galaxies and the relative contribution of
galaxies and QSOs to the total H i ionizing background (Khaire et al., 2016;
Khaire, 2017), (iv) and the effect of non-radiative processes like heating by
cosmic rays (Nath & Biermann, 1993; Samui et al., 2005), blazars (Chang et
al., 2012; Puchwein et al., 2012) and dark matter annihilation (Cirelli et
al., 2009; Liu et al., 2020).

The thermal parameters of the IGM have thereby been measured by (i)
decomposing absorption features into multi-component Voigt profiles and
identifying the lower envelope of the line width (b-parameter) vs H i column
density (Schaye et al., 1999, 2000; Bolton et al., 2014; Hiss et al., 2018;
Telikova et al., 2019), (ii) measuring the small scale suppression in the
transmitted flux power spectrum (McDonald, 2003; Kim et al., 2004; Walther et
al., 2019; Boera et al., 2019), (iii) using curvature (defined as,
$\kappa=F^{\prime\prime}/[1+(F^{\prime})^{2}]^{3/2}$ where $F^{\prime}$ and
$F^{\prime\prime}$ are the first and second derivatives of the normalised flux
with respect to the wavelength) statistics (Becker et al., 2011; Boera et al.,
2014; Padmanabhan et al., 2015) or (iv) characterising Fourier modes of the
flux distribution in the range sensitive to thermal parameters using wavelet
analysis (Zaldarriaga, 2002; Theuns et al., 2002; Lidz et al., 2010; Garzilli
et al., 2012). Even if the same observational data is used, the sensitivity to
small changes in thermal parameters and associated systematic uncertainties
are found to be different for different flux statistics (see §5 for a detailed
discussion of this). Ideally, one would thus like to calculate all these
different flux statistics simultaneously and consistently for the same data
set to obtain joint best fit values of the thermal parameters and the
associated errors. In order to obtain robust and consistent results it is also
important to perform such an analysis using a consistent set of model
parameters.

As simulated data are an integral part of the parameter estimation from
Ly$\alpha$ forest data, the reliability and accuracy of the extracted
parameters depend on the assumptions involved in the simulations and the
ability to perform a wide range of simulations that span a sufficiently wide
parameter space with good i.e. not too coarse sampling. Generating such a set
of simulations accounting for the relevant physical effects is challenging. At
any given epoch the pressure smoothing of the density field depends on the
past thermal history (Gnedin & Hui, 1998; Peeples et al., 2010; Kulkarni et
al., 2015; Nasir et al., 2016; Rorai et al., 2017a) and one needs self-
consistent high-resolution cosmological hydrodynamical simulations of the IGM
(Springel, 2005; Almgren et al., 2013; Bolton et al., 2017). Practically one
achieves a range of $T_{0}$ and $\gamma$ in simulations by varying photo-
heating rates of H i and He ii as a function of redshift (see for example,
Becker et al., 2011). Such simulations are computationally expensive which
limits the thermal parameter space that can be probed at any given redshift.

As noted before, most simulations used for this purpose are performed assuming
a uniform ionizing background and ionization equilibrium. While these are good
approximations after He ii reionization is complete such simulations do not
capture all relevant physics before and during He ii reionization, i.e.,
$2.7\leq z\leq 4.0$ (see Puchwein et al., 2015, 2019; Gaikwad et al., 2019).
For this one needs simulations that also account for the non-equilibrium
effects during He ii reionization. Ideally one would also like the simulations
to incorporate the spatially inhomogeneous nature of He ii reionization and
the corresponding fluctuations of the UV background. However, such simulations
will not only require to include radiative transfer, but also large box sizes
as well as high resolution and thus very high dynamic range to capture the
relevant physical processes accurately.

The main aim of this work is to obtain a consistent measurement of thermal
parameters using the larger data sets that have become publicly available
recently. Thanks to compilations like kodiaq dr2 (O’Meara et al., 2015, 2017,
from KECK/HIRES archival data) and UVES squad dr1 (Murphy et al., 2019, from
VLT/UVES archival data), large samples of QSO spectra with high resolution
($\sim 6$ km s-1 adequate to resolve the thermally broadened Ly$\alpha$
absorption lines) and high S/N are now available for analysis. These samples
have dramatically increased the number of available QSO spectra that have been
reduced and continuum normalised using uniform techniques.

Additional impetus to perform a consistent measurement of thermal parameters
for these large samples comes from two tools we have developed over the past
few years: (i) the Code for Ionization and Temperature Evolution (cite,
Gaikwad et al., 2017a) and (ii) the VoIgt profile Parameter Estimation Routine
(viper, Gaikwad et al., 2017b), an automatic Voigt profile fitting code that
decomposes Ly$\alpha$ absorption spectra into Voigt profile components. cite
not only allows us to construct models with a wide range of thermal and
ionization histories efficiently without running full hydrodynamical
simulations, but also enables us to calculate the non-equilibrium evolution of
the thermal and ionization state of the gas (see Gaikwad et al., 2019, for
details). cite has been shown to reproduce the results of full hydro-
simulations to well within 10 percent accuracy (Gaikwad et al., 2018). viper
runs on parallel architectures and allows us to perform a Voigt profile
analysis for large samples of observed and simulated spectra consistently (see
Maitra et al., 2019, 2020, for a clustering analysis using Voigt profile
components). Thanks to these two codes, we are in a position to measure the
physical parameters of the IGM from Ly$\alpha$ forest data by simultaneously
using the different statistics mentioned above.

We present here the measurements of thermal parameters ($T_{0}$ and $\gamma$)
over the redshift range $2\leq z\leq 4$ in 10 redshift bins of width $\Delta
z=0.2$ using 103 high resolution, high S/N QSO absorption spectra drawn from
the kodiaq dr2 sample. The paper is organized as follows, in §2 we present the
details of the observational data used in this work and compare our new
measurements of the mean flux as a function of redshift, H i column density
distribution, flux PDF and power-spectrum with measurements in the literature.
We describe our simulations and explain how we generate simulated spectra for
a wide range of finely spaced thermal parameters in §3. We provide details of
four different flux statistics of Ly$\alpha$ forest used to measure thermal
parameters in §4 and in Appendix E. In §5 we show our measurements of thermal
parameters, present our error analysis and compare our measurements with
measurements in the literature. In §6 we show predictions for the thermal
parameter evolution for different UVB models. Finally we summarize in §7. Note
that we make all the measurements from observational data available to the
community on request for future use.

The default cosmological parameters used here are
$(\Omega_{\Lambda},\Omega_{m},\Omega_{b},\sigma_{8},n_{s},h,Y)=(0.69,0.31,0.0486,0.83,0.96,0.674,0.24)$,
consistent with a flat $\Lambda$CDM cosmology (Planck Collaboration et al.,
2014). All distances are given in comoving units unless specified. We have
expressed $\Gamma_{\rm HI}$ in units of $10^{-12}\>{\rm s}^{-1}$ denoted by
$\Gamma_{\rm 12}$. We often refer to $T_{0},\gamma$ as thermal parameters in
this work.

##  2 Observations

Figure 1: Each horizontal line shows the Ly$\alpha$ redshift range corresponding to the wavelength range between Ly$\alpha$ and Ly$\beta$ emission for individual QSOs in our sample. We exclude QSO proximity regions of $10$ pMpc from the QSO towards us. The spectra are divided in ten redshift bins, each of size $\Delta=0.2$ centered on $z=2.0,2.2,\cdots,3.8$ as shown by the vertical dashed lines. The number of lines of sight in a given redshift bin are shown in the boxes.  Table 1: Observed mean Ly$\alpha$ transmitted flux $z\>\pm\>dz$ | ${\rm N_{los}}$ | Median SNR |  $\langle F\rangle$ $\pm$ $d\langle F\rangle$
---|---|---|---
$2.0\>\pm\>0.1$ | 32 | 14.7 | $0.8690\>\pm\>0.0214$
$2.2\>\pm\>0.1$ | 37 | 21.2 | $0.8261\>\pm\>0.0206$
$2.4\>\pm\>0.1$ | 33 | 19.2 | $0.7919\>\pm\>0.0210$
$2.6\>\pm\>0.1$ | 26 | 19.3 | $0.7665\>\pm\>0.0216$
$2.8\>\pm\>0.1$ | 28 | 19.7 | $0.7398\>\pm\>0.0212$
$3.0\>\pm\>0.1$ | 26 | 23.4 | $0.7105\>\pm\>0.0213$
$3.2\>\pm\>0.1$ | 16 | 22.0 | $0.6731\>\pm\>0.0223$
$3.4\>\pm\>0.1$ | 12 | 22.4 | $0.5927\>\pm\>0.0247$
$3.6\>\pm\>0.1$ | 5 | 22.8 | $0.5320\>\pm\>0.0280$
$3.8\>\pm\>0.1$ | 7 | 22.3 | $0.4695\>\pm\>0.0278$

We use observed spectra from the second data release of the Keck Observatory
Database of Ionized Absorption toward Quasars (kodiaq dr2) survey (O’Meara et
al., 2015, 2017)111https://www2.keck.hawaii.edu/koa/public/koa.php. The sample
consists of 300 QSO spectra with emission redshifts $z\leq 5.3$. All available
spectra are continuum normalised and the data product provides normalised flux
and the associated error as a function of wavelength. Many of these QSOs were
observed more than once with different exposure times. We co-added all the
spectra using the procedure described in Appendix A.

In total there are 214 QSO spectra that cover the Ly$\alpha$ forest in the
range $1.9\leq z\leq 4$. We manually checked all the coadded spectra and
excluded spectra if one or more of the following criteria are satisfied: (i)
sightline does not (partially or fully) contain Ly$\alpha$ forest in the
redshift range $1.9\leq z\leq 4$, (ii) sightline contains Damped Ly$\alpha$
(DLA) or sub-DLA systems, (iii) the sightline contains large spectral gaps, or
(iv) the median S/N per pixel along the sightlines is smaller than $5$. After
excluding the QSO spectra with above criteria, the resulting sample consists
of 103 QSO spectra. Fig. 1 shows the Ly$\alpha$ redshift range corresponding
to the wavelength range between the Ly$\alpha$ and Ly$\beta$ emission lines of
the 103 QSOs in our sample. The Ly$\alpha$ absorption close to the QSOs is
expected to be influenced by the enhanced ionizing flux due to the QSOs
(Carswell et al., 1982; Kulkarni & Fall, 1993; Lidz et al., 2007; Calverley et
al., 2011; Bolton et al., 2012). We exclude such biased QSO proximity regions
(i.e 10 pMpc from the QSO towards us) in our subsequent analysis.

In order to measure the evolution of thermal parameters, we divide our sample
into ten redshift bins centered on $z=2.0,2.2,2.4,\cdots,3.6~{}{\rm
and}~{}3.8$ with a bin width of $\Delta z=0.2$. Fig. 1 shows the number of QSO
spectra contributing to our sample in each of these ten redshift bins. The
properties of the observed spectra in these redshift bins are summarized in
Table 1. The median S/N of the observed sample is $>10$ in all redshift bins.
The observed Ly$\alpha$ forest regions are usually contaminated by metal lines
that produce narrow absorption features and potentially could bias our
measurements of $T_{0}$ and $\gamma$. To account for this we manually
identified contaminating metal line absorption using the list of known
intervening metal line systems along these sightlines. It is not always
possible to identify all the metal lines contaminating the Ly$\alpha$ forest.
However, as we fit all the observed spectra with our automated Voigt profile
fitting routine viper (Gaikwad et al., 2017b) we can mitigate this. We treat
all lines with $b\leq 8$ km s-1 as metal lines. Finally, we replace all metal
lines by continuum and add noise to the replaced regions. We have checked the
effect of residual metal line contamination on our $T_{0}-\gamma$ measurements
and found the effect to be marginal (see Appendix G).

There could be additional errors in the observed flux due to continuum fitting
uncertainties. Continuum fitting uncertainties depend not only on the number
of unabsorbed spectral regions used in the fit but also on the observed S/N
per pixel in these regions and the QSO spectral energy distribution (SED).
Since the true QSO continuum is unknown, an exact quantification of the
contribution of continuum fitting uncertainty to the error in the normalised
flux is not available for the kodiaq dr2 sample. The continuum fitting
uncertainty is expected to be of the order of a few percent for moderate S/N
data (O’Meara et al., 2015). We thus allow for the possibility of a systematic
error of $\pm 5$ per cent in the normalised flux for our final $T_{0}$ and
$\gamma$ measurements to account for continuum uncertainties.

###  2.1 Comparison of observed statistics with previous measurements

Figure 2:  The redshift evolution of the observed mean flux from this work
(red stars with errorbar) is compared with those from Kirkman et al. (2005);
Kim et al. (2007); Faucher-Giguère et al. (2008); Becker et al. (2013). The
$\langle F\rangle$ evolution in Becker et al. (2013) is obtained from a large
number of moderate resolution and low SNR SDSS spectra, while the $\langle
F\rangle$ evolution in Kirkman et al. (2005); Kim et al. (2007); Faucher-
Giguère et al. (2008) were obtained using high resolution QSO absorption
spectra. Our $\langle F\rangle$ measurements (based on kodiaq DR2 sample) are
in agreement with Kirkman et al. (2005, maximum deviation $\sim 2\sigma$);
Becker et al. (2013, maximum deviation $\sim 2\sigma$). However, there is a
slight increase in $\langle F\rangle$ at $2.9\leq z\leq 3.3$ compared to
Becker et al. (2013). The presence or absence of this possible excess in
$\langle F\rangle$ has only a marginal effect on our measured $T_{0}-\gamma$
but may be interesting in its own right as we will discuss later (see §2 for
details).  Figure 3:  Panels A1, B1 and C1 show a comparison of our
measurements (blue solid curves) of the flux PDF/power spectrum (FPDF/FPS) and
column density distribution (CDDF) with other measurements in the literature.
Panels A2, B2 and C2 show the differences between our measurements and those
from the literature. Kim et al. (2007); Calura et al. (2012); Rollinde et al.
(2013) calculated the FPDF at $2.71\leq z\leq 3.21$ using 8, 2 and 5 high-
resolution UVES spectra, respectively, while our sample consists of $\sim 25$
high-resolution spectra. Panels A1 and A2 show that the maximum difference
between our FPDF (in the range $0.1\leq F\leq 0.9$) and that of Kim et al.
(2007); Calura et al. (2012); Rollinde et al. (2013) is $\sim 13$ percent.
Panel B1 shows the comparison of our FPS with that from Walther et al. (2018)
at $2.9\leq z\leq 3.1$. Unlike Walther et al. (2018), we neither subtract
noise power from the FPS nor do we account for window effects in the FPS
calculation. Note that the effect of noise and finite resolution of the
spectrograph play a dominant role for the FPS measurement at $k>0.15\;{\rm
km^{-1}\>s}$ which we do not use for our measurements. We use the FPS at
$0.01\leq k\>{\rm(s\>km^{-1})}\leq 0.1$ to measure thermal parameters. In this
range our FPS is in good agreement (maximum difference is 19 percent) with
that from Walther et al. (2018). Panels C1 and C2 show that our CDDF is within
3.75 percent of that from Kim et al. (2013) at $1.9\leq z\leq 3.2$. In panel
C1, we define the CDDF similar to Kim et al. (2013). However, for rest of this
paper we define the CDDF in terms of $\log{\rm N_{HI}}$ and $dz$. A comparison
of errors of the observed flux statistics from this work with the literature
is discussed in detail in appendix C.

In this section, we derive the statistics of the transmitted Ly$\alpha$ flux
for our sample and compare them with results in the literature. In particular,
we compare the evolution of mean flux, flux probability distribution function
(FPDF), flux power-spectrum (FPS) and H i column density distribution function
(CDDF) from our sample with other measurements in the literature. In §4 we
describe in detail our method of calculating these flux statistics (along with
other statistics) from observations and simulations.

The last column in Table 1 provides the mean observed transmitted flux
measured in each redshift bin. The error on the mean flux given also accounts
for the systematic uncertainty due to continuum placement. As expected, the
mean flux decreases monotonically with redshift due to the increase in the
opacity of the IGM at higher redshifts. In Fig. 2, we compare the evolution of
the mean flux $\langle F\rangle$ of our sample (see Table 1) with that from
Kirkman et al. (2005); Kim et al. (2007); Faucher-Giguère et al. (2008) and
Becker et al. (2013). For clarity, we do not show the errors of $\langle
F\rangle$ for the measurements from the literature. The $\langle F\rangle$
evolution in Kirkman et al. (2005); Kim et al. (2007) and Faucher-Giguère et
al. (2008) is obtained using high-resolution QSO absorption spectra, while
Becker et al. (2013) obtained their $\langle F\rangle$ evolution using very
large number of SDSS spectra that have lower resolution and low SNR compared
to the typical high-resolution data used in the literature and our work here.

Our $\langle F\rangle$ evolution is broadly consistent with Kirkman et al.
(2005); Kim et al. (2007); Faucher-Giguère et al. (2008); Becker et al. (2013)
with some notable differences. $\langle F\rangle$ at $z=2.2,3.0$ and $3.2$ in
our sample is smaller by $\sim 2\sigma$ than that of Kirkman et al. (2005);
Kim et al. (2007). One possible reason could be the number of QSO lines of
sight used in their work, $\leq 10$ at those redshifts. Our $\langle F\rangle$
measurement at $2.9\leq z\leq 3.5$ ( $2.1\leq z\leq 2.5$) is systematically
larger (smaller) than that of Becker et al. (2013). The $\langle F\rangle$
evolution in our sample is in good agreement with that from Faucher-Giguère et
al. (2008) over the full redshift range. We see a slight enhancement in
$\langle F\rangle$ at $2.9\leq z\leq 3.5$ in the kodiaq dr2 sample (see
appendix B for more details). We have also analyzed the squad dr1 QSO sample
for the purpose of calculating the $\langle F\rangle$ evolution (Murphy et
al., 2019). We find that the enhancement in $\langle F\rangle$ is less
prominent in squad dr1, but the $\langle F\rangle$ evolution shows a change in
slope at $z>3.2$. It is noteworthy that the statistics sensitive to thermal
parameters are in good agreement with each other for the kodiaq dr2 and squad
dr1 samples. Furthermore, note again that we rescale simulated optical depths
to match the observed $\langle F\rangle$ which reduces the effect of
differences in $\langle F\rangle$ on thermal parameters. The presence or
absence of the slight excess of $\langle F\rangle$ discussed above has thus
only a marginal effect on measurements of $T_{0}-\gamma$ presented in this
work.

We compare the FPDF at $2.71\leq z\leq 3.21$ from our sample with those from
Kim et al. (2007); Calura et al. (2012); Rollinde et al. (2013) in Fig. 3
(panel A1 and A2). The number of spectra used by Calura et al. (2012);
Rollinde et al. (2013) and Kim et al. (2007) is 2, 5 and 8 respectively, while
our sample contains $\sim 25$ spectra. Despite this, our FPDF in the range
$0.1\leq F\leq 0.9$ is within $\sim 13$ percent of that from Kim et al.
(2007); Calura et al. (2012); Rollinde et al. (2013). The FPDF statistics
appears reasonably well converged even for small number of spectra. In this
work, we focus on the FPDF in the flux range $0.1\leq F\leq 0.9$ because the
flux at $F<0.1$ ($F>0.9$) could be dominated by sky subtraction (continuum
placement) uncertainty.

In panel B1 and B2 of Fig. 3, we show a comparison of the FPS from our sample
with that from Walther et al. (2018). Even though the sample used in Walther
et al. (2018) and our work here is the same (i.e, the kodiaq DR2 sample), the
number of QSOs per redshift bin is different because of our selection criteria
(see §2). Our method of calculating the FPS is also different. Since our
sample is selected such that the spectra do not contain spectral gaps, we
compute the power spectrum using FFT. We also forward model the simulated
Ly$\alpha$ forest spectra to mimic the noise and instrumental broadening
properties of the observed spectra. Hence when calculating the FPS, we neither
subtract the noise power nor deconvolve the instrumental broadening. We also
replace metal lines with continuum and add noise in the replaced regions,
while Walther et al. (2018) mask metal line regions. Despite these
differences, when we compare our FPS at $0.01\leq k\>{\rm(s\>km^{-1})}\leq
0.1$ with that from Walther et al. (2018), we find good agreement. Note that
we use the FPS only at $0.01\leq k\>{\rm(s\>km^{-1})}\leq 0.1$ to measure
$T_{0}$ and $\gamma$, as these scales are most sensitive to the variation of
thermal parameters.

In Fig. 3 (panel C1-C2), we compare the CDDF from Kim et al. (2013) with that
obtained from our sample using viper. For a fair comparison, we restrict the
observed spectra of our sample to the redshift range $1.9\leq z\leq 3.2$ and
we use the same definition of the CDDF. Kim et al. (2013) present the CDDF for
the range $13\leq$ $\log{\rm N_{HI}}$ $\leq 18$ since their sample is
incomplete at $\log{\rm N_{HI}}$ $<13$. We account for the incompleteness in
our measurement by calculating the sensitivity curve hence our CDDF
measurements are shown from $\log{\rm N_{HI}}$ $=12.6$ to $15.8$. Our observed
CDDF is in good agreement ($\sim 3.75$ percent) with that of Kim et al.
(2013). At $\log{\rm N_{HI}}$ $>15$, we find more systems compared to Kim et
al. (2013). This may be due to the fact that viper fits only Ly$\alpha$
absorption while Kim et al. (2013) fit Ly$\alpha$ and Ly$\beta$ absorption
simultaneously. Despite this, the two CDDF in the range $13\leq$ $\log{\rm
N_{HI}}$ $\leq 15$ are consistent with each other. We refer the reader to
appendix C for a discussion of the associated errors in these statistics for
our measurements and those in the literature.

##  3 Simulations

We evolve the cosmological density, velocity and temperature field using the
smooth particle hydrodynamic (SPH) gadget-3 222http://wwwmpa.mpa-
garching.mpg.de/gadget/ code. The initial conditions are generated at $z=99$
using 2lpt (Scoccimarro et al., 2012)333http://cosmo.nyu.edu/roman/2LPT/. Our
fiducial simulations have a box size of $10\>h^{-1}$ cMpc and have $512^{3}$
baryon particles and an equal number of dark matter particles corresponding to
a gas particle mass of $\sim 10^{5}\>{\rm M_{\odot}}$. The gravitational
softening length is set to $\sim 0.65\>h^{-1}$ ckpc. We found this to be the
best compromise in terms of resolution, dynamic range and ability to run a
sufficiently fine grid of thermal parameters. We have performed a resolution
study using the Sherwood simulation suite
444https://www.nottingham.ac.uk/astronomy/sherwood/ and show that our results
are sufficiently well converged for our choice of simulation (see §4.1, Bolton
et al., 2017). We also use the Sherwood simulation suite to demonstrate that
our measurements are not significantly affected by the somewhat limited box-
size and thus spatial dynamical range.

gadget-3 is a modified version of the publicly available gadget-2 code
(Springel, 2005) that incorporates the radiative heating and cooling by a time
varying, spatially uniform UV background. For our fiducial simulations we
incorporate the KS19 UVB model (QSO SED index $\alpha=-1.4$) by modifying the
TREECOOL file in gadget-3. We store the output of the gadget-3 code at equal
redshift intervals $z=6.0,5.9,\cdots,2.1,2.0$. We employ a simplified star
formation criteria that converts particles with $\Delta>1000$ and $T<10^{5}$
into stars (the so called quick Ly$\alpha$ option, Viel et al., 2004a) and do
not include AGN or stellar feedback.

To generate physically motivated thermal histories one needs to perform
computationally expensive cosmological simulations for a range of UVB models
(Becker et al., 2011; Walther et al., 2019). In this work, we follow the
approach laid out by Gaikwad et al. (2018) and explore the thermal parameter
space efficiently. Our procedure to simulate the thermal history is: (i) We
perform a gadget-3 simulation with the KS19 UVB model using equilibrium
ionization evolution equation. The typical thermal parameters of such a
simulation at $2.0<z<4.0$ are $T_{0}\sim 9500$ K and $\gamma\sim 1.5$. (ii) To
obtain the variation in thermal parameters, we modify the photo-heating rates
in the UVB code (see §3). We then solve the ionization and thermal evolution
equation for each particle on gadget-3 outputs using our thermal evolution
code Code for Ionization and Temperature Evolution (cite, Gaikwad et al.,
2017a) (iii) We apply pressure smoothing corrections while extracting the
$n_{\rm HI},T$ and $v$ fields along a sightline by convolving the SPH kernel
with the pressure smoothing Gaussian kernel (Gaikwad et al., 2018). (iv)
Finally, we compute the Ly$\alpha$ optical depth accounting for thermal
broadening, natural line width broadening and peculiar velocity effects. In
Gaikwad et al. (2017a, 2019), we showed that thermal parameters can be
constrained within $1\sigma$ uncertainty provided the mean flux is matched
between simulation and observations. However unlike Gaikwad et al. (2019), in
this work, we solve for the ionization evolution assuming equilibrium while
solving the thermal evolution using the non-equilibrium equations.

We vary the thermal parameter evolution by scaling the H i, He i and He ii
photo-heating rates of the KS19 UVB (see Becker et al., 2011, for a similar
approach). We generate $T_{0}-\gamma$ combinations for a finely sampled grid
such that $T_{0}$ is varied from $\sim 6000$ K to $\sim 24000$ K in steps of
500 K while $\gamma$ is varied from $\sim 0.7$ to $\sim 2.0$ in steps of 0.05
at $z=3$. The corresponding $T_{0}$ and $\gamma$ values at other redshifts
differ by $\leq 7$ percent. We have thus simulated the Ly$\alpha$ forest for
$37\times 27=999$ different thermal histories. The computational time required
to run and extract Ly$\alpha$ forest spectra for 999 UVB models is around
$\sim 2.25$ million cpu hours. Note that the number of thermal parameters
probed in this work is larger by a factor of $\sim 13$ and $\sim 50$ than
those by Walther et al. (2019, $N_{\rm thermal}=76$) and Becker et al. (2011,
$N_{\rm thermal}=18$), respectively. We refer the reader to appendix D for
more details.

We generate the overdensity ($\Delta$), neutral fraction ( $f_{\rm HI}$),
temperature ($T$) and peculiar velocity $(v)$ fields along 20000 random
skewers through our simulation box at $z=1.9,2.0,2.1,\cdots,3.7,3.8~{}{\rm
and}~{}3.9$. In each redshift bin we concatenate the fields to match the
observed redshift path length ($\Delta z=0.2$). However, while computing the
FPS, we do not concatenate spectra as there is no correlation beyond the
length of the simulation box. The Ly$\alpha$ optical depth is generated from
the $\Delta$, $f_{\rm HI}$, $T$ and $v$ fields accounting for Doppler
broadening, natural line broadening and peculiar velocity effects (Choudhury
et al., 2001; Padmanabhan et al., 2014). For each redshift bin, we construct a
single simulated mock sample with the same number of spectra as observed in
that redshift bin (see Gaikwad et al., 2017a, for a similar method at
$z<0.5$). The simulated spectra in a mock sample mimic the observed spectra by
(i) resampling the pixel distribution similar to the observed data, (ii)
convolving the flux field with Gaussian instrumental broadening (FWHM $\sim 6$
km s-1) and (iii) by adding Gaussian random noise generated from the observed
S/N per pixel array. We generate 100 such mock samples to constitute a mock
suite for each redshift bin. For example, in the redshift range $2.9\leq z\leq
3.1$, each mock sample consists of 26 simulated spectra and a mock suite
consists of $26\times 100=2600$ simulated spectra. We use these 100 mock
samples to estimate the errors on different flux statistics. In order to
compare the simulations with observations, we calculate the flux statistics
that we will show to be sensitive to thermal parameters in the next section.

##  4 Method

Figure 4: Effects of changing the temperature field on the properties of the
Ly$\alpha$ forest. Panel A shows the simulated Ly$\alpha$ flux from a hot
($T_{0}\sim 20000$ K ) and a cold ($T_{0}\sim 10000$ K) model accounting for
instrumental broadening and finite S/N. As expected, the absorption features
are systematically broader and shallower in the hot model than in the cold
model. Panel B shows the wavelet amplitude from the two models obtained by
convolving the Ly$\alpha$ flux field from panel A with a wavelet of scale 50
km s-1 (see appendix E.2). Panel C shows the corresponding curvature
amplitudes (see appendix E.3). The wavelet and curvature amplitudes are
systematically smaller for the hot model as compared to those from the cold
model. Panels D and E show the Voigt profile decomposition of spectra for the
cold and hot model obtained using viper respectively. The Doppler parameter
$b$ is systematically larger for the hot model compared to the cold model (see
appendix E.4). The black curve in panels D and E show the best-fit Voigt
profiles. The green curve near $F=0$ in these panels shows the residual
between input and fitted flux. Vertical ticks in panel D and E show the
locations of statistically significant Voigt components identified and fitted
by viper automatically.

The increase in temperature of the IGM due to He ii reionization has the
following main effects on the H i Ly$\alpha$ forest (i) broadening of the
absorption features (Schaye et al., 2000; Lidz et al., 2010; Becker et al.,
2011), (ii) decrease in neutral fraction due to the temperature dependence of
the recombination rate (Rauch et al., 1997; Bolton & Haehnelt, 2007; Becker &
Bolton, 2013; Viel et al., 2017; Khaire et al., 2019) and (iii) pressure
smoothing of the density and flux fields (Gnedin & Hui, 1998; Peeples et al.,
2010; Kulkarni et al., 2015; Lukić et al., 2015; Rorai et al., 2017a; Maitra
et al., 2019). In the literature, $T_{0}$ and $\gamma$ have been measured by
two different kind of statistics, namely, those derived from the transmitted
Ly$\alpha$ directly and those derived by fitting the transmitted Ly$\alpha$
flux with multi-component Voigt profiles. We use here four statistics, namely,
the flux power spectrum, the wavelet statistics, the curvature statistics and
the line width $(b)$ distribution (hereafter BPDF) of Voigt profile components
to measure the thermal parameters. In Fig. 4, we illustrate the effect of
increased temperature on the Ly$\alpha$ flux, wavelet field and curvature
field. We also show examples of viper fits (along with residuals) to spectra
of the two models. Our method of computing these statistics is described in
detail in appendix E.

###  4.1 Validation of our approach using mock data

Figure 5:  The left, middle and right panels show the recovery of thermal
parameters using FPS, BPDF, wavelet and curvature statistics at
$z=2,3.2~{}{\rm and}~{}S3.8$ respectively. The $1\sigma$ joint constraints
(using all statistics simultaneously) on the thermal parameters are shown by
black contours and the shaded regions. We here treat the Ly$\alpha$ forest
generated from the large dynamic range L40N2048 Sherwood simulation as mock
data. We post-process our fiducial L10N512 gadget-3 simulation with cite and
generate model Ly$\alpha$ forest spectra for 999 different thermal parameters
by scaling the KS19 photo-heating rates. The initial conditions of the L10N512
simulation used in this work are different from any of the ICs used in the
Sherwood simulation suite. The mock and model Ly$\alpha$ forest data both
mimic the observational properties at the respective redshifts. The true
values (cyan stars) of the Sherwood simulation used to create the mock data
sample lies within the gray shaded region at all three redshifts. Thus our
method successfully recovers the thermal parameters despite the smaller
spatial dynamic range of the simulations used and the approximate nature of
our modeling of the thermal history.

Before measuring the thermal parameters from observations, we demonstrate the
accuracy of our method in recovering the thermal parameters with an end-to-end
test on mock data created from our fiducial hydro-simulation(s) and from the
much larger dynamic range simulation from the Sherwood simulation suite
(computed with cosmological parameters similar to what we use in our models).
This is with the aim to (i) test the sensitivity of statistics to thermal
parameters, (ii) study the degeneracy between thermal parameters, (iii)
quantify the accuracy of the method, (iv) check if there are any systematic
effects between true and recovered thermal parameters, (v) check if the
simulations are sufficiently converged and (vi) test for the effect of Jeans
smoothing.

The effect of box size and resolution on statistics of the Ly$\alpha$ forest
using a range of Sherwood simulations are discussed in appendix F. As shown in
Fig. 25, the L10N512 simulation used in this work is sufficiently converged
for the corresponding Sherwood model. To test how well our models can recover
the thermal parameters, we generate Ly$\alpha$ forest statistics in three
redshift bins from mock data based on the L40N2048 simulations in the Sherwood
suite which has the same resolution and four times larger box size than our
fiducial hydro-simulation The flux statistics are generated from $\sim 20000$
skewers as explained in §3. The mock data and our model spectra both closely
mimic the properties of the observed sample e.g. number of QSOs, S/N,
resolution and redshift path length in the corresponding redshift bin. The
errors for each statistics are calculated as discussed in §E.5. When using
FPS, BPDF, wavelet and curvature PDF to characterise the thermal parameters we
compute the $\chi^{2}$ between model and mock data (see appendix E.5 for
details).

Fig. 5 shows the recovery of thermal parameters using our models and the mock
data generated from the L40N2048 Sherwood simulation. For all four statistics,
the true thermal parameters of the mock data are within the $1\sigma$ contours
and are thus close to the best fit measured thermal parameters. The true
thermal parameters are also well within the $1\sigma$ contours of the measured
values if the four statistics are combined (shown by the black contours and
the shaded region). As expected the joint analysis contours are narrower than
those from individual statistics. The best fit thermal parameters do not show
any systematic deviation from the true thermal parameters indicating that any
uncertainty in recovery is likely to be statistical in nature.

###  4.2 The effect of varying Jeans smoothing

The width of absorption features depends on the instantaneous thermal state as
well as on the entire past thermal history due to the smoothing effect of the
thermal pressure on the spatial distribution of the gas, an effect that has
been dubbed Jeans smoothing (Gnedin & Hui, 1998; Peeples et al., 2010;
Kulkarni et al., 2015). At the redshifts considered here, well past H i
reionization and with the additional heat input from He ii reionization
ongoing, the memory of different possible hydrogen reionization histories is
expected to be modest. We have verified this by repeating our end-to-end test
with mock data from an otherwise identical simulation of the Sherwood suite
which has the same instantaneous temperatures, but where hydrogen reionization
occurs considerably later (the zr9 Sherwood simulation). The difference in the
measured thermal parameters was hardly noticeable ($<0.8$ percent). The range
of possible He ii reionization histories consistent with the He ii opacity
data is also rather small (Worseck et al., 2019). Hence the effect of the
differences in Jeans smoothing for the same instantaneous temperature and
different He ii reionization histories is again expected to have a small
effect on our measurements of thermal parameters.

###  4.3 The effect of spatial fluctuation of the temperature-density relation
due to inhomogeneous He ii reionization

In reality, He ii reionization will be spatially inhomogeneous and instead of
a well defined temperature density relation there will be a range of
temperature density relations depending on when the He ii in a particular
region was reionized (Rorai et al., 2017b; Upton Sanderbeck & Bird, 2020). In
Gaikwad et al. (2020) we have shown that for the equivalent situation during H
i reionization assuming a homogeneous UVB model nevertheless recovers the
median thermal state reasonably well. We have verified that this should also
be the case for He ii reionization by producing flux statistics for mock
spectra obtained with a wide range of temperature density relations. As the
box size of our simulations is much smaller than the expected size of a region
where He ii is reionized simultaneously this should mimic the effect of
inhomogeneous He ii reionization reasonably well. The thermal parameters
measured from these samples of simulated mock spectra with a range of
temperature-density relations are indeed within 1$\sigma$ of the median
values. We are not sure if this can be interpreted as a systematic bias, but
we note that at and after the peak in temperature the inferred $T_{0}$ values
were about 1$\sigma$ lower than the median while at $z>3$ the difference was
hardly noticeable.

##  5 Results

In this section, we present our measurements of $T_{0}$ and $\gamma$ from the
kodiaq DR2 sample. In order to make a fair comparison, we (i) derive the
statistics from simulations and observations in the same way, (ii) calculate
the error on each statistics from simulations and/or observations, (iii)
minimize the $\chi^{2}$ between data and model and (iv) find the best fit
model that corresponds to the minimum $\chi^{2}$ and the $1\sigma$ uncertainty
on the parameters corresponding to $\chi^{2}_{\rm min}\pm\Delta\chi^{2}$
($\Delta\chi^{2}=3.50$ for 2 parameters + $\langle F\rangle$ normalization,
Avni, 1976). We refer the reader to appendix E.5 for more details on our error
estimation and $\chi^{2}$ minimization method.

###  5.1 $T_{0}-\gamma$ constraints in the redshift range $1.9\leq z\leq 3.9$

Figure 6: Panel A1 shows $1\sigma$ constraints for $T_{0}$ and $\gamma$ from
the wavelet PDF (blue dashed contour), curvature PDF (magenta dotted contour),
BPDF (red solid contour) and FPS (green solid contour) statistics individually
at $2.9\leq z\leq 3.1$. The joint constraints from all the four statistics is
shown by the gray shaded region. The constraints are obtained by minimizing
the $\chi^{2}$ between data and model. The wavelet constraints correspond to
joint constraints from wavelet PDFs for wavelet scales $s_{\rm
n}=30,40,\cdots,100$ km s-1. The BPDF constraints correspond to joint
constraints from the BPDF calculated in 9 different $\log{\rm N_{HI}}$ bins
[13.0,13.2], [13.2,13.4], $\cdots$, [14.4,14.6]. The curvature PDF gives
tighter constraints for $T_{0}$ as compared to $\gamma$. The $1\sigma$
contours from all the statistics are in good agreement with each other. Panels
A2 and A3 show the marginalized $\gamma$ and $T_{0}$ distribution from all
four statistics and our joint analysis. The best fit $T_{0}$ and $\gamma$
values (cyan star in panel A) lie in the overlapping region of constraints
from the the four flux statistics. We show the 1D reduced $\chi^{2}$
distribution in Fig. 26.  Figure 7: Panel A1 shows a comparison of the FPS
from the observations (blue squares with dashed line) with that of the best
fit model [$T_{0}=14500K$ and $\gamma=1.25$] (red solid line) at $2.9\leq
z\leq 3.1$ generated by cite. The measurement errors are calculated from the
observed data using the bootstrap method (gray shaded region). Panel A2 shows
the residuals between the observed distribution and the best fit model. Panel
B1, C1, D1, E1 and F1 are similar to panel A1 except the observations and
predictions of the best fit models are shown for wavelet PDF $(s_{\rm n}=50$
km s-1) , BPDF ($13.4\leq$ $\log{\rm N_{HI}}$ $\leq 13.6$), curvature PDF,
wavelet PDF $(s_{\rm n}=70$ km s-1) and BPDF ($13.8\leq$ $\log{\rm N_{HI}}$
$\leq 14.0$), respectively. Panel B2 to F2 are similar to panel A2 and show
again the residuals between observations and best fit model. For all the
statistics, the best fit model is in good agreement ($\sim 1\sigma$) with
observations. Note that the best fit $T_{0},\gamma$ given in Table 2 and Table
5 corresponds to interpolated best fit values.

Fig. 6 shows measurements of $T_{0}$ and $\gamma$ in redshift bin $2.9\leq
z\leq 3.1$. Panel A1 shows the measurements of thermal parameters from
individual as well as the combined statistics. The constraints from all four
statistics agree with each other within the $1\sigma$ uncertainties. The
constraints on $\gamma$ from the curvature statistics are poor compared to
those from the other statistics. This is because the curvature PDF is less
sensitive to $\gamma$ compared to $T_{0}$. In previous measurements using
curvature statistics, $T_{0}$ is measured at a characteristic density
$\overline{\Delta}$ for an assumed value of $\gamma$. This does not allow for
an independent measurement of $\gamma$. It is also evident that $T_{0}$ and
$\gamma$ are anti-correlated for the wavelet statistics, BPDF and FPS, while
there is little or no correlation between $T_{0}$ and $\gamma$ for the
curvature statistics. Panels A2 and A3 show measurements of $\gamma$ and
$T_{0}$ marginalizing over $T_{0}$ and $\gamma$, respectively. The
marginalized distributions also show that the $T_{0}$ and $\gamma$
measurements using different statistics are consistent with each other within
$1\sigma$. Note that the data and model Ly$\alpha$ forest spectra are the same
in all cases only the method/statistics used to measure $T_{0}-\gamma$ is
different. In order to get tighter and more robust measurements of $T_{0}$ and
$\gamma$, we also plot the joint constraints for $T_{0}$ and $\gamma$ from all
the statistics in panels A1-A3 of Fig. 6. We calculate the joint constraints
by adding the $\chi^{2}$ between model and data for all the four statistics.
For simplicity, we ignore the correlation among different statistics. Hence
the statistical uncertainty in joint constraints may be somewhat
underestimated. However, as shown in Fig. 5 when combining the statistics, the
fiducial $T_{0}-\gamma$ is recovered well within $1\sigma$. This suggests that
the statistical uncertainty for the measurement from the combined statistics
is realistic. The joint constrains for $T_{0}$ and $\gamma$ are indeed tighter
than the corresponding constraints from any individual statistics. The best
fit value ($T_{0}=14750$ K and $\gamma=1.23$) for the joint constraints from
all the four statistics is shown as the cyan star in panel A1 of Fig. 6 (also
see Fig. 26).

Fig. 7 shows the comparison of FPDF, BPDF, wavelet PDF and curvature PDF from
observations with the best fit model at $2.9\leq z\leq 3.1$. The corresponding
residuals are shown in the lower panels. The best fit models for all the
statistics are in $1\sigma$ agreement with the corresponding observations. The
$\chi^{2}$ per degree of freedom for all the statistics varies in the range
$0.7\leq\chi^{2}_{\rm dof}\leq 1.35$ (see Fig. 26). The predictions of our
best fit model obtained from the joint analysis are thus in reasonably good
agreement with observations for all the statistics.

Figure 8: Each panel is same as panel A1 in Fig. 6 except that the $1\sigma$
contours for $T_{0}$ and $\gamma$ from all the statistics are shown for
different redshift bins. The $1\sigma$ contours for different statistics are
in good agreement with each other at all redshifts. The joint constraints (the
gray shaded regions) are more robust and tighter than the constraints from the
individual statistics. The redshift evolution of $T_{0}$ and $\gamma$ is
clearly evident from this figure. $T_{0}$ and $\gamma$ are anti-correlated as
expected. The degree of correlation decreases with increasing redshift as most
of the Ly$\alpha$ forest probes gas typically at cosmic mean density
($\Delta\sim 1$) at $z\sim 3.8$.  Figure 9: The top and bottom panels show the
evolution of $T_{0}$ and $\gamma$, respectively, obtained using five methods
(i) Wavelet PDF (black squares), (ii) Curvature PDF (green diamonds), (iii)
BPDF (red stars), (iv) FPS (blue circles) and (v) joint constraints using all
the four statistics (magenta triangles). The errorbars show the $1\sigma$
statistical uncertainties on $T_{0}$ and $\gamma$ for all the 5 methods. All
$T_{0}$ and $\gamma$ measurements are consistent with each other in all
redshift bins. Note, however, that the errors for $T_{0}$ and $\gamma$ are
different for different methods due to different sensitivities of these
methods to variations of these parameters. The joint constraints minimize the
scatter of the best fit values and the resultant $1\sigma$ uncertainty on
$T_{0}$ and $\gamma$ is smaller than that for the corresponding individual
statistics.

Fig. 8 shows measurements of $T_{0}$ and $\gamma$ in the other redshift bins
spanning $1.9\leq z\leq 3.9$. In all the redshift bins, the $T_{0}$ and
$\gamma$ measurements from the individual statistics are in good agreement
with each other. The $\gamma$ measurements from the curvature statistics are
less tight than those from the other statistics. The uncertainty in
$T_{0}-\gamma$ using the joint constraints is significantly smaller than that
from the individual statistics at all redshifts. The joint $T_{0}-\gamma$
constraints are thereby in good agreement with the overlapping regions of all
the individual statistics. The variation of the $T_{0}$ and $\gamma$
measurements with redshift is also evident from Fig. 8. As expected the
$T_{0}$ and $\gamma$ constraints are anti-correlated in the lower redshift
bins where the absorption features probe densities above the mean. With
increase in redshift the $T_{0}$ and $\gamma$ correlation becomes weaker,
similar to the finding of Walther et al. (2019). This is because at $z\geq
3.5$ the Ly$\alpha$ forest is mostly sensitive to approximately mean cosmic
density $\Delta\sim 1$.

In Fig. 9, we compare the evolution of $T_{0}$ and $\gamma$ from the
individual statistics and joint analysis. The $T_{0}$ ($\gamma$) measurements
are marginalized over $\gamma$ ($T_{0}$). We can clearly see an evolution of
$T_{0}$ and $\gamma$ with redshift with a clearly identifiable peak in $T_{0}$
that coincides with a minimum in $\gamma$ at $z\sim 3$. The $T_{0}$ and
$\gamma$ measurements from each individual statistics as well as our joint
constraints are consistent within $1\sigma$ in all the redshift bins. As
expected, the $T_{0}-\gamma$ uncertainty for the joint analysis is smaller
than the corresponding uncertainty from the individual flux statistics. The
best fit $T_{0}$ and $\gamma$ values show some scatter between neighboring
redshift bins for the individual statistics due to the differences in
sensitivity of individual statistics to $T_{0}$ and $\gamma$ . The scatter of
the best fit values is considerably smaller for the joint measurements.

###  5.2 $T_{0}-\gamma$ Error Budget

Figure 10: The effect of continuum placement uncertainty and metal contamination on $T_{0}$ and $\gamma$ constraints from the joint analysis of data for $2.9\leq z\leq 3.1$. Comparison of the green contours (with metals) with the red contour (without metals) shows that $T_{0}$ and $\gamma$ are underpredicted by $\sim 2$ percent if metal contamination is not accounted for. The effect of continuum placement uncertainty (assumed to be $\pm 5$ percent) biases the measurements such that a low continuum predicts lower $T_{0}-\gamma$ measurements by 3% or less. We have accounted for both continuum placement and metal contamination uncertainties in our final measurements (see Table 2, Table 5 and appendix G for details).  Table 2: Measurements of $T_{0}$, $\gamma$ with total uncertainty (see Table 5 for error budget). $z\>\pm\>dz$ |  $T_{0}\pm\delta T_{0}$ (K) | $\gamma\pm\delta\gamma$
---|---|---
2.0 $\pm$ 0.1 | 9500 $\pm$ 1393 | 1.500 $\pm$ 0.096
2.2 $\pm$ 0.1 | 11000 $\pm$ 1028 | 1.425 $\pm$ 0.133
2.4 $\pm$ 0.1 | 12750 $\pm$ 1132 | 1.325 $\pm$ 0.122
2.6 $\pm$ 0.1 | 13500 $\pm$ 1390 | 1.275 $\pm$ 0.122
2.8 $\pm$ 0.1 | 14750 $\pm$ 1341 | 1.250 $\pm$ 0.109
3.0 $\pm$ 0.1 | 14750 $\pm$ 1322 | 1.225 $\pm$ 0.120
3.2 $\pm$ 0.1 | 12750 $\pm$ 1493 | 1.275 $\pm$ 0.129
3.4 $\pm$ 0.1 | 11250 $\pm$ 1125 | 1.350 $\pm$ 0.108
3.6 $\pm$ 0.1 | 10250 $\pm$ 1070 | 1.400 $\pm$ 0.101
3.8 $\pm$ 0.1 | 9250 $\pm$ 876 | 1.525 $\pm$ 0.140

We consider the following important uncertainties when measuring $T_{0}$ and
$\gamma$: (i) modeling uncertainty, (ii) continuum placement uncertainty,
(iii) uncertainty due to metal contamination and (iv) cosmological parameter
uncertainty. Fig. 10 shows the effect of continuum placement uncertainty and
metal contamination uncertainty on our measurements of $T_{0}$ and $\gamma$
from our joint analysis of the four different flux statistics (for $2.9\leq
z\leq 3.1)$. The uncertainty in $T_{0}$ and $\gamma$ due to continuum
placement uncertainty (assumed to be $\pm 5$ percent, see O’Meara et al.,
2017) is less than 3 percent. The continuum placement uncertainty mainly
affects the mean flux of the observed sample. However, when measuring $T_{0}$
and $\gamma$, we rescale the optical depth in the simulations to match the
observed mean flux. Contamination by narrow metal absorption lines in the
observed Ly$\alpha$ forest can also potentially bias our $T_{0}$ and $\gamma$
measurements. We thus remove the metal lines from the Ly$\alpha$ forest to
minimize their effect. However, it is still possible that some metal lines are
not identified due to insufficient wavelength coverage in the higher
wavelength side of the Ly$\alpha$ emission and/or blending effects. By
applying a cutoff in the fitted $b$ parameters, we show in Fig. 10 that the
contribution of such metal lines to the $T_{0}-\gamma$ uncertainty is not
larger than $\sim 2$ percent. Similarly, we find that the uncertainty in
cosmological parameters can lead to 0.5 percent uncertainty in $T_{0}$ and
$\gamma$. One of the main uncertainties in the $T_{0}$ and $\gamma$
measurements comes from our modeling of the Ly$\alpha$ forest. Since we vary
the thermal state of the gas by post-processing of gadget-3 simulations, the
dynamical impact of pressure smoothing may not have been captured accurately.
We find that the effect of pressure smoothing on the $T_{0}$ and $\gamma$
uncertainty is below $2.5$ percent. We refer the reader to Appendix G for a
detailed discussion of the effect of all the above uncertainties on our
$T_{0}-\gamma$ measurements.

Table 2 provides our $T_{0}$ and $\gamma$ measurements at different redshifts.
In Table 5 we give $T_{0}$ and $\gamma$ errors contributed by various above
listed uncertainties. The continuum placement uncertainty is systematic, while
the other uncertainties are statistical in nature. We add the statistical
uncertainties in quadrature while the systematic uncertainty is additive in
nature. In the rest of the paper, we show and discuss constraints on $T_{0}$
and $\gamma$ with total uncertainties as given in Table 2.

###  5.3 Evolution of $T_{0}$ and $\gamma$

Figure 11: The figure shows a comparison of our $T_{0}-\gamma$ measurements
with measurements from the literature for the four flux statistics. The BPDF,
Curvature, FPS and Wavelet statistics are shown in panel A1-A2, B1-B2, C1-C2
and D1-D2, respectively. The joint constraints on thermal parameters from this
work are shown by the black stars in each panel. For a fair comparison, we
also show the measurements of thermal parameters obtained using individual
statistics (i.e., same as Fig. 9, showing only statistical uncertainty). The
uncertainty of our measurements for individual statistics appears to be
consistent with being statistical in nature and is slightly bigger than that
from the joint analysis. For the BPDF statistics, our $T_{0}$ measurements are
typically smaller than that obtained from the literature (panel A1) while the
$\gamma$ measurements are similar to that in the literature. For the curvature
statistics we show the $T_{0}$ measurements obtained using the characteristic
overdensity method ( $\overline{\Delta}$ see §5.4.2, Becker et al., 2011). We
use our $\gamma$ measurements obtained from the joint analysis as a prior to
calculate the $T_{0}$ uncertainty for the $\overline{\Delta}$ method in panel
B1. The $T_{0}$ measurements using the $\overline{\Delta}$ method are
systematically larger than that from the PDF method. However, the $T_{0}$ in
both measurements are in good agreement (maximum difference $\sim 1.1\sigma$)
with each other. Our $T_{0}$ measurements at $z\geq 2.6$ using the
$\overline{\Delta}$ method are within $1\sigma$ of that from Becker et al.
(2011); Boera et al. (2014). Panel C1 and C2 shows that the $T_{0}$ ($\gamma$)
measurements for the FPS statistics are systematically larger (smaller) than
those from Walther et al. (2019). The systematic differences in $T_{0}$ and
$\gamma$ are likely due to constraints on pressure smoothing in Walther et al.
(2019, see §5.4.3). However, our $T_{0}$ measurement using the FPS statistics
is within $1\sigma$ of the corresponding measurement by Walther et al. (2019).
Panel D1 and D2 show the comparison of thermal parameters from this work with
that from Lidz et al. (2010); Garzilli et al. (2012) using wavelet statistics.
The improvement in the uncertainty of $T_{0}$ and $\gamma$ is due to our
larger observational sample and the simultaneous measurements from multiple
wavelet scales (see §5.4.4). Similar to the curvature statistics, we also show
the $T_{0}$ measurements using the $\overline{\Delta}$ method (see §5.4.4).

Fig. 11 shows the evolution of $T_{0}$ and $\gamma$ (as given in Table 2 and
Table 5) from the joint analysis with total uncertainty in the redshift range
$2\leq z\leq 4$. Initially $T_{0}$ ($\gamma$) is small (large) at $3.7\leq
z\leq 3.9$ then it increases (decreases) attaining a maximum (minimum) at
$2.7\leq z\leq 3.1$. Subsequently $T_{0}$ ($\gamma$) decreases (increases) at
$z\leq 2.7$. As we will discuss in §6, the peak in the evolution of $T_{0}$
and $\gamma$ is due to the additional heating from He ii reionization. Such a
peak in $T_{0}$ evolution has been suggested based on other published
measurements, but note the rather large scatter between different measurements
(see Fig. 11 and Fig. 31). Our measurements show for the first time a well
defined peak and the expected smooth evolution of $T_{0}$ and $\gamma$ in the
redshift range $2.0\leq z\leq 3.8$, consistently for all four different flux
statistics. We attribute this to the following main differences of our study:
(i) our observed sample contains a larger number of QSO sightlines than
previous studies, (ii) the simulated Ly$\alpha$ forest in our analysis is
generated for a finely sampled $T_{0}-\gamma$ grid than in previous studies,
(iii) we treat the data and simulations on an exact equal footing i.e., our
simulated mock spectra are mimicked to match the observed sample, (iv) we
apply the same procedure/algorithm to calculate the Ly$\alpha$ flux and Voigt
profile statistics to observed and simulated spectra (v) by combining the
$T_{0}-\gamma$ measurements from different statistics, we mitigate possible
biases of the individual flux statistics.

###  5.4 Comparison with other measurements

Table 3: Characteristic density probed by curvature and wavelet statistics in our simulations  $z$ | Curvature $\overline{\Delta}$ | Wavelet $\overline{\Delta}$ | Wavelet $\overline{\Delta}$
---|---|---|---
|  |  $s_{\rm n}=50$ km s-1 |  $s_{\rm n}=100$ km s-1
2.0 | 4.85 | 4.95 | 5.37
2.2 | 4.57 | 4.30 | 4.98
2.4 | 3.90 | 3.89 | 4.37
2.6 | 3.54 | 3.56 | 4.01
2.8 | 3.19 | 3.26 | 3.50
3.0 | 2.74 | 2.87 | 3.02
3.2 | 2.48 | 2.53 | 2.67
3.4 | 2.17 | 2.16 | 2.27
3.6 | 2.00 | 1.90 | 2.25
3.8 | 2.08 | 1.90 | 2.72

We now compare the $T_{0}$ and $\gamma$ measurements from this work with that
from the literature. Fig. 11 shows a compilation of $T_{0}$ and $\gamma$
measurements from various papers. $T_{0}$ and $\gamma$ have been measured in
the past (with variants) of the same four different statistics we have used
here. In order to do a fair comparison, we compare our measured $T_{0}$ and
$\gamma$ from the joint analysis and the individual statistics with the
corresponding statistics from the literature.

####  5.4.1 Comparison with measurments from $b-\log{\rm N_{HI}}$

In panel A1 and A2 of Fig. 11, we compare our thermal parameter evolution with
that from Schaye (2001); Bolton et al. (2014); Rorai et al. (2018); Hiss et
al. (2018, 2019); Telikova et al. (2019) obtained using the $b-\log{\rm
N_{HI}}$ distribution. This method is based on fitting the Ly$\alpha$ forest
spectra with Voigt profiles and defining a lower envelope of the 2D
$b-\log{\rm N_{HI}}$ distribution. The motivation behind this is that for the
absorption components with low $b$ the broadening should be dominated by
thermal broadening.

Our $T_{0}$ and $\gamma$ measurements from this statistics are consistent
within $1\sigma$ with that of Bolton et al. (2014); Rorai et al. (2018). The
$\gamma$ measurements by Hiss et al. (2018); Telikova et al. (2019) are in
good agreement with our $\gamma$ measurements. However, their $T_{0}$
measurements are significantly higher than our measurements at $2.6\leq z\leq
3.6$. Fitting a lower envelope to the $b-\log{\rm N_{HI}}$ distribution
although well motivated, is somewhat subjective. Finding a lower envelope in
an objective way is difficult (but see Telikova et al., 2019; Hiss et al.,
2019). The lower envelope is found by iteratively rejecting the lines with low
$b$ parameter until convergence is achieved. More (less) rejection of low $b$
parameter lines can lead to systematically higher (lower) inferred
temperatures. However, the high values of $T_{0}$ of Hiss et al. (2018);
Telikova et al. (2019) obtained from the $b-\log{\rm N_{HI}}$ distribution are
not consistent with the measurements obtained using e.g. the FPS by Walther et
al. (2019).

Manual Voigt profile fitting is labour intensive and the number of Voigt
components obtained from observations is rather limited in many analyses and
often the $b-\log{\rm N_{HI}}$ plane is sampled poorly. To avoid these
systematic biases, we use the 1D $b$ distribution (BPDF) computed in different
$\log{\rm N_{HI}}$ bins. The BPDF does not rely on rejecting low $b$
parameters to find a lower envelope. At the same time the BPDF statistics is
sensitive to both $T_{0}$ and $\gamma$ (see appendix E.4 for details). Our
method of measuring $T_{0}$ and $\gamma$ using Voigt profile parameters is
thus significantly different from that used in the literature. Unlike previous
work, our $T_{0}$ and $\gamma$ measurements using the BPDF statistics are
consistent with our measurements using the other statistics as shown in Fig.
9.

####  5.4.2 Comparison with measurments from curvature statistics

Figure 12:  The figure shows the existence of a characteristic overdensity (
$\overline{\Delta}$) in our simulations probed by the curvature statistics
(left panel) and wavelet statistics for two wavelet scales (50 km s-1 and 100
km s-1 for middle and right panels respectively). Simulations with different
$\gamma$ are shown by different colors. For visualisation purposes, we show
only points for selected $\gamma$ values. In the actual analysis, we use all
the simulated $\gamma$ values of our finely sampled $T_{0}-\gamma$ grids. The
red curve in each panel shows the one-to-one relation between mean curvature
(or wavelet) amplitude and temperature at $\overline{\Delta}$. The one-to-one
relation is obtained by fitting a power-law. The higher wavelet scale
corresponds to slightly higher characteristic density indicating that
$\overline{\Delta}$ depends on the scales probed by the statistics. The
$\overline{\Delta}$ probed by the curvature statistics is similar to that
probed by the wavelet statistics for $s_{\rm n}=50$ km s-1. The redshift
evolution of $\overline{\Delta}$ from this work (see Table 3) is in good
agreement with that obtained by Becker et al. (2011) for the curvature
statistics.

In panel B1 and B2 of Fig. 11, we compare our thermal parameter evolution with
that from Becker et al. (2011); Boera et al. (2014) obtained using the
curvature statistics. These authors measured $T_{0}$ by measuring the
temperature at a characteristic density $\overline{\Delta}$ (
$T(\overline{\Delta})$) that is a one-to-one empirical function of the mean of
the absolute curvature irrespective of $\gamma$. To map $T(\overline{\Delta})$
to the temperature at cosmic mean density $T_{0}$, one needs to however assume
a value of $\gamma$.

In Fig. 9 we show the measurement of thermal parameters using the PDF of the
curvature statistics. To make a more direct comparison between our measurement
with that from Becker et al. (2011); Boera et al. (2014), we also calculate
the $T_{0}$ using the characteristic (over-)density method. Fig. 12 and Table
3 shows the $\overline{\Delta}$ probed by the curvature statistics in our
simulations. The evolution of $\overline{\Delta}$ from this work is in good
agreement with that from Becker et al. (2011). We obtain the one-to-one
empirical relation between mean absolute curvature and temperature at
$\overline{\Delta}$ (as shown by the red curve in Fig. 12). For a given
observed mean absolute curvature, we measure the temperature at
$\overline{\Delta}$. We convert the temperature at $\overline{\Delta}$ into
$T_{0}$ using this empirical relation assuming a range of $\gamma$. For this
we use the $\gamma$ values (with 1$\sigma$ uncertainty) from our joint
constraints. The green squares in panel B1 of Fig. 11 show our $T_{0}$
measurements using the $\overline{\Delta}$ method. The errors on the $T_{0}$
measurements for the $\overline{\Delta}$ method account for the uncertainty in
$\gamma$ only.

In Fig. 11, we also show the $T_{0}$ measurements using the
$\overline{\Delta}$ method and that using the joint analysis. The best fit
$T_{0}$ values from the $\overline{\Delta}$ method are systematically higher
than that from the joint analysis as well as those obtained using the
curvature PDF (see Fig. 9) at $z<3.4$. The $\overline{\Delta}$ method assumes
that most of the absorption in the Ly$\alpha$ forest at any given redshift is
associated with densities $\sim$ $\overline{\Delta}$. When converting
temperature at $\overline{\Delta}$ to $T_{0}$, one assumes a single
$\overline{\Delta}$ value. However, in a realistic scenario the Ly$\alpha$
absorption will arise from a range in densities. Unlike the
$\overline{\Delta}$ method, the $T_{0}$ measured in our joint analysis (using
PDFs and PS) is contributed by Ly$\alpha$ absorption coming from all
densities. We suspect this to be the most likely reason for the systematically
larger $T_{0}$ obtained with the $\overline{\Delta}$ method. Another
interesting feature of the $T_{0}$ evolution obtained using the
$\overline{\Delta}$ method is that the uncertainty of the $T_{0}$ measurements
increases with decreasing redshift. This is mainly because of the uncertainty
in $\gamma$ and the increase of $\overline{\Delta}$ at lower redshifts.

At $z>3.4$, our $T_{0}$ measurements from the $\overline{\Delta}$ method are
consistent with that from Becker et al. (2011) irrespective of the assumed
value of $\gamma$. This is because the Ly$\alpha$ forest at these redshifts is
sensitive to densities close to the cosmic mean density $\Delta\sim 1$. This
is reflected in the smaller uncertainty of the $T_{0}$ measurements at $z>3$.
As $\gamma$ decreases ($1.2<\gamma<1.5$) in our joint measurements at $2.6\leq
z\leq 3.4$, $T_{0}$ increases and is in good agreement with that from Becker
et al. (2011); Boera et al. (2014) for $\gamma=1.3$. At $z<2.6$, $\gamma$
increases from 1.2 to 1.5, and our $T_{0}$ measurement using the
$\overline{\Delta}$ method is consistent with that from Becker et al. (2011);
Boera et al. (2014) for $\gamma=1.5$.

####  5.4.3 Comparison with measurements from the FPS

Panel C1 and C2 compares our $T_{0}$ measurements (FPS and joint analysis)
with those from Walther et al. (2019) using the FPS statistics. Similar to
this work, Walther et al. (2019) measured the $T_{0}$ and $\gamma$ evolution
using the FPS calculated from the kodiaq DR2 sample. However, their sample
selection criteria is different from ours. Despite this the FPS from our work
is in good agreement with that from Walther et al. (2019). Fig. 11 shows that
our $T_{0}$ measurements are also in good agreement (within $1\sigma$
uncertainty) with their measurements. Note, however, that our $T_{0}$
measurements are systematically larger (albeit within the errors) than their
measurements at $2\leq z\leq 3$. The difference in the $\gamma$ evolution is
more pronounced. Our $\gamma$ evolution shows a minimum at $z=3$ while their
measurement shows a nearly flat $\gamma\sim 1.7$ at $2\leq z\leq 3.5$.

Since the observed FPS are in good agreement with each other, the moderate
differences in the $T_{0}-\gamma$ measurements are likely to come from
differences in the simulations/analysis. Walther et al. (2019) vary $T_{0}$,
$\gamma$ and the pressure smoothing scale $\lambda_{p}$ as a free parameter
using the OH17 UVB. They derive a single pressure smoothing scale at a given
redshift for the entire simulation box using a cutoff in the 3D real space
flux power spectrum that is calculated without accounting for peculiar
velocity and thermal broadening. In reality, the pressure smoothing scale is
not an independent parameter. This is because the pressure smoothing scale is
set by temperature and density of the gas particles i.e.,
$\lambda_{p}\propto(T/\Delta)^{1/2}\propto[T_{0}\>\Delta^{\gamma-2}]^{1/2}$.
Simultaneous fitting of $T_{0},\gamma$ and $\lambda_{p}$ can result in
degeneracies of $\lambda_{p}-T_{0}$ (see Fig. 6, Fig. 6 in Rorai et al., 2018;
Walther et al., 2019, respectively). An overestimation of $\lambda_{p}$ can
result in a systematic decrease and increase in $T_{0}$ and $\gamma$,
respectively. In our simulations, the pressure smoothing scale is set by the
density and temperature of each SPH particle which varies for each particle.
Hence effectively we fit $T_{0}-\gamma$ as two free parameters to match the
simulations with observations. It is also important to note that all UVB
models predict $\gamma<1.55$ at $2\leq z\leq 4$ (for equilibrium as well as
non-equilibrium evolution of the ionization). As we show in §6 the $\gamma$
evolution in Walther et al. (2019) is very difficult to reproduce with any UVB
model.

####  5.4.4 Comparison of measurements from the wavelet statistics

We compare our $T_{0}$ and $\gamma$ measurements using the wavelet statistics
with those from Lidz et al. (2010) and Garzilli et al. (2012) in panel D1 and
D2 of Fig. 11 respectively. Our $T_{0}-\gamma$ measurements from the wavelet
analysis are more consistent with other flux statistics than those of Lidz et
al. (2010) and Garzilli et al. (2012). The main difference between previous
works and this work are the observations and the methodology we use. Our
observed sample size ($N_{\rm spectra}=103$) is larger than that of previous
studies using wavelet statistics (40 and 20 for Lidz et al. (2010) and
Garzilli et al. (2012) respectively). To mitigate noise effects, Lidz et al.
(2010) apply boxcar smoothing on the wavelet field which somewhat reduces the
sensitivity of the wavelets to $T_{0}$ and $\gamma$. In this work, we compute
the PDF from the unsmoothed wavelet field using kernel density estimation
which converges faster than using a simple histogram reducing the effect of
noise on the PDF. In addition we simultaneously measure $T_{0}$ and $\gamma$
at wavelet scales $s_{\rm n}=30,40,\cdots,90,100$ km s-1. Lidz et al. (2010)
and Garzilli et al. (2012) use only one wavelet scale at a time to measure
$T_{0}$ and $\gamma$. As a result, our $T_{0}$ and $\gamma$ constraints are
tighter than that of Lidz et al. (2010) and Garzilli et al. (2012).

We find the characteristic over-density ( $\overline{\Delta}$) probed by the
wavelet statistics for $s_{\rm n}=50$ km s-1 and $s_{\rm n}=100$ km s-1,
respectively (see Fig. 12). Table 3 also summarizes the redshift evolution of
$\overline{\Delta}$ for these two wavelet scales. The characteristic density
varies similarly with other wavelet scales. We find that the scatter in the
empirical relation increases with increasing wavelet scales reducing the
sensitivity to $T_{0}$. This is expected as on large scales the baryons trace
the dark matter density field and large scales are less sensitive to
variations in $T_{0}$. We also measure $T_{0}$ using the $\overline{\Delta}$
method for the wavelet statistics assuming the $\gamma$ evolution from our
joint analysis. The $T_{0}$ values for all the wavelet scales are very
similar. In Fig. 11, we show the $T_{0}$ values obtained using the
$\overline{\Delta}$ method (green squares, for $s_{\rm n}=50$ km s-1) and that
obtained using the joint analysis (black stars). The $T_{0}$ values from the
two methods are in good agreement with each other (within $1\sigma$).

##  6 Discussion: Implication for UVB models and HeII reionization

Figure 13:  The figure shows the evolution of thermal parameters obtained from
our joint analysis (black star points) and the predictions of the uniform UVB
models of ; ; ; ;  (different curves). Panels A1 and A2 show the $T_{0}$,
$\gamma$ evolution obtained for UVB models assuming equilibrium ionization
evolution. The observed $T_{0}$ and $\gamma$ evolution is significantly
different from those predicted by the uniform UVB models assuming ionization
equilibrium evolution. Panels B1 and B2 show the predicted $T_{0}$, $\gamma$
evolution of the UVB models assuming (physically very well motivated) non-
equilibrium ionization evolution. The $T_{0}$ evolution for the ; ;  non-
equilibrium models at $z<3$ are consistent with that from observations. As we
will discuss below the mismatch between observed and non-equilibrium UVB
models at $z>3$ is likely due to a somewhat wrong evolution of the assumed He
ii ionizing emissivity and/or (effective) He ii photo-heating rates. Note also
that on the rising side of the temperature peak He ii reionization will not
yet be complete and the UVB is not expected to be spatially homogeneous.The
$\gamma$ evolution of the ;  non-equilibrium UVB models is in good agreement
with the observed values. However $T_{0}$ is systematically higher in these
models. We have thus run the models of ;  also with reduced photo-heating
rates. The observed $T_{0}$ evolution is well reproduced if the photo-heating
rate in these two models is lowered by factors of 0.8 and 0.7, respectively
(see Fig. 28 and 29). The good match between observation and non-equilibrium ;
UVB models with rescaled photo-heating rates is due to the rapid evolution of
the He ii ion fraction in these two models.  Figure 14:  Panel A1 and A2 are
similar to panel B1 and B2 in Fig. 13 except that the photo-heating rates of
the ; ;  UVB models are scaled by a factors of 0.8, 0.9 and 0.7, respectively.
Note that such a scaling with a constant factor does not affect $\gamma$.
Panel B1 and B2 shows the evolution of the volume averaged He ii and H i
fraction in our simulations for all the UVB models. $f_{\rm HeII}$ evolves
rather rapidly at $2.9\leq z\leq 3.5$ in the ;  models (shaded region in panel
B1 and B2). Correspondingly $f_{\rm HI}$ in the same redshift range is similar
for all UVB models (except P19).

In this section, we compare our measurements of $T_{0}$ and $\gamma$ with
predictions from published spatially homogeneous UVB models for the effect of
He ii reionization on the thermal evolution of the IGM. In Fig. 13, we compare
the observed $T_{0}-\gamma$ evolution with that predicted by the UVB models of
; ; ;  and FG20. We use cite to obtain the evolution of $T_{0}$ and $\gamma$
for these UVB models, post-processing gadget-3 outputs.

###  6.1 Equilibrium vs Non-equilibrium UVB models

Panels A1 and A2 in Fig. 13 show that the different UVB models with
equilibrium ionization evolution show a very similar $T_{0}-\gamma$ evolution
in the redshift range $2\leq z\leq 4$ irrespective of the temperature after H
i reionization. The observed $T_{0}$ ($\gamma$) is systematically larger
(smaller) than that predicted from equilibrium models at $2.2\leq z\leq 3.4$ (
$2.2\leq z\leq 3.6$). It appears that it is not possible to reproduce the
higher $T_{0}$ (smaller $\gamma$) seen in the observed data with UVB models
assuming equilibrium ionization evolution. With equilibrium ionization
evolution, one assumes that the gas is ionized instantaneously as the ionizing
radiation field is turned on. Because of this, the He ii fraction in
equilibrium model is systematically underestimated (Puchwein et al., 2015;
Gaikwad et al., 2019). Since the photo-heating is proportional to the He ii
fraction, the temperature is underestimated. Similarly, the He ii fraction
depends on density in photo-ionization equilibrium and as a result the photo-
heating of the gas is density dependent. As shown in Gaikwad et al. (2019) the
density dependent photo-heating rates result in $\gamma\sim 1.5$ in
equilibrium ionization models.

The He ii reionization will in reality also be spatially inhomogeneous, and
spatially homogeneous UVB models can only give an average evolution of the
thermal history. However, in any given place the He iii ionization fronts will
move too fast for equilibrium ionization evolution to be a good approximation.
We have thus also implemented physically well motivated non-equilibrium
ionization evolution in cite. Panels B1 and B2 in Fig. 13 show the comparison
of our $T_{0}$ and $\gamma$ measurements with the prediction from UVB models
assuming non-equilibrium evolution. The $T_{0}$ ($\gamma$) values are
systematically larger (smaller) than in the corresponding equilibrium models
and the differences between the different UVB models are larger. The non-
equilibrium predictions by all the UVB models also show a more pronounced
temperature peak. This peak is generally somewhat higher than that in our
measurements from the observations and occurs in some of the models also
somewhat earlier.

###  6.2 The temperature peak/rise due to HeII reionization predicted by
spatially homogeneous UVB models

Table 4: Redshift and extent of He ii reionization in UVB models

UVB Model |  $z_{\rm mid}$[a] |  $\Delta z$[b] |  $f_{\rm scale}$[c] |  $u_{0}(z=3)$[d]
---|---|---|---|---
FG20 | 3.2 | 0.3 | 0.70 $\pm$ 0.10 | 2.97
P19 | 3.5 | 0.8 | 0.90 $\pm$ 0.20 | 3.05
KS19 | 3.6 | 0.8 | 1.00 $\pm$ 0.20 | 2.43
HM12 | 4.0 | 1.0 | 1.00 $\pm$ 0.20 | 3.14
OH17 | 3.2 | 0.4 | 0.75 $\pm$ 0.15 | 3.09

  * a

The redshift of reionization is defined as $z_{\rm mid}$ corresponding to
epoch when $f_{\rm HeII}=0.5$

  * b

The extent of reionization is defined as $\Delta z=z_{\rm start}-z_{\rm end}$.
where $z_{\rm start},\>z_{\rm end}$ corresponds to the redshift when $f_{\rm
HeII}=0.75,0.02$ respectively.

  * c

Factor by which photo-heating rates are scaled to better match the observed
$T_{0}$ evolution at $z<3$. Errors represent the error in the scaling factor
that corresponds to the $1\sigma$ uncertainty for $T_{0}$ (see Fig. 28 and
Fig. 29).

  * d

Cumulative energy deposited (in ${\rm eV}\>m^{-1}_{p}$) into the IGM at $z=3$
for non-equilibrium UVB models without scaling the photo-heating rates (see
appendix I).

The basic concept of the different UVB models is similar, but there are
considerable differences in the detailed assumptions. ;  compute photo-
ionization and photo-heating rates during ( H i and He ii) reionization
differently than ; ; . To calculate photo-ionization and photo-heating rates
during He ii reionization, ;  assume the He ii reionization history ( $f_{\rm
HeII}$ evolution) to be a free parametric function. The reionization history
is calibrated to match the observed $\tau_{\rm eff,HeII}$ evolution (Worseck
et al., 2019). To compute the photo-heating rates during reionization, ;  then
use another free parameter namely the total heat injected during He ii
reionization ($\Delta T_{\rm HeII}$). ; ;  on the other hand try to predict
photo-ionisation and photo-heating rates based on observed source luminosity
functions and SEDs and a model of the He ii opacity based on observations.
There is thus considerable variety in the assumptions of spatially homogeneous
UVB models. As a result the exact time and height of the temperature peak
predicted by these models should be considered somewhat uncertain. As can be
seen in panel B1/B2 in Fig. 13, the timing of the temperature peak in the ;
UVB models agrees very well with the timing of the peak in our measurements.
The evolution of $\gamma$ in these models is also in better agreement with our
measurements than those from ; ; . However, $T_{0}$ is systematically larger
than our measured values at all redshifts. Given the uncertainty in the
amplitude of the photo-heating rates in all the UVB models, we have rescaled
the H i, He i and He ii photo-heating rates in order to better match our
measured evolution of the thermal parameters. Such a rescaling changes $T_{0}$
systematically, while the $\gamma$ evolution remains largely unaffected.

Panels A1 and A2 in Fig. 14 show the evolution of $T_{0}$ and $\gamma$ for
scaled photo-heating rates in the different UVB models. In Table 4, we
summarize the photo-heating rate scale factor required in different UVB models
to match the observed $T_{0}$ evolution at $z<3$ (see also Fig. 28 and Fig. 29
for details). As expected, all the rescaled UVB models show a better agreement
with the $T_{0}$ evolution at $z<3$. However, the $T_{0}$ evolution at $z>3$
in the rescaled ;  UVB models is in much better agreement with our measured
values than the ; ;  UVB models. The shape of the $T_{0}-\gamma$ evolution is
closely related to the He ii reionization history ( $f_{\rm HeII}$ evolution
Gaikwad et al., 2019). It is noteworthy here that the reionization history
used in the ;  UVB models is calibrated to match the observed $\tau_{\rm
eff,HeII}$ evolution from Worseck et al. (2019). The good match between the
shape of the observed $T_{0}-\gamma$ evolution and that predicted by the UVB
models with scaled photo-heating rates at $z>3$ suggests that the
$T_{0}-\gamma$ evolution from our work is in good agreement with the observed
$\tau_{\rm eff,HeII}$ evolution from Worseck et al. (2019). Such consistency
is important because the evolution of $T_{0}$ and $\gamma$ at $2\leq z\leq 4$
is mainly driven by the $f_{\rm HeII}$ evolution (see Fig. 5 in Gaikwad et
al., 2019). We should here emphasize that the observed dataset, methodology
and parameters that we used here are completely complementary to those used to
measure the $\tau_{\rm eff,HeII}$ evolution in Worseck et al. (2019).

Despite the good match between observations and the predictions by non-
equilibrium ;  UVB models, we would like to caution the reader that these UVB
models do not account for inhomogeneous He ii reionization and the spatial
fluctuations in the He ii ionizing background which are expected to be
present. To capture the effects of such spatial fluctuations, will, however
require rather challenging high dynamic range radiative transfer simulations
(ideally coupled to the hydrodynamics). We hope to do this in future work but
caution that it will be very expensive to calibrate such simulations to match
the observed data (see Gaikwad et al., 2020).

###  6.3 Late vs Early HeII reionization

In Fig. 14, we compare the evolution of the He ii and H i fraction predicted
by the different UVB models with non-equilibrium ionisation evolution. As
already discussed the amount of heat injected into the IGM is coupled to the
$f_{\rm HeII}$ evolution and we have scaled the photo-heating rates to match
the observed $T_{0}$ evolution. To quantify the differences between the
various UVB models, we define the mid point of He ii reionization as the
redshift, where $f_{\rm HeII}=0.5$. We also define the extent of reionization
($\Delta z$) as corresponding to the difference between the redshifts when
$f_{\rm HeII}=0.75$ and $f_{\rm HeII}=0.02$. These definitions facilitate the
comparison of the UVB models quantitatively, but note that they may differ
from definitions used in the literature (see Table 4).

The midpoint of He ii reionization in the ;  models ($z_{\rm mid}\sim 3.2$) is
at lower redshift than that for the models by ; ; . The observed
$T_{0}-\gamma$ evolution from our measurements is consistent with the
relatively late He ii reionization in the ;  models. Fig. 14 also shows that
the evolution of $f_{\rm HeII}$ in the ; ;  models is more gradual than that
in the ;  UVB models. The extent of He ii reionization $\Delta z$ in the ;
UVB models $(<0.4)$ is smaller than that in the ; ;  UVB models (see Table 4).
We should note here that the actual duration and heating due to He ii
reionization may be different if the spatially inhomogeneous nature of He ii
reionization and the resulting fluctuations in the He ii ionizing background
are taken into account. With this caveat in mind we note that for the
spatially uniform UVB models considered in this work, our measured
$T_{0}-\gamma$ evolution is consistent with late and rapid He ii reionization.

We refer readers to appendix H and I for more details on the effect of the UVB
on the evolution of thermal parameters and cumulative energy deposited in to
the IGM. In particular, we show the effect of observational and modeling
uncertainty in the UVB model on thermal parameters in appendices H.1 and H.2.
We discuss the likely modification needed for the physically motivated non-
equilibrium UVB models to reproduce our measured thermal parameter evolution
in appendices H.3 and H.4. We discuss the evolution of the cumulative energy
deposited in IGM by various UVB models in appendix I.

Finally, we emphasize again that even though we use equilibrium KS19 UVB
models to vary the thermal parameters in the simulations, our final measured
thermal parameter evolution is consistent with the predictions by the
physically motivated non-equilibrium evolution for the ;  UVB models with
moderately reduced photo-heating rates. This gives us additional confidence
that our measurements of thermal parameters are robust with regard to
observational and modeling systematics.

##  7 Summary

We have measured the thermal parameters $T_{0}$ and $\gamma$ of the IGM at
$2\leq z\leq 4$ using 103 high-resolution, high signal-to-noise QSO absorption
spectra drawn from the kodiaq DR2 sample using four different Ly$\alpha$ flux
statistics. The measurements are calibrated with a suite of simulations based
on a high-resolution smoothed particle hydrodynamics simulation post-processed
with our IGM thermal history code “Code for Ionization and Temperature
Evolution” (cite) that finely samples the grid in $T_{0}$ and $\gamma$. The
salient points of our work are:

  * •

We use a sample of 103 QSO spectra satisfying the following criteria; (i)
spectral S/N $>5$, (ii) no large spectral gaps and (iii) no identified DLAs or
sub-DLAs to avoid systematic biases. We have divided our sample in redshift
bins having $\Delta z=0.2$ with centers at $z=2.0,2.2,\cdots,3.6~{}{\rm
and}~{}3.8$. The typical number of lines of sight used per redshift bin is
larger than what has been used in previous studies. We have identified metal
lines contaminating the Ly$\alpha$ forest and replaced them with continuum and
added noise. Our simulations, based on L10N512 gadget-3 SPH simulation, sample
$T_{0}$ from $\sim 6000$ K to $24000$ K in bins of $\Delta T_{0}=500$ K and
$\gamma$ from $\sim 0.7$ to $\sim 2.0$ in bins in of $\delta\gamma=0.05$
creating a densely sampled $T_{0}-\gamma$ grids with $37\times 27=999$ thermal
parameter variations. The mock spectra created from these simulations account
for the redshift path length, finite S/N and instrumental broadening of the
spectrograph of the observed sample.

  * •

Using our code “VoIgt profile Parameter Estimation Routine” (viper), we fit
the observed and simulated Ly$\alpha$ forest spectra with multi-component
Voigt profiles. We derive 6 flux statistics for observed and simulated
spectra: (i) flux power spectrum (FPS), (ii) wavelet statistics (iii)
curvature statistics, (iv) $b$ distribution function (BPDF), (v) transmitted
flux probability distribution function (FPDF) and (iv) H i column density
distribution function (CDDF). The observed FPS, FPDF, CDDF and mean flux
evolution from our sample is in good agreement with measurements in the
literature. We will make the observed statistics available to the community on
request. By varying one parameter at a time, we show the sensitivity of these
statistics to thermal parameters of the IGM.

  * •

We measure the thermal parameters $T_{0}$ and $\gamma$ in the redshift range
$2.0\leq z\leq 3.8$ by quantitatively comparing four of the six flux
statistic, namely FPS, wavelet, curvature and BPDF statistics. The $T_{0}$ and
$\gamma$ measurements using different statistics are consistent with each
other within $1\sigma$ at all redshifts. By combining the four statistics we
thus find joint constraints on $T_{0}$ and $\gamma$ that are more robust and
tighter than the measured values from the individual statistics. We estimate
the total uncertainty on $T_{0}$ and $\gamma$ by accounting for the
uncertainty due to continuum fitting, metal line contamination, cosmological
parameters and modeling uncertainties. We find that our measured $T_{0}$
varies from $9250$ K at $z=3.8$ to $14750$ K at $z=3$ to $9500$ K at $z=2$,
while $\gamma$ evolves from $1.525$ at $z=3.8$ to $1.225$ at $z=3$ to $1.500$
at $z=2$. We robustly detect a maximum in $T_{0}$ and minimum in $\gamma$ that
occur simultaneously at $z\simeq 3$.

  * •

We compare our measurements to the evolution of $T_{0}$ and $\gamma$ predicted
by spatially homogeneous UVB background models by ; ; ; ;  assuming
equilibrium and non-equilibrium ionization evolution. Our temperature
measurements are difficult to reproduce assuming equilibrium ionization
evolution. Our measurements are, however, in good agreement with the thermal
evolution predicted by the UVB models that have been calibrated to match the
observed He ii opacity evolution (; ) if we assume non-equilibrium ionization
evolution and reduce the photo-heating rates in these models by a moderate
factor 0.7-0.8. The timing of the peak in our measured temperatures due to the
extra heat input during He ii reionization appears thus to be nicely
consistent with the timing of He ii reionization inferred from the He ii
opacity evolution and appears to be rather rapid. Our measured temperatures
are also (to varying degree) systematically lower than those in the other
three models. These three models predict He ii reionization to be more
extended and either an onset of the temperature rise (P19) or a timing of the
peak that is too early (; ) to be consistent with our temperature
measurements.

The simulations used to calibrate our temperature measurements assume time
varying but spatially uniform UVB models. While spatial fluctuations in the He
ii ionizing background are clearly important during He ii reionization
($z>3$), simulating these self-consistently is beyond the scope of the current
work. Large suites of high dynamic range, large scale, multi-frequency
cosmological radiative transfer simulations will be required to capture this
effect accurately and to sample the observational uncertainties regarding
relevant properties of the ionizing sources. We hope to perform such
simulations in future. The tests we have performed and our experience with
hydrogen reionization suggest that modeling the thermal evolution with
homogeneous UVB models should nevertheless recover the evolution of the
(median) thermal state of the IGM not only after but also before overlap of He
iii regions reasonably well.

## Acknowledgments

We thank the DiRAC@Durham facility managed by the Institute for Computational
Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The
equipment was funded by BEIS capital funding via STFC capital grants
ST/P002293/1, ST/R002371/1 and ST/S002502/1, Durham University and STFC
operations grant ST/R000832/1. DiRAC is part of the National e-Infrastructure.
Computations in this work were also performed using the Perseus cluster at
IUCAA, the HPC cluster at NCRA and the CALX machines at IoA. Support by ERC
Advanced Grant 320596 ‘The Emergence of Structure During the Epoch of
reionization’ is gratefully acknowledged. PG and MGH acknowledge the support
of the UK Science and Technology Facilities Council (STFC). TRC acknowledges
support of the Department of Atomic Energy, Government of India, under project
no. 12-R&D-TFR-5.02-0700. We thank J. Bolton for making the Sherwood
simulation suite available to our work. The Sherwood simulations were
performed using the Curie supercomputer at the Tre Grand Centre de Calcul
(TGCC), and the DiRAC Data Analytic system at the University of Cambridge,
operated by the University of Cambridge High Performance Computing Service on
behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was
funded by BIS National E-infrastructure capital grant (ST/K001590/1), STFC
capital grants ST/H008861/1 and ST/H00887X/1, and STFC DiRAC Operations grant
ST/K00333X/1. DiRAC is part of the National E- Infrastructure. We thank Vid
Irsic, Ewald Puchwein and Jamie Bolton for useful comments on the manuscript.

## References

  * Almgren et al. (2013) Almgren A. S., Bell J. B., Lijewski M. J., Lukić Z., Van Andel E., 2013, ApJ, 765, 39
  * Avni (1976) Avni Y., 1976, ApJ, 210, 642
  * Becker & Bolton (2013) Becker G. D., Bolton J. S., 2013, MNRAS, 436, 1023
  * Becker et al. (2011) Becker G. D., Bolton J. S., Haehnelt M. G., Sargent W. L. W., 2011, MNRAS, 410, 1096
  * Becker et al. (2013) Becker G. D., Hewett P. C., Worseck G., Prochaska J. X., 2013, MNRAS, 430, 2067
  * Boera et al. (2014) Boera E., Murphy M. T., Becker G. D., Bolton J. S., 2014, MNRAS, 441, 1916
  * Boera et al. (2019) Boera E., Becker G. D., Bolton J. S., Nasir F., 2019, ApJ, 872, 101
  * Bolton & Haehnelt (2007) Bolton J. S., Haehnelt M. G., 2007, MNRAS, 382, 325
  * Bolton et al. (2012) Bolton J. S., Becker G. D., Raskutti S., Wyithe J. S. B., Haehnelt M. G., Sargent W. L. W., 2012, MNRAS, 419, 2880
  * Bolton et al. (2014) Bolton J. S., Becker G. D., Haehnelt M. G., Viel M., 2014, MNRAS, 438, 2499
  * Bolton et al. (2017) Bolton J. S., Puchwein E., Sijacki D., Haehnelt M. G., Kim T.-S., Meiksin A., Regan J. A., Viel M., 2017, MNRAS, 464, 897
  * Calura et al. (2012) Calura F., Tescari E., D’Odorico V., Viel M., Cristiani S., Kim T. S., Bolton J. S., 2012, MNRAS, 422, 3019
  * Calverley et al. (2011) Calverley A. P., Becker G. D., Haehnelt M. G., Bolton J. S., 2011, MNRAS, 412, 2543
  * Carswell et al. (1982) Carswell R. F., Whelan J. A. J., Smith M. G., Boksenberg A., Tytler D., 1982, MNRAS, 198, 91
  * Chang et al. (2012) Chang P., Broderick A. E., Pfrommer C., 2012, ApJ, 752, 23
  * Choudhury et al. (2001) Choudhury T. R., Srianand R., Padmanabhan T., 2001, ApJ, 559, 29
  * Cirelli et al. (2009) Cirelli M., Iocco F., Panci P., 2009, J. Cosmology Astropart. Phys., 2009, 009
  * D’Aloisio et al. (2020) D’Aloisio A., McQuinn M., Trac H., Cain C., Mesinger A., 2020, arXiv e-prints, p. arXiv:2002.02467
  * Faucher-Giguère (2020) Faucher-Giguère C.-A., 2020, MNRAS, 493, 1614
  * Faucher-Giguère et al. (2008) Faucher-Giguère C.-A., Lidz A., Hernquist L., 2008, Science, 319, 52
  * Gaikwad et al. (2017a) Gaikwad P., Khaire V., Choudhury T. R., Srianand R., 2017a, MNRAS, 466, 838
  * Gaikwad et al. (2017b) Gaikwad P., Srianand R., Choudhury T. R., Khaire V., 2017b, MNRAS, 467, 3172
  * Gaikwad et al. (2018) Gaikwad P., Choudhury T. R., Srianand R., Khaire V., 2018, MNRAS, 474, 2233
  * Gaikwad et al. (2019) Gaikwad P., Srianand R., Khaire V., Choudhury T. R., 2019, MNRAS, 490, 1588
  * Gaikwad et al. (2020) Gaikwad P., et al., 2020, arXiv e-prints, p. arXiv:2001.10018
  * Garzilli et al. (2012) Garzilli A., Bolton J. S., Kim T. S., Leach S., Viel M., 2012, MNRAS, 424, 1723
  * Garzilli et al. (2017) Garzilli A., Boyarsky A., Ruchayskiy O., 2017, Physics Letters B, 773, 258
  * Gnedin & Hui (1998) Gnedin N. Y., Hui L., 1998, MNRAS, 296, 44
  * Haardt & Madau (1996) Haardt F., Madau P., 1996, ApJ, 461, 20
  * Haardt & Madau (2012) Haardt F., Madau P., 2012, ApJ, 746, 125
  * Hiss et al. (2018) Hiss H., Walther M., Hennawi J. F., Oñorbe J., O’Meara J. M., Rorai A., Lukić Z., 2018, ApJ, 865, 42
  * Hiss et al. (2019) Hiss H., Walther M., Oñorbe J., Hennawi J. F., 2019, ApJ, 876, 71
  * Hui & Gnedin (1997) Hui L., Gnedin N. Y., 1997, MNRAS, 292, 27
  * Iršič et al. (2017a) Iršič V., et al., 2017a, Phys. Rev. D, 96, 023522
  * Iršič et al. (2017b) Iršič V., Viel M., Haehnelt M. G., Bolton J. S., Becker G. D., 2017b, Phys. Rev. Lett., 119, 031302
  * Khaire (2017) Khaire V., 2017, MNRAS, 471, 255
  * Khaire & Srianand (2015) Khaire V., Srianand R., 2015, ApJ, 805, 33
  * Khaire & Srianand (2019) Khaire V., Srianand R., 2019, MNRAS, 484, 4174
  * Khaire et al. (2016) Khaire V., Srianand R., Choudhury T. R., Gaikwad P., 2016, MNRAS, 457, 4051
  * Khaire et al. (2019) Khaire V., et al., 2019, MNRAS, 486, 769
  * Kim et al. (2004) Kim T. S., Viel M., Haehnelt M. G., Carswell R. F., Cristiani S., 2004, MNRAS, 347, 355
  * Kim et al. (2007) Kim T. S., Bolton J. S., Viel M., Haehnelt M. G., Carswell R. F., 2007, MNRAS, 382, 1657
  * Kim et al. (2013) Kim T. S., Partl A. M., Carswell R. F., Müller V., 2013, A&A, 552, A77
  * Kirkman et al. (2005) Kirkman D., et al., 2005, MNRAS, 360, 1373
  * Kulkarni & Fall (1993) Kulkarni V. P., Fall S. M., 1993, ApJ, 413, L63
  * Kulkarni et al. (2015) Kulkarni G., Hennawi J. F., Oñorbe J., Rorai A., Springel V., 2015, ApJ, 812, 30
  * Kulkarni et al. (2019) Kulkarni G., Worseck G., Hennawi J. F., 2019, MNRAS, 488, 1035
  * Lidz et al. (2007) Lidz A., McQuinn M., Zaldarriaga M., Hernquist L., Dutta S., 2007, ApJ, 670, 39
  * 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
  * Liu et al. (2020) Liu H., Qin W., Ridgway G. W., Slatyer T. R., 2020, arXiv e-prints, p. arXiv:2008.01084
  * Lukić et al. (2015) Lukić Z., Stark C. W., Nugent P., White M., Meiksin A. A., Almgren A., 2015, MNRAS, 446, 3697
  * Lusso et al. (2014) Lusso E., et al., 2014, ApJ, 784, 176
  * Maitra et al. (2019) Maitra S., Srianand R., Petitjean P., Rahmani H., Gaikwad P., Choudhury T. R., Pichon C., 2019, MNRAS, p. 2422
  * Maitra et al. (2020) Maitra S., Srianand R., Gaikwad P., Choudhury T. R., Paranjape A., Petitjean P., 2020, arXiv e-prints, p. arXiv:2005.05346
  * McDonald (2003) McDonald P., 2003, ApJ, 585, 34
  * Murphy et al. (2019) Murphy M. T., Kacprzak G. G., Savorgnan G. A. D., Carswell R. F., 2019, MNRAS, 482, 3458
  * Nasir et al. (2016) Nasir F., Bolton J. S., Becker G. D., 2016, MNRAS, 463, 2335
  * Nath & Biermann (1993) Nath B. B., Biermann P. L., 1993, MNRAS, 265, 241
  * Oñorbe et al. (2017) Oñorbe J., Hennawi J. F., Lukić Z., 2017, ApJ, 837, 106
  * Oñorbe et al. (2019) Oñorbe J., Davies F. B., Lukić Z., Hennawi J. F., Sorini D., 2019, MNRAS, 486, 4075
  * O’Meara et al. (2015) O’Meara J. M., et al., 2015, AJ, 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, AJ, 154, 114
  * Padmanabhan et al. (2014) Padmanabhan H., Choudhury T. R., Srianand R., 2014, MNRAS, 443, 3761
  * Padmanabhan et al. (2015) Padmanabhan H., Srianand R., Choudhury T. R., 2015, MNRAS, 450, L29
  * Peeples et al. (2010) Peeples M. S., Weinberg D. H., Davé R., Fardal M. A., Katz N., 2010, MNRAS, 404, 1281
  * Planck Collaboration et al. (2014) Planck Collaboration et al., 2014, A&A, 571, A16
  * Press et al. (1992) Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P., 1992, Numerical recipes in FORTRAN. The art of scientific computing
  * Puchwein et al. (2012) Puchwein E., Pfrommer C., Springel V., Broderick A. E., Chang P., 2012, MNRAS, 423, 149
  * Puchwein et al. (2015) Puchwein E., Bolton J. S., Haehnelt M. G., Madau P., Becker G. D., Haardt F., 2015, MNRAS, 450, 4081
  * Puchwein et al. (2019) Puchwein E., Haardt F., Haehnelt M. G., Madau P., 2019, MNRAS, 485, 47
  * Rauch et al. (1997) Rauch M., et al., 1997, ApJ, 489, 7
  * Rollinde et al. (2013) Rollinde E., Theuns T., Schaye J., Pâris I., Petitjean P., 2013, MNRAS, 428, 540
  * Rorai et al. (2017a) Rorai A., et al., 2017a, Science, 356, 418
  * Rorai et al. (2017b) Rorai A., et al., 2017b, MNRAS, 466, 2690
  * 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
  * Samui et al. (2005) Samui S., Subramanian K., Srianand R., 2005, in 29th International Cosmic Ray Conference (ICRC29), Volume 9. p. 215 (arXiv:astro-ph/0505590)
  * Schaye (2001) Schaye J., 2001, ApJ, 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
  * Scoccimarro et al. (2012) Scoccimarro R., Hui L., Manera M., Chan K. C., 2012, Phys. Rev. D, 85, 083002
  * Shull & Danforth (2020) Shull M., Danforth C., 2020, arXiv e-prints, p. arXiv:2007.02948
  * Springel (2005) Springel V., 2005, MNRAS, 364, 1105
  * Stevans et al. (2014) Stevans M. L., Shull J. M., Danforth C. W., Tilton E. M., 2014, ApJ, 794, 75
  * Telikova et al. (2019) Telikova K. N., Shternin P. S., Balashev S. A., 2019, ApJ, 887, 205
  * Theuns et al. (2002) Theuns T., Zaroubi S., Kim T.-S., Tzanavaris P., Carswell R. F., 2002, MNRAS, 332, 367
  * Upton Sanderbeck & Bird (2020) Upton Sanderbeck P., Bird S., 2020, arXiv e-prints, p. arXiv:2002.05733
  * Upton Sanderbeck et al. (2016) Upton Sanderbeck P. R., D’Aloisio A., McQuinn M. J., 2016, MNRAS, 460, 1885
  * Viel et al. (2004a) Viel M., Haehnelt M. G., Springel V., 2004a, MNRAS, 354, 684
  * Viel et al. (2004b) Viel M., Weller J., Haehnelt M. G., 2004b, MNRAS, 355, L23
  * Viel et al. (2017) Viel M., Haehnelt M. G., Bolton J. S., Kim T.-S., Puchwein E., Nasir F., Wakker B. P., 2017, MNRAS, 467, L86
  * Walther et al. (2018) Walther M., Hennawi J. F., Hiss H., Oñorbe J., Lee K.-G., Rorai A., O’Meara J., 2018, ApJ, 852, 22
  * Walther et al. (2019) Walther M., Oñorbe J., Hennawi J. F., Lukić Z., 2019, ApJ, 872, 13
  * Webb & Carswell (1991) Webb J. K., Carswell R. F., 1991, in Shaver P. A., Wampler E. J., Wolfe A. M., eds, Quasar Absorption Lines. p. 3
  * Worseck et al. (2011) Worseck G., et al., 2011, ApJ, 733, L24
  * Worseck et al. (2016) Worseck G., Prochaska J. X., Hennawi J. F., McQuinn M., 2016, ApJ, 825, 144
  * Worseck et al. (2019) Worseck G., Davies F. B., Hennawi J. F., Prochaska J. X., 2019, ApJ, 875, 111
  * Zaldarriaga (2002) Zaldarriaga M., 2002, ApJ, 564, 153

##  Appendix A Coaddition of observed spectra

Figure 15:  Panel A, B and C show the observed spectrum of QSO J095852+120245
at three different epochs having exposure times of 7200 s, 4800 s and 5100 s,
respectively. Many pixels in the individual exposures are bad pixels due to
observational systematics. Panel D shows the coadded spectrum obtained by our
method as discussed in appendix A. Our coaddition procedure takes care of bad
pixels due to observational systematics.

Most of the QSOs in the kodiaq dr2 sample were observed more than once with
different exposure times. In this section, we discuss our procedure to coadd
the spectra. Let a QSO be observed $N$ times with exposure times
$t_{1}$,$t_{2}$,$\cdots$,$t_{N}$. Let $F_{i}(\lambda)$, $\sigma_{i}(\lambda)$
be the value of the normalized flux and associated error at wavelength
$\lambda$ in the $i^{\rm th}$ observation. The coadded flux $F_{\rm
coadd}(\lambda)$ and coadded variance $\sigma^{2}_{\rm coadd}$ of the flux is
given by,

$\displaystyle F_{\rm coadd}(\lambda)$
$\displaystyle=\frac{\sum\limits_{i=1}^{N}\;t_{i}\>w_{i}\>F_{i}(\lambda)}{\sum\limits_{i=1}^{N}\;t_{i}\>w_{i}}$
(1) $\displaystyle\sigma^{2}_{\rm coadd}(\lambda)$
$\displaystyle=\frac{\sum\limits_{i=1}^{N}\;t_{i}\>w_{i}\>\sigma^{2}_{i}(\lambda)}{\sum\limits_{i=1}^{N}\;t_{i}\>w_{i}}$

where $w_{i}$ is weight for the $i^{\rm th}$ pixel. $w_{i}=0$ for bad pixels
and $w_{i}=1$ for good pixels. For the kodiaq dr2 sample, we define a pixel as
a bad pixel if it satisfies any of the following criteria, (i) pixels with
negative error i.e., $\sigma_{i}(\lambda)\leq 0$, (ii) pixels with
$\sigma_{i}(\lambda)\gg 1$ (i.e., SNR per pixel $\ll 1$) (iii)
$F_{i}(\lambda)=0$ (due to sky background $F_{i}(\lambda)\neq 0$ for good
pixel), (iv) negative outliers i.e., pixels with
$F_{i}(\lambda)\leq-3\sigma_{i}(\lambda)$ and (v) positive outliers i.e.,
pixels with $F_{i}(\lambda)\geq 1+3\sigma_{i}(\lambda)$,

Fig. 15 shows an example of coadded spectrum obtained using our method.
Spectrum obtained in individual exposures show many bad pixels due to
observational systematics. Our coaddition method accounts for these effects as
shown in panel D of Fig. 15.

##  Appendix B Excess mean flux $\langle F\rangle$ at $z=3.2$

Figure 16:  Panels A to D show the observed mean flux along individual QSO
sightlines in different redshift bins. The dashed line in each panel shows the
mean flux of the sample at the given redshift. We exclude the possibility of
larger mean flux at $3.1\leq z\leq 3.2$ due to any outliers. Panel E shows the
observed FPDF in four redshift bins. The vertical lines in panel E show the
observed mean flux in the corresponding redshift bins. Similar to panel A-D,
the FPDF comparison in panel E does not show any unusual evolution.

In this section we test if the observed excess in mean flux at $3.1\leq z\leq
3.3$ in the kodiaq dr2 sample is due to the presence of any unusual QSO
sightlines. To check if there are any outliers in this particular redshift
bin, we plot the mean flux along individual sightlines in Fig. 16 for four
redshift bins. Fig. 16 illustrates that the distribution of mean flux along
sightlines is random around the mean flux (dashed line in panel A to D) of the
entire sample for all redshift bins. We also do not see any outlier at
$3.1\leq z\leq 3.3$. Since the mean flux along a sightline is an average
quantity, we also check if there is any significant deviation in the observed
FPDF evolution in panel E of Fig. 16. The observed FPDF at $3.1\leq z\leq 3.3$
also does not show any unusual deviation from other redshift bins. The excess
in $\langle F\rangle$ at $3.1\leq z\leq 3.3$ may thus be real in the kodiaq
dr2 sample. Another reason for the excess $\langle F\rangle$ could be that the
number of sightlines at $3.1\leq z\leq 3.3$ are limited in our sample. Our
preliminary analysis using UVES squad dr1 (Murphy et al., 2019) while showing
distinct change in the slope at z$\sim 3.2$ does not show any prominent
excess. We plan to investigate this further in our future work. However, we
emphasize that the effect of excess in $\langle F\rangle$ at $3.1\leq z\leq
3.3$ has a marginal effect on our measurements of thermal parameters. This is
because for any thermal parameter variation we rescale the optical depth in
our simulations to match the mean flux. Thus, even if there are any
systematics in the $\langle F\rangle$ evolution in the kodiaq dr2 sample, this
does not significantly change the main results of this work.

##  Appendix C Comparison of observed flux statistics errors

Figure 17: Each panel is same as the corresponding panel in Fig. 3 except that
the errors of the observed FPDF, FPS and CDDF statistics are shown. For visual
purposes, the errors of the FPDF, FPS and CDDF statistics are inflated by a
factor of 2, 2 and 4 respectively. For a fair comparison, the errors from our
sample are calculated using the bootstrap method. In general the errors of the
FPDF, FPS and CDDF statistics from our sample are smaller because of the large
sample size. The errors of the FPDF in Rollinde et al. (2013) are computed
from simulated mock spectra and hence are larger than that of the other
measurements. The errors of the FPS statistics from this work are similar to
those in Walther et al. (2018, for the no metal case) due to similar sample
size and data quality. The errors of our CDDF at $\log{\rm N_{HI}}$ $<13.6$
are larger due to the somewhat smaller S/N of the spectra in our sample. At
$\log{\rm N_{HI}}$ $>13.6$, the errors of our CDDF are smaller than that in
Kim et al. (2013).

In §2.1 we discuss the consistency of observed flux statistics from this work
with that in the literature. In this section, we compare the errors of the
observed FPDF, FPS and CDDF from this work with measurements in the
literature. The errors for observed statistics are normally determined using
the bootstrap method (except in Rollinde et al., 2013). For a fair comparison,
we also compute the errors using the bootstrap method. Fig. 17 shows a
comparison of the errors for different statistics. In general, the errors in
different statistics from this work are either consistent or slightly smaller
than those in the literature.

Panel A1 in Fig. 17 shows that the error in our observed FPDF is smaller than
that in Kim et al. (2007); Calura et al. (2012); Rollinde et al. (2013). The
smaller errors are due to the large number of QSOs in our sample. The large
error in the FPDF of Rollinde et al. (2013) is due to the difference in their
method of calculating the errors. Rollinde et al. (2013) argue that the errors
computed using the bootstrap method are underestimated and may not account for
the cosmic variance accurately. Hence they derive the errors from large number
of simulated mocks. Furthermore, their sample size was limited to $<5$ at the
redshift of interest. We have calculated the errors on our other observed
statistics using both methods and accounted for cosmic variance (see appendix
E.5 for details).

Panel B1 in Fig. 17 shows a comparison of the error of the FPS from this work
with that from Walther et al. (2018). The errors (for the no metal case) are
similar in the two works. This is expected because the number of QSOs per
redshift bin are similar in both cases. Similarly, the data quality is the
same because the subsets of QSOs are both drawn from the kodiaq dr2 survey.
The consistency in FPS and its error of this work with that from Walther et
al. (2018) reflects the consistency of both the data and our method of
calculating the FPS. This consistency is important because the FPS is one of
the statistics used in this work to measure thermal parameters.

In panel C1 of Fig. 17, we compare the errors of our CDDF with that from Kim
et al. (2013). Kim et al. (2013) measured the CDDF from 18 high resolution
high S/N (usually $>50$) Ly$\alpha$ forest spectra. On the other hand we
measure the CDDF from a sample of $\sim 35$ high resolution but moderate S/N
spectra. Thus the data size and data quality is different in both the works.
These differences of the observed dataset are reflected in the CDDF errors.
For example the errors in our CDDF at $\log{\rm N_{HI}}$ $<13.6$ is larger
than that from Kim et al. (2013). This is because low ${\rm N_{HI}}$ systems
are affected more by the S/N of the spectra. Since S/N of our sample is
smaller than in Kim et al. (2013), our CDDF error at small $\log{\rm N_{HI}}$
is correspondingly large. At high $\log{\rm N_{HI}}$ $>13.6$, S/N is no longer
an issue since lines are detected with relatively high significance level. The
errors in our CDDF at $\log{\rm N_{HI}}$ $>13.6$ are smaller than that from
Kim et al. (2013) because the number of spectra in our sample are nearly twice
that in Kim et al. (2013). Even though there are differences in data quality
and sample size, the CDDF and the associated errors from the two works are in
good agreement with each other. The agreement between the CDDF demonstrates
the consistency of viper in fitting the Ly$\alpha$ forest with Voigt profiles.
This consistency is important because we use the $b$ parameter of Voigt
profile fits to constrain thermal parameters. In summary, the observed
statistics with associated errors from this work are in good agreement with
that from the literature.

##  Appendix D Variation in thermal parameters

Figure 18: The top and bottom panels show the variation in thermal parameters
$T_{0}$ and $\gamma$ obtained by rescaling the H i, He i and He ii photo-
heating rates of the KS19 UVB with our post-processing module cite. A simple
constant scaling of photo-heating rates (by factor $a$) leads to variation in
$T_{0}$ while $\gamma$ changes relatively little. A density dependent scaling
(by a factor $b$) leads to a variation in $\gamma$ for similar values of
$T_{0}$. cite evolves the thermal and ionization state of the IGM by using the
outputs of a gadget-3 simulation. A self-consistent gadget-3 simulation has
been performed with the KS19 UVB (the curves with $a=1$ and $b=0$.). For
generating the variation in thermal parameters, we assume ionization
equilibrium (but no thermal equilibrium) in cite. We generate 999 different
thermal parameter evolution histories.

In this section, we explain our method of generating different thermal
histories for a given UVB model. We vary the thermal parameter evolution by
modifying the H i, He i and He ii photo-heating rates of the KS19 UVB.
Following the approach of Becker et al. (2011), we scale the photo-heating
rates ($\epsilon_{\rm i}$) by two free parameters $a,b$ such that
$\epsilon_{\rm i}=a\>\Delta^{b}\>\epsilon_{\rm i}^{\rm KS19}$ where ${\rm
i}\equiv$ [ H i, He i, He ii ]. A change in factor $a$ ($b$) leads to
variation in $T_{0}$ ($\gamma$) while keeping the evolution of $\gamma$
($T_{0}$) relatively unchanged. It is important to note that the variation in
thermal parameters obtained by scaling the photo-heating rate is not physical.
Ideally, one would like to vary the properties of the ionizing sources such
as, luminosity function, QSO SED index, redshift of reionization etc., and
generate the UVB model from a cosmological 1D radiative transfer code.
However, such an approach cannot produce a large variation in $T_{0}$ and
$\gamma$ (but see Oñorbe et al., 2019). Our approach, even though not
physical, is useful in practice to probe a large range of the parameters
$T_{0}$ and $\gamma$ on a finely sampled grid. Fig. 18 shows examples of the
evolution in $T_{0}$ and $\gamma$ obtained by varying $a$ and $b$ for several
thermal histories. The case $a=1$ and $b=0$ shows the $T_{0}$ and $\gamma$
evolution obtained using the default photo-heating rates from the KS19 UVB.
Note that the shock heated gas ($T>10^{5}$ K) is mostly unaffected by the
photo-heating rate scaling, while the low density gas responds to this scaling
differently and changes normalization ($T_{0}$) and slope ($\gamma$) of the
TDR.

When measuring the thermal parameters, it is important to sample the
$T_{0}-\gamma$ plane densely enough. A coarse sampling of $T_{0}-\gamma$ plane
can lead to incorrect best fit values and excessive smoothing of the
$\chi^{2}$ field which, in turn, can result in underestimation of the
uncertainties of the thermal parameters. On the other hand a very dense
sampling of $T_{0}-\gamma$ would be computationally very expensive. We found
it to be the best compromise to choose the $T_{0}-\gamma$ grids such that
$T_{0}$ is varied from $\sim 6000$ K to $\sim 24000$ K in steps of 500 K while
$\gamma$ is varied from $\sim 0.7$ to $\sim 2.0$ in steps of 0.05 at
$z=3$555The corresponding $T_{0}$ and $\gamma$ at other redshifts differ by
$\sim 7$ percent. We have thus simulated the Ly$\alpha$ forest for $37\times
27=999$ different thermal histories. The computational time required to run
the thermal histories and extract mock spectra for 999 UVB models is $\sim
2.25$ million cpu hours.

##  Appendix E The sensitivity of Ly$\alpha$ forest statistics to $T_{0}$,
$\gamma$

###  E.1 Flux power spectrum (FPS)

Figure 19: The figure shows the sensitivity of the dimensionless FPS to
$T_{0}$ and $\gamma$ for three models (i) $T_{0}=10000$ K, $\gamma=1.50$ (blue
solid line), (ii) $T_{0}=20000$ K, $\gamma=1.50$ (green dashed line) and (iii)
$T_{0}=10000$ K, $\gamma=1.10$ (red dotted line) at $2.9\leq z\leq 3.1$. The
FPS is systematically lower for models with higher $T_{0}$ at $0.02\lesssim k$
( km-1 s) $\lesssim 0.2$, while the FPS is systematically higher for models
with smaller $\gamma$ at $0.02\lesssim k$ ( km-1 s) $\lesssim 0.2$. The FPS at
scales $k\lesssim 0.02$ is similar in all models because the underlying dark
matter density field is the same in all the models, while the FPS at $0.2\geq
k$ is dominated by observational systematics such as noise and instrumental
broadening.

The FPS is sensitive to a wide range of parameters: thermal parameters,
$\Gamma_{\rm HI}$, cosmological parameters ($\Omega_{b},\Omega_{m},n_{s}$),
the free-streaming of dark matter particles and neutrinos etc (Viel et al.,
2004b; Garzilli et al., 2017; Iršič et al., 2017a, b; Gaikwad et al., 2017a;
Khaire et al., 2019; Gaikwad et al., 2020). For a given set of cosmological
parameters in the $\Lambda$CDM framework, the FPS is sensitive to small scale
smoothing of the Ly$\alpha$ flux arising from heating of the IGM. Fig. 4
(panel A) illustrates the effect of increased temperature on the Ly$\alpha$
flux. The flux in the higher-$T_{0}$ model (for brevity, let us call it the
hot model) is somewhat smoother and shows less peaked absorption as compared
to the lower-$T_{0}$ (the cold model) model.

To calculate the FPS, we take the Fourier transform $\mathscr{F}$ of the flux
contrast $\delta_{f}=(F/\langle F\rangle)-1$ where $\langle F\rangle$ is the
mean flux in a given redshift bin. We bin the Fourier power
$|\mathscr{F}|^{2}$ in logarithmic bins centered on $k$ ( km-1 s)
$=0.005,\cdots,10.0$ and with bin width $\delta\log k=0.125$. This probes
scales from $\sim 1250$ km s-1 to $\sim 0.628$ km s-1. The lower and upper
limits of $k$ are motivated by the size of our simulation box and the
resolution of the observed spectra, respectively666While computing the FPS, we
do not concatenate spectra as there is no correlation beyond the length of the
simulation box..

Fig. 19 shows the sensitivity of the FPS to the thermal parameters $T_{0}$ and
$\gamma$. The power at $0.01\leq k\>{\rm(s\>km^{-1})}\leq 0.1$ for a model
with higher $T_{0}$ (smaller $\gamma$) is smaller due to more smoothing of the
flux. The FPS at scales $k<0.01$ km-1 s for all the models is the same because
the dark matter density field is the same for all the models, while the FPS at
$k>0.3$ km-1 s is dominated by observational systematic effects such as
instrumental broadening and finite S/N per pixel. The FPS at $0.02\leq
k\>{\rm(s\>km^{-1})}\leq 0.3$ ( $20\leq v\>{\rm(km\>s^{-1})}\leq 300$) is
clearly sensitive to thermal parameters.

###  E.2 Wavelet Statistics (WS)

Figure 20: The left, middle and right panels show the sensitivity of the
wavelet PDF to $T_{0}$ and $\gamma$ for a wavelet scale of $s_{n}=30,\>50$ and
$100$ km s-1, respectively (for $2.9\leq z\leq 3.1$). Each panel shows the
wavelet PDF for three models (i) $T_{0}=10000$ K, $\gamma=1.50$ (blue solid
line), (ii) $T_{0}=20000$ K, $\gamma=1.50$ (green dashed line) and (iii)
$T_{0}=10000$ K, $\gamma=1.10$ (red dotted line). The wavelet PDF is
systematically shifted to higher $\log{\rm A_{Sn}}$ for models with smaller
$T_{0}$, while the wavelet PDF is systematically shifted to smaller $\log{\rm
A_{Sn}}$ for model with larger $\gamma$. The effect of $T_{0}$ and $\gamma$
variation is smaller at small wavelet scales, i.e., $s_{n}=30$ km s-1 (panel
A). The effect of $T_{0}$ on the wavelet PDF is larger at intermediate wavelet
scale $s_{\rm n}=50$ km s-1, whereas the effect of $\gamma$ is larger at
higher wavelet scales $s_{\rm n}=100$ km s-1. We show in Fig. 6 that the
simultaneous constraints using wavelet PDFs at scales $30\leq s_{\rm n}$ ( km
s-1) $\leq 100$ yield better constraints on $T_{0}$ and $\gamma$. All the
wavelet PDFs are calculated using KDE (see appendix E.6).

Wavelets are wave-like oscillations that are localized in both the real space
and frequency domains. One can use a suitable wavelet to extract the scales in
transmitted Ly$\alpha$ flux that are sensitive to $T_{0}$ and $\gamma$
(Zaldarriaga, 2002; Theuns et al., 2002; Lidz et al., 2010; Garzilli et al.,
2012). Mathematically, this is equivalent to a convolution of the transmitted
Ly$\alpha$ flux with the wavelets. Throughout this work, we use a Morlet
wavelet that looks like a Gaussian in Fourier space777We have also performed
analysis with Mexican hat and Shannon wavelets and found that the results are
insensitive to the choice of wavelet.. The oscillation frequency of the
wavelet depends on the center of the Gaussian in Fourier space. While the
width of the wavelet is related to the width of the Gaussian in Fourier space.
The Morlet wavelet has the functional form,

$\displaystyle\Psi(v,s_{\rm n})$ $\displaystyle=C\>\exp\bigg{[}-\frac{2\pi
vi}{s_{\rm n}}\bigg{]}\>\exp\bigg{[}\frac{-v^{2}}{2s^{2}_{\rm n}}\bigg{]}$ (2)

where $v$ is velocity, $s_{\rm n}$ is wavelet scale, $C$ is a normalization
constant obtained by demanding that wavelets are square integrable i.e.,
$\int\>|\Psi(v,s_{\rm n})|^{2}\>dv=1$. The wavelet transform of the flux
(hereafter the wavelet field $A_{\rm Sn}$) is the convolution of the wavelet
$\Psi(v,s_{\rm n})$ with the transmitted Ly$\alpha$ flux ($F$).

$\displaystyle A_{S{\rm n}}$ $\displaystyle=\int\limits_{v_{\rm min}}^{v_{\rm
max}}\>dv^{\prime}\>F(v^{\prime})\>\Psi(v-v^{\prime},s_{\rm n}).$ (3)

Fig. 4 (panel B) demonstrates the sensitivity of the wavelet amplitude (for a
wavelet scale $s_{\rm n}=50$ km s-1) to the temperature of the IGM888Wavelet
amplitudes can be negative hence we plot logarithm of absolute values of
wavelet amplitude.. The wavelet amplitudes ($\log A_{\rm 50}$) in the hot
model are systematically smaller than those from the cold model. It is also
evident that the wavelet amplitude picks out small scale variations in the
flux field. To quantify the effect of $T_{0}$ and $\gamma$ on the wavelet
amplitude, we focus on the probability distribution function of the wavelet
field. Usually, the wavelet field needs to be smoothed (e.g., using a boxcar
filter of $\sim 1000$ km s-1, see Lidz et al., 2010) to reduce the effect of
flux noise on the wavelet PDF. However, we find that such a boxcar smoothing
reduces the sensitivity of the wavelet PDF to $T_{0}$ and $\gamma$. Instead we
use kernel density estimation to calculate the wavelet PDF. This reduces the
effect of noise while preserving the sensitivity of the wavelet PDF to $T_{0}$
and $\gamma$ (see Appendix E.6 for details).

Fig. 20 shows the comparison of the wavelet PDF for 3 models (i) $T_{0}=10000$
K, $\gamma=1.5$, (ii) $T_{0}=20000$ K, $\gamma=1.5$, and (iii) $T_{0}=10000$
K, $\gamma=1.1$. The wavelet PDF is systematically shifting to lower values
for higher $T_{0}$ at all the three wavelet scales $s_{\rm n}=30,50,100$ km
s-1. Fig. 20 also shows that the wavelet scale $s_{\rm n}=50$ km s-1 is most
sensitive to variation in $T_{0}$. At smaller wavelet scales $s_{\rm n}=30$ km
s-1 it is less sensitive to $T_{0}$ due to observational effects (noise and
resolution) while larger wavelet scales $s_{\rm n}=100$ km s-1 trace the
underlying dark matter density field which is the same for the three models.
The wavelet amplitudes are systematically lower for higher $\gamma$ values.
This is expected because for the same normalization $T_{0}$, a higher $\gamma$
value corresponds to larger temperature at $\Delta>1$. Since a significant
fraction of the Ly$\alpha$ absorption at $z\sim 3$ comes from $\Delta>1$, the
Ly$\alpha$ absorption lines in models with higher $\gamma$ are expected to be
broader. But it is interesting to note that the wavelet scale $s_{\rm n}=100$
km s-1 is slightly more sensitive to $\gamma$ than to $T_{0}$. Thus, if we
simultaneously use wavelet PDFs at scales $30\leq s_{\rm n}$ ( km s-1) $\leq
100$, the constrains on $T_{0}$ and $\gamma$ should be better than using the
wavelet PDF for the individual scales.

###  E.3 Curvature Statistics

Figure 21: The figure shows the sensitivity of the curvature PDF to $T_{0}$
and $\gamma$ at $2.9\leq z\leq 3.1$. The curvature PDF is systematically
shifted to higher (smaller) $\log\langle|\kappa|\rangle$ for smaller $T_{0}$
(larger $\gamma$). The curvature PDFs are calculated using KDE (see appendix
E.6). The curvature PDFs show a similar behaviour as the wavelet PDF except
that it is less sensitive to thermal parameters than the wavelet PDF (at
$s_{\rm n}=50$ km s-1).

The thermal broadening of the absorption lines can also be characterized by
the curvature of the absorption features (Becker et al., 2011). The curvature
($\kappa$) is essentially a measure of rate of change of direction of a point
that moves on a curve and is defined as

$\displaystyle\kappa$
$\displaystyle\equiv\frac{F^{\prime\prime}}{[1+(F^{\prime})^{2}]^{3/2}}\;\;,$
(4)

where $F^{\prime}=dF/dv$ and $F^{\prime\prime}=d^{2}F/dv^{2}$ are the first
and second derivatives of the flux field (with respect to velocity),
respectively. Similar to wavelet amplitudes, curvature can also be affected by
the finite S/N of the spectra. To circumvent this difficulty, Becker et al.
(2011) fitted a spline curve to the transmitted Ly$\alpha$ flux field and used
that to compute the curvature. We, however, use a simpler approach of
smoothing the flux field by a Gaussian filter of FWHM $\sim 10$ km s-1 (for a
similar approach, see Padmanabhan et al., 2015). Fig. 4 (panel C) shows the
comparison of the curvature field for the hot and cold models. Similar to the
wavelet field, the curvature amplitude is higher for the cold model as
compared to the hot model. In this work, we use the curvature PDF to measure
$T_{0}$ and $\gamma$. Fig. 21 shows the sensitivity of the curvature PDF to
$T_{0}$ and $\gamma$. The curvature PDF is shifted to lower values of
$\log\langle|\kappa|\rangle$ for hot models as compared to cold models. The
curvature PDF is also sensitive to variations in $\gamma$. However, note that
the curvature PDF shows less sensitivity to $T_{0}$ and $\gamma$ than the
wavelet PDF.

###  E.4 Doppler($b$)-parameter probability distribution function (BPDF)

Figure 22: Left, middle and right panels show the $b-\log{\rm N_{HI}}$
distribution for models (i) $T_{0}=10000$ K, $\gamma=1.50$, (ii) $T_{0}=20000$
K, $\gamma=1.50$ and (iii) $T_{0}=10000$ K, $\gamma=1.10$, respectively. The
$b-\log{\rm N_{HI}}$ distribution shows a systematic increase in $b$ values
for the higher $T_{0}$ model while the slope of the correlation remains
approximately the same. For higher values of $\gamma$, the $b-\log{\rm
N_{HI}}$ distribution becomes steeper. The above distributions are calculated
from $\sim 1300$ simulated spectra at $2.9\leq z\leq 3.1$. Due to the limited