# The forest at EndEoR: The effect of Lyman Limit Systems on the End of
Reionisation

Ivelin Georgiev 1 ID , Garrelt Mellema1 ID , and Sambit K. Giri2 ID
1Department of Astronomy and Oskar Klein Centre, AlbaNova, Stockholm
University, SE-10691 Stockholm, Sweden
2Nordita, KTH Royal Institute of Technology and Stockholm University, Hannes
Alfvéns väg 12, SE-106 91 Stockholm, Sweden
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The final stages of cosmic reionisation (EndEoR) are expected to be strongly
regulated by the residual neutral hydrogen in the already ionised regions of
the Universe. Its presence limits the mean distance that ionising photons can
travel and hence, the extent of the regions that sources of ionising photons
can affect. The structures containing most of this residual neutral hydrogen
are typically unresolved in large-scale simulations of reionisation. Here, we
investigate and compare a range of approaches for including the effect of
these small-scale absorbers, also known as Lyman limit systems (LLS), in such
simulations. We evaluate the impact of these different approaches on the
reionisation history, the evolution of the ultraviolet background, and its
fluctuations. We also compare to observational results on the distribution of
Lyman-$\alpha$ opacity towards the EndEoR and the measured mean free path of
ionising photons. We further consider their effect on the 21-cm power
spectrum. We find that although each of the different approaches can match
some of the observed probes of the final stages of reionisation, only the use
of a redshift-dependent and position-dependent LLS model is able to reproduce
all of them. We therefore recommend that large-scale reionisation simulations,
which aim to describe both the state of the ionised and neutral intergalactic
medium use such an approach, although the other, simpler approaches are
applicable depending on the science goal of the simulation.

###### keywords:

theory – large-scale structure of Universe – reionisation

††pubyear: 2024††pagerange: The forest at EndEoR: The effect of Lyman Limit
Systems on the End of Reionisation–The forest at EndEoR: The effect of Lyman
Limit Systems on the End of Reionisation

##  1 Introduction

The Epoch of Reionisation (EoR) was the period in the first billion years of
the Universe during which the neutral hydrogen (H i ) present in the
intergalactic medium (IGM) was ionised by ultraviolet (UV) radiation. The
ionising radiation is thought to originate from the first stars, galaxies, and
accreting black holes, which formed due to the gravitational collapse of
baryons into the dark matter halos. This phase transition spans a large range
of physical scales. In addition to the cosmological properties of the large-
scale IGM, the EoR is also defined by the astrophysical properties of galaxy
formation and evolution (see Pritchard & Loeb, 2012; Dayal & Ferrara, 2018,
for a review).

A key constraint on reionisation comes from the measurement of the Thomson
scattering optical depth ($\tau$) of Cosmic Microwave Background (CMB) photons
(Planck Collaboration et al., 2020) against free electrons produced during the
EoR. The scattering causes polarization of the CMB on scales larger than the
horizon size and suppresses temperature anisotropies on scales smaller than
the horizon size. The most recent constraint of $\tau=0.051\pm 0.006$ suggests
that the halfway point of the EoR was at $z_{\mathrm{mid}}\sim 8$ (assuming an
instantaneous reionisation model). Another key probe comes from the scattering
of Lyman-$\alpha$ (Ly-$\alpha$, rest-frame wavelength of 1216 Å), photons
passing through an ionised region, by residual neutral hydrogen along a line-
of-sight (LOS). Notable observations include the Ly$\alpha$ equivalent width
(Mason et al., 2018, 2019; Hayes & Scarlata, 2023; Hoag et al., 2019; Bolan et
al., 2022; Bruton et al., 2023; Morishita et al., 2023; Jones et al., 2023),
the Ly$\alpha$ dark pixel fraction (McGreer et al., 2015; Jin et al., 2023) as
well as the decreasing fraction of Ly$\alpha$ emitters among Lyman break
galaxies (Ouchi et al., 2010; Endsley & Stark, 2022) can be utilised to place
constraints on the global neutral fraction of hydrogen at a given redshift
(see the top panel of Fig. 1 in grey).

Particularly, spectra of Quasi-Stellar Objects (QSO or quasars) at high-$z$
have been found to contain a wealth of information about the state of the IGM.
Observations of the quasar damping wing (Gunn & Peterson, 1965) in high-$z$
quasar spectra indicate the scattering of Ly-$\alpha$ photons by neutral
hydrogen (H i) and can be used to estimate its fraction as a function of
redshift (e.g. Mortlock et al., 2011; Greig et al., 2017; Bañados et al.,
2018; Davies et al., 2018; Greig et al., 2019; Wang et al., 2020; Yang et al.,
2020; Greig et al., 2022). Therefore, QSOs are an excellent probe of the end
stages of reionisation (hereafter EndEoR) (see Davies et al., 2024, for a
comprehensive review). How the EndEoR occurred is essential for understanding
the ionising budget of the universe (Cain et al., 2023), and the role of the
last islands of neutral hydrogen (Giri et al., 2024). While observations from
Becker et al. (2001); Fan et al. (2006) initially indicated the endpoint of
reionisation to be around $z\sim 6$, the Ly$\alpha$ opacities such as those
derived from QSO spectra in Bosman et al. (2022) are difficult to reproduce
with reionisation simulations, requiring an extended reionisation down to
$z\approx 5.3$ (see, e.g. Kulkarni et al., 2019; Keating et al., 2020).

Moreover, once a region of the IGM is re-ionised, its internal state depends
on the ionised hydrogen (H ii ) recombination rate against the ionisation rate
set by the photon density or Ultraviolet Background (UVB) in the ionised IGM.
How the mean UVB within such ionised areas changed over time depends on the
growth and merging of such ionised structures. We will refer to the mean UVB
value in the later stages of the EoR and post-EoR as the UV Background.
Measurements of this quantity towards the EndEoR exist (Calverley et al.,
2011; Wyithe & Bolton, 2011; Gaikwad et al., 2023) but most simulations
struggle to reproduce this number, typically reporting UVB values higher than
the measurements (e.g. Bolton & Haehnelt, 2007; Finlator et al., 2009; Aubert
et al., 2015; Dixon et al., 2016; Shukla et al., 2016; Ocvirk et al., 2020).
The value of the UVB depends on the mean free path (MFP) of ionising (Lyman
continuum) photons, which is the mean distance an ionising photon travels
before being absorbed (e.g. Fan et al., 2024). The MFP can be set by the
absorption of ionising photons by neutral regions at the edge of the ionised
bubble or by absorptions within the ionised regions. Absorption within the H
ii bubble could occur due to residual neutral hydrogen or due to the presence
of discrete dense systems with an optical depth to ionising photons of
$\tau>1$, referred to as Lyman Limit Systems (LLS, see e.g. Alvarez & Abel,
2012; McQuinn et al., 2011; Shukla et al., 2016; D’Aloisio et al., 2020;
Theuns & Chan, 2024; Fan et al., 2024). The latter case is dominant in the
post-reionisation Universe as fully neutral islands have disappeared. The MFP
value towards the EndEoR has been measured by Becker et al. (2021); Zhu et al.
(2023); Gaikwad et al. (2023) and found to be significantly smaller compared
to the inferred value from post-EoR observations in Worseck et al. (2014)
(hereafter W14). Therefore, the evolution of the MFP of ionising photons is
particularly sensitive to the final stages of reionisation and vice versa,
while its value is influenced by the presence of LLS.

Modelling the EndEoR has not received much attention in the literature, partly
because the required simulation techniques are non-trivial but also because
the prospective probe of the 21-cm signal from the IGM will be weaker. In
reionisation simulations, the opacity of the IGM for ionising photons becomes
essential, and there exists quite a bit of uncertainty about this quantity.
For example, the excursion set approach used in 21CMFAST (Mesinger &
Furlanetto, 2007) is not applicable once most of the IGM has reionised. This
limitation of the approach had triggered the development of ISLANDFAST (Xu et
al., 2016). This code begins reionisation of the IGM using the same method as
in 21CMFAST and later switches to another formalism developed in Xu et al.
(2014) to follow the evolution of neutral hydrogen islands.

In this work, we use the state-of-the-art C2Ray simulation framework that does
not have this limitation as it numerically solves the three-dimensional
radiative transfer equation. The main objective of this work is centred on
studying the final stages of reionisation. Notably, we investigate the role of
unresolved LLS, which we will also refer to as small-scale absorbers, in
large-scale reionisation simulations. We construct and compare various
approaches to modelling the effect of these small-scale absorbers to examine
how their implementation impacts the properties of reionisation. To achieve
this, we model multiple observables of reionisation, from probing the state of
the ionised IGM in quasars to studying the neutral islands with the 21-cm
signal from the IGM during the EndEoR.

The paper is organised as follows. Sec. 2 describes the simulations used
within this study, focusing on the LLS implementations necessary for our
analysis. In Sec. 3, we examine the role of the small-scale absorbers on the
global parameters of reionisation and their effect on the statistics of the
UVB background. In Sec. 4, we examine the observables at the EndEoR, such as
the effective Ly$\alpha$ optical depth, the mean free path of ionising
photons, as well as the effect on the neutral islands at the end of
reionisation through the lens of the 21 cm power spectrum. We summarise our
findings and elaborate on future improvements in Sec. 5. Unless stated
otherwise, we address comoving megaparsec as Mpc.

##  2 Modelling

In this section, we introduce our simulation framework and present the methods
used to model the LLS.

###  2.1 Reionisation Simulations

The methodology for producing our reionisation simulations was described in
detail in Iliev et al. (2006); Mellema et al. (2006b) and Dixon et al. (2016).
Here we only provide an outline. First, the density and halo fields are
computed with the cosmological $N$-body code CUBEP3M (Harnois-Déraps et al.,
2013). Halos were found on-the-fly using the spherical overdensity halo-finder
described in Watson et al. (2013). Afterwards, these outputs are post-
processed using the $\mathrm{C}^{2}$Ray radiative transfer code (Mellema et
al., 2006a).

The sources of ionising photons are associated with the halos. For the
simulations used here, the ionising source model makes sure that the total
number of ionising photons at a given redshift is always proportional to the
fraction of matter which has collapsed into dark matter halos above the mass
threshold, $f_{\mathrm{coll}}$. To achieve this, the rate of ionising photons
escaping a halo of a given mass is calculated by multiplying the growth rate
of halos by a constant $\zeta$, an emissivity parameter. The resulting value
is then distributed over all halos above the minimal mass contributing
ionising photons by assuming a linear relationship between the ionising
luminosity and the mass of the halo (see Giri & Mellema, 2021, for more
details). This source model is described by the equation

$\dot{N}_{\gamma}=\zeta\frac{M}{m_{\mathrm{p}}}\frac{\Omega_{\mathrm{b}}}{\Omega_{\mathrm{m}}}\frac{1}{f_{\mathrm{coll}}}\frac{\mathrm{d}f_{\mathrm{coll}}}{\mathrm{d}t}\
,$ (1)

where $\dot{N}_{\gamma}$ is the emissivity of a single source in units of s-1,
$M$ is the mass of the dark matter halo, $m_{\rm p}$ is the mass of a proton,
and $\Omega_{\mathrm{m}}$ and $\Omega_{\mathrm{b}}$ are the density parameters
for dark and baryonic matter. The simulations presented in Tab. 1 have a
volume of $(L_{\mathrm{box}}\approx 349~{}\mathrm{cMpc})^{3}$ with $4000^{3}$
dark matter particles, resolving halos of mass down to $10^{9}M_{\odot}$. The
density field is gridded to a Cartesian 3D mesh with a size of
$(N_{\mathrm{grid}}=250)^{3}$ by applying a smooth-particle-hydrodynamics
kernel (Monaghan, 1992; Iliev et al., 2014). We note that unlike the
simulations described in Dixon et al. (2016), no subgrid haloes and radiative
feedback effects are used. We have adopted a flat $\Lambda$CDM cosmology with
parameters ($\Omega_{\mathrm{m}},\Omega_{\mathrm{b}},h,n_{s},\sigma_{8}$) =
(0.27, 0.044, 0.7, 0.96, 0.8) in agreement with WMAP (Hinshaw et al., 2013)
and Planck results (Planck Collaboration et al., 2014).

We can express the relationship between the ionisation rate and the global
emissivity ($\dot{N}_{\mathrm{glob}}$, given in units of s-1 cm-3) (see eq.10
of Becker et al., 2021) as

$\Gamma=\dot{N}_{\mathrm{glob}}\epsilon\sigma_{\text{H\,{i}}~{}}\lambda_{\mathrm{MFP}},$
(2)

where $\sigma_{\text{H\,{i}}~{}}=6.3\times 10^{-18}\mathrm{cm}^{-2}$ is the
hydrogen ionisation cross-section,
$\epsilon=\alpha_{s}/(\alpha_{\mathrm{bg}}+2.75)\approx 0.5$ is the efficiency
parameter and $\lambda_{\mathrm{MFP}}$ is the mean free path of ionising
photons. Given the constraints of the UVB value by Calverley et al. (2011) and
the MFP measurements of Becker et al. (2021); Zhu et al. (2023); Roth et al.
(2023); Satyavolu et al. (2023), we estimate an emissivity parameter of
$\zeta=50$ is required for the simulations used within this body of work to
have completed reionisation around $z=6$. We use this value for all
simulations in this paper.

###  2.2 Modelling the effect of small-scale absorbers.

All simulations in this work, listed in Tab. 1, are identical in terms of
density fields, halo lists and source model. However, they differ in the
implementation of the effect of small-scale absorbers. In this section, we
describe the five different ways in which we include their effect.

The r40 simulation models the MFP of ionising photos due to the presence of
small-scale absorbers by limiting the distance photons are allowed to travel
from a source. Such a spherical cut-off barrier is easy to implement, reduces
the cost of the radiative transfer by limiting it to a certain distance and is
also typically how semi-numerical codes based on the excursion set model of
Furlanetto et al. (2004) implement the effect. Examples of such codes are the
original 21CMFAST code (Mesinger & Furlanetto, 2007), simfast21 (Santos et
al., 2010) and Sem-Num (Majumdar et al., 2014). Once a ray reaches the
distance $R_{\mathrm{max}}$, its associated photons are assumed to have been
absorbed by an unresolved LLS and are removed. For the r40 simulation,
$R_{\mathrm{max}}=40$ cMpc was chosen. This value is used at all redshifts.

The C2 simulation instead models the effect of small-scale absorbers by
boosting the recombination rate through the introduction of a clumping factor,
$C=\langle\rho^{2}\rangle/\langle\rho\rangle^{2}$ where $\rho$ is the density
of matter. The choice of global clumping factor is motivated by the common
conception that it is sufficient to model small-scale absorbers (Wu et al.,
2021, and references within), as most ionising photons are absorbed in the
fully ionised (and hence recombining) outer layers of self-shielded systems.
Here we use a simple global clumping factor and set its value to $C=2$ at all
times. Note that this clumping comes in addition to the clumping already
resolved within the simulation from the baryon density fluctuations. Mao et
al. (2020) compared the use of global clumping factors and density-dependent
ones in simulations similar to the ones we use here and showed that global
clumping factors are a reasonable approximation for density-dependent
clumping. Bianco et al. (2021) showed that the impact of stochasticity on the
density-dependent clumping is relatively small. In order to limit the
computational cost and stop photons from crossing more than half a volume
size, the C2 simulation also includes a cut-off barrier of
R${}_{\mathrm{max}}=70$ cMpc. The choice of the cut-off is motivated by the
MFP of ionising photos at the post-EoR reported by W14 at $z=5$.

The three remaining simulations model the effect of small-scale absorbers by
adding additional, subgrid, opacity for ionising photons, characterized by a
parameter $\lambda_{\mathrm{MFP}}$, the distance at which this additional
opacity corresponds to an optical depth $\tau_{\mathrm{LL}}$ of 1 for Lyman
limit photons. This additional opacity is frequency-dependent following the
frequency-dependence of the cross-section of neutral hydrogen, approximately
$\nu^{-3}$. A value of $\tau_{\mathrm{LL}}=1$ implies that around 63 per cent
($1-e^{-\tau}$) of the photons at the Lyman limit will have been removed by
the subgrid opacity at $\lambda_{\mathrm{MFP}}$. For higher frequency photons
this fraction will be lower. Note that this is an additional opacity as the
neutral hydrogen on the grid will also remove ionising photons. Unlike the
cut-off barrier $R_{\mathrm{max}}$ model, absorption due to small-scale
structure is not instantaneous at a fixed distance from the source, but occurs
gradually at all distances from a source. We therefore will sometimes refer to
it as a ”diffuse barrier” model, although no real barrier exists in this
approach. As the small-scale absorbers are distributed throughout the IGM,
this model is more physical than the cut-off barrier approach. It was
previously described in Shukla et al. (2016) and in the appendix of Georgiev
et al. (2022) and is similar in philosophy to the approach described in Davies
& Furlanetto (2022).

Simulations LLS40, eLLS and pLLS all use this ”diffuse barrier” model for the
effect from small-scale absorbers. The difference is that in LLS40,
$\lambda_{\mathrm{MFP}}$ is uniform and fixed at 40 cMpc for all redshifts. In
eLLS, $\lambda_{\mathrm{MFP}}$ is also uniform but evolving according to

$\lambda^{\mathrm{W}}_{\mathrm{MFP}}=37\mathrm{\,pMpc}\left(\frac{1+z}{5}\right)^{-5.4}\,,$
(3)

a fit taken from Worseck et al. (2014). As
$\lambda^{\mathrm{W}}_{\mathrm{MFP}}\approx$ 40 cMpc at z=6, the LLS40 and
eLLS simulations have the same MFP at that redshift. The pLLS simulation uses
the same evolving $\lambda^{\mathrm{W}}_{\mathrm{MFP}}$ fit to determine the
average subgrid opacity, but distributes the opacities proportional to the sum
of the virial cross-sections of halos in a cell, which makes the subgrid
opacity position dependent. It assumes that the small-scale absorbers are
mostly associated with dark matter halos (see Theuns & Chan, 2024). This
approach was described in more detail in sec. 3.3 of Shukla et al. (2016).

Table 1: C2-Ray simulations with a minimal mass of ionising sources of $M_{\mathrm{min}}=10^{9}M_{\odot}$, an emissivity parameter of $\zeta=50$, and their respective models of the LLS. Label | MFP Method | $\lambda_{\mathrm{MFP}}$
---|---|---
r40 | spherical cut-off barrier | 40 cMpc
C2 | Global clumping of C = 2 | 70 cMpc
eLLS | uniform absorption | equation (3)
pLLS | position-dependent absorption | equation (3)
LLS40 | uniform absorption | 40 cMpc
Figure 1: Top Panel: Reionisation histories for the models in Table 1.
Constraints on the ionisation histories from Ly$\alpha$ equivalent width
(stars) (Mason et al., 2018, 2019; Hoag et al., 2019; Bolan et al., 2022;
Bruton et al., 2023; Morishita et al., 2023; Jones et al., 2023), Ly$\alpha$
dark pixel fraction (circles) (McGreer et al., 2015; Jin et al., 2023), and
Ly$\alpha$ emitters (diamonds) (Ouchi et al., 2010) as well as the damping
wings from quasars (squares) (Greig et al., 2019; Yang et al., 2020; Greig et
al., 2022) and galaxies (pentagons)(Hsiao et al., 2023; Umeda et al., 2023).
Middle panel: The Thomson optical depth constraints as reported by Planck
Collaboration et al. (2020). Lowest Panel: The evolution of the UVB
background, compared against the measurements from Calverley et al.
(2011),Wyithe & Bolton (2011), and Gaikwad et al. (2023). The inclusion of LLS
sub-grid models both delays the final stages of the EoR and regulates the
evolution of the UVB. Moreover, boosting the clumping factor throughout the
simulation (C2 model, in yellow) results in overshooting the UVB at EndEoR,
whilst a low cut-off of the MFP (r40 model, in green) limits its evolution.

##  3 Effect of small-scale absorbers on reionisation and the ultraviolet
background

In this section, we will compare the results of the five simulations in terms
of their reionisation history and the evolution and statistical properties of
the UVB.

###  3.1 The global parameters of reionisation

In Figure 1, we plot the global observables characterising our five
reionisation models. This figure also introduces the colour scheme for these
simulations which we will use throughout this paper. The top panel shows the
evolution of the volume-averaged neutral fraction $\bar{x}_{\rm\text{H\,{i}}}$
with observational constraints (see Bruton et al., 2023; Keating et al.,
2023). All simulations are consistent with these observational constraints.
The plot also shows that of the five simulations only r40 completes
reionisation by redshift 6. Due to limitations of the $N$-body results, we
were unable to continue the reionisation simulations beyond $z=6$.

The middle panel shows the CMB Thomson optical depths against the measurement
of $\tau=0.051\pm 0.006$ from Planck Collaboration et al. (2020). The
reionisation histories of all models are consistent with those measurements.

In the bottom panel of Fig. 1, we plot the average UVB as it evolves with
redshift during the EoR. We calculate the average UVB by normalising the
ionisation rate output from our simulations with respect to the ionisation
fraction ($x_{\rm\text{H\,{ii}}}$) in each cell following
$\sum_{i}\Gamma_{i}x_{\mathrm{H_{II}},i}/\sum_{i}x_{\mathrm{H_{II}},i}$. The
measurements of the UVB shown are as reported by Calverley et al. (2011);
Wyithe & Bolton (2011); Gaikwad et al. (2023). This plot shows considerable
differences in the evolution of the UVB, which we will now discuss in some
detail.

How LLS are modelled in our simulations directly dictates the evolution of the
average UVB. First, we consider the C2 model (shown in yellow), which has a
global clumping factor of $C=2$ and cut-off scale of $R_{\mathrm{max}}=70$
cMpc. The addition of a global clumping parameter results in a boost of the
recombination rate (see Bianco et al., 2021). A large value of the clumping
factor will, therefore, impact the UVB and be more pronounced in denser
regions of the simulation. This becomes evident by examining the evolution of
the UVB for the C2 model compared to the remaining models. The first ionising
sources are located in the regions of the highest density. The clumping factor
will increase the already large recombination rates there and thus limit the
growth of these initial H ii regions. This initially delays the EoR, evident
in the neutral faction’s slower evolution and the lower UVB amplitude at
$z>8$. Despite this, below $z=8$, the clumping model shows the strongest
evolution in the UVB and overshoots the observations at $z=6$. This is due to
two key factors. After $\bar{x}_{\rm\text{H\,{i}}}\approx 0.8$, the ionised
bubbles expand past their local over-densities and into the lower density IGM,
in addition to a considerable overlap of ionised regions (Lidz et al., 2008).
Regions previously in equilibrium with the ionisation output of their local
sources are affected by the radiation of adjacent sources, boosting the UVB.
Moreover, as $R_{\mathrm{max}}=70$ cMpc for the C2 model, the ionising
radiation can in principle propagate further than in our other models. This
behaviour becomes significant past the overlap phase at the end stages of
reionisation when most cells within the simulation host ionising sources and
the remaining neutral islands reside in the deepest voids.

The opposite effect on the UVB can be seen when considering the r40 model,
which has no additional clumping but a cut-off spherical barrier of
$R_{\mathrm{max}}=40$ cMpc. Compared to C2, the cut-off barrier model (seen in
green in the bottom panel of Fig. 1) initially has a higher UVB value, which
increases during reionisation and eventually asymptotes around $z\approx 6.3$.
During the early phases reionisation in r40 is more rapid than in C2 as there
is no additional clumping to impede the growth of the H ii bubbles. However,
towards the end stages of the EoR, the characteristic scale of H ii bubbles
approaches the scale of the cut-off barrier, which impedes the further growth
of these regions111Note that H ii bubbles will continue to grow in size due to
overlap, yet no single source can create an ionised structure larger than the
cut-off barrier (see appendix of Georgiev et al., 2022).. Moreover, towards
the end of reionisation, most neutral islands are far away from ionising
sources, as they are located in the voids (Giri et al., 2019, 2024). The sharp
cut-off of 40 cMpc results in fewer ionising photons reaching the outskirts of
the neutral islands. Therefore, once a neutral island becomes ionised, its
local ionisation rate is lower, and the residual neutral fraction is higher
than in models with larger $R_{\mathrm{max}}$ values. Because of this, most
recently ionised cells in these models will have a lower ionisation rate,
pushing down the global mean. In addition to r40, we have run models with
lower cut-off barriers of 10 and 20 cMpc to confirm this effect. While low
values of the cut-off barrier can be motivated by observations, the MFP is the
mean value around a broad distribution, and in this approach, the free paths
larger than the cut-off barrier are not modelled. We argue that this will also
impact the evolution of the UVB.

For all the diffuse barrier models, the UVB is slower to evolve than in r40,
and reionisation is delayed (also see fig.2 of Shukla et al., 2016). The
redshift-dependent opacity from the LLS in eLLS (in fuchsia) further slows the
reionisation process by impeding the growth of H ii bubbles compared to LLS40.
The addition of the position and density dependence in pLLS (in red) further
exacerbates this effect such that regions with more halos will have a higher
opacity due to LLS in the voids of the simulation. However, as the EoR
progresses and ionised bubbles overlap more and more, the evolution of the
mean UVB for LLS40, eLLS, and pLLS exhibit a steeper gradient compared to C2
and r40 at the EndEoR and converge at $z=6$. This behaviour follows the same
line of thinking as the role of the global clumping parameter. However,
because the additional opacity removes ionising photons, assumed to be
absorbed by unresolved LLS, it does not affect the recombination rate.
Therefore, it does not result in a boost of the UVB as seen in the C2 model.
This effect is most pronounced in the pLLS model since the additional position
and density-dependent opacity causes dense regions to be more optically thick
to ionising photons. On the other hand, low-density regions which host fewer
halos will have less additional opacity added for the pLLS model, compared to
LLS40 and eLLS. Because of this difference, the reionisation history seen at
the top panel of Fig. 1 is gradual for pLLS but becomes steeper at the end of
reionisation.

In summary, for the C2 model the clumping initially suppresses the evolution
of the UVB, however, it overshoots its amplitude compared to the measurements
at $z=6$. The role of the absorption due to clumping diminishes as the overlap
of H ii regions becomes significant and UV radiation is free to propagate up
to 70 cMpc. Meanwhile, the diffuse barrier models allow a fraction of the
ionising photons to permeate past the defined MFP while the r40 model imposes
a firm upper limit. For the r40 case, the sizes of H ii bubbles are sharply
constrained by the imposed cut-off, and the UVB exhibits an asymptotic trend
towards the EndEoR.

To conclude the discussion of the evolution of the mean UVB we would like to
point out that the observed trends in the diffuse barrier models are more
consistent with the UVB evolution reported in high-resolution state-of-the-art
hydrodynamic simulations. This can be seen by comparing to the results shown
in fig. 5 of Garaldi et al. (2022) (THESAN) and fig. 1 of Lewis et al. (2022)
(CODA III). However, it should be remembered that those simulations consider
volumes which are around $(4)^{3}$ times smaller than our simulations and
hence do not sample the same range of environments.

Figure 2: PDFs of the logarithm of the UVB for the set of models showcased in
Tab. 1. Note that UVB is normalised against the mean of the distribution at
the given redshift. Values at the right of each colour bar are the redshift,
whilst those to the left are the volume-averaged neutral fraction. The figure
at the bottom right row includes the UVB statistics for each model at a fixed
redshift of 6. A value of unity is plotted as a vertical line, to highlight
the peak of each PDF. The role of the LLS model directly influences the shape
of the UVB PDFs. The right tail of each PDF (values smaller than the mean) is
sensitive to small-scale fluctuation induced by the implementation of the
small-scale absorbers. The effect is most pronounced in the pLLS model (middle
row, right column, and seen in red), where the LLS are both position and
redshift-dependent. Conversely, the left tail of each PDF (values larger than
the mean) is sensitive to the cumulative effect of the LLS models, i.e. the
overall MFP. For example, the r40 model (middle row, left column, and seen in
green) imposes a fixed barrier of 40 cMpc, freezing the evolution of the
large-scale UVB fluctuations (note each PDF is centred around the mean value
at the given redshift).  Figure 3: Light cones of the reionisation models
presented in Tab. 1 during the EndEoR, where the $y$-axis for each panel
covers a redshifts period lower than $z=6.6$, while the $x$-axis is the length
of simulation. The neutral fraction is plotted in a logarithmic scale, the
colour scheme is used to highlight the fluctuations of the ionised and neutral
IGM. For example, the dark black regions are indicative of neutral hydrogen
islands, and the red regions indicate the ionised medium, particularly the
residual neutral hydrogen. Fluctuations in the ionised IGM are pronounced in
shades of yellow and green, highlighting the role of the small-scale
absorbers. Because in ionisation equilibrium the ionisation rate scales as
$\Gamma\sim 1/x_{\mathrm{HI}}$ these panels also provide a good indication of
the UVB fluctuations shown in Fig. 2.

###  3.2 Statistical analysis of the UVB

In this section, we expand our analysis of the role of small-scale absorbers
on the UVB from Sec. 3.1 by going beyond the single average value and instead
examine the Probability Density Function (PDF) of the field. In Fig. 2, we
present the PDF of the logarithm of the ionisation rate for all ionised cells
($x_{\rm\text{H\,{ii}}}\geq 0.5$). The PDF has been normalised against the
logarithm of the mean UVB value at the given redshift to highlight the
evolution of the fluctuations (i.e.
$\log_{10}(\Gamma)/\log_{10}(\overline{\Gamma})$). Hence a ratio below unity
(leftward for each panel in Fig. 2) is indicative of a UVB larger than the
global mean and vice versa. As we saw above, while all the source models used
in this paper are identical, the various implementations of the effect of
small-scale absorbers result in a diverse range of reionisation histories. We
present the PDFs for each model for the final stages of reionisation (roughly
$z\leq 6.6$ and $\bar{x}_{\rm\text{H\,{ii}}}\geq 0.8$). Moreover, we include
light cone slices of the neutral fraction of hydrogen for our models in Figure
3. Compared to coeval slices of the simulation, where both axes are dimensions
of length, the light cones show the evolution of the neutral fraction across a
redshift period over the length of the simulation, as the $y$ and $x$-axis,
respectively. This serves to better visualise the role of small-scale
absorbers over time, specifically in the case of the redshift-dependent models
eLLS and pLLS.

All models exhibit similar PDF distributions at high redshifts (and,
therefore, high global neutral fractions), as seen in the lightest shades in
Fig. 2. The distributions are initially skewed to values higher than the mean
(left), while there are significant outliers of the PDFs for values lower than
the mean (right). In this period, while reionisation is significantly
underway, the remaining neutral islands in the cosmic voids significantly
affect the ionisation rates in low and high-density regions alike and,
therefore, the fluctuations in the UVB.

For example, we can imagine two adjacent regions, one of high and one of low
density, separated by an area of neutral hydrogen. The high-density region
within the simulation will likely host multiple clusters of ionising sources.
Its local UVB value will be higher than that of the region with fewer sources
(for example, see the lower left area of all panels in Fig. 3). The lower-
density region can be in a state of local ionisation equilibrium if it is
sufficiently shielded by the neutral islands. Naturally, as the neutral island
in our example disappears with the progress of EoR, and the two regions become
interconnected, radiation from the cluster of sources interacts with the IGM
in the low-density region, pushing the local UVB to evolve to the global mean,
reducing the fluctuations in the UVB distribution. In this regime, the role of
small-scale absorbers (and how they are modelled) is most evident in the UVB
fluctuations in the PDFs. The LLS model not only limits the global limit of
how far ionising photons can travel (the MFP), but it also affects the UVB
fluctuations inside the ionised region (i.e. small scales).

We can observe multiple trends due to the nature of each method of modelling
the small-scale absorbers. We first discuss the sides of the PDFs of our
models throughout the EndEoR. The shift of the peak of the PDFs with a
redshift to their mean value is indicative of the evolution of the ionisation
rate to a state of equilibrium. Moreover, we can initially associate the left
side of the PDF with high-density regions inside ionised bubbles, occupied by
numerous sources, which in the case of the pLLS model are sensitive to the
presence of small-scale absorbers. Meanwhile, the right side of the PDF is
more sensitive to the low-density regions, which interact with the ionisation
font, which is modulated by the MFP.

For the models r40,LLS40,C2 in the left column of Fig. 2, we note that the
outliers on the right side of each PDF decrease with redshift, indicating the
reduction of large-scale fluctuations in the UVB. The common factor between
the diffuse barrier (blue model), the cut-off barrier (green model), and the
clumping model (orange model) is that they do not directly model the small-
scale absorbers, resulting in a more homogeneous local UVB field, which
adjusts to the global UVB value as overlap occurs. This is especially true for
the r40 and LLS40, which model the LLS by placing a non-evolving upper limit
on the MFP. On the other hand, the clumping factor in the C2 model accounts
for the position-dependence of small-scale absorbers. Rather than being
absorbed by the self-shielded regions photons are free to propagate across the
simulation, leading to an increased UVB value as discussed in Sec. 3.1.
Conversely, the models on the right column of Fig. 2 (eLLS,pLLS), exhibit
large-scale fluctuations, which do not decrease with the EndEoR. For the
redshift and density-dependent model pLLS, this is a natural consequence of
adding the small-scale absorbers, which modulate the large-scale UVB. Adding a
redshift dependence of the LLS evolution in the eLLS model induces this effect
compared to the LLS40 model, as its increasing MFP value limits the large-
scale UVB fluctuations throughout reionisation.

We now turn our attention to the fluctuations seen on the right side from the
mean value of each PDF222Remember that the PDFs are centred around their
evolving mean value., which is sensitive to the small-scale fluctuations of
the UVB and the presence of the small-scale absorbers. The previous point is
now more apparent by comparing the r40 model against the LLS40, eLLS, and
pLLS. The cut-off barrier model places an upper bound on the MFP. Hence, all
small-scale fluctuations are smoothed over within a given region. Thus, the
PDF’s right hand side does not significantly evolve compared to the mean
value. A diffuse barrier model, such as in LLS40, removes more local photons
due to the absorption by self-shielded regions. Additionally, the diffuse
barrier allows photons from distant sources to propagate further than those in
the cut-off barrier in r40. Moreover, the evolving diffuse barrier in eLLS
reproduces this effect and induces large-scale fluctuations in the UVB.

The most interesting effect can be observed in the pLLS model, as compared to
the clumping model C2. By accounting for the density/position and redshift
dependence of the LLS, these systems induce small-scale UVB fluctuations and
regulate the large-scale UVB (as seen in the top middle panel of Fig. 3,
considering that in equilibrium $\Gamma\sim 1/x_{\text{H\,{i}}~{}}$). The
resulting PDF closely follows a log-normal distribution. While the clumping
model induces small-scale fluctuations, as seen in the lower left of Fig. 2,
it cannot modulate the large-scale UVB fluctuations. The C2 simulation results
in a PDF which is strongly peaked around the mean value (or in other words has
a high kurtosis), compared to the other models.

Examining the lower right panel of Fig. 2, we can summarise the role of the
MFP on the UVB distribution. The role of small-scale absorbers is two-fold as
they will profoundly affect the small-scale absorptions as well as modulate
the MFP of ionising photons and, hence, the large-scale UVB. By modelling the
presence of small-scale absorbers, we can reproduce the large-scale UVB
fluctuations. However, the PDF converges to a log-normal distribution for the
pLLS model, highlighting that the UVB is sensitive to the outliers induced by
the small-scale absorbers. Furthermore, the C2 model cannot accurately
reproduce the large-scale UVB fluctuations. In the following sections, we show
how these effects on the UVB will also fundamentally impact the observables of
reionisation.

##  4 Observables of the End of Reionisation

After having analysed the reionisation histories and the UVB properties of our
set of simulations, we will now consider three key observables for the EndEoR,
namely the fluctuations in the Ly$\alpha$ opacity, the MFP for ionising
photons and the 21-cm power spectrum from the EoR. Measurements for the first
two already exist, the 21-cm power spectrum will hopefully be measured in the
near future.

###  4.1 The Gunn-Peterson optical depth

A probe for the state of the IGM along the line of sight can be found by
analysing the Ly$\alpha$ transmission in quasar spectra (Fan et al., 2006).
The observed large-scale fluctuations in high-$z$ quasar spectra are an active
area of exploration (see Gaikwad et al., 2023, and references within).
Interpretation of the results coupled with numerical modelling indicates that
the fluctuations cannot solely be explained by a uniform UVB in a varying
density field (see Ishimoto et al., 2022, and references within). Possible
explanations of the fluctuations which have been proposed are the effect of
density fluctuations in the IGM and temperature fluctuations in the IGM (eg.
D’Aloisio et al., 2015; Keating et al., 2018). Davies & Furlanetto (2016)
found that a fluctuating UVB shaped by the presence of self-shielding systems
plays a significant role in their models. Including the fluctuating MFP
substantially improves agreement with observations. Therefore, modelling the
effective Ly$\alpha$ optical depth ($\tau_{\mathrm{eff}}$) can not only
constrain models of reionisation but can also be used as a probe of the role
of small-scale absorbers. In this section, we simulate mock observations of
the Ly$\alpha$ effective optical depth and study how the choice of MFP
implementation and the presence of neutral islands impact these observations.

The Ly$\alpha$ optical depth can be written as follows (e.g. Pagano & Liu,
2021):

$\tau_{\mathrm{Ly}\alpha}=\int_{z_{\mathrm{emi}}}^{z_{\mathrm{obs}}}n_{\mathrm{HI}}\sigma_{\mathrm{Ly}\alpha}(\lambda,T)\frac{cdt}{dz}dz,$
(4)

where $n_{\mathrm{HI}}$ is the column density of neutral hydrogen,
$\sigma_{\mathrm{Ly}\alpha}$ is the Ly$\alpha$ cross section, which depends on
the wavelength $\lambda$ and the temperature of the system $T$. To properly
describe the Ly$\alpha$ absorption profile and its temperature dependence, we
model it with a Voigt profile using the methodology from Tepper-García (2006)

$\sigma_{\mathrm{Ly}\alpha}(\lambda,T)=\frac{e^{2}\lambda_{i}^{2}f_{i}\Gamma_{i}}{4\sqrt{\pi}m_{e}cb^{2}}H(a,x),\quad
a=\lambda^{2}\Gamma_{i}/4\pi cb.$ (5)

$\Gamma_{i}$ and $\lambda_{i}$ are the damping constant and the wavelength of
Ly$\alpha$ while $b=\sqrt{2kT/m}$ is the Doppler parameter. In this study, we
assign the temperature of the ionised cells following the methodology
described in Qin et al. (2021) (see eq.8) as follows

$T_{g}^{\gamma}=T_{\mathrm{ion,I}}^{\gamma}\bigg{[}\bigg{(}\frac{\mathcal{Z}}{\mathcal{Z_{\mathrm{ion}}}}\bigg{)}^{3}\frac{\rho_{b}}{\rho_{\mathrm{b,ion}}}\bigg{]}^{\frac{2\gamma}{3}}\frac{\mathrm{exp}(\mathcal{Z}^{2.5})}{\mathrm{exp}(\mathcal{Z}^{2.5}_{\mathrm{ion}})}+T_{\mathrm{lim}}^{\gamma}\frac{\rho_{b}}{\overline{\rho_{b}}},$
(6)

where $\mathcal{Z}=(1+z)/7.1$ and equation of state index is $\gamma=1.7$
while $\rho_{b}$ is the baryon density (see McQuinn & Upton Sanderbeck, 2016,
for a detailed analysis). The subscript ’ion’ indicates the quantity at the
ionisation redshift $z_{\mathrm{ion}}$, $T_{\mathrm{ion,I}}$ is the
temperature after the I-front has propagated, which is assumed to be constant
$T_{\mathrm{ion,I}}=2\times 10^{4}$ K. The final relaxation temperature is
$T_{\mathrm{lim}}^{\gamma}=1.775\mathcal{Z}\times 10^{4}$ K.333Neutral cells
are assumed to have a temperature of around 200K. Additionally, $H(a,x)$ is
the Voigt-Hjerting function:

$H(a,x)=\frac{a}{\pi}\int_{-\infty}^{\infty}\frac{\mathrm{exp}[-y^{2}]}{(x-y)^{2}+a^{2}},\quad\mathrm{where}\quad
y=v/b.$ (7)

We measure the transmitted flux in our simulations averaged over $N$ cells.
See sec. 5 of B22 for more details. The effective optical depth444Note that we
calculate the Ly$\alpha$ optical depth from the neutral hydrogen on the grid
and do not include the effect of the subgrid neutral hydrogen assumed to be
present and absorbing ionising photons in the LLS40, eLLS and pLLS
simulations. This allows us to make a direct comparison with the other two
simulations, r40 and C2, which do not introduce subgrid H i . However, it does
mean we may be underestimating the actual Ly$\alpha$ opacity implied by these
models. is given as

$\tau_{\mathrm{eff}}=-\ln\bigg{(}\frac{1}{N}\sum_{i}\mathrm{exp}\big{[}-\tau_{\mathrm{Ly}\alpha}(i)\big{]}\bigg{)},$
(8)

Our results are computed for transmission over 30 cMpc skewers, to avoid the
evolution of the signal along the LOS in the observable, while we also
generate 70 cMpc skewers to calculate the fit parameters from Tab. 2 to match
the fit values from B22. We sample $5\times 10^{3}$ rays for each simulation
in this study, taking a random direction on the $x$,$y$, and $z$ planes to
construct the skewer. Note we have also conducted checks with different sample
sizes as well as sampling the rays at random locations versus at the position
of the ionising sources. In both cases, we find our results to be consistent.

####  4.1.1 Cumulative Distribution Functions of the effective Lyman-$\alpha$
opacity (CDFs)

Figure 4: The Cumulative Distribution Functions (CDFs) of the effective
Ly$\alpha$ opacity across 30 cMpc during the final stages of reionisation. The
models and their nomenclature are described in Tab. 1. The dark shaded regions
are the 1$\sigma$ uncertainties of the measurements reported by B22 (for
$z=5.433,5.525,5.620,5.716,5.814,5.915,6.017,6.122,6.230$). The skewing of
each CDF with increasing redshift is indicative of lines of sight with a
higher Ly$\alpha$ opacity, hinting that the EoR is still underway. All
simulated skewers chosen within this analysis are averaged over 30 cMpc, to
match the reported data. While each thick line is the mean distribution of
$\tau_{\mathrm{eff}}$, the 10 thin lines represent random LOS, which have been
chosen to highlight the variation of the distribution. The bottom right panel
shows the simulated and real data at redshift 6, normalised against the mean
value of each distribution. Including redshift and position dependence in the
LLS implementation (pLLS model) leads to a more accurate model of
observations. The r40 model (left middle row in green) is quantitatively
excluded as it fails to reproduce high effective Ly$\alpha$ opacity at the end
of reionisation. Meanwhile, models such as eLLS and C2 (right top row in
purple and left bottom row in yellow, respectively) can marginally reproduce
the observables.

We present the Cumulative Distribution Function (CDF) for our set of models
against the data reported by B22 (seen in grey) in Fig. 4. We focus our
discussion on the final stages of reionisation, visualising output for
redshifts where the global neutral fraction is $\bar{x}_{\rm\text{H\,{i}}}\leq
0.20$, for the first three rows of Fig. 4 (thick lines). To highlight the
effect of sample variance in our analysis we visualise randomly selected CDFs
(thin lines), choosing 10 for clarity. The analysis in Gaikwad et al. (2023)
reports a neutral fraction of $\bar{x}_{\rm\text{H\,{i}}}\approx 0.15$ at a
redshift 6. While we do not directly match our models with the observations,
both data sets describe a universe in the final stages of reionisation.
Therefore, we can study the shape of the functions to gain insight into the
distributions and study the different implementations of the effects of small-
scale absorbers. Hence, in the right panel of the third row of Figure 4, we
include each output of our simulation closest to the observed CDF, all of
which are normalised against their mean value.

Lastly, we also perform a best-fit analysis of the data following the
parametrisation555The original form of the equation includes a third parameter
$C$, which we do not consider as we equate
$\tau_{\mathrm{eff}}(1+z_{0})=\tau_{0}$. The $C$ parameter is more relevant
for data with implicit noise, which is not included in our simulations. of
$\tau_{\mathrm{eff}}$ reported in eq.3 of B22 that is

$\tau_{\mathrm{eff}}(1+z)=\tau_{0}\bigg{(}\frac{1+z}{1+z_{0}}\bigg{)}^{\beta},$
(9)

where $z$ and $\beta$ are free parameters, and $z_{0}=6$. Our best-fit values
are listed in Tab. 2. The CDF results of B22 at redshifts of 5 rise rapidly
from zero to unity, indicative of Gaussian distribution, and become slanted
with increasing redshift due to the positive skewness of distribution
attributed to the presence of neutral hydrogen. Each $\tau_{\mathrm{eff}}$ is
the cumulative value of a spatial skewer, which itself captures diverse
environments of different density, temperature, and ionisation states. Hence,
the effective value will capture some average quantity of the IGM across the
line of sight. A low value is indicative of a low-density region, while high
values are thought to indicate a higher concentration of residual neutral
hydrogen. How the peak of the distribution evolves is therefore related to
both the evolution of the UVB as well as the disappearance of H i islands. For
example, at high-$z$, the 1$\sigma$ distribution reported by B22 (rightmost
grey CDFs of Fig. 4) are positively skewed, hinting at an IGM with a higher
fraction of residual neutral hydrogen as well as the existence of neutral
islands.

The mean of the distribution is shifted due to the evolving residual neutral
fraction, while the presence of neutral islands results in a higher
$\tau_{\mathrm{eff}}$, which positively skews the data. The r40 model (in
green, middle row, left column) does not exhibit the described skewness, with
a visibly slight variance from sight-line to sight-line. Following Section
3.2, the implementation of the cut-off barrier results in ionised regions with
a high-amplitude, ubiquitous, and slowly evolving UVB value across redshift.
The fluctuations in the UVB, and hence the Ly$\alpha$ optical depth, are then
set either by the density field, the 40 cMpc barrier to ionising photons, or
the neutral islands. The slow evolution of the UVB, noted in Fig. 1, results
in a weaker redshift revolution of $\tau_{\mathrm{eff}}$ seen in the low value
of $\beta=4.2\pm 0.3$ Tab. 2, compared to the observations (where
$\beta=13.7\pm 1.5$). Moreover, the relatively small scatter in the random
sight lines is due to the suppression of the fluctuations in the IGM, which
is, on average, more irradiated compared to the other models.

Conversely, the pLLS model is marked by a rapid evolution of the CDF against
redshift, which is consistent with the values reported by B22. Although the
simulation has not been completely reionised in Fig. 4 (middle row, right
column, in red), we can observe this rapid evolution in tandem with the
positive skewness of the distribution measured in observations. As discussed
in Section 3.2, introducing the effect of small-scale absorbers results in
fluctuations of the UVB in dense regions. While we cannot compare directly to
observations, we note that the shape of the CDF at redshift 6 in the middle
right panel of Fig. 4, is most similar in shape and form to that of the data.
Lastly, the rapid evolution of the lower values of $\tau_{\mathrm{eff}}$ seen
in the pLLS model is more consistent with those observed in the data shown in
grey.

The eLLS and LLS40 simulations exhibit similar characteristics to the pLLS
model. However, the evolution of the CDFs is impeded by excluding the density-
dependence of the LLS (which results in small-scale fluctuations in the UVB
for pLLS as discussed in Section 3.2), resulting in a slower evolution of the
low values of $\tau_{\mathrm{eff}}$. Nevertheless, eLLS is broadly consistent
in shape and evolution with the global features of the Ly$\alpha$ CDFs, as
evidenced by the redshift evolution of $\beta=11.5\pm 0.7$, closest of the
remaining models to the value derived from observations.

Lastly, the clumping simulation C2 reproduces the broad features of the
measurements for high Ly$\alpha$ opacities, similar to the data. However, the
evolution of the mean value is slower than those noted in observations
($\beta=10.4\pm 0.8$,), and the evolution of the lowest values of
$\tau_{\mathrm{eff}}$ appears to be impeded.

Our conclusion here follows the one from Section 3.2; the model in which the
effect of small-scale absorbers depends on both redshift and position
reproduces the observations better. The cut-off barrier model fails to
reproduce the high effective Ly$\alpha$ opacity across EndEoR. Models such as
eLLS and C2 can marginally reproduce the observables while pLLS is most
successful in reproducing the features inferred from the observational data.

Table 2: Best fit parameters of equation (9) of the effective Ly$\alpha$ opacities ($\tau_{\mathrm{eff}}$), derived for 70 cMpc skewers to match the reported values from B22. Note that we choose $z_{0}=6$ for our fit, which is our lowest available redshift. Label | $\tau_{o}$ | $\beta$
---|---|---
B22 | 0.3 $\pm$ 0.08 | 13.7 $\pm$ 1.5
r40, cut-off fixed barrier | 2.3 $\pm$ 0.1 | 4.2 $\pm$ 0.3
C2, global clumping parameter of two | 2.5 $\pm$ 0.1 | 10.4 $\pm$ 0.8
LLS40, diffuse fixed barrier | 2.5 $\pm$ 0.1 | 8.5 $\pm$ 0.6
eLLS, diffuse evolving barrier | 2.7 $\pm$ 0.1 | 11.5 $\pm$ 0.7
pLLS, density-dependent diffuse barrier | 2.8 $\pm$ 0.1 | 13.3 $\pm$ 0.6

###  4.2 The mean free path of ionising photons from the EoR.

Figure 5: Upper Panel: PDF distribution of the free path of ionising photons
within the simulations reported in Tab. 1. The thick solid lines are the PDF
functions at redshift 6, while the thin dashed lines are normalised to a fixed
neutral fraction of 15 %. The colour scheme and naming notation are reported
across this body of work. Bottom Panel: Distribution for the cut-off barrier
model r40 as well as the PDFs of the free paths of the diffuse barrier models,
excluding the role of the additional opacity added to model the LLS.
Conversely, the thin red line is the same analysis conducted only on the
small-scale absorbers used in the pLLS model, and is consistent with the thin
black line is the reported MFP value by W14. On the one hand, the r40 and C2
models cannot reproduce the small values of the free paths reported in Becker
et al. (2021); Bosman (2021); Gaikwad et al. (2023). On the other hand, the
diffuse barrier models, LLS40 and eLLS, can match the MFP. It is only the pLLS
model seen at the top panel in red, which reproduces a PDF of the MFP, which
both matches the mean value and has a large variation, favouring small free
paths. Figure 6: Evolution of the median value of the free path evolution over
the end stages of reionisation as well as the standard deviation around the
mean value (shaded areas). The labels and colouring scheme follow that in Fig.
5. The black dashed line is the MFP evolution depicted in W14, whilst the red
dashed line is the MFP of the LLS model implemented for the pLLS simulation.
It is only the free path distributions from the diffuse barrier models (LLS40,
eLLS,pLLS) in blue, fuchsia, and red, respectively, which agree with the
measurements from Gaikwad et al. (2023).

We recreate the ray-tracing approach described in Section 4.1 and follow the
Lyman continuum photons in the ionised medium during the EndEoR for a given
LOS. We calculate the opacity to Lyman continuum radiation as

$\tau_{\mathrm{LyC}}=\int_{0}^{1}n_{\text{H\,{i}}~{}}\epsilon\sigma_{\text{H\,{i}}~{}}dx,$
(10)

where $n_{\text{H\,{i}}~{}}$ and $\sigma_{\text{H\,{i}}~{}}$ are the number
density of hydrogen (see eq. (2)). Following the definition of the MFP of
ionising photons, we trace the cumulative opacity to ionising radiation up to
a value of unity ($\tau_{\mathrm{LyC}}\approx 1$) from the positions of the
ionising sources in our simulation. This allows us to probe multiple free
paths (FP), which probe the IGM at the end of reionisation, exploring the role
of residual neutral hydrogen along the lines-of-sight (LOS) and the remaining
extended neutral islands. For the sake of consistency we use exactly the same
random $2.5\times 10^{5}$ rays when analysing each of our five simulations.
This way we can better contextualise the direct effect that the implementation
of the effect of small-scale absorbers plays in determining the MFP
measurements. Additionally, because the diffuse barrier approach model the
effect as an additional subgrid component within our simulation (see Section
3.1), simulations LLS40, eLLS and pLLS allow us to perform our analysis either
by removing or including the contribution of the subgrid small-scale
absorbers. While the absorbers have a time-dependent effect on the opacity in
the IGM, performing this analysis allows us to directly probe how the LLS
model can affect the MFP observable.

The resulting PDFs are presented in the upper panel of Fig. 5, for a fixed
redshift of 6 (solid) or a fixed neutral fraction of approximately 15 %
(dashed). The bottom panel includes the same analysis of the diffuse barrier
simulations excluding the subgrid opacity from the small-scale absorbers.

What is immediately apparent for the free path PDFs is that the clumping
simulation C2 at redshift 6 produces a distribution, which is Gaussian and has
a mean of around 80 cMpc ($\approx 11$ pMpc). Initially, it would be tempting
to relate the large value of the MFP to the low neutral fraction of the
simulation ($\approx 2$ %), however, the PDF for the same model at a higher
redshift and neutral reaction (yellow dashed line), remains biased towards
large free paths in contrast to current measurements (e.g. Becker et al.,
2021). Rather, the behaviour of the distribution is related to the residual
neutral hydrogen in the IGM. The additional opacity to the IGM, indirectly
introduced through the boosted recombination rate, is insufficient in reducing
the length of the free paths of the ionising photons (also see Mao et al.,
2020, where much higher clumping values than here are used). Because
$R_{\mathrm{max}}=70$ cMpc, the majority of ionising sources during EndEoR are
free to interact with each other, further ionising the IGM between them and
boosting the value of the MFP.

A similar distribution is also observed for the r40 model seen in green in the
bottom panel of Fig. 5. Limiting the propagation of the ionising photons to 40
cMpc for r40 when conducting the radiative transfer does not directly
reproduce lower free path values when sampled in this manner. The cut-off
barrier, which serves as an upper limit of ionising photons, is ineffective at
the end of EoR, during which the majority of the cells in the simulation have
active sources of radiation, resulting in a highly ionised IGM. Note that the
shape of this PDF and most distributions are skewed towards low free path
values. We can intuitively imagine the overall PDF of the free paths as a
convolution of two different distributions. On the one hand, for small values
of the free paths, the function follows a Poisson distribution, sensitive to
the rays absorbed by neutral islands in the vicinity of the ionising source.
On the other hand, for large values of the free paths of ionising photons, the
distribution resembles a normal distribution, as it probes the ionised
intergalactic medium, which is typically in a state of equilibrium during the
final stages of reionisation.

The same result is observed to a varying degree for the diffuse barrier models
LLS40, eLLS, pLLS in the bottom panel of Fig. 5 (solid lines), for the case
when the additional subgrid absorption due to small-scale absorbers is removed
from the analysis. Despite this, we can see that by including redshift and
density dependence in the model of the LLS, the resulting distribution (even
with the aforementioned absorbers being removed) converges to distributions
skewed to lower values and a lower MFP value. We can, therefore, infer that by
including the additional opacity from small-scale absorbers, we not only
remove the excess photons which should have interacted with the self-shielded
systems but also regulate the overall state of the ionised IGM. Naturally, the
more physically motivated the implementation of the effect of these absorbers,
the stronger the effect on the simulation.

As an additional result, we can sample the free paths between the density-
dependent absorbers included in the pLLS simulation in seen at the bottom
panel of Fig. 5 as a dashed red line. We confirm that the produced
distribution peaks at a $\lambda_{\mathrm{MFP}}=48$ cMpc as it has been
normalised to match the fit from W14. In the case of the additional opacity
for the LLS40 and eLLS, the additional opacity would peak as a delta function
(as $\tau=1$ at $\lambda_{\mathrm{MFP}}=40$ cMpc by definition).

It is only when we convolve both of the red functions (dashed and solid) for
the pLLS model from the bottom panel of Fig. 5, that we arrive at the true
distribution in the top panel of the same figure (similarly for LLS40 and
eLLS). Unlike the other diffuse barrier models, whose distributions become
shifted and centred at the predefined MFP value, the position-dependent model
produces a log-normal distribution skewed towards smaller free paths for
ionising photons. At the same time, some rays can travel to as much as 60
cMpc. This model is more natural than LLS40 and eLLS, both of which, due to
the scaling, disfavour low and high values of the free paths alike. The
complete PDF produced by the pLLS run naturally reproduces a large variety of
possible sight-lines such as those reported in fig.2 of Bosman (2021), some
more consistent with the ionised universe, while others are more sensitive to
the presence of small-scale absorbers or neutral islands. Crucially, this is
confirmed when examining the distributions for a fixed neutral fraction
(dashed lines in the top panel), highlighting the conclusion that a position-
dependent model of the LLS affects both the small and large-scale physics of
the IGM. This is in concordance with the UVB results in sec. 3.2.

Lastly, in Fig. 6, we plot the evolution of the mean free path of ionising
photons for each of our simulations. We choose to represent the evolution of
the median value of each distribution but also include the $1\sigma$ around
the mean value to better capture the full range of the statistics. In
addition, we compare the MFP evolution from W14 (from equation (3), seen in
black) to that derived from our small-scale absorbers within the simulation.
The inferred evolution slightly differs in amplitude and is biased towards
larger values of the MFP. This can be attributed to the normalisation, which
assumes that the global opacity must be normalised against the reported mean
free path. As expected, the mean values derived for each of the diffuse models
are significantly lower than for the remaining models. Most notably, the
reported measurements of the MFP from Bosman (2021) and Gaikwad et al. (2023)
are consistent within the distribution of free paths from the eLLS and pLLS
models.

###  4.3 The 21-cm signal from the epoch of reionisation

In this section, we turn our attention to the 21-cm power spectrum from
reionisation. Compared to the observables discussed in Sec. 8 and 4.2, which
probe the ionised IGM, the 21-cm signal is sensitive to the neutral hydrogen
present in the IGM. We examine the impact LLS modelling has on the shape and
amplitude of power spectrum.

####  4.3.1 Differential brightness temperature

The 21-cm line is the result of the spin-flip transition, occurring in the
ground state of neutral hydrogen. The 21-cm radiation observed by radio
interferometry telescopes is defined as the differential surface temperature
against the CMB. In the Rayleigh-Jeans limit it can be written as follows
(e.g. Mellema et al., 2013)

$\delta
T_{\mathrm{b}}(\mathbf{r},z)=\frac{1-e^{-\tau(\mathbf{r},z)}}{1+z}(T_{s}(\mathbf{r},z)-T_{\mathrm{CMB}}(z))\,,$
(11)

where $\tau$ is the optical depth of the gas and $T_{s}$ is the spin
temperature. The spin temperature characterises the level populations of the
21-cm transition and has to be different from $T_{\mathrm{CMB}}$ to produce an
observable signal. At the redshifts of reionisation $T_{s}$ is expected to be
similar to the kinetic temperature of the gas due to the frequent interaction
of the intergalactic H i with Ly$\alpha$ photons produced by star-forming
galaxies (Field, 1959; Wouthuysen, 1952). X-ray heating caused by the same
galaxies is expected to increase this kinetic temperature (Pritchard &
Furlanetto, 2007).

Assuming an optically thin gas $\tau\ll 1$ at high redshifts, Eq. 11 can be
rewritten as

$\displaystyle\delta
T_{\mathrm{b}}(\mathbf{r},z)=T_{0}(\mathbf{r},z)\bar{x}_{\rm\text{H\,{i}}}(1+\delta_{\text{H\,{i}}~{}}(\mathbf{r}))(1+\delta_{\rho}(\mathbf{r}))\
,$ (12) $\displaystyle T_{0}(\mathbf{r},z)\approx
27\bigg{(}\frac{1+z}{10}\bigg{)}^{1/2}\bigg{(}\frac{T_{\mathrm{s}}(\mathbf{r},z)-T_{\mathrm{CMB}}(z)}{T_{\mathrm{s}}(\mathbf{r},z)}\bigg{)}$
$\displaystyle\bigg{(}\frac{\Omega_{\mathrm{b}}}{0.044}\frac{h}{0.7}\bigg{)}\bigg{(}\frac{\Omega_{\mathrm{m}}}{0.27}\bigg{)}^{-1/2}\bigg{(}\frac{1-Y_{\mathrm{p}}}{1-0.248}\bigg{)}~{}\mathrm{mK},$
(13)

Here the variable $\rho$ is the mass density, whilst
$\delta_{\rho}=\rho/\bar{\rho}-1$ is the corresponding fluctuation. Similarly,
$\delta_{\text{H\,{i}}~{}}$ denotes the fluctuation of the neutral fraction
$x_{\rm\text{H\,{i}}}$ field. The term $T_{0}$ collects all the cosmological
parameters as well as the spin temperature dependence. Note that within this
work we do not consider the effects of redshift space distortions. We analyse
the 1D spherically averaged power spectra of the 21 cm signal shown in Eq. 12.
We also define the dimensionless power spectrum at a certain wave-number $k$
as $\Delta^{2}(k)=k^{3}P(k)/(2\pi^{2})$.

A large variety of upper limit values on the 21 cm power spectrum have been
reported by several low-frequency radio interferometers (at various redshifts
and scales), e.g., the Low-Frequency Array (LOFAR)666www.lofar.org (Mertens et
al., 2020), the Hydrogen Epoch of Reionization Array
(HERA)777https://reionization.org (Abdurashidova et al., 2022; HERA
Collaboration et al., 2023), the Murchison Widefield Array
(MWA)888www.mwatelescope.org (Trott et al., 2020; Yoshiura et al., 2021), the
Giant Metrewave Radio Telescope (GMRT)999www.gmrt.ncra.tifr.res.in (Paciga et
al., 2013), as well as the forthcoming Square Kilometre Array
(SKA)101010www.skatelescope.org (Koopmans et al., 2015). However, no detection
has been made and the derived upper limits set only weak constraints on
astrophysical quantities (Abdurashidova et al., 2022; HERA Collaboration et
al., 2023).

####  4.3.2 The 21 cm power spectrum as a probe of the neutral IGM during the
final stages of reionisation

Figure 7: The 21 cm power spectrum plots of the models are presented in Tab. 1
againts $k$-scale. The upper panel show the power spectrum at a fixed redshift
of 6 for all models. The implementation of the LLS for the diffuse barrier
models (LLS40, eLLS, pLLS) results in a delay of the reionisation process. The
higher fraction of neutral hydrogen present in each simulation affects the
amplitude of the power spectra, accordingly. Note that for the r40, where a
cut-off barrier of $R_{\mathrm{max}}=40$ cMpc is employed across the EoR,
results in a steeper slope of the 21 cm power spectrum at low-$k$. The lower
panel instead compares each 21 cm power spectrum at a fixed global neutral
fraction of 9 %. The resulting 21 cm power spectra have more or less similar
amplitudes due to the fixed global neutral fraction but are sensitive to the
distribution of the neutral islands. A large amplitude of the power spectrum
at low-$k$ (hence, high physical scale), such as in the C2 model, indicates
that the majority of the neutral islands in the simulation are fewer and
concentrated in the furthest voids (see the lower right panel of Fig. 3).

In Fig. 7, we present the 21 cm power spectrum, both at a fixed redshift of 6
(upper panel) and as the closest possible global neutral fraction of around 9
% (lower panel). The amplitude of the 21 cm power spectrum is notably higher
for the diffuse barrier models. Because the diffuse barrier models produce a
delayed EoR, more neutral hydrogen is available at redshift 6 (as seen at the
bottom $y$-axes of the panels in Fig. 3). This is most apparent when comparing
the power spectra at a fixed neutral fraction in the bottom panel of Fig. 7,
where the overall quantity of H i is similar across all simulations. In doing
so, we can reduce the role of the global neutral fraction when comparing our
results, and see a quantitative decrease in the ratios. Despite this, the
spatial distribution of the neutral islands also plays an important role in
the 21 cm power spectrum.

Taking the r40 model, for example, the bubble growth is restricted up to 40
cMpc across all of the EoR. The effect of this is most pronounced when the H
ii bubble size reaches the cut-off scale, which in the case for r40 is at
redshifts $z<6.3$ (seen as the flattening of the UVB amplitude at Sec. 4). At
a redshift of 6 (the upper panel of Fig. 7), the low-$k$ 21-cm power spectrum
has a steeper slope compared to the other models, as traces the density field,
while larger $k$-values are influenced by the ionisation field (see sec. 5 of
Georgiev et al., 2022). In comparison, in the LLS40 model, where the barrier
is also fixed yet diffuse, the transition scale is smoother since 37 % of the
photons (at the Lyman limit) can travel further than the imposed 40 cMpc (see
Section 1). Note that, in general, both models produce similar results when
compared at a similar global neutral fraction (bottom panel of Fig. 7).
Conversely, the C2 model does not display a similar feature and indeed has
more power on large scales, when compared in the bottom panel of Fig. 7.

Hence, we have also analysed the same outputs using the MFP-BSD methodology
from Mesinger & Furlanetto (2007) to probe the H i island sizes at a fixed
global neutral fraction. Note that, we confirm that the C2 simulation produces
slightly larger ionised regions, however, we find the method cannot
distinguish between the distribution of the r40, LLS40, eLLS, pLLS models.
Therefore, by allowing the photons to freely propagate across the simulations
and instead adjusting the recombination factor, the resulting neutral islands
at the end of reionisation are mainly larger in size and group at the furthest
voids of the simulation, resulting in more power in large scales. The eLLS and
pLLS models look largely consistent with one another, as their main difference
is centred on how the small-scale absorbers are modelled within the ionised
regions. Compared to the remaining models, these diffuse barrier models
produce neutral islands which are smaller in size and more scattered across
the simulation due to the fluctuations introduced by the LLS.

##  5 Conclusion and discussion

This body of work examines the end stages of the Epoch of Reionisation
(EndEoR) using a set of 3D radiative transfer simulations. The final stages of
reionisation are particularly difficult to model numerically, due to the
significant amount of ionising sources and the complex topology of its ionised
overlapping structures. Particularly, modelling the UVB within such
simulations often overestimates the measurements reported by Calverley et al.
(2011). Despite this, the EndEoR is a treasure trove of information about the
process of reionisation, and the first stars and galaxies. Utilising
information from high redshift quasars, the reported low values of the mean
free path of ionising photons at a redshift of 6 (Becker et al., 2021; Zhu et
al., 2023) as well as higher values of the effective Ly$\alpha$ opacity along
the line-of-sight (Bosman, 2021) indicate at an extended reionisation period
(Keating et al., 2018; Kulkarni et al., 2019). In addition, the large
Gunn–Peterson trough reported by Becker et al. (2015) hints at the presence of
islands of neutral hydrogen, which can be utilised to study reionisation
through the 21 cm signal (Giri et al., 2019).

This work aims to explore the role of small-scale optically thick self-
shielded systems, which are often not included in large-scale simulations of
reionisation. Referred to as Lyman Limit systems, these systems are often
associated with hydrogen present in the outskirts of galaxies and/or in the
filaments between them (Theuns & Chan, 2024; Fan et al., 2024). LLS have an
effect on the observables of reionisation, such as the effective Ly$\alpha$
opacity, the mean free path of ionising (Lyman continuum) photons, as well as
the 21 cm signal from the intergalactic medium. We investigate several methods
to model these small-scale absorbers with varying assumptions (see Section 2.2
and Table 1). All four C2-Ray models are identical in set-up save for the
modelling of the LLS, described below:

  1. 1.

In the r40 model ionising photons are cut off at a fixed distance of 40 cMpc
for all redshifts. This can be seen as a model which only considers the large-
scale effect of the LLS on reionisation.

  2. 2.

The LLS40 model uses an additional subgrid opacity to impede the propagation
of ionising photons. Following the definition of the mean free path of
ionising photons, at 40 cMpc the cumulative additional opacity is equal to
unity. Hence, this barrier is uniformly distributed and porous, as roughly 37
% of the photons (at the Lyman limit) can permeate past the boundary
condition. Unlike the r40 model, the LLS40 will also impact the ionisation
rate within the H ii regions.

  3. 3.

The eLLS model uses the same additional opacity approach as LLS40 but evolves
the value of the mean free path according to the fit from W14, see eq. 3. It
thus captures the likely redshift dependence of the number of LLS within the
simulation, although calibrated on post-EoR observables.

  4. 4.

In the pLLS we use the same MFP fit as in eLLS but distribute the additional
opacity according to the a recipe based on the cross-section of dark matter
halos, rather than uniformly (see Shukla et al., 2016, for a description).
This model captures both the likely evolution of the number of LLS as well as
the probable connection between LLS and collapsed structures.

  5. 5.

Lastly, model C2 attempts to model the presence of small-scale absorbers by
boosting the global recombination rate by a factor of two. We include this
method, as it is an often-used solution within the community.

The conclusions of our study can be formulated as follows:

  1. 1.

We find that by including the presence of LLS in our simulations, we can
naturally reproduce delayed reionisation histories. Including redshift and/or
density dependence in our LLS models further extends the EoR. Moreover, we
quantify how the role of clumping in our C2 model results in an overshooting
of the UVB as well as how the choice of a cut-off barrier for the r40 case can
produce nonphysical features in our observables. Conversely, the diffuse
barrier models LLS40, eLLS, pLLS are more consistent with the evolution of the
UVB reported by high-resolution small-scale simulations of reionisation such
as THESAN (Garaldi et al., 2022) and CODA III (Lewis et al., 2022).

  2. 2.

By conducting a statistical analysis of the UVB for our set of models in Sec.
3.2, we highlight that the role of Lyman Limit Systems and how they are
modelled will affect reionisation globally as well as locally. Models such as
r40, LLS40, which only model the impact of the value of the MFP on ionising
photons, affect the UVB on the larges-scales (captured by the right side of
the PDFs in Fig. 2) but not the small-scales, where LLS are not directly
modelled, and the ionisation rate is higher. Including the redshift evolution
of the LLS in eLLS marginally improves on this. Crucially, the density and
redshift-dependent pLLS model, which models the LLS on the cell level, results
in small-scale fluctuation in the ionised regions, resulting in a distinct
log-normal shape of the PDF.

We can summarise that the properties of the Lyman Limit systems affect the
global parameters of reionisation. In the last section of the paper, we study
whether these findings are consistent with observations at the EndEoR.

  1. 3.

We conduct mock observations of the effective Ly$\alpha$ optical depth
($\tau_{\mathrm{eff}}$) and compare the features of our results against the
measurements from the QSO sample from B22. In Section. 4.1 (see Figure 4), we
confirm our findings in points 1. and 2. necessitating the role of the proper
modelling of small-scale absorbers during reionisation. While the eLLS,C2
models marginally reproduce the features of the observable, only the pLLS
model can match the redshift evolution of the effective Ly$\alpha$ optical
depth ($\beta=13.3\pm 0.6$) from the measurements ($\beta=13.7\pm 1.3$).

  2. 4.

We further our study focusing on direct observations of the mean free path of
ionising photons, which as reported by Becker et al. (2021) and Gaikwad et al.
(2023) is particularly low at a redshift of 6 and another indicator of a
delayed reionisation. In Section 6 we sample lines of sight from our
simulations and construct a PDF of the free paths of ionising photons at a
redshift of 6 or a fixed neutral fraction of 15 % (see Fig. 5). We find that
models such r40 and C2 overshoot the value of the mean free path as their
implementations and cannot induce the required opacity to limit Lyman limit
photons. And while all diffuse barrier models (LLS40, eLLS, pLLS) reproduce
low values of the MFP at a redshift of 6 $\lambda_{\mathrm{MFP}}\approx 25$
cMpc, their distributions differ massively due to the nature of their
implementations. Notably, the pLLS model follows a log-normal distribution
skewed towards low values of the MFP, and both redshift-dependent absorption
models (pLLS and eLLS) have distributions with standard deviations overlapping
within the MFP observations (see Fig. 6).

  3. 5.

Lastly, we probe the effect of Lyman Limit Systems on the neutral medium
during EndEoR by examining the 21 cm signal, particularly the 21 cm power
spectrum as shown in Fig. 7. At redshift 6, we find differences in amplitude
by an order of magnitude between models such r40 and pLLS, which can be
attributed to the delay of the reionisation process discussed in 1. Note, the
flat shape of the r40 is consistent with the discussion in Georgiev et al.
(2022). When comparing the 21 cm power spectra at a fixed neutral fraction of
9 %, we can further peer into key differences between models. We find that the
clumping model has a higher large-scale amplitude of the 21 cm signal,
indicative of larger neutral islands. Since the clumping model has no imposed
MFP, the majority of neutral islands trace the furthest voids from the
clusters of ionised sources. And while the LLS40 model imposes a barrier of 40
cMpc, its permissivity allows for ionising photons to traverse past it,
resulting in a similar effect on the neutral islands. Conversely, the r40 is
seen to have power across all $k$-scales, because most neutral islands will be
located in regions 40 cMpc away from ionising sources.

Lastly, both eLLS and pLLS have a similar shape and amplitude. This occurs
because the H i islands are less sensitive to the modelling differences of the
small-scale absorbers within the ionised bubbles, but rather on the large-
scale mean free path imposed by the LLS and its redshift evolution.

In summary, we find that modelling the LLS during the epoch of reionisation
has a crucial role in regulating the duration of reionisation as well its
observables such as the UVB background, Ly$\alpha$ effective optical depth,
mean free path of ionising photons as well as the shape of neutral islands
detentions via the 21 cm signal. In concordance with Georgiev et al. (2022) we
conclude that cut-off barrier models such as r40, which only place an upper
limit on the MFP value, produce inconsistent results with the current
knowledge within the field. Meanwhile, the additional clumping added in the C2
allows the model to marginally reproduce some of the features from the
$\tau_{\mathrm{eff}}$ observations. However, it is insufficient in reducing
the average distance of UV photons, resulting in higher values of the UVB and
MFP at redshift 6, compared to observations. This also impacts the 21-cm power
spectrum, which when compared to all models at a fixed neutral fraction, has
more power on low-$k$ scales. Introducing a uniform opacity diffuse barrier
like in LLS40 resolves many of the aforementioned problems, but is not fully
consistent with observations of the evolution of the Ly$\alpha$ optical depth
observations. Lastly, the evolving diffuse barrier model eLLS improves on the
latter and reproduces the observables adequately, but it is only when we
account for the density-dependence of the absorbers such as in pLLS is
sufficient to reproduce the measurements of $\tau_{\mathrm{eff}}$ and the MFP.
When comparing the 21-cm power spectrum at $\bar{x}_{\rm\text{H\,{i}}}\approx
0.09$, both eLLS and pLLS are identical, while LLS40 is shown to match the r40
model and have more power on low-$k$ scales. The addition of a redshift and
position dependence to the model of small-scale absorbers is shown to delay
the EndEoR (comparing LLS40,eLLS,pLLS respectively).

Although we covered a wide range of approaches for including the effect of
unresolved absorbers, there remain further models, to investigate. For
example, the evolving hard barrier model from Giri et al. (2024) can be
combined with the halo model in Theuns & Chan (2024), which reproduces the MFP
evolution from W14, although the disadvantages of the hard barrier approach as
shown in the current paper will also be present in this model. A more
promising approach is the density-dependent clumping models reported in
Sobacchi & Mesinger (2014); Bianco et al. (2021). We will investigate the
sensitivity of the parameters of this fit for models such as pLLS. Moreover,
we will extend our simulations to a redshift of $z\approx 5.5$ to further
explore the EndEoR, however, and explore the role of our source model on the
duration of reionisation (e.g. Cain et al., 2023).

## Acknowledgements

IG acknowledges the support from the Gustaf och Ellen Kobbs stipendiestiftelse
and the Alva and Lennart Dahlmark Research Grant. GM’s research is supported
by the Swedish Research Council project grant 2020-04691_VR. Nordita is
supported in part by NordForsk. Other results were obtained on resources
provided by the National Academic Infrastructure for Supercomputing in Sweden
(NAISS) at PDC (DARDEL), Royal Institute of Technology, Stockholm. We have
utilised the following python packages for manipulating the simulation outputs
and plotting results: numpy (Harris et al., 2020), scipy (Virtanen et al.,
2020), matplotlib (Hunter, 2007) and tools21cm111111https://github.com/sambit-
giri/tools21cm package (Giri et al., 2020).

## Data Availability

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

## References

  * Abdurashidova et al. (2022) Abdurashidova Z., et al., 2022, ApJ, 925, 221
  * Alvarez & Abel (2012) Alvarez M. A., Abel T., 2012, ApJ, 747, 126
  * Aubert et al. (2015) Aubert D., Deparis N., Ocvirk P., 2015, MNRAS, 454, 1012
  * Bañados et al. (2018) Bañados E., et al., 2018, Nature, 553, 473
  * Becker et al. (2001) Becker R. H., et al., 2001, AJ, 122, 2850
  * Becker et al. (2015) Becker G. D., Bolton J. S., Madau P., Pettini M., Ryan-Weber E. V., Venemans B. P., 2015, MNRAS, 447, 3402
  * Becker et al. (2021) Becker G. D., D’Aloisio A., Christenson H. M., Zhu Y., Worseck G., Bolton J. S., 2021, MNRAS, 508, 1853
  * Bianco et al. (2021) Bianco M., Iliev I. T., Ahn K., Giri S. K., Mao Y., Park H., Shapiro P. R., 2021, MNRAS, 504, 2443
  * Bolan et al. (2022) Bolan P., et al., 2022, MNRAS, 517, 3263
  * Bolton & Haehnelt (2007) Bolton J. S., Haehnelt M. G., 2007, MNRAS, 382, 325
  * Bosman (2021) Bosman S. E. I., 2021, arXiv e-prints, p. arXiv:2108.12446
  * Bosman et al. (2022) Bosman S. E. I., et al., 2022, MNRAS, 514, 55
  * Bruton et al. (2023) Bruton S., Lin Y.-H., Scarlata C., Hayes M. J., 2023, ApJ, 949, L40
  * Cain et al. (2023) Cain C., D’Aloisio A., Lopez G., Gangolli N., Roth J. T., 2023, arXiv e-prints, p. arXiv:2311.13638
  * Calverley et al. (2011) Calverley A. P., Becker G. D., Haehnelt M. G., Bolton J. S., 2011, MNRAS, 412, 2543
  * D’Aloisio et al. (2015) D’Aloisio A., McQuinn M., Trac H., 2015, ApJ, 813, L38
  * D’Aloisio et al. (2020) D’Aloisio A., McQuinn M., Trac H., Cain C., Mesinger A., 2020, ApJ, 898, 149
  * Davies & Furlanetto (2016) Davies F. B., Furlanetto S. R., 2016, MNRAS, 460, 1328
  * Davies & Furlanetto (2022) Davies F. B., Furlanetto S. R., 2022, MNRAS, 514, 1302
  * Davies et al. (2018) Davies F. B., et al., 2018, ApJ, 864, 143
  * Davies et al. (2024) Davies F. B., et al., 2024, ApJ, 965, 134
  * Dayal & Ferrara (2018) Dayal P., Ferrara A., 2018, Physics Reports, 780, 1
  * Dixon et al. (2016) Dixon K. L., Iliev I. T., Mellema G., Ahn K., Shapiro P. R., 2016, MNRAS, 456, 3011
  * Endsley & Stark (2022) Endsley R., Stark D. P., 2022, MNRAS, 511, 6042
  * Fan et al. (2006) Fan X., Carilli C. L., Keating B., 2006, ARA&A, 44, 415
  * Fan et al. (2024) Fan J., Zhu H., Avestruz C., Gnedin N. Y., 2024, ApJ, 963, 45
  * Field (1959) Field G. B., 1959, ApJ, 129, 536
  * Finlator et al. (2009) Finlator K., Özel F., Davé R., 2009, MNRAS, 393, 1090
  * Furlanetto et al. (2004) Furlanetto S. R., Zaldarriaga M., Hernquist L., 2004, ApJ, 613, 1
  * Gaikwad et al. (2023) Gaikwad P., et al., 2023, MNRAS, 525, 4093
  * Garaldi et al. (2022) Garaldi E., Kannan R., Smith A., Springel V., Pakmor R., Vogelsberger M., Hernquist L., 2022, MNRAS, 512, 4909
  * Georgiev et al. (2022) Georgiev I., Mellema G., Giri S. K., Mondal R., 2022, MNRAS, 513, 5109
  * Giri & Mellema (2021) Giri S. K., Mellema G., 2021, Monthly Notices of the Royal Astronomical Society, 505, 1863
  * Giri et al. (2019) Giri S. K., Mellema G., Aldheimer T., Dixon K. L., Iliev I. T., 2019, MNRAS, 489, 1590
  * Giri et al. (2020) Giri S., Mellema G., Jensen H., 2020, The Journal of Open Source Software, 5, 2363
  * Giri et al. (2024) Giri S. K., Bianco M., Schaeffer T., Iliev I. T., Mellema G., Schneider A., 2024, arXiv preprint arXiv:2403.04838
  * Greig et al. (2017) Greig B., Mesinger A., Haiman Z., Simcoe R. A., 2017, MNRAS, 466, 4239
  * Greig et al. (2019) Greig B., Mesinger A., Bañados E., 2019, MNRAS, 484, 5094
  * Greig et al. (2022) Greig B., Mesinger A., Davies F. B., Wang F., Yang J., Hennawi J. F., 2022, MNRAS, 512, 5390
  * Gunn & Peterson (1965) Gunn J. E., Peterson B. A., 1965, ApJ, 142, 1633
  * HERA Collaboration et al. (2023) HERA Collaboration et al., 2023, ApJ, 945, 124
  * Harnois-Déraps et al. (2013) Harnois-Déraps J., Pen U.-L., Iliev I. T., Merz H., Emberson J. D., Desjacques V., 2013, MNRAS, 436, 540
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357
  * Hayes & Scarlata (2023) Hayes M. J., Scarlata C., 2023, ApJ, 954, L14
  * Hinshaw et al. (2013) Hinshaw G., Larson D., Komatsu E., et al. 2013, ApJS, 208, 19
  * Hoag et al. (2019) Hoag A., et al., 2019, ApJ, 878, 12
  * Hsiao et al. (2023) Hsiao T. Y.-Y., et al., 2023, arXiv e-prints, p. arXiv:2305.03042
  * Hunter (2007) Hunter J. D., 2007, Computing in Science & Engineering, 9, 90
  * Iliev et al. (2006) Iliev I. T., et al., 2006, MNRAS, 371, 1057
  * Iliev et al. (2014) Iliev I. T., Mellema G., Ahn K., Shapiro P. R., Mao Y., Pen U.-L., 2014, MNRAS, 439, 725
  * Ishimoto et al. (2022) Ishimoto R., et al., 2022, MNRAS, 515, 5914
  * Jin et al. (2023) Jin X., et al., 2023, ApJ, 942, 59
  * Jones et al. (2023) Jones G. C., et al., 2023, arXiv e-prints, p. arXiv:2306.02471
  * Keating et al. (2018) Keating L. C., Puchwein E., Haehnelt M. G., 2018, MNRAS, 477, 5501
  * Keating et al. (2020) Keating L. C., Weinberger L. H., Kulkarni G., Haehnelt M. G., Chardin J., Aubert D., 2020, MNRAS, 491, 1736
  * Keating et al. (2023) Keating L. C., Bolton J. S., Cullen F., Haehnelt M. G., Puchwein E., Kulkarni G., 2023, arXiv e-prints, p. arXiv:2308.05800
  * Koopmans et al. (2015) Koopmans L., et al., 2015, in Advancing Astrophysics with the Square Kilometre Array (AASKA14). p. 1 (arXiv:1505.07568), doi:10.22323/1.215.0001
  * Kulkarni et al. (2019) Kulkarni G., Keating L. C., Haehnelt M. G., Bosman S. E. I., Puchwein E., Chardin J., Aubert D., 2019, MNRAS, 485, L24
  * Lewis et al. (2022) Lewis J. S. W., et al., 2022, MNRAS, 516, 3389
  * Lidz et al. (2008) Lidz A., Zahn O., McQuinn M., Zaldarriaga M., Hernquist L., 2008, ApJ, 680, 962
  * Majumdar et al. (2014) Majumdar S., Mellema G., Datta K. K., Jensen H., Choudhury T. R., Bharadwaj S., Friedrich M. M., 2014, Monthly Notices of the Royal Astronomical Society, 443, 2843
  * Mao et al. (2020) Mao Y., Koda J., Shapiro P. R., Iliev I. T., Mellema G., Park H., Ahn K., Bianco M., 2020, MNRAS, 491, 1600
  * Mason et al. (2018) Mason C. A., et al., 2018, ApJ, 857, L11
  * Mason et al. (2019) Mason C. A., et al., 2019, MNRAS, 485, 3947
  * McGreer et al. (2015) McGreer I. D., Mesinger A., D’Odorico V., 2015, MNRAS, 447, 499
  * McQuinn & Upton Sanderbeck (2016) McQuinn M., Upton Sanderbeck P. R., 2016, MNRAS, 456, 47
  * McQuinn et al. (2011) McQuinn M., Oh S. P., Faucher-Giguère C.-A., 2011, ApJ, 743, 82
  * Mellema et al. (2006a) Mellema G., Iliev I. T., Alvarez M. A., Shapiro P. R., 2006a, New Astron., 11, 374
  * Mellema et al. (2006b) Mellema G., Iliev I. T., Pen U.-L., Shapiro P. R., 2006b, MNRAS, 372, 679
  * Mellema et al. (2013) Mellema G., et al., 2013, Experimental Astronomy, 36, 235
  * Mertens et al. (2020) Mertens F. G., et al., 2020, MNRAS, 493, 1662
  * Mesinger & Furlanetto (2007) Mesinger A., Furlanetto S., 2007, ApJ, 669, 663
  * Monaghan (1992) Monaghan J. J., 1992, ARA&A, 30, 543
  * Morishita et al. (2023) Morishita T., et al., 2023, ApJ, 947, L24
  * Mortlock et al. (2011) Mortlock D. J., et al., 2011, Nature, 474, 616
  * Ocvirk et al. (2020) Ocvirk P., et al., 2020, MNRAS, 496, 4087
  * Ouchi et al. (2010) Ouchi M., et al., 2010, ApJ, 723, 869
  * Paciga et al. (2013) Paciga G., et al., 2013, MNRAS, 433, 639
  * Pagano & Liu (2021) Pagano M., Liu A., 2021, MNRAS, 508, 1915
  * Planck Collaboration et al. (2014) Planck Collaboration Ade P. A. R., Aghanim N., Armitage-Caplan C., Arnaud M., et al. 2014, A&A, 571, A16
  * Planck Collaboration et al. (2020) Planck Collaboration et al., 2020, A&A, 643, A42
  * Pritchard & Furlanetto (2007) Pritchard J. R., Furlanetto S. R., 2007, MNRAS, 376, 1680
  * Pritchard & Loeb (2012) Pritchard J. R., Loeb A., 2012, Reports on Progress in Physics, 75, 086901
  * Qin et al. (2021) Qin Y., Mesinger A., Bosman S. E. I., Viel M., 2021, MNRAS, 506, 2390
  * Roth et al. (2023) Roth J. T., D’Aloisio A., Cain C., Wilson B., Zhu Y., Becker G. D., 2023, arXiv e-prints, p. arXiv:2311.06348
  * Santos et al. (2010) Santos M. G., Ferramacho L., Silva M., Amblard A., Cooray A., 2010, Monthly Notices of the Royal Astronomical Society, 406, 2421
  * Satyavolu et al. (2023) Satyavolu S., Kulkarni G., Keating L. C., Haehnelt M. G., 2023, arXiv e-prints, p. arXiv:2311.06344
  * Shukla et al. (2016) Shukla H., Mellema G., Iliev I. T., Shapiro P. R., 2016, MNRAS, 458, 135
  * Sobacchi & Mesinger (2014) Sobacchi E., Mesinger A., 2014, MNRAS, 440, 1662
  * Tepper-García (2006) Tepper-García T., 2006, MNRAS, 369, 2025
  * Theuns & Chan (2024) Theuns T., Chan T. K., 2024, MNRAS, 527, 689
  * Trott et al. (2020) Trott C. M., et al., 2020, MNRAS, 493, 4711
  * Umeda et al. (2023) Umeda H., Ouchi M., Nakajima K., Harikane Y., Ono Y., Xu Y., Isobe Y., Zhang Y., 2023, arXiv e-prints, p. arXiv:2306.00487
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Wang et al. (2020) Wang F., et al., 2020, ApJ, 896, 23
  * Watson et al. (2013) Watson W. A., Iliev I. T., D’Aloisio A., Knebe A., Shapiro P. R., Yepes G., 2013, MNRAS, 433, 1230
  * Worseck et al. (2014) Worseck G., et al., 2014, MNRAS, 445, 1745
  * Wouthuysen (1952) Wouthuysen S. A., 1952, AJ, 57, 31
  * Wu et al. (2021) Wu X., McQuinn M., Eisenstein D., 2021, J. Cosmology Astropart. Phys., 2021, 042
  * Wyithe & Bolton (2011) Wyithe J. S. B., Bolton J. S., 2011, MNRAS, 412, 1926
  * Xu et al. (2014) Xu Y., Yue B., Su M., Fan Z., Chen X., 2014, The Astrophysical Journal, 781, 97
  * Xu et al. (2016) Xu H., Wise J. H., Norman M. L., Ahn K., O’Shea B. W., 2016, ApJ, 833, 84
  * Yang et al. (2020) Yang J., et al., 2020, ApJ, 904, 26
  * Yoshiura et al. (2021) Yoshiura S., et al., 2021, MNRAS, 505, 4775
  * Zhu et al. (2023) Zhu Y., et al., 2023, ApJ, 955, 115