The detectability of strong 21 centimetre forest absorbers from the diffuse intergalactic medium in late reionisation models

A late end to reionisation at redshift $z\simeq 5.3$ is consistent with observed spatial variations in the Ly$\alpha$ forest transmission and the deficit of Ly$\alpha$ emitting galaxies around extended Ly$\alpha$ absorption troughs at $z=5.5$. In this model, large islands of neutral hydrogen should persist in the diffuse intergalactic medium (IGM) until $z\simeq 6$. We use a novel, hybrid approach that combines high resolution cosmological hydrodynamical simulations with radiative transfer to predict the incidence of strong 21 cm forest absorbers with optical depths $\tau_{21}>10^{-2}$ from the diffuse IGM in these late reionisation models. We include the effect of redshift space distortions on the simulated 21 cm forest spectra, and treat the highly uncertain heating of the pre-reionisation IGM by soft X-rays as a free parameter. For a model with only modest IGM pre-heating, such that average gas kinetic temperatures in the diffuse IGM remain below $T_{\rm K}\simeq 10^{2} \rm\, K$, we find that strong 21 cm forest absorption lines should persist until $z=6$. For a sample of $\sim 10$ sufficiently radio loud background sources, a null-detection of 21 cm forest absorbers at $z\simeq 6$ with SKA1-low or possibly LOFAR should provide an informative lower limit on the still largely unconstrained soft X-ray background at high redshift and the temperature of the pre-reionisation IGM.


INTRODUCTION
At present, the premier technique for examining the small-scale structure of intergalactic neutral hydrogen approaching the reionisation era is Lyman series absorption in the spectra of luminous quasars Eilers et al. 2017;Bosman et al. 2018;Yang et al. 2020b). However, it is challenging to probe the intergalactic medium (IGM) much beyond redshift 6 with this approach. The large cross-section for Ly scattering means the IGM becomes opaque to Ly photons at neutral hydrogen fractions as low as HI 10 −4 . An alternative transition that overcomes this limitation is the hyperfine 21 cm line, which has a cross-section that is a factor ★ E-mail: tomas.soltinsky@nottingham.ac.uk ∼ 10 7 smaller than the Ly transition 1 . If radio bright sources such as high redshift quasars (Bañados et al. 2021) or gamma-ray bursts (e.g. Ioka & Mészáros 2005;Ciardi et al. 2015b) can be identified during the reionisation era, the intervening neutral IGM may be observed as a 21 cm forest of absorption lines in their spectra. This can be achieved either through the direct identification of individual absorption features (Carilli et al. 2002;Furlanetto & Loeb 2002;Meiksin 2011;Xu et al. 2011;Ciardi et al. 2013;Semelin 2016;Villanueva-Domingo & Ichiki 2021), or by the statistical detection of the average 21 cm forest absorption (Mack & Wyithe 2012;Ewall-Wice et al. 2014;Thyagarajan 2020). This approach is highly complementary to proposed tomographic studies of the redshifted 21 cm line and measurements of the 21 cm power spectrum during reionisation (e.g. Mertens et al. 2020;Trott et al. 2020), as it is subject to a different set of systematic uncertainties Pritchard & Loeb 2012).
However, any detection of the 21 cm forest relies on the identification of sufficient numbers of radio-loud sources and the existence of cold, neutral gas in the IGM at 6. While neither of these criteria are guaranteed, the prospects for both have improved somewhat in the last few years. Approximately ∼ 10 radio-loud active galactic nuclei are now known at 5.5 < < 6.5 (e.g. Bañados et al. 2018bBañados et al. , 2021Liu et al. 2021;Ighina et al. 2021), including the = 6.1 blazar PSO J0309+27 with a flux density 147 MHz = 64.2±6.2 mJy (Belladitta et al. 2020). The Low Frequency Array (LOFAR) Twometre Sky Survey (LoTSS, Shimwell et al. 2017;Kondapally et al. 2020), the Giant Metrewave Radio Telescope (GMRT) all sky radio survey at 150 MHz (Intema et al. 2017), and the Galactic and Extragalactic All-sky Murchison Widefield Array survey (GLEAM, Wayth et al. 2015) are also projected to detect hundreds of bright > 6 radio sources, particularly if coupled with large spectroscopic follow-up programmes such as the William Herschel Telescope Enhanced Area Velocity Explorer (WEAVE)-LOFAR survey (Smith et al. 2016).
Furthermore, there is now growing evidence that reionisation ended rather late, and possibly even extended to redshifts as late as 5.3 Nasir & D'Aloisio 2020;Qin et al. 2021). This picture is motivated by the large spatial fluctuations observed in the Ly forest transmission at 5.5 Eilers et al. 2018). A late end to reionisation is also consistent with the electron scattering optical depth inferred from the cosmic microwave background (Pagano et al. 2020), the observed deficit of Ly emitting galaxies around extended Ly absorption troughs Kashino et al. 2020;Keating et al. 2020), the clustering of Ly emitters (Weinberger et al. 2019), the thermal widths of Ly forest transmission spikes at > 5 ( Gaikwad et al. 2020), and the mean free path of ionising photons at = 6 (Becker et al. 2021). If this interpretation proves to be correct (but see D'Aloisio et al. 2015;Davies & Furlanetto 2016;Chardin et al. 2017;Meiksin 2020, for alternative explanations), then there should still be large islands of neutral hydrogen in the IGM as late as = 6 (e.g. Lidz et al. 2007;Mesinger 2010). If this neutral gas has not already been heated to kinetic temperatures K 10 3 K by the soft X-ray background, then it may be possible to detect 21 cm absorbers in the pre-reionisation IGM at 6. Alternatively, a null-detection could provide a useful lower limit on the soft X-ray background at high redshift.
The goal of this work is to investigate this possibility further. We use a set of high resolution hydrodynamical cosmological simulations drawn from the Sherwood-Relics simulation programme (Puchwein et al. in prep). Using a novel hybrid approach, these are combined with the ATON radiative transfer code (Aubert & Teyssier 2008) to model the small-scale structure of the IGM. Following Kulkarni et al. (2019), we consider a model with late reionisation ending at = 5.3, and contrast this with a simulation that has an earlier end to reionisation at = 6.7. We pay particular attention to some of the common assumptions adopted in previous models of the 21 cm forest that affect the absorption signature on small scales. This includes the treatment of gas peculiar motions and thermal broadening, the coupling of the spin temperature to the Ly background, and the effect of pressure (or Jeans) smoothing on the IGM. Our approach is therefore closest to the earlier work by Semelin (2016), although we do not follow spatial variations in the X-ray and Ly backgrounds. Instead, we attempt to explore a broader range of parameter space for spatially uniform X-ray heating using hydrodynamical simulations that use several different reionisation histories and have an improved mass resolution (by a factor ∼ 27). Note, however, that even the high resolution cosmological simulations considered here will still only capture the 21 cm absorption that arises from the diffuse IGM. We therefore do not model the (uncertain amount of) absorption from neutral gas in haloes below the atomic cooling threshold, or from the cold interstellar medium in much rarer, more massive haloes that host high redshift galaxies (see e.g. Furlanetto & Loeb 2002;Meiksin 2011). This paper is structured as follows. We start by describing our numerical model of the IGM in Section 2, and our calculation of the 21 cm optical depths in Section 3. We examine how different modelling assumptions affect the observability of strong 21 cm forest absorbers in Section 4.1 and 4.2, and estimate how a null-detection of strong 21 cm absorbers at redshift 6 with LOFAR or the Square Kilometre Array (SKA) could constrain the high redshift soft X-ray background in Section 4.3. Finally, we conclude in Section 5. The appendix contains further technical details regarding our methodology and modelling assumptions.

Hydrodynamical simulations with radiative transfer
We model the 21 cm forest during inhomogeneous reionisation using a sub-set of the high resolution cosmological hydrodynamical simulations drawn from the Sherwood-Relics simulation programme (see Gaikwad et al. 2020, for an initial application of these models). The Sherwood-Relics simulations were performed with a modified version of the P-GADGET-3 code -which is itself an updated version of the GADGET-2 code described in (Springel 2005) -and uses the same initial conditions as the earlier Sherwood simulation suite (Bolton et al. 2017). In this work we adopt a flat ΛCDM cosmology with Ω Λ = 0.692, Ω m = 0.308, Ω b = 0.0482, 8 = 0.829, s = 0.961, ℎ = 0.678, consistent with Planck Collaboration et al. (2014), and a primordial helium fraction by mass of p = 0.24 (Hsyu et al. 2020). The simulations have a volume (40ℎ −1 cMpc) 3 and track 2 × 2048 3 dark matter and gas particles. This yields a dark matter particle mass of 7.9 × 10 5 and resolves dark matter haloes with masses greater than ∼ 2.5 × 10 7 . This high mass resolution is necessary for capturing the small-scale intergalactic structure probed by the 21 cm forest (cf. Semelin 2016). We furthermore adopt a simple but computationally efficient scheme for converting high density gas into collisionless particles that robustly predicts the properties of the IGM. If a gas particle has an overdensity Δ = 1 + > 1000 and kinetic temperature K < 10 5 K, it is converted into a collisionless star particle (Viel et al. 2004). We have verified this simplified approach is sufficient for modelling the 21cm forest in the diffuse IGM by direct comparison to a full sub-grid star formation model (see Appendix A for further details). The main effect of this approximation is the removal of dense gas from haloes, which slightly reduces the number of strong 21 cm absorbers in models with no X-ray heating.
In order to include the effect of inhomogeneous reionisation by UV photons on the IGM, the Sherwood-Relics simulations are combined with the moment-based, M1-closure radiative transfer code ATON (Aubert & Teyssier 2008). We adopt a novel hybrid approach that captures the small-scale hydrodynamical response of the gas in the simulations to patchy heating during reionisation Table 1. Hydrodynamical simulations used in this work. From left to right, the columns give the name of the simulation, the nature of the reionisation model, the redshift when the IGM is fully reionised, and the prescription for converting dense gas into collisionless star particles, which follows either Viel et al. (2004) or Puchwein & Springel (2013). The first three simulations are part of the Sherwood-Relics suite. The final two models are optically thin simulations with rapid reionisation at 15, and are described in Appendix A. All models have a volume of (40ℎ −1 cMpc) 3 and include 2 × 2048 3 dark matter and gas particles. (see also Oñorbe et al. 2019, for a related approach). Our hybrid RT/hydrodynamical simulations use inputs in the form of 3D maps of the reionisation redshift and H I photo-ionisation rate, produced by ATON simulations performed on the P-GADGET-3 outputs in post-processing. These maps are then fed back into a re-run of the P-GADGET-3 model, where they are called within a non-equilibrium thermo-chemistry solver (Puchwein et al. 2015). Following Kulkarni et al. (2019), the ionising sources in the ATON simulations have luminosities proportional to the halo mass with a redshiftdependent normalization, and the mean energy of ionising photons is assumed to be 18.6 eV. Further details can be found in Gaikwad et al. (2020) and Puchwein et al. (in prep). The main advantage of this approach is that since the post-processing step using the ATON radiative transfer simulations is computationally cheap compared to the hydrodynamical simulations, we may empirically calibrate the source model to yield a reionisation history that is consistent with a wide range of observational constraints. This avoids many of the uncertainties associated with direct hydrodynamical modelling of the source population. Note, however, the relatively small volume of our simulations means the patchy structure of reionisation will not be fully captured on scales larger than the box size. This will lead to smaller neutral islands and an earlier percolation of ionised regions relative to simulations performed in a larger volume (Iliev et al. 2014;Kaur et al. 2020). We therefore adjust the ionising emissivity in the models by hand to achieve a given reionisation history; this scaling is equivalent to varying the escape fraction of ionising photons. In addition, while our ATON simulations self-consistently follow the propagation of ionising photons using a 2048 3 Cartesian grid, self-shielded regions below the size of the grid cells ( 20ℎ −1 ckpc) will not be resolved. We attempt to partially correct for this by implementing a correction for the self-shielding of dense gas in all our simulations in post-processing, using the results of Chardin et al. (2018). We find, however, that this correction makes almost no difference to our final results, as the majority of the strong 21 cm absorbers in our simulations arise from the diffuse IGM.
We consider two different reionisation histories in this work, where the IGM becomes fully ionised at either = 5.3 (model zr53) or = 6.7 (model zr67). The filling fraction of ionised gas and the Thomson scattering optical depth predicted by these models are displayed in Fig. 1. Both of these models are consistent with current observational constraints on the timing of reionisation. As already discussed, the reionisation model that ends at = 5.3 is furthermore consistent with the large fluctuations in the Ly forest transmission The filling factor of ionised gas in the zr53 (blue solid curve) and zr67 (fuchsia solid curve) simulations, compared to observational constraints from dark gaps in the Ly forest (McGreer et al. 2015), the damping wing in high redshift quasar spectra (Bañados et al. 2018a;Davies et al. 2018;Greig et al. 2017;Greig et al. 2019;Wang et al. 2020;Yang et al. 2020a) and Ly emitting galaxies (Mason et al. 2018(Mason et al. , 2019. For clarity of presentation a small offset has been added to the redshifts of some of the data points. Bottom: The Thomson scattering optical depth to cosmic microwave background photons. The black line with the shaded region shows the Planck Collaboration et al. (2020) measurement. observed at 5.5 Kulkarni et al. 2019;Keating et al. 2020). Finally, we also use a simulation (zr53-homog) that has been constructed to give exactly the same globally averaged reionisation history as the zr53 model, but using a spatially uniform ionising background. A comparison between the zr53 and zr53homog simulations therefore allows us to estimate the uncertain effect that pressure smoothing (from e.g. UV photo-heating) may have on the gas in the pre-reionisation IGM (see Section 4.2). All the simulations used in this work are listed in Table 1, where the final two models listed are used in Appendix A only.

Heating of neutral gas by the X-ray and Ly backgrounds
Absorption features in the 21 cm forest arise from neutral hydrogen in the IGM. In addition to modelling the inhomogeneous reionisation of the IGM by UV photons, we must therefore also consider the temperature and ionisation state of gas that is optically thick to Lyman continuum photons. This heating is attributable to adiabatic compression and shocks -which are already included within our hydrodynamical simulations -and the X-ray and (to a lesser extent) Ly radiation backgrounds at high redshift (Ciardi et al. 2010), which are not. Hence, we now describe the procedure we use to include spatially uniform X-ray and Ly heating in our simulations, by recalculating the density dependent temperature and ionisation state of the neutral gas in our hybrid simulations in post-processing. Further details on the model we use are also provided in Appendix B.
As we do not directly model the star formation rate in our simulations, rather than using a detailed model for the number and spectral energy distribution of X-ray sources at high redshift, for simplicity and ease of comparison to the existing literature we instead follow the approach introduced by Furlanetto (2006b) for parameterising the comoving X-ray background emissivity. This uses the observed relationship between the star formation rate, SFR, and hard X-ray band luminosity (2-10 keV) for star-forming galaxies at = 0 (Gilfanov et al. 2004;Lehmer et al. 2016). Furlanetto (2006b) adopt the normalisation X = 3.4 × 10 40 erg s −1 for the total X-ray luminosity at photon energies > 0.2 keV, assuming a power-law spectral index X = 1.5. The X-ray efficiency, X , parameterises the large uncertainty in the extrapolation of Eq. (1) toward higher redshift. Using the conversion X,0.2 keV = X ( X − 1)/ 0.2keV , the corresponding comoving Xray emissivity is We assume a power-law spectrum with X = 1.5, and use the fit to the observed comoving star formation rate density from We assume that SFR = 0 at redshifts > ★ = 14, and have verified that adopting ★ > 14 does not change our predictions for 21 cm absorption at ≤ 12. The UV background emissivity at the Ly wavelength from stars in our model is instead given by where we have used the conversion between SFR and UV luminosity at 1500 Å from Madau & Dickinson (2014) and assumed a flat spectrum ( ∝ 0 ) in the UV, where the Ly efficiency parameterises the uncertain amplitude. We adopt = 1 as the fiducial value in this work, but note that this parameter is uncertain and the Ly emissivity should furthermore vary spatially (see e.g. Fig. 4 in Semelin 2016). For illustrative purposes we therefore also show some results for the much smaller value of = 0.01 (but note that in practise and the reionisation history are not fully decoupled). The primary effect of increasing (decreasing) the Ly efficiency is to produce a tighter (weaker) coupling of the H I spin and kinetic temperatures. A smaller value of may be more appropriate for absorbers that are distant from the sources of Ly background photons. Instead of a flat UV spectrum we also considered the = 1 (solid black curve) and the specific intensity of the X-ray background for photon energies 0.2 keV (blue curve), 1 keV (fuchsia curve) and 2 keV (orange curve), assuming an X-ray efficiency of X = 1. The X-ray specific intensities have been multiplied by a factor of 10 5 for presentation purposes. For comparison, the dashed curves show the X-ray specific intensities evaluated in the optically thin limit.
power-law population II and III spectra used by , but the strength of the Ly coupling in our model is not very sensitive to this choice at the redshifts we consider.
With these emissivities in hand, we may evaluate the solution to the cosmological radiative transfer equation (see Eq. (B1) in Appendix B) to obtain the X-ray specific intensity at photon energies 0.2-30 keV (Pritchard & Loeb 2012). Similarly, we obtain the specific intensity of the Ly background by evaluating Eq. (B7), following . Fig. 2 shows the redshift evolution of the specific intensity of the Ly background from stellar emission, ,★ ( ), and the specific intensity of the X-ray background at three different energies, 0.2 keV, 1 keV and 2 keV. The dashed curves show the X-ray specific intensities in the optically thin limit, i.e. when the optical depth of the intervening IGM to X-ray photons is set to zero in Eq. (B1). Note that 2.0 keV remains almost unchanged in the optically thin limit, implying the IGM is transparent to photons emitted with energies ≥ 2 keV at < ∼ 10 (cf.

McQuinn 2012).
The unresolved soft X-ray background at = 0 places an upper limit on the contribution of high redshift sources to the hard X-ray background, since these photons may redshift without significant absorption to = 0 (Dĳkstra et al. 2004;Mc-Quinn 2012). When assuming X = 1.8, integrating our model specific intensity in the soft X-ray band (0.5-2 keV) at = 0 yields 0.5−2keV = 2.9 × 10 −12 erg s −1 cm −2 deg −2 . This value is consistent with the unresolved soft X-ray background obtained from Chandra observations of the COSMOS legacy field, 0.5−2keV = 2.9±0.16×10 −12 erg s −1 cm −2 deg −2 (Cappelluti et al. 2017). Note, however, the = 0 soft X-ray background does not provide a direct constraint on the very uncertain soft X-ray background at high redshift (see e.g. Dĳkstra et al. 2012;Fialkov et al. 2017). Recently, Greig et al. (2021a) have presented the first weak, model . The redshift evolution of the gas kinetic temperature, K , (solid curves) and spin temperature, S , (dashed curves) at mean density following X-ray background heating by photons with = 0.2-30 keV. The different coloured curves correspond to efficiency parameter X = 10 (green curves), 1 (orange curves), 0.1 (fuchsia curves) and 0.01 (blue curves). For comparison, the CMB temperature, CMB = 2.73 K(1 + ) corresponds to the dot-dashed curve, and the kinetic temperature for adiabatic heating and cooling only, ad,0 = 2.73 K(1 + z) 2 /(1 + dec ), is shown by the dotted curve. We assume the gas thermally decouples from the CMB at dec = 147.8 ). The filled diamonds and grey shading correspond to the gas kinetic temperature measurements from Ly transmission spikes in quasar spectra (Gaikwad et al. 2020) and Ly absorption lines in quasar proximity zones (Bolton et al. 2012), respectively. The filled circles show the model dependent lower limits on the H I spin temperature obtained from LOFAR (Greig et al. 2021b) and MWA (Greig et al. 2021a). dependent lower limits on the soft X-ray background emissivity at 6.5 ≤ ≤ 8.7 using the Murchison Widefield Array (MWA) upper limits on the 21 cm power spectrum (Trott et al. 2020), where X,0.5−2 keV 10 34.5 erg s −1 cMpc −3 . For comparison, for an Xray efficiency of X = 0.01, our X-ray background model gives X,0.5−2 keV = 10 36.0 erg s −1 cMpc −3 at = 8.1, which is well above the Greig et al. (2021a) lower limit.
Given the specific intensities of the X-ray and Ly radiation backgrounds, we next compute the thermal evolution of the IGM that remains optically thick to UV photons, but is heated by X-ray and Ly backgrounds that are assumed to be spatially uniform on the scale of our simulated volume. 2 We follow the procedure described in Appendix B for this purpose. Fig. 3 displays the temperature evolution of a gas parcel at mean density for four different values of the X-ray efficiency parameter X . An approximate lower limit on X is provided by the recent constraints on the spin temperature from upper limits on the 21 cm power spectrum at 9.1 obtained with LOFAR (Mertens et al. 2020), and at = 6.5−8.7 from MWA (Trott et al. 2020). These data disfavour very cold reionisation models with 2 The mean free path to X-ray photons is X = 5 cMpc x −1 HI (1 + ) −1 (E/0.2 keV) 3 [ (1 + z)/10] −2 . Fluctuations in the temperature of soft X-ray heated gas on ∼ 10 cMpc scales are thus expected (Pritchard & Furlanetto 2007;Ross et al. 2017;Eide et al. 2018). These fluctuations would not, however, be adequately captured in our small simulation volume. no X-ray heating (Greig et al. 2021b;Ghara et al. 2020;Greig et al. 2021a). An approximate upper limit on X at > 6 is provided by Ly absorption measurements of the kinetic temperature at 5-6, after the IGM has been photo-ionised and heated by UV photons (Bolton et al. 2012;Gaikwad et al. 2020). These data are consistent with X 10. Adopting larger X-ray efficiencies in our model would overheat the low density IGM by = 6.

THE 21 CM FOREST OPTICAL DEPTH
We now turn to the calculation of the 21 cm optical depth. The 21 cm line arises from the hyperfine structure of the hydrogen atom, and is determined by the relative orientation of the proton and electron spin, where the ground state energy level is split into a singlet and triplet state. A photon with rest-frame wavelength 21 = 21.11 cm, or equivalently frequency 21 = 1420.41 MHz, can induce a transition between these two states.
In the absence of redshift space distortions, the optical depth to 21 cm photons at redshift is = 0.27 HI 1 + 10 where HI is the H I number density, S is the spin temperature, 10 = 2.85 × 10 −15 s −1 is the Einstein spontaneous emission coefficient for the hyperfine transition, is the gas overdensity and ( ) is the Hubble parameter (Madau et al. 1997). Note the factor of 0.27 in the second equality is cosmology dependent. Absorption will therefore be most readily observable for dense, cold and significantly neutral hydrogen gas. The H I spin temperature, a measure of the relative occupation numbers of the singlet and triplet states, is (Field 1958 where CMB = 2.73(1 + ) K is the temperature of the cosmic microwave background (CMB, Fixsen 2009), is the Ly colour temperature and c , are the coupling coefficients for collisions and Ly photon scattering, respectively. If c + 1, the H I spin temperature is coupled to the gas kinetic temperature, and if c + 1 it is coupled to the CMB temperature. The collisional coupling coefficient is where ★ = ℎ p 21 / B , and HH 10 , eH 10 , pH 10 are the temperature dependent de-excitation rates for collisions between hydrogen atoms, electrons and hydrogen atoms, and protons and hydrogen atoms, respectively. We use the convenient fitting functions to the deexcitation rates from Kuhlen et al. (2006) and Liszt (2001), modified to better agree with tabulated values for HH 10 , eH 10 (Furlanetto & Furlanetto 2007a), and pH 10 (Furlanetto & Furlanetto 2007b) over the range 1 K ≤ T K ≤ 10 4 K.
The coupling coefficient for Ly scattering is (Wouthuysen 1952;Field 1958;Madau et al. 1997 where = 1215.67 Å, Λ = 6.265 × 10 8 s −1 is the Einstein spontaneous emission coefficient for the Ly transition, is a factor of order unity that corrects for the spectral distortions in the Ly spectrum, and is the proper Ly specific intensity in units erg s −1 cm −2 Hz −1 sr −1 . We use the fits provided by Hirata (2006) to calculate and , where S , and must be solved for iteratively.
In this work we also include the effect of redshift space distortions on the 21 cm forest absorption features. In our calculation of the 21 cm optical depth, we therefore include a convolution with the Gaussian line profile and incorporate the gas peculiar velocities from our hybrid RT/hydrodynamical simulations. The optical depth in Eq. (5) may then be calculated in discrete form as (e.g. for pixel with Hubble velocity H, and velocity width 3 . Here = (2 B K / H ) 1/2 is the Doppler parameter, K is the gas kinetic temperature, = H, + pec, , and pec is the peculiar velocity of the gas. We evaluate Eq. (9) in our simulations by extracting a total of 5000 periodic lines of sight, drawn parallel to the simulation box axes at redshift intervals of Δ = 0.1 over the range 5 ≤ ≤ 12. The total path length we use to make our mock 21 cm forest spectra at each output redshift is therefore 200ℎ −1 cGpc.
The redshift evolution of the transmission, = − 21 , for a random selection of 21 cm forest spectra drawn from the zr53 simulation is shown in Fig. 4, for three different X-ray efficiencies. No instrumental features have been added to the simulated data. The detailed small-scale structure of the 21 cm absorption is displayed in the insets. One can see the strong effect that X-ray heating has on the the 21 cm absorption as the X-ray efficiency parameter is increased from X = 0.01 in the top panel, to X = 1 in the bottom panel (cf. Xu et al. 2011;Mack & Wyithe 2012). The redshift evolution Figure 5. The optical depth weighted temperature-density plane for gas in the zr53 simulation at redshift = 9 (top), 7.5 (middle) and 6 (bottom), for an X-ray efficiency X = 0.1 and Ly efficiency = 1. The colour scale shows the average 21 cm optical depth at each point in the plane. The number density of points increases by 1 dex for each contour level. due to the increasing filling factor of warm ( K ∼ 10 4 K), photoionised gas is also apparent. In particular, the occurrence of gaps in the 21 cm forest absorption due to extended regions of ionised gas increases toward lower redshift.
In order to better identify the gas associated with the absorption, we calculate the optical depth weighted density, Δ w = 1 + w , and optical depth weighted kinetic temperature, K,w , for each pixel in our zr53 mock spectra for X = 0.1. This is analogous to the approach used to study the properties of gas responsible for absorption in the Ly forest (Schaye et al. 1999); peculiar motions (and to a much lesser extent, line broadening) would otherwise distort the mapping between 21 cm optical depth, temperature and gas density. The results are shown in Fig. 5, where the temperature-density plane is displayed for the zr53 simulation at three different redshifts: = 9 (top), 7.5 (middle) and 6 (bottom). The colour bar and con-tours show the average 21 cm optical depth and the relative number density of the pixels, respectively. The gas distribution in Fig. 5 is bimodal, with the bulk of the pixels associated with either warm ( K ∼ 10 4 K), photo-ionized gas or cold ( K ≤ 10 2 K), significantly neutral regions (see also Ciardi et al. 2013;Semelin 2016). The plume of gas at intermediate temperatures has been heated by shocks from structure formation. Note, furthermore, that in this very late reionisation model the IGM is still not fully ionised by = 6. The largest optical depths in the model arise not from the highest density gas, but the cold, diffuse IGM with 3 < Δ < 10. This is because gas at higher densities is typically reionised early due to proximity to the ionising sources, and also because gas around haloes (with Δ 100) is shock-heated and partially collisionally ionised. Note again, however, there is no cold, star forming gas in this simulation -for further discussion of this point see Appendix A. Toward lower redshift, the increase in the minimum kinetic temperature of the neutral gas due to X-ray heating, the partial ionisation of the H I by secondary electrons and collisions, and the decrease in the proper number density of gas at fixed overdensity, all conspire to lower the maximum optical depth. The contours furthermore show that the regions with the largest optical depths are at least 100 times rarer than the bulk of the cold, neutral gas. Nevertheless, in this very late reionisation model, it remains possible that some detectable 21 cm absorption may persist as late as 6. We now explore this possibility in more detail.

The volume averaged 21 cm optical depth
We first consider the redshift evolution of the volume averaged 21 cm optical depth, 21 , in the zr53 simulation, displayed as the solid curves in Fig. 6 for our fiducial model with = 1. In the upper panel, we test the common assumption that, as a result of the Wouthuysen-Field effect, the spin temperature becomes strongly coupled to the gas kinetic temperature during the later stages of reionisation, such that S = K (e.g. Xu et al. 2009;Mack & Wyithe 2012;Ciardi et al. 2013). This is shown by the dashed curves in the upper panel of Fig. 6. As also noted by Semelin (2016), a full calculation of S using Eq. (6) can either reduce or enhance 21 cm optical depths relative to the value obtained assuming strong coupling. This is caused by a partial coupling of the spin temperature to the CMB temperature; if K < CMB , the full calculation will result in a higher spin temperature and smaller 21 cm optical depth, and vice versa.
This can be observed in Fig. 6 for X = 0.1 (fuchsia curves), where 21 for the full calculation assuming = 1 (solid curves) is smaller than the S = K case (dashed curves) at 8, but is greater at lower redshifts. This coincides with the temperature evolution shown in Fig. 3, particularly the transition from K < CMB (and S > K ) at > 8 to K > CMB (and < ) at < 8. Similarly, in the case of a weaker ( X = 0.01, blue curves) or stronger ( X = 1, orange curves) X-ray background, the full S calculation respectively decreases or increases 21 relative the the strong coupling approximation. The dotted curves furthermore show the 21 redshift evolution for significantly weaker Ly coupling, with = 0.01. In this case S is now decoupled from K and has a value similar to CMB . The weak coupling means 21 is significantly increased in the models with efficient X-ray heating. Hence, while the assumption of strong coupling, S = K , remains a reasonable approximation if = 1, this will not be the case if the background Ly emissivity is significantly overestimated in our fiducial model (i.e. 1). The lower panel of Fig. 6 instead shows 21 for the two different reionisation histories in Fig. 1. Both of these reionisation models are broadly consistent with existing constraints on the timing of reionisation, and the zr53 model furthermore successfully reproduces the large fluctuations in the Ly forest opacity at = 5.5 . For comparison, we also show 21 from Ciardi et al. (2013) as the dotted curve. This includes X-ray and Ly heating following Ciardi et al. (2010), and is most similar to our zr67 simulation with X 1. The differences between this work and Ciardi et al. (2013) are due to different assumptions for the X-ray emissivity and the reionisation history. A later end to reioni-sation means 21 in Fig. 6 remains significantly larger than earlier reionisation models at redshifts 6 7. If reionisation does indeed complete late, such that large neutral islands persist in the IGM at 6 (e.g. Lidz et al. 2007;Mesinger 2010), this suggests 21 cm forest absorption lines may be more readily observable than previously thought at these redshifts.

The differential number density of 21 cm absorption lines
We now consider the number density of individual absorption lines in our high resolution mock spectra. We present this as the total number of lines, , within a given optical depth bin, per unit redshift (see also Furlanetto 2006a;Shimabukuro et al. 2014), where The absorption lines in our simulated 21 cm forest spectra are identified following a similar method to Garzilli et al. (2015), who identify absorption lines in mock Ly forest spectra as local optical depth maxima located between two minima. In this work, we require that the local maxima must have a prominence (i.e. be higher by a certain value than the minima) that corresponds to a factor of 1.001 difference in the transmitted flux, = − 21 , between the line base and peak. We then define the optical depth for each identified line as being equal to the local maximum. We find this method is robust for lines with 21 ≥ 10 −2 (i.e. = − 21 ≤ 0.99), but for optical depths below this threshold the number of lines is sensitive to the choice for the prominence, and is thus unreliable. The number density distributions, 21 ( 21 , ), for different model parameters at three different redshifts, = 9, 7.5 and 6, are displayed for our fiducial model with = 1 in Fig. 7 (for an illustration of the effect of these model parameter variations on individual absorbers, see Appendix C). Each column corresponds to a different model parameter choice, each row shows a different redshift, and in each panel we show the distribution for three X-ray efficiencies: X = 0.01 (blue curves), X = 0.1 (fuchsia curves) and X = 1 (orange curves). The peak of the distribution is at 21 ≤ 0.1, and it shifts to lower amplitudes and smaller optical depths as the IGM reionises and the spin temperature of the Xray heated gas increases. The distribution also has an extended tail toward higher optical depths. While strong 21 cm absorbers will be rare, this suggests that for X ∼ 0.1, features with a transmission of = − 21 0.9 should still be present at = 7.5 in the late reionisation model (see also Fig. 4).
In the first column of Fig. 7 we re-examine the effect of strong Ly coupling on the distribution of 21 cm optical depths. As was the case for the volume averaged optical depth in Fig. 6, the impact is relatively modest for low X-ray efficiencies: for X = 0.01 at = 6, the two cases are almost identical. For X = 1, however, the abundance of features with 21 ≥ 0.01 for S = K is more than 50 per cent smaller than the full calculation at = 9. In either case, however, by = 7.5 most gas in the X = 1 model has 21 < 10 −2 , and will therefore be challenging to detect directly. However, the dotted curves also demonstrate that if = 0.01, the weak coupling of S to K allows strong 21 cm absorbers to still be observable at = 6, even for X = 1. We consider the effect of gas peculiar velocities on the 21 cm forest in the second column of Fig. 7. Redshift space distortions are well known to impact on the observability of the high redshift 21 cm signal (Bharadwaj & Ali 2004;Mao et al. 2012;Majumdar et al. 2020). We do this by creating mock 21 cm spectra that ignore the effect of gas peculiar motions, such that pec = 0 in Eq. (9). 7. The differential number density of absorption lines in synthetic 21 cm forest spectra. Each row shows the distribution at redshift = 9 (top row), 7.5 (middle row) and 6 (bottom row) for our fiducial model with Ly efficiency = 1, and in each panel the distribution is shown for three X-ray efficiencies, X = 0.01 (blue curves), X = 0.1 (fuchsia curves) and X = 1 (orange curves). Each column displays the zr53 simulation (solid curves) compared to models where one of the parameter choices is varied (dashed curves). These parameters are, from left to right, the assumption of strong Ly coupling (i.e. S = K ), neglecting the effect of peculiar velocities (i.e pec = 0), pressure smoothing due to a uniform rather than patchy UV photo-heating rate (i.e the zr53-homog model) and an earlier end to reionisation (the zr67 model). In the first column, we also show the number density distribution for very weak Ly coupling (i.e. = 0.01, dotted curves). The black dotted curves in the last column show the case of no reionisation or X-ray heating (i.e. S = K = ad = 2.73 K (1 + ) 2/3 (1 + z) 2 /(1 + z dec ), where dec = 147.8 , and HI = 1).
The results are shown by the dashed curves. While the position of the peak in the number density distribution is unchanged, the high optical depth tail is strongly affected, particularly for inefficient Xray heating. Ignoring peculiar velocities within 21 cm forest models can therefore significantly reduce the incidence of the strongest absorbers, and this will have a negative impact on the predicted observability of the 21 cm forest. Qualitatively, this agrees with the assessment of Semelin (2016), who also included the effect of gas peculiar motions in their models.
As our hybrid simulations self-consistently model the hydrodynamical response of gas to photo-heating by the inhomogeneous UV radiation field, we may also estimate the effect of (the lack of) pressure (Jeans) smoothing on the 21 cm forest. Inhomogeneous reionisation introduces large scale gas temperature fluctuations in the IGM (Keating et al. 2018), and these lead to differences in the local gas pressure that smooth the structure of the IGM on different scales (e.g. Gnedin & Hui 1998;Kulkarni et al. 2015;Nasir et al. 2016;D'Aloisio et al. 2020). In the absence of significant Xray heating, the neutral gas responsible for the 21 cm forest should therefore experience minimal pressure smoothing compared to the photo-ionised IGM. We therefore compare the results of our zr53 model to the zr53-homog simulation in the third column of Fig. 7. The latter model has exactly the same initial conditions and volume averaged reionisation history as zr53, but all the gas in the simulation volume is instead heated simultaneously (i.e. we do not follow the radiative transfer for UV photons).
The dashed curves in the third column of Fig. 7 show the line density distribution obtained from the density and peculiar velocity fields in the zr53-homog model (differences due to HI , K and S in the two models have been removed). We observe that there is a small, redshift dependent difference between the two distributions, such that the simulation with the homogeneous UV background exhibits fewer strong absorption lines. This is because the gas responsible for the highest optical depths in the 21 cm forest (see Fig. 5) is still cold within the hybrid model, and hence has slightly higher density due to the smaller pressure smoothing scale.
We caution, however, that this comparison will still not fully capture the effect of pressure smoothing on 21 cm forest absorbers.
where J is the Jeans scale, is the mean molecular weight of hydrogen and helium assuming primordial composition ( = 1.22 for fully neutral gas, = 0.59 for fully ionised), and J = p / J is a factor of order unity that accounts for the finite time required for gas to dynamically respond to a change in pressure. For comparison, the mean interparticle separation and gravitational softening length in our simulations are 19.5ℎ −1 ckpc and 0.78ℎ −1 ckpc, respectively. Eq. (11) thus implies that the pressure smoothing scale for typical 21 cm forest absorbers is not fully resolved in our simulations (see also Emberson et al. 2013). We furthermore do not capture the 21 cm absorption from minihaloes with < 2.5 × 10 7 (Furlanetto 2006a). Larger differences could then be observed in Fig. 7 for fully resolved gas. On the other hand, although we follow the dynamical response of gas to heating by UV photons, the X-ray heating of the neutral gas in our hybrid simulation is applied in post-processing. It is therefore decoupled from the hydrodynamics, and this may then underestimate the impact of pressure smoothing on cold gas for high X-ray efficiencies. Regardless of these modelling uncertainties, however, this suggests that the effect of the pressure smoothing scale on the 21 cm forest in the diffuse IGM remains small compared to the substantial impact of X-ray heating on the spin temperature at ≤ 10. Finally, in the fourth column of Fig. 7 the effect of the reionisation history is displayed for the zr53 (solid curves) and zr67 (dashed curves) simulations. For comparison, the dotted curves also show the line number density distribution under the assumption of no reionisation or X-ray heating (i.e. S = K = ad and HI = 1). As expected, the two reionisation models are significantly different at = 6; there are no strong absorption features with 21 > 10 −2 in zr67 model, as reionisation has already completed by this time. At = 7.5 one can see that there are also fewer absorption features in the zr67 model due to the larger volume of ionised gas. However, the differences between the two models become smaller with increasing redshift. This again demonstrates that for reionisation models that complete at < 6, the 21 cm forest may remain observable if sufficiently bright radio sources exist at 6 < < 7. Alternatively, a null-detection could place an interesting limit on the very uncertain X-ray background (e.g. Mack & Wyithe 2012). We now investigate this possibility further.

Detectability of strong 21 cm forest absorbers at redshift
= 6 for late reionisation and X-ray heating A detection of the 21 cm forest relies on the identification of objects at high redshift that are sufficiently radio bright to act as background sources. Based on a model for the radio galaxy luminosity function at > 6, Saxena et al. (2017) predict around one radio source per 400 square degrees at a flux density limit of 150 MHz = 3.5 mJy, and at least ∼ 30 bright sources with 150 MHz > 15 mJy (see also Bolgar et al. 2018). Ongoing observational programmes such as the LOFAR Two-metre Sky Survey (LoTSS, Shimwell et al. 2017;Kondapally et al. 2020), the Giant Metrewave Radio Telescope (GMRT) all sky radio survey at 150 MHz (Intema et al. 2017), and the Galactic and Extragalactic All-sky Murchison Widefield Array survey (GLEAM, Wayth et al. 2015) should furthermore detect hundreds of bright > 6 radio sources. Encouragingly, a small number of radio-loud sources have already been identified at > 5.5 (e.g. Bañados et al. 2018b), including the = 6.1 blazar PSO J0309+27 with a flux density 147 MHz = 64.2 ± 6.2 mJy (Belladitta et al. 2020) We now use our hydrodynamical simulations to assess the feasibility of detecting the 21 cm forest in late reionisation models, assuming = 1. We shall calculate the minimum redshift path length, Δ min , necessary for detecting a single, strong (i.e. 21 > 0.01) absorption line with a minimum transmission at some arbitrary threshold th = − 21, th . For a signal-to-noise ratio S/N, the minimum flux density contrast, Δ min , detectable by an interferometric radio array is then (e.g. Ciardi et al. 2015a), where sys is the system temperature, Δ is the bandwidth, eff is the effective area of the telescope, int is the integration time, and min is the minimum intrinsic flux density a radio source must have to allow detection of a 21 cm absorption feature with a minimum at a flux density of abs = min In what follows, we shall adopt values for the sensitivity, eff / sys , in Eq. (12) appropriate for LOFAR, SKA1-low and SKA2, where eff / sys 80 m 2 K −1 , 600 m 2 K −1 and 5500 m 2 K −1 , respectively 4 (Braun et al. 2019). Additionally, to approximately model the effect of spectral resolution on the data we convolve our mock spectra with a boxcar function. Following the bandwidths adopted in Ciardi et al. (2015b), we assume boxcar widths of 10 kHz and 5 kHz for LOFAR and SKA1-low, respectively. For a more futuristic measurement with SKA2, we assume a smaller bandwidth and adopt a boxcar width of 1 kHz.
First, in Fig. 8, we show the minimum redshift path length Δ min required to detect a single 21 cm absorption line in the minimum transmission threshold th -redshift plane for three different X-ray efficiencies X (upper panels), or in the X -redshift plane for three different transmission thresholds th (lower panels). Note that for now we assume a sufficient number of background radio sources exists for such a measurement; we consider the issue of detectability at = 6 further in Fig. 9. The mock spectra used in Fig. 8 are drawn from the zr53 simulation and have been convolved with a boxcar of width 5 kHz (i.e. our assumed SKA1-low bandwidth). Unshaded white regions indicate where no absorbers are present over our total simulated path length of 200ℎ −1 cGpc. Fig. 8 shows that no absorption features with th ≤ 0.77 should be present at 8 for even a very low X-ray efficiency of X = 0.01 in the late reionisation model. Similarly, almost no strong 21 cm absorption with th 0.99 will exist at < 7 for X ≥ 1. This highlights the challenging nature of 21 cm forest measurements from the diffuse IGM, even if reionisation ends very late, and also how sensitive the 21 cm forest absorption is to X-ray heating. Proposals to use the 21 cm forest as a sensitive probe for distinguishing between different cosmological or dark matter models using the diffuse IGM are therefore likely to be restricted to very high redshifts, prior to any substantial X-ray heating of the IGM.
As a reference, the black curves in Fig. 8 correspond to the redshift path length obtainable by a hypothetical observation of 1, 4 Note that in reality the sensitivity eff / sys is frequency dependent. However, over the frequency range we consider, 142 MHz ≤ 21 /(1 + z) ≤ 203 MHz, this dependence is reasonably weak. See fig. 8 in Braun et al. (2019) for further details. Figure 8. The minimum redshift path length, Δ min , required to observe a single 21 cm absorption feature in the zr53 simulation assuming that a sufficient number of background radio sources exist. The mock spectra have been convolved with a boxcar of width 5 kHz to approximately model the effect of spectral resolution on the lines. In the upper panels we show Δ min in the th -z plane for an X-ray efficiency factor of X = 1 (left), 0.1 (middle) and 0.01 (right). In the lower panels we instead show Δ min in the X -z plane for a 21 cm absorption feature with minimum transmission ≤ th = 0.8 (left), 0.9 (middle) and 0.99 (right). Here we also note the minimum intrinsic flux density, min , that a background radio source must have such that an absorption line with a minimum at ≤ th is detectable with SKA1-low at a signal-to-noise of S/N = 5 and integration time of int = 1000 hr (see Eq. 12). The unshaded white regions are where no absorbers are present over our total simulated path length of 200ℎ −1 cGpc. The thick black curves in each panel track the redshift path length that would be covered by the observation of 1 (dotted), 10 (dashed) and 100 (solid) radio sources assuming redshift bins of width Δ = 0.2.
10 or 100 radio sources of sufficient brightness in redshift bins of width Δ = 0.2 (i.e. an observation of radio sources provides a total redshift path length of 0.2 ). 5 A null-detection over this path length would provide a model dependent lower limit on the X-ray background emissivity, such that X ≥ X,max , where X,max is the maximum X-ray efficiency that retains at least one strong absorption feature with ≤ th . From the lower middle panel in Fig. 8, the null-detection of a feature with th < 0.9 at = 9 in 1 (10) radio source(s) implies X,max 0.04 ( X.max 0.07). The parameter space that lies below the black curves would then be disfavoured. In practice, however, radio telescope sensitivity, spectral resolution and the availability of sufficiently bright background radio sources will impact upon the detectability of strong lines. We quantify this in Fig. 9, where similarly to Fig. 8 we show Δ min , but now in the Xth plane at redshift = 6. This is shown for our LOFAR (left), SKA1-low (middle) and SKA2 (right) model assumptions, where we have convolved the synthetic spectra with a 5 The choice of Δ = 0.2 is somewhat arbitrary -we require a bin that is small enough that redshift evolution is not significant, but large enough to probe a reasonable path length. For reference, increasing the bin size to Δ = 0.4 would approximately halve the number of background sources required to detect a single absorber with th , assuming minimal redshift evolution across the bin. Table 2. The maximum X-ray background efficiency, X,max , that retains at least one strong 21 cm absorption feature with transmission ≤ th in our synthetic 21 cm forest spectra, for a redshift path length corresponding to bright radio sources covering a redshift bin of width Δ = 0.2, centred at redshift = 6. The mock spectra have been convolved with a boxcar of width 5kHz to approximately model the effect of observed bandwidth on the lines. The minimum intrinsic flux density of the background source, min , required to detect a line with th at a signal-to-noise of S/N = 5 with SKA1-low is calculated using Eq. (12), assuming a bandwidth Δ = 5 kHz, sensitivity eff / sys = 600 m 2 K −1 and integration time of int = 1000 hr. The expected number of radio sources in the sky with min at = 6, N S17 , are estimated from Saxena et al. (2017) (their fig. 11). In the event of a null-detection of an absorption feature with th , the X,max values give a (model dependent) lower limit on the X-ray efficiency.  boxcar of width 10 kHz, 5 kHz and 1 kHz, respectively. The minimum intrinsic source flux density, min , required to detected a line with th has also been calculated using Eq. (12) and is displayed Figure 9. As for Fig. 8, but now the minimum redshift path length Δ min is shown in the Xth plane at redshift = 6. The mock spectra have been convolved with a boxcar of width 10 kHz (left), 5 kHz (middle) and 1 kHz (right) to approximately model the effect of our assumed bandwidths for LOFAR, SKA1-low and SKA2, respectively. Note the scale on the horizontal axis is different in each panel. The upper horizontal axis now also shows the minimum intrinsic flux density, min , required for a background radio source, such that a line with minimum transmission ≤ th is detectable at a signal-to-noise of S/N = 5 by LOFAR with an integration time of int = 1000 hr (left), by SKA1-low with int = 1000 hr (middle) and by SKA2 with int = 100 hr (right). In the left panel (LOFAR) the thick dashed curve is truncated where the number of available background radio sources predicted by Saxena et al. (2017) with min at 6 falls below 10. The curve for 100 sources (solid) is not shown, as this exceeds the expected radio source number from Saxena et al. (2017) at the required min . Table 3. As for Table 2, except the mock 21 cm forest spectra are now smoothed with a boxcar of width 10 kHz and the minimum source flux densities, min , have been computed for LOFAR using a bandwidth Δ = 10 kHz, sensitivity eff / sys = 80 m 2 K −1 and integration time of int = 1000 hr. Dashes mean that X,max is not measurable due to the lack of expected sources.  on the horizontal top axis. Here we assume a strong absorption line with minimum transmission th is detected with S/N = 5 for an integration time of int = 1000 hrs with LOFAR and SKA1-low, and int = 100 hrs with SKA2. First, one can see that if using a more sensitive telescope with higher spectral resolution it is possible to detect deeper, narrower absorption features. Moreover, tighter constraints on the X-ray efficiency X may also be obtained. For example, at = 6, there are no absorption features with ≤ 0.95 for X > 0.01 if observed by LOFAR. However, this increases to X > 0.025 for SKA1-low and X > 0.05 for SKA2. The minimum source flux density required to detect an absorption feature at fixed S/N = 5 also decreases significantly, thus increasing the number of potentially suitable background radio sources.
We quantify this in more detail in Tables 2, 3 and 4, where we list the maximum X-ray efficiency, X,max , that retains at least one 21 cm absorption feature at = 6 with a transmission minimum ≤ th over a path length of Δ = 0.2, Δ = 2 or Δ = 20 in the zr53 simulation. This corresponds to = 1, 10 and 100 sources, respectively, for redshift bins of width Δ = 0.2. We also give the minimum flux density, min , required to detect an absorption line with th at S/N = 5. Additionally, we give the expected number of background sources in the sky at 6 with min reported by Saxena et al. (2017) for an observing time of 100 hrs with the standard LOFAR configuration (see their fig. 11). As a quantitative example, using Table 2 (SKA1-low), for a ten background sources with 203 MHz = 3.4 mJy, on average we would expect to detect at least one 21 cm absorption line with < 0.95 at = 6.0 ± 0.1 if X ≤ 0.007. A null-detection would instead imply a lower limit of X > 0.007. Within our model, this X-ray efficiency may be converted to an estimate of the soft X-ray band emissivity at 0.5-2 keV, where X,0.5−2 keV = 10 38.3 X erg s −1 cMpc −3 at = 6. Hence X > 0.007 corresponds to X,0.5−2 keV > 10 36.1 erg s −1 cMpc −3 . Alternatively, from Table 3 (LOFAR), the null-detection of an absorption line with < 0.95 at = 6.0 ± 0.1 in the spectra of 10 radio bright sources with 203 MHz = 18.2 mJy would imply a slightly weaker constraint of X > 0.001 and X,0.5−2 keV > 10 35.3 erg s −1 cMpc −3 . This suggests that lower limits on the soft X-ray background emissivity at high redshift from a null-detection of the 21 cm forest may complement existing constraints from upper limits on the 21 cm power spectrum (Greig et al. 2021a). We note, however, these results are highly model dependent. If the Ly coupling is very weak (i.e. if 1), or there is a significant contribution to the 21 cm forest absorption from unresolved small scale structure, the X,max values in Table 2-4 will translate to lower limits on X that are conservative.

CONCLUSIONS
We have used very high resolution hydrodynamical simulations combined with a novel approach for modelling patchy reionisation to model the 21 cm forest during the epoch of reionisation. Our simulations have been performed as part of the Sherwood-Relics simulation programme (Puchwein et al. in prep). In particular, we have considered the observability of strong ( 21 > 10 −2 ) 21 cm absorbers in a late reionisation model consistent with the large Ly forest transmission fluctuations observed at = 5.5 , where large neutral islands of intergalactic gas persist until 6 Keating et al. 2020). We also explore a wide range of assumptions for X-ray heating in the pre-reionisation intergalactic medium (IGM), and have assessed the importance of several common modelling assumptions for the predicted incidence of strong 21 cm absorbers. Our key results are summarised as follows: • In a model of late reionisation ending at = 5.3, for an X-ray efficiency parameter X 0.1 (i.e. for relatively modest X-ray pre-heating of neutral hydrogen gas, such that the gas kinetic temperature K 10 2 K) strong 21 cm absorption lines with optical depths 21 ≥ 0.01 situated in neutral islands of intergalactic gas should persist until = 6. In this case, the 21 cm absorbers with the largest optical depths should arise from cold, diffuse gas with overdensities 3 < Δ < 10 and kinetic temperatures K < 10 2 K. A null-detection of 21 cm forest absorbers at = 6 may therefore place a valuable lower limit on the high redshift soft X-ray background and/or the kinetic temperature of the diffuse pre-reionisation IGM in the neutral islands. With ∼ 10 radio-loud active galactic nuclei now known at 5.5 < < 6.5 (e.g. Bañados et al. 2018b;Liu et al. 2021) and the prospect of more radio-loud sources being identified in the next few years, this possibility merits further investigation.
• By far the largest uncertainty in models of the 21 cm forest is the heating of the pre-reionisation IGM by the soft X-ray background (see also Mack & Wyithe 2012). In the absence of strong constraints on the soft X-ray background at ≥ 6, proposals to use the 21 cm forest to distinguish between cosmological models (where differences between competing models are small compared to the effect of X-ray heating) will likely be restricted to redshifts prior to the build-up of the soft X-ray background. Uncertainties in the strength of the Wouthuysen-Field coupling will also be important to consider if the Ly background is significantly weaker than expected from extrapolating the observed star formation rate density to > 6. By contrast, we find the effect of uncertain pressure/Jeans smoothing on the 21 cm absorption from the diffuse IGM should remain comparatively small.
• Models of the 21 cm forest must include the effect of gas peculiar motions on absorption line formation to accurately predict the incidence of strong absorption features (see also Semelin 2016). Ignoring redshift space distortions reduces the incidence of the strongest 21 cm forest absorbers, and results in a maximum optical depth in the 21 cm forest that is up to a factor of ∼ 10 smaller compared to a model that correctly incorporates gas peculiar velocities.
• We present model dependent estimates for the minimum redshift path length required to detect a single, strong 21 cm forest absorption feature as a function of redshift and X-ray efficiency parameter, X within a late reionisation model that ends at redshift = 5.3. At = 6.0 ± 0.1 for an integration time of int = 1000 hrs per background radio source, a null-detection of 21 cm forest absorbers with < 0.95 at a signal-to-noise of S/N = 5 in the spectra of 10 radio sources with 203 MHz > 3.4 mJy (> 18.2 mJy) using SKA1-low (LOFAR) implies a soft X-ray background emissivity X,0.5−2 keV > 10 36.1(35.3) erg s −1 cMpc −3 . As the soft X-ray background at high redshift is still largely unconstrained, this suggests lower limits on the X-ray emissivity from a null-detection of the 21 cm forest could provide a valuable alternative constraint that complements existing and forthcoming constraints from upper limits on the 21 cm power spectrum.
While the calculation we present in this work is illustrative, a more careful forward modelling of the 21 cm absorption data is still required. We have not considered how to recover absorption features from noisy data beyond the simple signal-to-noise calculation adopted here, or how an imperfect knowledge of the radio source continuum and/or radio background might impact upon the detectability of 21 cm absorbers. Uncertainties in other parameters such as the reionisation history and the Ly background emissivity should furthermore be marginalised over to obtain a robust lower limit on the soft X-ray background. Our simulations do not account for the absorption from unresolved mini-haloes with masses < 2.5 × 10 7 , and will lack coherent regions of neutral gas on scales greater than our box size of 40ℎ −1 cMpc. On the other hand, even a modest amount of feedback, either in the form of photoevaporation (Park et al. 2016;Nakatani et al. 2020) or feedback from star formation (Meiksin 2011) will substantially reduce the absorption signature from minihaloes. These feedback effects may be particularly important during the final stages of reionisation at 6, where any remaining 21 cm absorption should arise from neutral islands in the diffuse IGM.
More detailed models the 21 cm forest will require either radiation-hydrodynamical simulations that encompass a formidable dynamic range, and/or multi-scale, hybrid approaches that adopt sub-grid models for unresolved absorbers and their response to feedback. Both must furthermore cover a very large and uncertain parameter space. Nevertheless, we conclude that if reionisation completes at < 6, the prospects for using SKA1-low or possibly LOFAR to place an independent constraint on the soft X-ray background using strong absorbers in the 21 cm forest are encouraging. managed by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility. The equipment was funded by BEIS capital funding via STFC capital grants ST/P002293/1 and ST/R002371/1, Durham University and STFC operations grant ST/R000832/1. DiRAC is part of the National e-Infrastructure. This work has made use of matplotlib (Hunter 2007), astropy (Astropy Collaboration et al. 2013), numpy (Harris et al. 2020) and scipy (Virtanen et al. 2020). TŠ is supported by the University of Nottingham Vice Chancellor's Scholarship for Research Excellence (EU). JSB acknowledges the support of a Royal Society University Research Fellowship. JSB and NH are also supported by STFC consolidated grant ST/T000171/1. MGH acknowledges support from the UKRI STFC (grant numbers ST/N000927/1 and ST/S000623/1). We thank Volker Springel for making P-GADGET-3 available.

DATA AVAILABILITY
All data and analysis code used in this work are available from the first author on reasonable request. An open access preprint of the manuscript will be made available at arXiv.org.

APPENDIX A: TEST OF THE PRESCRIPTION FOR CONVERTING DENSE GAS INTO COLLISIONLESS PARTICLES
As discussed in Section 2.1, we adopt a simplified scheme for the treatment of dense, star forming gas in the Sherwood-Relics simulations, where all gas particles with density Δ > 1000 and temperature K < 10 5 K are converted to collisionless star particles (Viel et al. 2004). As a consequence, very dense, cold halo gas is not included in the Sherwood-Relics models. We test whether this approximation affects our results for the 21 cm forest in Fig. A1. Here we compare two models drawn from the Sherwood simulation suite (Bolton et al. 2017) that use the same box size, mass resolution and initial conditions as the other simulations used in this work. These two additional simulations use either the simplified scheme used in this study (QLy ) or the star formation and energy driven winds prescription of Puchwein & Springel (2013, PS13). The only difference between these two models is the incorporation of dense, star forming gas within the PS13 simulation.
In Fig. A1, we show the differential line number distribution obtained after applying the neutral fraction, gas kinetic and spin temperature from the patchy zr53 simulation to the native density and peculiar velocity fields from the QLy and PS13 models. As before, we consider three different X-ray efficiencies. We observe little to no difference in the statistics of the 21 cm forest computed from these two simulations. This is because the highest density gas is usually located close to ionising sources, and so is often too hot, ionised or rare to exhibit significant amounts of strong absorption in the hyperfine line. This is further illustrated by the cyan curves in Fig. A1, where the mock 21 cm forest spectra are instead computed assuming a fully neutral, isothermal gas with S = K = ad,0 , where ad,0 = 2.73 K(1 + ) 2 /148.8 is the gas temperature assuming adiabatic cooling at the mean density. Small differences due to the presence of the high density gas in the PS13 simulation are now apparent in the tail of the distribution at 21 10. However, if we also include the adiabatic heating of the gas by compression, such that ad = ad,0 (1 + ) 2/3 , these models become almost identical. We conclude that the approximate treatment of dense, star forming gas we adopt in this work should not significantly change our key results. The relative rarity of 21 cm absorption from cold gas within massive haloes suggests this population will in any case be completely dominated by 21 cm absorbers from the diffuse IGM and/or minihaloes during reionisation.

APPENDIX B: CALCULATION OF THE X-RAY AND Ly SPECIFIC INTENSITIES, GAS KINETIC TEMPERATURE AND IONISATION STATE
In this section we describe the model for the X-ray heated IGM introduced in Section 2.2. The X-ray background is primarily responsible Figure A1. The differential line number density distribution for 21 cm forest absorption features in simulations that use two different implementations for the treatment of dense gas. The dashed curves displays the simplified approach used in this work (the QLy simulation), whereas the solid curves use the sub-grid star formation and feedback model from Puchwein & Springel (2013) (the PS13 simulation, solid line). The results are shown at three different redshifts: = 9 (top), 7.5 (middle) and 6 (bottom). The X-ray efficiencies are X = 0.01 (blue curves), 0.1 (fuchsia curves) and 1 (orange curves). For comparison, the cyan curves show the distribution for fully neutral, unheated gas with temperature equal to the adiabatic temperature at mean density (i.e. HI = 1 and S = K = ad,0 ). for ionising and heating the intergalactic medium prior to reionisation. The proper specific intensity of the X-ray background, X, [erg s −1 cm −2 Hz −1 sr −1 ], is given by the solution to the cosmological radiative transfer equation (Haardt & Madau 1996 where , is the comoving X-ray emissivity, ★ is the redshift when X-ray emitting sources first form, and the emission frequency, of a photon emitted at redshift and observed at frequency and redshift is = (1 + ) (1 + ) .
The optical depth encountered by a photon observed at frequency is where the sum is over the species = H I, He I, He II, and , are the photo-ionisation cross-sections .
The photo-ionisation rates for species = H I, He I, He II are where is the frequency of the ionisation threshold for species . The second term in Eq. (B4) arises from secondary ionisations due to collisions with energetic photo-electrons, where Φ i is the number of secondary ionisations per primary photo-electron of energy ℎ p ( − i ) for a free electron fraction of e (Shull & van Steenberg 1985). The corresponding photo-heating rates are where heat is the fraction of the primary photo-electron energy that contributes to the heating of the gas. We use the tabulated results from Furlanetto & Stoever (2010) for Φ i and heat . The Compton scattering of free electrons off X-ray background photons will also heat the IGM (Madau & Efstathiou 1999). The Compton heating rate is where T = 6.65 × 10 −25 cm 2 is the Thomson cross-section, appropriate for X-ray photons with energy < ∼ 100 keV (i.e. relativistic effects may be ignored).
The Ly background has two contributions: emission from stars, and Ly photons produced by the excitation of H I atoms by X-ray photons. The proper Ly specific intensity from stars, ,★ , requires consideration of both Ly and higher order Lyman series photons. This is because the Ly photons redshift into resonance at redshift and generate Ly photons via a series of radiative cascades to lower energies , such that where n is the probability of producing a Ly photon from a Figure C1. Left: An example line of sight drawn from our mock 21 cm forest spectra at ∼ 7.5 for an X-ray efficiency of X = 0.1 and Ly efficiency = 1. From top to bottom, we show the zr53 model with the solid black curve, and compare this to several model parameter variations (orange dotted curves): gas peculiar velocities set to zero (top), pressure smoothing under the assumption of homogeneous heating in the zr53-homog simulation (middle) and an earlier end to reionisation in the zr67 model (bottom). Right: The quantities responsible for the observed differences between the spectra displayed in the left column. From top to bottom, these are the gas peculiar velocity, pec , normalised gas density Δ = (1 + ) = / , and hydrogen neutral fraction, HI . In the middle panels we also present a zoomed-in view of an absorption feature (left) and the associated density peak (right), which has been broadened by pressure smoothing in the simulation with homogeneous heating.