# The impact of varying inhomogeneous reionization histories on metrics of
Ly$\alpha$ opacity

Caitlin C. Doughty,1 Joseph F. Hennawi,2,1, Jose Oñorbe,3, Frederick B.
Davies,4 Zarija Lukić 5
1Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA, Leiden, the
Netherlands
2 Department of Physics, University of California, Santa Barbara, CA 93106,
USA
3 Facultad de Físicas, Universidad de Sevilla, Avda. Reina Mercedes s/n,
Campus Reina Mercedes, E-41012 Seville, Spain
4 Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg,
Germany
5 Lawrence Berkeley National Laboratory, CA 94720-8139, USA E-mail:
<EMAIL_ADDRESS>(CCD)

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The epoch of hydrogen reionization is complete by $z=5$, but its progression
at higher redshifts is uncertain. Measurements of Ly$\alpha$ forest opacity
show large scatter at $z<6$, suggestive of spatial fluctuations in neutral
fraction ($x_{\mathrm{HI}}$), temperature, or ionizing background, either
individually or in combination. However, these effects are degenerate,
necessitating modeling these physics in tandem in order to properly interpret
the observations. We begin this process by developing a framework for modeling
the reionization history and associated temperature fluctuations, with the
intention of incorporating ionizing background fluctuations at a later time.
To do this, we generate several reionization histories using semi-numerical
code AMBER, selecting histories with volume-weighted neutral fractions that
adhere to the observed CMB optical depth and dark pixel fractions.
Implementing these histories in the Nyx cosmological hydrodynamics code, we
examine the evolution of gas within the simulation, and the associated metrics
of the Ly$\alpha$ forest opacity. We find that the pressure smoothing scale
within the IGM is strongly correlated with the adiabatic index of the
temperature-density relation. We find that while models with 20,000 K
photoheating at reionization are better able to reproduce the shape of the
observed $z=5$ 1D flux power spectrum than those with 10,000 K, they fail to
match the highest wavenumbers. The simulated autocorrelation function and
optical depth distributions are systematically low and narrow, respectively,
compared to the observed values, but are in better agreement when the
reionization history is longer in duration, more symmetric in its distribution
of reionization redshifts, or if there are remaining neutral regions at $z<6$.
The systematically low variance likely requires the addition of a fluctuating
UVB.

###### keywords:

intergalactic medium – dark ages, reionization, first stars – quasars:
absorption lines

††pubyear: 2023††pagerange: The impact of varying inhomogeneous reionization
histories on metrics of Ly$\alpha$ opacity–The impact of varying inhomogeneous
reionization histories on metrics of Ly$\alpha$ opacity

##  1 Introduction

The epoch of reionization (EoR) is the last major phase transition in the
Universe’s history, when the neutral intergalactic medium (IGM) was ionized
for the first time since before recombination (e.g. Loeb & Barkana, 2001).
Initiating when the first stars began to form, it could not progress in
earnest until a suitable star formation rate density was achieved, since the
production rate of ionizing photons had to outstrip the rate of recombination.

While the amount of neutral material could be measured through observations of
hydrogen Lyman series absorption, the high mean density of the IGM at $z>3$
combined with the high oscillator strength of the Ly$\alpha$ transition at
$z=5$ means saturation occurs even for sub-percent neutral fractions,
preventing direct determination of its timeline. The best constraints on the
EoR timeline have come from measurements of the Cosmic Microwave Background
(CMB) from the Planck experiment, giving $\tau_{\mathrm{CMB}}=0.054\pm 0.007$
(Planck Collaboration et al., 2020). This was combined with assumed
reionization histories using the FlexKnot model (Millea & Bouchet, 2018),
which treats the timeline as a series of “knots” in $z$ and $x_{\mathrm{HI}}$
space, and used to establish an approximate timeline with a midpoint at
$z=7.7\pm 0.7$. The next best constraints will eventually come from 21 cm
experiments like the Hydrogen Epoch of Reionization Array (DeBoer et al.,
2017), which will be able to tomographically map neutral hydrogen throughout
the reionization epoch, allowing for a thorough assessment of the history.
However, while some experiments have preliminary results placing constraints
on pre-reionization X-ray heating (HERA Collaboration et al., 2023),
tomographic results are still a ways off, with much work currently going
towards diminishing significant foreground contamination from interlying
celestial structures as well as Earth-based signals.

Currently, there is a growing body of high-$z$ quasar and galaxy observations
that are also being used to chart the progression of reionization. QSO damping
wings at $7<z<7.5$ are producing a wide range of neutral fractions, from
almost 0 to 0.9 including the uncertainties, depending on the quasar and
analysis method (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). The large variation in these measurements, especially towards very low
neutral fractions, may indicate issues with the inference procedure, or that
the quasars are located in biased regions that are more ionized at early times
(e.g. Costa et al., 2014). Using galaxies in a similar way may hold promise,
though there is additional uncertainty with determining their intrinsic
spectra, which are more variable than those of quasars (Keating et al., 2023).
Additionally, IGM damping wings observed in galaxy spectra may be contaminated
by gas in the background galaxy’s intergalactic or circumgalactic media
(Heintz et al., 2023).

For the last stages of reionization occurring at $z<6$, it is possible to
measure a nonzero Ly$\alpha$ flux $F$, and learn about the neutral fraction by
considering the amount of “dark” $F\sim 0$ regions compared to those with
$F>0$. The most certain results are from McGreer et al. (2015), which by
quantification of the dark pixel fraction established upper limits on the
volume averaged abundance of neutral hydrogen to 0.09 at $z=5.6$. Other
studies have achieved model-dependent results through evaluation of Ly$\alpha$
and $\beta$ dark gaps in tandem with simulations (Zhu et al., 2022; Jin et
al., 2023).

The full distribution of Ly$\alpha$ opacities contains more information than
dark pixels alone, the values being impacted by the field values of neutral
fraction, gas density, temperature, and the ultraviolet background (UVB). For
a fully ionized IGM, the Ly$\alpha$ opacity ($-\ln{F}$) will vary with these
values according to

$\tau\propto n_{\mathrm{HI}}\propto T^{-0.7}\propto\Gamma_{\mathrm{HI}}^{-1}$
(1)

These measurements have been used to constrain the neutral fraction, mean free
path, and photoionization rate (Bosman et al., 2022; Gaikwad et al., 2023; Zhu
et al., 2023). The fluctuations in these fields also contribute to scale-
dependent variations in power, glimpsed through clustering measures like the
Ly$\alpha$ power spectrum and flux autocorrelation function (Viel et al.,
2013a; Walther et al., 2018; Boera et al., 2019; Karaçaylı et al., 2022).

In order to properly interpret observations, it is essential to compare to
predictions, and for such a complex process as reionization, this necessitates
the use of numerical simulations. These simulations must be able to accurately
capture _(1)_ the structure at all relevant scales down to the pressure
smoothing scale that characterizes the small-scale cutoff in power due to the
IGM thermal history (Kulkarni et al., 2015; Puchwein et al., 2023; Doughty et
al., 2023) and _(2)_ the inhomogeneous reionization and post-reionization
thermal state (including the temperature at mean density and the slope of the
temperature density relation). However, given the known degeneracies between
reionization timing, the amount of photoheating to expect, and the UVB
amplitude compounded with the overall uncertainty regarding the reionization
timeline, it is important to consider a wider range of reionization histories,
simulated using an appropriate code, and at an appropriate resolution to
capture the pressure smoothing.

Here, we begin this work by using the semi-numerical code AMBER to create a
variety of inhomogeneous reionization histories that are compatible with the
best available constraints on reionization: the _Planck_ CMB results and
model-independent dark pixel fractions. Using a small set of sample histories,
and varying the heat injection within the cosmological hydrodynamics code Nyx,
we demonstrate how the parameters of the reionization histories impact gas
within the IGM and lead to degenerate conclusions results for a preferred
ionization history when examined using a simplistic “$\chi$-by-eye” evaluation
technique. In Section 2 we begin by describing our methods, including the
hydrodynamical simulation setup, modeling of inhomogeneous reionization, and
justification for our chosen test $x_{\mathrm{HI,v}}(z)$ histories. In Section
3 we first explore how the histories affect the temperature and density
evolution of gas in the IGM, including the pressure smoothing. We then delve
into the effect on the observable Ly$\alpha$ opacities and compare them to a
subset of existing observations. We conclude the work by synthesizing the
results in Section 4 and presenting a summary in Section 5.

##  2 Simulations

In order to model reionization, we must include separate components to evolve
the formation of structure as well as the actual ionization process. The
reionization of neutral hydrogen is a prolonged process, with the bulk
occurring below $z=9$. It is induced in overdense regions earlier due to the
greater abundance of ionizing photons there, resulting in a patchy field of H
II. This generates fluctuations in temperature as well, both impacting the
opacity in the Ly$\alpha$ forest. In this section we describe the method for
modeling this patchy reionization, and the details of the histories we choose.

###  2.1 Modeling an inhomogeneous reionization

Reionization can be modeled through the inclusion of _in-situ_ radiative
transfer (e.g. Rosdahl et al., 2018, 2022; Kannan et al., 2022), or by post-
processing of the volume (e.g. Puchwein et al., 2023). More cheaply, semi-
analytical methods may be used to generate realistic ionization fields, as in
Battaglia et al. (2013), or tools like 21cmfast (Mesinger et al., 2011) and
FASTPM (Feng et al., 2016).

To generate the reionization field for our simulations we elect to use the
abundance-matching model AMBER (Trac et al., 2022), which creates a gridded
reionization field by matching a selected mass-weighted reionization history
to a rank-ordered pseudo-radiation field. This matching is based on the
assumption that grid cells encountering a stronger radiation field will
reionize first. The first step is determining a reionization history, which is
parametrized in terms of a midpoint $z_{\mathrm{mid}}$, duration $\Delta z$,
and asymmetry $A_{z}$. The midpoint is the redshift where the mass-weighted
ionized fraction $x_{\mathrm{HII,m}}$ equals 0.5. “Early” and “late” points in
the history are defined as $x_{\mathrm{HII,m}}=0.05$ and
$x_{\mathrm{HII,m}}=0.95$, respectively, and calculated as

$z_{\mathrm{early}}=z_{\mathrm{mid}}+\frac{\Delta
zA_{z}}{1+A_{z}}\>\>;\>\>z_{\mathrm{late}}=z_{\mathrm{mid}}-\frac{\Delta
zA_{z}}{1+A_{z}}$ (2)

Using these definitions, the duration is equal to the difference between
$z_{\mathrm{early}}$ and $z_{\mathrm{late}}$, and the asymmetry is the ratio
$(z_{\mathrm{early}}-z_{\mathrm{mid}})/(z_{\mathrm{mid}}-z_{\mathrm{late}})$.
A Weibull function (Weibull, 1951) is used to interpolate between the three
fixed points $z_{\mathrm{mid}}$, $z_{\mathrm{early}}$, and $z_{\mathrm{late}}$
and arrive at a fully defined history.

Figure 1: The volume-weighted neutral fraction (_top_) and PDFs of the
reionization redshift (_bottom_) for four $z_{\mathrm{mid}}=6.9$ example
reionization histories from AMBER, demonstrating the impact of varying the
duration ($\Delta z$) and asymmetry ($A_{z}$) parameters. With all other
parameters held constant, the (1) duration dictates the peak fraction being
ionized, and (2) asymmetry dictates whether the distribution of
$z_{\mathrm{re}}$ is symmetric about the peak. In practice, an asymmetric
model leads to a “tail” in $x_{\mathrm{HI,v}}$ out to lower $z$. While not
shown here, the mass-weighted midpoint shifts the ($x_{\mathrm{HI,v}}$,$z$)
intersection point horizontally.

In the base AMBER code, the density field is acquired by first evolving a
Gaussian random field (adjusted to match the cosmology-determined power
spectrum) using linear perturbation theory (LPT; Bouchet et al., 1995;
Scoccimarro, 1998), with the option to use either the Zel’dovich approximation
(Zel’dovich, 1970) or the more accurate second order LPT (2LPT). Assignment of
the particles to a grid is performed using two staggered, interlaced meshes to
improve the accuracy of the overdensity fields, and avoid effects such as
aliasing (Hockney & Eastwood, 1988). We do not use the built-in functionality
of AMBER to generate initial conditions from a Gaussian random field. Instead
we extract the $z=0$ linear density field from our initial conditions
generator (see Section 2.2.1 for more details) and use AMBER to evolve it
using 2LPT to $z=z_{\mathrm{mid}}$.

From the overdensity field, the excursion set formalism is used to model the
collapsed mass fraction within the simulation volume, and extended Press-
Schechter theory is used to derive a halo mass function and the resultant halo
density field. To extend this to the radiation field, the background is
defined as usual

$\Gamma_{\mathrm{HI}}(\textbf{x)}=4\pi\int^{\infty}_{\nu_{\mathrm{HI}}}\frac{J_{\nu}(\nu,\textbf{x})}{h\nu}\sigma_{\mathrm{HI}}(\nu)\mathrm{d}\nu$
(3)

where $\sigma_{\mathrm{HI}}$ is the cross-section for hydrogen
photoionization. $J_{\mathrm{\nu}}$ is the specific intensity, itself given by

$J_{\mathrm{\nu}}(\nu,\textbf{x})=\int\frac{S_{\nu}(\nu,\textbf{x'})}{4\pi|\textbf{x}-\textbf{x}^{\prime}|^{2}}\exp\left(-\frac{|\textbf{x}-\textbf{x}^{\prime}|}{\lambda_{\mathrm{mfp}}}\right)\mathrm{d^{3}x^{\prime}}$
(4)

with source function $S_{\nu}$, and the mean free path of ionizing photons
$\lambda_{\mathrm{mfp}}$. Recent studies of the mean free path suggest
significant evolution between $5<z<6$ (Becker et al., 2021; Gaikwad et al.,
2023; Zhu et al., 2023), but for the majority of the reionization process that
occurs at $z>6$, this value is uncertain. The mean free path
$\lambda_{\mathrm{mfp}}$ is defined using the standard theoretical definition,
the distance over which the transmitted flux fraction drops to $1/e$. Given
that the mean free path mainly affects the topology of reionization rather
than the history, we simply use a constant 3 $h^{-1}$Mpc for all our
histories, and leave a more detailed consideration for future work.

The source function,

$S_{\nu}(\nu,\textbf{x})=f_{\mathrm{esc}}\epsilon_{\nu}(\nu)\rho_{\mathrm{SFR}}(\textbf{x}),$
(5)

is dependent on several astrophysical variables such as the radiation escape
fraction $f_{\mathrm{esc}}$, the radiative energy per unit star formation rate
(SFR) per $\nu$. The star formation rate density is in turn given by

$\rho_{\mathrm{SFR}}(\textbf{x})=\frac{f_{\mathrm{star}}\rho_{\mathrm{halo}}(\textbf{x})}{\tau_{\mathrm{star}}}$
(6)

where $\rho_{\mathrm{halo}}$ is the halo mass density field, and
$f_{\mathrm{star}}$ and $\tau_{\mathrm{star}}$ are the star formation
efficiency and timescale, respectively. Equations 5 and 6 are not themselves a
part of the reionization parametrization, and the variables therein are simply
tuned to produce a history with the desired values.

All $N$ cells are assigned a rank order $n$ based on the values of the
radiation field, and the mass fraction is calculated for each rank

$\frac{\sum^{n}_{m=1}1+\delta_{m}(z)}{\sum^{N}_{m=1}1+\delta_{m}(z)}=\bar{x}_{i}(z)$
(7)

where $\bar{x}_{i}(z)$ is the ionized fraction at redshift $z$ and
$\delta_{m}$ is the overdensity of the $m$th cell at $z$. For each mass
fraction, the reionization redshift can be determined by passing the fraction
$\bar{x_{i}}$ to the inverted Weibull function. This uses the density field at
$z_{\mathrm{mid}}$, but a few iterative steps are performed to adjust these
values to the field at $z=z_{\mathrm{re}}$ assuming linear growth of the
overdensities from the reionization midpoint.

We show some example histories in the top panel of Figure 1, for a constant
reionization midpoint of $z=6.9$ and two choices each for asymmetry and
duration. The midpoint plays the primary role in positioning the reionization
history, here producing a volume-weighted neutral fraction of approximately 75
per cent at $z=6.9$. For a constant asymmetry parameter of $A_{\mathrm{z}}=7$
(dashed lines), a longer duration of $\Delta z=4$ compared to 2 increases the
time between the starting and ending redshifts, so reduces the steepness of
the history. For a constant duration $\Delta z=2$ (thin lines), a greater
asymmetry steepens the history profile and in particular removes any lingering
neutral regions at late times.

To visualize this another way, we show the normalized PDFs of the redshifts at
which the cells reionize in these example models (bottom panel Figure 1).
Here, the distributions for the two symmetric models ($A_{\mathrm{z}}=1$) peak
at the center of an approximately normal distribution, with a roughly equal
numbers of cells reionizing above and below the peak. Their peaks also contain
only a small fraction of the total cells, $<10$ percent. The asymmetric
distributions on the other hand have many more occurring at
$z>z_{\mathrm{peak}}$ and $\sim 30$ percent in the peak bin itself. The peaks
do not occur at the mass-weighted midpoint, but rather at a lower redshift. It
is apparent that the models with asymmetric histories will experience an
abrupt ionizing event, coupled with a heat injection that may result in a
relatively uniform temperature field and a more coherent pressure smoothing
scale.

Table 1: A summary of the simulations used in Section 3. The $z_{\mathrm{re}}$ fields have a coarse grid of $128^{3}$ and are generated using AMBER to abundance match the radiation field to an arbitrary reionization history, as described in Section 2.1. $z_{\mathrm{early}}$, $z_{\mathrm{mid}}$, and $z_{\mathrm{late}}$ are the points that define the arbitrary history in terms of mass-weighted H I fractions. All of the runs use $L_{\mathrm{box}}=20$ $h^{-1}$Mpc with a fixed grid resolution $2048^{3}$ cells, corresponding to spatial resolution $\Delta x=10$ $h^{-1}$kpc. The naming convention distinguishes between midpoints (Late $z\sim 6$, Middle $z\sim 7$, and Early $z\sim 8$, durations (Short $\Delta z\leq 3$ or Long $\Delta z>3$), and asymmetry (Asymmetric $A_{z}>1$ or Symmetric $A_{z}=1$). The heat injection is further indicated using Cold (10000 K) or Hot (20000 K). name |  $z_{\mathrm{early}}$ ($x_{\mathrm{HI,m}}=0.95$) |  $z_{\mathrm{mid}}$ ($x_{\mathrm{HI,m}}=0.5$) |  $z_{\mathrm{late}}$ ($x_{\mathrm{HI,m}}=0.05$) |  $z_{\mathrm{end}}$ ($x_{\mathrm{HI,m}}=0.0$) | $\Delta z$ | $A_{\mathrm{z}}$ |  $\Delta T_{\mathrm{re}}$ (K) | $\tau_{\mathrm{e}}$
---|---|---|---|---|---|---|---|---
LateShortAsym | 7.9 | 6.1 | 5.9 | 5.9 | 2 | 8 | 1e4, 2e4 | 0.0398
MidShortAsym | 8.7 | 6.9 | 6.7 | 6.7 | 2 | 8 | 1e4, 2e4 | 0.0475
MidLongAsym | 21.2 | 7.3 | 6.2 | 6.2 | 15 | 13 | 1e4, 2e4 | 0.0596
EarlyLongAsym | 21.8 | 7.9 | 6.8 | 6.8 | 15 | 13 | 1e4, 2e4 | 0.0661
EarlyShortSym | 9.7 | 8.2 | 6.7 | 5.4 | 3 | 1 | 1e4, 2e4 | 0.0549
EarlyLongSym | 10.2 | 8.2 | 6.2 | 4.4 | 4 | 1 | 1e4, 2e4 | 0.0529
Figure 2: _Top:_ Reionization histories for our selected models generated
using AMBER. The dark and light gray shaded regions indicate the constraints
from Planck Collaboration et al. (2020) using the FlexKnot model. Overplotted
are observationally-inferred values from McGreer et al. (2015) (purple), Zhu
et al. (2022) (red), Jin et al. (2023) (salmon), Davies et al. (2018) (blue),
Wang et al. (2020) (orange), and Yang et al. (2020) (green). _Middle:_ PDFs of
the reionization redshifts for each model. Visually, the most major
differences are between symmetric Sym and asymmetric Asym models, with the
latter having high peaks. _Bottom:_ The Thompson scattering optical depths for
the reionization models in Figure 2, calculated according to Equation 8. The
dark and light gray shaded regions indicate the $1$ and $2\sigma$ limits from
Planck Collaboration et al. (2020), and the colored ticks on the right y-axis
show the full predicted $\tau_{e}$ when integrated out to the surface of last
scattering. Broadly, a model that is more reionized earlier has a higher
$\tau_{e}$, since there are more free electrons in the IGM to contribute to
scattering.

###  2.2 Selecting reionization histories

Figure 2 shows our selected reionization histories as the volume-weighted H I
neutral fraction ($x_{\mathrm{HI,v}}$), overplotted with select observational
constraints,111We note again that the model parameters describe the _mass-
weighted_ reionization history, while the history in Figure 2 is volume-
weighted. including Planck Collaboration et al. (2020) 2$\sigma$ limits
indicated in gray. Given our goal of probing the parameter space suggested by
the strongest constraints (taken here to be _Planck_ FlexKnot and McGreer et
al., 2015), we select several combinations of ($z_{\mathrm{mid}}$, $\Delta z$,
$A_{z}$) to span most of this swath; the details of these models are presented
in Table 1. For the naming convention, we use Early, Mid, and Late to indicate
early ($z\sim 8$), middling ($z\sim 7$), and late ($z\sim 6$) mass-weighted
midpoints. For duration, we use Short ($\Delta z<4$) or Long ($\Delta z\geq
4$), a somewhat arbitrary designation. For asymmetry we use Asym ($A_{z}>1$)
and Sym ($A_{z}=1$). Lastly, we use hot ($\Delta T_{\mathrm{re}}=20,000$ K)
and cold ($\Delta T_{\mathrm{re}}=10,000$ K) to refer to the amount of heat
injection as a region reionizes (see Section 2.2.1 for details of this
implementation).

We examine values within a reasonable parameter space
$6.1<z_{\mathrm{mid}}<8.2$, $2<\Delta z<15$, and $1<A_{\mathrm{z}}<13$. Within
this space, there are many possible histories, but we manually select six that
extend from extremely late, short histories (upper left of panel _a_) to
earlier, and more symmetric ones (lower right of panel _a_). Most of our
selected histories do not violate the $2\sigma$ limits of FlexKnot with the
exception of the EarlyLongSym model, which approximately matches the McGreer
et al. (2015) limits. However, the true endpoints $z_{\mathrm{end}}$ where the
volume-weighted neutral fraction reaches zero sometimes occur quite late,
$z_{\mathrm{end}}=5.4$ and 4.4 for EarlyShortSym and EarlyLongSym,
respectively. EarlyLongSym in particular is quite late to reionize fully, but
we retain it as an extreme example history.

Besides McGreer et al. (2015) in purple, the other colored points in Figure 2
indicate inferences from Ly$\alpha$ and Ly$\beta$ dark gaps (Zhu et al., 2022;
Jin et al., 2023). For $z\geq 7$, we show select results from analyses of
quasar damping wings (Davies et al., 2018; Bañados et al., 2018; Greig et al.,
2019). Other methods of inferring this history, whose results we do not
include here, include the Ly$\alpha$ equivalent width distribution (e.g. Mason
et al., 2018; Hoag et al., 2019; Bruton et al., 2023) and Ly$\alpha$
clustering (e.g. Ouchi et al., 2010). Many of these $z>7$ observationally-
inferred $x_{\mathrm{HI,v}}$ values are in some contention with one another,
related to systematic discrepancies between the methods but also potentially
due to the patchy topology of reionization and/or the biased nature of certain
environments (e.g. around quasars). However, these other results are broadly
in support of the _Planck_ limits, and the histories we have chosen.

In the center panel of Figure 2 we show the normalized PDFs of the redshifts
at which the cells reionized. The four asymmetric models (Asym) are highly
peaked, so it is expected that these models will resemble a “flash” or
instantaneous reionization history, where a large amount of heat is injected
to the IGM at a single redshift. The symmetric models (Sym) show a “normal”
distribution of $z_{\mathrm{re}}$, which will lead to a gradual injection of
heat in the simulations and a wide distribution of pressure smoothing lengths
in the IGM.

Using the reionization fields generated in AMBER, we calculate the optical
depth to Thompson scattering of the cosmic microwave background, $\tau_{e}$.
To accomplish this, we use

$\tau_{\mathrm{e}}=\int_{0}^{z_{\mathrm{ls}}}\frac{c\;\sigma_{\mathrm{T}}\;n_{e}\left(z\right)}{\left(1+z\right)\;H\left(z\right)}\mathrm{d}z$
(8)

where $z_{\mathrm{ls}}$ is the redshift of last scattering ($z\sim 1100$),
$\sigma_{T}$ is the Thompson scattering cross-section, $n_{e}$ is the electron
density, and $H(z)$ is the Hubble parameter at redshift $z$. For the electron
density, we assume contributions from both H II and He II for $z>3$, and add
contributions from He III for $z\leq 3$, assuming that quasars have reionized
all the He II by this point in time.

With this approximation we arrive at the curves in the bottom panel of Figure
2. The lowest and highest $\tau_{\mathrm{e}}$ correspond to models
LateShortAsym and EarlyLongAsym, whose largest differences are in their
midpoints and durations. LateShortAsym initiates later and is much shorter,
leading to a smaller integrated electron density along the line of sight
compared to the EarlyLongAsym model, which initiates quite early. The two
models closest to the _Planck_ result are EarlyShortSym and EarlyLongSym,
which bracket the measurement of $\tau_{e}\approx 0.054$.

####  2.2.1 Nyx

For our simulations we use the open source code222http://github.com/AMReX-
Astro/Nyx Nyx, a highly parallel, adaptive mesh, finite-volume N-body
compressible hydrodynamics solver for cosmological simulations (Almgren et
al., 2013; Sexton et al., 2021). Nyx scales well on CPU- and GPU-based
machines, and has been used extensively in studies of the IGM and Ly$\alpha$
forest (see e.g. Lukić et al., 2015; Oñorbe et al., 2019; Chabanier et al.,
2023; Wolfson et al., 2023b; Jacobus et al., 2023). We use a domain of 20
$h^{-1}$Mpc with spatial resolution 10 $h^{-1}$kpc , which is higher than
typically used in IGM studies, as we have found this to be necessary for
convergence at $z>5$ given our numerical setup (Doughty et al., 2023).

We generate the transfer functions using CAMB (Lewis et al., 2000), and the
initial conditions using the code CICASS (O’Leary & McQuinn, 2012), including
a $\langle\Delta v\rangle=30$ km/s streaming velocity between baryons and dark
matter at recombination. We have found that inclusion of a non-zero streaming
velocity is found to create differences on the few per cent level in
Ly$\alpha$ power, in particular leading to slightly reduced small scale power
with respect to a $\langle\Delta v\rangle=0$ km/s simulation, also generated
via CICASS. This difference increases with resolution, but down to at least
$\sim\Delta x=2$ $h^{-1}$kpc is subdominant to other effects. Initial
conditions are generated for $z=200$, and used to inform the overdensity field
in AMBER as described in Section 2.1. We assume a $\Lambda$CDM cosmology
consistent with results from Planck (Planck Collaboration et al., 2020),
$\left(\Omega_{M},\Omega_{\Lambda},\Omega_{b},h,X_{H}\right)=$ (0.315, 0.685,
0.049, 0.675, 0.76).

Nyx excludes galaxy formation physics and we omit any type of feedback
implementation, which should have a minimal effect on the high-$z$ Ly$\alpha$
forest (Viel et al., 2013b), so the main impacts on the gas come from
structure formation and the associated shock heating, and the reionization
event itself. To model inhomogeneous reionization in Nyx, we use the coarse-
grid reionization fields generated from AMBER, which are mapped onto the finer
hydrodynamics grid. Once a hydrodynamics grid cell reaches its reionization
redshift, a user-defined amount of heat is injected dependent on the current
gas temperature and the self-consistently calculated hydrogen neutral
fraction.

$T^{\mathrm{post-reion}}=\left(x_{\mathrm{HI}}^{\mathrm{pre-
reion}}\right)\left(\Delta T_{\mathrm{re}}-T^{\mathrm{pre-
reion}}\right)+T^{\mathrm{pre-reion}}$ (9)

where $T^{\mathrm{pre-reion}}$ and $x_{\mathrm{HI}}^{\mathrm{pre-reion}}$ are
the pre-reionization temperature and neutral fraction for a given cell. It can
be heated by a possible maximum of $\Delta T_{\mathrm{re}}$, either 10000
(cold) or 20000 K (hot) in this work. We output simulation snapshots in bins
of $\Delta z=0.1$ from $5.0\leq z\leq 6.0$, and use these for all the analyses
presented here.

We generate physical and Ly$\alpha$ optical depth skewers cast along all axes
of the simulation, and in the positive or negative direction. The optical
depths are calculated assuming that they are optically thin and in
thermodynamic equilibrium, and do not use the self-consistent neutral fraction
from Nyx. In theory, this could mean we would retrieve a non-zero flux from a
cell that has not yet been reionized, which is not realistic. To address this,
we set the H I photoionization rate to an extremely low value in cells that
have not reionized by the snapshot redshift, thereby preventing any undue
transmission.

Figure 3: 50 km/s slices of a 20 $h^{-1}$Mpc, $2048^{3}$ simulation with a
$128^{3}$ reionization field. The left column shows the gas overdensity field,
while the center and right columns show the temperature and reionization
redshift fields, respectively. Since the simulations have initial conditions
with the same seed, the overdensity structures are the same; the differences
come from the reionization histories. The top row shows the $\Delta
T_{\mathrm{re}}$ versions of LongShortAsym and EarlyShortSym (red and orange
curves in Figure 2), with LongShortAsym being rapid with a late midpoint, and
EarlyShortSym being prolonged with an early one. The rapid LongShortAsym
history gives a homogeneous-looking $z_{\mathrm{re}}$ field, and a
correspondingly featureless temperature field.

##  3 Results

In this work, we are attempting to examine the trends in several Ly$\alpha$
opacity-based metrics of reionization, the 1D flux power spectrum,
autocorrelation function, and effective optical depth distribution. To
facilitate an understanding of those trends, we first characterize the impact
of the differing reionization histories on the density and temperature fields,
phase-space distribution of gas, and the temperature at mean density.

###  3.1 Measures of gas in the simulation

####  3.1.1 Slices

In Figure 3 we show 50 km/s slices in overdensity and temperature at $z=5.5$,
and the reionization field for two reionization models: on the top row is the
late, rapid, asymmetric, and hot (LateShortAsym) model, while on the bottom is
an early, rapid, symmetric history, hot model (EarlyShortSym). For both
models, overdense structures are visible in the left column as lighter colors,
revealing the filamentary structure of the cosmic web. The temperature panels
(center column) show hotter temperatures ($>10^{4.5}$ K) co-located with the
filaments, which trace collapsed structures as they have been shock heated to
their virial temperatures. Due to the identical seeds used in the initial
conditions, the same overdensity structure is naturally present in both
models.

Beyond these commonalities, the differing natures of the two histories become
apparent. The abrupt LateShortAsym model shows a very uniform reionization
field (right column), with only the _most_ overdense regions ionizing at
barely over $z=6$, while EarlyShortSym has a much wider dynamic range in
$z_{\mathrm{re}}$ values. Accordingly, LateShortAsym has a relatively
featureless temperature field in the underdense IGM, nearly all cells hovering
around $\log\;(T\mathrm{/K})\sim 4.2$. There are minor excursions to lower
temperatures ($\sim 10^{4.0}$ K) around the largest filaments, for example
near (8, 8) $h^{-1}$Mpc , but they are quite minimal and difficult to pick
out. For EarlyShortSym on the other hand, the regions surrounding the
filaments are those that reionized first, indicated by yellow coloring in the
reionization field in the rightmost column. Since those regions became
reionized first, they have had the most time to cool, and have achieved much
lower temperatures than the regions which were more recently ionized. This
slice also captures a few completely un-ionized regions in the EarlyShortSym
simulation, which doesn’t fully complete (in a volume-weighted sense) until
$z=5.4$; these are the glaringly cold purple/blue cells visible in the lower
center panel of Figure 3. Since they have not been photoheated by reionization
and are physically distant from structure formation-induced shock heating,
they have been permitted to cool adiabatically since recombination, reaching
temperatures of $<100$ K.

####  3.1.2 Temperature density relation

Figure 4:  Temperature density relations of all 12 models, grouped based on
visual similarity into rows, from top to bottom: Short asymmetric (ShortAsym),
long asymmetric (LongAsym), and symmetric (Sym). From this visualization, it
is apparent that an earlier reionization midpoint allows the gas to further
relax, and move away from a flat isothermal relation (e.g. LateShortAsym
versus MidLongAsym). Asymmetry causes a large injection of heat at later
times, more resembling an instantaneous reionization, while a symmetric early
history permits more gas to have relaxed into the expected temperature-density
relation.

It is expected that gas in the IGM will settle into a relatively tight
relation between temperature and density after reionization has completed (Hui
& Gnedin, 1997)

$T(\Delta)=T_{0}\Delta^{\gamma-1}$ (10)

where $\Delta$ is the overdensity, $\rho/\langle\rho\rangle$, and $\gamma$ is
the adiabatic index, found to be $\sim 1.6$ in a relaxed medium. For an
inhomogeneous reionization, however, the amount of time since reionization
occurred will vary between different regions in the simulation, meaning that
only a fraction of the gas will obey the relation, and instead have a wider
overall spread in properties (Oñorbe et al., 2019). To examine the impact of
the reionization parameters on the distribution, we plot the temperature-
density distributions for all models at $z=5.5$ in Figure 4.

For the short asymmetric histories (ShortAsym, top row), the bulk of the gas
occupies a relatively narrow range of temperatures for a given overdensity,
although this is not yet the relaxed temperature-density relation of Equation
10. All else being equal, models with higher heat injection have a flatter
distribution, resembling that obtained when assuming an instantaneous
reionization, whereas the colder models have moved closer to their eventual
log-log distribution. Long, asymmetric histories (LongAsym; second row) tend
to have gas split between the final $\Delta-T$ relation with steeper $\gamma$
and a flatter, unsettled distribution more resembling the ShortAsym models.
These are visible as two dark red lines at $\log\Delta<1$.

Symmetric models (Sym; bottom row) create two distinguishing features in their
distributions. First, there is gas remaining below the temperature density
relation, whereas for all the other models this region of the phase space is
essentially empty. This gas, still mostly at IGM overdensities, is removed by
the reionization process, but these reionization histories end the latest,
both after $z=5.5$.333There is a “hook” of gas that first increases
logarithmically from $(\log\Delta,\log(T/\mathrm{K})=(1.0,3.5)$ before
becoming flat out to $\log\Delta>1.4$. This is likely shock heated gas that in
the other simulations has already been removed from this phase by additional
heating reionization process. Further, these models resemble the LongAsym ones
in that they have a wider range of temperatures below $\log\Delta=1$ than the
short histories. However, more of their gas has settled into the final log-log
temperature-density relation, and they are lacking material in the upper dark
red line present in the LongAsym histories. As seen in Figure 2, it seems that
asymmetric histories contribute to the majority of the gas being reionized in
a single, impulsive event, causing a higher temperature and flatter $\rho-T$
relation prior to the IGM’s thermodynamic relaxation.

####  3.1.3 Temperature at mean density

Figure 5: The temperature at mean density for our simulation runs, overplotted
with values derived from $z\geq 5$ observations (Garzilli et al., 2017; Boera
et al., 2019; Walther et al., 2019; Gaikwad et al., 2020). The hot models are
indicated with solid lines, while dashed lines are for the cold ones. The
model trajectories span a wide range of temperatures, but are not well-
constrained by the observations, which are themselves in contention with one
another.

An important part of Equation 10 is the temperature at mean density, $T_{0}$,
as it is reflective a combination of heating and cooling events impacting the
majority of gas in the Universe. We calculate the temperature at mean density
from $5<z<6$ for each of our simulation runs, and show them in Figure 5. These
are calculated from the full snapshot data outputted at each $\Delta z=0.1$
interval redshift, within the overdensity range $0.95<\Delta<1.05$. As a point
of comparison, we overplot observationally-derived temperatures from several
studies (Garzilli et al., 2017; Boera et al., 2019; Walther et al., 2019;
Gaikwad et al., 2020). Since a nonzero Ly$\alpha$ signal can’t be obtained
from the un-ionized IGM, values derived from observations are probing
reionized regions, so we also restrict contributions to cells that have been
ionized by masking out those with $z_{\mathrm{re}}<z_{i}$. In practice, this
is only relevant for models with $z_{\mathrm{end}}<6.0$.

First, we find that all of the cold models (dashed lines) have very similar
temperatures, $\sim$8000 K, and do not show much cooling over $5<z<6$, instead
being relatively constant. The hot models all show greater temperature
evolution, and the inter-model variation covers a much wider range of
temperature. We find that our extreme LateShortAsym model creates a $T_{0}$
that is too high by a factor of $\sim 1.5$, lying above even the relatively
high derived temperatures in Gaikwad et al. (2020). The other hot models
MidShortAsym, MidLongAsym, and EarlyLongAsym all cut through the Gaikwad et
al. (2020) datapoints, with MidShortAsym and MidLongAsym overshooting the
$T_{0}<10000$ K measurements from other studies. EarlyShortSym and
EarlyLongSym lie the closest to the cold model trajectories, approximately
passing between the low $T_{0}$’s from Garzilli et al. (2017) and Walther et
al. (2019) and the high values from Gaikwad et al. (2020), and matching the
$z\sim 5$ data point from Boera et al. (2019).

It is interesting to note the cases where different reionization histories
have induced very similar $T_{0}$ trajectories. The MidShortAsym and
MidLongAsym histories have nearly identical $T_{0}$ evolution, with the
MidLongAsym history being more ionized at $z>6.6$ and completing about $\Delta
z=0.5$ later, while having a similar midpoint. The two Sym histories are even
more similar. Evidently, despite EarlyLongSym’s extremely late end, the
continued heat injection at late times is from such a small number of cells
($x_{\mathrm{HI,v}}(z\sim 5.4)=0.05$), that it doesn’t matter for the average
temperature. With the exception of the hot EarlyLongAsym model, the similar
histories are ones with similar mass-weighted midpoints.

####  3.1.4 Pressure smoothing scale

The instantaneous temperature of the IGM is important for understanding its
structure, but in actuality it takes time for the gas to thermodynamically
respond to a heating event like reionization. To quantify the physical
response at a given redshift, we examine the pressure smoothing scale, also
called the filtering scale (Gnedin & Hui, 1998), describing the physical scale
below which there is a cutoff in smaller structures. This is due to pressure
smoothing induced by the instantaneous temperature and also the thermal
history.

To do this, we first calculate the 3D matter density power spectra in all our
redshift bins, restricting it to the gas at overdensities $0.2<\Delta<0.4$ to
isolate the underdense gas contributing to the Ly$\alpha$ forest. We then
compare the 3D power $P_{\mathrm{rad}}(z)$ to that in an otherwise identical
simulation that has been allowed to adiabatically cool with no reionization
event, and thus excludes any pressure smoothing effect, $P_{\mathrm{ad}}(z)$.
The expected relationship is characterized as

$P_{\mathrm{rad}}(k)=N\cdot
P_{\mathrm{ad}}(k)\cdot\exp(-k^{2}/k_{\mathrm{ps}}^{2})^{2}$ (11)

as defined in Puchwein et al. (2023) (and matching our previous treatment in
Doughty et al., 2023). We then fit for a normalization parameter $N$ and the
pressure smoothing wavenumber $k_{\mathrm{ps}}$.

We plot the normalized ratios $P_{\mathrm{rad}}/P_{\mathrm{ad}}$ for both of
the ShortAsym models as an example in Figure 6. We include both the hot and
cold versions (solid and dashed lines) and at $z=5.5$ and $z=5.0$ (less and
more saturated in color, respectively). Models that are more removed from the
$\mathrm{ratio}=1$ line have more pressure smoothing, i.e. they have less
power on small scales and so fall off farther to the left.

Figure 6: The normalized ratios of Equation 11 showing the loss of small scale
power due to heat-induced pressure smoothing, using the ShortAsym models as an
example. The lines with lower (higher) saturation are taken at $z=5.5$
($z=5.0$), showing how the smoothing effect is lower at higher $z$.
LateShortAsym and MidShortAsym are identical except for their reionization
midpoint, and the earlier midpoint in MidShortAsym leads to increased
smoothing when considered at the same redshift. Hot models (solid lines) also
are more smoothed than cold ones (dashed). Figure 7: _Top and middle_ : The
evolution of the pressure smoothing scale (top) and the measured adiabatic
index (middle) with time since the mass-weighted reionization midpoint. The
points indicate times corresponding to each of our redshift snapshots. Since
both are related to thermodynamic relaxation, they both increase with time.
_Bottom_ : The correlation between the measured pressure smoothing scale and
the adiabatic index, measured in the same simulation snapshots. Both values
increase with time, and so are strongly correlated with one another for a
given model. Many of the model trajectories overlie one another, but there is
an obvious separation between the hot and cold models (solid and dashed lines,
respectively). Figure 8: Power spectra of 1D Ly$\alpha$ flux, with data from
Boera et al. (2019) at $z=5$. All hot models (left panel) have extremely
similar values and generally lie quite close to the observed values, while the
cold models (right panel) differ from each other, and demonstrate excessive
small scale power. For the cold models, the order of increasing high $k$ power
matches the order of increasing temperature at mean density. This is likely
because the temperature is an indicator of time since the reionization
midpoint, and a more recent midpoint results in a smaller pressure smoothing
scale.

There are trends here which can serve as a reality check on the method, by
confirming our intuition about how the pressure smoothing scale should evolve
with time and temperature. First, the hot versions of each reionization
history peel away from the $\mathrm{ratio}=1$ line at larger physical scales
than their cold counterparts, which makes sense given the increased speed of
sound in hotter gas, and the consequently increased ability for the gas to
expand post-reionization. The history with the $z\sim 6$ midpoint
LateShortAsym shows less pressure smoothing than the earlier one, which is
sensible given that it has had less time to physically respond to the heat
injection. Lastly, of course the pressure smoothing is larger for any given
model at $z=5$ than at $z=5.5$, again since more time has passed since the
thermal injection event.

Given that both the pressure smoothing scale and the slope of the $\rho-T$
relation $\gamma$ are related to the time elapsed post-reionization, we
examine the relationship between these two quantities. Fitting the $\rho-T$
distributions in Figure 4 we extract the adiabatic index $\gamma$ and plot
both it and our measured pressure smoothing scale $\lambda_{\mathrm{ps}}$
(equal to $2\pi/k_{\mathrm{ps}}$) against the time elapsed since the mass-
weighted reionization midpoint, $\Delta t_{\mathrm{re}}$. The results are
plotted in Figure 7.

The fitted pressure smoothing scale curves (top panel) start where $\Delta
t_{\mathrm{re}}$ is the time difference between the mass-weighted midpoint and
the time at $z=6$, and increase approximately linearly with time.
Extrapolating the curves to $\Delta t_{\mathrm{re}}=0$, they radiate from some
initial smoothing scale $\lesssim 5$ $h^{-1}$kpc. Comparing between hot and
cold versions of the same reionization histories, the hot models have steeper
slopes and have larger values for a given $\Delta t_{\mathrm{re}}$. The
general linear trend is supportive of the simple expansion model shown in
Puchwein et al. (2023).

The time evolution of the adiabatic indices is a bit messier, as while they
increase with time the rate of change slows and there is a turnover at higher
$\gamma$. This is particularly evident in the cold models, which also start at
higher values than their hot counterparts, a logical result given that it does
not take as long for a model with low heat injection to thermodynamically
relax into the expected power-law relation. Unlike the pressure smoothing
scale, these curves do not all emanate from the same point; however, we should
expect that the “initial” $\gamma$ would be close to one, particularly for the
Asym models with their quasi-instantaneous reionizations and resultant nearly
isothermal temperature distributions.

We also directly compare $\lambda_{\mathrm{ps}}$ and $\gamma$ measured at the
same redshifts in the bottom panel of Figure 7, with the hot models indicated
with circles and the cold models with squares. Of course, the values both
increase with decreasing redshift (lower left to upper right within the
panel), and the hot and cold models occupy different regions of the plot, with
cold models having higher $\gamma$’s for a given $\lambda_{\mathrm{ps}}$. Many
of the $\lambda_{\mathrm{ps}}-\gamma$ curves even lie on top of one another,
particularly for the hot models, even for significantly different reionization
histories. Based on these results, in the Ly$\alpha$ metrics we should expect
to see evidence of more small scale structure in the cold models, even if they
are more thermodynamically relaxed.

###  3.2 Metrics of the Ly$\alpha$ forest

####  3.2.1 1D Ly$\alpha$ forest power spectrum

The 1D Ly$\alpha$ flux power spectrum characterizes the fluctuations in the
flux field as a function of the wavenumber $k$

$P_{\mathrm{Ly\alpha}}(k)=\int d^{3}\bm{x}\,\xi(x)\,e^{-i\,\bm{k}\,\bm{x}}$
(12)

where $\xi(\bm{x})$ is the two-point correlation function, or autocorrelation
function

$\xi(x)=\langle\delta_{F}(\bm{x}_{1})\,\delta_{F}(\bm{x}_{2})\rangle$ (13)

with $x=\lvert\,\bm{x_{1}}-\bm{x_{2}}\,\rvert$ and $\delta_{F}=F/\langle
F\rangle-1$. Greater power indicates greater fluctuations in the flux contrast
on that spatial scale, and thus a larger amount of structure with size
$\lambda=2\pi/k$. We calculate the power for our models at $z=5$, and plot
them in Figure 8 alongside the measurements from Boera et al. (2019) who were
able to probe $-2.2\leq\log$$(k/(\mathrm{s\;km}^{-1}))$$\leq-0.7$, reaching an
unprecedented high $k$ mode achievable through use of high S/N spectra. All
models have had their optical depths adjusted to produce a mean flux of
$\langle F\rangle=0.18$, matching the observations.

We find that all of the hot models (left panel) have extremely similar power
spectra, with the largest inter-model differences not exceeding the
observational uncertainties. Additionally, they are all within _approximately_
one $\sigma$ of the observed values for the entire $\log k$ range, except for
the bins at $\log$$(k/(\mathrm{s\;km}^{-1}))$$=-1.2$ and
$\log$$(k/(\mathrm{s\;km}^{-1}))$$\geq-0.8$. The cold models (right panel)
differ more from one another than the hot ones, especially for higher $k$,
with the largest difference amounting to $\sim 3\sigma$ at
$\log$$(k/(\mathrm{s\;km}^{-1}))$$=-0.9$. As noted earlier, all of the cold
models have relatively low $T_{0}$, but the ordering of the trend in small
scale power is interesting, because the models with higher power are those
with slightly higher temperature at mean density. Thus, for cooler
reionization models it doesn’t seem that the temperature difference is the
driver of the differences. It may instead be that a slightly higher
temperature here is more indicative of a later or ongoing reionization,
meaning the small structure of the IGM has had less time to react to
reionization-induced heating, resulting in a smaller pressure smoothing scale.

Also interesting is the shape of the cold models, which compared to the hot
models show a steeper average slope for
$\log$$(k/(\mathrm{s\;km}^{-1}))$$<-1.6$ and a flatter slope for
$\log$$(k/(\mathrm{s\;km}^{-1}))$$>-1.4$. These models clearly do not create
the correct shape in the power spectrum to match Boera et al. (2019). However,
at the same time, the hot models don’t produce sufficient
$\log$$(k/(\mathrm{s\;km}^{-1}))$$>-0.9$ power to match the observations, and
both the hot and cold models may be too steep at
$\log$$(k/(\mathrm{s\;km}^{-1}))$$<-1.8$. The latter could potentially be
alleviated by introducing UVB fluctuations, which would preferentially impact
large scales (Oñorbe et al., 2019). For the mismatch at small scales, perhaps
an intermediate injection temperature, $\approx$15000 K, is necessary.

Figure 9: The flux autocorrelation functions from $z=5.1$ to $z=6.0$
overplotted with the results from Wolfson et al. (2023a). The mean fluxes of
each skewer set have been adjusted to match the observed quantities in each
$\Delta z=0.1$ redshift bin. The simulations generally lie low compared to the
observed points especially at smaller lags, $<700$ km/s. The exceptions are
models that still have significant neutral material, such as LateShortAsym at
$z=6.0$ and EarlyLongSym in all $z$ bins, although the latter is still
systematically low. Figure 10: The same as Figure 9, but for the cold models.
The model differences here are smaller than for the hot models, likely as a
result of reduced temperature fluctuations given the lower maximum heat
injection at reionization. They are nearly identical to the hot models at lags
$\Delta v>500$ km/s, but lower at smaller lags.

####  3.2.2 Ly$\alpha$ autocorrelation function

The flux autocorrelation function $\xi(x)$ goes into the calculation of the
power spectrum, and so contains another representation of the same
information. However, it is easier to measure than the power spectrum, and is
less impacted by observational necessities such as masking. Therefore, it is
useful to understand its sensitivity to the reionization history. Typically,
the autocorrelation is given as the function of the velocity lag between
pixels and calculated on the flux contrast. However, at $z>5$ when the mean
flux is quite low, this can result in large vertical shifts in the function
from a relatively small uncertainty in the mean flux, so we elect to calculate
it on the raw flux field instead. For each raw model skewer, we first apply
smoothing to account for the approximate instrumental resolution of X-Shooter,
$R\sim 8800$, and then calculate the function as

$\xi_{F}\left(\Delta v\right)=\langle F(v)F(v+\Delta v)\rangle$ (14)

where $F(v)$ is the flux in the pixel at velocity $v$ and $\Delta v$ is the
applied velocity lag. Since the simulation reinforces a periodic boundary
condition at the domain edges, there is a close physical connection between
pixels at opposite ends of each skewer, so we account for this in our
calculation. We plot the results for the hot models from $5.1\leq z\leq 6.0$
in Figure 9 and the cold models in Figure 10. Plotted alongside are the
results from the Wolfson et al. (2023a) analysis of XQR-30 data, where we have
matched our skewers to their measured mean fluxes for each $\Delta z$ bin.

With the exception of LateShortAsym, and the two Sym models, the values lie
very close together in all redshift bins, showing the greatest inter-model
deviations mostly at smaller lags. At $\Delta v>700$ km/s or so, nearly all
the models converge to $\langle F\rangle^{2}$, which is expected for
$\xi_{F}$, whereas for the more standard $\xi_{\delta_{F}}$ it would converge
to $\sim 0.0$. LateShortAsym lies above the other models at $z=6.0$, but has
settled by $z=5.9$. This rapid change can be explained by the continuing
reionization of this model at these redshifts: it is still 60 per cent neutral
at $z=6$ and does not finish reionizing until $z=5.9$, so large fluctuations
in the field of neutral gas cause an elevated $\xi_{F}$ at $z=6$. The two Sym
models both have higher $\xi_{F}$ than the others, and the more neutral of the
two (EarlyLongSym) is consistently the highest for $z\leq 5.9$. This persists
down to $z=5.3$, at which point EarlyShortSym finishes reionizing and rises to
match EarlyLongSym at $\Delta v<500$ km/s.

With regards to the data from Wolfson et al. (2023b), while the modeled
autocorrelation values often fall within the 1$\sigma$ limits, they fall
systematically low, particularly at $5.4\leq z\leq 5.8$. The agreement is a
bit better for $z\leq 5.2$ and $z\geq 5.9$, at the higher $z$ especially
because the measurements are quite noisy. They also match the observations
better at larger lags, where everything beings to converge to $\langle
F\rangle^{2}$.

The general trends in the cold models (Figure 10) are the same, but there are
some differences with respect to the hot models. First, the cold $\xi_{F}$
values are slightly lower, except for LateShortAsym at $z\geq 5.9$. Second,
they are lower specifically at $\Delta v<500$ km/s, but converge to the same
values as the hot models at larger velocity lags. It seems that the hot models
are generating greater temperature fluctuations to elevate the overall
variance and by extension the zero-lag autocorrelation function.

####  3.2.3 Cumulative distribution of effective optical depths

While the mean flux and the associated effective optical depth
$\tau_{\mathrm{eff}}$ measured using quasar sightlines is a useful global
statistic, there is more information contained within the distribution of mean
fluxes along each individual line of sight, $\langle F\rangle_{\mathrm{los}}$.
In particular, it is another way of characterizing the amount of scatter in
the opacity, with a broader distribution reflecting greater variation in the
opacities, and the position of the center of the distribution being reflective
of the typical amount of saturation (Fan et al., 2023).

We calculate the mean flux and the corresponding effective optical depth for
each model skewer, and plot their cumulative distribution for the hot and cold
models in Figures 11 and 12. We overplot the measurements from Bosman et al.
(2022), measured along 20 $h^{-1}$Mpc segments of the XQR-30 sightlines. We
adjust the ensemble mean flux to the average flux of all their sightlines,
since this does not exactly match the reported mean flux per $\Delta z=0.1$
bin. For their data, we adhere to their convention for lower and upper limits,
and these are indicated by the gray shading in the figures.444Bosman et al.
(2022) limits are based on non-detections in their dataset, defined as lines
of sight where the local mean flux is less than $2\sigma_{\langle F\rangle}$.
For these non-detections, the bounds are determined by setting $\langle
F\rangle_{\mathrm{los}}=0.0$ and $2\sigma_{\langle F\rangle}$.

First considering the hot models, MidAsym and EarlyAsym are fairly similar to
one another across the entire redshift range, differing mostly in the upper 40
percent of the CDF at the higher $\tau_{\mathrm{eff}}$ end of the
distribution. Within this subset of models, the two Long ones tend to have a
slightly wider distribution. The Late history is extremely wide at
$5.9<z<6.0$, but has settled to match MidShortAsym by $z=5.8$, once it has
finished reionizing and the Ly$\alpha$ opacity is no longer impacted by
remaining opaque neutral regions. It is interesting that these two histories
are identical in $A_{z}$ and $\Delta z$, and result in virtually identical
CDFs despite their different $z_{\mathrm{mid}}$.

The Sym histories are consistently the most different from the other models,
since they only complete reionization at $z<6$ they retain cells that are
wholly opaque. EarlyLongSym is the most impacted by this, since it has not
even reached one percent neutral by $z=5$. However, EarlyShortSym finishes
reionizing at $z=5.4$, and remains a bit wider compared to the other models
even by $z=5.0$. It seems that the longer the duration and the more uniform
the distribution of reionization redshifts (center panel of Figure 2), the
wider the $\tau_{\mathrm{eff}}$ CDF, likely as a byproduct of the temperature
fluctuations introduced by a prolonged reionization process. Of course, there
is also an impact from the ongoing reionization in the Sym models, and this is
why it will be necessary in the future to do a more thorough exploration of
the reionization history parameter space.

The cold models in Figure 12 have narrower distributions than the hot ones,
with the exception of the Late history. Here, it is apparent that the breadth
of the $\tau_{\mathrm{eff}}$ CDF is related to the amplitude of the
autocorrelation function, where a wider CDF is associated with a higher
$\xi_{F}$, which is to be expected given that these are both probes of the
opacity fluctuations.

The models approximately overlap Bosman et al. (2022) best at $z\leq 5.2$ and
$z\geq 5.8$, when the observed CDF is quite narrow and wide, respectively. In
most of the other redshift bins, the models tend to be too narrow, missing
sightlines with more extreme low and high opacities. The exception is the
EarlyLongSym model, which is decently overlapping the upper limits, where the
flux has been “optimistically” set to $2\sigma_{\mathrm{F}}$, in all $z$ bins.

Figure 11: Effective optical depth distributions for all models in $\Delta
z=0.1$ intervals, plotted with observations using the data from Bosman et al.
(2022) reprocessed for segment lengths matching our domain size, obtained via
private communication. The mean flux for each bin is set to the mean of all
the observed lines of sight, because this does not match the reported mean
flux in each redshift bin. The models are generally very similar to one
another, except in the upper 40 per cent of the CDF, and show a narrower
distribution than the data. The exception is the late-finishing EarlyLongSym,
which generally matches the upper limits of the observed distribution. Figure
12: Same as Figure 11, but for the cold models. As with the autocorrelation
functions, there are reduced inter-model differences in the case of cold
models. Additionally, the cold models have narrower CDFs than the hot ones.

The excessive narrowness of the majority of the model CDFs, especially when
combined with the systematically low $\xi_{F}$ values in the previous section,
is indicative that in the future our models must introduce greater
fluctuations in the temperature, neutral fraction, or UVB. These models
intrinsically include fluctuations in temperature because of the inclusion of
inhomogeneous reionization. However, given the trends in the
$\tau_{\mathrm{eff}}$ CDF it seems likely that our use of highly asymmetric
histories in the Asym models is suppressing temperature fluctuations,
narrowing the distributions, causing a greater mismatch with Bosman et al.
(2022). Additionally, while the Sym models approximately span the range of
“permitted” reionization histories ending at $z<6$, there are intermediate
histories that could also be considered. Lastly, these simulations do not
contain a varying UVB, although this will be added in the future, and will
certainly help to alleviate the discrepancy.

##  4 Discussion

In this work, we have used the semi-numerical code AMBER to implement
inhomogeneous reionization histories into the cosmological hydrodynamics code
Nyx. Running at a higher resolution ($\sim 10$ $h^{-1}$kpc) than most
comparable investigations, and using a variety of reionization timings, we
have shown how certain histories can better reproduce various metrics of
Ly$\alpha$ in the $5<z<6$ range, and how certain metrics are systematically
discrepant regardless of the model.

It is common for numerical simulations to fail to reproduce the relatively
high values of the $\log$$(k/(\mathrm{s\;km}^{-1}))$$>-0.9$ end of the
Ly$\alpha$ flux power spectrum of Boera et al. (2019), and our hot models also
reproduce this issue. However, we were able to overshoot the power on this
scale by implementing histories with a lower heat injection than is commonly
assumed while still decently reproducing the observationally-derived $T_{0}$.
Thus it seems that some intermediate temperature could be found to match the
power precisely. However, these results only further confirm those of previous
studies, to emphasize the degeneracies inherent in these metrics (e.g. Wu et
al., 2019), since despite the wildly different pressure smoothing scales in
our hot models, their $z=5.0$ power spectra are essentially indistinguishable.
This motivates the need for experiments involving careful inference over many
reionization models covering a wider range of the parameters considered here.

Related to this, Khan et al. (2023) recently investigated the effect of
initial conditions on Ly$\alpha$ forest power and found that use of a glass
particle initialization scheme can produce a large excess in the power at
$\log$$(k/(\mathrm{s\;km}^{-1}))$$>-1.3$, up to 50 percent at higher $z$. This
appears to result from the impact of particle coupling on both the resultant
density field and peculiar velocities, and in an SPH simulation may be reduced
by using the gravitational softening length of SPH particles (adaptive
softening) rather than using the scale determined by dark matter particles. We
use glass initialization for both dark matter and gas particles, with the gas
particles having a staggered displacement with respect to the dark matter.
Although our setup is significantly different than in Khan et al. (2023), it
is possible that our high $k$ power could arise from this unphysical coupling
effect. If this is the case, then our cold models with their currently
excessive high $k$ power would likely provide a decent match to the Boera et
al. (2019) $z=5$ power spectrum, disfavoring the commonly assumed $\Delta
T_{\mathrm{re}}$ of 20000 K.

The temperature at mean density is not directly observable, but still provides
a useful point of comparison for the models, especially since the amount of
heat injection is not certain. The majority of observations find quite low
$T_{0}$ from $5.0<z<5.8$, less than 9000 K, which several of our models match
at $z\sim 5.0$ but overshoot at higher redshifts. These have varying heat
injection and reionization histories. The hotter derived $T_{0}$ values of
Gaikwad et al. (2020) are matched by all hot models except for the very late
and rapid LateShortAsym history, which is far _too_ hot. It doesn’t seem
possible to achieve the very low temperatures of Garzilli et al. (2017) or
Walther et al. (2019), unless using a very low heat injection, leading to a
flat evolution at $5<z<6$. The Gaikwad et al. (2020) measurements are often
interpreted as suggesting a decrease in temperature and thus a period of
significant cooling following the end of reionization, which is only a feature
of the hot histories, and is also stronger for Asym histories. Generally
though, it seems impossible to choose a preferred reionization history on the
basis of the observationally-derived temperature at mean density, until they
begin to converge to some distinctive trend.

Considering the more direct measures of gas phase, Ly$\alpha$ opacity metrics,
we find that the metrics appear to be sensitive to different details of the
reionization history. Histories where the reionization midpoint occurs later
have less relaxed $\rho-T$ distributions and smaller pressure smoothing
scales, which leaves the potential for them to have increased small scale
fluctuations. As a result, we can see that the small scale
$P_{\mathrm{Ly\alpha}}$ in $\Delta T_{\mathrm{re}}=10,000$ K reionization
histories can give an indication of the reionization midpoint, and the
pressure smoothing scale. However, in practice these fluctuations are only
observable if $\Delta T_{\mathrm{re}}$ is sufficiently low, less than $20,000$
K, otherwise they are completely washed out in the 1D flux power spectrum.
$\xi_{F}$ and the $\tau_{\mathrm{eff}}$ CDF on the other hand appear to be
sensitive to the duration $\Delta z$ and asymmetry $A_{z}$ of the histories,
with a longer or more symmetric history leading to higher $\xi_{F}$ and wider
CDFs. This may be because the more varied timing in heat injections leads to
greater temperature fluctuations.

Our simulations can span the 1D flux power spectrum measurements easily, but
we are less able to match the autocorrelation functions and effective optical
depth distributions. Given the trends we see with regard to $\Delta z$ and
$A_{z}$, exploration of more Sym-like histories could partially alleviate
this. However, the main issue is likely the lack of UVB fluctuations in this
model, as Wolfson et al. (2023b) demonstrated that more fluctuations through a
smaller mean free path would increase the autocorrelation function values
across the range of lags considered here. Further, this change would primarily
impact large scales (Oñorbe et al., 2019), thus minimally affecting the $k$
modes covered in our 1D flux power spectra. This supports results from both
Bosman et al. (2022) and Zhu et al. (2022) that rule out a homogeneous UVB for
$z\gtrapprox 5.4$.

Most of our reionization histories were chosen on the basis that they adhered
to the _Planck_ FlexKnot (originating in Millea & Bouchet, 2018) limits in
volume-weighted neutral fraction, which only permits “knots” at $z>6$ in order
to match the reionization redshift range implied by observations of Gunn-
Peterson troughs. It thus enforces $x_{\mathrm{HI,v}}<0.02$ by $z=5.8$, which
is arguably not the best assumption given the current Ly$\alpha$ forest
evidence for a late reionization. We consider only one model that violates
this result, and instead adheres to the dark pixel fraction measurements from
McGreer et al. (2015) and the directly measured $\tau_{\mathrm{e}}$. In the
future, a better choice would be to explore a larger number of EarlyLongSym-
like models, that permit $x_{\mathrm{HI,v}}>0$ for $z<6$.

##  5 Conclusions

In this work we have evaluated how a variety of reionization histories affect
the IGM in Nyx cosmological simulations and several resultant Ly$\alpha$
forest metrics: the 1D flux power spectrum, autocorrelation function, and
effective optical depth distribution. We accomplished this using the semi-
numerical code AMBER, which allowed us to parametrize the evolution of the
neutral fraction with redshift with the shape of the reionization histories
rather than the details of the physics. We have found that:

  * •

Reionization histories with longer durations and greater symmetry produce gas
that spans a wide range of properties in $\rho-T$ space. This naturally will
lead to more fluctuations in the temperature and density fields.

  * •

The ensemble pressure smoothing scale measured in the IGM grows linearly with
time since the mass-weighted midpoint, and hotter heat injection (20,000
versus 10,000 K) leads to a faster growth. This scale is also correlated with
the adiabatic index $\gamma$ as measured from the $\rho-T$ distribution. Low
heat injection models also have a higher adiabatic index for a given pressure
smoothing scale.

  * •

Models with lower heat injection at reionization (10,000 K) have shorter
pressure smoothing scales in the IGM, and their average temperature post-
reionization is low enough that the remaining small scale power is still
visible in the 1D Ly$\alpha$ flux power spectrum.

  * •

Hot ($\Delta T_{\mathrm{re}}=20000$ K) reionization histories better reproduce
the shape of the observed 1D Ly$\alpha$ flux power spectrum than cold ($\Delta
T_{\mathrm{re}}=10000$ K) histories. However, hot histories also underproduce
observed $\log$$(k/(\mathrm{s\;km}^{-1}))$$>-0.9$ power, while the cold ones
overproduce it.

  * •

All of the models systematically underestimate the observed autocorrelation
function, and the disagreement is worse at $5.3\leq z\leq 5.8$ than at $z\leq
5.2$ and $z\geq 5.9$. The models that tend to match the best are those that
don’t finish reionizing until $z\leq 5.4$ and that have more symmetric
reionization histories, or an approximately normal distribution of
reionization redshifts. Hotter $\Delta T_{\mathrm{re}}$ models produce higher
autocorrelation values than cold ones, but only at $\Delta v<500$ km/s.

  * •

Similar to the autocorrelation function, the models underpredict the observed
width of the distribution of effective optical depth measurements, with the
exception of our extreme history, that does not even reach
$x_{\mathrm{HI,v}}=1$ per cent by $z=5.0$. Simulations with similar durations
and asymmetries produce similar $\tau_{\mathrm{eff}}$ CDFs, suggesting that
the differences here arise from temperature or density fluctuations in the
post-reionization IGM.

Notably, our experimental setup includes fluctuations in temperature and
hydrogen neutral fraction, but excludes any kind of ionizing background
fluctuation. Greater spatial variation in the UVB is known to contribute to a
higher $\xi_{F}$ and wider CDFs, and has already been deemed necessary for
matching observations down to at least $z=5.3$ (Bosman et al., 2022; Zhu et
al., 2022). Thus, in order to properly compare to observations, it seems
absolutely necessary to include UVB fluctuations.

This work demonstrates the sensitivity of the Ly$\alpha$ forest metrics to the
reionization history, but also shows the situations where the effects of model
parameters become degenerate. In order to place quantitative constraints on
the history, it will be necessary to use a larger model grid that captures all
of the competing physical effects, and with fewer arbitrary restrictions
placed on the histories, such as forcing an end at $z\geq 6$.

## Acknowledgements

CCD thanks Hyunbae Park and Jean Sexton for assistance in implementing
inhomogeneous reionizations in Nyx. CCD thanks Hy Trac for assistance in
working with AMBER, and Sarah Bosman for providing re-calculated effective
optical depth distributions. This research used resources of the National
Energy Research Scientific Computing Center, which is supported by the Office
of Science of the U.S. Department of Energy under Contract No. DE-
AC02-05CH11231. This research further used resources of the Oak Ridge
Leadership Computing Facility at the Oak Ridge National Laboratory, which is
supported by the Office of Science of the U.S. Department of Energy under
Contract No. DE-AC05-00OR22725.

## Data Availability

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

## References

  * Almgren et al. (2013) Almgren A. S., Bell J. B., Lijewski M. J., Lukić Z., Van Andel E., 2013, ApJ, 765, 39
  * Bañados et al. (2018) Bañados E., et al., 2018, Nature, 553, 473
  * Battaglia et al. (2013) Battaglia N., Trac H., Cen R., Loeb A., 2013, ApJ, 776, 81
  * Becker et al. (2021) Becker G. D., D’Aloisio A., Christenson H. M., Zhu Y., Worseck G., Bolton J. S., 2021, MNRAS, 508, 1853
  * Boera et al. (2019) Boera E., Becker G. D., Bolton J. S., Nasir F., 2019, ApJ, 872, 101
  * Bosman et al. (2022) Bosman S. E. I., et al., 2022, MNRAS, 514, 55
  * Bouchet et al. (1995) Bouchet F. R., Colombi S., Hivon E., Juszkiewicz R., 1995, A&A, 296, 575
  * Bruton et al. (2023) Bruton S., Lin Y.-H., Scarlata C., Hayes M. J., 2023, arXiv e-prints, p. arXiv:2303.03419
  * Chabanier et al. (2023) Chabanier S., et al., 2023, MNRAS, 518, 3754
  * Costa et al. (2014) Costa T., Sijacki D., Trenti M., Haehnelt M. G., 2014, MNRAS, 439, 2146
  * Davies et al. (2018) Davies F. B., et al., 2018, ApJ, 864, 142
  * DeBoer et al. (2017) DeBoer D. R., et al., 2017, PASP, 129, 045001
  * Doughty et al. (2023) Doughty C. C., Hennawi J. F., Davies F. B., Lukić Z., Oñorbe J., 2023, arXiv e-prints, p. arXiv:2305.16200
  * Fan et al. (2023) Fan X., Bañados E., Simcoe R. A., 2023, ARA&A, 61, 373
  * Feng et al. (2016) Feng Y., Chu M.-Y., Seljak U., McDonald P., 2016, MNRAS, 463, 2273
  * Gaikwad et al. (2020) Gaikwad P., et al., 2020, MNRAS, 494, 5091
  * Gaikwad et al. (2023) Gaikwad P., et al., 2023, arXiv e-prints, p. arXiv:2304.02038
  * Garzilli et al. (2017) Garzilli A., Boyarsky A., Ruchayskiy O., 2017, Physics Letters B, 773, 258
  * Gnedin & Hui (1998) Gnedin N. Y., Hui L., 1998, MNRAS, 296, 44
  * 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
  * HERA Collaboration et al. (2023) HERA Collaboration et al., 2023, ApJ, 945, 124
  * Heintz et al. (2023) Heintz K. E., et al., 2023, arXiv e-prints, p. arXiv:2306.00647
  * Hoag et al. (2019) Hoag A., et al., 2019, ApJ, 878, 12
  * Hockney & Eastwood (1988) Hockney R. W., Eastwood J. W., 1988, Computer simulation using particles
  * Hui & Gnedin (1997) Hui L., Gnedin N. Y., 1997, MNRAS, 292, 27
  * Jacobus et al. (2023) Jacobus C., Harrington P., Lukić Z., 2023, arXiv e-prints, p. arXiv:2308.02637
  * Jin et al. (2023) Jin X., et al., 2023, ApJ, 942, 59
  * Kannan et al. (2022) Kannan R., Garaldi E., Smith A., Pakmor R., Springel V., Vogelsberger M., Hernquist L., 2022, MNRAS, 511, 4005
  * Karaçaylı et al. (2022) Karaçaylı N. G., et al., 2022, MNRAS, 509, 2842
  * 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
  * Khan et al. (2023) Khan N. K., Kulkarni G., Bolton J. S., Haehnelt M. G., Iršič V., Puchwein E., Asthana S., 2023, arXiv e-prints, p. arXiv:2310.07767
  * Kulkarni et al. (2015) Kulkarni G., Hennawi J. F., Oñorbe J., Rorai A., Springel V., 2015, ApJ, 812, 30
  * Lewis et al. (2000) Lewis A., Challinor A., Lasenby A., 2000, ApJ, 538, 473
  * Loeb & Barkana (2001) Loeb A., Barkana R., 2001, ARA&A, 39, 19
  * Lukić et al. (2015) Lukić Z., Stark C. W., Nugent P., White M., Meiksin A. A., Almgren A., 2015, MNRAS, 446, 3697
  * Mason et al. (2018) Mason C. A., Treu T., Dijkstra M., Mesinger A., Trenti M., Pentericci L., de Barros S., Vanzella E., 2018, ApJ, 856, 2
  * McGreer et al. (2015) McGreer I. D., Mesinger A., D’Odorico V., 2015, MNRAS, 447, 499
  * Mesinger et al. (2011) Mesinger A., Furlanetto S., Cen R., 2011, MNRAS, 411, 955
  * Millea & Bouchet (2018) Millea M., Bouchet F., 2018, A&A, 617, A96
  * Oñorbe et al. (2019) Oñorbe J., Davies F. B., Lukić Z., Hennawi J. F., Sorini D., 2019, MNRAS, 486, 4075
  * O’Leary & McQuinn (2012) O’Leary R. M., McQuinn M., 2012, ApJ, 760, 4
  * Ouchi et al. (2010) Ouchi M., et al., 2010, ApJ, 723, 869
  * Planck Collaboration et al. (2020) Planck Collaboration et al., 2020, A&A, 641, A6
  * Puchwein et al. (2023) Puchwein E., et al., 2023, MNRAS, 519, 6162
  * Rosdahl et al. (2018) Rosdahl J., et al., 2018, MNRAS, 479, 994
  * Rosdahl et al. (2022) Rosdahl J., et al., 2022, MNRAS, 515, 2386
  * Scoccimarro (1998) Scoccimarro R., 1998, MNRAS, 299, 1097
  * Sexton et al. (2021) Sexton J., Lukic Z., Almgren A., Daley C., Friesen B., Myers A., Zhang W., 2021, The Journal of Open Source Software, 6, 3068
  * Trac et al. (2022) Trac H., Chen N., Holst I., Alvarez M. A., Cen R., 2022, ApJ, 927, 186
  * Viel et al. (2013a) Viel M., Becker G. D., Bolton J. S., Haehnelt M. G., 2013a, Phys. Rev. D, 88, 043502
  * Viel et al. (2013b) Viel M., Schaye J., Booth C. M., 2013b, MNRAS, 429, 1734
  * Walther et al. (2018) Walther M., Hennawi J. F., Hiss H., Oñorbe J., Lee K.-G., Rorai A., O’Meara J., 2018, ApJ, 852, 22
  * Walther et al. (2019) Walther M., Oñorbe J., Hennawi J. F., Lukić Z., 2019, ApJ, 872, 13
  * Wang et al. (2020) Wang F., et al., 2020, ApJ, 896, 23
  * Weibull (1951) Weibull W., 1951, Journal of Applied Mechanics, 18, 293
  * Wolfson et al. (2023a) Wolfson M., et al., 2023a, arXiv e-prints, p. arXiv:2309.03341
  * Wolfson et al. (2023b) Wolfson M., Hennawi J. F., Davies F. B., Oñorbe J., 2023b, MNRAS, 521, 4056
  * Wu et al. (2019) Wu X., McQuinn M., Kannan R., D’Aloisio A., Bird S., Marinacci F., Davé R., Hernquist L., 2019, MNRAS, 490, 3177
  * Yang et al. (2020) Yang J., et al., 2020, ApJ, 897, L14
  * Zel’dovich (1970) Zel’dovich Y. B., 1970, A&A, 5, 84
  * Zhu et al. (2022) Zhu Y., et al., 2022, ApJ, 932, 76
  * Zhu et al. (2023) Zhu Y., et al., 2023, arXiv e-prints, p. arXiv:2308.04614