Probing the thermal state of the intergalactic medium at z > 5 with the transmission spikes in high-resolution Ly α forest spectra

We compare a sample of ﬁve high-resolution, high S/N Ly α forest spectra of bright 6 < z < ∼ 6 . 5 QSOs aimed at spectrally resolving the last remaining transmission spikes at z > 5 with those obtained from mock absorption spectra from the Sherwood and Sherwood-Relics suites of hydrodynamical simulations of the intergalactic medium (IGM). We use a proﬁle ﬁtting procedure for the inverted transmitted ﬂux, 1 − F , similar to the widely used Voigt proﬁle ﬁtting of the transmitted ﬂux F at lower redshifts, to characterise the transmission spikes that probe predominately underdense regions of the IGM. We are able to reproduce the width and height distributions of the transmission spikes, both with optically thin simulations of the post-reionization Universe using a homogeneous UV background and full radiative transfer simulations of a late reionization model. We ﬁnd that the width of the ﬁtted components of the simulated transmission spikes is very sensitive to the instantaneous temperature of the reionized IGM. The internal structures of the spikes are more prominent in low temperature models of the IGM. The width distribution of the observed transmission spikes, which require high spectral resolution ( ≤ 8 km s − 1 ) to be resolved, is reproduced for optically thin simulations with a temperature at mean density of T 0 = (11000 ± 1600 , 10500 ± 2100 , 12000 ± 2200) K at z = (5 . 4 , 5 . 6 , 5 . 8). This is weakly dependent on the slope of the temperature-density relation, which is favored to be moderately steeper than isothermal. In the inhomogeneous, late reionization, full radiative transfer simulations where islands of neutral hydrogen persist to z ∼ 5 . 3, the width distribution of the observed transmission spikes is consistent with the range of T 0 caused by spatial ﬂuctuations in the temperature-density relation. ∼ 0 . 33) than those in the optically thin simulations (∆ ∼ 0 . 25). This is due to the spatial ﬂuctuations in the photo-ionization rate and temperature that are present in the radiative transfer simulations. We discuss the ∆ τ vs log (cid:101) N HI correlation for the optically thin and radiative transfer simulations in § 4.2.1 and § 6.2.1, respectively.


INTRODUCTION
The reionization of intergalactic H i and He i by ultraviolet photons from stars and black holes in the first galaxies is one of the major phase transitions of the universe (Fan Downloaded from https://academic.oup.com/mnras/article-abstract/doi/10.1093/mnras/staa907/5818338 by University of Nottingham user on 17 April 202017 April et al. 200117 April , 2006Robertson et al. 2010;Bolton et al. 2011;Planck Collaboration et al. 2014. Photo-heating during the reionization increases the temperature of the intergalactic medium (IGM) (Hui & Gnedin 1997;Trac et al. 2008;McQuinn et al. 2009;McQuinn & Upton Sanderbeck 2016). The thermal and ionization histories of the Universe are thus interlinked, and can be used in tandem to understand the process of reionization and the nature of ionizing sources (Haehnelt & Steinmetz 1998;Furlanetto & Oh 2009;McQuinn et al. 2011;Becker et al. 2015;Chardin et al. 2017;Davies et al. 2018;Keating et al. 2018;Kulkarni et al. 2019b).
The H i Lyα forest observed along sightlines towards bright QSOs is frequently used to constrain the thermal and ionization state of the IGM at z < 5. For underdense or moderately overdense regions, the thermal state of the IGM is often approximated as a power law, parameterized as the temperature at mean density and the slope of the power law (T0 and γ Lidz et al. 2010;Becker et al. 2011;Rudie et al. 2012;Garzilli et al. 2012;Bolton et al. 2014;Boera et al. 2014;Hiss et al. 2018;Telikova et al. 2018Telikova et al. , 2019Walther et al. 2019;Boera et al. 2019). Other relevant parameters are the H i photo-ionization rate (ΓHI Rauch et al. 1997;Bolton & Haehnelt 2007;Faucher-Giguère et al. 2008;Calverley et al. 2011;Becker & Bolton 2013;Gaikwad et al. 2017a,b;Viel et al. 2017;Khaire et al. 2019) and the pressure smoothing scale (Gnedin & Hui 1998;Peeples et al. 2010;Kulkarni et al. 2015;Lukić et al. 2015;Rorai et al. 2017a). Due to the rapid increase of the Lyα opacity with redshift, the transmitted Lyα flux at z > 5 is close to zero, with the occurrence of a few transmission spikes indicating that the reionization process is inhomogeneous Bosman et al. 2018;Eilers et al. 2018). Analysis of the transmission spikes towards ULAS J1120+0641 QSO (zem = 7.084, Becker et al. 2015; Barnett et al. 2017) with numerical simulations suggests that these spikes correspond to underdense, highly ionized regions of gas (Gnedin et al. 2017;Chardin et al. 2018;Kakiichi et al. 2018;Garaldi et al. 2019;Nasir & D'Aloisio 2019). These studies find that the number and height of the spikes are sensitive to the ionization fraction ( xHII) of the IGM, which in turn depends on the photo-ionization rate and the temperature of the gas in the ionized regions. Perhaps somewhat surprisingly, however, Garaldi et al. (2019) found that the spike shape (especially the widths of the spikes) appeared to be only weakly correlated with the temperature of the IGM. We revisit this question here with a larger sample of higher resolution, higher quality Lyα forest spectra which we compare to high-resolution hydrodynamical simulations of the IGM using the Sherwood and Sherwood-Relics simulation suites. These incorporate simulations with a homogeneous UV background as well as full radiative transfer simulations of inhomogeneous reionization Kulkarni et al. 2019b, Puchwein et.al. 2020. Note that in the main analysis we used the Sherwood-Relics simulation suite, but complemented these with simulations from the Sherwood simulation suite for additional tests performed in the Appendices.
There are several difficulties in probing the effect of the thermal state of the IGM on Lyα transmission spikes. High-redshift QSOs, while intrinsically luminous, are nevertheless faint and observed spectra are normally taken at moderate resolution as e.g. offered by VLT/X-Shooter, i.e. with 35 km s −1 or worse. As a result most spectra of highredshift transmission spikes are of rather modest quality and the number of well observed transmission spikes is necessarily still small due to the rarity and faintness of high-z background QSOs (Kulkarni et al. 2019a). We improve this situation here and present a sample of 5 high resolution (FWHM ∼ 6 km s −1 ) and high S/N (∼ 10) QSO absorption spectra obtained using the Magellan Inamori Kyocera Echelle (MIKE) spectrograph on the Magellan II telescope (Bernstein et al. 2003), and the High-Resolution Echelle Spectrograph (HIRES) on the Keck I telescope (Vogt et al. 1994).
Similarly, simulated transmission spikes may be drawn from simulations with moderate mass or spatial resolution that is only sufficient to produce mock spectra mimicking moderate resolution spectrographs like X-Shooter (35 km s −1 , but see Garaldi et al. 2019). The thermal smoothing scale of the Lyα transmitted flux due to the IGM temperature (T ∼ 10 4 K, b ∼ 13 km s −1 ) is, however, significantly smaller than this. As discussed in detail by Bolton & Becker (2009), rather high mass or spatial resolution is required to resolve the small scale structure in the underdense regions of the IGM probed by Lyα forest spectra at high redshift. Previous theoretical work has focused on analyzing the spikes in radiative transfer simulations (Gnedin et al. 2017;Chardin et al. 2018;Garaldi et al. 2019). While these simulations are physically well motivated, they are computationally expensive and it is hard to disentangle the effect of the thermal state of the IGM on the transmission spikes from the reionization history and numerical limitations.
We have thus chosen to analyze the transmission spikes first in very high resolution, high dynamic range optically thin simulations with different thermal and reionization histories where a single parameter is varied at a time while keeping other parameters fixed. We therefore use simulations from the Sherwood and Sherwood-Relics simulation suite , Puchwein et.al. 2020 to show how spike properties depend on the IGM thermal state, the H i photo-ionization rate ΓHI, and (in the appendices) the pressure smoothing scale, as well as the mass resolution and box size of the simulations. Once we have established this we investigate the effect of inhomogeneous reionization in more physically motivated, spatially inhomogeneous H i reionization simulations including radiative transfer (see Kulkarni et al. 2019b;Keating et al. 2019).
Another problem when comparing simulated and observed Lyα forest spectra is the accurate characterisation of the transmission spike properties. Transmission spikes are often asymmetric and transmission features appear "blended". The often used simple definition of height and width of spikes based on the maximum and FWHM of the transmitted flux will thus not capture the detailed information contained in the complex spike shapes, and could be the reason that the analyses performed so far show little or no correlations with astrophysical parameters (see e.g. the width vs temperature correlation in Garaldi et al. 2019). Characterizing the shape of transmission spikes becomes even more crucial for high S/N, high resolution QSO absorption spectra. In practice the problem is very similar to that of the characterisation of Lyα absorption lines at lower redshift. To utilise the practical experience gained Downloaded from https://academic.oup.com/mnras/article-abstract/doi/10.1093/mnras/staa907/5818338 by University of Nottingham user on 17 April 2020 in this area with existing software packages (e.g. Gaikwad et al. 2017b) we characterize the transmission spikes by fitting Voigt profiles to the "inverted" transmitted flux, 1 − F . We will later show that the fitted parameters obtained in this way are well correlated with physical properties (e.g. the density and temperature) of the gas associated with the transmission spikes. We will also study how the statistics of the fitted parameters depend on the astrophysical parameters ΓHI, T0 and γ. Unlike for absorption lines, there is no direct physical motivation for fitting Voigt profiles to 1-F, so this should be considered as a purely heuristic approach to comparing simulated and observed spectra. However, as we will see this does not mean that the fit parameters obtained do not correlate with physical properties.
The main goal of this paper is to constrain the thermal state of the IGM at 5.3 ≤ z ≤ 5.9. The paper is organized as follows: In §2 we present the high resolution spectra of 5 z > 6 QSOs and discuss qualitatively the physical origin of Lyα transmission spikes. We present the properties of transmission spikes in optically thin simulations in §3 and §4. We demonstrate the sensitivity of spike statistics to the IGM thermal state in §5. The main results of the paper are presented in §5.3 and §6 by comparing the observed spike statistics with those from optically thin and radiative transfer simulations. We summarize our findings in §7. We assume a flat ΛCDM cosmological model (ΩΛ, Ωm, Ω b , σ8, ns, h, Y ) ≡ (0.692, 0.308, 0.0482, 0.829, 0.961, 0.678, 0.24) consistent with Planck Collaboration et al. (2014,2018). All distances are given in comoving units unless specified. ΓHI expressed in units of 10 −12 s −1 is denoted by Γ12.

TRANSMISSION SPIKES IN
HIGH-RESOLUTION, HIGH-REDSHIFT Lyα FOREST SPECTRA 2.1 Observations: high-resolution spectra of transmission spikes The data consist of the high resolution echelle spectra of five recently discovered z > 6 QSOs. The objects were chosen for their brightness, and individual targets were further selected to maximize the exposure time during a given observing run. Table 1 gives the name of each object, the emission redshift zem, the J AB band magnitude (from the compilation of Ross & Cross 2019), the total on-source exposure time T in hours, and a typical signal-to-noise ratio per pixel. The objects were observed under mostly photometric conditions in sub-arcsec seeing. The first four objects were observed with the MIKE instrument (Bernstein et al. 2003) on the Magellan II telescope at Las Campanas Observatory. A 0.5" wide slit gave a measured spectral resolution of 5 km s −1 (FWHM). The spectra were binned onto 2 km s −1 wide pixels. The spectrum of SDSSJ010013.02+280225.8 was obtained with the HIRES instrument (Vogt et al. 1994) on the Keck I telescope, and a 0.861" wide slit, giving a resolution of 6.1 km s −1 (FWHM), sampled by 2.5 km s −1 wide bins. The data were reduced with a custom pipeline .
Optimal sky-subtraction on the individual, un-rectified exposures was performed according to the prescription by Kelson (2003).
For the continuum model, a power law was assumed. As the high resolution echelle spectra only had limited coverage of the unabsorbed region redward of the Lyα emission, the continuum was derived from flux-calibrated, low resolution spectra of the same QSOs which extended further to the red. For objects 1-3 (in table 1), the continua were determined from the discovery spectra, using fits to regions redward of Lyα avoiding broad emission lines, whereas the power law slopes for objects 4 and 5 were taken from the literature cited. The continuum was then scaled uniformly to match the flux-calibrated high resolution spectrum in the overlap region with the low-resolution spectrum redward of Lyα, and divided into the spectrum. To correct for the rapidly variable region of the spectrum near the Lyα emission line, the emission line region was fitted with a higher order polynomial in the previously continuum divided spectrum, which then was multiplied into the previous continuum fit. The final continuum thus obtained was divided into the data. As we show later, the width of the fitted components of the transmission spikes are relatively robust to continuum fitting uncertainty.
During data reduction a certain degree of smoothing is introduced into the data that appears as a discrepancy between the RMS fluctuations in the final data and the propagated error array, and generally results in an underestimate of the reduced χ 2 ν when fitting line profiles. To counteract this problem, a correction factor (generally a number close to unity varying slowly with wavelength) was derived from the observed ratio between the RMS fluctuations and the error array, as determined from wavelength windows in the spectrum with zero flux. A linear fit to the correction factor as a function of wavelength was divided into the error array to obtain a reduced χ 2 ν ∼ 1 during profile fitting.

Characterising width and height of individual components
Fig. 1 compares transmission spikes in a high-resolution Lyα forest observation of the QSO J043947.08+163415.7 with the MIKE spectrograph with simulated spikes drawn from the Sherwood-Relics simulation suite (see Section 3.1, Puchwein et.al. 2020 in prep), for cold and hot models with a spatially uniform UV background. The transmission spikes have complex shapes composed of many asymmetric and blended features. Even isolated transmission spikes are often highly asymmetric and consist of two or more "components". In order to facilitate a more quantitative discussion of the transmission spikes, we focus on two properties, namely the height and width of individually identifiable components. We quantify these (see §4) by fitting the "inverted" transmitted flux, 1 − F , with multi-component Voigt profiles, similar to the fitting of Voigt profiles of the transmitted flux, F , often employed at lower redshift to characterise absorption lines. The best fits to observed and simulated spectra are shown in Fig. 1 by the red curve. The location of individual identified components are marked by black vertical lines. Fig. 1 show simulated spectra for the hot and cold model and illustrate the sensitivity of spike features to the thermal state of the IGM. The main effect of increase in IGM temperature is to reduce the internal structure of the spikes (i.e., more blending). As a result, fewer and broader components are required to fit the transmission spikes in the hot model. In general, the spikes in the hot model are (i) Downloaded from https://academic.oup.com/mnras/article-abstract/doi/10.1093/mnras/staa907/5818338 by University of Nottingham user on 17 April 2020  Figure 1. Examples of transmission spikes from observed spectra (panel A1) and simulated spectra from cold (panel B1) and hot (panel C1) optically thin simulations drawn from the Sherwood-Relics simulation suite at z ∼ 5.59. The widths and heights of the spikes are sensitive to the thermal and ionization state of the IGM. The shape and number of spikes in the observed spectra are similar to the simulated spectra from the hot model. The resolution and noise properties of the simulated spectra are chosen to match the observed spectra. As described in the main text we fit the "inverted" transmitted flux, 1 − F , with multi-component Voigt profiles using the software package viper described in detail in Gaikwad et al. (2017b). In viper, the number of components to be fitted in given region is decided automatically by minimizing the Akaike information criteria with correction (see Gaikwad et al. 2017b, for details). The top panels show the input spectra (F , blue solid curve) and fitted spectra (F ,in red solid curve). The black solid lines mark the location of the centres of Voigt components identified and fitted by viper. The bottom panels show that the residuals between the input and fitted spectra are random and less than 11 per cent. The number of components identified in the cold model is larger than in the hot model due to the smoother transmitted flux distribution.
broader, (ii) larger in height, (iii) more asymmetric and (iv) more blended (hence fewer in number) than those in the corresponding cold model. It is interesting to note that the number and shape of the spikes in the hot model are qualitatively very similar to those in the observed spectra. We will analyse this more quantitatively below.

The origins of transmission spikes
The complex shapes of the transmission spikes shown in the last section will be a superposition of features in the line-ofsight distribution of density, temperature, photo-ionization rate and peculiar velocity. To get a better feel for this in Fig. 2 we show mock spectra where we isolate the effect of varying these parameters, one at a time. Fig. 2 illustrates how the variation in any of these physical quantities along a sightline can lead to regions of smaller H i optical depth and hence transmission spikes. Panels A1-A6, B1-B6 and C1-C6 show the effect of underdensity, enhanced ΓHI and enhanced temperature along a sightline, respectively. All these effects can result in a smaller number density of neutral hydrogen, nHI, along the sightline. Panels D1-D6 show that a diverging peculiar velocity field along the sightline can also produce transmission spikes. Hot, ionized underdense regions subjected to high photo-ionization rates and diverging peculiar velocities will thus produce the most prominent transmission spikes in high−z QSO absorption spectra (Gnedin et al. 2017;Chardin et al. 2018;Garaldi et al. 2019). The transmission spikes shown in Fig. 2 are by construction isolated, symmetric and simple. As discussed in the previous section, transmission spikes in both observed spectra and spectra drawn from cosmological hydrodynamical simulations have complicated shapes, as generally more than one of the above effects contribute.

The Sherwood and Sherwood-Relics simulation suites
We use cosmological hydrodynamical simulations from the Sherwood-Relics simulation suite to investigate transmission spikes in the high-redshift Lyα forest; see Table 2 for an overview. The simulations were performed with a modified version of the p-gadget-3 code (itself an updated In each column, a different physical parameter is varied. In the first three columns (i.e. for varying ∆, Γ 12 and T ) the occurrence of spikes is due to a change in n HI . Panels D1-D6 instead shows the formation of a transmission spike due to a diverging velocity flow along the sightline while n HI remains constant. Realistic transmission spikes (see Fig. 4), have more complicated shapes and will be due to a combination of these effects.
version of the gadget-2 code presented in Springel 2005). The code uses a tree-particle mesh gravity solver for following cosmic structure formation and a manifestly energy and entropy-conserving smoothed particle hydrodynamics scheme (Springel & Hernquist 2002) for following the hydrodynamics. The Sherwood-Relics simulations build upon the original Sherwood simulation suite (which is used in Appendix ?? to study numerical convergence) in that the initial conditions were generated in the same way, and much of the modeling of the IGM and Lyα forest is based on similar methods ) 1 . Our main production runs follow 2 × 2048 3 particles in a (40 h −1 cMpc) 3 volume, corresponding to a gas mass resolution of 9.97 × 10 4 h −1 M . For the gravitational softening we adopt 0.78 h −1 ckpc. Star formation is treated with a rather simplistic but numerically efficient scheme in which all gas particles with densities larger than 1000 times the mean cosmic baryon density and temperatures smaller than 10 5 K are converted to collisionless star particles. While this does not produce realistic galaxies, it allows robust predic-1 https://www.nottingham.ac.uk/astronomy/sherwood/ tions of the properties of the IGM (Viel et al. 2004a). Photoheating and photo-ionization are followed based on external UV background models. In a departure from the Sherwood simulations, we use a non-equilibrium ionization and cooling/heating solver (Puchwein et al. 2015;Gaikwad et al. 2019) for following the thermochemistry of hydrogen and helium. This ensures that no artificial delay between the reionization of gas and its photo-heating is present. We have also replaced the slightly modified Haardt & Madau (2012) UV background used in Sherwood with the fiducial UV background model from Puchwein et al. (2019) in our default run. This results in a more realistic reionization history with hydrogen reionization finishing at z ≈ 6.2. The cold /hot models were obtained by decreasing/increasing the H i and He i photo-heating rates in the fiducial UV background model by a factor of 2 and the He ii photo-heating rate by a factor of 1.7, while keeping all photo-ionization rates fixed. Mock Lyα forest spectra were extracted from all simulations as described in Bolton et al. (2017) (see also Gaikwad et al. 2018Gaikwad et al. , 2019. Through out this work, we use simulations from the Sherwood-Relics simulation suite for the main analysis. Ad-ditional simulations from the Sherwood simulation suite are used for convergence tests as described in Appendix ??. 3.2 The temperature density relation (TDR) in optically thin simulations Fig. 3 shows the TDR in the hot and cold models of the Sherwood-Relics simulation suite. We shall compare the hot and cold models to study the effect of temperature on transmission spikes. Note that the hot model is not only hotter than the cold and the default model but also has a flatter temperature density relation, and that the ionization state of the gas is also different in the models. This is because the recombination coefficient is temperature dependent. For the same photo-ionization rate the H i fraction is therefore smaller in the hot model. The models aton and patchy incorporate the effect of inhomogeneous UV background that we describe in 6.

Examples of transmission spikes in the hot and cold optically thin simulations
Fig . 4 shows the relevant physical properties along the two sightlines in the hot and cold models. As expected, the hot model shows (i) a smoother density (∆) field in real space, (ii) a smaller H i fraction ( xHI), (iii) a larger temperature (T ) and (iv) a smoother velocity (v) field compared to the corresponding cold model. The smoothing of the real-space density field and the velocity field can be attributed to the increased effect of pressure smoothing in the hot model. The resultant Lyα optical depth and transmitted flux calculated from the ∆, xHI, T and v fields are shown in panels A6-B6 and panels A7-B7, respectively. Note that the location of spikes in the green and yellow shaded regions in redshift space differs in real space due to the effect of peculiar velocities. The transmission spikes in the hot and cold models have complicated shapes, qualitatively similar to that in the observed spectra ( Fig. 1). The smoother transmission features in the simulated spectra of the hot model are more similar to those in the observed spectra than those in the more "spiky" spectra in the cold model. In both models, the transmission spikes correspond to regions of low H i optical depth ( τHI ≤ 4, black dashed line). The peaks in transmission are well correlated with those in the optical depth weighted overdensities ( ∆τ < 0.8). 2 It is clear from Fig. 4 that the spikes in the hot model are smoother, broader and more prominent than in the cold model. Furthermore the number of individual transmission components is significantly smaller in the hot model than in the cold model due to the "thermal blending" of transmission features. The shape and number of transmission spikes in our simulated spectra are clearly sensitive to the thermal and ionization state of the IGM in a manner that we will quantitatively discuss below. Table 2 shows the physical effects responsible for the occurrence of transmission spikes in the optically thin hot and cold simulations. Most of the spikes (∼ 89 percent) in the hot and cold models occur in underdense regions. Around 15 percent of spikes show a diverging velocity field along the sightline. The effect of enhanced temperature on the occurrence of transmission spikes is marginal in both hot and cold optically thin simulations.
In summary, Fig. 4 and Table 2 show that underdense ( ∆τ < 0.8), more highly ionized (minimum in xHI) and hotter regions along a sightline produce more prominent spikes. This motivates us to quantify the shape of spikes and to introduce statistics that are sensitive to the thermal and ionization parameters of the IGM. Most of the absorption features in the high-z (z ∼ 6) Lyα forest are saturated (F ∼ 0) and strongly blended. It is well known that Voigt profile decompositions are highly degenerate for saturated lines and very sensitive to systematic errors due to continuum fitting and treatment of noise properties (Webb & Carswell 1991;Fernández-Soto et al. 1996). The inverted transmitted flux, 1−F , however, becomes similar in appearance to the absorption features in the Lyα forest at lower redshift, where Voigt profile fitting is much less problematic. For convenience, we have thus fitted Voigt profiles to 1 − F , building on existing experience with Voigt profile decomposition of complex blended spectral profiles. At the redshift we consider here the transmission spikes are observed to be unsaturated, so effectively we use our VoIgt profile Parameter Estimation Routine (viper, Gaikwad et al. 2017b) to fit multi-component Voigt profiles to the transmission profiles. Similar to absorption lines, a simple, isolated and symmetric spike is fitted by 3 parameters: (i) a spike centre (λc), (ii) the logarithm of the pseudo-column density (denoted by log NHI) and (iii) a spike width (b). The pseudo-column density is thereby a measure of the deficiency of H i along the sightline where the spike occurs. For example, a larger value of log NHI means a large H i deficit hence a more prominent spike. Our measured log NHI are sensitive to the H i photoionization rate ΓHI. The interpretation of the other two parameters i.e., λc and width of the spikes remains unchanged when we fit 1 − F . As we show below, the distribution of spike widths is sensitive to the thermal state of the IGM and can be used to constrain the temperature of the IGM. 3 Our main aim is to use the distribution of spike widths to constrain the thermal state of the IGM. As we will show later, these are much less sensitive to IGM ionization state and continuum fitting uncertainties than the heights of the spikes. Note that we rescale the optical depth in different The aton simulation is post-processed with the radiative transfer code ATON, while the patchy simulation includes the effect of pressure/Jeans smoothing as well as the shock heating of the gas. For optically thin simulations the TDR in underdense and moderately overdense regions can be approximated as a power-law relation (T = T 0 ∆ γ−1 ) at ∆ ≤ 10. The best-fitting relation in the cold and hot models is shown by the black dashed lines. By construction, the temperature of the gas with ∆ ≤ 10 in the hot model is consistently larger than that in the cold model. As a result, the heights and widths of the components to be fitted to the transmission features are expected to be different in the hot and cold models. Unlike the optically thin simulations, the radiative transfer simulations (aton and patchy) do not exhibit a single power-law TDR at ∆ ≤ 10. For visual purposes, we show two power-laws that cover the range in temperature for the radiative transfer runs (panel C and D). The absence of gas with T > 30000 K in aton is because the shock heating is not captured self-consistently in the aton runs. We plot mass weighted (SPH particle) temperature and density in panel A, B and D, while we plot volume weighted temperature and density (calculated on grids) for the aton model (panel C). As a result, the number of gas elements (at ∆ ≤ 10) between the two straight power law TDR lines is larger in the aton simulations than in the patchy simulations. Note that gas with ∆ > 10 3 and T < 10 5 K has been converted into stars in all these simulations (Viel et al. 2004a). We discuss the TDR for optically thin and radiative transfer simulations in §3.2 and §6.2, respectively. Table 2. Origin of transmission spikes due to variation in physical parameters in our simulations at 5.5 < z < 5.7 (as shown in Fig. 2). a Fraction of all spikes with ∆τ ≤ 1. b Fraction of all spikes with Γ HI ≥ Γ HI,median where Γ HI,median is the optical depth weighted median Γ HI calculated from all the sightlines (i.e., regions with and without spikes). c Fraction of all spikes with Tτ ≥ T τ,median where T τ,median is the optical depth weighted median temperature calculated from all the sightlines (i.e., regions with and without spikes). d Fraction of all spikes with ∆v ≥ 6 km s −1 where ∆v is difference between mean velocity on redward and blueward side of spike center. For diverging velocity flow ∆v > 0, whereas for converging velocity flow ∆v < 0. The limit of ∆v = 6 km s −1 corresponds to the spectral resolution of the instrument. e L40N2048 DEFAULT L40N2048 COLD and L40N2048 HOT are optically thin simulations that do not include fluctuations in Γ HI . models to match observations as to account for the uncertainty in ΓHI at 5.3 < z < 5.9 (rescaling is not applied in Fig. 4). We show such rescaling does not significantly affect the widths of the line in Appendix ??. We now turn to connecting the observed spike shapes to optical depth weighted physical properties of the IGM in the optically thin simulations (see Eq. 12 in Gaikwad et al. 2017a). In Fig. 5, we show the optical depth weighted overdensity against spike height as measured by (pseudo) column density log NHI for the cold and hot models. In both mod- els the corresponding overdensity of spikes is ∆τ < 1 i.e., most of the spikes occur in underdense regions. Note that applying a rescaling of optical depth to correct for the uncertain amplitude of the UV background results in a systematic increase or decrease of log NHI. However, the overdensities corresponding to spikes are relatively robust. This is also evident from Fig. 4, where the spikes in the hot and cold models correspond to similar overdensities (rescaling is not applied in Fig. 4).

Fig
. 5 also illustrates that log NHI and ∆τ are anticorrelated i.e. a larger spike height and higher H i pseudo-column density corresponds to smaller overdensity ∆τ . We quantify the degree of anti-correlation by fitting a straight line of the form ∆τ = ∆0 [ NHI / 10 12.8 ] β . The normalisation (∆0 =0.24 for cold and 0.26 for hot) and slope (β = −0.23 for cold and −0.20 for hot) of the correlation in the two models are in good agreement with each other. This suggests that the thermal state of the gas does not have a strong effect on the ∆τ − log NHI correlation when optical depths are rescaled to match observations. Note, however, that the scatter in ∆τ for a given log NHI (e.g., at log NHI ∼ 12.6) log ( N HI /cm −2 ) . Panels A, B, C and D show the correlation of optical depth weighted overdensity ( ∆τ ) with the pseudo-column density ( log N HI ) of transmission spikes in, respectively, the cold, hot, aton and patchy models at 5.5 < z < 5.7. Irrespective of the model, the spikes correspond to underdensities i.e., ∆ ∼ 0.25 for log N HI = 12.8 in the optically thin simulations and ∆ ∼ 0.33 for log N HI = 12.8 in the radiative transfer simulations. The correlation can be fitted with a straight line (cyan dotted line) and shows good agreement between the various models. The transmission spikes in the radiative transfer simulations are produced by regions with slightly larger densities (∆ ∼ 0.33) than those in the optically thin simulations (∆ ∼ 0.25). This is due to the spatial fluctuations in the photo-ionization rate and temperature that are present in the radiative transfer simulations. We discuss the ∆τ vs log N HI correlation for the optically thin and radiative transfer simulations in §4.2.1 and §6.2.1, respectively. correlation is smaller in the hot model. This is because the spikes are smoother and hence there is less variation in ∆τ .

4.2.2
The dependence of spike width on temperature: Tτ vs log b The widths of the components fitted to the spikes are sensitive to the instantaneous temperature along the sightline (see Fig. 1 and Fig. 4). Fig. 6 shows the correlation of optical depth weighted temperature ( Tτ ) with spike width (b) for the cold and hot model. The range in temperature associated with spikes is small for both models. This is expected as the slope of the TDR (Fig. 3) for both models is relatively flat and the temperature associated with ∆ < 1 is relatively constant. Fig. 6 illustrates that the spike widths are well correlated with temperature. The spike widths are systematically larger in the hot model (log b ∼ 1.3) compared to the cold model (log b ∼ 1.05). Even though the range in temperature is small (δ log Tτ ∼ 0.1) for both models, the scatter in log b is relatively large (δ log b ∼ 0.5). In §5.4 we will use the spike width distribution to constrain the temperature of the IGM.

The correlation of spike width and height: log b vs log NHI
The relation of absorption line widths with column density and its relation with the thermal state of the gas at z < 4 has been widely discussed in the literature (Schaye et al. 1999;Bolton et al. 2014;Gaikwad et al. 2017b;Rorai et al. 2017bRorai et al. , 2018Hiss et al. 2018Hiss et al. , 2019. The equivalent relation for the Voigt profile parameters of the transmission spikes in our simulated spectra is compared in Fig. 7 to the observed spectra (white crosses). Fig. 7 shows a strong positive correlation between log b and log NHI for both models. We fit this correlation with a straight line of the form b = b0 [ NHI / 10 12.8 ] α with b0 = (16.65, 10.93) km s −1 ) and α = (0.32, 0.41)) for the hot and cold model, respectively. The hot model is in significantly better agreement with the observations than the cold model. The somewhat flatter slope in the hot model is likely due the flatter TDR in the hot simulation. Note further that the scatter in log NHI is slightly larger in the hot model. This is because the spikes are more blended and hence less distinctive ( Fig. 1 and Fig. 4). The scatter in log b is similar.
In summary, we find strong correlations of the Voigt profile parameters with physical quantities in the optically thin simulations, where : (i) the gas probed by the transmission spikes is typically underdense (∆ ∼ 0.3) and (ii) the spike widths (heights) are strongly (anti-) correlated with temperature (density).

Characterising the flux distribution in transmission spikes
Constraints on cosmological and astrophysical parameters from Lyα forest data have been obtained typically by using either a variety of statistical measures of the Lyα transmitted flux and/or Voigt profile decomposition (Storrie-Lombardi et al. 1996;Penton et al. 2000;McDonald et al. 2000McDonald et al. , 2005Viel et al. 2004bViel et al. , 2009Becker et al. 2011;Shull et al. 2012). Here we briefly consider both approaches before focusing on constraints on the thermal state of the IGM from the width distribution of transmission spikes. The transmitted flux based statistics (such as the probability distribution function and power spectrum) are straightforward to derive from simulations and observations. For the Voigt profile parameter based statistics, we have fitted Voigt profiles to the inverted transmitted flux, 1 − F , using viper 4 . The simulated spectra mimic the observed spectra in terms of S/N and instrumental resolution. of spikes in, respectively, the cold, hot, aton and patchy models at 5.5 < z < 5.7. The b parameter and Tτ are systematically larger in the hot model compared to the cold model. As shown in Fig. 3, Tτ is larger in the aton model compared to the patchy model due to the presence of more gas at ∆ < 1 and T < 30000 K in the aton model. The scatter in temperature in the RT simulation is much larger than that in optically thin simulations due to the presence of UV background and temperature fluctuations. We discuss the Tτ vs (log b) correlation for optically thin and radiative transfer simulations in §4.2.2 and §6.2.2 respectively. with the pseudo-column density ( log N HI ) of the transmission spikes in, respectively, the cold, hot, aton and patchy models at 5.5 < z < 5.7. The correlation is fitted with a straight line (cyan dotted line). The log b parameter at fixed log N HI is systematically larger in the hot model (b ∼ 16.65 km s −1 at log N HI = 12.8) compared to the cold model (b ∼ 10.93 km s −1 at log N HI = 12.8). The slope of the correlation is steeper for the cold (∼ 0.41) model compared to the hot model (∼ 0.32). The b parameters in the aton model (b ∼ 14.96 km s −1 at log N HI = 12.8) are systematically larger than in the patchy model (b ∼ 12.95 km s −1 at log N HI = 12.8). This is because temperatures in the aton model are larger than in the patchy models for ∆ < 1 (see Fig. 6). The white crosses show the scatter in log b and log N HI in the observed spectra. We discuss the (log b) vs log N HI correlation for optically thin and radiative transfer simulations in §4.2.3 and §6.2.3 respectively.
for these effects when determining the best fit parameters, the 1σ statistical uncertainty on best fit parameters and a significance level for each Voigt component. For deriving the spike statistics, we chose only those Voigt components with relative error on parameters ≤ 0.5 and with a significance level ≥ 3. We consider three statistics of the flux distribution in the transmission spikes that are sensitive to astrophysical parameters i.e., ΓHI, T0 and γ (for a given cosmology) which are discussed below.

Spike width (b-parameter) distribution function
For absorption features the line width distribution is frequently used as a diagnostic for the gas temperature, turbulence, and the impact of stellar and AGN feedback on the IGM at z <  Nasir et al. 2017). By contrast, the width of the individual spikes is not a direct measure of the temperature of the low density gas (i.e., b spike = 2kBT /m). Nevertheless we find that the widths of transmission spike components is systematically larger if the temperature of the IGM is larger. 5

Pseudo Column Density Distribution Function (pCDDF)
Similar to the H i column density distribution function (CDDF) at low redshift, we define the pseudo-CDDF (pCDDF) as the number of spikes with a pseudo column density in the range log NHI to δ log NHI in the redshift interval z to z+δz Shull et al. 2012). We calculate the pCDDF in 7 log NHI bins centered at 12.7, 13.1, · · · , 13.9 with dlog NHI = 0.2. This choice of bins is motivated by the S/N and resolution of the observed spectra. The pCDDF characterizes the height and number of spikes in a given redshift bin, and is sensitive to the thermal and ionization parameters of the IGM.

Transmitted flux power spectrum (FPS)
The

Error estimation
We estimate the error for the spike and flux statistics from the simulation using an approach similar to that in Rollinde et al. (2013) and Gaikwad et al. (2018). For this, we generate samples of 5 simulated Lyα forest spectra corresponding to 5 observed QSO spectra with the same observational property i.e., redshift path length, noise property, number of pixels etc. A collection of 5 spectra constitutes a single mock sample. We generate 80 such mock sample and compute the covariance matrix for each statistics. We find that the covariance matrix is converged and dominated by diagonal terms for all the statistics. We have also estimated the covariance matrix from observations using a bootstrap method. Here we find that the off-diagonal terms of the covariance matrix estimated with the bootstrap method are not converged. The bootstrap method (diagonal terms) underestimates the error by ∼ 15 percent as compared to those estimated from the simulations. Throughout this work, we use bootstrap errors increased by a corresponding factor for each statistics. 6 The smallest k mode corresponds to a scale of ∼ 10h −1 cMpc. Fig. 8 compares the spike width distribution, pCDDF and FPS from observations with that of the optically thin simulations for three different redshift bins centered at z = (5.4, 5.6, 5.8). 7 The corresponding residuals suggest a good level of agreement between the default model and observations for all three statistics. The residuals for the hot model with T0 ∼ 11000 K and γ ∼ 1.2 are somewhat larger. Overall the three statistics for the default and hot model are in noticeably better agreement than those for the cold model, with differences in the Doppler parameter distribution being most pronounced where agreement with the default model is significant better than with both the hot and cold model. The cold model also predicts more power on small scales, as well as a larger number of spikes with small log NHI than observed, and is clearly the model that is least consistent with the observations. We further quantify the degree of agreement between models and observation in appendix ??.

Constraining thermal parameters with optically thin simulations
In the optically thin simulations there is a well defined TDR in underdense and moderately overdense regions. It is thus common practice to study the effect of the thermal state on flux statistics by imposing different TDRs by rescaling temperatures (Hui & Gnedin 1997;Rorai et al. 2017bRorai et al. , 2018. In Appendix ?? we use such rescaling to show how the three flux statistics we consider here depend on thermal parameters T0 and γ and the photo-ionization rate, and demonstrate that it is the width distribution of the transmission spikes which is most sensitive to the thermal state of the gas. We further show that, for the reionization and thermal history models we consider in this work, the spike statistics are only very weakly affected by pressure (Jeans) smoothing at the typical gas densities probed by the Lyα transmission spikes (Gnedin & Hui 1998;Theuns et al. 2000;Peeples et al. 2010;Kulkarni et al. 2015;Lukić et al. 2015;Nasir et al. 2016;Maitra et al. 2019;Wu et al. 2019). We use the spike width distribution to obtain an estimate of the best fit values and uncertainty for the thermal parameters 8 . We vary T0 and γ by assuming a TDR of the form T = T0 ∆ γ−1 for ∆ < 10 and T = T0 10 γ−1 for ∆ ≥ 10. We use the ∆ and v fields from the optically thin default model, which falls in between the cold and hot models. We vary T0 between 6000 K to 20000 K in steps of 500 K and γ between 0.4 to 2.0 in steps of 0.05. For each model we (i) compute the Lyα transmitted flux, (ii) post-process the transmitted flux to match observations, (iii) fit the inverted transmitted flux with Voigt profiles and (iv) compute spike statistics. We use 2000 sightlines for each model. with those from cold (red stars), default (orange triangles) and hot (blue squares) optically thin simulations at 5.3 < z < 5.5. The 1σ uncertainties are estimated from the simulated spectra and are shown by the grey shaded regions. Panels A4, A5 and A6 show the corresponding residuals between the models and observations. The gray shaded region in panels A4-A6 corresponds to the gray shaded region in panels A1-A3. Panels B1-B6 and C1-C6 are similar to panels A1-A6 except for the different redshift range, 5.5 < z < 5.7 and 5.7 < z < 5.9, respectively. Note that the hot and default models are in better agreement with observations than the cold model at all redshifts. The agreement between model and observed data is discussed quantitatively in appendix ??. The spike statistics in the default model are very close to those of the best fit model obtained by varying T 0 and γ in §5.4.
Downloaded from https://academic.oup.com/mnras/article-abstract/doi/10.1093/mnras/staa907/5818338 by University of Nottingham user on 17 April 2020 9 shows the 1σ constraints on T0 − γ in three redshift bins. 9 Fig. 9 also shows the marginal distributions of T0 and γ. Table 3 summarizes the best fit values and uncertainty on T0 and γ values under the assumption that the IGM is optically thin. The uncertainty on the T0 and γ measurement accounts for the uncertainty due to continuum fitting, mean flux and Jeans smoothing effects (see Appendix ?? and ??). The best fit values and uncertainty on T0 and γ are consistent with each other within 1σ for the three redshift bins. As the spikes are probing predominantly gas at densities lower than the mean, the inferred values for T0 and γ are strongly correlated. We will come back to this in Appendix ??. Fig. 10 compares the evolution of T0 and γ from this work with a variety of measurements in the literature Bolton et al. 2012Bolton et al. , 2014Boera et al. 2014;Rorai et al. 2017b;Hiss et al. 2018;Walther et al. 2019;Boera et al. 2019). 10 Theoretical models for the evolution of T0 and γ from Puchwein et al. (2019) and Haardt & Madau (2012) are shown by the dashed and dotted curves, respectively 11 . The T0 and γ evolution in these models is obtained for a uniform but time evolving UVB and assuming nonequilibrium ionization evolution. Our T0 and γ constraints in the redshift range 5.3 < z < 5.9 are consistent with the corresponding evolution of the default model in Puchwein et al. (2019) within 1σ. Fig. 10 also shows that the T0 evolution in the hot (cold ) model is systematically larger (smaller) than in the corresponding default model. The errors on T0 and γ account for the statistical and systematic uncertainty (mainly due to continuum fitting). It is interesting to note that the uncertainty on our T0 measurement is smaller than the T0 evolution spanned by the hot and cold models. Our T0 (γ) measurements are higher (lower) than the measurement of Walther et al. (2019). Note, however, that the T0 and γ constraints in Walther et al. (2019) are obtained using the FPS whereas we obtained the T0 and γ constraints from the spike width distribution that is less sensitive to continuum placement and ΓHI uncertainty (see Appendix ??). Further note that the best fit T0 and γ values obtained for the three redshift bins are close to those in the default optically thin simulation. Fig. 8 also shows that the FPS and pCDDF statistics of the default model are consistent with the observations within 1.5 σ. The best fit model obtained by matching the spike width distribution with observations should therefore also have FPS and pCDDF statistics in good agreement with observations. Our γ constraints (γ ∼ 1.2) correspond to a TDR that is moderately steeper than isothermal. However, as we show later, there is no single power law TDR at the redshift of our analysis since the reionization is a patchy inhomogeneous process with different regions reionising at different times. At 5.3 < z < 5.9, γ is therefore not well Table 3. Constraints on T 0 − γ from optically thin simulations

Redshift
T 0 ± δT 0 γ ± δγ 5.3 < z < 5.5 11000 ± 1600 1.20 ± 0.18 5.5 < z < 5.7 10500 ± 2100 1.28 ± 0.19 5.7 < z < 5.9 12000 ± 2200 1.04 ± 0.22 defined in our RT simulations (Keating et al. 2018). In summary, the T0 (γ) constraints obtained in this work are larger (smaller) than those obtained by ). We do not see a significant evolution of T0 and γ in the redshift range 5.3 < z < 5.9. Transmission spikes in optically thin simulations are mainly produced by fluctuations in the density field and the effect of peculiar velocities. Furthermore, the temperature is strongly and tightly correlated with density in optically thin simulations. However, at the redshifts considered here this almost is certainly not realistic. One would expect a large scatter in temperature for a given density as the reionization process will be inhomogeneous, with different regions ionized at different times (Abel & Haehnelt 1999;Miralda-Escudé et al. 2000;Trac et al. 2008;Choudhury et al. 2009). The resulting spatial fluctuations in the amplitude of the UVB and the TDR are not present in optically thin simulations 12 . However, as shown in §2.3, transmission spikes can also be produced by fluctuations in the UVB amplitude and/or temperature. Including these radiative transfer (RT) effects is therefore particularly relevant if reionization ends as late as suggested by the large spatial fluctuations in the Lyα forest opacity (Keating et al. 2019;Kulkarni et al. 2019b).

The radiative transfer simulations in the Sherwood-Relics simulation suite
In addition to the optically thin simulations discussed in §3.1, we have also performed post-processed radiative transfer simulations and hybrid radiative transfer/hydrodynamical simulations as part of the Sherwood-Relics simulation suite. The former simulations model patchy reionization by performing the radiative transfer in post processing on optically thin simulations. This captures many aspects of patchy reionization such as large spatial fluctuations in the photo-ionization rate and temperature, but misses the hydrodynamic response of the IGM to the heating and can thus not accurately predict spatial variations in the pressure smoothing or the distribution of shock heated gas. The hybrid simulations aim to capture these aspects as well.
The post-processed radiative transfer simulations were performed with the GPU-accelerated aton code Aubert & C C C 5.3 < z < 5.5 5.5 < z < 5.7 5.7 < z < 5.9 Figure 9. Panel A shows 1σ constraints on T 0 and γ obtained by comparing the transmission spike width distribution from optically thin simulations with observations at 5.3 < z < 5.5 (red dashed curve), 5.5 < z < 5.7 (green solid curve) and 5.7 < z < 5.9 (blue dotted curve). T 0 and γ are varied in post-processing assuming a power-law TDR (the effect of Jeans smoothing is small, see Fig. ?? in the appendix). Panel B and C shows the marginal distributions for T 0 and γ respectively. The best fit T 0 and γ for 5.3 < z < 5.5 are shown by the red square in panel A and red dashed lines in panel B and C, respectively. Corresponding best fit values for 5.5 < z < 5.7 and 5.7 < z < 5.9 are shown by the green star and solid line and blue circle and dotted line, respectively.
Teyssier (2008), which uses a moment based radiative transfer scheme along with the M1 closure relation. The advection of the radiation was performed using the full speed of light and a single frequency bin for all ionizing photons. Ionizing sources were inserted into dark matter halos as in Kulkarni et al. (2019b). The (mean) energy of ionizing photons was assumed to be 18.6 eV. In the following, we will refer to the post-processed radiative transfer simulation performed on top of our default optically thin simulation by the term aton simulation. It used 2048 3 cells in the (40 h −1 Mpc) 3 box and hydrogen reionization completes at z ≈ 5.2, consistent with the late reionization history found to be favoured by large scale Lyα forest fluctuations (Kulkarni et al. 2019b). The hybrid radiative transfer/hydrodynamical simulation, referred to as the patchy simulation, takes the reionization redshift and H i photo-ionization rate maps produced in the aton simulation as inputs. These are fed to our modified version of p-gadget-3, where they are used in the nonequilibrium thermochemistry solver instead of an external homogeneous UV background model. To obtain consistent density and radiation fields, we use the same initial conditions as in the default optically thin simulation on which the aton run is based. At each timestep, for each SPH particle we check whether it resides in a region in which reionization has already begun. This is assumed to be the case if the ionized fraction in the corresponding cell of the aton simu-lation has exceeded 3 percent. All particles located in such regions are assumed to be exposed to an ionizing radiation field, which is obtained by interpolating the ΓHI maps produced by aton in redshift and reading out the value of the cell containing the particle. This value is then adopted for ΓHI in the non-equilibrium thermochemistry solver. The H i photo-heating rate is computed from ΓHI assuming the same mean ionizing photon energy, 18.6 eV, as in aton. As we do not follow He i and He ii ionizing radiation separately, we use a few simple assumptions to set their photo-ionization and heating rates. For He i we use the same photo-ionization rate as for H i, but we adopt a photo-heating rate that is 30 percent larger than that of H i. For He ii, we use the rates of the fiducial UV background model of Puchwein et al. (2019). This hybrid method results in ionized regions and inhomogeneous photo-heating that closely match those in the parent radiative transfer run, while at the same time following the hydrodynamics and hence including consistent pressure smoothing, as well as shock heating.

The thermal state of the gas in full radiative transfer simulations
The main difference in the radiative transfer simulation is that there are large spatial variations in the TDR (Trac et al. 2008;Keating et al. 2018;Kulkarni et al. 2019b). In panels C and D of Fig. 3 we can compare the TDRs from the aton and patchy RT simulations to the optically thin simulations in panels A and B. Unlike the optically thin simulations, the TDR in the radiative transfer simulations at ∆ < 10 can not be described by a single power law (Bolton et al. 2004). The regions that have been ionized most recently have a flat TDR while regions that have been ionized earlier have progressively steeper TDRs. In the redshift range considered here there are also still significant spatial fluctuations in the amplitude of the UVB. As a result, the variations in temperature for a given ∆ < 10 is large. A crucial difference between the aton and patchy simulations is also evident in Fig.  3. The aton simulation does not account self-consistently for shock heating of the gas. There is thus much less gas with T > 30000 K in the aton simulation than in the patchy simulation. As a consequence there is less gas with ∆ < 10 in the temperature range described by two straight lines (see Fig. 3) in the patchy simulation. As we will show later this has the effect of producing slightly larger spike widths in the aton simulations. Table 2 also shows the physical effects responsible for the occurrence of spikes in the aton and patchy radiative transfer simulations. Similar to the optically thin simulations, most of the spikes (∼ 89 percent) in the aton and patchy simulations occur in underdense regions and for ∼ 18 percent of spikes the gas shows a diverging velocity field along the sightline. However, in contrast to optically thin simulations, around 50 percent of spikes show an enhancement of ΓHI and temperature. Note here that the recent observations of high redshift Lyman alpha emitters and Lyman break galaxies suggests that the transmission spikes are spatially correlated with the ionizing radiation escaping these galaxies (Meyer et al. 2019). Thus, for the transmission spikes in the full radiative transfer simulations, all the physical processes discussed in §2.3 and Fig. 2 contribute to the occurrence of spikes.
Bolton+2014 Telikova+2019 This Work (Gaikwad+2020) Figure 10. The evolution of thermal parameters T 0 and γ from the literature and in this work (red stars with error bars) are shown in the top and bottom panels, respectively Bolton et al. 2012Bolton et al. , 2014Boera et al. 2014;Rorai et al. 2017b;Hiss et al. 2018;Walther et al. 2019;Boera et al. 2019;Telikova et al. 2019). Note that the temperature constraints from Becker et al. (2011) and Boera et al. (2014) are not measured at the mean density, and have therefore been scaled to a T 0 value assuming a TDR with slope γ = 1.3. The T 0 and γ evolution in the hot, default and cold Sherwood-Relics simulations is shown by the blue dash-dotted, black dashed and red dotted curves respectively. The corresponding T 0 and γ evolution in the Haardt & Madau (2012) UVB synthesis model is shown by a green dotted curve. This T 0 and γ evolution is obtained by assuming a uniform UVB and solving for the non-equilibrium ionization evolution (Haardt & Madau 2012;Puchwein et al. 2019). Our constraints on T 0 and γ are consistent within 1σ with Puchwein et al.  Fig. 5 we compare the dependence of ∆τ on log NHI for the aton and patchy simulations to the optically thin cold and hot simulations. Similar to the optically thin simulations (Fig. 5), ∆τ and log NHI are anti-correlated. The normalization (∆0) and slope (β) are similar in both simulations. The ∆0 (β) in the RT simulations are slightly larger (smaller) compared to the optically thin simulations. The spikes (irrespective of log NHI) in the RT simulations are produced from underdensities somewhat larger than in the optically thin simulations. This is expected, as spikes (at a given log NHI) in the RT simulations are produced by all four physical effects we discussed previously, i.e., fluctuations in density, peculiar velocity, UVB and temperature. Fig. 5 also shows that the scatter in ∆τ (at a given log NHI) is larger in the radiative transfer simulations due to fluctuations in UVB, temperature and pressure smoothing effects. Note, however, that in all simulations the spikes occur in underdense regions with ∆ < 1.
6.2.2 Dependence of temperature on spike width: Tτ vs b Fig. 6 compares the dependence of optical depth weighted temperature ( Tτ ) on spike widths (b-parameter) for the RT simulations with that in the optically thin simulations. Unlike the optically thin simulations, the RT simulations show a large scatter in temperature for a given spike width due to fluctuations in the UVB amplitude and temperature. Furthermore, the temperature in the patchy simulation is smaller than in the corresponding aton simulation. This is because (i) the amount of gas with ∆ < 1 and T < 30000 K is larger in the patchy simulation and (ii) due to the post-processed nature of the aton simulation the (adiabatic) change in the temperature due to changes in density (the d∆/dt term) is not accounted. As a result, the spike widths are also slightly smaller in the patchy simulations.
6.2.3 The relation of spike width and height: b vs log NHI Fig. 7 compares the b− log NHI correlation for the aton and patchy simulations to that in the optically thin simulations. Similar to the optically thin simulations, log b and log NHI are strongly anti-correlated in the RT simulations. The nor-malization of the correlation b0 is smaller in the patchy (∼ 12.95 km s −1 ) simulation than in the aton (∼ 14.96 km s −1 ) simulation due to the smaller temperature of the gas probed by the transmission spikes in the latter. The slope of the correlation (α = 0.32 for aton and α = 0.34 for patchy) and the scatter in the TDR (at ∆ < 1 in Fig. 3) is relatively similar for both simulations. It is interesting to note here that α in the RT simulations is similar to that in the hot optically thin simulation, while b0 in the RT simulations is smaller than in the hot optically thin simulations. The smaller value of b0 in the RT simulation is a consequence of fluctuations in UVB amplitude and temperature. In summary, the RT simulations include the effects of fluctuations in the UVB amplitude and temperature which are missing in the optically thin simulations. Due to these effects: (i) the TDR in RT simulations cannot be described by a single power-law, (ii) the typical densities responsible for transmission spikes in the RT simulations (∆0 ∼ 0.33) is slightly larger than in optically thin simulations (∆0 ∼ 0.25) and (iii) the scatter in temperature and spike width are larger.

Transmission spike properties: Full radiative transfer simulations vs observations
We now compare the three statistics from the RT simulations with observations in Fig. 11. Each panel is similar to Fig. 8, except we now use the aton and patchy simulations. Similar to the optically thin simulations, the mean transmitted flux in the aton and patchy models is matched to that in the observed spectra to account for the uncertainty in the continuum placement and UV background amplitude. Note that this rescaling has surprisingly little effect on the width of the transmission spikes (see appendix ??). Comparison of Fig. 11 with Fig. 8 shows that the RT simulations are in perhaps even better agreement with the observations than the optically thin simulations at all redshifts. Note, however, that there is still considerable freedom to adjust the overall temperature in both the RT and optically thin simulation. Unfortunately, we therefore don't think that this (marginally) better agreement should be interpreted as (hard) evidence for the spatial variations of the TDR predicted by our RT simulations. The three statistics are furthermore very similar in the aton and patchy simulations at all redshifts. The spike widths in the aton simulations are somewhat larger than in the patchy simulation and are in marginally better agreement with observations than in the patchy simulations (see appendix ?? for a comparison of the goodness of fit for the different models). As explained in the previous section, this is due to the effect of slightly higher gas temperatures in the aton simulations. Note, however, also that both RT simulations are mono-frequency and the normalisation of the temperature distribution predicted by the RT simulations is still somewhat uncertain. Fig. 12 shows the comparison of the T0 and γ evolution with that from T0 − γ measured assuming a uniform UVB. Since there is no single power-law TDR in the RT simulations (see Fig. 3), we show a range in T0 and γ evolution. To obtain this range in T0 and γ, we find the 16 th and 84 th percentile temperature in four ∆ bins. We then fit a power-law TDR to the 16 th and 84 th percentile temperature values (see Fig. 3 13 ). Fig. 12 shows that the T0 − γ constraints obtained here are consistent with the range in the T0 − γ evolution seen in the patchy and aton RT simulations.
Thus, quite remarkably the patchy simulation based on the self-consistent reionization model of Kulkarni et al. (2019b) that (i) simulates cosmological density and velocity fields, (ii) includes spatial fluctuations in the UV background and TDR, (iii) accounts for the pressure smoothing of gas and (iv) matches the Thomson scattering optical depth (Planck Collaboration et al. 2018), also produces transmission spike properties consistent with those in observed high-resolution, high-z QSO absorption spectra. Note, however, that there is still some uncertainty in the post-reionization temperatures in RT simulations (see D'Aloisio et al. 2019, Puchwein et.al. 2020

CONCLUSIONS
We have explored here for the first time the use of the transmission spikes observed in high-z QSO absorption spectra as a tool to probe the physical state of the IGM near the tail end of hydrogen reionization. We constrain the thermal state of the IGM at 5.3 < z < 5.9 by comparing the properties of Lyα transmission spikes from a sample of 5 high resolution (vFWHM ∼ 6 km s −1 ) and high S/N (∼ 10) QSO absorption spectra with that from state-of-the-art, high resolution optically thin simulations run with gadget-3 (Springel 2005) from the Sherwood and Sherwood-Relics suites, as well as a simulation post-processed with the radiative transfer code aton (Aubert & Teyssier 2008). The main results of this work are as follows.
• In full radiative transfer simulations regions with low density, enhancement in the photo-ionization rate ΓHI, enhancement in temperature and a diverging peculiar velocity field along the line of sight can all contribute to the occurrence of transmission spikes at high redshift. Most of the spikes (∼ 90 per cent) in optically thin and radiative transfer simulations occur in regions with density ∆ < 1. Optically thin simulations do not account for the effect of enhanced temperatures in recently ionized regions and the resulting spatial fluctuations in the temperature-density relation. Due to the assumed spatially homogeneous UV background amplitude they also do not account for the occurrence of transmission spikes due to enhancements in the photo-ionization rate. About 50 per cent of the transmission spikes in our RT simulations show the effect of an enhanced ΓHI and enhanced temperature. In the RT simulation the transmission spikes are often due to either hot, recently ionized, very underdense regions with a diverging line of sight peculiar velocity field or due to somewhat less underdense and colder regions with an enhanced photo-ionization rate.
• The width of the components fitted to the asymmetric and blended transmission spikes are very sensitive to This Work (Gaikwad+2020) Figure 12. Same as Fig. 10, except the T 0 and γ evolution is shown from the patchy RT simulation (shaded region) at 4 ≤ z ≤ 8. Since there is no single power-law TDR in the RT simulations (see Fig. 3), the shaded region displays the 16 th and 84 th percentiles of T 0 and γ (see §6.3). For comparison, we also show the T 0 and γ evolution from RT models (50 th percentile, blue solid line), the uniform UVB default (red dashed line) and Haardt & Madau (2012) (green dotted line) model. The T 0 and γ measured by comparing the observed spike width distribution with that from optically thin simulations are in good agreement with the median T 0 and γ evolution from RT simulations.
the instantaneous temperature of the gas and are significantly broader in the optically thin hot simulation than in the corresponding cold simulation. To quantify this, we have fitted multi-component Voigt profiles to the inverted transmitted flux 1 − F in both simulated and observed spectra with our automated code viper (Gaikwad et al. 2017b). We derive the transmitted flux power spectrum, pseudo column density distribution function (pCDDF) and spike width (b-parameter) distribution functions for simulated and observed spectra. We show that the spike width distribution is the statistic that is sensitive to the thermal state of the IGM. The dependence of the shape of the FPS and pCDDF on the temperature of the absorbing gas is somewhat weaker while their normalisation is more sensitive to the ionization state/neutral fraction of the IGM.
• We associate the observable properties of spikes with the physical properties of gas in simulations by studying the ∆τ − log NHI, Tτ −b and b− log NHI correlations. These correlations show that the underdensity of gas associated with spikes is similar i.e., ∆τ ∼ 0.3 at log NHI ∼12.8 in both optically thin and radiative transfer simulations. The spike widths in both simulations are sensitive to the temperature of the gas (b ∼ 10.9 km s −1 for T0 ∼ 7500 K and b ∼ 16.6 km s −1 for T0 ∼ 14000 K). However, the temperature scatter for a given density is larger in the radiative transfer simula-tion compared to the optically thin simulation. As a result, a significant fraction of spikes in the radiative transfer simulations are due to hotter temperatures in recently ionized regions of the Universe.
• We have compared the three observed flux statistics with those derived from our optically thin hot and cold simulations in three redshift bins. The statistics for the hot and default model are in significantly better agreement with those from observations in all the 3 redshift bins. The Doppler parameter distribution which is most sensitive to the instantaneous temperature of the gas is in significantly better agreement for the default model than for both the cold and hot model. We constrain thermal parameters by varying T0 and γ in post-processed optically thin simulations. The best fit values at 5.3 ≤ z ≤ 5.5, 5.5 ≤ z ≤ 5.7, 5.7 < z < 5.9 are T0 ∼ 11000 ± 1600, 10500 ± 2100, 12000 ± 2200 K and γ ∼ 1.20 ± 0.18,1.28 ± 0.19, 1.05 ± 0.22 respectively. We do not find significant evolution in T0 and γ over 5.3 < z < 5.9.
• We have also compared the three statistics in physically motivated radiative transfer simulations with those from observations. Unlike optically thin simulations, radiative transfer simulations incorporate spatial fluctuations in the amplitude of the UVB and the TDR. As a result the scatter in the TDR is large and a single power-law cannot describe the TDR in radiative transfer simulations. The observed spike Downloaded from https://academic.oup.com/mnras/article-abstract/doi/10.1093/mnras/staa907/5818338 by University of Nottingham user on 17 April 2020 statistics in our radiative transfer simulation of late reionization with neutral islands persisting to z ∼ 5.3 are in good agreement (see appendix ?? for details) with observations in all three redshift bins .
Our work shows the potential of transmission spike shapes (heights and widths) for constraining the thermal history of the IGM near the tail-end of H i reionization, complimentary to other transmitted flux based methods. In future, a much larger sample of high resolution, high S/N and high-z QSO absorption spectra should become available thanks to 30-40 m class optical telescopes. These larger data sets, complemented by further improved radiative transfer simulations, promise to put tight constraints on the nature and the exact timing of H i reionization.