# The short ionizing photon mean free path at z=6 in Cosmic Dawn III, a new
fully-coupled radiation-hydrodynamical simulation of the Epoch of Reionization

Joseph S. W. Lewis1,2, Pierre Ocvirk3, Jenny G. Sorce4,5, Yohan Dubois6,
Dominique Aubert3, Luke Conaboy7, Paul R. Shapiro8, Taha Dawoodbhoy8, Romain
Teyssier9, Gustavo Yepes10,11, Stefan Gottlöber4, Yann Rasera2, Kyungjin
Ahn13, Ilian T. Iliev7, Hyunbae Park14,15,16, Émilie Thélie3
1Zentrum für Astronomie der Universität Heidelberg, Institut für Theoretische
Astrophysik, Albert-Ueberle-Straße 2, 69120 Heidelberg, Germany
2Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg,
Germany
3Observatoire Astronomique de Strasbourg, Université de Strasbourg, CNRS UMR
7550, 11 rue de l’Université, 67000 Strasbourg, France
4Leibniz-Institut für Astrophysik Potsdam (AIP), An der Sternwarte 16, D-14482
Potsdam, Germany
5Université Paris-Saclay, CNRS, Institut d’Astrophysique Spatiale, 91405,
Orsay, France
6Institut d’Astrophysique de Paris, UMR 7095, CNRS, UPMC Univ. Paris VI, 98
bis boulevard Arago, 75014 Paris, France
7Astronomy Center, Department of Physics & Astronomy, Pevensey II Building,
University of Sussex, Falmer, Brighton BN1 9QH, United Kingdom
8Department of Astronomy, University Texas, Austin, TX 78712-1083, USA
9Institute for Theoretical Physics, University of Zurich, Winterthurerstrasse
190, CH-8057 Zürich, Switzerland
10 Departamento de Física Teórica M-8, Universidad Autónoma de Madrid,
Cantoblanco, 28049, Madrid, Spain
11 Centro de Investigación Avanzada en Física Fundamental (CIAFF), Universidad
Autónoma de Madrid, 28049 Madrid, Spain
12Laboratoire Univers et Théories, Université de Paris,Observatoire de Paris,
Université PSL, CNRS, F-92190 Meudon, France
13Chosun University, 375 Seosuk-dong, Dong-gu, Gwangjiu 501-759, Korea
14Lawrence Berkeley National Laboratory, CA 94720-8139, USA
15Berkeley Center for Cosmological Physics, UC Berkeley, CA 94720, USA
16Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583,
Japan

(Accepted 2022 August 08. Received 2022 August 05; in original form 2022
February 11)

###### Abstract

Recent determinations of the mean free path of ionising photons (mfp) in the
intergalactic medium (IGM) at $\rm z=6$ are lower than many theoretical
predictions. In order to gain insight, we investigate the evolution of the mfp
in our new massive fully coupled radiation hydrodynamics cosmological
simulation of reionization: Cosmic Dawn III (CoDa III). CoDa III’s scale ($\rm
94^{3}\,cMpc^{3}$) and resolution ($\rm 8192^{3}$ grid) make it particularly
suitable to study the IGM during reionization. The simulation was performed
with RAMSES-CUDATON on Summit, and used 131072 processors coupled to 24576
GPUs, making it the largest reionization simulation, and largest ever RAMSES
simulation. A superior agreement with global constraints on reionization is
obtained in CoDa III over CoDa II, especially for the evolution of the neutral
hydrogen fraction and the cosmic photo-ionization rate, thanks to an improved
calibration, later end of reionization ($\rm z=5.6$), and higher spatial
resolution. Analyzing the mfp, we find that CoDa III reproduces the most
recent observations very well, from $\rm z=6$ to $\rm z=4.6$. We show that the
distribution of the mfp in CoDa III is bimodal, with short (neutral) and long
(ionized) mfp modes, due to the patchiness of reionization and the co-
existence of neutral versus ionized regions during reionization. The neutral
mode peaks at sub-kpc to kpc scales of mfp, while the ionized mode peak
evolves from $\rm 0.1Mpc/h$ at $\rm z=7$ to $\sim 10$ Mpc/h at $\rm z=5.2$.
Computing the mfp as the average of the ionized mode provides the best match
to the recent observational determinations. The distribution reduces to a
single neutral (ionized) mode at $\rm z>13$ ($\rm z<5$).

###### keywords:

cosmology: reionization – galaxies: high-redshift

††pubyear: 2022††pagerange: The short ionizing photon mean free path at z=6 in
Cosmic Dawn III, a new fully-coupled radiation-hydrodynamical simulation of
the Epoch of Reionization–A

##  1 Introduction

The Epoch of Reionization (EoR hereafter) began as the first ionising galactic
sources started to ionize their immediate surroundings in the inter-galactic
medium (IGM) some few hundred million years after the Big Bang. This epoch
lasted approximately $600\text{\,}\mathrm{Myr}$, until the reionization of the
IGM was complete at some point between $\rm z\approx 6$ and $\rm z\approx 5.3$
(Bosman et al., 2022, Bo22 hereafter). Constraining this distant epoch offers
unique insight into the formation of the first galaxies and stars, and is an
important step in bridging the gap between our knowledge of ’modern’ galaxies,
and the first moments after the Big Bang. Indeed, the rising and expanding UV
background during the EoR may suppress star formation in low mass galaxies
(Ocvirk et al., 2016, 2020; Dawoodbhoy et al., 2018, the former two will be
called Oc16 and Oc20 from hereon). This in turn motivated efforts to
investigate and predict possible signatures of the EoR in populations of
galaxies of the lowest mass, which are best observed in the Local Group
(Ocvirk & Aubert, 2011; Iliev et al., 2011; Ocvirk et al., 2013, 2014; Dixon
et al., 2018).

Due to the very nature of the EoR and the extreme distances involved, direct
observations of the reionizing sources as they drive reionization are
extremely challenging. Despite this, over the past decades innovative indirect
constraints on both the sources of reionization and on the IGM have emerged.
Further, the future is bright as upcoming space and ground based observatories
(e.g. the James Webb Space Telescope, the Nancy Grace Roman Space Telescope,
or the Extremely Large Telescope) will improve our capability to target the
drivers of reionization. At the same time, current and planned 21cm radio
observations (e.g. EDGES , HERA, LOFAR, NenuFAR, SKA) may bring new
constraints on the ionisation of the IGM during the EoR and at cosmic dawn.
Thus, the study of the EoR appears to be set for rapid progress over the
coming years. In preparation of this new era of EoR observations, the
simulation community is hard at work making predictions for these new
observational campaigns, and investigating the physical processes that drive
the sources of reionization. However, the predictive power of these
simulations relies on the quality of their agreement with existing
observational constraints, some of which have proved difficult to match.

For instance, recent determinations of the mean free path of ionising photons
(mfp hereinafter) , at a higher redshift than ever before, have shown that
models and simulations of reionization overestimate the mfp when $\rm z\approx
6$ (and possibly when $\rm 5.24\lesssim z\lesssim 6$, see Rahmati & Schaye,
2018; D’Aloisio et al., 2020; Keating et al., 2020), casting doubt on the
realism of their portrayal of the ionisation state of the IGM and thus cosmic
reionization as a whole (See Becker et al., 2021, Be21 from hereon). Already,
new work has investigated the implications of these new data, finding that
they require late and rapid reionization scenarios (Davies et al., 2021; Cain
et al., 2021, Da21 and Ca21 from hereon), which are also favoured by
constraints on the $\rm Lyman\mathchar 45\relax\alpha$ forest (Kulkarni et
al., 2019; Bosman et al., 2022). In fact, Ca21 have been able to obtain a
reasonable agreement with both the new and old mfp constraints. We consider
below these two recent works attempting to match the new very low mfp
measurement of Be21 and show how using Cosmic Dawn III (CoDa III )111More
information can be found about the whole simulation suite and the CoDa project
here: https://coda-simulation.github.io/. may improve upon them. Both studies
above use a fairly low resolution simulation of the IGM ($\rm\sim 1cMpc$ cell
size), and must recourse to complex sub-grid models, to:

  * •

model the source population: Da21’s description of sources involves a range of
parameters which fold in a number of complex and coupled aspects, such as the
sensitivity of sources to SN and AGN feedback, their star formation
efficiency, their escape fractions and the global dependence of these
processes with respect to mass. However, as we show in Ocvirk et al.
(2021)(Oc21 hereafter), radiative suppression of star formation in low mass
haloes is crucial to obtaining a realistic Lyman-alpha forest after overlap.
This regulation, for instance, is not included in Da21, and CoDa III improves
on their approach by allowing this process to occur self-consistently thanks
to the fully coupled radiation hydrodynamics (RHD) framework employed.

  * •

set the transmission properties of the IGM sub-grid: Ca21 uses a sub-grid
parameterization of the mean free path, which relies on simulations of
irradiated patches (D’Aloisio et al., 2020). Such simulations do not include
star formation, and one can therefore estimate that their densest patch
structures would at some point form stars and therefore, both decrease their
own mean free path and the mean free path around them.

Our approach circumvents the need for using such a sub-grid model of the mean
free path, as the whole volume of CoDa III is eligible to star formation
provided the relevant criteria in density and temperature are met. Moreover,
the mfp is self-consistently set in CoDa III by the physical state of the IGM,
and is self-consistently affected by radiation received by sources near and
far and the high-resolution hydrodynamical description of the IGM.

Beyond these two studies, it is desirable, in all cases, to deploy a fully
coupled RHD formalism as we do here, in order to break free from the
limitations imposed by the parameterization of the physical processes at hand
as in ; , and check that the results they claim can be confirmed in more self-
consistent simulations where more physics are resolved and therefore do not
require such a level of sub-grid modelling. In this direction, the THESAN
project (Kannan et al., 2021; Garaldi et al., 2021) has recently performed a
series of AREPO-RT (Kannan et al., 2019) simulations of the EoR. While they
reach high spatial resolution inside of galaxies, CoDa III improves on their
approach in at least 2 aspects: their fiducial run, THESAN-1, stops at z=5.5,
whereas we were able to go beyond this redshift, down to z=4.6, and as a
result compare CoDa III mfp with observations in a larger redshift range,
which, most importantly, includes post-reionization measurements as well.
Also, THESAN-1 uses a typical speed of light reduction of 0.2, which can have
a strong impact on the post-reionization residual neutral hydrogen fraction
(Ocvirk et al., 2019), and to some extent change the timing of reionization
(Deparis et al., 2019), complicating the interpretation of their simulation
results. To avoid such difficulties, Cosmic Dawn III uses the full speed of
light thanks to the GPU-optimized radiative transfer module CUDATON ( ; , ).

First, we present our methodology, including the particularities of our
simulation code RAMSES-CUDATION, the simulation setup, and the computational
steps we take to derive the mfp. Then, we present the global properties of
reionization in the CoDa III simulation, and the evolution of the mfp in CoDa
III. Finally we conclude our work.

##  2 Methods

###  2.1 Simulation

#### Overview of RAMSES-CUDATON

The CoDa III (Cosmic Dawn III) simulation was performed using the fully
coupled radiation and hydrodynamics code RAMSES-CUDATON (presented in ; , ).
Here we briefly recall its main properties as well as details concerning its
deployment.

RAMSES-CUDATON results from the coupling of the RAMSES code for N-body dark
matter dynamics, hydrodynamics, star formation and stellar feedback (Teyssier,
2002), with the radiative transfer module ATON (Aubert & Teyssier, 2008),
which handles the radiative transfer of ionising photons via the M1 method
(Levermore, 1984), hydrogen photo-ionisation and thermo-chemistry. Motivated
by the important computational cost of radiative transfer, we accelerate the
ATON module using GPUs (Aubert & Teyssier, 2010), whilst the remaining RAMSES
physics run on CPUs. The performance boost thus obtained allows us to run
fully coupled simulations with a full speed of light, circumventing pitfalls
encountered by certain reduced speed of light approaches (Gnedin, 2016;
Deparis et al., 2019), in particular on the residual neutral fraction after
overlap (Ocvirk et al., 2019).

#### Improvements over our previous simulation CoDa II

For the most part, the code and simulation setup are as described in Oc20, but
with higher resolution. Other features of CoDa III which differ from past CoDa
simulations were fully tested and calibrated by a large suite of smaller
simulations of the same numerical resolution as CoDa III to ensure a good
match to observational constraints on the EoR. While the chief purpose of
these calibration runs was to prepare for CoDa III, the calibration runs
themselves were novel enough to warrant new science as described in Ocvirk et
al. (2021). Since the new features of CoDa III were already described there,
we refer the reader to that paper for more detail. We quickly enumerate the
novelties of CoDa III with respect to CoDa II hereafter.

  * •

CoDa III uses a $8192^{3}$ grid while CoDa II had a $4096^{3}$ grid, for the
same volume, meaning its physical resolution is a factor of two finer, and its
mass resolution is eight times higher.

  * •

We take a more detailed approach with respect to our source model than in the
previous CoDa simulations. In CoDa III, each stellar particle receives an
ionising emissivity value based on its mass, age, and metallicity. We pre-
compute emissivities using BPASS V2.2.1 Eldridge & Stanway (2020), considering
that each stellar particle is a star cluster of about
${10}^{4}\text{\,}\mathrm{M_{\odot}}$ forming over a few Myr as in Oc21. As a
result, in CoDa III individual stellar particles’ ionising luminosities
decrease more gradually and realistically than in CoDa II. 222In the previous
CoDa simulations, stellar particles had a fixed stellar luminosity until their
massive stars went supernova, at which point their ionising luminosities were
set to 0.

  * •

One of the main model differences with respect to CoDa II lies within the star
formation sub-grid model. In CoDa III, for a gas cell to form stars, we not
only require that it be dense enough, but also that it be cooler than
$2\text{\times}{10}^{4}\text{\,}\mathrm{K}$. This additional temperature
threshold for star formation has been widely used in similar simulations
before CoDa III including CoDa I (Ocvirk et al., 2016) and others (e.g. Dubois
& Teyssier, 2008). Its physical interpretation is as follows: gas cells hotter
than the threshold are so because of shocks, usually caused by supernovae,
which are not favoured sites for star formation until they manage to cool back
down to below the $2\text{\times}{10}^{4}\text{\,}\mathrm{K}$ threshold.

We also found that including the temperature threshold allowed our simulations
to achieve better agreement with the residual neutral fraction at the end of
the EoR and with the ionisation rate of ionized regions, thus advocating its
use in CoDa III. For more information and a comparison between star formation
with, and without the temperature threshold, we refer the reader to Oc21.

  * •

We use a sub-grid stellar escape fraction of $\rm f_{esc}^{sub}=1$, whereas we
used $\rm f_{esc}^{sub}=0.42$ in CoDa II, i.e. all photons produced by the
stellar particles in CoDa III are released in the host cell. We shall fully
investigate the escape of ionising photons in an upcoming paper, however we
can already tentatively explain why we require a higher $\rm f_{esc}^{sub}$ in
CoDa III: First, and due to the changes to the star formation model, star
formation in CoDa III is not permitted in very hot ionised cells. For the same
amount of star formation it is then simple to imagine that one requires more
photons to complete Reionization in CoDa III with respect to CoDa II. Second,
CoDa III’s higher resolution likely gives rises to more and denser clumps of
gas in the IGM, potentially further raising the required number of photons to
complete Reionization. Third, the ionising emissivity of stellar particles
evolves in a smoother fashion in CoDa III, plausibly altering the balance
between ionization and recombination (although we found this effect to be
relatively moderate in smaller test volumes).

  * •

Just as in CoDa I & CoDa II (the two previous Cosmic Dawn simulation
projects), upon reaching $10\text{\,}\mathrm{Myr}$ of age, a mass fraction
$\eta_{\rm SN}=0.2$ of stellar particles detonates as supernovae. Each
supernova event injects $51\text{\,}\mathrm{ergs}$ of energy for every
$10\text{\,}\mathrm{M_{\odot}}$ of progenitor mass into its host cell (using
the kinetic feedback of Dubois & Teyssier, 2008). After the supernova event, a
long lived particle of mass $\rm(1-\eta_{SN})M_{birth}$ remains. Unlike CoDa I
and CoDa II, CoDa III features the standard RAMSES chemical enrichment333for
clarification: this is not a new code feature, but was not activated at run
time in the previous simulations, for simplicity.. Metals are produced by
supernovae. A fraction $\rm y=0.05$ of the total supernova mass is re-injected
as metals. Metallicity is treated as a passive scalar, and metals are thus
advected along with the gas.

  * •

CoDa III features a physically motivated model for the formation and
destruction of dust in galaxies on the fly, that we coupled to the ionizing
radiative transfer of ATON (see Lewis et al., 2022, for details). An
independent paper on the dust model is also in preparation (Dubois et al.).
The dust aspects of CoDa III will be analysed in a forthcoming paper, and can
be left aside for the present study.

For quick reference, Table 1 recaps the essential physical and numerical
parameters of CoDa III.

#### Initial conditions & Cosmology

As in previous CoDa simulations, CoDa III ’s initial conditions were
constrained so as to produce facsimiles of the progenitors of the Local Group
(See Sorce et al., 2016; Sorce & Tempel, 2018) within the Constrained Local
Universe Simulations project (CLUES444https://www.clues-project.org/cms/). We
will use this unique opportunity to investigate the reionization of the Local
Group in a forthcoming paper. For the present study and for the purpose of
this paper, we can forgo the Local Universe aspect and treat CoDa III as a
random realization of the Universe, as its mass function, for instance, is
compatible with that of a random patch of Universe. The cosmology used
($\Omega_{\Lambda}=0.693$, $\rm\Omega_{m}=0.307$, $\rm\Omega_{b}=0.0482$, $\rm
H_{0}=$
$67.77\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{Mpc}}^{-1}$,
$\rm\sigma_{s}=0.829$, $\rm n=0.963$) here is taken from the Planck 2014
constraints ( Planck14, )555these constraints are compatible with subsequent
Planck publications, as summarised in Table 1. The initial conditions were
generated at $\rm z=150$.

#### Deployment

CoDa III was run on Summit at the Oak Ridge National Laboratory / Oak Ridge
Leadership Computing Facility. It was deployed on 131072 processors coupled
with 24576 GPUs, distributed on 4096 nodes, and ran for $\sim$ 10 days,
allowing CoDa III to reach a final redshift of z=4.6, i.e. sufficiently late
after overlap that we can consider reionization complete and sample well into
the post-overlap phase. Snapshots were written every $\sim 10$ Myr, we use a
subset of them to derive the properties of the IGM in the following section.
Subsequent processing of the more than 20 PB of data outputs took place at the
same facility on Andes and Summit.

Setup
---
Grid size & Dark matter particle number | $8192^{3}$
Box size | $94.43\text{\,}\mathrm{cMpc}$
Force resolution |
…comoving | $11.53\text{\,}\mathrm{ckpc}$
…physical ($\rm z=6$) | $1.65\text{\,}\mathrm{kpc}$
Dark matter particle mass | $5.09\text{\times}{10}^{4}\text{\,}\mathrm{M_{\odot}}$
Average cell gas mass | $9.375\text{\times}{10}^{3}\text{\,}\mathrm{M_{\odot}}$
Stellar particle mass | $11\,732\text{\,}\mathrm{M_{\odot}}$
Star formation & feedback
Density threshold for star formation | 50<$\rm\rho_{gas}$>
Temperature threshold for star formation | $2\text{\times}{10}^{4}\text{\,}\mathrm{K}$
Star formation efficiency $\rm\epsilon_{\star}$ | 0.03
Massive star lifetime | $10\text{\,}\mathrm{Myr}$
Supernova energy | ${10}^{51}\text{\,}\mathrm{ergs}$
Supernova mass fraction, $\rm\eta_{SN}$ | 0.2
Supernova ejecta metal mass fraction | 0.05
Radiation
Stellar ionising emissivity model | BPASS V2.2.1 binary
(from Eldridge & Stanway, 2020) |
Stellar particle sub-grid escape fraction, $\rm f_{esc}^{sub}$ | $1.0$
Effective photon energy | $20.28\text{\,}\mathrm{eV}$
Effective HI cross-section (at $20.28\text{\,}\mathrm{eV}$) | $2.493\text{\times}{10}^{-22}\text{\,}{\mathrm{m}}^{2}$
Reionization
CMB electron scattering optical depth, $\rm\tau_{es}(z=14)$ | 0.0497
reionization mid-point, $\rm z_{re,\,50\%}$ | 6.81
reionization complete, $\rm z_{re,\,99.99\%}$ | 5.53
Table 1: Overview of the CoDa III simulation.

###  2.2 Computing the mean free path of ionising photons

To compute the mfp, we employ four approaches:

  * •

First, we employ a method (called average transmission method from hereon)
inspired by Kulkarni et al. (2016) . For a given simulation snapshot, we draw
1024 lines of sight (LoS) across the whole simulation box in one direction.
Our LoS are evenly separated so as to evenly sample the box.

For each LoS, we then compute the transmitted flux at
$912\text{\,}\mathrm{\textup{\AA}}$ as a function of position along the LoS x.
These can then be averaged into one average transmitted flux F(x), from which
we obtain $\rm\lambda_{\rm AT}$ by fitting it with Eq. 1 for the parameter
$\lambda$:

$\begin{array}[]{lll}\rm F(x)&=&\rm
F_{0}exp\Big{(}-\frac{x}{\lambda}\Big{)}\end{array}$ (1)

where $\rm F_{0}$ is a normalisation parameter, equivalent to a flux incident
where $\rm x=0$, that is free during the fit. For key redshifts, we check that
our results are converged with respect to the number of lines of sight, and
find no significant difference in our result with 8x or 64x more lines of
sight. In order to mimic the small number statistics regime of the
observations, we draw samples of only 13 LoS (the same LoS sample size as in
Be21, at z=6) from our full sample of 1024. We repeat this 105 times, and our
final $\rm\lambda_{\rm AT}$ result is then the median of these $10^{5}$
realizations. Also of interest, this approach yields the standard deviation of
this set of mfp.

  * •

Second, in order to access the distribution of mfp belonging to individual
lines of sight ($\rm\lambda_{\rm los}$ hereafter), we fit Eq. 1 to the
transmitted flux of every LoS individually. We can then also compute the mean,
median, and spread of the resulting distribution at every redshift. The
results from bad fits (reduced $\rm\chi^{2}>1.0$) are discarded. This occurs
only 5 times for the $1024\times 14$ LoS that we process, and correspond to
lines of sight that cross a large overdensity. In Appendix A, we compare this
method to that of Rahmati & Schaye (2018) (also found in Garaldi et al.,
2021), and show that both methods for computing the mean free path along
individual LoS are remarkably consistent.

  * •

Third, as we will see in Sec. 3.2, our $\rm\lambda_{\rm los}$ are divided into
two families that reside in neutral, and ionized areas of the IGM. We can
easily isolate the evolution of the $\rm\lambda_{\rm los}$ from the ionized
regions so as to more closely mimic the observational context of the
measurement of mfp, since such measurements cannot be performed in still
neutral medium. To this end, for every redshift we divide the distribution of
$\rm\lambda_{\rm los}$ in two modes, at the approximate position of the
minimum of the distributions and therefore the transition between the modes,
and consider only the long mode of the distribution, corresponding to ionized
LoSs. We refer to this approach as the ionised method. We use the notation
$\rm\lambda_{\rm i}$ to refer to mfps computed in this manner.

##  3 Results

###  3.1 Global properties of the IGM during the EoR

The top left panel of Fig. 1 shows the evolution of the average volume
weighted neutral fraction of hydrogen ($\rm\langle x_{HI}\rangle$) in CoDa III
and CoDa II, along with a range of constraints. In both simulations,
$\rm\langle x_{HI}\rangle$ starts to drop noticeably between $\rm z=10$ and
$\rm z=9$. However, in CoDa III, reionization ends later and less abruptly
than in CoDa II, only finishing by $\rm z\approx 5.6$ ($\rm\approx 6$ in CoDa
II). Further, in CoDa II the final average value of the neutral fraction is
far lower than in CoDa III. This latter difference is due to the combination
of the new star formation sub-grid model and improved resolution, as explained
in Oc21. Considering this, along with the rather late reionization obtained
(compared to the canonical value of z=6 for the end of reionization in use
since Fan et al., 2006), CoDa III is in much better agreement with the
constraints from the transmission in the $\rm Lyman\mathchar 45\relax\alpha$
forest (and modelling) on $\rm\langle x_{HI}\rangle$ after reionization. In
particular, the agreement with the data from Bo22 is excellent. The match with
constraints from Becker et al. (2015) is also remarkable in CoDa III after the
ankle of the $\rm\langle x_{HI}\rangle$ curve, although Becker et al. (2015)
favours slightly earlier reionization than in CoDa III. Finally we find that
the older constraints from Fan et al. (2006) are disfavoured both by the more
recent determinations of Bo22 and Becker et al. (2015), and our own
simulations. The bottom left panel of Fig. 1 shows the evolution of the
average ionised fraction of Hydrogen ($\rm x_{HII}$) on a linear scale, along
with the same constraints. CoDa III is in better agreement with a greater
fraction of the collection of constraints for $\rm x_{HII}\gtrsim 0.1$ than
CoDa II. Overall the agreement of CoDa III with modern constraints on the
evolution of the progress of reionization is very good.

Figure 1: _Top left_ : Volume weighted average neutral fraction along with a
collection of constraints including Fan et al. (2006); Becker et al. (2015);
Bosman et al. (2022); Ouchi et al. (2010); McGreer et al. (2015); Totani et
al. (2016); Wang et al. (2021); D̆urovŭíková et al. (2020); Bañados et al.
(2018); Schroeder et al. (2013); Mortlock et al. (2011); Tilvi et al. (2014);
Pentericci et al. (2014); Schenker et al. (2014); Ono et al. (2012); Mason et
al. (2018); Hoag et al. (2019); Greig & Mesinger (2017); Jung et al. (2020).
_Top right_ : Ionisation rate in ionized regions when compared to constraints
from Calverley et al. (2011); Becker & Bolton (2013); D’Aloisio et al. (2018);
Becker et al. (2021). _Bottom Left_ : Volume weighted ionised fraction.
_Bottom Right_ : Optical depth to CMB photons due to scattering by free
electrons, compared to the Planck collaboration results Planck18. CoDa III is
in very good agreement with constraints on the history of ionisation of the
IGM as well as on the ionization rate, whilst remaining compatible with the
optical depth due to electron scattering measured by Planck, despite a rather
late end of reionization. Overall the improvement over CoDa II is very
significant.

The top right panel of Fig. 1 shows the evolution of the ionisation rate in
ionized regions ($\rm\langle\Gamma_{HII}\rangle$) in both CoDa III and CoDa
II. In the former, $\rm\langle\Gamma_{HII}\rangle$ rises rapidly by roughly a
factor of 10 between $\rm z=7$ and $\rm z=5.5$, after which the slope becomes
shallower. CoDa III ’s $\rm\langle\Gamma_{HII}\rangle$ is a very good fit to
the mix of observational constraints available, and a great improvement on the
$\rm\langle\Gamma_{HII}\rangle$ from CoDa II, which overshoots all constraints
for the presented redshift range. This excess is symptomatic of CoDa II ’s
excessive ionisation of the IGM (see ; , for discussions about this). As with
$\rm\langle x_{HI}\rangle$, the origin of this improvement lies both in the
switch to a late reionization calibration and in the improved star formation
sub-grid model, as well as increased resolution (as discussed in Oc21, ).

The bottom right panel of Fig. 1 shows that although reionization ends later
in CoDa III, the optical depth due to scattering on free electrons seen by the
Cosmic Microwave Background remains compatible with the constraints given by
the latest Planck results ( Planck18, ), demonstrating that a late
reionization scenario may still yield a high enough electron scattering
optical depth666An even higher optical depth may still be achieved by adding
the contribution from mini-halos at $\rm z\geq 15$, adding an extra
$\rm\Delta\tau_{\rm z\geq 15}\approx 0.008$, that are not treated in this work
but suggested by Ahn & Shapiro (2021)..

Figure 2: Average evolution of the gas temperature in CoDa III (blue) and CoDa
II (orange). For comparison, constraints from Boera et al. (2019); Walther et
al. (2019); Gaikwad et al. (2020) are also shown.

Fig. 2 shows the average evolution of the gas temperature in CoDa III (blue)
and CoDa II (orange). In both simulations, the average temperature rises as
Reionization progresses and a larger and larger fraction of the IGM is photo-
heated. Eventually, the temperatures reach a characteristic peak as Hydrogen
Reionization completes, before cooling slightly. We find that the average
temperature in CoDa III (and to some extent in CoDa II) is consistent with
recent constraints when available ($\rm 6\geq z\geq 4.6$). CoDa II heats up
much more rapidly, which is coherent with the earlier Reionization and higher
photo-ionisation rates shown in Fig. 1.

Overall reionization in CoDa III reaches a better agreement with observational
constraints on the IGM and other works from the literature than CoDa II did,
signalling a marked improvement in our description of reionization.

###  3.2 Results: the mean free path of ionising photons

We now move to our titular subject: the mean free path of ionising photons.
Fig 3 shows the evolution of mfp over time in CoDa III. Observations from
Worseck et al. (2014) and Be21 show that the mfp increases over time for $\rm
z<6$. The full blue line with errorbars shows that in CoDa III,
$\rm\lambda_{\rm AT}$ increases from $\approx\,$$1\text{\,}\mathrm{pMpc}$ at
$\rm z\approx 6$ to $\approx\,$$20\text{\,}\mathrm{pMpc}$ at $\rm z\approx
4.6$. In fact, $\rm\lambda_{\rm AT}$ grows faster between $\rm z\approx 5.75$
and $\rm z\approx 5.2$, allowing CoDa III to match the newest low mfp from
observations at $\rm z=6$ (Becker et al., 2021) and the high observed mfp at
$\rm z\lesssim 5.25$, whereas most other theoretical work has tended to
overshoot the observed mfp at z=6. Because of the stronger evolution of the
mfp in CoDa III, we report values that are significantly lower than in other
theoretical work for $\rm z\gtrsim 6$, where observations are very
challenging. This points to substantially different progressions of
reionization in the IGM of these other simulations, feasibly implying
differences in ionising emissivities and source populations.

Already, previous work (; ) has highlighted that the latest determinations of
mfp by Be21 point to a late and rapid reionization, which CoDa III matches.
However, unlike in Cain et al. (2021) we do not need to invoke unresolved
photon sinks in order to obtain our agreement with the observed mfp.

The THESAN simulation (Garaldi et al., 2021) also finds a good agreement with
the $\rm z=6$ data point from Be21, however, unfortunately, the simulation
ended too soon to comment on their agreement with lower redshift constraints,
and it is therefore difficult to gauge the overall success of their simulation
in reproducing the evolution of the properties of the IGM over the whole EoR,
down to the post-overlap regime. In any case, it is possible that fully
coupled RHD simulations are better able to match the mfp constraint from Be21
at $\rm z=6$ than post-processing or modelling methods, perhaps because of the
more consistent treatment of the radiative-hydrodynamical interactions in the
IGM and possible radiative feedback in galaxies, and/or the superior
resolution.

Figure 3: Evolution of the mean free path with redshift in CoDa III, and
comparison to observational constraints from (Worseck et al., 2014; Becker et
al., 2021) and previous simulation results from the literature (Keating et
al., 2020; D’Aloisio et al., 2020; Cain et al., 2021; Garaldi et al., 2021).
The blue graphics show results from CoDa III. The full blue line with circles
and errorbars show the median and $\rm 25^{th}-75^{th}$ percentile region for
$\rm\lambda_{\rm AT}$. The full blue line with triangles shows the mean
$\rm\lambda_{\rm i}$, the blue dashed line and blue region show the median and
$\rm 25^{th}-75^{th}$ percentile region $\rm\lambda_{\rm i}$. The median mean
free path in CoDa III is in good agreement with existing observational
constraints from $\rm z=4.5$ till $\rm z=6$.

In order to gain further insight into what drives the mfp measured in CoDa
III, we show the distribution of individual LoS mfp (i.e. $\rm\lambda_{\rm
los}$) at several redshifts in Fig. 4. At $\rm z=12.9$, before reionization is
fully underway, the distribution of $\rm\lambda_{\rm los}$ has a single mode
near ${10}^{-4}\text{\,}\mathrm{pMpc}$. Although this is shorter than the
simulation’s cell size, it is physical: indeed, the general formula for the
mean free path of ionising photons reads $1/\sigma\rho_{\rm HI}$, where
$\sigma$ is the interaction cross-section and $\rho_{\rm HI}$ the number
density of neutral hydrogen atoms. In the dense, neutral high-redshift
Universe before reionization, mean free paths shorter than our cell size
naturally occur and are fully consistent with the measurements presented here.

By $\rm z=7.3$, reionization is in full swing, and the distribution of
$\rm\lambda_{\rm los}$ is now bi-modal, displaying a short and a long mode.
The short mode peaks around $8\text{\times}{10}^{-4}\text{\,}\mathrm{pMpc}$
and seems to descend from the $\rm z=12.9$ mode. However, approximately 30% of
LoS are now scattered around a secondary, longer mode peaking near
$\rm\lambda_{\rm los}$=${10}^{-1}\text{\,}\mathrm{pMpc}$. This bi-modality
reflects the co-existence of large fully ionized regions already by $\rm
z=7.3$, when reionization is almost complete in half of CoDa III ’s volume
(about 40% ionized reading from the bottom left panel of Fig. 1). At $\rm
z=6.1$, both modes have moved to slightly higher $\rm\lambda_{\rm los}$ values
(${10}^{-3}\text{\,}\mathrm{pMpc}$ for the short $\rm\lambda_{\rm los}$ peak,
and $0.64\text{\,}\mathrm{pMpc}$ for the long $\rm\lambda_{\rm los}$ peak),
and the relative weight between modes has swapped: now most LoS belong to the
long $\rm\lambda_{\rm los}$ peak.

By $\rm z=5.2$, reionization is complete, all of the LoS have been ionized and
the short mode has disappeared in favour of the long mode. The evolution of
the $\rm\lambda_{\rm los}$ distribution reveals that the CoDa III
$\rm\lambda_{\rm los}$ evolves due to both gradual changes along every LoS
(due to cosmic expansion driving the average physical density to decrease over
time), and to rapid changes along short $\rm\lambda_{\rm los}$ LoS due to
reionization. The rapid increase of $\rm\lambda_{\rm AT}$ in Fig. 3 between
$\rm z=6.1$ and $\rm z=5.2$ could be driven by the disappearance of the
shortest $\rm\lambda_{\rm los}$ values as more and more LoS experience
reionization.

Figure 4: Distribution of $\rm\lambda_{\rm los}$. Dashed vertical lines show
the distribution mean for each redshift. For readability, only a subset of
redshifts roughly 200 Myr apart between z=12.9 and z=5.2 is shown. The
$\rm\lambda_{\rm los}$ are distributed in a neutral opaque mode of short mean
free paths, and an ionised transparent mode of long mean free paths. The
relative weight of these modes changes as reionization progresses.

Fig. 4 shows that the transition between the ionised mode of the
$\rm\lambda_{\rm los}$ distribution and the neutral mode occurs near
${10}^{-2}\text{\,}\mathrm{pMpc}$). Accordingly we select only those
$\rm\lambda_{\rm los}$ longer than this threshold for the ionised method. The
blue solid curve with triangles in Fig. 3 shows the evolution of
$\rm\lambda_{\rm i}$ in CoDa III. For $\rm z>5$, both the mean and the median
$\rm\lambda_{\rm i}$ are lower than $\rm\lambda_{\rm AT}$. This makes our
prediction based on the mean $\rm\lambda_{\rm i}$ a better match to the $\rm
z=6$ constraint from Be21, and also shows that the average transmission method
(measuring the mean free path from an average simulated transmission spectrum)
is biased towards high transmissions and overestimates the median
$\rm\lambda_{\rm los}$ at the redshifts considered in Fig. 3. For $\rm z\leq
5$, the mean $\rm\lambda_{\rm i}$ is slightly above the constraints from
Worseck et al. (2014), however the median remains in very good agreement and
there is a broad scatter, comparable to the error bar of Be21.

##  4 Conclusions

In this article, we have presented first results from our latest large scale
simulation of the Epoch of Reionization, Cosmic Dawn III (CoDa III ), the
largest fully coupled RHD simulation of Reionization ever performed. It
improves on our previous CoDa II simulation in many aspects, such as using a
two times higher spatial resolution and 8 times higher mass resolution, as
well as updates of the physical modules and sub-grid models, in particular
star formation and stellar source modelling. We show that, thanks to these
improvements, CoDa III compares much better than CoDa II with observational
constraints of the high-redshift Universe, in particular the evolution of the
neutral hydrogen fraction which sees a dramatic improvement, as well as the
ionizing rate in ionized regions, while the evolution of the ionized fraction
and the electron scattering optical depth are also in good agreement with the
constraints available. Moreover, we showed that CoDa III is the first fully
coupled radiative hydrodynamics simulation to reproduce the latest constraints
on the ionising photon mean free path all the way from $\rm z=6$ down to $\rm
z=4.6$ thanks to a late reionization ($\rm z_{rei}=5.6$) and a rapid evolution
of the ionising mean free path. In particular, CoDa III shows that the
surprisingly low mean free path measured by Be21 can naturally occur in a
well-calibrated simulation of the EoR such as ours. The distribution of mean
free paths during reionization is generally bi-modal, displaying a short,
neutral mode and a long, ionized mode. Both modes shift to higher mfp as
reionization progresses: The neutral mode peak shifts from sub-kpc to kpc
scales of mfp, and the ionized mode peak from 0.1 Mpc/h (at z=7) to $\sim 10$
Mpc/h at ($\rm z=5.2$). The balance of these modes reflects the degree of
ionisation of the simulated volume. The bi-modal distribution reduces to a
single short (long) mode before (after) reionization, respectively. Finally,
by using 3 different methods for computing the mean free path from CoDa III,
we show that the choice of the method is important. In particular, measuring
the mean free path from an average transmitted flux stacking a number of lines
of sight tends to overestimate the mean free path otherwise obtained as the
average or the median of the mfp of individual lines of sight. Overall, CoDa
III’s description of reionization in the IGM appears to be a marked
improvement over CoDa II, providing an adequate basis with which to perform a
broad range of studies of the EoR, such as $\rm Lyman\mathchar 45\relax\alpha$
transmission in the IGM, but also of star formation suppression, and the
ionising photons budget of galaxies.

## Acknowledgements

The authors thank F. Davies and H. Katz for fruitful discussions. Cosmic Dawn
III was performed on Summit at Oak Ridge National Laboratory / Oak Ridge
Leadership Computing Facility (Project AST031) using an INCITE 2020
allocation. JL acknowledges support from the DFG via the Heidelberg Cluster of
Excellence STRUCTURES in the framework of Germany’s Excellence Strategy (grant
EXC-2181/1 - 390900948). JS thanks Sergey Pilipenko for sharing the latest
version of the ginnungagap code as well as for his guidance with using it. JS
acknowledges support from the ANR LOCALIZATION project, grant ANR-21-CE31-0019
of the French Agence Nationale de la Recherche. JS gratefully acknowledges the
Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for providing
computing time to build the initial conditions of the simulations for this
project on the GCS Supercomputer SuperMUC-NG at Leibniz Supercomputing Centre
(www.lrz.de). GY acknowledges financial support from MICIU/FEDER under project
grant PGC2018-094975-C21. KA was supported by NRF-2016R1D1A1B04935414,
2021R1A2C1095136, and 2016R1A5A1013277. ITI was supported by the Science and
Technology Facilities Council [grant numbers ST/I000976/1 and ST/T000473/1]
and the Southeast Physics Network (SEP-Net). HP was supported by the World
Premier International Research Center Initiative (WPI), MEXT, Japan and JSPS
KAKENHI grant No. 19K23455.

## Data availability

The authors will share the data used for this article upon reasonable request.

## References

  * Ahn & Shapiro (2021) Ahn K., Shapiro P. R., 2021, arXiv:2011.03582 [astro-ph]
  * Aubert & Teyssier (2008) Aubert D., Teyssier R., 2008, Mon Not R Astron Soc, 387, 295
  * Aubert & Teyssier (2010) Aubert D., Teyssier R., 2010, The Astrophysical Journal, 724, 244
  * Bañados et al. (2018) Bañados E., et al., 2018, Nature, 553, 473
  * Becker & Bolton (2013) Becker G. D., Bolton J. S., 2013, Monthly Notices of the Royal Astronomical Society, 436, 1023
  * Becker et al. (2015) Becker G. D., Bolton J. S., Madau P., Pettini M., Ryan-Weber E. V., Venemans B. P., 2015, Monthly Notices of the Royal Astronomical Society, 447, 3402
  * Becker et al. (2021) Becker G. D., D’Aloisio A., Christenson H. M., Zhu Y., Worseck G., Bolton J. S., 2021, arXiv e-prints, 2103, arXiv:2103.16610
  * Boera et al. (2019) Boera E., Becker G. D., Bolton J. S., Nasir F., 2019, The Astrophysical Journal, 872, 101
  * Bosman et al. (2022) Bosman S. E. I., et al., 2022, MNRAS, 514, 55
  * Cain et al. (2021) Cain C., D’Aloisio A., Gangolli N., Becker G. D., 2021, arXiv e-prints, 2105, arXiv:2105.10511
  * Calverley et al. (2011) Calverley A. P., Becker G. D., Haehnelt M. G., Bolton J. S., 2011, Monthly Notices of the Royal Astronomical Society, 412, 2543
  * D’Aloisio et al. (2018) D’Aloisio A., McQuinn M., Davies F. B., Furlanetto S. R., 2018, Monthly Notices of the Royal Astronomical Society, 473, 560
  * D’Aloisio et al. (2020) D’Aloisio A., McQuinn M., Trac H., Cain C., Mesinger A., 2020, arXiv e-prints, 2002, arXiv:2002.02467
  * Davies et al. (2021) Davies F. B., Bosman S. E. I., Furlanetto S. R., Becker G. D., D’Aloisio A., 2021, arXiv e-prints, p. arXiv:2105.10518
  * Dawoodbhoy et al. (2018) Dawoodbhoy T., et al., 2018, MNRAS, 480, 1740
  * Deparis et al. (2019) Deparis N., Aubert D., Ocvirk P., Chardin J., Lewis J., 2019, Astronomy and Astrophysics, 622, A142
  * Dixon et al. (2018) Dixon K. L., Iliev I. T., Gottlöber S., Yepes G., Knebe A., Libeskind N., Hoffman Y., 2018, Monthly Notices of the Royal Astronomical Society, 477, 867
  * Dubois & Teyssier (2008) Dubois Y., Teyssier R., 2008, Astronomy and Astrophysics, 477, 79
  * D̆urovŭíková et al. (2020) D̆urovŭíková D., Katz H., Bosman S. E. I., Davies F. B., Devriendt J., Slyz A., 2020, Monthly Notices of the Royal Astronomical Society, 493, 4256
  * Eldridge & Stanway (2020) Eldridge J. J., Stanway E. R., 2020, arXiv e-prints, 2005, arXiv:2005.11883
  * Fan et al. (2006) Fan X., et al., 2006, The Astronomical Journal, 132, 117
  * Gaikwad et al. (2020) Gaikwad P., et al., 2020, Monthly Notices of the Royal Astronomical Society, 494, 5091
  * Garaldi et al. (2021) Garaldi E., Kannan R., Smith A., Springel V., Pakmor R., Vogelsberger M., Hernquist L., 2021, arXiv e-prints, p. arXiv:2110.01628
  * Gnedin (2016) Gnedin N. Y., 2016, The Astrophysical Journal, 833, 66
  * Greig & Mesinger (2017) Greig B., Mesinger A., 2017, MNRAS, 465, 4838
  * Hoag et al. (2019) Hoag A., et al., 2019, The Astrophysical Journal, 878, 12
  * Iliev et al. (2011) Iliev I. T., Moore B., Gottlöber S., Yepes G., Hoffman Y., Mellema G., 2011, Monthly Notices of the Royal Astronomical Society, 413, 2093
  * Jung et al. (2020) Jung I., et al., 2020, ApJ, 904, 144
  * Kannan et al. (2019) Kannan R., Vogelsberger M., Marinacci F., McKinnon R., Pakmor R., Springel V., 2019, MNRAS, 485, 117
  * Kannan et al. (2021) Kannan R., Garaldi E., Smith A., Pakmor R., Springel V., Vogelsberger M., Hernquist L., 2021, arXiv e-prints, p. arXiv:2110.00584
  * Keating et al. (2020) Keating L. C., Weinberger L. H., Kulkarni G., Haehnelt M. G., Chardin J., Aubert D., 2020, Monthly Notices of the Royal Astronomical Society, 491, 1736
  * Kulkarni et al. (2016) Kulkarni G., Choudhury T. R., Puchwein E., Haehnelt M. G., 2016, Mon. Not. R. Astron. Soc., 463, 2583
  * Kulkarni et al. (2019) Kulkarni G., Keating L. C., Haehnelt M. G., Bosman S. E. I., Puchwein E., Chardin J., Aubert D., 2019, Monthly Notices of the Royal Astronomical Society: Letters, 485, L24
  * Levermore (1984) Levermore C. D., 1984, Journal of Quantitative Spectroscopy and Radiative Transfer, 31, 149
  * Lewis et al. (2022) Lewis J. S. W., Ocvirk P., Dubois Y., Aubert D., Chardin J., Gillet N., Thélie E., 2022
  * Mason et al. (2018) Mason C. A., Treu T., Dijkstra M., Mesinger A., Trenti M., Pentericci L., de Barros S., Vanzella E., 2018, The Astrophysical Journal, 856, 2
  * McGreer et al. (2015) McGreer I. D., Mesinger A., D’Odorico V., 2015, Monthly Notices of the Royal Astronomical Society, 447, 499
  * Mortlock et al. (2011) Mortlock D. J., et al., 2011, Nature, 474, 616
  * Ocvirk & Aubert (2011) Ocvirk P., Aubert D., 2011, MNRAS, 417, L93
  * Ocvirk et al. (2013) Ocvirk P., Aubert D., Chardin J., Knebe A., Libeskind N., Gottlöber S., Yepes G., Hoffman Y., 2013, ApJ, 777, 51
  * Ocvirk et al. (2014) Ocvirk P., et al., 2014, ApJ, 794, 20
  * Ocvirk et al. (2016) Ocvirk P., et al., 2016, Monthly Notices of the Royal Astronomical Society, 463, 1462
  * Ocvirk et al. (2019) Ocvirk P., Aubert D., Chardin J., Deparis N., Lewis J., 2019, Astronomy and Astrophysics, 626, A77
  * Ocvirk et al. (2020) Ocvirk P., et al., 2020, Monthly Notices of the Royal Astronomical Society
  * Ocvirk et al. (2021) Ocvirk P., Lewis J. S. W., Gillet N., Chardin J., Aubert D., Deparis N., Thelie E., 2021, arXiv e-prints, 2105, arXiv:2105.01663
  * Ono et al. (2012) Ono Y., et al., 2012, The Astrophysical Journal, 744, 83
  * Ouchi et al. (2010) Ouchi M., et al., 2010, The Astrophysical Journal, 723, 869
  * Pentericci et al. (2014) Pentericci L., et al., 2014, The Astrophysical Journal, 793, 113
  * Planck Collaboration et al. (2014) Planck Collaboration et al., 2014, Astronomy & Astrophysics, Volume 571, id.A16, <66> pp., 571, A16
  * Rahmati & Schaye (2018) Rahmati A., Schaye J., 2018, Monthly Notices of the Royal Astronomical Society, 478, 5123
  * Schenker et al. (2014) Schenker M. A., Ellis R. S., Konidaris N. P., Stark D. P., 2014, The Astrophysical Journal, 795, 20
  * Schroeder et al. (2013) Schroeder J., Mesinger A., Haiman Z., 2013, Monthly Notices of the Royal Astronomical Society, 428, 3058
  * Sorce & Tempel (2018) Sorce J. G., Tempel E., 2018, Monthly Notices of the Royal Astronomical Society, 476, 4362
  * Sorce et al. (2016) Sorce J. G., et al., 2016, Monthly Notices of the Royal Astronomical Society, 455, 2078
  * Teyssier (2002) Teyssier R., 2002, A&A, 385, 337
  * The Planck Collaboration (2020) The Planck Collaboration 2020, arXiv:2007.04997 [astro-ph]
  * Tilvi et al. (2014) Tilvi V., et al., 2014, The Astrophysical Journal, 794, 5
  * Totani et al. (2016) Totani T., Aoki K., Hattori T., Kawai N., 2016, Publications of the Astronomical Society of Japan, 68, 15
  * Walther et al. (2019) Walther M., Oñorbe J., Hennawi J. F., Lukić Z., 2019, The Astrophysical Journal, 872, 13
  * Wang et al. (2021) Wang B., et al., 2021, arXiv e-prints, 2104, arXiv:2104.03432
  * Worseck et al. (2014) Worseck G., et al., 2014, Monthly Notices of the Royal Astronomical Society, 445, 1745

##  Appendix A Consistency check of mean free paths for individual lines of
sight

Figure 5: _Left:_ Distribution of $\rm\lambda_{\rm los}$. Dashed vertical
lines show the distribution mean for each redshift. For readability, only a
subset of redshifts roughly 200 Myr apart between z=12.9 and z=5.2 is shown.
The $\rm\lambda_{\rm los}$ are distributed in a neutral opaque mode of short
mean free paths, and an ionised transparent mode of long mean free paths. The
relative weight of these modes changes as reionization progresses. _Right:_ We
compute the evolution of the mean free path along ionised lines of sight, by
applying the same methodology to the intercept distribution as used to compute
the $\rm\lambda_{\rm los}$. The dotted lines show the result when compared to
the other techniques we employ.

In order to verify the $\rm\lambda_{\rm los}$ and $\rm\lambda_{\rm i}$ values
for individual lines of sight, we also compute these using another method
found in the literature (e.g. Rahmati & Schaye, 2018). Using Eq. 1, one finds
that:

$\begin{array}[]{l}\rm\frac{F(x)}{F_{0}}=exp(-1),\,if\,x=\lambda\end{array}$
(2)

Applying this to each LoS, we search for the cells whose transmission $\rm
T_{a}(x_{a})=F_{a}(x_{a})/{F_{0}}$ and $\rm T_{b}(x_{b})$ bracket $\rm
exp(-1)$. We then interpolate between the cells to find x so that $\rm
T(x)=exp(-1)$. There are two limit cases that we must deal with. First, we
must treat the case of very transparent LoS for which $\lambda$ is longer than
CoDa III’s box length. In order to estimate the mean free path along these
lines of sight, we equate their high transmissions with the absence of large
overdensities and assume that the density of neutral Hydrogen is constant and
small along the LoS. Thus, we can obtain an estimate of the full mean free
path by extrapolating the additional distance necessary for the integrated
optical depth to reach one. Second, there are the cases of very opaque LoS (or
at least starting in an opaque cell) where $\lambda$ is smaller than the cell
size of CoDa III. Here, assuming that this cell does not contain an unresolved
absorber, we use the cell’s optical depth to deduce $\lambda$.

The left panel of Fig. 5 shows the distributions for both the fitting method
presented in 2.2, and the method described above (that we call intercept
method from hereon). Both methods are very consistent for short, long, and
intermediate mean free paths, and throughout Reionization in CoDa III.
However, there are some discrepancies, most notably for the shortest and
longest mean free paths. This isn’t too surprising as the exponential fit is
harder in these cases, and at the same time we must make some strong
assumptions for the intercept method. The right panel of Fig. 5, shows the
average evolution of the ionised lines of sight as computed with all tested
techniques, including the intercept method. As expected based on the
distributions, the result is very similar to the evolution of the mean free
path as derived using fitting. In fact, for $\rm z\leq 6$, the two methods are
almost indistinguishable. Though at higher redshifts, the intercept method
seems to predict a slightly higher mean free path. That these two methods
independently yield such similar mean free paths, shows the robustness of our
technique and results.


# A multistage framework for studying the evolution of jets and high-$p_{T}$
probes in small collision systems

Abhijit<EMAIL_ADDRESS>for The JETSCAPE Collaboration,
Department of Physics and Astronomy, Wayne State University,
666 W. Hancock St., Detroit, Michigan, USA

###### Abstract

Understanding the modification of jets and high-$p_{T}$ probes in small
systems requires the integration of soft and hard physics. We present recent
developments in extending the JETSCAPE framework to build an event generator,
which includes correlations between soft and hard partons, to study jet
observables in small systems. The multi-scale physics of the collision is
separated into different stages. Hard scatterings are first sampled at binary
collision positions provided by the Glauber geometry. They are then propagated
backward in space-time following an initial-state shower to obtain the
initiating partons’ energies and momenta before the collision. These energies
and momenta are then subtracted from the incoming colliding nucleons for soft-
particle production, modeled by the 3D-Glauber + hydrodynamics + hadronic
transport framework. This new hybrid approach (X-SCAPE) includes non-trivial
correlations between jet and soft particle productions in small systems. We
calibrate this framework with the final state hadron’s $p_{T}$-spectra from
low to high $p_{T}$ in $p$-$p$, and and then compare with the spectra in
$p$-$Pb$ collisions from the LHC. We also present results for additional
observables such as the distributions of event activity as a function of the
hardest jet $p_{T}$ in forward and mid-rapidity for both $p$-$p$ and $p$-$Pb$
collisions.

##  1 Introduction

Over the last decade there has been extensive interest in the possibility that
a quark-gluon plasma (QGP) could have been produced in small collision systems
such as $p$-$p$ and $p$-$A$, especially in events with very high multiplicity
at top RHIC and LHC energies. The first indication of this was the observation
of a ridge in two particle correlations over long ranges in rapidity in both
$p$-$p$ and $p$-$A$ collisions
ridge,*ALICE:2012eyl,*ATLAS:2012cix,*CMS:2010ifv,*CMS:2012qk,*PHENIX:2013ktj,*PHENIX:2014fnc,*STAR:2015kak.
The observation of collective flow in $p$-$Pb$ collisions at the LHC
v2,*ATLAS:2013jmi,*ALICE:2014dwt,*CMS:2015yux, reinforced via multi-particle
correlations, along with successful theory comparisons, strongly suggests the
formation of expanding droplets of QGP in these collisions.

The _standard model_ of heavy-ion collisions
reviews,*Muller:2012zq,*Jacobs:2004qv,*Gyulassy:2007kla,*Arslandok:2023utm
includes the modification of hard jets as a signature of the QGP. This raises
the question of the effect on jets (or hard processes in general) from these
small droplets of QGP formed in $p$-$p$ and $p$-$Pb$ collisions. A study of
the nuclear modification factor $R_{pA}$ for both jets and leading hadrons, in
minimum bias events, contains no observable modification, i.e., $R_{pA}\simeq
1$ [21, 22]. However, if events with jets are binned in centrality, different
measures of centrality yield differing variations of $R_{pA}$ with centrality,
with some measures yielding a large suppression $R_{pA}\lesssim 1$ for
nominally central events, coupled with an enhancement $R_{pA}\gtrsim 1$ for
peripheral events [21]. Along with this, one also observes a residual
azimuthal anisotropy at high $p_{T}$, even above $20$ GeV, indicative of jet
modification highptv2,*ATLAS:2014qaj,*ATLAS:2019vcm.

A possible reason for these conflicting results could simply be consequences
of small system size and the much tighter constraints of energy momentum
conservation in the initial stages of such systems. In the typical heavy-ion
simulation
jetscapejets,*JETSCAPE:2021ehl,*JETSCAPE:2022jer,*JETSCAPE:2022hcb,*JETSCAPE:2023hqn,
one simulates the bulk medium evolution first
jetscapebulk,*JETSCAPE:2020mzn,*JETSCAPE:2020shq. Given the fluctuating
configuration of the initial state, a certain amount of energy and momentum is
deposited over a range of space-time, and then evolved using a viscous fluid
dynamical simulation. After this, one introduces hard jets whose location
depends on the binary collision profile and these then propagate and interact
with the pre-calculated fluid medium. However, no attempt is made to balance
the energy-momentum required for the production of the jet with the energy
momentum deposited towards the production of the bulk medium, as the energy
within a jet is a tiny fraction of the energy within the bulk medium. This is
no longer the case for small collision systems [37].

In spite of the flexibility of the JETSCAPE framework jetscape-
framework,*JETSCAPE:2017eso,*Putschke:2019yrg, it cannot be used to accurately
simulate the collisions of small systems where exact initial state
conservation of energy-momentum may be required, i.e., the framework does not
provide tools to carry this out easily. This situation is ameliorated by the
new X-SCAPE framework, which provides various infrastructure enhancements to
aid in the incorporation of exact energy-momentum conservation and other
initial state effects required in small systems. In Sec. 2, we will describe
the salient features of the X-SCAPE project. Preliminary results using the new
coupled modules of i-MATTER for the hard sector and 3D MC-Glauber for the soft
sector are described in Sec. 3. An outlook is presented in Sec. 4.

##  2 The X-SCAPE framework

The X-ion collisions with a Statistically and Computationally Advanced Program
Envelope (X-SCAPE) is the second project of the JETSCAPE collaboration [42].
The main goals of the new framework are to create a JETSCAPE like framework,
which is fully backwards compatible, can be extended to encompass the physics
of small collision systems, low energy heavy-ion collisions relevant to the
beam energy scan, and make the first extension to electron-ion collisions.

In Fig. 1, we present the flow chart of the X-SCAPE framework. Two important
differences with the older JETSCAPE framework should immediately become
obvious. The input systems are not limited to be $p$-$p$ or $A$-$A$ and can
accommodate any system. Also, there is an intentional correlation between the
hard scattering, soft stopping and energy deposition modules, connected to the
initial state (leftmost red rectangle in Fig. 1). These will require new
modules for both the soft energy deposition and the hard scattering, which can
correlate the energy-momentum being deposited in these systems, described in
the next section.

Figure 1: The flow chart of the X-SCAPE framework

The other new aspects of the new framework are the concurrent running of the
hadronic cascade (bottom yellow arrow in Fig. 1) in tandem with the fluid
dynamical simulation (middle green arrow). This aspect is required for the
future simulation of lower energy heavy-ion collisions. The other new aspect
is the study of the spectator (top red arrow), which is a required aspect of
future simulations of electron ion collisions. In the remainder of these
proceedings, we will focus on the simulations of small systems using this new
framework.

##  3 Results with the new modules: i-MATTER and 3D MC-Glauber

In order to simulate the physics of jet and bulk correlations in $p$-$p$ and
$p$-$A$ collisions, we introduce two new modules: 3D MC-Glauber, which
generates the bulk or soft energy and Baryon number deposition, setting up the
initial state for the subsequent fluid simulation, and i-MATTER which
generates the configuration of hard partons prior to the hard scattering.

Figure 2: Particle production from the bulk simulation with and without a hard
scattering in the event.

Using the X-SCAPE framework, we have introduced a new methodology for
simulating small systems in $p$-$Pb$ and $p$-$p$ collisions: The soft sector
of all colliding systems, ranging from low multiplicity $p$-$p$ to the highest
multiplicity $p$-$A$ (and obviously extending to all $A$-$A$) collisions will
be simulated by a viscous fluid dynamical simulation. While fluid dynamics
represents actual physics in dense systems such as $A$-$A$ or high
multiplicity $p$-$A$, its custodial energy-momentum conservation conditions
should provide a reliable model for a majority of soft observables in $p$-$p$.

Figure 3: Left: Charged pion spectra from $p$-$p$ and $p$-$A$ collisions at
$\sqrt{s}=5.02$ TeV as measured by the ALICE experiment compared with
preliminary calculations using i-MATTER, PYTHIA and MATTER (final state),
along with 3D MC-Glauber and MUSIC within the X-SCAPE framework. Right: Jet
spectra measured in the same experiment and compared to the output from the
hard sector from X-SCAPE.

The above paradigm successfully reproduces a variety of observables in $p$-$p$
and $p$-$A$ collisions, when combined with a parametrized initial state such
as 3D MC-Glauber, as used in this work. The 3D MC-Glauber model
MCG,*Shen:2017bsr,*Shen:2022oyg simulates each incoming proton in terms of
several (typically 4) hot-spots. Interactions between these opposing hot
spots, calculated using sub-nucleon cross sections generates flux tube-like
strings whose spatial extent simulates the deposited energy and baryon number.
These energy and baryon number distributions become inputs to the subsequent
fluid dynamical simulation simulated using the MUSIC generator. Hadronization
of the soft sector is carried out using an event-by-event Cooper-Frye routine
called iSS.

Figure 4: Left: Nuclear modification factor for charged hadrons for minimum
bias $p$-$Pb$ collisions at $\sqrt{s}=5.02$ TeV compared to preliminary
results from the X-SCAPE framework. Right: Nuclear modification factor for
jets in the same experiment compared to resutls from X-SCAPE. Spectra are
presented in Fig. 3.

Hard processes are simulated in 4 steps: Hard scatterings are sampled using
PYTHIA at locations obtained from sampling the binary hot-spot collision
profile from 3D MC-Glauber. Then the initial state radiation is done with
i-MATTER, which maintains exact energy-momentum conservation in the backwards
shower, along with locations of each splits and the correct color tags for all
partons. The initial partons at a large negative time are matched with a
remnant anti-parton with the requisite anti-color so that the entire hard
sector is a color singlet.

Figure 5: The calculated combined (soft + hard process) transverse energy
produced in $p$-$p$ collisions at $\sqrt{s}=5.02$ TeV, with $|\eta|<1$ (left),
and $-4.9<\eta<-3.2$ (right), as a function of the leading jet $p_{T}$.

Both partons and remnants are then conveyed to a second instance of 3D MC-
Glauber where their energy and momentum are subtracted from one of the hot
spots. The depleted nucleons are then allowed to undergo collision leading to
the formation of a bulk medium that lacks the energies of the hard partons
constituting the initial state of the hard scattering. The final hadronization
has two sources, the hadronization of the bulk medium using Cooper-Frye and
the hadronization of the hard sector using the Lund hadronization model. The
effect of the energy-momentum depletion from the soft sector can be seen in
Fig. 2, where only the hadrons from the bulk sector, over the entire rapidity
range, are plotted for different energies within the hard sector, controlled
by $\hat{p}_{T}$, the hardest transverse momentum exchange within the hard
sector.

Figure 6: Same as Fig. 5 above, but for $p$-$Pb$ collisions.

The value of $\hat{p}_{T}$ varies over a range from a minimum
$\hat{p}_{T}^{0}$ to the maximum allowed by kinematics. The value of
$\hat{p}_{T}^{0}$ controls the energy in the hard sector and thus the energy
depletion of the soft sector. We vary this parameter to obtain the best fit
with the minimum bias charged pion spectra in $p$-$p$ collisions. This is
presented in the left panel of Fig. 3. The multiple binary collision enhanced
spectra for $p$-$Pb$ is plotted in the same figure. Jets are obtained by
clustering hadrons from the hard sector and are presented in the right plot of
Fig. 3. Nuclear modification factors for both these spectra are presented in
Fig. 4. This is the first time a consistent description of the $R_{AA}$ that
explains hadron data from $0\lesssim p_{T}\lesssim 20$ GeV, for $p$-$Pb$
collisions has been obtained from one simulation.

With a successful description of the hadron and jet spectrum in $p$-$p$ and
$p$-$Pb$ collsions we ask how much does jet production affect multiplicity or
transverse energy $E_{T}$ production at both mid and forward rapidity. This is
plotted in Fig. 5 for $p$-$p$, and in Fig. 6 for $p$-$Pb$. Both these plots
show qualitatively similar results. While the particle multiplicity at mid-
rapidity increases steadily with the energy of the hardest jet. The
multiplicity first rises and then drops at very forward rapidity where the
effects of energy conservation between hard and soft sectors are clearly
observed.

##  4 Summary

Following these successful preliminary results, the upcoming simulations will
feature a Bayesian calibration of the parameters determining the extent of
hard and soft sectors. Nuclear modification of the incoming parton
distribution functions as well as any possible energy loss of the hard partons
within the droplets of the QGP are currently being explored. Also, the effect
of initial state fluctuations leading to azimuthal anisotropies in both hard
and soft sectors is being studied.

_ACK_ : A.M. is supported in part by the NSF (OAC-2004571) and by the DOE (DE-
SC0013460).

ALICE:2012eyl,ATLAS:2012cix,CMS:2010ifv,CMS:2012qk,PHENIX:2013ktj,PHENIX:2014fnc,STAR:2015kak

ALICE:2012eyl,ATLAS:2012cix,CMS:2010ifv,CMS:2012qk,PHENIX:2013ktj,PHENIX:2014fnc,STAR:2015kak

ATLAS:2013jmi,ALICE:2014dwt,CMS:2015yux

ATLAS:2013jmi,ALICE:2014dwt,CMS:2015yux

Muller:2012zq,Jacobs:2004qv,Gyulassy:2007kla,Arslandok:2023utm

Muller:2012zq,Jacobs:2004qv,Gyulassy:2007kla,Arslandok:2023utm

ATLAS:2014qaj,ATLAS:2019vcm

ATLAS:2014qaj,ATLAS:2019vcm

JETSCAPE:2021ehl,JETSCAPE:2022jer,JETSCAPE:2022hcb,JETSCAPE:2023hqn

JETSCAPE:2021ehl,JETSCAPE:2022jer,JETSCAPE:2022hcb,JETSCAPE:2023hqn

JETSCAPE:2020mzn,JETSCAPE:2020shq

JETSCAPE:2020mzn,JETSCAPE:2020shq

JETSCAPE:2017eso,Putschke:2019yrg

JETSCAPE:2017eso,Putschke:2019yrg

Shen:2017bsr,Shen:2022oyg

Shen:2017bsr,Shen:2022oyg

## References

  * [1]
  * [2] Betty Abelev  “Long-range angular correlations on the near and away side in $p$-Pb collisions at $\sqrt{s_{NN}}=5.02$ TeV”  In _Phys. Lett. B_ 719, 2013, pp. 29–41  DOI: 10.1016/j.physletb.2013.01.012
  * [3] Betty Abelev  “Long-range angular correlations on the near and away side in $p$-Pb collisions at $\sqrt{s_{NN}}=5.02$ TeV”  In _Phys. Lett. B_ 719, 2013, pp. 29–41  DOI: 10.1016/j.physletb.2013.01.012
  * [4] Georges Aad  “Observation of Associated Near-Side and Away-Side Long-Range Correlations in $\sqrt{s_{NN}}$=5.02 TeV Proton-Lead Collisions with the ATLAS Detector”  In _Phys. Rev. Lett._ 110.18, 2013, pp. 182302  DOI: 10.1103/PhysRevLett.110.182302
  * [5] Vardan Khachatryan  “Observation of Long-Range Near-Side Angular Correlations in Proton-Proton Collisions at the LHC”  In _JHEP_ 09, 2010, pp. 091  DOI: 10.1007/JHEP09(2010)091
  * [6] Serguei Chatrchyan  “Observation of Long-Range Near-Side Angular Correlations in Proton-Lead Collisions at the LHC”  In _Phys. Lett. B_ 718, 2013, pp. 795–814  DOI: 10.1016/j.physletb.2012.11.025
  * [7] A. Adare  “Quadrupole Anisotropy in Dihadron Azimuthal Correlations in Central $d+$Au Collisions at $\sqrt{s_{{}_{NN}}}$=200 GeV”  In _Phys. Rev. Lett._ 111.21, 2013, pp. 212301  DOI: 10.1103/PhysRevLett.111.212301
  * [8] A. Adare  “Measurement of long-range angular correlation and quadrupole anisotropy of pions and (anti)protons in central $d+$Au collisions at $\sqrt{s_{{}_{NN}}}$=200 GeV”  In _Phys. Rev. Lett._ 114.19, 2015, pp. 192301  DOI: 10.1103/PhysRevLett.114.192301
  * [9] L. Adamczyk  “Long-range pseudorapidity dihadron correlations in $d$+Au collisions at $\sqrt{s_{\rm NN}}=200$ GeV”  In _Phys. Lett. B_ 747, 2015, pp. 265–271  DOI: 10.1016/j.physletb.2015.05.075
  * [10]
  * [11] Georges Aad  “Measurement with the ATLAS detector of multi-particle azimuthal correlations in p+Pb collisions at $\sqrt{s_{NN}}$ =5.02 TeV”  In _Phys. Lett. B_ 725, 2013, pp. 60–78  DOI: 10.1016/j.physletb.2013.06.057
  * [12] Georges Aad  “Measurement with the ATLAS detector of multi-particle azimuthal correlations in p+Pb collisions at $\sqrt{s_{NN}}$ =5.02 TeV”  In _Phys. Lett. B_ 725, 2013, pp. 60–78  DOI: 10.1016/j.physletb.2013.06.057
  * [13] Betty Bezverkhny Abelev  “Multiparticle azimuthal correlations in p -Pb and Pb-Pb collisions at the CERN Large Hadron Collider”  In _Phys. Rev. C_ 90.5, 2014, pp. 054901  DOI: 10.1103/PhysRevC.90.054901
  * [14] Vardan Khachatryan  “Evidence for Collective Multiparticle Correlations in p-Pb Collisions”  In _Phys. Rev. Lett._ 115.1, 2015, pp. 012301  DOI: 10.1103/PhysRevLett.115.012301
  * [15]
  * [16] Berndt Muller, Jurgen Schukraft and Boleslaw Wyslouch  “First Results from Pb+Pb collisions at the LHC”  In _Ann. Rev. Nucl. Part. Sci._ 62, 2012, pp. 361–386  DOI: 10.1146/annurev-nucl-102711-094910
  * [17] Berndt Muller, Jurgen Schukraft and Boleslaw Wyslouch  “First Results from Pb+Pb collisions at the LHC”  In _Ann. Rev. Nucl. Part. Sci._ 62, 2012, pp. 361–386  DOI: 10.1146/annurev-nucl-102711-094910
  * [18] Peter Jacobs and Xin-Nian Wang  “Matter in extremis: Ultrarelativistic nuclear collisions at RHIC”  In _Prog. Part. Nucl. Phys._ 54, 2005, pp. 443–534  DOI: 10.1016/j.ppnp.2004.09.001
  * [19] Miklos Gyulassy and Larry McLerran  “How and Where to go Within the Standard Model?: sQGP/CGC @ RHIC/LHC”  In _Subnucl. Ser._ 42, 2007, pp. 109–159  DOI: 10.1142/9789812708427˙0005
  * [20] M. Arslandok  “Hot QCD White Paper”, 2023  arXiv:2303.17254 [nucl-ex]
  * [21] Georges Aad  “Centrality and rapidity dependence of inclusive jet production in $\sqrt{s_{\mathrm{NN}}}=5.02$ TeV proton-lead collisions with the ATLAS detector”  In _Phys. Lett. B_ 748, 2015, pp. 392–413  DOI: 10.1016/j.physletb.2015.07.023
  * [22] Vardan Khachatryan  “Charged-particle nuclear modification factors in PbPb and pPb collisions at $\sqrt{s_{\mathrm{N}\;\mathrm{N}}}=5.02$ TeV”  In _JHEP_ 04, 2017, pp. 039  DOI: 10.1007/JHEP04(2017)039
  * [23]
  * [24] Georges Aad  “Measurement of long-range pseudorapidity correlations and azimuthal harmonics in $\sqrt{s_{NN}}=5.02$ TeV proton-lead collisions with the ATLAS detector”  In _Phys. Rev. C_ 90.4, 2014, pp. 044906  DOI: 10.1103/PhysRevC.90.044906
  * [25] Georges Aad  “Measurement of long-range pseudorapidity correlations and azimuthal harmonics in $\sqrt{s_{NN}}=5.02$ TeV proton-lead collisions with the ATLAS detector”  In _Phys. Rev. C_ 90.4, 2014, pp. 044906  DOI: 10.1103/PhysRevC.90.044906
  * [26] Georges Aad  “Transverse momentum and process dependent azimuthal anisotropies in $\sqrt{s_{\mathrm{NN}}}=8.16$ TeV $p$+Pb collisions with the ATLAS detector”  In _Eur. Phys. J. C_ 80.1, 2020, pp. 73  DOI: 10.1140/epjc/s10052-020-7624-4
  * [27]
  * [28] S. Cao  “Determining the jet transport coefficient q̂ from inclusive hadron suppression measurements using Bayesian parameter estimation”  In _Phys. Rev. C_ 104.2, 2021, pp. 024905  DOI: 10.1103/PhysRevC.104.024905
  * [29] S. Cao  “Determining the jet transport coefficient q̂ from inclusive hadron suppression measurements using Bayesian parameter estimation”  In _Phys. Rev. C_ 104.2, 2021, pp. 024905  DOI: 10.1103/PhysRevC.104.024905
  * [30] A. Kumar  “Inclusive jet and hadron suppression in a multistage approach”  In _Phys. Rev. C_ 107.3, 2023, pp. 034911  DOI: 10.1103/PhysRevC.107.034911
  * [31] W. Fan  “Multiscale evolution of charmed particles in a nuclear medium”  In _Phys. Rev. C_ 107.5, 2023, pp. 054901  DOI: 10.1103/PhysRevC.107.054901
  * [32] Y. Tachibana  “Hard Jet Substructure in a Multi-stage Approach”, 2023  arXiv:2301.02485 [hep-ph]
  * [33]
  * [34] D. Everett  “Multisystem Bayesian constraints on the transport coefficients of QCD matter”  In _Phys. Rev. C_ 103.5, 2021, pp. 054904  DOI: 10.1103/PhysRevC.103.054904
  * [35] D. Everett  “Multisystem Bayesian constraints on the transport coefficients of QCD matter”  In _Phys. Rev. C_ 103.5, 2021, pp. 054904  DOI: 10.1103/PhysRevC.103.054904
  * [36] D. Everett  “Phenomenological constraints on the transport properties of QCD matter with data-driven model averaging”  In _Phys. Rev. Lett._ 126.24, 2021, pp. 242301  DOI: 10.1103/PhysRevLett.126.242301
  * [37] Michael Kordell and Abhijit Majumder  “Jets in d(p)-A Collisions: Color Transparency or Energy Conservation”  In _Phys. Rev. C_ 97.5, 2018, pp. 054904  DOI: 10.1103/PhysRevC.97.054904
  * [38]
  * [39] S. Cao  “Multistage Monte-Carlo simulation of jet modification in a static medium”  In _Phys. Rev. C_ 96.2, 2017, pp. 024909  DOI: 10.1103/PhysRevC.96.024909
  * [40] S. Cao  “Multistage Monte-Carlo simulation of jet modification in a static medium”  In _Phys. Rev. C_ 96.2, 2017, pp. 024909  DOI: 10.1103/PhysRevC.96.024909
  * [41] J. H. Putschke  “The JETSCAPE framework”, 2019  arXiv:1903.07706 [nucl-th]
  * [42] https://github.com/JETSCAPE/X-SCAPE-COMP
  * [43]
  * [44] Chun Shen and Björn Schenke  “Dynamical initial state model for relativistic heavy-ion collisions”  In _Phys. Rev. C_ 97.2, 2018, pp. 024907  DOI: 10.1103/PhysRevC.97.024907
  * [45] Chun Shen and Björn Schenke  “Dynamical initial state model for relativistic heavy-ion collisions”  In _Phys. Rev. C_ 97.2, 2018, pp. 024907  DOI: 10.1103/PhysRevC.97.024907
  * [46] Chun Shen and Björn Schenke  “Longitudinal dynamics and particle production in relativistic nuclear collisions”  In _Phys. Rev. C_ 105.6, 2022, pp. 064905  DOI: 10.1103/PhysRevC.105.064905