# Photon and dilepton production in high energy heavy ion collisions

Takao Sakaguchi<EMAIL_ADDRESS>Physics Department, Brookhaven National
Laboratory, Upton, NY 11973-5000, USA

###### Abstract

The recent results on direct photons and dileptons in high energy heavy ion
collisions, obtained particularly at RHIC and LHC are reviewed. The results
are new not only in terms of the probes, but also in terms of the precision.
We will discuss the physics learned from the results.

###### keywords:

Photons, Dileptons, RHIC, LHC, QGP

###### pacs:

25.75.-q, 25.75.Bh, 25.75.Cj, 25.75.Nq

##  1 Introduction

Electromagnetic radiations are an excellent probe for extracting
thermodynamical information of a matter produced in nucleus-nucleus
collisions, as they are emitted from all the stages of collisions with wide
$Q^{2}$, and don’t interact strongly with medium once produced [1]. They
appear in two figures, namely, photons ($\gamma$) that have zero-mass and
virtual photons ($\gamma^{*}$) that have finite mass. Experimentalists often
refer virtual photons as dileptons since they have been measuring virtual
photons through lepton-pair channels ($\gamma^{*}\rightarrow ee,\mu\mu$), by
which a wide kinematic range in invariant mass and transverse momentum
($p_{T}$) can be explored. Photons are produced through a Compton scattering
of quarks and gluons ($qg\rightarrow q\gamma$) or an annihilation of quarks
and anti-quarks ($q\overline{q}\rightarrow g\gamma$) as leading-order (LO)
processes, and the next-to-leading order (NLO) process is dominated by
Bremsstrahlung and fragmentation ($qg\rightarrow qg\gamma$), as depicted in
Figure 1.

Figure 1: Production processes of direct photons.

Their yields are proportional to $\alpha\alpha_{s}$, which are $\sim$40 times
lower than those of hadrons that are produced in strong interactions. Virtual
photons are produced by annihilation of quarks and anti-quarks. The initial
hard scattering process that takes place in the beginning of the collisions
produces relatively high $p_{T}$ photons, often referred as hard photons. They
have same LO and NLO processes, and are called prompt photons and fragment
photons, respectively. The production rate of these photons are well described
by a NLO pQCD calculation [2]. Photons will be emitted from the hot and dense
medium (quark-gluon plasma: QGP) in high energy nucleus collisions and
manifest at moderate $p_{T}$(1$<p_{T}<$3 GeV/$c$) if the QGP is formed [3]. We
often call these photons as thermal photons. The thermal photons are of
interest in exploring thermodynamical nature of the QGP, such as temperature.
One can also obtain the degree of freedom of the system by combining the
temperature with measurement of the energy density of the system as
$g\propto\epsilon/T^{4}$. For $p_{T}<$1 GeV, the photons are predominantly
contributed from hadron gas state via the processes of
$\pi\pi(\rho)\rightarrow\gamma\rho(\pi)$, $\pi K^{*}\rightarrow K\gamma$ and
etc., which are no longer the quark-gluon level interaction. We often refer
these photons as hadron-gas photons. Photons from Compton scattering of hard-
scattered partons and partons in the medium (jet-photon conversion), or
Bremsstrahlung of the hard scattered partons in the medium have been predicted
to contribute in the $p_{T}$ range of $p_{T}>$2.5 GeV/$c$ if QGP is created
[4]. Figure 2 shows a landscape of photon sources as a function of formation
time and $p_{T}$.

Figure 2: Manifestation of photons from various sources as a function of
formation time and $p_{T}$.

As described above, photons have rich information on the state they are
emitted. However these photons are overwhelmed by the large background coming
from the decay of hadrons. $\pi^{0}\rightarrow\gamma\gamma$ and
$\eta\rightarrow\gamma\gamma$ are the major contributors to the background
photons ($\sim$95 %). The signal to background ratio at 1-3 GeV/$c$ is of the
order $\sim$10 %. One can understand the difficulty of the measurement by
comparing the uncertainty of the best $\pi^{0}$ measurement at RHIC [5] which
directly relates to the precision of background determination. There is a
certain probability that photons produced by the same process acquire virtual
mass and decay into lepton-pairs (shown as one another degree of freedom in
Figure 1(b)). This process is called internal conversion process and is
different from the virtual photon production described above. We will explain
these photons in detail in a later section. There have been several attempts
on measuring thermal photons in high energy nuclear collisions as an evidence
of QGP formation. The first sizable signal was reported by the WA98
experiments at SPS in 1$<p_{T}<$3 GeV/$c$ where a calculation predicts that
QGP photons manifest [6]. At that time, the hard photons were not measured
because the statistics ran out at the $p_{T}$ where the hard photons start
arising. Estimating the hard photon contribution had to rely on a theoretical
guidance that had large ambiguity, therefore the measurement was not able to
exhibit the thermal photon contribution [3].

The source of dileptons depend on their mass and $p_{T}$ of the measurement.
The low mass region ($m_{ee}<$1 GeV/$c^{2}$) at low $p_{T}$ is predominantly
from the in-medium decay of $\rho$ mesons and/or thermal radiation as shown in
Figure 3(a) [7].

Figure 3: (a, left) A calculation of dilepton source in low mass region. (b,
right) Outline of background components and resonances in dilepton invariant
mass spectra.

The same mass region at higher $p_{T}$ or higher mass region at lower $p_{T}$
are from mainly from the thermal radiation. The dileptons of interest are
obtained after subtracting large combinatoric background arising from Dalitz
decays of $\pi^{0}$ or $\eta$ (e.g., $\pi^{0},\eta\rightarrow
e^{+}e^{-}\gamma$). The outline of the background components and resonances in
dilepton mass spectra is shown in Figure 3(b) [8].

In this paper, we review the recent measurements of photons and dileptons from
RHIC and LHC experiments and discuss what we have learned from the data.

##  2 Hard production of photons and dileptons

One of the big success in electromagnetic radiation measurements in
relativistic heavy ion collisions is the observation of high $p_{T}$ direct
photons that are produced in initial hard scattering [9]. Figure 4(a) and (b)
show the latest direct photon $p_{T}$ spectra in Au+Au collisions at
$\sqrt{s_{NN}}$=200 GeV for various centralities, and the nuclear modification
factor ($R_{AA}$) for the 0-10 % centrality, respectively [10].

Figure 4: (a, left) Direct photon $p_{T}$ spectra in Au+Au collisions at
$\sqrt{s_{NN}}$=200 GeV measured by the PHENIX experiment. (b, right) Nuclear
modification factor ($R_{AA}$) for direct photons in 10 % central Au+Au
collisions [10].

The $R_{AA}$ is consistent with unity within quoted uncertainty, implying that
the photons from Au+Au collisions are consistent with ones expected from $p+p$
collisions. This result is not very trivial since photons from hard scattering
include both prompt and fragment components. The fragment photons may be
reduced in central Au+Au collisions due to the parton energy loss. The data is
compared with several model predictions. It is interesting to note that a
model that includes reduction of fragment photon due to energy loss of hard
scattered partons and increase of jet-photon conversion photons (coherent +
conversion +$\Delta E$ ) is not consistent with the data. The small
suppression in $R_{AA}$ seen in the highest $p_{T}$ is likely due to the fact
that the ratio of the yields in Au+Au to $p+p$ was computed without taking the
isospin dependence of direct photon production into account [11]. The LHC
heavy ion runs have the cms energy of 2.76 TeV where the hard photon
production is copious. One can make a photon isolation cut to enrich the
prompt photon component even in heavy ion collisions. Figure 5(a) shows the
isolated photon $p_{T}$ distributions in Pb+Pb collisions at
$\sqrt{s_{NN}}$=2.76 TeV measured by the ATLAS experiment at LHC, together
with ones by the CMS experiments [12, 13]. The lines show the expected values
from JetPHOX and PYTHIA event generators.

Figure 5: (a, left) Prompt photon $p_{T}$ spectra measured in Pb+Pb collisions
at $\sqrt{s_{NN}}$=2.76 TeV by the ATLAS and CMS experiments together with
JETPHOX and PYTHIA simulation [12, 13]. (b, right) Yield of $Z$-bosons scaled
by the nuclear thickness function ($T_{AB}$) as a function of centrality
measured by the CMS experiment [15].

The prompt photon component should be unaffected by the medium created in
heavy ion collisions, and was experimentally confirmed by these measurements.
Considering the fact that neither prompt photons from LHC nor prompt+fragment
photons from RHIC are suppressed, the contribution of photons from
fragmentation and jet-photon conversion could be smaller than predicted.
Several measurements benefited from the fact that these hard electromagnetic
radiations are well under control. For instance, direct photons are utilized
to quantify jet energy loss as they carry the initial momenta of jets [14,
12]. A new hard probe that became available at the LHC energy is $Z$-bosons.
The $Z$-bosons have a peak mass of $\sim$80 GeV/$c^{2}$, and are produced in
the medium with a lifetime of 0.1 fm/$c$. Therefore, they carry information on
the initial states of the collisions, and decay before they are affected by
the medium. The $Z$-boson yields scaled by nuclear thickness function
($T_{AB}$) as a function of $N_{part}$ (centrality) published by the CMS
experiments are shown in Figure 5(b) [15]. The similar results are recently
published by the ATLAS experiments [16]. The $Z$-boson yields are found to
follow the scaling of the initial hard scattering process. $Z$-boson is also
ideal for serving as a reference in many similar measurements, such as $Z$-jet
correlations [17]. We point out that $Z$-bosons can serve as one of the most
reliable tools to normalize the dilepton spectra between $p+p$ and A+A
collisions.

##  3 Thermal photon production

###  3.1 Spectra

In the conventional real photon measurement, single photons can be observed
after a huge amount of background photons coming from hadron decays
($\pi^{0}$, $\eta$, $\eta^{\prime}$ and $\omega$, etc.) are subtracted off
from the inclusive photon distributions. This fact makes it very difficult to
look at the signal at low $p_{T}$, where thermal photons from QGP manifest. A
breakthrough was made by utilizing internal conversion of photons [18].
Because these photons decay into $e^{+}e^{-}$, one can use measurement
technique for dileptons. The relation between real photon production and the
$e^{+}e^{-}$ pairs decaying from associated internal conversion photon
production can be described as follow [19]:

$\frac{d^{2}n_{ee}}{dm_{ee}}=\frac{2\alpha}{3\pi}\frac{1}{m_{ee}}\sqrt{1-\frac{4m_{e}^{2}}{m_{ee}^{2}}}\Bigl{(}1+\frac{2m_{e}^{2}}{m_{ee}^{2}}\Bigr{)}Sdn_{\gamma}$

where $\alpha$ is the fine structure constant, $m_{e}$ and $m_{ee}$ are the
masses of the electron and the $e^{+}e^{-}$ pair respectively, and $S$ is a
process dependent factor that goes to 1 when $m_{ee}\rightarrow 0$ or
$m_{ee}\ll p_{T}$. This equation also applies to the relation between the
photons from hadron decays (e.g. $\pi^{0}\rightarrow\gamma\gamma$) and the
$e^{+}e^{-}$ pairs from Dalitz decays ($\pi^{0}\rightarrow e^{+}e^{-}\gamma$).
For $\pi^{0}$ and $\eta$, the factor $S$ is given as
$S=|F(m_{ee}^{2})|^{2}(1-m_{ee}^{2}/M_{h}^{2})^{3}$, where $M_{h}$ is the
meson mass and $F(m_{ee}^{2})$ is the form factor. The analysis assumes that
the form factor for direct photons is $F(m_{ee}^{2})=1$ similar to a purely
point-like process. We can select the higher $e+e-$ invariant mass region
where $\pi^{0}$ contribution becomes off. This will eliminate the large
background coming from $\pi^{0}$ Dalitz decay. Figure 6 shows the $e^{+}e^{-}$
invariant mass distribution in minimum bias Au+Au collisions for
1.0$<p_{T}<$1.5 GeV/$c$ measured by the PHENIX experiment [18].

Figure 6: $e^{+}e^{-}$ invariant mass distribution at 1.0$<p_{T}<$1.5 GeV/$c$
together with the cocktail calculation of known sources and direct photon
internal conversion [18].

The $e^{+}e^{-}$ mass spectra were fit with the function that have terms of
the cocktail calculation of known sources ($e^{+}e^{-}$ from various hadron
Dalitz decays) and the direct photon internal conversion:

$f(m_{ee})=(1-r)f_{c}(m_{ee})+r~{}f_{\rm dir}(m_{ee})$

where $f_{c}(m_{ee})$ is the shape of the cocktail mass distribution, $f_{\rm
dir}(m_{ee})$ is the expected shape of the direct photon internal conversion.
One can obtain the signal to background ratio for a given mass window ($r$ in
the plot) at a given $p_{T}$. Using the Kroll-Wada formula [20], $r$ is
associated with the ratio at zero-mass and thus is converted to the ratio of
direct to inclusive photons:

$r=\frac{\gamma^{*}_{\rm dir}(m_{ee}>0.15)}{\gamma^{*}_{\rm
inc}(m_{ee}>0.15)}\propto\ \frac{\gamma^{*}_{\rm dir}(m_{ee}\approx
0)}{\gamma^{*}_{\rm inc}(m_{ee}\approx 0)}\ =\frac{\gamma_{\rm
dir}}{\gamma_{\rm inc}}\equiv r_{\gamma}$

Finally, the direct photon $p_{T}$ spectra is calculated as $\gamma_{\rm
inc}\times r_{\gamma}$. Figure 7(a) shows the direct photon $p_{T}$ spectra in
Au+Au and $p+p$ collisions as obtained by this procedure.

Figure 7: (a, left) PHENIX results on invariant yields of direct photons in
Au+Au collisions at various centralities together with that in $p+p$
collisions. Lines show the $N_{coll}$-scaled $p+p$ fit functions with
exponential functions [18]. (b, right) ALICE results on invariant yields of
direct photons in 2.76 TeV Pb+Pb collisions [21]

The spectra were fit with $N_{coll}$-scaled $p+p$ fit function with
exponential function, and the slope parameters of $\sim$220 MeV, which is
almost independent of centrality, were obtained. The ALICE experiment recently
measured real photons using external conversion technique, namely, by looking
at the conversion of real photons into $e^{+}e^{-}$ in the material of the
inner detectors. Figure 7(b) shows the direct photon $p_{T}$ spectrum measured
by the ALICE experiment in 0-40 % Pb+Pb collisions at $\sqrt{s_{NN}}$=2.76
TeV, together with the exponential fit to the low $p_{T}$ region [21]. The
ratio of the slope parameter from PHENIX and ALICE measurements is 1.38. From
the published result on the Bjorken energy density at Pb+Pb top centrality
from the CMS experiment [22] and one at Au+Au from the PHENIX experiment [23],
one can find that the ratio of Bjorken energy density of the LHC Pb+Pb
collisions ($\sim$14 GeV/fm3) to that of RHIC Au+Au ($\sim$5.7 GeV/fm3) is
2.6, which is smaller than the one expected the ratio of slope parameters
(1.384 = 3.65). This is because the photons measured experimentally are a sum
of the photons from all the stages from the initial to the final state of
collisions and their slope parameters reflect ”average” temperature, while the
energy density is measured at the thermalization. In order to obtain the
temperature at all the stages as well as the initial energy density, one has
to run a hydrodynamical simulation that has a realistic time profile of the
system. One might ask if the excess of photons over the initial hard
scattering process is due to cold nuclear matter effect (CNM) such as
$p_{T}$-broadening. PHENIX has recently measured the direct photon production
in $d$+Au collisions to quantify the CNM effect, using the same internal
conversion technique applied to Au+Au collisions. Figure 8 shows the $R_{AA}$
of the direct photons in $d$+Au and Au+Au collisions along with a model on the
CNM and parton energy loss effect [24]. It was found that the CNM effect to
the direct photon production is negligible compared to the large excess seen
in Au+Au collisions.

Figure 8: PHENIX results of direct photon $R_{AA}$ in $d$+Au and Au+Au
collisions at $\sqrt{s_{NN}}$=200 GeV [24].

###  3.2 Flow

After the system reached to local equilibrium, the system proceeds to
hydrodynamic expansion and its interaction with hard scattered partons, which
results in the anisotropic emission of particles. The magnitude of the
anisotropy can be studied from the azimuthal angle distribution of particles
relative to the second order event plane angle ($v_{2}$, called as elliptic
flow). For particles with low $p_{T}$ ($p_{T}<3$ GeV/$c$), the $v_{2}$ are
understood in terms of pressure-gradient anisotropy in an initial ”almond-
shaped” collision zone produced in non-central collisions. Recently, a large
$v_{2}$ of particles and its scaling in terms of kinetic energy have been
found for identified charged hadrons at RHIC [25]. It suggested that the
system is locally in equilibrium as early as 0.4 fm/$c$, and the flow occurs
at the partonic level [26]. There are predictions that photons also have a
collective motion and their $v_{2}$ show different signs and/or magnitudes
depending on the production processes [27, 28, 29]. The observable is useful
to disentangle the contributions from various photon sources in the $p_{T}$
region where they intermix. The photons from hadron-gas interaction and
thermal radiation may follow the collective expansion of a system, and give a
positive $v_{2}$. Those photons produced by jet-photon conversion or in-medium
Bremsstrahlung will increase as the size of the medium to traverse increases
and result in a negative $v_{2}$. The fragment photons will give positive
$v_{2}$ since larger energy loss of jets is expected in the orthogonal
direction to the event plane. PHENIX has measured the $v_{2}$ of direct
photons by subtracting the $v_{2}$ of hadron decay photons off from that of
the inclusive photons, following the formula below:

${v_{2}}^{dir.}=(R\times{v_{2}}^{incl.}-{v_{2}}^{bkgd.})/(R-1),\ \ \
R=(\gamma/\pi^{0})_{meas}/(\gamma/\pi^{0})_{bkgd}$

Here, $R$ is obtained either from internal conversion or external conversion
method, and $v_{2}^{incl.}$ is obtained from real photons or their external
conversions. Figure 9(a) shows the direct photon $v_{2}$ in minimum bias Au+Au
collisions at $\sqrt{s_{NN}}$=200 GeV measured by the PHENIX experiment using
both internal conversion and external conversion technique [30, 31].

Figure 9:  (a, left) Direct photon elliptic flow ($v_{2}$) in minimum bias 200
GeV Au+Au collisions measured by the PHENIX experiment, using internal
conversion (marked as arXiv:1105.4126) and external conversion technique [30,
31]. (b, right) Direct photon $v_{2}$ in 0-40 % 2.76 TeV Pb+Pb collisions
measured by the ALICE experiment using external conversion technique [32].

The $v_{2}$ of direct photons is sizable and positive, and comparable to the
flow of hadrons for $p_{T}<$3 GeV/$c$. The ALICE experiments also obtained the
direct photon $v_{2}$ in 0-40% Pb+Pb collisions, using external conversion
technique recently [32]. Although the energy density and the temperature of
the two systems are very different, the $v_{2}$ at LHC is surprisingly similar
to what PHENIX has found at RHIC.

There are many models that tried explaining the RHIC result. Several models
predicted the positive flow of the photons assuming the photons are boosted
with hydrodynamic expansion of the system, but the magnitudes from these
models are significantly lower than the measurement [33]. There are two models
that give relatively large flows.

Figure 10: Direct photon $v_{2}$ from a theoretical model increasing hadron-
gas interaction contribution [34]. (a, left) $v_{2}$ from the model
calculation, and (b, right) corresponding $p_{T}$ spectra.

Figure 10 shows a model calculation with increasing photon contribution from
hadron-gas interaction [34]. Since the hadrons have a large positive flow as
we observed, the photons produced by the interaction with these hadrons result
in a large flow. The models (a) and (b) correspond to two ways of
incorporating hard photon contribution, namely, a pQCD parametrization and a
fit to PHENIX p+p data, respectively. The spectra in low $p_{T}$ (1$<p_{T}<$3
GeV/$c$) region in this model, where QGP photons are said to dominate, is
overwhelmed by the hadron-gas interaction, and QGP contribution is hardly
seen. Figure 11 shows a model calculation for direct photon higher order flow
($v_{n}$) for two initial conditions, Glauber-based (MCGlb) and CGC-based
(MCKLN) conditions, and the corresponding $p_{T}$ spectrum for the MCGlb case
[35, 36].

Figure 11: Direct photon $v_{2}$ from a theoretical model with including
initial state fluctuation and shear viscosity [35, 36]. (a, left) $v_{2}$ from
the model calculation, and (b, right) corresponding $p_{T}$ spectra for MCGlb
case.

The shear viscosity is increased from $\eta/s$=0.08 to 0.20 when switching
from MCGlb to MCKLN, in order that the model still describes the flow of
hadrons. The dashed and sold lines are before and after viscous corrections
are applied on the rates. The reason that this model gives higher values of
$v_{n}$ for MCKLN is that the rate of QGP photons are reduced in order to
compensate the viscous entropy production. For a reference, a recent 3+1D
hydrodynamic calculation with a new CGC-inspired initial state (IP-Glasma)
gives a good description of $v_{n}$ for charged hadrons at RHIC with shear
viscosity of $\eta/s$=0.08 [37]. This implies that the determination of
$\eta/s$ is significantly affected by initial conditions. Both models
effectively assume a reduction of QGP photons and call for the hadron-gas
photon contribution. Further development at the theory side to explain the
data is clearly deserved. An another model tries to explain the large flow by
the interaction of photons with the strong magnetic field existing in the non-
central collisions [38]. The measurement of triangular flow ($v_{3}$) may be
useful to discriminate the models; for instance, the strong magnetic field
scenario gives $v_{3}\sim 0$ while hydrodynamical expansion scenario gives
sizable positive $v_{3}$ [39].

##  4 Low mass dileptons

The PHENIX experiment has performed the first measurement of the $e^{+}e^{-}$
invariant mass spectra in $p+p$ and Au+Au collisions at $\sqrt{s_{NN}}$=200
GeV at RHIC as shown in Figure 12 [19].

Figure 12: (a, left) $e^{+}e^{-}$ mass spectra in minimum bias 200 GeV Au+Au
collisions measured by the PHENIX experiments together with cocktail
calculation of known sources. (b, right) Centrality dependence of $e^{+}e^{-}$
mass spectra and corresponding cocktail calculations for the same dataset
[19]. Figure 13: (a, left) $e^{+}e^{-}$ mass spectra in various centralities
in 200 GeV Au+Au collisions measured by the STAR experiment. (b, right) Ratios
of data to cocktail calculation for each centralities [40].

Figure 14:  (a, left) $e^{+}e^{-}$ mass spectra in 20-40 % Au+Au collisions
measured by the PHENIX experiment in 2010, with a hadron blind detector (HBD)
installed. (b, right) Comparison of the ratio of the integrated yields in the
LMR to the cocktail calculations in 2004 and 2010 dataset [31].

The left panel shows the minimum bias mass spectra with hadronic cocktail
components. The right panel shows the centrality dependence of mass spectra
with cocktail calculations. The $p+p$ results are well reproduced by the
cocktail of known sources of $e^{+}e^{-}$, whereas the Au+Au data show a
strong enhancement at low mass region (LMR, $0.15\\!<M_{ee}\\!<\\!0.75$
GeV/$c^{2}$) compared to cocktail calculations. The enhancement for the 0-10 %
central collisions is a factor of $7.6\pm 0.5^{stat}\pm 1.3^{syst}\pm
1.5^{model}$ and that for minimum bias collisions is $4.7\pm 0.4^{stat}\pm
1.5^{syst}\pm 0.9^{model}$, respectively. The STAR experiment also recently
obtained $e^{+}e^{-}$ mass spectra in $p+p$ and Au+Au collisions at
$\sqrt{s_{NN}}=200$ GeV as shown in Figure 13 [40, 41]. In the same LMR, the
enhancement for 0-10 % central collisions is $1.72\pm 0.10\pm 0.50$, and that
for minimum bias collisions is $1.53\pm 0.07\pm 0.41$, respectively. Clearly,
there is a discrepancy in the magnitudes of the excess by a factor of 3-4
between two experiments. PHENIX has installed a hadron blind detector (HBD) in
2010 run in order to reduce systematic uncertainty of the measurement and also
to confirm the previous result. The HBD is a Cherenkov detector with CF4 gas
and rejects the $e^{+}e^{-}$ tracks from photon conversions and Dalitz decay
of $\pi^{0}$’s, which are major background in LMR measurement, by looking at
the opening angle of pair tracks in a magnetic field free region [42]. Figure
14(a) shows the $e^{+}e^{-}$ mass spectra measured in 20-40 % Au+Au collisions
in the 2010 run with the HBD [31]. The analysis for the most central
collisions are still in progress. PHENIX measurements of 2004 and 2010 have
several differences in both data (magnetic field and detector material
budget), and cocktail calculations (MC@NLO is used for charm contribution
estimate in 2010). PHENIX has made a comparison of the data/cocktail ratio for
the integrated yield in LMR obtained in the 2004 and 2010 runs for the three
centrality bins as shown in Figure 14(b). It is seen that the two runs give
consistent results within uncertainties. It should be noted that the previous
PHENIX data showed most of the excess is seen in most central events (0-20 %),
whose analysis is still in progress. The discrepancy between PHENIX and STAR
results persists until PHENIX comes up with the new result with HBD. As for
the STAR data, some concern is put on the fact that the cocktail calculation
in $p+p$ collisions over-predicts data points in the LMR, though they are
still consistent within systematic errors. If one correct for the over-
prediction, it is possible STAR and PHENIX see the same amount of excess.
Further effort to understand the discrepancy is deserved.

As measured for photons, the measurement of dilepton flow is useful. There is
a radial flow measurement of dileptons by the NA60 experiment [43], and an
elliptic flow measurement by the STAR experiment [44, 45] that could help
disentangling the source of the dileptons that contribute to this particular
mass region. A theory study is also in progress [46].

###  4.1 Energy dependence of LMR dilepton production

Figure 15: $e^{+}e^{-}$ mass spectra measured by the STAR experiment in 19.6,
62.4 and 200 GeV Au+Au collisions, together with cocktail calculation and a
theoretical model assuming in-medium modification of $\rho$ spectral function
[40].

The excess in LMR has been explained by various models including in-medium
broadening of the $\rho$ mesons or their mass shift [47]. The measurement of
the energy dependence of the excess may provide a discrimination of these
models. The STAR experiment recently came up with the $e^{+}e^{-}$ invariant
mass spectra for minimum bias Au+Au collisions at $\sqrt{s_{NN}}=$ 19.6, 62.4,
and 200 GeV as shown in Figure 15 [40]. STAR observed a qualitatively similar
excess as observed by the CERES measurements in the Pb$+$Au at
$\sqrt{s_{\mathrm{NN}}}=$17.2 GeV [48]. In each of the three panels the hadron
cocktail simulation includes contributions from Dalitz decays, photon
conversions (19.6 GeV only), and the dielectron decay of the $\omega$ and
$\phi$ vector mesons. The cocktail simulations purposely exclude contributions
from $\rho$ mesons. Instead, these are explicitly included in the model
calculation by Rapp [47] which involve in-medium modifications of the $\rho$
meson spectral shape in the isentropic fireball evolution. The LMR enhancement
measured by STAR are consistently agreeing with these model calculations
within the quoted errors. One should be careful on taking the absolute
magnitude of the excess in all energies, provided that there is still the
issue of inconsistency between STAR and PHENIX. Nonetheless, it is interesting
to see that the relative change as a function cms energy is well described by
this model calculation. If it is the case, the flow of this mass region may
exhibit the $KE_{T}$ scaling with assuming they are all $\rho$.

##  5 Summary

The recent results on direct photons and dileptons in high energy heavy ion
collisions, obtained particularly at RHIC and LHC are reviewed. The results
are new not only in terms of the probes, but also in terms of the precision.
Much progress has been made in understanding high $p_{T}$ direct photons as
well as $Z$-bosons with the latest RHIC and LHC results. The soft single
photons have been studied down to lower $p_{T}$ using internal conversion and
external conversion technique at RHIC and LHC, and exhibit the average
temperature of the system. A large flow of soft photons was also observed both
at RHIC and LHC, which are not explained by models considering QGP only so
far. The low mass dilepton excess was observed at the PHENIX and STAR
experiments at RHIC in Au+Au collisions, but are not quantitatively agreeing
each other. The flow of dilepton is as important as the one of photons, and
experiments should make an effort to improve the measurement. The STAR results
on dileptons at several cms energies shows that the excess in LMR is
consistent with in-medium modification of the $\rho$ meson spectral function.
If it is the case, the flow of this mass region may exhibit the $KE_{T}$
scaling with assuming they are all $\rho$. This is an interesting topic to
explore.

## References

  * [1] See for example, J-e. Alam, S. Raha, B. Sinha, Phys. Rep. 273 (1996) 243.
  * [2] L. E. Gordon and W. Vogelsang, Phys. Rev. D 48, 3136 (1993).
  * [3] S. Turbide, R. Rapp and C. Gale, Phys. Rev. C 69, 014903 (2004).
  * [4] R. J. Fries, B. Muller and D. K. Srivastava, Phys. Rev. C 72, 041902 (2005).
  * [5] A. Adare et al. [PHENIX Collaboration], Phys. Rev. C 87, 034911 (2013).
  * [6] M. M. Aggarwal et al. [WA98 Collaboration], Phys. Rev. Lett. 85, 3595 (2000).
  * [7] R. Rapp, Phys. Rev. C 63, 054907 (2001) [Phys. Rev. A 83, 043833 (2011)]
  * [8] A. Drees, Nucl. Phys. A 830, 435C (2009).
  * [9] S. S. Adler et al. [PHENIX Collaboration], Phys. Rev. Lett. 94, 232301 (2005).
  * [10] S. Afanasiev et al. [PHENIX Collaboration], Phys. Rev. Lett. 109, 152302 (2012).
  * [11] F. Arleo, JHEP 0609, 015 (2006).
  * [12] P. Steinberg [ATLAS Collaboration], Nucl. Phys. A 904-905, 233c (2013).
  * [13] S. Chatrchyan et al. [CMS Collaboration], Phys. Lett. B 710, 256 (2012).
  * [14] A. Adare et al. [PHENIX Collaboration], Phys. Rev. Lett. 111, 032301 (2013).
  * [15] S. Chatrchyan et al. [CMS Collaboration], Phys. Rev. Lett. 106, 212301 (2011).
  * [16] G. Aad et al. [ATLAS Collaboration], Phys. Rev. Lett. 110, 022301 (2013).
  * [17] S. Chatrchyan et al. [CMS Collaboration], JHEP 1305, 090 (2013).
  * [18] A. Adare et al. [PHENIX Collaboration], Phys. Rev. Lett. 104, 132301 (2010).
  * [19] A. Adare et al. [PHENIX Collaboration], Phys. Rev. C 81, 034911 (2010).
  * [20] N. M. Kroll and W. Wada, Phys. Rev. 98, 1355 (1955).
  * [21] M. Wilde [ALICE Collaboration], Nucl. Phys. A 904-905, 573c (2013).
  * [22] S. Chatrchyan et al. [CMS Collaboration], Phys. Rev. Lett. 109, 152303 (2012).
  * [23] S. S. Adler et al. [PHENIX Collaboration], Phys. Rev. C 71, 034908 (2005) [Erratum-ibid. C 71, 049901 (2005)].
  * [24] A. Adare, et al. [PHENIX Collaboration], Phys. Rev. C 87, 054907 (2013).
  * [25] See for example, A. Adare et al. [PHENIX Collaboration], Phys. Rev. Lett. 98, 162301 (2007).
  * [26] B. Schenke, S. Jeon and C. Gale, Phys. Rev. C 82, 014903 (2010).
  * [27] S. Turbide, C. Gale and R. J. Fries, Phys. Rev. Lett. 96, 032303 (2006).
  * [28] R. Chatterjee, E. S. Frodermann, U. W. Heinz and D. K. Srivastava, Phys. Rev. Lett. 96, 202302 (2006).
  * [29] S. Turbide, C. Gale, E. Frodermann and U. Heinz, Phys. Rev. C 77, 024909 (2008).
  * [30] A. Adare et al. [PHENIX Collaboration], Phys. Rev. Lett. 109, 122302 (2012).
  * [31] I. Tserruya [PHENIX Collaboration], Nucl. Phys. A 904-905, 225c (2013).
  * [32] D. Lohner [ALICE Collaboration], J. Phys. Conf. Ser. 446, 012028 (2013).
  * [33] R. Chatterjee and D. K. Srivastava, Nucl. Phys. A 830, 503C (2009).
  * [34] H. van Hees, C. Gale and R. Rapp, Phys. Rev. C 84, 054906 (2011).
  * [35] C. Shen, U. W. Heinz, J. -F. Paquet, I. Kozlov and C. Gale, arXiv:1308.2111 [nucl-th].
  * [36] C. Shen, U. WHeinz, J. -F. Paquet and C. Gale, arXiv:1308.2440 [nucl-th].
  * [37] C. Gale, S. Jeon, B. Schenke, P. Tribedy and R. Venugopalan, Phys. Rev. Lett. 110, 012302 (2013).
  * [38] See for example, A. Bzdak and V. Skokov, Phys. Rev. Lett. 110, 192301 (2013).
  * [39] See for example, R. Chatterjee, D. K. Srivastava and T. Renk, arXiv:1401.7464 [hep-ph].
  * [40] F. Geurts [STAR Collaboration], Nucl. Phys. A904-905 2013, 217c (2013).
  * [41] L. Adamczyk et al. [STAR Collaboration], arXiv:1312.7397 [hep-ex].
  * [42] W. Anderson, B. Azmoun, A. Cherlin, C. Y. Chi, Z. Citron, M. Connors, A. Dubey and J. M. Durham et al., Nucl. Instrum. Meth. A 646, 35 (2011).
  * [43] R. Arnaldi et al. [NA60 Collaboration], Phys. Rev. Lett. 100, 022302 (2008).
  * [44] X. Cui [STAR Collaboration], J. Phys. Conf. Ser. 446, 012047 (2013).
  * [45] L. Adamczyk et al. [STAR Collaboration], arXiv:1402.1791 [nucl-ex].
  * [46] G. Vujanovic, C. Young, B. Schenke, R. Rapp, S. Jeon and C. Gale, arXiv:1312.0676 [nucl-th].
  * [47] See for example, R. Rapp and J. Wambach, Eur. Phys. J. A 6, 415 (1999).
  * [48] G. Agakichiev et al. [CERES Collaboration], Eur. Phys. J. C 41, 475 (2005).


11footnotetext: Corresponding author

# Efficient hybrid technique for generating sub-grid haloes in reionization
simulations

Ankur Barsode  and Tirthankar Roy Choudhury

###### Abstract

Simulating the distribution of cosmological neutral hydrogen (HI) during the
epoch of reionization requires a high dynamic range and is hence
computationally expensive. The size of the simulation is dictated by the
largest scales one aims to probe, while the resolution is determined by the
smallest dark matter haloes capable of hosting the first stars. We present a
hybrid approach where the density and tidal fields of a large-volume, low-
resolution simulation are combined with small haloes from a small-volume,
high-resolution box. By merging these two boxes of relatively lower dynamic
range, we achieve an effective high-dynamic range simulation using only 13% of
the computational resources required for a full high-dynamic range simulation.
Our method accurately reproduces the one- and two-point statistics of the halo
field, its cross-correlation with the dark matter density field, and the two-
point statistics of the HI field computed using a semi-numerical code, all
within 10% accuracy at large scales and across different redshifts. Our
technique, combined with semi-numerical models of reionization, provides a
resource-efficient tool for modeling the HI distribution at high redshifts.

##  1 Introduction

Numerical simulations of reionization, particularly those designed for
inference with upcoming 21 cm surveys, require a high dynamic range (for
recent reviews, see [1, 2]). At one end, simulation volumes need to be
$\gtrsim 200$ cMpc to model the largest scales probed by these surveys [3, 4,
5, 6]. At the other end, it is necessary to resolve relatively small mass
haloes, given that the first galaxies formed in atomically cooled haloes with
masses $\sim 10^{8}\,\mathrm{M}_{\odot}$. Typically, these simulations would
require $\gtrsim 10^{11}$ matter particles which, in turn, would demand
$\gtrsim 10$ TB of memory. The requirements become even more demanding if one
aims to resolve haloes where molecular cooling is efficient. A high number of
particles is also necessary for resolving small spatial scales, which is
crucial for accurately modeling the sinks of ionizing photons.

In semi-numerical simulations of reionization, where the radiative transfer
physics is replaced by an appropriate photon counting algorithm, a common
compromise is to use large-volume simulation boxes with relatively coarse
resolution [7, 8, 9, 10, 11, 12, 13, 14]. In some extensions to these models,
the effect of recombinations is accounted for through sub-grid models [15,
16]. The coarse resolution of these simulations also implies that they are not
able to resolve the small mass haloes. Instead, they use the gridded density
fields to compute the conditional halo mass functions within each grid cell to
populate them with sub-grid haloes. This approach is typically sufficient
because the primary interest lies in modeling the large-scale 21 cm signal
pertinent to upcoming experiments, rather than in the detailed properties of
the haloes (e.g., their density profiles). Therefore, only the location and
mass of the haloes are required for these simulations.

Such sub-grid prescriptions for populating small mass haloes in coarse
resolution simulations have been widely used for modeling the epoch of
reionization. One of the early methods was based on using conditional Press-
Schechter mass function [17] scaled to match the Sheth-Tormen mass function
[18], which was found to provide a good match to the halo bias from $N$-body
simulations [19]. This method has been widely used in semi-numerical models of
reionization for populating haloes that are otherwise unresolved in the box
[20, 21, 8, 22, 23]. A slightly alternate approach has been to use the
analytical conditional mass function based on the ellipsoidal collapse to
populate the unresolved haloes [14, 24]. Alternatively, a fit to the
conditional mass function, obtained from a high-resolution simulation capable
of resolving the smallest haloes forming stars [25, 26] has also been used.
Obvious extensions to these techniques involve introducing scatter in the halo
numbers [8, 24].

The conditional halo mass functions based on analytical calculations or
functional fits encounter several challenges. For instance, the analytical
relations require knowledge of the Lagrangian density as opposed to the
Eulerian ones obtained from the simulations. Usually, the Eulerian to
Lagrangian mapping is done using spherical approximation [27], which may not
be accurate when the grid cell sizes become small. To bypass this reliance on
analytical relations, an alternate approach could be to combine information
from two independent simulation boxes, one high-resolution but small volume,
and another coarse resolution with large volume. The motivation behind such
hybrid approaches comes from the fact that computational cost increases with
the dynamic range, so two simulations with narrow dynamic ranges can be quite
inexpensive as compared to simulating with the full dynamic range. Such
methods have been used to model the density and velocity fields during the
late stages of reionization [28, 29], where they combined a large-scale dark-
matter-only simulation with a high-resolution hydrodynamical simulation to
obtain a large-scale, high-resolution HI density field.

In this work, we refine the hybrid method based on multiple simulation volumes
to populate low-mass haloes in the coarse resolution box that are otherwise
unresolved. Our method does not require any analytical fits to the conditional
mass functions or any mapping to the Lagrangian field, rather depending on
matching similar grid cells from the two boxes. We also improve on the
existing methods in the sense that we use not only the density to condition
the mass function but also the tidal field. As it turns out, using the tidal
information leads to more robust conditional mass functions in agreement with
the actual one. We restrict our analysis to high redshifts $z\gtrsim 6$ and
evaluate the success of our method with respect to how accurately it can
produce the large-scale ionized hydrogen field.

The outline of the paper is as follows: The following section illustrates our
methodology for constructing the hybrid box containing the sub-grid haloes. We
demonstrate the accuracy of our method by comparing it with a full high-
dynamic range simulation in section §3. Finally, we conclude in §4 with a
brief discussion of the advancements our method can bring to the studies of
reionization.

##  2 Methodology

Figure 1:  Illustration of the methodology to construct the hybrid box: (1)
high-mass haloes from the large low-resolution box (LB) are directly copied to
the hybrid box (HB), (2) for each grid cell in LB, best matching tiles are
found in a small high-resolution box (SB) using either the density field or
the eigenvalues of the tidal tensor field, (3) the low-mass haloes and dark
matter density field smoothed at the grid scale are picked up from these
locations and substituted in the hybrid box, ensuring that the large scale
density field in HB is identical to LB.

In this section, we describe our algorithm for constructing the hybrid
simulation boxes in detail.

Let us assume that we want to construct a hybrid box equivalent to a (cubical)
simulation box of volume $V_{\mathrm{box}}=L_{\mathrm{box}}^{3}$ with
$N_{\mathrm{part}}$ particles. The smallest halo that can be resolved in this
simulation would depend on the particle mass resolution
$M_{\mathrm{part}}=(V_{\mathrm{box}}/N_{\mathrm{part}})~{}\bar{\rho}_{m}$,
where $\bar{\rho}_{m}$ is the mean density of matter. Our procedure for
constructing the hybrid box requires two simulations of relatively coarser
resolution, one having a size as large as the target simulation and another
having a mass resolution the same as the target one. The details are as
described in the following:

  * •

Large Box (LB): This simulation box has a volume
$V_{\mathrm{box,large}}=V_{\mathrm{box}}$ corresponding to the ultimate volume
of the simulation we want to construct. Let $N_{\mathrm{part,large}}$ be the
number of particles in this simulation. Clearly, if we want this simulation to
be less expensive than the target simulation, it must contain fewer particles
than the target one, i.e., $N_{\mathrm{part,large}}<N_{\mathrm{part}}$.

  * •

Small Box (SB): This simulation box has a volume $V_{\mathrm{box,small}}$ with
$N_{\mathrm{part,small}}$ number of particles, chosen such that the mass
resolution is equal to the target simulation, i.e.,

$\dfrac{V_{\mathrm{box,small}}}{N_{\mathrm{part,small}}}=\dfrac{V_{\mathrm{box}}}{N_{\mathrm{part}}}.$
(2.1)

As with the LB, we want $N_{\mathrm{part,small}}<N_{\mathrm{part}}$, hence
$V_{\mathrm{box,small}}<V_{\mathrm{box}}$. This also implies that the mass
resolution for SB is finer than the LB, i.e.,
$M_{\mathrm{part,small}}<M_{\mathrm{part,large}}$.

The simulation outputs are used to compute the cloud-in-cell (CIC)-smoothed
density contrast
$\delta_{m}(\mathbf{x})\equiv\rho_{m}(\mathbf{x})/\bar{\rho}_{m}$ on a uniform
grid. The grid size $\Delta x_{\mathrm{CIC}}$ for each box is chosen to be
equal to the mean inter-particle distance in that box.

However, in order to make the semi-numerical algorithms more efficient, it is
often the case that the ionization fields are generated at a relatively
coarser resolution. Denoting this target grid size as $\Delta
x_{\mathrm{grid}}$ ($\geq\Delta x_{\mathrm{CIC}}$ for either box), we smooth
the high-resolution density field to this resolution using the boxcar filter.
The number of grid cells for the two simulation boxes are given by

$N_{\mathrm{grid,small}}=\frac{V_{\mathrm{box,small}}}{\Delta
x_{\mathrm{grid}}^{3}}<\frac{V_{\mathrm{box,large}}}{\Delta
x_{\mathrm{grid}}^{3}}=N_{\mathrm{grid,large}}.$ (2.2)

This smoothed density field is used for computing the tidal tensor field

$\nabla^{2}T_{ij}(\mathbf{x})=\frac{\partial^{2}\delta_{m}(\mathbf{x})}{\partial
x_{i}\partial x_{j}},~{}~{}~{}i,j=1,2,3.$ (2.3)

In practice, we compute the derivative with respect to the spatial coordinates
$x_{i}$ in Fourier space and then transform it back to the configuration
space. The main quantities of our interest are the three eigenvalues
$\lambda_{1}(\mathbf{x}),\lambda_{2}(\mathbf{x}),\lambda_{3}(\mathbf{x})$ of
$T_{ij}(\mathbf{x})$, ordered such that
$\lambda_{1}(\mathbf{x})\leq\lambda_{2}(\mathbf{x})\leq\lambda_{3}(\mathbf{x})$.
Note that the density contrast
$\delta_{m}(\mathbf{x})=\lambda_{1}(\mathbf{x})+\lambda_{2}(\mathbf{x})+\lambda_{3}(\mathbf{x})$
is the sum of these eigenvalues of the tidal tensor.

We further identify and generate the halo catalogue for each of the two boxes,
using the Friend-of-Friends (FoF) algorithm [30]. The smallest halo
$M_{h,\mathrm{min}}$ used in this work is assumed to have 20 particles. Note
that this minimum halo mass would be smaller in SB than in LB
$M_{h,\mathrm{min,small}}<M_{h,\mathrm{min,large}}$, since SB has a resolution
finer than LB. Also, since the resolution for SB is chosen to be identical to
that of the target box, we have $M_{h,\mathrm{min,small}}=M_{h,\mathrm{min}}$.
Once identified, we treat the haloes as point masses having a mass equal to
that of the halo, and located at the center of mass of its constituent
particles.

To compute the power spectrum of the halo field, its cross-spectrum with the
density field, and ultimately the ionization field (to be discussed later), we
also need to transform these point mass haloes into a gridded halo field
$\delta_{h}(\mathbf{x})$. This is done by first CIC-depositing the haloes for
the two boxes at the fine resolution of $\Delta x_{\mathrm{CIC}}$, and then
boxcar-smoothing them to the target resolution $\Delta x_{\mathrm{grid}}$.

We also run another simulation, called Reference Box (RB), which is used for
comparing the results obtained from our hybrid algorithm (discussed in the
next section §3). The specifications of this simulation are identical to the
ones we wish to achieve using the hybrid algorithm. Specifically, it has a
volume $V_{\mathrm{box}}=L_{\mathrm{box}}^{3}$ with $N_{\mathrm{part}}$
particles. The particle positions for this box are CIC-smoothed to obtain the
density field, which is then boxcar-smoothed to a coarse resolution $\Delta
x_{\mathrm{grid}}$ for later use. The haloes are identified using the same FoF
algorithm as above and the halo field at the desired resolution is computed in
a manner identical to those for the LB and SB. Note that this box does _not_
play any role in the implementation of the hybrid algorithm, it is only used
for validating the results obtained from the hybrid algorithm. We also control
the seed for the initial realization of the cosmological fields for the RB so
that the large-scale fields are identical to the LB. This ensures that the
effects of cosmic variance are minimized at scales comparable to
$L_{\mathrm{box}}$.

In this work, we choose to work with the case where,

$N_{\mathrm{part,small}}=N_{\mathrm{part,large}}=\frac{N_{\mathrm{part}}}{8},$
(2.4)

which implies that both LB and SB contain 8 times fewer particles than the
target simulation. This also implies the following relations between the box
lengths

$L_{\mathrm{box,small}}=\frac{L_{\mathrm{box,large}}}{2}=\frac{L_{\mathrm{box}}}{2},$
(2.5)

the number of grid cells along each direction

$N_{\mathrm{grid,small}}=\frac{N_{\mathrm{grid,large}}}{8},$ (2.6)

and the minimum halo mass

$M_{h,\mathrm{min,small}}=\frac{M_{h,\mathrm{min,large}}}{8}=M_{h,\mathrm{min}}.$
(2.7)

For the analyses in this paper, dark-matter-only simulations were run using
GADGET-2222https://wwwmpa.mpa-garching.mpg.de/gadget/ [31]. A flat
$\Lambda$CDM cosmology was assumed with
$h=0.678,\Omega_{m}=0.308,\Omega_{b}=0.0,\sigma_{8}=0.829,n_{s}=0.961$. The
initial conditions were generated at $z=60$ using N-GenIC [32]. The initial
power spectrum was generated using the analytical fit of Eisenstein & Hu [33].
The details of the simulations used in this work are mentioned in Table 1. Our
RB contains $1024^{3}$ particles, the maximum we can use given the memory
available to us. We choose the box size $80h^{-1}\mathrm{cMpc}$ which is
slightly smaller than what is ideal for modeling the 21 cm signal, and it
allows us to probe halo masses an order of magnitude larger than the
atomically cooled haloes. The SB and LB each contain $512^{3}$ particles,
hence requiring almost an order of magnitude less memory than RB. The total
combined run time for SB and LB is $\sim 7$ times less than RB.

Name | $L_{\mathrm{box}}$ | $N_{\mathrm{part}}$ | $M_{\mathrm{min}}$ | CPU timea | Max RAM
---|---|---|---|---|---
| ($h^{-1}\mathrm{cMpc}$) |  | ($h^{-1}\mathrm{M}_{\odot}$) | (hours) | (GB)
Large Box (LB) | 80 | $512^{3}$ | $6.52\times 10^{9}$ | $\approx 2.1\times 10^{2}$ | $\approx 20$
Small Box (SB) | 40 | $512^{3}$ | $8.15\times 10^{8}$ | $\approx 2.2\times 10^{2}$ | $\approx 20$
Reference Box (RB) | 80 | $1024^{3}$ | $8.15\times 10^{8}$ | $\approx 2.9\times 10^{3}$ | $\approx 160$

asimulations were run only up to $z=5$.

Table 1:

###  2.1 Constructing the hybrid box

The LB contains the large-scale density field on a uniform grid of size
$\Delta x_{\mathrm{grid}}$, which is to be used for generating the ionization
field. However, this simulation can resolve haloes only above a mass
$M_{h,\mathrm{min,large}}=6.52\times 10^{9}h^{-1}\mathrm{M}_{\odot}$, while
our aim is to resolve haloes as small as $M_{h,\mathrm{min}}=8.15\times
10^{8}h^{-1}\mathrm{M}_{\odot}$, as obtained in the RB. It is possible to
resolve these haloes in the SB which, however, is unable to access the largest
scales we are interested in. To obtain the conditional halo mass function for
$M_{h}\geq M_{h,\mathrm{min}}$ at every grid cell in LB, we use the grid cells
from the SB. This implicitly assumes that the conditional mass function is
determined primarily by local conditions. Under such an assumption, for a
given grid cell in LB, one only needs to find cells in the SB that have
similar properties.

Now, regarding which properties could be appropriate, we consider two options:

  1. 1.

Density: The first option is that we find a grid cell in the SB which have a
matter density contrast $\delta_{m}$ similar to the cell in LB under
consideration. Once such a match is identified, we simply populate all the
haloes with masses between $M_{\mathrm{min,small}}$ and
$M_{\mathrm{min,large}}$ in the cell in LB. This avoids any double-counting,
as all haloes with $M_{h}\geq M_{\mathrm{min,large}}$ are already present in
the LB. It is likely that the match between the $\delta_{m}$’s from the LB and
SB would not be perfect, so to conserve the mass within haloes, we scale the
masses of all these sub-grid haloes by the ratio of the densities. This leads
to a hybrid box with a population of haloes with $M_{h}\geq
M_{h,\mathrm{min,small}}$.

  2. 2.

Tidal tensor: We also explore another possibility that the conditional halo
properties may not be determined solely by the local density, but also by the
tidal environment. To account for this possibility, we use the eigenvalues
$\lambda_{i}$ of the tidal tensor to find the matching grid cells between LB
and SB. While trying to identify the cell with the best match, we simply find
one having the nearest value of
$\sum_{i=1}^{3}(\lambda_{i}^{\mathrm{LB}}-\lambda_{i}^{\mathrm{SB}})^{2}$.
Following this identification of the best-matched cell, we populate the sub-
grid haloes in a manner identical to the density-based method described above.
In particular, in this case, too, we scale the halo masses by the ratio of the
cell densities.

A schematic representation of the above algorithm can be found in Figure 1.

The above method provides sub-grid haloes in the density field of the LB,
which forms our Hybrid Box (HB). Note that the above methods would lead to
discontinuities in the halo field at the boundaries of the grid cells, hence
the correlation functions would show a lack of correlation at scales $k\gtrsim
2\pi/\Delta x_{\mathrm{grid}}$. However, both the above methods are expected
to preserve the density-halo correlation at larger scales $k\lesssim
2\pi/\Delta x_{\mathrm{grid}}$. We will discuss these in more detail in the
next section §3 by comparing the results from the HB with the RB.

###  2.2 Ionization field

Our main quantity of interest is the ionization field. We use the semi-
numerical model SCRIPT333https://bitbucket.org/rctirthankar/script [14, 34] to
generate the ionization fields for both HB and RB. The algorithm requires a
density field in a rectangular grid, along with haloes that are capable of
producing ionizing photons. We assume the ionizing photons produced in a halo
to be proportional to the halo mass $M_{h}$, with the constant of
proportionality being denoted by the ionization efficiency $\zeta$. We also
assume that all haloes with $M_{h}\geq M_{h,\mathrm{min}}$ are capable of
forming stars and contributing to the ionizing photon budget.

The algorithm of SCRIPT conserves the number of photons explicitly and hence
provides large-scale ionization maps that are independent of the grid size
$\Delta x_{\mathrm{grid}}$ used for generating the ionization fields. Given
the density field on a grid and the number of ionizing photons produced by
haloes in a grid cell, the algorithm assigns ionized regions of appropriate
volumes around the ‘source’ cells. This process is repeated independently for
all cells containing sources in the box, and it may result in some grid cells
that receive photons from multiple source cells being ‘overionized’. These
excess ionizing photons in these unphysical overionized cells are subsequently
distributed among the surrounding neighboring cells which are yet to be fully
ionized. The process is continued until all the overionized cells are properly
accounted for.

##  3 Results and Discussion

In this section, we compare our hybrid method for generating sub-grid haloes
with the reference box. We compare our results for three choices of the grid
size, namely, $\Delta x_{\mathrm{grid}}=1.25$, $0.625$, and
$0.3125~{}h^{-1}\mathrm{cMpc}$, which correspond to
$N_{\mathrm{grid}}=64^{3}$, $128^{3}$ and $256^{3}$, respectively. Note that
smaller grid sizes will allow us to probe the power spectra at higher
$k$-values. However, such grids would contain relatively fewer haloes, and
hence the variance would be larger. Consequently, this may lead to more
discrepancies between the distribution of haloes in the cells in HB compared
to the RB. We choose the reference redshift as $z=7$ for most of the results,
but we also show the evolutionary trends of our scheme in the latter part.

###  3.1 Matter fields

Let us first compare the gridded density fields of HB and RB. _Note that the
density field for the HB is simply that of LB_ , and the details of the
hybridization scheme do not play any role here. The power spectra $P_{m}(k)$
of the density fields for the two boxes are shown in Figure 2, and their
relative discrepancy is shown in the lower panel of the same figure. The
lowest $k$ we consider for the plot is $\sim 2\pi/(40h^{-1}\mathrm{cMpc})$ as
smaller $k$-modes show effects of aliasing arising from the box length of
$80h^{-1}\mathrm{cMpc}$. The maximum $k$ used depends on the grid size $\Delta
x_{\mathrm{grid}}$.

It is clear that the power spectra are identical at large scales $k\lesssim
1h/$ cMpc. This is because we have controlled the seeds for generating the
initial density field of the RB and LB such that they are identical at large
scales. At smaller scales, the power of LB is lower than the RB, a direct
consequence of the lower force resolution of LB. At $k\gtrsim 4h/$ cMpc, the
difference exceeds $\sim 5\%$.

Figure 2: Comparison of hybrid dark matter density power spectra at $z=7$ with
the reference simulation (thin black lines) for different grid sizes as
mentioned in the legend. Solid lines (scattered stars) show the spectra of
boxes generated by matching tidal-tensor eigenvalues (density) between grid
cells of LB and SB. Also shown is the power spectrum of the LB (green
crosses). The lower panel shows the % error in each power spectrum with
respect to the reference box. Note that our hybridization scheme ensures that
the density field of the hybrid box (HB) is identical to the large low-
resolution box (LB). The suppression in the power at large-$k$ is due to lower
force resolution in the LB.

###  3.2 Hybrid halo catalogue

Figure 3: Comparison of hybrid halo mass functions at $z=7$ with the reference
simulation. The red dotted line with shaded errors show the results from the
reference box (RB), while the solid lines are for the hybrid boxes built using
the tidal tensor matching criteria, and scattered stars are the hybrid boxes
built using the density matching criteria. There is little difference between
the latter two. The mass function corresponding to LB and SB are shown in
green crosses and teal circles respectively. The lower panel shows the % error
in each mass function with respect to the reference box. The error bars on
each mass function are calculated assuming the number of haloes in each mass
bin follows a Poisson distribution.

The most important outcome of our method is the high dynamic range hybrid halo
catalogue. Figure 3 shows the mass function of the haloes in the HB, compared
against that from RB, at a redshift of 7. We also show the mass functions
obtained from LB and SB. The RB covers a wide range of haloes in the range
$10^{9}h^{-1}\mathrm{M}_{\odot}\lesssim M_{h}\lesssim 5\times
10^{11}h^{-1}\mathrm{M}_{\odot}$. In contrast, the LB haloes are limited to
$M_{h}\gtrsim 10^{10}h^{-1}\mathrm{M}_{\odot}$, while the number of high mass
haloes ($M_{h}\gtrsim 10^{11}h^{-1}\mathrm{M}_{\odot}$) in the SB become
rarer, as is clear from the size of the error-bars. The HB, which is
constructed using both LB and SB, provides a mass function that is very
similar to the RB. At low masses $M_{h}\lesssim
10^{11}h^{-1}\mathrm{M}_{\odot}$, the agreement is within $10\%$. At large
masses, the match is less accurate due to Poisson fluctuations caused by the
limited number of haloes. It is important to note that this discrepancy is
solely due to the differences in the LB and RB and does not directly affect
our hybridization scheme. The agreement between the halo mass functions of HB
and RB has little to no dependence on the matching criteria of the cells
(density/tidal tensor) or the grid sizes.

Figure 4: Joint probability distribution of the collapsed fraction
$f_{\mathrm{coll}}$ and density $1+\delta_{m}$ for hybrid boxes (solid/dotted
lines corresponding to tidal-tensor/density matching criteria) as well as
reference boxes smoothed at the same resolution as the hybrid box (light
shaded regions surrounded by dotted lines). The central panel shows the joint
probability distribution, while the upper and side panels show a marginalized
probability distribution of $\log_{10}{f_{\mathrm{coll}}}$ and $1+\delta_{m}$
respectively.

Given that the HB produces a good match to the total halo mass function, we
next check if the haloes are properly correlated with the underlying density
field. We compute the joint distribution of the fraction of mass collapsed to
form haloes
$f_{\mathrm{coll}}=\rho_{m}^{-1}\int_{M_{\mathrm{min}}}^{\infty}{\text{d}}M_{h}~{}M_{h}~{}{\text{d}}n/{\text{d}}M_{h}$,
and the density $\delta_{m}$. The results are shown in Figure 4. As is clear,
the HB captures qualitative features of the correlation between $\delta_{m}$
and $f_{\mathrm{coll}}$ as found in RB. Furthermore, the quantitative
agreement between the two for grid sizes $\Delta x_{\mathrm{grid}}\geq
0.625h^{-1}\mathrm{cMpc}$ too is excellent. There are some minor disagreements
between the HB and RB at low $\delta_{m}$ and low $f_{\mathrm{coll}}$, which
could probably arise from fewer number of haloes in low-density cells that are
not well captured in the SB. However, these cells contribute negligibly to the
ionization photon budget and hence are expected to play a relatively minor
role in the ionization field. The disagreement in the joint distribution for
$\Delta x_{\mathrm{grid}}=0.3125h^{-1}\mathrm{cMpc}$ indicates that the method
tends to perform less reliably for smaller grid sizes, possibly due to the
large variance in halo numbers. We also find that the matching criteria of the
grid cells (tidal tensor/density) do not affect the joint distribution in any
significant way.

Figure 5: Power spectrum $P_{hh}(k)$ of the halo overdensity field scaled by
the density power spectrum $P_{m}(k)$ at $z=7$ for three values of the grid
size $\Delta x_{\mathrm{grid}}$. Thin black lines represent the halo power
spectrum for the reference box (RB), while solid lines (scattered stars) show
the same for the hybrid box (HB) generated using the tidal tensor (density)
matching method. The lower panel shows the % error in each scaled power
spectrum with respect to the reference box.

We now turn to the two-point function of the haloes. Figure 5 shows the auto
power spectrum $P_{hh}(k)$ of the halo overdensity field scaled by the density
power spectrum $P_{m}(k)$. We find that unlike the halo mass function or the
$f_{\mathrm{coll}}-\delta_{m}$ joint distribution, the power spectra of HB
built using the closeness of the tidal tensor eigenvalues as matching criteria
agree better with the RB than those built using the closeness of the density.
This indicates the importance of the tidal environment of haloes in shaping
their large-scale distribution. We also find that the match between HB and RB
is worse for the smallest grid size $\Delta
x_{\mathrm{grid}}=0.3125h^{-1}\mathrm{cMpc}$ we use, disagreements being
$\gtrsim 10\%$ at $k\gtrsim 0.5h/$ cMpc. This is consistent with the results
of the $\delta_{m}-f_{\mathrm{coll}}$ joint distribution. For the other two
grid sizes, the accuracy is within 5% at large scales, but disagreements
increase for small scales. This is because of both the lack of power in the
density field for the LB, and also the possible discontinuities in the halo
field across the grid cells boundaries arising due to the haloes being
substituted from independent grids of the SB.

Figure 6: Cross power spectrum $P_{mh}(k)$ of the halo overdensity with matter
density, scaled by the density power spectrum $P_{m}(k)$, at $z=7$ for three
values of the grid size $\Delta x_{\mathrm{grid}}$. Thin black lines represent
the cross power spectrum for the reference box (RB), while solid lines
(scattered stars) show the same for the hybrid box (HB) generated using the
tidal tensor (density) matching method. The lower panel shows the % error in
each scaled power spectrum with respect to the reference box.

Similar conclusions can be drawn from the cross power spectra $P_{mh}(k)$ of
the density and halo fields. We show the ratio $P_{mh}(k)/P_{m}(k)$ in Figure
6. The figure shows clear agreement between HB and RB that is within $5\%$ at
small $k$ for the larger grid sizes $\Delta x_{\mathrm{grid}}=0.625$ and
$1.25h^{-1}\mathrm{cMpc}$. These results related to the power spectra imply
that the HB not only has the right number of haloes, but they are also in the
right places with respect to the density field.

###  3.3 Hybrid HI fields

We next compare the HI density field obtained from the HB with that from the
RB. It is worth emphasizing that to get the correct HI field, the halo
catalogue needs to be complete down to the lowest mass that can produce
ionizing photons, and the haloes must be correctly correlated with the
underlying dark matter density field. As indicated by our results so far, we
indeed have reasonably accurate density and halo fields to be able to build HI
maps. While comparing the results at $z=7$, we tune the ionization efficiency
$\zeta$ in HB such that the mass-averaged global neutral fraction
$x_{\mathrm{HI}}$ is equal to that in RB.

Figure 7: Power spectrum $P_{HI}(k)$ of the HI density field scaled by the
density power spectrum $P_{m}(k)$ at $z=7$. Results are shown for three mass-
averaged neutral fractions, obtained by appropriately tuning the ionization
efficiency $\zeta$ for three values of the grid size $\Delta
x_{\mathrm{grid}}$. Thin black lines represent the HI power spectrum for the
reference box (RB), while solid lines (scattered stars) show the same for the
hybrid box (HB) generated using the tidal tensor (density) matching method.
The lower panel shows the % error in each scaled power spectrum with respect
to the reference box. Figure 8: Evolution of HI and dark matter density with
redshift/neutral fraction ($x_{HI}$). The ionization efficiency $\zeta$ was
chosen such that the reference box (RB) had a mass-averaged neutral fraction
$x_{HI}=0.5$ at $z=7$. The leftmost column shows the slices of HI and dark
matter density fields for hybrid boxes generated with grid size $\Delta
x_{grid}=0.625h^{-1}$cMpc ($N_{\mathrm{grid}}=128^{3}$) using a tidal tensor-
based matching criterion. The rightmost column shows the same for a reference
simulation smoothed at the same resolution. Certain regions of interest are
zoomed in and shown in the middle two columns. Figure 9: Power spectrum
$P_{HI}(k)$ of the HI density field scaled by the density power spectrum
$P_{m}(k)$as a function of redshift. Results are shown for three wavenumbers
$k$ and for three values of $\Delta x_{\mathrm{grid}}$. The ionization
efficiency $\zeta$ was chosen such that the reference box (RB) had a mass-
averaged neutral fraction $x_{HI}=0.5$ at $z=7$. Thin black lines represent
the HI power spectrum for the reference box (RB), while solid lines (scattered
stars) show the same for the hybrid box (HB) generated using the tidal tensor
(density) matching method. The lower panel shows the % error in each scaled
power spectrum with respect to the reference box.

Figure 7 shows the HI power spectra $P_{HI}(k)$, scaled by $P_{m}(k)$, at
$z=7$, obtained using our method for three values of mass-averaged neutral
fraction. As in the case of power spectra of the haloes, we find that the HB
built using the closeness of tidal tensor eigenvalues as the matching criteria
match better with the RB. So let us focus on the results obtained using the
match in tidal tensor eigenvalues (solid lines). For values of
$x_{\mathrm{HI}}<0.75$, the agreement between the HB and RB is within $\sim
10\%$ for $k\lesssim 2.5h/$ cMpc, irrespective of the grid size used. For
$x_{\mathrm{HI}}\gtrsim 0.75$, the disagreement exceeds $10\%$ for the
smallest grid size $\Delta x_{\mathrm{grid}}=0.3125h^{-1}\mathrm{cMpc}$. The
disagreement seems to be noisy and could arise due to fewer ionized bubbles at
these stages. This shows that there is some scope for improving the method for
early stages of reionization and for smaller grid sizes. However, for the
larger two grid sizes, the match is within $10\%$. We also find that for the
largest grid size $\Delta x_{\mathrm{grid}}=1.25h^{-1}\mathrm{cMpc}$ used in
this work, the match is within $5\%$ if we limit to relatively later stages of
reionization $x_{\mathrm{HI}}\lesssim 0.75$. In summary, we find that it can
provide HI power spectra accurate to within $10\%$ for $k\lesssim 2.5h/$cMpc
for grid size $\Delta x_{\mathrm{grid}}\geq 0.625h^{-1}\mathrm{cMpc}$.

We next check the performance of our method at other redshifts. We choose five
redshifts $z=9,8,7,6,5$. The efficiency $\zeta$ is fixed at a value that leads
to $50\%$ ionization at $z=7$ in the RB. This leads to a reionization history
with $x_{\mathrm{HI}}=0.86,0.73,0.5,0.14,0.0$ at $z=9,8,7,6,5$, respectively.
The evolution of the HI density field is shown in Figure 8 for $\Delta
x_{\mathrm{grid}}=0.625h^{-1}\mathrm{cMpc}$ for both HB and RB. As is clear
from the comparison, the HB captures the large-scale distribution quite
faithfully. There are some differences at smaller scales, which is expected as
the halo populations in HB and RB are not identical.

Finally, we show the evolution of $P_{HI}(k)/P_{m}(k)$ with redshift for three
values of $k$ in Figure 9. We again find that their values for the HB and RB
agree to within $10\%$ for all the redshifts for $\Delta x_{\mathrm{grid}}\geq
0.625h^{-1}\mathrm{cMpc}$ when we use the closeness of tidal-tensor as the
matching criterion for grid cells in LB and SB.

###  3.4 Computational Advantage

The most significant advantage of our hybrid method is that we can obtain
high-dynamic range information from low-dynamic range simulation boxes having
much fewer particles, therefore requiring much less memory to run. The peak
memory requirement for building the hybrid box from parent boxes is designed
to never exceed that of loading one of those parent boxes in memory.
Therefore, we can build the entire hybrid box using a computer that can only
simulate

$\frac{N_{\mathrm{part}}}{N_{\mathrm{part,large}}}=\frac{N_{\mathrm{part}}}{N_{\mathrm{part,small}}}$
(3.1)

($=8$ in our case) times fewer particles.

Currently, our hybrid method requires parent simulations that are run with 8
times fewer particles than what a full simulation would require, thus reducing
the memory requirement of GADGET-2 simulations by 87%. One could get away with
an even smaller amount of memory at the cost of increasing the errors in the
HI power spectra.

Since the LB and SB both contain fewer particles than RB, the computational
time taken for constructing the HB is significantly smaller. As mentioned in
Table 1, we can generate the LB and SB in about $7$ times less time compared
to the RB. Compared to the amount of time it takes to run the parent
simulations, our hybrid method is nearly instantaneous, taking only a few
minutes to run on a single CPU core.

##  4 Summary

Numerical studies of reionization can be quite expensive due to the high
dynamic range demanded by such simulations [1, 2]. It is therefore challenging
to explore the space of unknown parameters when fitting theoretical models to
data. While this challenge has been addressed to a large extent via
computationally efficient semi-numerical simulations, they still have the
shortcoming that small haloes, potentially hosting the first stars, cannot be
resolved directly [7, 8, 9, 10, 11, 12, 13, 14]. This problem is usually
addressed using some sub-grid prescriptions to populate the simulation volume
with unresolved haloes, e.g., those based on an analytical conditional halo
mass function [15, 16, 21, 20, 8, 22, 23, 24].

In this work, we take an alternate hybrid approach to populate a low-
resolution simulation volume with sub-grid haloes by combining information
from two simulation boxes having relatively narrow dynamic ranges. Our hybrid
scheme uses the density field from a coarse resolution simulation, having a
size that can access the large scales required in the problem. The low-mass
haloes cannot be resolved in this simulation, so another high-resolution
simulation of smaller volume is run to resolve these haloes. We then match
grid cells from the large box and the small box based on their physical
properties (e.g., density and/or tidal environment) and populate haloes from
the small box to the large one. Our method can provide an effective high
dynamic range simulation by combining two boxes having relatively narrow
dynamic ranges. An advantage of our method is that it does not require
assuming any functional form for the underlying halo mass function. Instead,
it is derived self-consistently from the simulations used, minimizing any
uncertainty arising from the choice of the functional form [35].

While finding the matching cells between the large and small boxes, we explore
two different criteria. The first is to find a matching cell that is closest
in the density contrast $\delta_{m}$, and the second method uses the match
between the eigenvalues of the tidal tensor. Our analyses show that the second
method provides a better match between the hybrid and reference box, thus
highlighting the importance of the tidal environment in determining the large-
scale distribution of haloes.

We find that, compared to a reference simulation done using the cosmological
simulation GADGET-2 and semi-numerical code SCRIPT, our hybrid method
accurately reproduces the HI maps as well as summary statistics relevant to
the studies of reionization. The hybrid halo mass function, large-scale halo
auto-power spectra and halo-matter cross-power spectra, and the HI power
spectra are all reproduced within $\sim 10\%$ at scales $k\lesssim 2.5h/$cMpc.
Our method gives accurate results at a range of redshifts relevant to the
studies of reionization, enabling detailed studies of the temporal evolution
of the 21 cm signal.

Our method is at least 7 times faster than a full GADGET-2 simulation having
the equivalent dynamic range. More importantly, the hybrid box can be built
with only $13\%$ of the memory required for a full simulation. In principle,
our method can be even less expensive if one is only interested in the really
large scales, and can afford some loss of accuracy at smaller scales. In such
a case, boxes with 4 times narrower dynamic range (as opposed to 2 times
narrower boxes that we have analyzed in this paper) can be explored. The
caveat to consider in such a case is that the small box may not sample the
extreme densities that are present in the large box, thus finding
appropriately matched grid cells could be difficult. In addition, the largest
halo mass reliably sampled in the small box might be smaller than the smallest
halo resolved in the large box, thus leading to a gap in the hybrid halo mass
function. Further analyses are required to quantify these effects.

Encouraged by our success in accurately reproducing the two-point statistics
of dark matter, halo, and HI density fields, we wish to extend our approach to
incorporate the effects of redshift-space distortions on the power spectra.
This would require building a hybrid velocity field that is well correlated
with the density, halo, and HI fields, and we leave it for future work. We are
also interested in measuring how well our hybrid method reproduces higher-
order statistics such as three-point correlations of these fields.

The efficiency of our method enables the creation of reasonably accurate halo
catalogues with relatively fewer computational resources, which can be used to
generate large-scale HI maps. This method can be particularly useful for
varying cosmological parameters with semi-numerical reionization simulations,
a traditionally challenging task [36]. With recent advancements in techniques
like simulation-based inference for constraining reionization models [37, 38,
39, 40], it is worth exploring how upcoming reionization observations can be
leveraged to probe cosmological parameters effectively using our hybridization
method.

## Acknowledgments

TRC acknowledges the insightful discussions with Aditya Chowdhury and Girish
Kulkarni during the initial stages of this work.

## Data Availability

The data generated during this work will be made available upon reasonable
request to the corresponding author.

## References

  * [1] N.Y. Gnedin and P. Madau, _Modeling cosmic reionization_ , _Living Reviews in Computational Astrophysics_ 8 (2022) 3 [2208.02260].
  * [2] T.R. Choudhury, _A short introduction to reionization physics_ , _General Relativity and Gravitation_ 54 (2022) 102 [2209.08558].
  * [3] I.T. Iliev, G. Mellema, K. Ahn, P.R. Shapiro, Y. Mao and U.-L. Pen, _Simulating cosmic reionization: how large a volume is large enough?_ , _MNRAS_ 439 (2014) 725 [1310.7463].
  * [4] H.D. Kaur, N. Gillet and A. Mesinger, _Minimum size of 21-cm simulations_ , _MNRAS_ 495 (2020) 2354 [2004.06709].
  * [5] S.K. Giri, A. Schneider, F. Maion and R.E. Angulo, _Suppressing variance in 21 cm signal simulations during reionization_ , _A &A_ 669 (2023) A6 [2209.01225].
  * [6] A. Acharya, E. Garaldi, B. Ciardi and Q.-b. Ma, _Cosmic variance suppression in radiation-hydrodynamic modelling of the reionization-era 21-cm signal_ , _MNRAS_ 529 (2024) 3793 [2310.13401].
  * [7] A. Mesinger and S. Furlanetto, _Efficient Simulations of Early Structure Formation and Reionization_ , _ApJ_ 669 (2007) 663 [0704.0946].
  * [8] A. Mesinger, S. Furlanetto and R. Cen, _21CMFAST: a fast, seminumerical simulation of the high-redshift 21-cm signal_ , _MNRAS_ 411 (2011) 955 [1003.3878].
  * [9] M.G. Santos, L. Ferramacho, M.B. Silva, A. Amblard and A. Cooray, _Fast large volume simulations of the 21-cm signal from the reionization and pre-reionization epochs_ , _MNRAS_ 406 (2010) 2421 [0911.2219].
  * [10] S. Majumdar, G. Mellema, K.K. Datta, H. Jensen, T.R. Choudhury, S. Bharadwaj et al., _On the use of seminumerical simulations in predicting the 21-cm signal from the epoch of reionization_ , _MNRAS_ 443 (2014) 2843 [1403.0941].
  * [11] S. Hassan, R. Davé, K. Finlator and M.G. Santos, _Simulating the 21 cm signal from reionization including non-linear ionizations and inhomogeneous recombinations_ , _MNRAS_ 457 (2016) 1550 [1510.04280].
  * [12] R. Mondal, S. Bharadwaj and S. Majumdar, _Statistics of the epoch of reionization (EoR) 21-cm signal - II. The evolution of the power-spectrum error-covariance_ , _MNRAS_ 464 (2017) 2992 [1606.03874].
  * [13] A. Hutter, _The accuracy of seminumerical reionization models in comparison with radiative transfer simulations_ , _MNRAS_ 477 (2018) 1549 [1803.00088].
  * [14] T.R. Choudhury and A. Paranjape, _Photon number conservation and the large-scale 21 cm power spectrum in seminumerical models of reionization_ , _MNRAS_ 481 (2018) 3821 [1807.00836].
  * [15] T.R. Choudhury, M.G. Haehnelt and J. Regan, _Inside-out or outside-in: the topology of reionization in the photon-starved regime suggested by Ly $\alpha$ forest data_, _MNRAS_ 394 (2009) 960 [0806.1524].
  * [16] E. Sobacchi and A. Mesinger, _Inhomogeneous recombinations during cosmic reionization_ , _MNRAS_ 440 (2014) 1662 [1402.2298].
  * [17] J.R. Bond, S. Cole, G. Efstathiou and N. Kaiser, _Excursion Set Mass Functions for Hierarchical Gaussian Fluctuations_ , _ApJ_ 379 (1991) 440.
  * [18] R.K. Sheth and G. Tormen, _Large-scale bias and the peak background split_ , _MNRAS_ 308 (1999) 119 [astro-ph/9901122].
  * [19] R. Barkana and A. Loeb, _Unusually Large Fluctuations in the Statistics of Galaxy Formation at High Redshift_ , _ApJ_ 609 (2004) 474 [astro-ph/0310338].
  * [20] O. Zahn, A. Lidz, M. McQuinn, S. Dutta, L. Hernquist, M. Zaldarriaga et al., _Simulations and Analytic Calculations of Bubble Growth during Hydrogen Reionization_ , _ApJ_ 654 (2007) 12 [astro-ph/0604177].
  * [21] R. Barkana and A. Loeb, _The difference PDF of 21-cm fluctuations: a powerful statistical tool for probing cosmic reionization_ , _MNRAS_ 384 (2008) 1069 [0705.3246].
  * [22] M.G. Santos, M.B. Silva, J.R. Pritchard, R. Cen and A. Cooray, _Probing the first galaxies with the Square Kilometer Array_ , _A &A_ 527 (2011) A93 [1009.0950].
  * [23] R. Ghara, T.R. Choudhury and K.K. Datta, _21 cm signal from cosmic dawn: imprints of spin temperature fluctuations and peculiar velocities_ , _MNRAS_ 447 (2015) 1806 [1406.4157].
  * [24] A. Doussot and B. Semelin, _A bubble size distribution model for the Epoch of Reionization_ , _A &A_ 667 (2022) A118 [2208.14044].
  * [25] K. Ahn, I.T. Iliev, P.R. Shapiro, G. Mellema, J. Koda and Y. Mao, _Detecting the Rise and Fall of the First Stars by Their Impact on Cosmic Reionization_ , _ApJ_ 756 (2012) L16 [1206.5007].
  * [26] K. Ahn, I.T. Iliev, P.R. Shapiro and C. Srisawat, _Non-linear bias of cosmological halo formation in the early universe_ , _MNRAS_ 450 (2015) 1486 [1407.2637].
  * [27] R.K. Sheth and G. Tormen, _An excursion set model of hierarchical clustering: ellipsoidal collapse and the moving barrier_ , _MNRAS_ 329 (2002) 61 [astro-ph/0105113].
  * [28] A. Mesinger, A. Aykutalp, E. Vanzella, L. Pentericci, A. Ferrara and M. Dijkstra, _Can the intergalactic medium cause a rapid drop in Ly $\alpha$ emission at z > 6?_, _MNRAS_ 446 (2015) 566 [1406.6373].
  * [29] T.R. Choudhury, E. Puchwein, M.G. Haehnelt and J.S. Bolton, _Lyman $\alpha$ emitters gone missing: evidence for late reionization?_, _MNRAS_ 452 (2015) 261 [1412.4790].
  * [30] M. Davis, G. Efstathiou, C.S. Frenk and S.D.M. White, _The evolution of large-scale structure in a universe dominated by cold dark matter_ , _ApJ_ 292 (1985) 371.
  * [31] V. Springel, _The cosmological simulation code GADGET-2_ , _MNRAS_ 364 (2005) 1105 [astro-ph/0505010].
  * [32] V. Springel, “N-GenIC: Cosmological structure initial conditions.” Astrophysics Source Code Library, record ascl:1502.003, feb, 2015.
  * [33] D.J. Eisenstein and W. Hu, _Power Spectra for Cold Dark Matter and Its Variants_ , _ApJ_ 511 (1999) 5 [astro-ph/9710252].
  * [34] T.R. Choudhury and A. Paranjape, “SCRIPT: Semi-numerical Code for ReIonization with PhoTon-conservation.” Astrophysics Source Code Library, record ascl:2204.013, apr, 2022.
  * [35] B. Greig, D. Prelogović, J. Mirocha, Y. Qin, Y.-S. Ting and A. Mesinger, _Exploring the role of the halo mass function for inferring astrophysical parameters during reionisation_ , _arXiv e-prints_ (2024) arXiv:2403.14061 [2403.14061].
  * [36] N.S. Kern, A. Liu, A.R. Parsons, A. Mesinger and B. Greig, _Emulating simulations of cosmic dawn for 21 cm power spectrum constraints on cosmology, reionization, and x-ray heating_ , _The Astrophysical Journal_ 848 (2017) 23.
  * [37] X. Zhao, Y. Mao, C. Cheng and B.D. Wandelt, _Simulation-based Inference of Reionization Parameters from 3D Tomographic 21 cm Light-cone Images_ , _ApJ_ 926 (2022) 151 [2105.03344].
  * [38] X. Zhao, Y. Mao and B.D. Wandelt, _Implicit Likelihood Inference of Reionization Parameters from the 21 cm Power Spectrum_ , _ApJ_ 933 (2022) 236 [2203.15734].
  * [39] D. Prelogović and A. Mesinger, _Exploring the likelihood of the 21-cm power spectrum with simulation-based inference_ , _MNRAS_ 524 (2023) 4239 [2305.03074].
  * [40] A. Saxena, A. Cole, S. Gazagnes, P.D. Meerburg, C. Weniger and S.J. Witte, _Constraining the X-ray heating and reionization using 21-cm power spectra with Marginal Neural Ratio Estimation_ , _MNRAS_ 525 (2023) 6097 [2303.07339].