11institutetext:  Observatoire de Paris, PSL Research University, Sorbonne
Université, CNRS, LERMA, 75014 Paris, France

# The LoReLi database: 21 cm signal inference with 3D radiative hydrodynamics
simulations

R.Meriot 11 B.Semelin 11

The Square Kilometer array is expected to measure the 21cm signal from the
Epoch of Reionization in the coming decade, and its pathfinders may provide a
statistical detection even earlier. Current reported upper limits start
putting constraints on the astrophysical parameters of the models of the Epoch
of Reionization (EoR). In order to interpret such data with 3D radiative
hydrodynamics simulations using Bayesian inference, we present the latest
developments of the Licorice code. Relying mainly on an implementation of the
halo Conditional Mass Function to account for unresolved star formation, they
allow for accurate simulations of the EoR at $256^{3}$ resolution. We use this
version of Licorice to produce the first iteration of LoReLi, a public dataset
now containing hundreds of 21cm signals computed from radiative hydrodynamics
simulations. We train a neural network on LoReLi to provide a fast emulator of
the Licorice powerspectra, LorEMU, that has $\sim 5\%$ rms error relative to
the simulated signals. LorEMU is used in an MCMC framework to perform Bayesian
inference, first on a mock observation composed of a simulated signal and a
thermal noise corresponding to 100h observations with the SKA. We then apply
our inference pipeline to the latest measurements from the HERA
interferometer. We report constraints on the X-ray emissivity, and confirm
that cold reionization scenarios are unlikely to accurately represent our
Universe.

###### Key Words.:

Cosmology: dark ages, reionization, first stars – Radiative transfer– Early
Universe – Methods: numerical

##  1 Introduction

The first billion years are key in the evolution of the Universe on large
scales. Following the hierarchical theory of structure formation and the
physics regulating star formation, the gas contained in massive halos cools
down, subsequently fragmenting, collapsing and giving birth to the first
stars. This marks the end of the dark ages. There are several processes
responsible for this cooling. The atomic hydrogen cooling channel allows
significant cooling down to $\sim 10^{4}$ K, that is below the virial
temperature of halos with mass larger than $M_{min}\approx
10^{8}\,\rm{M_{\odot}}$, allowing gravothermal collapse and star formation.
Molecular hydrogen cooling brings this minimal mass for star formation down to
$M_{min}\approx 10^{5}\,\rm{M_{\odot}}$, in so-called ”mini halos”. The
emergence of the first stars and galaxies, in a period called the Cosmic Dawn
(CD, $z\sim 15-30$), greatly impacted the neutral gas in the intergalactic
medium (IGM). UV radiation of the first generation of sources ionized their
local environment, while the still neutral IGM, farther away from the sources,
became heated by X-ray emissions. During the Epoch of Reionization (EoR,
$z\sim 6-15$), the ionized regions quickly grew until they covered the entire
IGM. The 21 cm line of neutral hydrogen is emitted by the neutral hydrogen in
the IGM and its intensity fluctuations are shaped by the astrophysical
properties of the first sources, making it an extremely promising window to
the first billion years. Furlanetto et al. (2006) provide a review of the
topic.

While valuable, the detection of this signal remains an observational
challenge. Foreground emissions, roughly 10 000 times stronger than the
expected intensity of the 21 cm signal, must be subtracted from observed data
or avoided in order to extract the cosmological signal. The subtraction
requires a very good handling of radio frequency interference, a deep and
exact sky model and accurate direction dependent calibration (see e.g. Mertens
et al., 2020). Many instrumental programs are tackling this challenge, the
most ambitious being the Square Kilometer Array (SKA). The unprecedented
sensitivity of the SKA will enable us to probe this epoch with sufficient
signal-to-noise ratio to build a full tomography of the signal between $z\sim
27$ and $z\sim 6$, providing a comprehensive view of how reionization
unfolded. The observation of the 21 cm signal by the SKA will start by the end
of this decade, but several instruments are already trying to measure summary
statistics of the signal. A first detection of the global (sky-averaged)
signal at $z\sim 16$ has been announced by Bowman et al. (2018) using the
EDGES instrument. The features of this detection, in particular the strong
intensity of the signal in absorption, are not compatible with standard models
of reionization and this results has not been confirmed by other global signal
experiment, such as SARAS-3 (Singh et al., 2022). Interferometers such as
LOFAR, MWA, NenuFAR, and HERA aim at measuring the power spectrum of the 21 cm
signal. While no detection is claimed yet, upper limits of the spectrum have
been devised ( e.g. Mertens et al. (2020); Trott et al. (2020); Abdurashidova
et al. (2023); Munshi et al. (in prep.), and references therein). In the
coming years, these instruments are expected to detect the signal or produce
upper limits informative enough to constrain non-exotic models of the EoR at
various scales and redshifts. Given these observational prospects, in order to
extract astrophysical information, the community must provide accurate models
of the first billion years as well as inference methods linking raw
observational data to the astrophysical processes encoded in the models.

On the one hand, sustained effort has been put into developing theoretical and
numerical approaches to model the EoR, ranging from fast semi-analytical codes
(e.g. Mesinger et al., 2011; Santos et al., 2010; Cohen et al., 2017; Murray
et al., 2020; Reis et al., 2021) to expensive 3D radiative transfer
simulations (e.g. Ciardi et al., 2003; Mellema et al., 2006; Baek et al.,
2010; Semelin et al., 2017; Ocvirk et al., 2020; Lewis et al., 2022; Doussot &
Semelin, 2022). Fast codes compute the dynamics of dark matter in the linear
regime and typically assume a constant gas to dark matter density ratio. On
large scales (more than a few cMpc), this yields similar results to full
numerical simulations that compute the non-linear dynamics of gravity and
hydrodynamics. However, star formation, a crucial process when modelling the
EoR, occurs on small scales. This leads to very different modelling
approaches. In high-resolution N-body simulations, individual dark matter
halos are resolved. Star formation within the halos can then be computed from
the local properties of the gas (typically density and temperature), though
the actual equations linking gas properties to star formation may differ
between simulation codes. 1D (spherically symmetric) radiative transfer
simulations assign baryons to halos formed in pre-computed dark-matter-only
simulation (Krause et al., 2018; Ghara et al., 2018; Schaeffer et al., 2023)
using simple recipes, and compute star formation from there. Semi analytical
codes usually do not resolve individual halos but estimate populations of
halos and stars in a given region, through for instance the CMF formalism
(Lacey & Cole, 1993). These different modelling approaches can potentially
induce different features in the reionization field. The second critical step
is the process through which UV and X-rays are accounted for. Fully numerical
simulations typically include 3D radiative transfer through the use of M1 or
various ray-tracing methods, while semi-analytical codes identify ionized and
heated region through simple (and easy to compute) photon budget arguments. A
comprehensive evaluation of the impact of these different modelling approaches
on the interpretation of the 21-cm signal remains to be performed.

On the other hand, applications of inference methods to the 21 cm signal have
been developing in the last 10 years. The most widely employed technique is
undoubtedly the Markov Chain Monte Carlo (MCMC) algorithm (see e.g. Greig &
Mesinger (2015, 2018); Schmit & Pritchard (2018)). In this approach, the
posterior distribution of the parameters is explored in a guided random walk.
This method is guaranteed to converge towards the true posterior probability
(under some reasonable assumptions). However, it requires that an explicit
likelihood can be formulated and it might require hundreds of thousands or
millions of steps (and thus of realisations of the model) to produce a good
estimate of the posterior, making its use only suited to fast numerical codes.
This led the community to explore other inference techniques, such as training
artificial neural networks to do the inference (Shimabukuro & Semelin, 2017;
Doussot et al., 2019; Kendall & Gal, 2017; Neutsch et al., 2022). While
requiring fewer simulations in the training sample, convergence towards the
true maximumlikelihood or true posterior is in this case not guaranteed, as
the training of the neural network can produce a systematic bias or a residual
additional variance. Another noteworthy family of methods is called
”simulation-based inference” or ”likelihood-free inference” (Alsing et al.,
2017; Alsing & Wandelt, 2018; Zhao et al., 2022; Prelogović & Mesinger, 2023).
It relies on the fact that a simulation is in and of itself a draw in the
joint probability of the observable data and the model parameter, a joint
probability that need not necessarily be explicitly writable. The joint
probability distribution is then fitted using a neural network acting as a
parametric function, allowing for a straightforward computation of the
posterior.

Due to computation time constraints, most inference methods have until now
only been applied to signals generated using semi-numerical approaches, as
they typically are a hundred thousand times faster than fully numerical
simulations. Indeed, performing a single MCMC inference with e.g. the Licorice
code running at high resolution would require $10^{11}$ CPU hours of
computation, well beyond current capabilities. Performing inference with such
a code at a reasonable cost therefore requires the use of machine learning
based approaches, that can reduce the number of simulations to a few
thousands. But even in this case, a decrease in resolution is necessary to
reach an affordable computing cost. This decrease must be made up for by a
significant amount of subgrid modelling. This constitutes the approach at the
heart of this work, where we use of the order of $10^{3}$ simulations to chart
a 4D astrophysical parameters space with $256^{3}$ simulations in $300$ cMpc
boxes. When then train neural network emulators for the power spectrum of the
signal and use them to perform MCMC inference.

In this paper, we present LoReLi, a collection of numerical simulations of the
EoR using the Licorice code. In section 2, we explain the recent modifications
of the code that make such a database possible. In section 3, we give details
about the LoReLi database, and in section 4, we show the results of parameter
inference applied to both mock data and the most recent data from the HERA
collaboration (Abdurashidova et al., 2022), thanks to an emulator of the code
trained on LoReLi.

##  2 The Licorice code

The simulations of the LoReLi dataset were run using the Licorice simulation
code (Semelin et al., 2007; Baek et al., 2009; Semelin, 2015; Semelin et al.,
2017). We give here a brief overview of the code and details about recent
developments. Licorice is an N-body simulation code dealing with baryons and
dark matter particles. Dynamics is solved using a TREE+SPH method (Springel et
al., 2001). Radiative transfer in UV and X-ray continuum is coupled to
dynamics through a Monte-Carlo scheme on an adaptive grid : photon packets are
emitted in random directions by source particles (baryon particles with a non-
zero stellar fraction), propagate at the speed of light, redshift, and deposit
their energy in gas particles lying in the encountered cells depending on
local density. They represent UV ionizing radiation and X-rays. ”Hard” X-rays,
after propagating over a distance larger than the size of the simulation box,
are added to a uniform background. The ionization and temperature states of
each gas particle are then updated several times per dynamical time step and
the latter affect dynamics through the pressure and artificial viscosity terms
in the SPH scheme (Monaghan, 1992; Springel et al., 2001).

###  2.1 Star formation model

In previous iterations of Licorice, gas particles with an SPH-calculated
overdensity $\delta$ greater than a fixed threshold $\delta_{thresh}$
transform at each dynamical time step a fraction of their gas mass into
stellar mass following a Schmidt-Kennicutt law of exponent one, using

$df_{*}=(f_{coll}-f_{*})\dfrac{dt}{\tau}$ (1)

where $f_{*}$ is the current stellar fraction of the particle, $f_{coll}=1$ if
$\delta>\delta_{thresh}$ and $f_{*}$ otherwise, $df_{*}$ the newly created
stellar fraction, $dt$ the duration of a time step, and $\tau$ an
astrophysical parameter of the simulation representing the typical gas-to-
stars conversion timescale, thus controlling the star formation rate (SFR).
However, this star formation numerical method is sensitive to an insufficient
resolution : not resolving the small scale modes of the density field will
smooth out the density peaks and prevent particles from reaching the density
threshold for star formation, causing Eq. (1) to underestimate the SFR.
Essentially, a simulation will lack the star formation of unresolved star-
forming halos. The mass of the lightest of these halos is estimated to be
around $\sim 10^{8}\,\rm{M_{\odot}}$, as it corresponds to a virial
temperature of $10^{4}$ K, the lowest temperature reached through the cooling
of atomic hydrogen. In boxes of hundreds of cMpc, large enough to capture
relevant large scale fluctuations of the 21 cm signal, resolving such halos
requires more than $4096^{3}$ particles and $10^{7\sim 8}$ cpuh. This
motivates a new implementation in the code of a subgrid model that aims at
giving a better estimate of star formation in unresolved halos, without
affecting resolved ones.

####  2.1.1 The Conditional Mass Function subgrid model

The chosen subgrid model relies on the Conditional Mass Function (CMF)
formalism (Lacey & Cole, 1993) to estimate the halo mass functions of regions
below the spatial resolution limit. In the CMF theory, the density of spheres
of decreasing radii follows a random walk, and statistically estimating the
radius at which the threshold $\delta_{c}$ is crossed for the first time gives
the CMF of regions parameterized by their radius $R_{0}$ and overdensity
$\delta_{0}$.

In Extended Press-Schechter (EPS), the simplest formulation of the CMF theory,
the fraction of the mass of a region (again parameterized by a radius $R_{0}$
and overdensity $\delta_{0}$) that lies in collapsed halos more massive than
$M_{min}$ is given by

$f_{coll}=\rm{erfc}\left[\frac{\delta_{c}-\delta_{0}}{\sqrt{2(\sigma(M_{min})^{2}-\sigma_{0}^{2}})}\right]$
(2)

Where $\sigma(M)$ is the density variance on mass scale $M$ and $\delta_{c}$
is the linear overdensity of collapse. The classical theory of structure
formation predicts $\delta_{c}\approx 1.68$ .

Sheth et al. (2001) devised a more accurate formalism (ST) to calculate the
mass function. During the random walk, instead of crossing a constant
threshold of size $\delta_{c}$, a ”moving barrier”

$B(\sigma^{2},z)=\sqrt{a}\delta_{c}[1+\beta(a\frac{\delta_{c}^{2}}{\sigma^{2}})^{-\alpha}]$

is adopted. This extends the EPS formalism from the physics of spherical
collapse to that of ellipsoidal collapse. $\alpha\approx 0.615,\beta\approx
0.485$ are parameters whose value come from ellipsoidal dynamics and $a\approx
0.707$ is fitted to accurately model simulations. Rubiño-Martín et al. (2008);
Tramonte et al. (2017) provided an expression for the CMF associated with ST :

$n_{c}(M)=-\sqrt{\frac{2}{\pi}}\frac{d\sigma}{dM}\frac{\rho_{0}}{M}\frac{|T(\sigma^{2}|\sigma_{0}^{2})|\sigma}{(\sigma^{2}-\sigma_{0}^{2})^{3/2}}\rm{exp}\left[-\frac{(B(\sigma^{2},z)-\delta_{0})^{2}}{2(\sigma^{2}-\sigma_{0}^{2})}\right]$
(3)

where

$T(\sigma^{2}|\sigma_{0}^{2})=\sum^{5}_{n=0}\frac{(\sigma_{0}^{2}-\sigma^{2})}{n!}\frac{\partial^{n}[B(\sigma^{2},z)-\delta_{0}]}{\partial(\sigma^{2})^{n}}.$
(4)

One can then easily calculate $M_{coll}$, the mass of the region contained in
halos more massive than $M_{min}$ :

$M_{coll}=V_{region}\int_{M_{min}}^{M_{region}}n_{c}(M)M\rm{dM}$ (5)

where $V_{region}$ is the lagrangian volume of the region and $M_{region}$ its
mass. We then calculate the collapsed fraction as

$f_{coll}=\frac{M_{coll}}{M_{region}}.$ (6)

Both EPS and ST CMFs are implemented in the version of Licorice used in the
present work. In practice, we calculate $f_{coll}$ for each gas particle. To
do so, we calculate the overdensity
$\delta_{0}=\frac{\rho_{0}-\rho_{m}}{\rho_{m}}$ of the region they represent,
where $\rho_{0}$ is the SPH-calculated density of the particle multiplied by
the cosmological ratio $\Omega_{m}/\Omega_{b}$. The volume of the region is
given by $V=\gamma\frac{4}{3}\pi R_{0}^{3}$, where $R_{0}$ is the SPH
smoothing radius. While the formulation of the CMFs relies on tophat filters,
the SPH smoothing kernel is not a tophat, so the $\gamma$ parameter is tuned
to represent the volume of the region of identical density smoothed with a
tophat kernel. Motivated by the fact that the integral of the SPH kernel
reaches 99% of its maximum value past $\sim 0.7R_{0}$, $\gamma=0.4$ is a
reasonable choice, that is validated by the resolution study presented below.

We use the resulting value of $f_{coll}$ in Eq. 1 instead of setting
$f_{coll}$ depending on some density threshold. It is done either with Eq. (2)
when using the EPS formalism, or with Eq. (3), (5) and (6) when using ST.

A final subtlety lies in the fact that in Licorice, all gas particles with a
non-zero stellar fraction are considered sources of photons. However, this CMF
formalism assigns a strictly positive $f_{coll}$ to all particles. This causes
particles to behave in an unphysical manner, as they instantly turn into
sources at the beginning of the simulation with infinitesimal stellar
fractions. To avoid that, a method to stochastically assign $f_{coll}$ was
adopted : when $df_{coll}$, the difference in $f_{coll}$ of a particle between
one time step and the next, is smaller than
$df_{coll,min}=\frac{M_{min}}{M_{region}}$, $df_{coll,min}$ is assigned to the
particle with a probability $\frac{df_{coll}}{df_{coll,min}}$. This
effectively prevents very weak, unphysical star formation at very high
redshifts ($z\gtrsim 25-30$ depending on $M_{min}$, approximately when star
formation starts in high-resolution simulations) while ensuring that the
average star formation density (SFRD) remains identical to that of the non-
stochastic model.

Figure 1: Evolution with redshift of the halo mass function of HIRRAH-21,
calculated using a FoF algorithm. We also show the history of theoretical PS
and ST HMFs for the same cosmology.

Figure 1 shows how the HMF of HIRRAH-21, a high-resolution Licorice simulation
($2048^{3}$ particles, with no subgrid model) presented in Doussot & Semelin
(2022) compares with the PS and ST models at different redshifts. The HMF of
HIRRAH-21 was calculated by finding halos in the simulation using a Friends-
of-Friends (FoF) algorithm. We find an excellent agreement with the ST HMF,
which lies a couple orders of magnitudes above the PS HMF. This is not
unexpected, given the studies in Reed et al. (2007), and this difference in
the HMFs necessarily appears in the CMFs. This justifies the implementation
and use of the ST CMF. However, for practical reasons, in the presented and
current version of the LoReLi database we compute star formation using EPS and
a different choice of an effective radius (10% smaller, which implies a $\sim
5\%$ larger $\sigma$) that produces a very similar SFRD. Future versions of
the database will use the ST CMF. A comparison of the conditional mass
functions of HIRRAH-21 at $z=9.48$ with the prediction of ST and EPS (with a
$\sigma$ scaled up by 5%), at different overdensities $\delta$, can be seen on
figure 2. A good agreement between simulation and models is observed, with
errors roughly within the fluctuations in the simulated CMF due to cosmic
variance. The resulting $f_{coll}$ are within $\sim 30\%$ of each other.

Figure 2: Evolution of the simulated CMF with the overdensity $\delta$ and of
the theoretical CMFs as predicted by the ST and EPS models. In order to match
the high resolution simulation, the $\sigma$ in the EPS model has been scaled
up by $5\%$ on average.

The connection between the simulated HMFs and star formation in the simulation
is indirect however, as star formation in Licorice does not necessarily occur
in identified dark matter halos : only the local density of the gas appears in
the star formation equation and not the properties of the dark matter that may
or may not be present nearby. In order to evaluate the results of this subgrid
model on star formation, we compare in figure 3 the star formation rate
density of HIRRAH-21 and LORRAH-21, a low-resolution simulation ($256^{3}$
particles with the subgrid modelling of star formation). The $M_{min}$
parameter of the subgrid model was set at $4\times 10^{9}\,\rm{M_{\odot}}$ to
match the mass of the lightest halos identified by FoF in HIRRAH-21, and all
other astrophysical and numerical parameters are identical in the two setups
except for the dynamical time step : $0.5\,\rm{Myr}$ in HIRRAH-21 and
$7\,\rm{Myr}$ in LORRAH-21. We observe a very good agreement between the SFRDs
of LORRAH-21 and HIRRAH-21 over the whole redshift range (with an RMS of the
relative error of $0.16$ between $z\sim 20$ and $z\sim 7$, and $0.09$ between
$z\sim 15$ and $z\sim 7$). The largest mismatch ($\sim 30\%$) occurs at $z\sim
18-20$, when the stochastic attribution of $f_{coll}$ is the most relevant. In
order to test this subgrid model at a different resolution, figure 3 also
shows the SFRD in a $1024^{3}$ simulation (from the 21SSD database Semelin et
al. (2017)) with the same parameters as HIRRAH-21 (but a minimal mass of DM
halos naturally set to $3.2\times 10^{10}\,\rm{M_{\odot}}$) and in LORRAH-21,
this time with $M_{min}$ set to the same value. In this case as well, an
excellent agreement is observed, implying that this CMF formalism is robust to
resolution and can accurately describe the behaviour of simulations of higher
resolution.

Figure 3: History of the star formation rates of HIRRAH-21 and LORRAH-21. As a
reference, observational data from J. Bouwens et al. (2016); Oesch et al.
(2018); McLeod et al. (2016) is included. The SFRD of a $1024^{3}$ particles
simulation is also plotted, and is consistent with that of a LORRAH-21
simulation using a $M_{min}$ parameter set at the value of the smallest
resolved halos in the $1024^{3}$ simulation.

####  2.1.2 Two-phase model of gas particles

Figure 4: History of the average temperature of all particles in both
LORRAH-21 and HIRRAH-21, as well as the temperature of weakly ($<2\%$) ionized
particles. In the case of LORRAH-21, the temperature of the neutral phase of
the particles is also shown.  Figure 5:  Global 21 cm signal in HIRRAH-21 and
LORRAH-21, the latter calculated using either the phase-averaged temperature
$\langle T\rangle$ of particles or the temperature of their neutral phase
$T_{HI}$. This shows that properly taking into account the temperature of the
neutral phase is critical to correctly model the 21 cm signal in low-
resolution numerical simulations.

The other significant adaptation of the code is the computation of a specific
temperature for the neutral phase within each gas particle. In previous
versions of the code, only the average temperature of particles was considered
(see Baek et al. (2009) for details) and used in 21 cm computations. This may
have no consequence in high-resolution simulations, as ionization fronts are
resolved, but becomes an issue as the resolution decreases. Low resolution
implies a high number of partially ionized particles that actually represent a
fully neutral and a fully ionized phase, and using the phase-averaged
temperature then leads to a poor estimation of the intensity of the 21 cm
signal. Indeed, the average temperature between the neutral and ionized gas is
greater than the temperature of the neutral phase alone, often by orders of
magnitudes. It does not affect codes that do not allow for partial ionization
such as 21cmFAST (Mesinger & Furlanetto, 2007), and post-processing solutions
have been designed to correct for this effect with the significant drawback of
having to run additional simulations without X-ray (Ross et al., 2016). To
compute the temperature of the neutral phase in Licorice, only the coupling to
the dynamics, the cosmological adiabatic expansion of the gas, and the heating
by X-ray were considered. This means computing the temperature evolution
equation in Baek et al. (2009) with only the HI atoms :

$\frac{dT_{HI}}{dt}=\frac{2}{3k_{B}n_{HI}}\left[-\frac{3}{2}k_{B}T_{HI}\frac{dn_{HI}}{dt}+\Lambda\right]$
(7)

where $\Lambda$ contains the X-ray heating and the adiabatic temperature
evolution of the gas due to the dynamics, as well as the viscosity term of the
SPH algorithm, from Monaghan (1992) and computed independently for the neutral
phase. It describes the evolution of the internal energy of the neutral phase
of a particle when interacting with the other particles in the SPH scheme.

The resulting evolution of temperatures is shown in figure 4. It confirms that
after the start of reionization ($z\lesssim 13$) the temperature of the
neutral phase of particles in LORRAH-21 is more than an order of magnitude
below their average temperature and on average significantly lower than the
temperature of weakly ionized ($\lesssim 2\%$) particles. It is much closer to
the temperature of weakly ionized particles of HIRRAH-21. This is used as an
imperfect proxy for the temperature of the neutral phase in HIRRAH-21, which
was not implemented at the time.

A decrease in resolution also affects the recombination rate $R$ of the
ionized gas, as it depends on the square of the HII density. The main
theoretical approach to counter this effect is to use a clumping factor.
Different versions of the clumping factor (inspired by different works :
Kaurov & Gnedin (2014); Mao et al. (2020); Chen et al. (2020); Bianco et al.
(2021)) have been implemented in Licorice. However, none of them managed to
reconcile the reionization timings of LORRAH-21 and HIRRAH-21.

The main cause of this failure is that any error in the calibration of the
parameters of the implemented clumping factor formula will cause an error in
the average ionized fraction in the simulated box that increases with time.
For instance, too small a clumping factor at a time step $n$ will cause too
few HII to recombine, leaving too many photons per HI at time step $n+1$,
further increasing the HII density. In conclusion, none of the theoretical
approaches to the clumping factor or fit to high-resolution simulations
available in the literature were accurate enough and matching the reionization
timings of our low and high-resolution simulations required a calibration of
the photon budget instead. To do so, a ”step” model of the UV escape fraction
$f_{esc}$ was implemented in LORRAH-21, depending on the local ionization
fraction $x_{ion}$ :

$\displaystyle f_{esc}=\cases{f}_{esc,post}$ ifxion>xthresh (8) (9)

$f_{esc,post}=0.1$, $f_{esc,pre}=0.01$, and $x_{thresh}=0.01$ are parameters
of this model, calibrated on HIRRAH-21 for LORRAH-21. This is justified by the
fact that the recombinations missing in LORRAH-21 occur in unresolved dense
regions that are likely to surround the unresolved sources and which must be
ionized before letting the UV photons escape the $\sim 1\,\rm{Mpc}$ region
represented by a gas particle. The ionization history of LORRAH-21 using this
model is consistent with the one of HIRRAH-21, something which was not
observed when the simulation code relied on a clumping factor implementation.
An approach using the clumping factor that leads to equivalent or more
accurate results might however be found in the future.

The global 21 cm signals of HIRRAH-21 and LORRAH-21 (calculated with the
average temperature and with the HI temperature) are shown in figure 5. As
expected, the 21 cm signal is strongly modified when calculated with the
phase-averaged temperature when compared to the result obtained using the HI
temperature. In addition, the 21 cm signal of LORRAH-21 calculated using the
HI temperature is in good agreement with the signal in HIRRAH-21. A small
mismatch occurs at $z\lesssim 11$, as heating occurs in LORRAH-21
approximately $0.15$ redshift units before HIRRAH-21, mainly due to the slight
difference in SFRD. Over the whole redshift range, LORRAH-21 is a good
approximation of HIRRAH-21, especially given the $10^{4}$ increase in
computation speed in LORRAH-21 caused by the drop in resolution : HIRRAH-21
required $\sim 3\times 10^{6}$ cpuh (hence why it was not run again with the
newly implemented subgrid models) while running LORRAH-21 only takes $\sim
3\times 10^{2}$ cpuh.

##  3 The LoReLi Database

We now present the LoReLi database111Available at https://21ssd.obspm.fr/,
which consists of 760 simulations with $256^{3}$ resolution elements run in
$300\,\rm{Mpc}$ boxes using the code presented in the previous section,
including the subgrid models for unresolved sources and the computation of the
neutral phase temperature for each particle. Initial conditions, generated
using the MUSIC code (Hahn & Abel, 2011), vary between simulations. The
4-parameters space sampled by the database was designed to loosely constrain
parameters value according to various probes of reionization. Here are the
varied parameters and explored ranges:

  * •

The gas-to-star conversion timescale $\tau\in[10\,Gyr,100\,Gyr]$ and minimum
halo mass $M_{min}\in[10^{8}\,\,\rm{M_{\odot}},4\times
10^{9}\,\,\rm{M_{\odot}}]$ are the parameters of the star formation model.

  * •

$f_{esc,post}\in\\{0.05,0.1,0.2,0.3,0.4,0.5\\}$ is the escape fraction of UV
radiation in particles with an ionized fraction higher than $x_{thresh}$.

  * •

The X-ray production efficiency $f_{x}\in[0.1,10]$, 9 logarithmically spaced
values. It drives the X-ray emissivity according to $L_{x}=3.4\times
10^{40}f_{x}\left(\frac{SFR}{1\,\rm{M_{\odot}}yr^{-}1}\right)\rm{ergs^{-1}}$
(Furlanetto et al., 2006).

Only astrophysical parameters were varied: cosmology was kept constant across
the database. Due to computation time constraints, $f_{esc,pre}$ and
$x_{thresh}$ were not varied and were kept at $0.01$, the same as in
LORRAH-21. We expect these parameters to have a smaller impact on the global
ionization field (except at high redshift) as the local ionization rises above
the threshold early in the source life. Numerical and cosmological parameters
are detailed in table 1.

$H_{0}$ | $\Omega_{0}$ | $\Omega_{b}$ | $\Omega_{\Lambda}$ | $\sigma_{8}$
---|---|---|---|---
67.8 km/s/Mpc | 0.308 | 0.0484 | 0.692 | 0.815

$L_{box}(\,\rm{Mpc})$ | dt (Myr) | $N_{part}$
---|---|---
294 | 7 | $256^{3}$

Table 1: Numerical and cosmological parameters of LoReLi simulations.

###  3.1 Observational constraints

Here we describe what regions of the parameters spaces were explored, the
observational constraints we used to restrict these regions, and present the
evolution of different relevant physical quantities in the LoReLi simulations.

####  3.1.1 Star formation rate

Figure 6: Star formation density history for the 17 ($\tau$, $M_{min}$)
parameters couples used in LoReLi simulations (solid). The parameters were
selected from a grid after a $\chi^{2}$ goodness-of-fit test. Only the
parameter couples yielding a p-value $p$ such that $1-p<0.95$ for at least one
of the plotted observational data sets were kept from the initial grid. For
reference, the SFR in HIRRAH-21 is also plotted (dashed).

The star formation parameters $\tau$ and $M_{min}$ were constrained using
recent observations of high redshift galaxies and SFR estimates by J. Bouwens
et al. (2016); McLeod et al. (2016); Oesch et al. (2018). A $\chi^{2}$ test
was performed to exclude models that fit none of the considered observational
data sets with a probability higher than 5 $\%$. . 17 $\\{\tau,M_{min}\\}$
couples were selected and the range of SFRD spanned in our database can be
seen in figure 6. JWST measurements are likely to give tighter and higher
redshift constraints in the near future.

####  3.1.2 The Thomson scattering optical depth

Figure 7: Thomson optical depth of the LoReLi signals. The ionization
histories of the simulations confine the optical depths to within 1$\sigma$ of
Planck’s 2015 measurement (dotted : central value, shaded region : 1$\sigma$).
As in 21SSD (Semelin et al., 2017), second ionization of He is assumed to
occur at $z=3$.

As late reionization scenarios are disfavoured by several reionization probes
(Fan et al., 2006; Mitra et al., 2011; Greig et al., 2019), $f_{esc}$ values
were chosen so that reionization in LoReLi ends between $z\sim 5$ and $z\sim
8$. Another characterization of reionization history is the Thomson optical
depth $\tau_{T}$. The values of $\tau_{T}$ of the LoReLi simulations are shown
in figure 7, as well as the results from the Planck collaboration (Ade et al.,
2016; White et al., 2018). While we did not use this observation to constraint
the parameter space of LoReLi, nearly all models are within $3\sigma$ of the
mean value in Planck 2018 cosmology and within $1\sigma$ of that of Planck
2015. We note that the value obtained in Planck 2018 assumes simple models for
reionization, and we do not exclude the couple LoReLi simulations that display
a $\gtrsim 3\sigma$ tension from the dataset.

Figure 8: Average ionization fraction of the LoReLi models, as well as
observations from Bouwens et al. (2015); Davies et al. (2018); Wang et al.
(2020). These observations were not used to calibrate the database, and
reionization end in all models between $z\sim 8$ and $z\sim 5$

No other constraints on reionization history were applied. In particular, no
attempt was made to match observational data on the evolution of the global
average ionized fraction $\langle x_{HII}\rangle$, as calibrating the escape
fraction parameters for each $\\{\tau,M_{min}\\}$ couple would require many
simulations. For completeness, we plot the evolution of $\langle
x_{HII}\rangle$ for all LoReLi simulations in figure 8.

###  3.2 21 cm signals in LoReLi

The differential brightness $\delta T_{b}$ of each 21 cm signal in LoReLi was
computed according to e.g. Furlanetto et al. (2006)

$\begin{split}\delta
T_{b}&=27x_{HI}(1+\delta)\left[\frac{T_{s}-T_{CMB}}{T_{s}}\right]\left[1+\frac{dv_{||}/dr_{||}}{H(z)}\right]^{-1}\\\
&\times\left[\frac{1+z}{10}\right]^{1/2}\left[\frac{\Omega_{b}}{0.044}\frac{h}{0.7}\right]\left[\frac{\Omega_{m}}{0.27}\right]^{1/2}\,\rm{mK}\end{split}$
(15)

here $x_{HI}$ is the local neutral fraction, $\delta$ the local overdensity,
$T_{s}$ the spin temperature of neutral hydrogen, $T_{CMB}$ the CMB
temperature at redshift $z$, $H(z)$ the Hubble parameter, and
$dv_{||}/dr_{||}$ the velocity gradient along the line of sight. The spin
temperature was calculated from simulated data according to the classical
equation (see e.g. Furlanetto et al. (2006)) :

$T_{s}^{-1}=\frac{T_{CMB}^{-1}+x_{c}T_{kin}^{-1}+x_{\alpha}T_{kin}^{-1}}{1+x_{c}+x_{\alpha}}$
(16)

where $T_{kin}$ is the kinetic temperature of hydrogen, $x_{c}$ and
$x_{\alpha}$ the collisional and Wouthuysen-Field couplings, respectively.
Collisional coupling is negligible in the considered redshift range, therefore
coupling to $T_{kin}$ is driven by the Wouthuysen-Field effect (Wouthuysen,
1952). $x_{\alpha}$ was computed on all saved snapshots for all simulations in
a post-processing step using the semi-analytical code SPINTER (Semelin et al.,
2023). The luminosities of Licorice particles were arranged on a $256^{3}$
grid, which SPINTER takes as input. It computes spherically symmetric
propagation kernels using MCMC ray-tracing that accounts for scattering in the
wings of Lyman-$\alpha$ line. These kernels are then convolved with the
emissivity field at previous redshifts to output grids of $x_{\alpha}$. The
outputs of SPINTER are close to those of full radiative transfer in the Lyman
bands, and calculating $x_{\alpha}$ for a single simulation takes a few cpuh,
making the cost of this post-processing step negligible. We show every global
21 cm signal in figure 9. As is, the database contains a wide range of models
that can be considered fiducial : they are shaped by conservative
observational constraints and no ”exotic” physics, such as non-standard dark
matter (Barkana, 2018) or extra radio background (Fialkov & Barkana, 2019),
was considered. In particular, the EDGES measurements are not compatible with
any of the LoReLi signals, as the absorption peaks of LoReLi signals are
shallower, at lower redshifts, and have milder slopes than EDGES-compatible
signals. However, some of these exotic physics, such as the extra radio
background, can easily be computed in post-processing steps from LoReLi
snapshots.

Figure 9: Global 21 cm brightness temperature of the LoReLi simulations. As
opposed to 21SSD, global brightness temperatures are calculated from full 3D
snapshots, not single slices of lightcones and therefore do not exhibit small
redshift-scale fluctuations. Color represents the value of $f_{X}$, the
parameter which impacts the depth of the signal the most.

##  4 Inference on power spectra

Figure 10: Power spectra of the LoReLi simulations at $k=0.1\,h/cMpc$ as a
function of $z$. They are compared to recent upper limits of various
instruments (Patil et al., 2017; Mertens et al., 2020; Trott et al., 2020;
Abdurashidova et al., 2023) as well as theoretical thermal noise power of
1000h observations with NenuFAR and LOFAR.

Figure 11: Power spectra of LoReLi models as a function of $k$ at redshift
9,11,14, from top to bottom. The value of the SKA thermal noise (1000h of
observation) is plotted to showcase that the LoReLi signals should be
detectable by the instrument.

In the previous section, we detailed the set of simulations at our disposal.
However, simulating many scenarios of the EoR is but the first step to
understand how reionization took place in our Universe. Indeed, one must then
be able to extract information from observations to determine which model is
the likeliest. In this section, we present how we perform this inference step
using the LoReLi dataset, first on mock data and then on actual observations
from the HERA instrument. LoReLi contains raw snapshots and lightcones
containing full particle information at 55 redshifts between 53.6 and 4.97, as
well as various physical quantities on 3D grids (data cubes), but in order to
compare models to future observational data, it is convenient to compress this
high dimensional data into summary statistics. The most standard choice is the
3D power spectrum of the 21 cm signal, which various instruments are currently
trying to measure (Patil et al., 2017; Mertens et al., 2020; Trott et al.,
2020; Abdurashidova et al., 2023). While it does not contain non-gaussian
information and summaries representing complementary information do exist, we
will focus on the power spectrum in the following. We show the power spectra
of the LoReLi simulations in figure 10 and in figure 11, as well as upper
limits from recent observations, and expected thermal noises of various
instruments. The power spectra were calculated using the coeval cubes of the
signal and not the lightcones. This approximation can have a significant (up
to $\sim 50\%$) effect on the spectrum (Datta et al., 2012, 2014) and taking
it into account will be at the centre of future improvements of our method. As
expected from classical 21cm theory, at higher redshift, the amplitude of the
spectra is mostly driven by the Lyman-$\alpha$ coupling, and becomes primarily
dictated by the value of $f_{X}$ at lower redshift. We anticipate that
improvements of these upper limits in the coming years will soon allow tight
constraints on our models.

Figure 12: Diagnostics of the performance of the emulator. Top left : at each
$k,z$, median across the dataset of the emulator error to thermal noise ratio.
Top right : at each $k,z$, median across the dataset of the emulator error to
cosmic variance ratio.. The error increases with cosmic variance, as the
deterministic network cannot replicate the noise induced by cosmic variance.
Bottom left :  100x mean squared relative error across the dataset between the
emulator prediction and the training sample. Bottom right :  loss at each k,
z. Faint signal at high z causes low loss but high relative error. Cosmic
variance causes high loss and relative error at low k. Additionally, we
checked that the mean square relative error exhibits no systematic trend
depending on the values of the astrophysical parameters.

Figure 13: Randomly selected examples of true (top) and emulated (bottom)
LoReLi signals qualitatively showing the resemblance between emulated and
simulated spectra

.

Figure 14: Results of multiple inferences done using several versions of
LorEMU, trained with different weights initializations (dashed). On the 2D
panels, only the 1$\sigma$ contours are shown, for clarity. The blue dotted
lines are the posteriors obtained using each emulator. The thick red line show
the posterior obtained using the average emulator on the noised signal, while
the inference on the noiseless signal is shown in black.  Figure 15: Results
of inference on HERA data using LorEMU. On the 1D panels, the priors are in
red and the posteriors in blue. Contours represent $68\%$ and $95\%$ levels of
confidence.

Classical Bayesian inference, for example through the use of MCMC, requires
large numbers of forward modelling instances ($\gtrsim 10^{5}$), often
sequentially. Given the computation cost of a single Licorice simulation, this
is unreasonable even at $256^{3}$ resolution. In order to perform parameter
inference on this power spectra dataset, our method consists in using LoReLi
as a training sample for a neural network that will function as an emulator of
Licorice, and then perform classical MCMC inference using the emulator as the
model. Indeed, once the emulator is trained, the computation cost of producing
a single signal drops to a few ms, which allows for a sufficiently large
number of steps in the Markov chain in a few hours of runtime.

###  4.1 Data preprocessing

####  4.1.1 Data properties

In order to efficiently train the network, it is necessary to preprocess the
data. The power spectra were computed for a list of 32 redshifts, with power
spectra being set to zero for $z$ between the redshift of full reionization
and the last redshift bin. However, in order to reduce the dimensionality of
the problem, we restrict our analysis to
$k=[0.23,0.33,0.46,0.66,0.93,1.31,1.86,2.64\rm{h/cMpc}]$ and
$z=[15.67,14.05,12.77,11.73,10.87,10.14,9.51,8.96]$. The $k$ bins width
defined such that the bins are non-overlapping. 21 cm power spectra span many
orders of magnitudes, which tends to hinder the training of networks.
Therefore, the logarithm of the power spectra was used before normalizing by
dividing by the maximum value across the dataset.

####  4.1.2 Noise in the signals

For an inference framework to have any practical application, it must be
performed on signals affected by instrumental effects. We will here neglect
foreground residuals and other systematics and focus on the SKA thermal noise
for an observation time of 100h. Like in Doussot et al. (2019), and following
McQuinn et al. (2006), the standard deviation of the power spectrum caused by
SKA thermal noise was computed as such :

$\delta
P(k,z)=\left[\sum_{|k|=k}\left(\frac{1}{\frac{Ax^{2}y}{\lambda(z)^{2}B^{2}}C(\textbf{k},\textbf{k})}\right)^{2}\right]^{-1/2}$
(17)

where $\lambda$ is the observed wavelength, $A$ the area of a 256-antenna
station, $B=10\,\rm{MHz}$ is the bandwidth, $x$ the comoving distance to the
observed redshift, $y$ the depth of field, and $C$ the detector covariance
matrix. In order to noise a signal, we simply add for each $k,z$ a realization
of a Gaussian random variable of mean 0 and standard deviation $\delta
P(k,z)$.

This noise is added to another form of stochasticity, the cosmic variance of
the simulation box (CV), that is implicit to our setup. We recall the equation
found in McQuinn et al. (2006) to compute the covariance matrix of CV, which
is added to the detector noise covariance matrix in the previous equation :

$C_{CV}(k_{i},k_{j})\approx
P(k_{i})\frac{\lambda^{2}B^{2}}{Ax^{2}y}\delta_{ij}$ (18)

Marginalizing over CV for each point in parameter space requires hundreds of
simulations at each point and is not feasible using LoReLi. However, as CV
mainly affects the largest scales, its effect can be at least partially
avoided by focusing analysis on $k\gtrsim 0.1\,\rm{h/cMpc}$. This motivates us
to target the analysis towards a region in $k$-space more affected by thermal
noise, which is far less computationally expensive to generate. It is also
likely that a machine learning method that tends to interpolate between the
LoReLi simulations, such as the one presented in the very next section,
smooths over the effects of CV altogether.

###  4.2  LorEMU

Layer type | Activation function | Regularization
---|---|---
Dense , 512 neurons | Leaky ReLu ($\alpha=0.05$) | L1 ($\lambda=5\times 10^{-6}$)
Dense, 512 neurons | Leaky ReLu ($\alpha=0.05$) | L1 ($\lambda=5\times 10^{-8}$)
Dense, 512 neurons | Leaky ReLu ($\alpha=0.05$) | None
Dense, 1 neuron | Sigmoid | None
Optimizer | Loss function | Batch size
Adam (learning rate = $5\times 10^{-4})$ | Mean Squared Error | 32
Table 2: Architecture and hyper parameters of LorEMU.

LorEMU, an emulator of the Licorice power spectra, was trained using the
LoReLi simulations. The network is a multilayer perceptron that takes the 4
astrophysical parameters and a (k, z) pair as inputs and outputs a noiseless
power spectrum value. The full architecture and parameters are shown in
table222The detailed signification of the technical terms can be found at
https://keras.io/ 2 and was implemented using the Keras framework (Chollet,
2015). It was trained in a supervised manner to predict the value of the
(noiseless) power spectrum given the inputs $k,z$, and the astrophysical
parameters values. Training took place over 200 epochs, using a batch size of
32, and minimizing a Mean Squared Error loss function using the Adam optimizer
with a learning rate of $5\times 10^{-4}$. The training set was composed of
$75\%$ of the 760 spectra, randomly selected, while the remaining $25\%$ of
the spectra were set aside to constitute the test set.

This network is deterministic. In particular, this means the output power
spectra do not include a proper contribution by CV. It is hypothesized that
the emulator learns to average over the CV affecting the large scales, as the
spectra in the training set that are close in parameter space have been
generated using different initial conditions of the density field and
therefore have different CV realizations. Various quantities describing the
accuracy of the emulator are depicted in figure 12. As indicated by the bottom
left plot, the emulator performs especially well at high $k$ and low
redshifts. At large scales, the error is caused by the fact that the emulator
cannot predict CV whereas spectra in the test set include a CV realization
(see top right plot). At high redshift, the signal is faint, resulting in low
loss but high relative error (see bottom plots). The rise in relative error at
the very lowest redshift may partially be due to simulations approaching the
end of reionization, but is in this case accompanied by a rise in the loss
function that is more difficult to explain. A similar issue occurred in
Jennings et al. (2019), and this region in (k,z) may be the focus of future
improvements of the emulator. Overall, the emulator outperforms the one in
Jennings et al. (2019) by a factor of $\sim 2$. We show examples of simulated
and emulated signals in figure 13.

####  4.2.1 Inference on mock data

LorEMU allows for a generation of a $8\times 8$ power spectrum in milliseconds
and is therefore a suitable method of forward modelling in an MCMC inference
pipeline. We show examples of Bayesian inference in figure 14, on a simulated
spectrum 333 calculated using the following astrophysical parameters :
$f_{X}=0.56$, $\tau=27\,Gyr$, $M_{min}=10^{9.06}\,\rm{M_{\odot}}$,
$f_{esc,post}=0.2$ assuming perfect foreground removal and noise corresponding
to a 100h SKA observation, as explained previously. The inferences were
performed using the emcee python package (Foreman-Mackey et al., 2013). The
log-likelihood of a model $y$ with astrophysical parameter set $\theta$, with
respect to data $x$, and with total variance $\sigma_{tot}$, can be explicitly
written :

$logL(y|\theta)=-\sum_{k,z}\frac{1}{2}\left(\frac{x(k,z)-y(\theta,k,z)}{\sigma_{tot}(k,z)}\right)^{2}$
(26)

It is worth noting that this form of the likelihood assumes that Fourier modes
are independent variables, an assumption discussed in e.g Prelogović &
Mesinger (2023). While the authors of this paper show that using a diagonal
covariance for the cosmic variance can lead to biased inference results, they
also show that the non-diagonal terms are mostly erased by the level of their
thermal noise, which is an order of magnitude below the one used in our work.
Therefore, we expect no significant contribution from these non-diagonal terms
in our case. In our case, $\sigma_{tot}$ is

$\sigma_{tot}=\sqrt{\sigma_{thermal}^{2}+\sigma_{CV}^{2}+\sigma_{train}^{2}}.$
(27)

It includes at each bin of $k,z$ the cosmic variance $\sigma_{CV}$, the
variance of the SKA thermal noise $\sigma_{thermal}$, and the ”training
variance” $\sigma_{train}$, which represents the variance in the training of
the network and is detailed thereafter.

According to Bayes theorem, the posterior distribution of the parameters can
be written :

$p(\theta|y)\propto L(y|\theta)\pi(\theta)$ (28)

Where $\pi(\theta)$ is the prior over the astrophysical parameters.

In order to estimate the error induced by an imperfect training of the
emulator, 9 different versions of LorEMU were trained. Each version uses the
same architecture, and was trained using different weights initialization and
with different random splits of the data to constitute the test and training
sets. The rest of the training pipeline was kept unchanged. Each one of them
was subsequently used as the forward simulator in the MCMC pipeline, producing
the different posteriors of figure 14. The inferences were done using 160
walkers, randomly initialized in the prior, and $\sim 800000$ total steps. The
prior on the astrophysical parameters is flat within the region in parameter
space explored in LoReLi and zero outside this region.

The inferences were all performed on the same simulated signal noised with the
same 100h SKA noise realisation, therefore the differences between inferences
are due to the different weight initializations and training stochasticity
between the different versions of LorEMU. This training variance can be
defined at each $k,z$, and astrophysical parameters set $\theta$ as

$\sigma_{train}^{2}=\frac{1}{N}\sum_{N}\left(P_{predicted}(k,z,\theta)-P_{target}(k,z,\theta)\right)^{2}$
(29)

where $P_{target}$ and $P_{predicted}$ are the training spectra and the
outputs of the emulator respectively, and the sum is taken over a large number
$N$ of emulators trained with different weight initializations. Rigorously
evaluating this variance would not only require a large number of emulators,
but also a target signal for any value of $\theta$, while we have signals only
for the parameter values present in the LoReLi database. To address those
difficulties we chose to average this quantity over $\theta$. The result is an
approximation of the training variance. However, it is justified by the fact
that the error shows no clear dependency on parameters. This variance adds a
bias in the inference with single emulators, and putting the training variance
in the likelihood widens the confidence contours. While negligible for
$f_{X}$, the effect of this training variance compares to that of the SKA 100h
thermal noise for the SFR parameters. As an attempt to marginalize over the
random initializations, figure 14 shows a posterior obtained using the average
predictions of the 9 emulators in the MCMC framework, following the bagging
method common in machine learning. In order to demonstrate that, for the
average emulator, the biases of the predictions are negligible compared to the
uncertainty induced by the noise, figure 14 also shows the posterior applied
to the noiseless signal (but with the SKA 100h thermal noise variance still
included in the likelihood). We see that in this case the confidence contours
are well centred in the target parameters values. This also indicates that our
inference results are only weakly affected by our choice of using a diagonal
covariance in the likelihood.

$\tau$ and $M_{min}$ are unsurprisingly strongly anticorrelated, as their
effects on star formation are degenerate. Similarly, $f_{X}$ is correlated
with $M_{min}$ and anticorrelated with $\tau$, as X-ray emissivity is
proportional to SFR. We note that these effects are unlikely to stem from our
choice of priors, since they are visible in confidence contours far away from
the flat prior boundaries. Additionally, $f_{esc}$ is the least strongly
constrained parameter, with large contours and non-negligible likelihood over
the whole prior. Except for $f_{esc}$, the true values of the astrophysical
parameters are within the $1\sigma$ contours obtained after inference with the
average emulator.

####  4.2.2 Inference on recent HERA data

In the previous section, we validated our method by applying it to mock
observations generated from our simulated data. We now turn to real data, and
provide an analysis of the most recent HERA observations (Abdurashidova et
al., 2023), at $z=10.4$ and $z=7.9$. As indicated on figure 11, the most
recent HERA upper limits are currently the only ones susceptible to constrain
our data set. We follow the procedure detailed in Abdurashidova et al. (2022)
(thereafter HERA2022), which we briefly summarize here. The observed data is
assumed to be the 21 cm signal, to which systematics (typically foreground
residuals and radio frequency interference) are added. The systematics power
is supposed to be positive, and marginalizing over the systematics yields the
likelihood :

$L(y|\theta)=\prod_{k,z}\frac{1}{2}\left[1+\rm{erf}\left(\frac{x(k,z)-y(\theta,k,z)}{\sqrt{2}\sigma(k,z)}\right)\right].$
(30)

The emulator was not retrained on spectra calculated with the $k$ bins of HERA
data with its window function. The posterior distribution of the astrophysical
parameters obtained through MCMC inference is shown in figure 15. The only
parameter for which the posterior distribution is different from its prior is
$f_{X}$. Qualitatively, our conclusions match those of HERA2022 : a cold
reionization of the Universe is an unlikely scenario. Quantitatively, however,
the constraints we put on $f_{X}$ appear looser than in HERA2022. In figure 6
of HERA2022, which shows the inference results as obtained using 21cmFAST, the
likelihood decreases for $f_{x}<1$, while we only observe such decrease for
$f_{X}<0.5$. This may stem from a difference in the prior on the SFR
parameters. Indeed, the observational constraints on SFR applied to LoReLi and
21cmFAST only exist at $z\leq 10$, leaving the functional form of the SFR
evolution completely free at higher redshift. Different posteriors can result
from this, as seminumerical codes typically assume that stars are formed
instantaneously as some fraction of the collapsed gas. As detailed in section
2, in Licorice, the stellar fraction evolves with time. If the stellar
fraction in Licorice is typically lower than in seminumerical codes at high
redshift and higher at lower redshift in order to be consistent with the same
observational constraints on SFRD, then this difference in the inference
result is expected. Unfortunately, as the codes function in significantly
different ways with different parametrizations, understanding the difference
in the predicted posteriors of both inferences would require an in-depth
comparison of Licorice and 21cmFAST (as well as with other codes) that does
not exist yet and lies beyond the scope of this work. Regarding figure 14 of
HERA2022, obtained using the seminumerical code described in HERA2022, an
additional explanation may be that the authors marginalize over $f_{r}$, a
parameter that controls the intensity of an exotic radio background, while we
implicitly fix $f_{r}$ to 0 as the only background we consider is the CMB.

As a final note, the inference was also attempted on upper limits from LOFAR
(Mertens et al., 2020) using the exact same methodology. As could be expected
from figure 11, no difference between priors and posteriors could be
established in that case.

##  5 Conclusions

In this work, we have presented new functionalities of the Licorice simulation
code that allow low resolution simulations of the EoR to reproduce many
features of high resolution simulations. They are twofold. One is an
implementation of the conditional mass function formalism as a way to
statistically estimate the mass of unresolved dark matter halos. We then
include the gas contained in these halos to compute star formation in low
resolution simulations, which would otherwise be severely underestimated.
Using the Sheth-Tormen formulation of the CMF, we find excellent agreement
with $1024^{3}$ and $2048^{3}$ simulations.

The other modification is a calculation of the temperature of the neutral
phase of each particle. Due to poorly resolved ionization fronts in low
resolution simulations, we found that using the phase-averaged temperature of
the particles leads to a significant overestimate of the intensity of the
21-cm signal. Using the neutral phase temperature produces a good agreement
between results at low and high resolution.

Together, these improvements allow for physically reasonable Licorice
simulations at $256^{3}$ resolution, that run in approximately $300$ cpuh
each. Consequently, this makes running many Licorice simulations
computationally feasible. Therefore, we present LoReLi, a growing dataset of
760 Licorice simulations with $256^{3}$ particles in $300\,\rm{Mpc}$ boxes.
Full particle data is saved at 32 redshift values, and LoReLi spans a
4-parameter space : the escape fraction $f_{esc,post}$, the X-ray emissivity
$f_{X}$, the minimum mass of star forming halos $M_{min}$, and the gas-to-star
conversion timescale $\tau$. The latter two parameters were calibrated using
constraints on the SFRD, and the first one was chosen so that reionization
ends in all models between $z\sim 5$ and $z\sim 8$. LoReLi therefore contains
a variety of standard models of the Epoch of Reionization.

In the first application of this dataset, we summarize our data into
independent power spectra values at 8 $k$ and 8 $z$ values. We then present
LorEMU, a neural network trained on LoReLi to produce accurate Licorice power
spectra given $k,z$, and a set of the 4 astrophysical parameters. During
training, LorEMU reaches $5\sim 10\%$ relative mean squared error averaged
over the dataset. We then perform Bayesian inference using LorEMU and MCMC on
a power spectrum generated using Licorice, onto which a noise realization
corresponding to 100h of SKA observations was added. Because different random
weights initializations affect the result of the network training, 9 different
versions of LorEMU were trained, and training variance was included into the
inference likelihood. By averaging the outputs of the different networks
during inference, we obtain accurate posterior distributions. We find that
this approach produces a posterior with very little bias by performing
inference on a noiseless simulated signal. Finally, we apply the same
inference pipeline to the most recent HERA upper limits and obtain constraints
on lower values of $f_{X}$, indicating that our Universe is unlikely to have
followed a cold reionization scenario.

The new implementation of Licorice suffers from some limitations. While star
formation in resolved halos is sensitive to the local environment through the
dynamics of the gas, unresolved halos of a given mass have similar star
formation rate, that only vary if they are located in particles with different
stellar fractions. Moreover, the CMF formalism provides the cosmic average of
the halo population in a given region with fixed overdensity. In reality, such
regions may present fluctuations around this average. A future improvement
will be to include stochasticity in the unresolved star formation rate, both
from the halo population fluctuations and from the star formation efficiency
of each unresolved halo. The latter will have to be parameterized with a last
one additional astrophysical parameter. Future versions of the LoReLi database
will include both a denser sampling of the parameter space and variations of
additional astrophysical parameters.

In this work, we have only begun to tap the potential of the LoReLi database.
Since we are modelling the 21 cm signal using full non-linear dynamics and 3D
radiative transfer, we can hope that the non-Gaussian properties of the signal
are well accounted for. It will be relevant to explore inference based on non-
Gaussian summary statistics of the signal. An explicit likelihood is usually
not available in this case. One solution is to choose summary statistics that,
while encoding non-Gaussian properties of the signal, have themselves near-
Gaussian behaviours (due to the central limit theorem) and model the
likelihood as a multivariate Gaussian. The other promising avenue is to use
Simulation Based Inference (e.g Zhao et al. (2022)). We will be working toward
both goals.

###### Acknowledgements.

This project was provided with computer and storage resources by GENCI at TGCC
thanks to the grant 2023-A0130413759 on the supercomputer Joliot Curie’s ROME
partition. This study was granted access to the HPC resources of MesoPSL
financed by the Région Île-de-France and the project EquipMeso (reference
ANR-10-EQPX-29-01) of the programme Investissements d’Avenir supervised by the
Agence Nationale pour la Recherche.

## References

  * Abdurashidova et al. (2023) Abdurashidova, T. H. C. Z., Adams, T., Aguirre, J. E., et al. 2023, The Astrophysical Journal, 945, 124
  * Abdurashidova et al. (2022) Abdurashidova, Z., Aguirre, J. E., Alexander, P., et al. 2022, The Astrophysical Journal, 924, 51
  * Ade et al. (2016) Ade, P. A., Aghanim, N., Arnaud, M., et al. 2016, Astronomy and Astrophysics, 594 [arXiv:1502.01589]
  * Alsing & Wandelt (2018) Alsing, J. & Wandelt, B. 2018, Monthly Notices of the Royal Astronomical Society: Letters, 476, L60
  * Alsing et al. (2017) Alsing, J., Wandelt, B., & Feeney, S. 2017, Massive optimal data compression and density estimation for scalable, likelihood-free inference in cosmology, Tech. rep.
  * Baek et al. (2009) Baek, S., Di Matteo, P., Semelin, B., Combes, F., & Revaz, Y. 2009, Astronomy and Astrophysics, 495, 389
  * Baek et al. (2010) Baek, S., Semelin, B., Di Matteo, P., Revaz, Y., & Combes, F. 2010, Astronomy and Astrophysics, 523 [arXiv:1003.0834]
  * Barkana (2018) Barkana, R. 2018, Nature, 555, 71
  * Bianco et al. (2021) Bianco, M., Iliev, I. T., Ahn, K., et al. 2021, Monthly Notices of the Royal Astronomical Society, 504, 2443
  * Bouwens et al. (2015) Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2015, Astrophysical Journal, 803, 1
  * Bowman et al. (2018) Bowman, J. D., Rogers, A. E., Monsalve, R. A., Mozdzen, T. J., & Mahesh, N. 2018, Nature, 555, 67
  * Chen et al. (2020) Chen, N., Doussot, A., Trac, H., & Cen, R. 2020, The Astrophysical Journal, 905, 132
  * Chollet (2015) Chollet, F. 2015, Keras
  * Ciardi et al. (2003) Ciardi, B., Stoehr, F., & White, S. D. M. 2003, Simulating intergalactic medium reionization, Tech. rep.
  * Cohen et al. (2017) Cohen, A., Fialkov, A., Barkana, R., & Lotem, M. 2017, Monthly Notices of the Royal Astronomical Society, 472, 1915
  * Datta et al. (2014) Datta, K. K., Jensen, H., Majumdar, S., et al. 2014, Monthly Notices of the Royal Astronomical Society, 442, 1491
  * Datta et al. (2012) Datta, K. K., Mellema, G., Mao, Y., et al. 2012, Monthly Notices of the Royal Astronomical Society, 424, 1877
  * Davies et al. (2018) Davies, F. B., Hennawi, J. F., Bañados, E., et al. 2018, The Astrophysical Journal, 864, 142
  * Doussot et al. (2019) Doussot, A., Eames, E., & Semelin, B. 2019, Monthly Notices of the Royal Astronomical Society, 490, 371
  * Doussot & Semelin (2022) Doussot, A. & Semelin, B. 2022, Astronomy and Astrophysics, 667, A118
  * Fan et al. (2006) Fan, X., Carilli, C. L., & Keating, B. 2006, Annual Review of Astronomy and Astrophysics, 44, 415
  * Fialkov & Barkana (2019) Fialkov, A. & Barkana, R. 2019, Monthly Notices of the Royal Astronomical Society, 486, 1763
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, Publications of the Astronomical Society of the Pacific, 125, 306
  * Furlanetto et al. (2006) Furlanetto, S. R., Peng Oh, S., & Briggs, F. H. 2006, Cosmology at low frequencies: The 21 cm transition and the high-redshift Universe
  * Ghara et al. (2018) Ghara, R., Mellema, G., Giri, S. K., et al. 2018, MNRAS, 476, 1741
  * Greig & Mesinger (2015) Greig, B. & Mesinger, A. 2015, Monthly Notices of the Royal Astronomical Society, 449, 4246
  * Greig & Mesinger (2018) Greig, B. & Mesinger, A. 2018, Monthly Notices of the Royal Astronomical Society, 477, 3217
  * Greig et al. (2019) Greig, B., Mesinger, A., & Bañados, E. 2019, Monthly Notices of the Royal Astronomical Society, 484, 5094
  * Hahn & Abel (2011) Hahn, O. & Abel, T. 2011, Mon. Not. R. Astron. Soc, 415, 2101
  * J. Bouwens et al. (2016) J. Bouwens, R., Aravena, M., Decarli, R., et al. 2016, The Astrophysical Journal, 833, 72
  * Jennings et al. (2019) Jennings, W. D., Watkinson, C. A., Abdalla, F. B., & Mcewen, J. D. 2019, MNRAS, 483, 2907
  * Kaurov & Gnedin (2014) Kaurov, A. A. & Gnedin, N. Y. 2014, The Astrophysical Journal, 787, 146
  * Kendall & Gal (2017) Kendall, A. & Gal, Y. 2017, in Advances in Neural Information Processing Systems, Vol. 2017-Decem, 5575–5585
  * Krause et al. (2018) Krause, F., Thomas, R. M., Zaroubi, S., & Abdalla, F. B. 2018, New Astronomy, 64, 9
  * Lacey & Cole (1993) Lacey, C. & Cole, S. 1993, Monthly Notices of the Royal Astronomical Society, 262
  * Lewis et al. (2022) Lewis, J. S., Ocvirk, P., Sorce, J. G., et al. 2022, Monthly Notices of the Royal Astronomical Society, 516, 3389
  * Mao et al. (2020) Mao, Y., Koda, J., Shapiro, P. R., et al. 2020, MNRAS, 491, 1600
  * McLeod et al. (2016) McLeod, D. J., McLure, R. J., & Dunlop, J. S. 2016, Monthly Notices of the Royal Astronomical Society, 459, 3812
  * McQuinn et al. (2006) McQuinn, M., Zahn, O., Zaldarriaga, M., Hernquist, L., & Furlanetto, S. R. 2006, The Astrophysical Journal, 653, 815
  * Mellema et al. (2006) Mellema, G., Iliev, I. T., Alvarez, M. A., & Shapiro, P. R. 2006, New Astronomy, 11, 374
  * Mertens et al. (2020) Mertens, F. G., Mevius, M., Koopmans, L. V. E., et al. 2020, MNRAS, 493, 1662
  * Mesinger & Furlanetto (2007) Mesinger, A. & Furlanetto, S. 2007, The Astrophysical Journal, 669, 663
  * Mesinger et al. (2011) Mesinger, A., Furlanetto, S., & Cen, R. 2011, Monthly Notices of the Royal Astronomical Society, 411, 955
  * Mitra et al. (2011) Mitra, S., Choudhury, T. R., & Ferrara, A. 2011, Mon. Not. R. Astron. Soc, 413, 1569
  * Monaghan (1992) Monaghan, J. 1992, in Annu.Rev.Astron.Astrophys, Vol. 97, 101–123
  * Murray et al. (2020) Murray, S., Greig, B., Mesinger, A., et al. 2020, Journal of Open Source Software, 5, 2582
  * Neutsch et al. (2022) Neutsch, S., Heneka, C., & Brüggen, M. 2022, Monthly Notices of the Royal Astronomical Society, 511, 3446
  * Ocvirk et al. (2020) Ocvirk, P., Aubert, D., Sorce, J. G., et al. 2020, Monthly Notices of the Royal Astronomical Society, 496, 4087
  * Oesch et al. (2018) Oesch, P. A., Bouwens, R. J., Illingworth, G. D., Labbé, I., & Stefanon, M. 2018, The Astrophysical Journal, 855, 105
  * Patil et al. (2017) Patil, A. H., Yatawatta, S., Koopmans, L. V. E., et al. 2017, The Astrophysical Journal, 838, 65
  * Prelogović & Mesinger (2023) Prelogović, D. & Mesinger, A. 2023, Monthly Notices of the Royal Astronomical Society
  * Reed et al. (2007) Reed, D. S., Bower, R., Frenk, C. S., Jenkins, A., & Theuns, T. 2007, Mon. Not. R. Astron. Soc, 374, 2
  * Reis et al. (2021) Reis, I., Fialkov, A., & Barkana, R. 2021, The subtlety of Ly-$\alpha$ photons: changing the expected range of the 21-cm signal, Tech. rep.
  * Ross et al. (2016) Ross, H. E., Dixon, K. L., Iliev, I. T., & Mellema, G. 2016, Simulating the Impact of X-ray Heating during the Cosmic Dawn, Tech. rep.
  * Rubiño-Martín et al. (2008) Rubiño-Martín, J. A., Betancort-Rijo, J., & Patiri, S. G. 2008, Monthly Notices of the Royal Astronomical Society, 386, 2181
  * Santos et al. (2010) Santos, M. G., Ferramacho, L., Silva, M. B., Amblard, A., & Cooray, A. 2010, Monthly Notices of the Royal Astronomical Society, 406, 2421
  * Schaeffer et al. (2023) Schaeffer, T., Giri, S. K., & Schneider, A. 2023, MNRAS, 000, 1
  * Schmit & Pritchard (2018) Schmit, C. J. & Pritchard, J. R. 2018, Monthly Notices of the Royal Astronomical Society, 475, 1213
  * Semelin (2015) Semelin, B. 2015, Detailed modelling of the 21-cm Forest, Tech. rep.
  * Semelin et al. (2007) Semelin, B., Combes, F., & Baek, S. 2007, Astronomy and Astrophysics, 474, 365
  * Semelin et al. (2017) Semelin, B., Eames, E., Bolgar, F., & Caillat, M. 2017, Monthly Notices of the Royal Astronomical Society, 472, 4508
  * Semelin et al. (2023) Semelin, B., Mériot, R., Mertens, F., et al. 2023, Astronomy and Astrophysics, 672, 162
  * Sheth et al. (2001) Sheth, R. K., Mo, H. J., & Tormen, G. 2001, Monthly Notices of the Royal Astronomical Society, 323, 1
  * Shimabukuro & Semelin (2017) Shimabukuro, H. & Semelin, B. 2017, MNRAS, 468, 3869
  * Singh et al. (2022) Singh, S., Nambissan T, J., Subrahmanyan, R., et al. 2022, Nature Astronomy, 6, 607
  * Springel et al. (2001) Springel, V., Yoshida, N., & White, S. D. 2001, New Astronomy, 6, 79
  * Tramonte et al. (2017) Tramonte, D., Rubiño-Martín, J. A., Betancort-Rijo, J., & Dalla Vecchia, C. 2017, Monthly Notices of the Royal Astronomical Society, 467, 3424
  * Trott et al. (2020) Trott, C. M., Jordan, C. H., Midgley, S., et al. 2020, Monthly Notices of the Royal Astronomical Society, 493, 4711
  * Wang et al. (2020) Wang, F., Davies, F. B., Yang, J., et al. 2020, The Astrophysical Journal, 896, 23
  * White et al. (2018) White, D. M., Zacchei, A., & Zonca, A. 2018, Planck 2018 : Diego 60 , O. Doré 62,10 , M. Douspis 53 , A. Ducout 54,52 , X. Dupac 35 , S. Dusini 61 , G. Efstathiou 65,57 * , F. Elsner 73 , T. A. Enßlin 73 , H. K. Eriksen 58 , Y. Fantaye 3,18 , M. Farhang 77 , J. Fergusson 11 , R. Fernandez-Cobos 60 ,, Tech. rep.
  * Wouthuysen (1952) Wouthuysen, S. A. 1952, Physica, 18, 75
  * Zhao et al. (2022) Zhao, X., Mao, Y., Cheng, C., & Wandelt, B. D. 2022, The Astrophysical Journal, 926, 151