# Employing online quantum random number generators for generating truly
random quantum states in _Mathematica_

Jarosław Adam Miszczak
<EMAIL_ADDRESS>
Institute of Theoretical and Applied Informatics,
Polish Academy of Sciences, Bałtycka 5, 44-100 Gliwice, Poland

###### Abstract

We present a new version of TRQS package for _Mathematica_ computing system.
The package allows harnessing quantum random number generators (QRNG) for
investigating the statistical properties of quantum states. It implements a
number of functions for generating random states. The new version of the
package adds the ability to use the on-line quantum random number generator
service and implements new functions for retrieving lists of random numbers.
Thanks to the introduced improvements, the new version provides faster access
to high-quality sources of random numbers and can be used in simulations
requiring large amount of random data.

Keywords: random density matrices ;quantum information ;quantum random number
generator ;on-line quantum random number generator service

PACS numbers: 03.67.-a ;02.70.Wz ;07.05.Tp

## New Version Program Summary

_Program title:_ TRQS
_Program author:_ Jarosław A. Miszczak
_Program web page:_ http://www.iitis.pl/ miszczak/trqs
_Distribution format:_ tar.gz
_No. of bytes in distributed program, including test data, etc.:_ $218\ 618$
_No. of lines in distributed program, including test data, etc.:_ $16\ 479$
_Programming language:_ Mathematica, C
_Computer:_ any supporting Mathematica in version 7 or higher
_Operating system:_ any platform supporting Mathematica; tested with GNU/Linux
(32 and 64 bit)
_RAM:_ case-dependent
_Classification:_ 4.15
_Catalogue identifier of previous version:_ AEKA_v1_0
_Reference in CPC:_ Comput. Phys. Commun. 183 (2012) 118 [1]
_Does the new version supersede the previous version?:_ Yes.
_Nature of problem:_ Generation of random density matrices and utilization of
high-quality random numbers for the purpose of computer simulation.
_Method of solution:_ Use of a physical quantum random number generator and an
on-line service providing access to the source of true random numbers
generated by quantum real number generator.
_Running time:_ Depends on the used source of randomness and the amount of
random data used in the experiment.
_Reasons for new version:_ Added support for the high-speed on-line quantum
random number generator and improved methods for retrieving lists of random
numbers.

##  1 Summary of revisions

The presented version provides two significant improvements.

The first one is the ability to use on-line Quantum Random Number Generation
service developed by PicoQuant GmbH and the Nano-Optics groups at the
Department of Physics of Humboldt University. The on-line service supported in
the version 2.0 of the TRQS package provides faster access to true randomness
sources constructed using the laws of quantum physics. The service is freely
available at [arng]. The use of this service allows using the presented
package with the need of a physical quantum random number generator.

The second improvement introduced in this version is the ability to retrieve
arrays of random data directly for the used source. This increases the speed
of the random number generation, especially in the case of an on-line service,
where it reduces the time necessary to establish the connection. Thanks to the
speed improvement of the presented version, the package can now be used in
simulations requiring larger amount of random data. Moreover, the functions
for generating random numbers provided by current version of the package more
closely follow the pattern of functions for generating pseudo-random numbers
provided in _Mathematica_.

##  2 Speed comparison

The implementation of the support for QRNG on-line service provides a
noticeable improvement in the speed of random number generation. For the
samples of real numbers of size $10^{1},10^{2},\dots,10^{7}$ the times
required to generate these samples using Quantis USB device and QRNG service
are compared in Fig. 1.

Figure 1: Comparison of the times required to generate a sample of random
numbers using Quantis device (red squares) and QRNG service (blue dots) in
log-log scale. In the case of the QRNG service the plot presents time averaged
over three experiments. In both cases the results were obtained by generating
matrices of random real numbers from $[0,1]$ using TrueRandomReal function.

The presented results show that the use of on-line service provides faster
access to random numbers. One should note, however, that the speed gain can
increase or decrease depending on the connection speed between the computer
and the server providing random numbers.

## Acknowledgements

Author would like to thank R. Heinen and M. Wahl for bringing his attention to
Ref. [3] and motivating this work and to C. Phillips and helping during the
development of the package. This work was supported by the Polish Ministry of
Science and Higher Education under the grant number IP2011 036371.

## References

  * [1] J.A. Miszczak, _Generating and using truly random quantum states in Mathematica_ , Computer Physics Communications, Vol. 183, No. 1 (2012), pp. 118-124. arXiv:1102.4598 DOI:10.1016/j.cpc.2011.08.002
  * [2] QRNG Service available on-line at https://qrng.physik.hu-berlin.de/.
  * [3] M. Wahl, M. Leifgen, M. Berlin, T. Röhlicke, H.-J. Rahn, O. Benson O., _An ultrafast quantum random number generator with provably bounded output bias based on photon arrival time measurements_ , Applied Physics Letters, Vol. 098, 171105 (2011). DOI:10.1063/1.3578456


# BEAM–BEAM EFFECTS IN THE HIGH-PILE-UP TESTS OF THE LHC

G. Trad  CERN  Geneva  Switzerland

###### Abstract

Investigating the beam–beam limit in the LHC is of great importance, since
identifying its source is crucial for the luminosity optimization scenario.
Several experiments were carried out to search for this limit and check
whether it is dominated by the head-on (HO) or the long-range (LR)
interactions. In this paper only the HO collision effects will be considered,
tracking the evolution of the maximum tune shift achieved during the dedicated
machine developments and the special high pile-up fills.

##  1 INTRODUCTION

Contrarily to electron machines, hadron machines are mostly limited by non-
linear effects and lifetime problems, in particular for beams with high
intensities [1]. The limitation appears as a slow emittance increase (over
hours) or beam losses (tails and dynamic aperture), bad lifetime or other
effects such as coherent beam–beam oscillations.

It is important to note that the head-on (HO) tune shift depends only on the
bunch intensity and the normalized emittance, i.e. is independent of the beta
function at the interaction point $\beta^{*}$ and the energy as shown in Eq.
(1). For the long-range (LR) interactions, the tune shift depends on the beam
separation $d_{\mathrm{sep}}$, thus on the crossing angle $\alpha$ and the
beam energy $\gamma$, as shown in Eq. (3).

$\Delta Q_{\mathrm{HO}}\propto\frac{N}{\varepsilon_{n}},$ (1)

$L\propto\frac{N^{2}}{\varepsilon_{n}}=\Delta Q_{\mathrm{HO}}\cdot N,$ (2)

$\Delta
Q_{\mathrm{LR}}\propto\frac{N}{d_{\mathrm{sep}}^{2}}=\frac{N\cdot\varepsilon_{n}}{\alpha^{2}\cdot\beta^{*}\cdot\gamma}.$
(3)

Therefore, identifying the source of these limits has vital significance for
the luminosity optimization strategy: if the HO collisions are dominating the
beam–beam limit, it is then advantageous to increase the bunch intensity
together with the transverse emittance since this would keep the tune shift
unaffected and increases the luminosity proportionally to the intensity (see
Eq. (2)). The luminosity is further increased by pushing more the focusing at
the interaction point (IP) without affecting the beam–beam parameter $\xi$:

$\xi_{21x,y}=\frac{N_{1}r_{p}}{2\pi\sqrt{\varepsilon_{x,y_{n1}}}\left(\sqrt{\varepsilon_{x,y_{n1}}}+\sqrt{\varepsilon_{y,x_{n1}}}\right)}.$
(4)

On the other hand, when the machine is limited by the LR interactions, a large
number of bunches with a moderate $\beta^{*}$ is preferred. Since $\Delta
Q_{LR}$ depends on $\beta^{*}$, any change of the optical function at the IP
requires one to adjust the crossing angle $\alpha$ to keep the LR tune shift
constant,

$\alpha\propto\sqrt{\frac{N\cdot\varepsilon_{n}}{\Delta
Q_{\mathrm{max}}\cdot\beta^{*}\cdot\gamma}}.$ (5)

Some confusion is related to the maximum expected HO beam–beam tune shift for
the LHC. The nominal HO tune shift was derived from the experience of the
Super Proton Synchrotron (SPS), taking into account possible contributions
from the lattice non-linearities and significant LR contributions. The nominal
value of $\xi=0.0037$ was defined to provide a coherent set of parameters to
reach the target luminosity, $L=10^{34}$ cm-2 s-1 [2, 3]. Therefore, it should
be considered as a conservative value and not as a real upper limit, in
particular in the absence of strong LR interactions.

##  2 Is the ‘nominal’ beam–beam HO tune shift reachable?

Early in 2010, even if no show-stoppers were expected, dedicated experiments
and observations during normal operation were planned to check the feasibility
of colliding beams with the nominal tune shift estimated in the design report.
In the following, we briefly list a few of these fills including the important
results observed [4].

###  2.1 Fill 1069 (May 2010)

####  2.1.1 Experimental setup

Beam 1B1-b1IP1 vs.  B2-b1IP5 vs.  B2-b1IP8 vs.  B2-b892B1-b1786IP8 vs.  B2-b1

Beam 2B2-b1IP1 vs.  B1-b1IP5 vs.  B1-b1B2-b892IP2 vs.  B1-b1IP8 vs.  B1-b1786

Figure 1: Collision scheme of the injected bunches in Fill 1069. Shown are the
colliding bunch pairs in the two beams in the four interaction points. The
slot numbers are given to specify a bunch position.

For the LHC fill in study (Fill 1069), two ‘almost’ nominal bunches of
intensity $\sim$1$\times$ $10^{11}$ ppb and normalized emittance of 3 $\mu$m
were injected in each ring.
The bunches were put in collision according to the standard operational cycle,
with $\beta*=10$ m in all IPs and the collision tunes ($Q_{x}=0.31$,
$Q_{y}=0.32$), at the injection energy of 450 GeV. The collision scheme is
shown in Fig. 1.

Figure 2: Total beam intensity evolution during Fill 1069. Beam 1 is shown in
blue while Beam 2 is shown in yellow.

####  2.1.2 Observations

When the separation bumps were ‘collapsed’ simultaneously to bring the beams
into collision, the lifetime dropped, especially in Beam 1, but a small tune
adjustment ($\Delta Q_{h}\sim+0.005$) was sufficient to stabilize it at 25 h
(Fig. 2). The resulting beam–beam tune shift exceeded the nominal value
reaching $\xi=0.004$ per IP. An increase of the vertical emittance in Beam 2
was observable and thought to be originating from an external excitation not
fully understood (the hump). Since the results were encouraging, and no
limitations were expected for the HO tune shift, it was decided to increase
the intensity of the colliding bunches exploring regions with a higher tune
shift.

###  2.2 Fill 1765 (May 2011)

####  2.2.1 Experimental setup

One high-intensity bunch was injected per beam ($\sim$1.6$\times$ $10^{11}$
ppb); the measured normalized emittance by the wire scanners (WSs) was
$\varepsilon_{{}_{x,y}}=1.2$ $\mu$m. With the collision tunes setting, the
beams were brought into collision at injection energy. The IP settings were
the following: $\beta*=11$ m and nominal crossing angles in all IPs, while the
spectrometers were off in IP2 and IP8. The transverse damper (ADT) was turned
on only at injection and set off afterwards [5, 6].

####  2.2.2 Observations

An increase in the vertical emittance once in collision was observed,
resulting in a 2.2 $\mu$m emittance and $\xi=0.009$/IP. Since the bunches were
colliding only in IP1 and IP5, a total $\xi=0.018$ was reached. A small tune
scan was tried as well at the end of the fill, to search for a better working
point. No lifetime effects were observed, just a minor emittance blow-up for
Beam 2.

###  2.3 Fill 1766: Part 1 (May 2011)

Following the encouraging beam–beam parameter reached in Fill 1765, it was
decided to maintain the high-intensity bunches and increase the number of
collisions. Therefore, two bunches per beam were injected, allowing collisions
in all interaction points but with a different collision pattern, as shown in
Fig. 3.

Beam 1B1-b100IP1 vs.  B2-b100IP2 vs.  B2-b991IP5 vs.  B2-b100B1-b1885IP8 vs.
B2-b991

Beam 2B2-b100IP1 vs.  B1-b100IP5 vs.  B1-b100B2-b991IP2 vs.  B1-b100IP8 vs.
B1-b1885

Figure 3: Collision scheme of the injected bunches in Fill 1766.

####  2.3.1 Observations

A signiﬁcant emittance increase was observed when beams started to collide in
IP8. This emittance increase is observed on the bunches colliding only in IP8,
i.e. with a single collision. The bunch with the largest number of collisions
had the smallest emittance increase. The emittances as a function of time
measured with the synchrotron light monitor (BSRT) are shown in Figs. 4 and 5.

Figure 4: Fill 1766: Beam 1 horizontal (upper plot) and vertical (lower plot)
normalized emittance evolution for the bunches $B1-b100$ (black) and $B1-b991$
(green), as measured by the BSRT. Figure 5: Fill 1766: Beam 2 horizontal
(upper plot) and vertical (lower plot) normalized emittance evolution for the
bunches $B2-b100$ (black) and $B2-b991$ (green), as measured by the BSRT.

This emittance increase did not start when the bunches were put into
collision, but during the optimization process, i.e. luminosity scan. It was
done manually since the automatic procedure did not find an optimum.

###  2.4 Fill 1766: Part 2 (May 2011)

With the same machine configuration as in Fill 1765, again two single bunches
colliding only in IP1 and 5 were injected to study in more detail the
parameter space, in particular the effect of the working point on the initial
emittance growth observed once the beams are put into collision. Two fillings
were done using the same filling scheme with roughly the same injected
intensities ($\sim$1.8$\times$ $10^{11}$ ppb) and normalized emittance
($\sim$1.32 $\mu$m), but changing the tunes starting working point. As shown
in Figs. 6 and 7, it is seen that moving the tune point from (0.31, 0.32) to
(0.31, 0.31) helped to increase the beams’ lifetime. The initial low lifetime
observed was thought to originate from a lattice resonance where the core of
the bunch was getting close to the 10th-order resonance; a small negative
shift for the vertical tune ($\Delta Q_{v}=-0.01$) improved the losses and
reduced the emittance blow-up since the large tune shift from the collisions
will bring the core particles below the resonance.

Figure 6: Fill 1766: horizontal (blue) and vertical (red) tune and the
lifetime (black) evolution with the (0.31, 0.32) tune working point for Beam 1
(upper plot) and Beam 2 (lower plot). Figure 7: Fill 1766: horizontal (blue)
and vertical (red) tune and the lifetime (black) evolution with the (0.31,
0.31) tune working point for Beam 1 (upper plot) and Beam 2 (lower plot).

Once this emittance blow-up was cured, a beam–beam linear parameter of
$\xi=0.017$/IP was achieved resulting in a total value $\xi=0.034$.

##  3 High pile-up tests: single bunches

To provide a calibration test for the LHC experiments, tests were performed to
generate a high pile-up in the detectors, i.e. a large number of events per
crossing [7, 8]. For this purpose a large number of bunches is not needed and
was done with single bunches per beam.

###  3.1 Experimental Setup

For the LHC fill in study (Fill 2201), one high-intensity bunch
($\sim$2.4$\times$ $10^{11}$ ppb) was injected in the bucket 1001 of both
rings. The bunches were accelerated to the flat-top energy (3.5 TeV) and put
in collision in ATLAS and CMS according to the standard operational cycle,
where the machine was in the standard configuration:
$\beta*=1$ m in IP1/5 and $\pm$120 $\mu$rad crossing angle,
$\beta*=10$ m in IP2 and $\pm$80 $\mu$rad crossing angle,
$\beta*=3$ m in IP8 and $\pm$250 $\mu$rad crossing angle.

###  3.2 Observations

The resulting luminosity was $\sim 4.7\times 10^{-30}$ cm-2 s-1 corresponding
to a pile-up (number of inelastic interactions per crossing) of $\mu\sim 31$
encounters per turn. The nominal number is $\mu\sim 20$ [2]. Observations of
the beam intensity and emittance showed the following:

####  3.2.1 Emittance

The bunch transverse emittances were frequently measured with the WSs
throughout the fill. At the injection energy an emittance growth was observed
only in the horizontal plane (12% for Beam 1 and 8% for Beam 2). Through the
ramp the beams’ emittances in both planes steadily grew reaching $\sim$3.1
$\mu$m at the beginning of the collisions (a growth of $44\%$ B1 and $27\%$ B2
in the horizontal plane and $30\%$ B1 and $7\%$ B2 in the vertical plane). In
about 3 h in collision the horizontal emittances grew by about 20% for both
beams while the vertical ones increased by $\sim 13\%$.

####  3.2.2 Intensity

Figure 8: Fill 2201: measured vs. fitted bunch intensity for Beam 1 (upper
plot) and Beam 2 (lower plot). The points refer to the fBCT measurement and
the lines are the linear fits representing the loss regimes observed. The
slopes of the linear fits corresponding to the constant loss rate are also
reported.

The beam intensity decay was studied and two loss regimes were observed for
Beam 1: the first regime ended after 1.5 h in collision with a loss rate of
$\sim 4.1\times 10^{9}$ p/h followed until the end of the fill with a constant
loss regime of $\sim 6.6\times 10^{9}$ p/h. For Beam 2, a constant loss rate
was recorded throughout all the fill $\sim 6\times 10^{9}$ p/h. A correlation
was observed between the transition between the two loss regimes in Beam 1 and
its bunch length growth [9], where, once the value of $\sim 8.95$ cm was
reached, corresponding to the initial bunch length of Beam 2, the losses
settled to their maximum value (Fig. 8).

####  3.2.3 Luminosity

A validation of the measured beam parameters was made through a comparison
between the published instantaneous luminosity from the experiments and the
analytically calculated one. The observed discrepancy was compatible with the
uncertainties in the values of the machine parameters ($\beta^{*}$) and the
monitors’ accuracy.

In addition, the luminosity evolution model developed for the TEVATRON was
applied using the LHC parameters, trying to predict the evolution of the beam
parameters. An agreement was only found for the period corresponding to the
first loss regime for Beam 1, while for the second loss regime and for the
Beam 2 intensity evolution, the longitudinal losses were underestimated almost
by a factor of 3.

##  4 2011 High-pile-up tests: multiple bunches

After the single-bunch high-pile-up tests, the experimenters requested a new
test to check whether multiple high-pile-up collisions can be effectively
processed.

###  4.1 Experimental Setup

Fill 2252 was dedicated to providing ATLAS and CMS with high-pile-up
collisions. A peak pile-up $\mu$ of almost 35 was achieved. Part of the test
was also to bring the pile-up to values $<0.5$ and this was achieved by
separating the beams in steps. The head-on collisions were restored after each
IP separation to study its effects. The collision scheme during this
experiment is shown in Fig. 9. The number of collisions is shown for every
bunch.

The bunch parameters (intensity, emittance and length) of high-pile-up
colliding bunches were observed while putting beams in collision with the same
machine setup as in Fill 2201 but with the following filling scheme:

Figure 9: Collision scheme of the injected bunches in Fill 2252. All bunches
are colliding HO in IP1 and 5, except for $B1-b995$ and $B2-b992$ colliding
respectively in IP8 with $B2-b101$ and in IP2 with $B1-b101$. The number of LR
encounters is shown as well.

###  4.2 Observations

After $\sim$ 100 min in single-bunch mode, the beams were re-separated into
three different phases, as shown in Fig. 10:

Figure 10: ATLAS and CMS instantaneous luminosity during the separation phases
shown respectively in blue and green.

  * –

Five separation steps in the horizontal plane of IP5 reducing the luminosity
from $\sim$ 44 Hz/$\mu$b to $\sim$ 10 Hz/$\mu$b,

  * –

Eleven separation steps in the vertical plane of IP1 reducing the luminosity
from $\sim$ 43 Hz/$\mu$b to $\sim$ 0.08 Hz/$\mu$b,

  * –

Two separation steps in the horizontal plane of IP5 reducing the luminosity
from $\sim$ 41 Hz/$\mu$b to $\sim$ 1.23 Hz/$\mu$b.

After each separation, the instantaneous luminosity was re-optimized to its
maximum value in the considered IP.

####  4.2.1 Emittance

The bunches can be classified into three categories:

  * –

bunches presenting a higher growth rate in the vertical plane during the
separation phases: $B2-b101$, $B2-b421$ and $B1-b1061$ (see Fig. 11),

Figure 11: Vertical emittances of bunches 101–421 (Beam 2) showing a
continuous smooth increase in the emittance while separating the beams in IP5.
The vertical coloured bar delimits the separation steps of every separation
phase in both IPs.

  * –

bunches presenting sudden blow-up throughout the separation steps: $B1-b101$,
$B1-b421$, $B2-b992$ and $B2-b1061$ (see Fig. 12).

Figure 12: Vertical blow-up of the vertical emittance of bunches 101–421 (Beam
1) and bunches 992–1061 (Beam 2) while separating the beams in IP5. The
vertical coloured bar delimits the separation steps of every separation phase
in both IPs.

  * –

bunches not affected by the separation steps.

The sudden increases in emittance were observed mainly in the vertical plane
during the steps of the first CMS separation (recall that in IP5 the vertical
plane is the separation plane).

####  4.2.2 Lifetime

Normalizing the bunch intensity curves to the initial value at the start of
collisions allows one to identify three groups of behaviours according to the
number of HO collisions the bunches experience (three HO, two HO, one HO).

Analysing the intensity data for the first 1.8 h, it is possible to
disentangle the separation effects from the initial evolution of the bunch
parameters.

Once in collision, for most bunches a linear decrease in intensity is observed
during the considered time period; a linear fit is applied to the measured
intensities and the slopes corresponding to the loss rates are shown in Fig.
13.

Figure 13: Slopes from linear fits of the bunch intensity decays for the first
100 min.

The highest losses are observed for bunches experiencing three HO collisions
while the lowest losses are for bunches with one HO collision only.

We observe the losses for each of the three separation phases mentioned
before; the loss trend can be organized into families (according to the
collision scheme) for all the steps.

  * –

_Separation phase 1 (IP5)_
The first separation step in IP5 had an effect on the luminosity in IP1, where
a drop of 3.5% was observed. No significant change in the intensity evolution
of all the bunches in both beams except for:

    * .

bunch 101 in Beam 1 (having three HO: IP1 and 5 with bunch 101 in Beam 2, IP2
with bunch 992 in Beam 2),

    * .

bunch 992 in Beam 2 (having one HO: IP2 with bunch 101 in Beam 1).

An increase in the loss rate is observed for $B1-b101$ after the first
separation step as shown in Fig. 14. For $B2-b992$, a sudden increase in
losses is seen in the first and last steps of the separation. Note that at the
end of the fifth step, once the beams are put back in HO collisions in IP5,
the loss rate for $B2-b992$ returned to the value it had before the beams’
manipulation as seen in Fig. 14. In order to explain and understand the link
between the behaviour of $B2-b992$, colliding only in IP2, and the separations
in IP5 it is important to consider the LR collisions this bunch is
experiencing in IP5 with $B1-b995B1$.

Figure 14: The evolution of the normalized intensity of $B2-b992$ (red) during
the first separation phase in IP5. The instantaneous luminosity in IP1 and IP5
is shown respectively in black and green.

  * –

_Separation phase 2 (IP1)_
$B1-b995$ and $B2-b992$ having one HO collision in IP8 and 2, respectively,
were not affected by the steps of this separation. For the three HO collisions
family, it can be observed that $B1-b101B1$ experienced higher loss rate at
half total separation and kept this loss regime until the beams were brought
back into HO collision; bunch $B2$-$b101$ entered a higher loss regime after
the first separation step and maintained this rate once back in the initial
situation of fully HO collisions. The slopes of the intensity evolution for
all the other bunches (two HO families) depend on the value of the separation
between the beams: lower losses were observed for higher separations. It is
worth observing that the bunch losses for both families in Beam 1 and Beam 2
colliding in IP1 and 5 at zero separation were higher at the end of this phase
(especially for Beam 1).

  * –

_Separation phase 3 (IP5)_
No clear variation in the intensity decay was observed for the bunches having
one HO and three HO in both beams, while for all the other bunches, the losses
were decreasing for an increasing beam separation. Again for some bunches in
Beam 1 having two HO, the initial loss rate was not restored at the end when
the beams were brought into HO collisions: higher losses were still observed
at the end.

####  4.2.3 Bunch length

We only highlight the bunch length evolution of Beam 2. It is worth pointing
out the behaviour of $B2-b992$ (see Fig. 15): a bunch shortening is observed
simultaneously to the intensity loss, corresponding to the first separation
step in IP5 (first separation phase). For all the other bunches, the bunch
length did not bear signs of the losses and the emittance blow-up (horizontal
and vertical).

Figure 15: The observed parameters for $B2-b992$ (intensity, transverse
emittances and bunch length) from the start of the fill until the end of the
collisions. The vertical coloured bar delimits the separation steps of every
separation phase in both IPs.

##  5 2012 high-pile-up

Profiting from the new optics in the SPS (Q20 optics), where the fractional
part of the tune was moved to 0.2, very high brightness bunches were put into
collision in the LHC aiming to establish a new record pile-up possibly up to
100. High-brightness bunches with intensity of 3 $\times$ $10^{11}$ ppb and
normalized emittance of 2 $\mu$m were used for this test. The energy ramp was
troublesome; the controlled longitudinal blow-up needed for the beam stability
was faulty, as shown in Fig. 16. The bunch length of Beam 2 was brought to the
target bunch length of 1.3 ns while Beam 1 bunch length went down to $\sim
0.4$ ns before reaching a value close to its target value of 1.2 ns. This
caused a deterioration in the beam quality throughout the ramp. The bunches
were put into collision with more than 10% of losses (end of ramp and at flat-
top energy) and a doubled emittance, as shown in Fig. 17.

Figure 16: Fill 2824: bunch length evolution for Beam 1 (blue) and Beam 2
(red) during the fill as measured by the beam-quality monitors in the LHC.
Figure 17: Fill 2824: intensity evolution of Beam 1 (magenta) and Beam 2
(green) along with energy ramp (red). The value of $\beta^{*}$ in IP1 is shown
as well to indicate the squeeze phase, and the instantaneous luminosity (blue)
in IP1 is shown as well as an indication for the start of collisions.

The resulting pile-up of 58 was the best reached in two tries. Once the issue
with the longitudinal blow-up was solved, two bunches per beam colliding only
in IP1 and IP5 were injected and with difficulty brought into collision at 4
TeV with a $\beta^{*}$ of 60 cm in both IPs. An instability arising during the
squeeze was observed leading again to some losses along with an emittance
increase. The maximum pile-up obtained was 70 in IP1 and 65 in IP5.

##  6 Summary

In this paper, we reported the main studies of the HO beam–beam effects in the
LHC with nominal, high-intensity, high-brightness single and multiple bunches.
It is worth noting that small contributions of the lattice non-linearities as
well as good settings of the machine allowed one to quickly reach the nominal
head-on beam–beam tune shift; it was also shown that the LHC allows very large
head-on tune shifts above nominal.
It has to be seen whether this can be maintained in the presence of many LR
interactions. Yet there is no reason to assume that an HO limit has been
reached for the moment in the LHC.

##  7 ACKNOWLEDGMENTS

I would like to thank all the CERN beam–beam working group members for their
contribution, help and advice, especially W. Herr and G. Papotti.

## References

  * [1] W. Herr et al., “Observations of Beam–Beam Effects at high intensity at the LHC,” Proceedings of International Particle Accelerator Conference 2011, San Sebastian Spain (2011).
  * [2] LHC Design Report, CERN-2004-003 (2004).
  * [3] W. Herr, “LHC Beam–Beam Effects – Review and Outlook”, Proc. LHC Operations Workshop, Evian, December 2010 (2010).
  * [4] G. Arduini, W. Herr and J.M. Jowett, “First Observation of Beam–Beam Interactions in High Intensity Collisions at the LHC,” CERN-ATS-Note-2010-039 PERF.
  * [5] W. Herr et al., “Head-on Beam–Beam Tune Shift with High Brightness Beams in the LHC,” CERN-ATS-Note-2011-029 MD.
  * [6] M Albert et al., “Head-on Beam–Beam Collisions with High Intensities and Long Range Beam–Beam Studies in the LHC,” CERN-ATS-Note-2011-058 MD.
  * [7] G. Trad et al., “Beam Parameters Observations During a High Pile-up Collisions fill,” CERN-ATS-Note-2011-105 MD.
  * [8] G. Trad et al., “Beam Parameters Observations During the Second High Pile-up Collisions fill in 2011,” CERN-ATS-Note-2012-010 MD.
  * [9] G. Papotti, T. Bohl, F. Follin, and U. Wehrle, “Longitudinal Beam Measurements at the LHC: the LHC Beam Quality Monitor,” IPAC11, San Sebastian, September 2011.


# Constraining the state of the intergalactic medium during the Epoch of
Reionization using MWA 21-cm signal observations

Raghunath Ghara1,2, Sambit K. Giri3, Benedetta Ciardi4, Garrelt Mellema5, and
Saleem Zaroubi1,6,2
1ARCO (Astrophysics Research Center), Department of Natural Sciences, The Open
University of Israel, 1 University Road, PO Box 808, Ra’anana 4353701, Israel
2Department of Physics, Technion, Haifa 32000, Israel
3Institute for Computational Science, University of Zurich,
Winterthurerstrasse 190, 8057 Zurich, Switzerland
4Max-Planck Institute for Astrophysics, Karl-Schwarzschild-Straße 1, 85748
Garching, Germany
5The Oskar Klein Centre, Department of Astronomy, Stockholm University,
AlbaNova, SE-10691 Stockholm, Sweden
6Kapteyn Astronomical Institute, University of Groningen, PO Box 800, 9700AV
Groningen, the Netherlands
E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The Murchison Widefield Array (MWA) team has derived new upper limits on the
spherically averaged power spectrum of the 21-cm signal at six redshifts in
the range $z\approx 6.5-8.7$. We use these upper limits and a Bayesian
inference framework to derive constraints on the ionization and thermal state
of the intergalactic medium (IGM) as well as on the strength of a possible
additional radio background. We do not find any constraints on the state of
the IGM for $z\gtrsim 7.8$ if no additional radio background is present. In
the presence of such a radio background, the 95 per cent credible intervals of
the disfavoured models at redshift $\gtrsim 6.5$ correspond to an IGM with a
volume averaged fraction of ionized regions below 0.6 and an average gas
temperature $\lesssim 10^{3}$ K. In these models, the heated regions are
characterised by a temperature larger than that of the radio background, and
by a distribution with characteristic size $\lesssim 10$ $h^{-1}$ Mpc and a
full width at half maximum (FWHM) of $\lesssim 30$ $h^{-1}$ Mpc. Within the
same credible interval limits, we exclude an additional radio background of at
least $0.008\%$ of the CMB at 1.42 GHz.

###### keywords:

radiative transfer - galaxies: formation - intergalactic medium - cosmology:
theory - dark ages, reionization, first stars - X-rays: galaxies

††pubyear: ††pagerange: Constraining the state of the intergalactic medium
during the Epoch of Reionization using MWA 21-cm signal
observations–References

##  1 INTRODUCTION

During the Epoch of Reionization (EoR) ultra-violet (UV) radiation from the
earliest generations of galaxies ionized the primordial neutral hydrogen (H i
) in the intergalactic medium (IGM). This global phase transition permanently
changed the state of the IGM and had far reaching consequences for the
subsequent formation and evolution of galaxies. Observations of the absorption
spectra of quasars at redshift $z\gtrsim 6$ (e.g. Fan et al., 2006; Bañados et
al., 2018; Davies et al., 2018; Greig et al., 2017; Greig et al., 2019; Wang
et al., 2020; Yang et al., 2020) and the measurement of the Thomson scattering
optical depth from observations of the Cosmic Microwave Background (CMB)
(Planck Collaboration et al., 2020) suggest that this transition took place
sometime between redshift $\approx$ 6 and 10 (e.g, Mitra et al., 2015).
However, many of its details, such as the exact timing of the process, the
ionization and thermal state of the IGM, the morphology of the ionized
regions, remain poorly understood (see e.g. Pritchard & Loeb, 2012; Zaroubi,
2013, for reviews).

The 21-cm signal produced by H i in the IGM is the most direct probe of the
reionization process as it depends directly on the H i density. It is
therefore expected to answer many of the outstanding questions, including the
nature of the sources responsible for reionization. Thus, substantial efforts
have been made to measure this signal. Observations of the 21-cm signal fall
into two categories. The first type uses large interferometric radio
telescopes to measure the statistical fluctuations of the signal. Examples of
such telescopes are the Low Frequency Array (LOFAR)111http://www.lofar.org/
(van Haarlem et al., 2013; Patil et al., 2017), the Precision Array for
Probing the Epoch of Reionization (PAPER)222http://eor.berkeley.edu/ (Parsons
et al., 2014; Kolopanis et al., 2019), the Murchison Widefield Array
(MWA)333http://www.mwatelescope.org/ (e.g. Tingay et al., 2013; Wayth et al.,
2018), the New extension in Nançay Upgrading loFAR
(NenuFAR)444https://nenufar.obs-nancay.fr/ (Zarka et al., 2012), the Upgraded
Giant Metrewave Radio Telescope (uGMRT) (Gupta et al., 2017), the Owens Valley
Radio Observatory Long Wavelength Array (OVRO-
LWA)555http://www.tauceti.caltech.edu/LWA/ and the Hydrogen Epoch of
Reionization Array (HERA)666https://reionization.org/ (DeBoer et al., 2017). A
more sensitive radio interferometer, the upcoming Square Kilometre Array
(SKA)777http://www.skatelescope.org/, should be able to also produce
tomographic images of the distribution of the signal on the sky (Mellema et
al., 2015; Ghara et al., 2017; Giri et al., 2018a, b; Giri, ; Ghara &
Choudhury, 2020; Giri & Mellema, 2020; Kapahtia et al., 2021). The second type
of observations use single antennae and attempt to measure the sky-averaged
21-cm signal and its evolution with time. Examples are EDGES (Bowman & Rogers,
2010), EDGES2 (Monsalve et al., 2017; Bowman et al., 2018), SARAS (Patra et
al., 2015), SARAS2 (Singh et al., 2017), SARAS3 (Nambissan et al., 2021),
PRIZM (Philip et al., 2019), BigHorns (Sokolowski et al., 2015), SciHi (Voytek
et al., 2014) and LEDA (Price et al., 2018).

One of the main challenges of these experiments is to remove the contribution
of the galactic and extra-galactic foregrounds, which is stronger than the
expected 21-cm signal by several orders of magnitude. Their smooth frequency
dependence should allow a separation from the fluctuating 21-cm signal,
allowing them to be either subtracted (Harker et al., 2009; Bonaldi & Brown,
2015; Chapman et al., 2016; Mertens et al., 2018; Hothi et al., 2021),
suppressed (Datta et al., 2007; Majumdar et al., 2012; Ghara et al., 2016) or
avoided (Datta et al., 2010; Liu et al., 2014). Moreover, an exquisite
calibration of the system is required to minimize the artefacts from strong
sources (Barry et al., 2016; Patil et al., 2017), and a long observation time
is necessary to reduce the instrumental noise. Further calibration challenges
may be posed by the hardware components of the telescope (see e.g., Kern et
al., 2019) and by the temporally and spatially varying ionosphere (see e.g.,
Mevius et al., 2016), to mention a few.

So far, no undisputed detection of the 21-cm signal from the EoR has been
achieved. Bowman et al. (2018) analysed observations with the EDGES2 low-band
antenna and claimed a detection of the averaged brightness temperature of the
21-cm signal ($\overline{\delta T_{\rm b}}$) at various redshifts. The
redshift evolution of $\overline{\delta T_{\rm b}}$ has a minimum at $z\approx
17$ with an amplitude much stronger than predicted by standard theories.
Therefore the validity of these claims have been debated (e.g. in Hills et
al., 2018; Draine & Miralda-Escudé, 2018; Singh & Subrahmanyan, 2019; Bradley
et al., 2019), and unconventional physical processes such as an unknown
cooling mechanism (see e.g., Tashiro et al., 2014; Barkana, 2018; Fialkov et
al., 2018; Muñoz & Loeb, 2018; Ghara & Mellema, 2020) or the presence of an
additional radio background on top of the CMB (Feng & Holder, 2018; Ewall-Wice
et al., 2018; Fialkov & Barkana, 2019; Mebane et al., 2020) have been invoked
to explain these results.

While none of the other attempts have to date claimed a detection of the
signal, several of them have provided upper limits. A 2$\sigma$ value of
$\Delta^{2}(k)<(248)^{2}~{}\rm mK^{2}$ for a $k$-scale $0.5~{}h~{}\rm
Mpc^{-1}$ was the very first upper limit on the power spectrum888Here and for
the rest of this paper we use the so-called dimensionless form of the power
spectrum, $\Delta^{2}(k)$, which is related to the standard spherically
averaged power spectrum $P(k)$ through $\Delta^{2}(k)=k^{3}P(k)/2\pi^{2}$. of
the 21-cm signal from redshift $8.6$, obtained from observations with the GMRT
(Paciga et al., 2013). Barry et al. (2019), Cheng et al. (2018), and Kolopanis
et al. (2019) provided upper limits using the observations with MWA and PAPER.
Patil et al. (2017) published the first upper limit on the power spectrum of
the signal from redshift between 9.6 and 10.6 using LOFAR observations. Later,
Gehlot et al. (2019, 2020) and Eastwood et al. (2019) probed a higher redshift
range and also provided upper limits on the power spectrum using observations
with the LOFAR-Low Band Antenna array and the OVRO-LWA, respectively.
Recently, Mertens et al. (2020) analysed 10 nights of LOFAR observations of
the North Celestial Pole and obtained a 2$\sigma$ upper limit on the 21-cm
signal power spectrum at $z\approx 9.1$ of $\Delta^{2}(k=0.075~{}h~{}\rm
Mpc^{-1})=(73)^{2}\rm~{}mK^{2}$.

The focus of this paper is on the upper limits in the redshift range 6.5 – 8.7
published by the MWA team (Trott et al., 2020). Their best 2$\sigma$ upper
limit at redshift 6.5 is $\Delta^{2}(k=0.14~{}h~{}\rm
Mpc^{-1})\approx(43)^{2}\rm~{}mK^{2}$, derived from 110 hours of integration
with the MWA high band on the EoR0 field. The upper limit values increase with
redshift. At $z=8.7$, the upper limit of $\Delta^{2}(k=0.14~{}h~{}\rm
Mpc^{-1})\approx(250)^{2}\rm~{}mK^{2}$ comes from 51 hours of observations of
the EoR2 field.

The emitted 21-cm signal from H i depends on the number density of the sources
present at those epochs, as well as on the emissivity of UV, X-ray, radio and
$\rm{Ly{\alpha}}$ photons from those sources (see e.g., Ross et al., 2019;
Reis et al., 2020a). Due to the complex dependence of the signal on the
sources, it is not straightforward to extract the astrophysical as well as the
cosmological information from 21-cm signal observations, and exploration of
many theoretical models of the signal is required to interpret the results.
Often, the signal prediction algorithms are used in a Bayesian inference
framework to explore the parameter space of the reionization models (e.g.
Greig & Mesinger, 2015, 2017; Park et al., 2019; Cohen et al., 2020; Ghara et
al., 2020; Mondal et al., 2020; Greig et al., 2021a; Greig et al., 2021b).

In this paper, we consider the new upper limits on the 21-cm signal power
spectrum from MWA observations (Trott et al., 2020) in the redshift range
spanning from 6.5 to 8.7 to constrain the IGM conditions and explore scenarios
for the EoR that are disfavoured by such limits. Our framework is based on
simulations of the power spectrum and IGM parameters and a Bayesian inference
framework. Just as in our work on the LOFAR upper limits for $z\approx 9.1$
(Ghara et al., 2020), our main focus is on deriving constraints on the
condition of the IGM. In addition, we also aim to put constraints on the
additional radio background which has been proposed to explain the EDGES2
results. This work differs from the interpretation of the MWA upper limits by
Greig et al. (2021a), who instead focused on the properties of the sources and
did not consider any additional radio background component to the CMB.

This paper is structured as follows. In Section 2, we describe the basic
methodology of our framework. We present our results in Section 3. We discuss
the outcome of this study from the point of view of other observations in
Section 4, before concluding in Section 5. The cosmological parameters used
throughout this study are the same of the $N$-body simulations employed here,
i.e. $\Omega_{\rm m}=0.27$, $\Omega_{\Lambda}=0.73$, $\Omega_{\rm B}=0.044$,
$h=0.7$ (Wilkinson Microwave Anisotropy Probe (WMAP); Hinshaw et al., 2013).

##  2 Framework

Here we describe the framework employed to constrain the IGM properties using
the observational upper limits from the MWA.

Table 1: Overview of the simulation parameters used in grizzly and their explored ranges. Simulation Parameters | Description | Explored range |
---|---|---|---
$\zeta$ | Ionization efficiency |  [$10^{-3},10^{3}$] |
$M_{\rm min}$ | Minimum mass of the UV emitting halos |  [$10^{9}~{}{\rm M}_{\odot},10^{12}~{}{\rm M}_{\odot}$] |
$M_{\rm min,X}$ | Minimum mass of the X-ray emitting halos |  [$10^{9}~{}{\rm M}_{\odot},10^{12}~{}{\rm M}_{\odot}$] |
$f_{X}$ | X-ray heating efficiency |  [$10^{-3}$, $10^{2}$] |
$A_{r}$ | Efficiency of the excess radio background | [0, 416] |

###  2.1 Modelling the 21-cm signal using grizzly

Our framework uses the grizzly code (Ghara et al., 2015a; Ghara et al., 2018)
to generate 21-cm signal maps for the redshift range $z\approx 6.5-8.7$ for
different combinations of astrophysical parameters. The code requires gridded
versions of cosmological density fields and halo catalogues as input. We use
density fields of length 500 $h^{-1}$ comoving megaparsec (Mpc) smoothed onto
$300^{3}$ grids (see e.g, Giri et al., 2019a, b; Kamran et al., 2021), which
were retrieved from the results of the PRACE999Partnership for Advanced
Computing in Europe: http://www.prace-ri.eu/ project PRACE4LOFAR. This box
size corresponds to a field-of-view of $4.3^{\circ}\times 4.3^{\circ}$ at
redshift $\approx$ 9\. The minimum mass of the dark matter halos as identified
on the fly with a spherical overdensity halo finder (Watson et al., 2013) is
$\approx 10^{9}$ M⊙. These are the same cosmological $N$-body simulation
previously used in Ghara et al. (2020).

The grizzly code is an independent implementation of the previously developed
bears algorithm (Thomas & Zaroubi, 2008; Thomas et al., 2009; Thomas &
Zaroubi, 2011; Krause et al., 2018). It is based on a one-dimensional
radiative transfer scheme which approximates the transfer of photons (UV,
X-rays and $\rm{Ly{\alpha}}$) by assuming that the effect from individual
sources is isotropic. We refer the reader to the original papers (Ghara et
al., 2015a; Ghara et al., 2018) for a more detailed and complete description
of the algorithm.

Given a set of astrophysical source parameters, grizzly generates ionization
and gas temperature maps at different redshifts which are then used to produce
differential brightness temperature ($\delta T_{\rm b}$) maps. The $\delta
T_{\rm b}$ of the 21-cm signal can be expressed as (see e.g, Madau et al.,
1997; Furlanetto et al., 2006),

$\displaystyle\delta T_{\rm b}(\mathbf{x},z)\\!$ $\displaystyle=$
$\displaystyle\\!27~{}x_{\rm HI}(\mathbf{x},z)[1+\delta_{\rm
B}(\mathbf{x},z)]\left(\frac{\Omega_{\rm B}h^{2}}{0.023}\right)$
$\displaystyle\times$ $\displaystyle\\!\left(\frac{0.15}{\Omega_{\rm
m}h^{2}}\frac{1+z}{10}\right)^{1/2}\left(1-\frac{T_{\gamma,\rm eff}}{T_{\rm
S}(\mathbf{x},z)}\right)\,\rm{mK},$

where $x_{\rm HI}$ and $\delta_{\rm B}$ denote the neutral hydrogen fraction
and baryonic density contrast respectively, each at position $\mathbf{x}$ and
redshift $z$. The quantity $T_{\gamma,\rm eff}$ is the temperature of the
radio background, which reduces to the usual CMB temperature, $T_{\gamma}(z)$
= $2.725~{}(1+z)$ K, in the absence of additional contributions. We assume
that the spin temperature of hydrogen in the IGM, $T_{\rm S}$, is equal to the
gas temperature. This condition occurs in the presence of a strong
$\rm{Ly{\alpha}}$ background which is expected at the redshifts of interest
(e.g. Barkana & Loeb, 2005; Pritchard & Loeb, 2012; Schneider et al., 2020).

Afterwards we calculate the dimensionless power spectrum of $\delta T_{\rm
b}$, denoted as $\Delta^{2}(k)$. Note that we have ignored the presence of
redshift space distortions (RSD) while evaluating the power spectrum for
different model parameters, as their impact remains rather small during the
epoch of our interest101010The effects of RSD can enhance the large-scale
power spectrum by a factor of $\approx 2$ in case it is driven by density
fluctuations. In our case, we underestimate the power spectra values for those
scenarios. , when ionization fluctuations dominate the power spectrum of
$\delta T_{\rm b}$ (Jensen et al., 2013; Majumdar et al., 2016; Ghara et al.,
2015a, b; Ross et al., 2020).

###  2.2 Simulation parameters in grizzly

grizzly uses the following set of parameters to generate the maps of
ionization fraction and gas temperature in the IGM. These are later used to
generate $\delta T_{\rm b}$ maps and from these the power spectra. We describe
the parameters below and summarise the range of values explored in Table 1.

  * •

Ionization efficiency ($\zeta$): This scales the rate of emission of ionizing
photons per unit stellar mass from a halo $\dot{N}_{i}$ as
$\dot{N}_{i}=\zeta\times 2.85\times 10^{45}~{}{\rm s^{-1}}~{}{\rm
M}_{\odot}^{-1}$. $\zeta=1$ corresponds to a model galaxy spectrum which has
been produced with the publicly available population synthesis code pegase2
111111http://www2.iap.fr/pegase/ (Fioc & Rocca-Volmerange, 1997) and employed
in the 1D radiative transfer. We assume that the stellar mass $M_{\star}$ of a
halo is related to the dark matter halo mass $M_{\rm halo}$ as
$M_{\star}=f_{\star}~{}\frac{\Omega_{\rm B}}{\Omega_{\rm m}}~{}M_{\rm halo}$,
where the star formation efficiency $f_{\star}$ is fixed at 0.02 (Behroozi &
Silk, 2015; Sun & Furlanetto, 2016). Note that the effect of $\zeta$ is
degenerate with that of $f_{\star}$, the emission rate of ionizing photons
from the sources, and their escape fraction into the IGM. Here we vary
Log10($\zeta$) between [-3, 3].

  * •

Minimum mass of the UV emitting halos ($M_{\rm min}$): We assume no emission
of ionizing photons from halo with mass less than $M_{\rm min}$. This is to
incorporate the fact that various mechanism such as the radiative and
mechanical feedback can reduce the star formation efficiency significantly
(see e.g., Hasegawa & Semelin, 2013; Dawoodbhoy et al., 2018). In general, one
expects that the radiative feedback suppresses star formation in halos with
mass $\lesssim 10^{9}~{}{\rm M}_{\odot}$ (see e.g., Wise et al., 2014; Dixon
et al., 2016). In our study the lowest value for $M_{\rm min}$ is $10^{9}$
${\rm M}_{\odot}$. This is limited by the mass resolution of our N-body
simulation (see Sect. 2.1). We vary Log10($\frac{M_{\rm min}}{{\rm
M}_{\odot}}$) between [9, 12].

  * •

Minimum mass of X-ray emitting halo ($M_{\rm min,X}$): In addition to the UV
emitting sources, grizzly also includes heating and ionization from X-ray
emitting sources, such as quasars and high-mass X-ray binaries, as not
necessarily all halos hosting stellar sources also host X-ray sources. We
assume no emission of X-ray photons from halos with mass less than $M_{\rm
min,X}$, while all other halos host X-ray emitting sources. We vary
Log10($\frac{M_{\rm min,X}}{{\rm M}_{\odot}}$) between [9, 12].

  * •

X-ray heating efficiency ($f_{X}$): This scales the emission rate of X-ray
photons per unit stellar mass from a halo as $\dot{N}_{X}=f_{X}\times
10^{42}~{}\rm s^{-1}~{}{\rm M}_{\odot}^{-1}$. We model the spectrum of an
X-ray source as a power-law of the energy $E$, i.e. $I_{X}(E)\propto
E^{-\alpha}$, where we fix the spectral index $\alpha$ to 1.2. The value of
$\dot{N}_{X}$ for $f_{X}=1$ is consistent with the measurements of high-mass
X-ray binaries in local star forming galaxies in the 0.5-8 keV band (Mineo et
al., 2012; Islam et al., 2019). We assume that the X-ray band spans from 100
eV to 10 keV. We vary Log10($f_{X}$) between [-3, 2].

  * •

Radio background efficiency ($A_{r}$): This parameter quantifies the
contribution of a potential radio background existing in addition to the CMB.
For the LOFAR upper limits this parameter was not considered in Ghara et al.
(2020) but was instead studied in a separate paper by Mondal et al. (2020). We
model the effective brightness temperature of the radio background as (Fialkov
& Barkana, 2019; Mondal et al., 2020),

$T_{\gamma,\rm eff}=T_{\gamma}\left[1+A_{r}\left(\frac{\nu_{\rm obs}}{78~{}\rm
MHz}\right)^{-2.6}\right],$ (2)

where $\nu_{\rm obs}$ is the frequency of observation. $T_{\gamma,\rm eff}$
can differ from $T_{\gamma}$ in presence of a radio background of any
astrophysical or cosmological origin. Its existence is motivated by the
evidence of an excess radio background above the CMB towards the Rayleigh-
Jeans part of the spectrum, as detected by ARCADE2 (Fixsen et al., 2011) and
LWA1 (Dowell & Taylor, 2018). The LWA1 measurement of the excess radio
background at frequency 40-80 MHz fits well with a power-law with spectral
index -2.58$\pm$0.05 with a temperature of $603^{+102}_{-92}$ mK at 1.42 GHz.
Varying $T_{\gamma,\rm eff}$ changes the brightness temperature maps although
the ionization and gas temperature maps from grizzly remain unchanged. We vary
$A_{r}$ from 0 to 416, where the lower bound stands for no excess radio
background and the upper bound represents the upper limits of 0.603 K as
measured by LWA1. Note that the radio background used in this study is
uniform, while in general an inhomogeneous radio background may enhance the
power spectrum even more than what we estimate here (Reis et al., 2020b). A
study with inhomogeneous radio background is beyond the scope of this paper.

Table 2: An overview of the IGM parameters considered in this paper. Note that the ‘heated region’ denotes the regions in the IGM with gas temperature larger than $T_{\gamma,\rm eff}$.  IGM Parameters | Description |
---|---|---
$\overline{x_{\rm HII}}$ | Volume averaged ionized fraction |
$\overline{T_{\rm K}}$ (K) |  Volume averaged gas temperature in the partially ionized IGM with $x_{\rm HII}<0.5$ |
$\overline{\delta T_{\rm b}}$ (mK) |  Mass averaged differential brightness temperature |
$f_{\rm heat}$ |  Volume fraction of regions with temperature larger than $T_{\gamma,\rm eff}$ |
$R^{\rm heat}_{\rm peak}$ (Mpc/h) |  Size at which the PDF of the size distribution of the heated regions peaks |
$\Delta R^{\rm heat}_{\rm FWHM}$ (Mpc/h) |  FWHM of the PDF of the size distribution of the heated regions |

###  2.3 Derived IGM parameters

One should realise that the 21-cm signal observations do not contain any
direct information about the properties of the sources at the redshift of
interest. Instead, the 21-cm signal characterises the state of the IGM at
those redshifts. Thus, the main focus of this study is on the IGM properties
rather than on the source properties, i.e. we aim to constrain the physical
properties of IGM, such as the ionization and thermal state, rather than the
astrophysical simulation parameters used to simulate the brightness
temperature maps. Below we list the IGM parameters that we consider here at
each redshift (see also in Table 2):

  * •

$\overline{x_{\rm HII}}$: Volume averaged ionized fraction of the IGM.

  * •

$\overline{T_{\rm K}}$ (K): Volume averaged gas temperature of the regions in
the IGM with ionization fraction $x_{\rm HII}<0.5$.

  * •

$\overline{\delta T_{\rm b}}$ (mK): Mass averaged differential brightness
temperature in the IGM.

  * •

$f_{\rm heat}$: Volume fraction of the IGM with gas temperature $T_{\rm
K}>T_{\gamma,\rm eff}$.

  * •

$R^{\rm heat}_{\rm peak}$ (cMpc): Size of the ‘heated regions’, i.e. regions
with $T_{\rm K}>T_{\gamma,\rm eff}$ at which the probability distribution
function (PDF) of the sizes peaks. To determine the size distribution of the
heated regions we use the mean free path (MFP) method, which is a Monte Carlo
based approach (Mesinger & Furlanetto, 2007; Giri et al., 2018a).

  * •

$\Delta R^{\rm heat}_{\rm FWHM}$ (cMpc): Full width at half maximum (FWHM) of
the PDF of the sizes of the heated regions.

Unlike in Ghara et al. (2020), here we have not considered separate IGM
parameters for the H ii regions. In the absence of X-ray heating, the ‘heated
regions’ in the IGM are equivalent to the H ii regions, as the kinetic
temperature of the gas within such regions is $\sim 10^{4}$K. Thus, for
insignificant heating in the neutral IGM, $f_{\rm heat}$ is equivalent to
$\overline{x_{\rm HII}}$, while $R^{\rm heat}_{\rm peak}$ and $\Delta R^{\rm
heat}_{\rm FWHM}$ characterise the size distribution of the H ii regions. Note
that we estimate these quantities from grizzly simulations at our simulation
resolution.

###  2.4 Bayesian inference framework

We generate $10^{5}$ power spectra and sets of IGM parameters at each redshift
point using grizzly, while employing the gridded sampling method to sample the
5D parameter space. We then use a linear interpolation scheme121212Note that
in our previous study (Ghara et al., 2020) we used an emulator based on the
Gaussian Process Regression (GPR) algorithm. However, GPR does not scale well
with the size of the training set (e.g. Rasmussen & Williams, 2006). As the
parameter space used here has a higher dimension than the one explored in
Ghara et al. (2020), a larger training set is needed. Tests showed that this
would lead the emulator to achieve a size of several tens of gigabytes,
slowing it down too much to be useful. on a regular grid to determine the
power spectrum for a set of parameters $\theta$. To avoid extrapolation, the
range of $\theta$ used is the same as the range of the gridded parameter space
covered by the grizzly simulations performed for this study. We use the
regular grid interpolation method implemented in scipy python package (Weiser
& Zarantonello, 1988; Virtanen et al., 2020). See Monsalve et al. (2018) for a
similar Bayesian framework that is based on interpolation in a regular grid of
4D parameter space. This linear interpolator is used in the emcee python MCMC
module (Foreman-Mackey et al., 2013) to explore the parameter space as
described in section 2.2. Note that we use only the smallest three $k$-bin
values from Trott et al. (2020), i.e. 0.142, 0.212 and 0.283 $h~{}\rm
Mpc^{-1}$, as the upper limits at other wavenumbers are too high to give any
additional constraints.

Table 3: The $\Delta^{2}_{\mathrm{up}}\equiv\Delta^{2}_{\mathrm{21}}+2\sigma_{21,\rm err}$ upper limits from MWA observations at different redshifts for different $k$-bins. See Trott et al. (2020) for details. We obtained $\Delta^{2}_{\mathrm{21}}$ and the corresponding error from private communication with the MWA team. k ($h~{}\rm Mpc^{-1}$) |  $\Delta^{2}_{\mathrm{up}}(k)$ (mK2) at $z$= 6.5 |  $\Delta^{2}_{\mathrm{up}}(k)$ (mK2) at $z$= 6.8 |  $\Delta^{2}_{\mathrm{up}}(k)$ (mK2) at $z$= 7.1 |  $\Delta^{2}_{\mathrm{up}}(k)$ (mK2) at $z$= 7.8 |  $\Delta^{2}_{\mathrm{up}}(k)$ (mK2) at $z$= 8.2 |  $\Delta^{2}_{\mathrm{up}}(k)$ (mK2) at $z$= 8.7 |
---|---|---|---|---|---|---|---
0.142 | $43.1^{2}$ | $60.1^{2}$ | $77.7^{2}$ | $154.2^{2}$ | $167.7^{2}$ | $249.6^{2}$ |
0.212 | $70.2^{2}$ | $90.0^{2}$ | $117.4^{2}$ | $247.5^{2}$ | $430.3^{2}$ | $569.9^{2}$ |
0.283 | $93.3^{2}$ | $114.1^{2}$ | $152.3^{2}$ | $314.5^{2}$ | $422.2^{2}$ | $562.5^{2}$ |
0.354 | $209.5^{2}$ | $243.9^{2}$ | $281.5^{2}$ | $460.1^{2}$ | $540.9^{2}$ | $688.1^{2}$ |
0.425 | $183.5^{2}$ | $221.3^{2}$ | $263.3^{2}$ | $804.4^{2}$ | $772.8^{2}$ | $963.2^{2}$ |
0.495 | $125.5^{2}$ | $169.0^{2}$ | $231.9^{2}$ | $466.8^{2}$ | $1402.6^{2}$ | $1854.5^{2}$ |
0.566 | $210.1^{2}$ | $255.4^{2}$ | $310.9^{2}$ | $484.4^{2}$ | $1109.9^{2}$ | $1546.0^{2}$ |
0.637 | $214.1^{2}$ | $260.3^{2}$ | $333.8^{2}$ | $501.0^{2}$ | $739.1^{2}$ | $962.3^{2}$ |
0.708 | $384.6^{2}$ | $383.1^{2}$ | $437.9^{2}$ | $613.4^{2}$ | $781.1^{2}$ | $947.6^{2}$ |

The MCMC analysis requires the calculation of a likelihood. As we are working
with upper limits we cannot use the formal $\chi^{2}$ which only applies for
actual measurements of a signal. Also, as the upper limits have quite high
values, we want our framework to find the models which are excluded, rather
than the ones which are allowed. Given a measured upper limit of the power
spectrum at redshift $z$ as
$\Delta^{2}_{\mathrm{21}}(k_{i},z)\pm\Delta^{2}_{\mathrm{21,err}}(k_{i},z)$,
the likelihood of a set of parameters $\bm{\theta}$ to be excluded is defined
as (see Ghara et al., 2020, for details)

$\mathcal{L}_{\rm
ex}(\theta,z)=1-\prod_{i}\frac{1}{2}\left[1+\mathrm{erf}\left(\frac{\Delta^{2}_{\mathrm{21}}(k_{i},z)-\Delta^{2}_{\mathrm{m}}(k_{i},\bm{\theta},z)}{\sqrt{2}\sigma(k_{i},z)}\right)\right],$
(3)

where $k_{i}$ denotes the 3 $k$-bins of the measured power spectrum, and
$\Delta^{2}_{\mathrm{m}}(k,\bm{\theta},z)$ the model power spectrum for a set
of parameters $\theta$ for $k$-scale at redshift $z$. The quantity $\sigma$
accounts for errors from the measurement, as well as the modelling error, i.e.
$\sigma^{2}=\Delta^{4}_{\mathrm{21,err}}+\Delta^{4}_{\mathrm{m,err}}$. We
adopt 30 per cent modelling error
$\Delta^{2}_{\mathrm{m,err}}(k_{i},z)=0.3\times\Delta^{2}_{\mathrm{m}}(k_{i},\bm{\theta},z)$
to incorporate errors from the modelling with the code grizzly and from the
interpolation scheme. In Ghara et al. (2018) we performed a detailed
comparison of grizzly with the 3D radiative transfer scheme C2ray, and found
that the error from grizzly is less than 10 per cent. We estimate the error
due to the interpolation scheme on a test set of the power spectra and found
it to be less than 30 per cent for the scales of interest. Note that the
sample variances on the power spectrum at the scales of interest are
negligible compared to the adopted modelling errors.

Figure 1: A subset of the power spectra of the 21-cm signal brightness
temperature generated using grizzly over the gridded parameter space. These
correspond to $A_{r}=0$ or, equivalently, $T_{\gamma,\rm
eff}=2.725\times(1+z)$ K. Note that we have also varied $A_{r}$ to get the
full power spectra set. Left and right panels correspond to $z\approx 6.5$ and
$z\approx 7.8$, respectively. The red points refer to the recent MWA 2$\sigma$
upper limits from Trott et al. (2020).

Unlike Greig et al. (2021a), we consider the upper limits at each redshift
independently to obtain constraints on the simulation parameters and
thereafter on the IGM parameters. One requires a background source model to
derive constraints on the simulation parameters using multi-redshift
constraints. Thus, the derived constraints on the IGM parameters at individual
redshifts from such a combined upper limits analysis would be source model
dependent. Moreover, those constraints would be biased as the likelihood in
that case would be mostly dominated by the best upper limit at $z\approx 6.5$
(see Table 3). Thus, to avoid dependence of the quantities of our interest on
the background source model, we consider the upper limits at each redshift
independently.

Once the MCMC chains that contain the lists of the simulation parameters used
in the MCMC analysis is built, we generate the list of corresponding IGM
parameters using the interpolation scheme. Finally, these are used to
determine the posterior distribution of the IGM parameters.

##  3 Results

In this section, we apply the Bayesian inference framework described in the
previous section to the measured upper limits from Trott et al. (2020). These
2$\sigma$ upper limits are also shown in Table 3. Clearly, the upper limits
are stronger at the larger (spatial) scales (lower ‘$k$’ values) and lower
redshifts. Figure 1 shows the set of the power spectra generated with grizzly
for $A_{r}=0$ at redshift 6.5 (left panel) and 7.8 (right panel),
respectively. The red points in the panels denote the measured MWA 2$\sigma$
upper limits on the power spectrum. Clearly, the power spectrum amplitude of a
significant fraction of the models is above the upper limit at scales
$\lesssim 0.2~{}h~{}\rm Mpc^{-1}$ at $z=6.5$, even in the absence of an excess
radio background. This is also true for redshifts 6.8 and 7.2. On the other
hand, the right panel shows that the upper limits at $z=7.8$ (and also for
$z=8.2$ and 8.7) remain above the power spectra values from grizzly for
$A_{r}=0$.

As explained in Section 2.4 the MCMC provides the probability of the models to
be excluded by the MWA upper limits. Each MCMC run consists of 10 walkers and
$10^{6}$ steps. We find that the MCMC chains converge well before the final
step. The explored parameters ranges are listed in Table 1. Note that all the
priors are flat or uniform in log-scale.

Table 4: Constraints on the excluded part of the IGM parameter space at different redshifts for the scenario where $A_{r}=0$, i.e. no excess radio background is considered in addition to the CMB. These are obtained from the MCMC analysis using the multi-redshift MWA upper limits from Trott et al. (2020). The priors are derived from the grizzly simulations. IGM Parameters |  Credible intervals | $z=6.5$ | $z=6.8$ | $z=7.2$ |
---|---|---|---|---|---
$\overline{x_{\rm HII}}$ | Prior | [0, 1] | [0, 1] | [0, 1] |
68$\%$ | [0, 0.15] | [0, 0.15] | [0, 0.26] |
|  95$\%$ | [0, 0.56] | [0, 0.53] | [0, 0.56] |
$\overline{T_{\rm K}}$ (K) | Prior |  [1.17, $10^{4}$] |  [1.28, $10^{4}$] |  [1.38, $10^{4}$] |
68$\%$ | [1.17, 40.1] | [1.28, 11] | [1.38, 15.7] |
|  95$\%$ | [1.17, 251] | [1.28, 158] | [1.38, 90] |
$f_{\rm heat}$ | Prior | [0, 1] | [0, 1] | [0, 1] |
68$\%$ | [0, 0.28] | [0, 0.18] | [0, 0.28] |
|  95$\%$ | [0, 0.63] | [0, 0.55] | [0, 0.58] |
$\overline{\delta T_{\rm b}}$ (mK) | Prior | [-376, 23.3) | [-367, 23.8) | [-360, 24.4) |
68$\%$ | [-376, -140] | [-367,-247] | [-360, -229] |
|  95$\%$ | [-376., -89] | [-367, -125] | [-360, -121] |
$R^{\rm heat}_{\rm peak}$ ($~{}h^{-1}\rm Mpc$) | Prior | [0, 198] | [0, 203] | [0, 201] |
68$\%$ | [0, 4.7] | [0, 8.5] | [0, 18.2] |
|  95$\%$ | [0, 18.] | [0, 17.5] | [0, 30] |
$\Delta R^{\rm heat}_{\rm FWHM}$ ($~{}h^{-1}\rm Mpc$) | Prior | [0, 440] | [0, 445] | [0, 452] |
68$\%$ | [0, 12.5] | [0, 16.6] | [0,35] |
|  95$\%$ | [0, 47] | [0,39] | [0, 63] |

To constrain the IGM parameters using our inference framework we consider a
case with $A_{r}=0$, i.e. no excess radio background contribution to the CMB,
and one in which we also vary the $A_{r}$ parameter. We will first explore the
scenario without any excess radio background in Section 3.1, while in Section
3.2 we present the constraints allowing $A_{r}$ parameter to vary.

Figure 2: Posterior distribution of the simulation parameters at $z\approx
6.5$ from the MCMC analysis. Here $A_{r}=0$. The color-bar shows the
marginalized probability of the models to be ruled out. The solid and dashed
curves corresponds to the 68 and 95 per cent credible intervals of the models
ruled out by the MWA upper limit at that redshift. The diagonal panels refer
to the marginalized probability distributions for ruled out simulation
parameters.  Figure 3: Posterior distribution of the IGM parameters at
$z\approx 6.5$ from the MCMC analysis. Here $A_{r}=0$. The color-bar shows the
marginalized probability of the models to be ruled out. The solid and dashed
curves corresponds to the 68 and 95 per cent credible intervals of the models
ruled out by the MWA upper limit at that redshift. The diagonal panels refer
to the marginalized probability distributions for ruled out IGM parameters.

###  3.1 Analysis without excess radio background

First, we explore the four-dimensional source parameter space, i.e. $\zeta$,
$M_{\rm min}$, $M_{\rm min,X}$ and $f_{X}$, using MCMC at six different
redshifts. We fix $A_{r}=0$, i.e. we initially assume that the radio
background is contributed by the CMB alone.

The posterior distribution of the source parameter space of the excluded
models at $z=6.5$, our fiducial redshift, is shown in Figure 2. Three main
scenarios can be envisioned that produce a high amplitude of the large-scale
power spectrum. In the first scenario the large-scale power spectrum is
dominated by $x_{\rm HII}$ fluctuations, which are enhanced by a distribution
of rare but large H ii regions in the IGM. Note that the large-scale power-
spectrum decreases for a highly ionized IGM with $\overline{x_{\rm
HII}}\gtrsim 0.5$, and almost vanishes for $\overline{x_{\rm HII}}\sim 1$.
This suggests that high values of $\zeta$ and low values of $M_{\rm min}$
should be disfavored among the excluded models, as confirmed in the figure.
From this, we deduce that $\zeta\lesssim 2.4$ and $M_{\rm min}\gtrsim 3\times
10^{10}~{}{\rm M}_{\odot}$ are the most likely to be excluded within the 68
per cent credible intervals by the MWA upper limit at redshift 6.5 alone. Note
also that a colder IGM, which is obtained with a smaller value of $f_{X}$,
will also enhance the large-scale power spectrum in this scenario. An
alternative scenario that produces a high amplitude of the large-scale power
spectrum is when the latter is dominated by $T_{\rm S}$/$T_{\rm K}$
fluctuations. This scenario can happen for a small number density of X-ray
emitting sources, as numerous X-ray sources and high X-ray emissivity can
easily heat the IGM above $T_{\gamma}$ and suppress the large-scale power
spectrum. Indeed, $f_{X}\lesssim 0.1$ and $M_{\rm min,X}\gtrsim 10^{11}~{}{\rm
M}_{\odot}$ are the most likely to be excluded within the 68 per cent credible
intervals by the MWA upper limit at redshift 6.5 alone. Alternatively, a
scenario with neutral and cold IGM (as indicated by low values of $\zeta$ and
$f_{X}$) where the large-scale power spectrum is dominated mainly by the
density fluctuations is also important.

While the parameters the $M_{\rm min}$, $M_{\rm min,X}$ and $f_{X}$ show a
smooth distribution over the explored ranges, one can see a distinct peak
around $\zeta\approx 1$ in the posterior distribution of $\zeta$ in Figure 2.
The combination of $\zeta\approx 1$ and a high value of $M_{\rm min}\gtrsim
10^{11}~{}{\rm M}_{\odot}$ maximises the contribution from ionisation
fluctuations. These ionization states correspond to $\overline{x_{\rm
HII}}\sim 0.5$. On the other hand, the other two scenarios, i.e. density
fluctuations in a cold and neutral IGM and $T_{\rm S}$ fluctuations in a
neutral IGM, maximise the PDF for $\overline{x_{\rm HII}}\sim 0$, which
requires small values for $\zeta$. This bimodal feature in the PDF of $\zeta$
is consistent with the results of Ghara et al. (2020), where we perform a
similar analysis for LOFAR upper limits at redshift 9.1. Note that we
considered a prior on the ionization fraction at redshift 9.1 in Ghara et al.
(2020).

Table 5: Constraints on the excluded part of the IGM parameter space at different redshifts for a varying radio background scenario. These are obtained from the MCMC analysis using the multi-redshift upper limits from MWA (Trott et al., 2020). The priors are derived from the simulations run with grizzly. IGM Parameters |  Credible intervals | $z=6.5$ | $z=6.8$ | $z=7.2$ | $z=7.8$ | $z=8.2$ | $z=8.7$ |
---|---|---|---|---|---|---|---|---
Log10($A_{r}$) | Prior | [-2, 2.6 ] | [-2, 2.6 ] | [-2, 2.6 ] | [-2, 2.6 ] | [-2, 2.6 ] | [-2, 2.6 ] |
68$\%$ | [0.4, 2.6 ] | [1, 2.6 ] | [1.4, 2.6 ] | [1.87, 2.6 ] | [1.88, 2.6 ] | [1.9, 2.6 ] |
|  95$\%$ | [-0.8, 2.6 ] | [-0.76, 2.6 ] | [-0.65, 2.6 ] | [1.36, 2.6 ] | [1.38, 2.6 ] | [1.43, 2.6 ] |
$\overline{x_{\rm HII}}$ | Prior | [0, 1] | [0, 1] | [0, 1] | [0, 1] | [0, 1] | [0, 1] |
68$\%$ | [0, 0.15] | [0, 0.15] | [0, 0.17] | [0, 0.12] | [0, 0.15] | [0, 0.12] |
|  95$\%$ | [0, 0.58] | [0, 0.53] | [0, 0.56] | [0, 0.56] | [0, 0.55] | [0, 0.52] |
$\overline{T_{\rm K}}$ (K) | Prior |  [1.17, $10^{4}$] |  [1.28, $10^{4}$] |  [1.38, $10^{4}$] |  [1.63, $10^{4}$] |  [1.75, $10^{4}$] |  [2, $10^{4}$] |
68$\%$ | [1.17, 208] | [1.28, 80.7] | [1.38, 69] | [1.63, 69] | [1.75, 68] | [2, 42] |
|  95$\%$ | [1.17, 1025] | [1.28, 795] | [1.38, 790] | [1.63, 602] | [1.75, 592] | [2, 400.] |
$f_{\rm heat}$ | Prior | [0, 1] | [0, 1] | [0, 1] | [0, 1] | [0, 1] | [0, 1] |
68$\%$ | [0, 0.33] | [0, 0.28] | [0, 0.28] | [0, 0.28] | [0, 0.23] | [0, 0.18] |
|  95$\%$ | [0, 0.64] | [0, 0.6] | [0, 0.59] | [0, 0.56] | [0, 0.57] | [0, 0.55] |
Log10($\frac{|\overline{\delta T_{\rm b}}|}{\rm mK}$) | Prior | [1.37, 4.6] | [1.38, 4.59] | [1.39, 4.58] | [1.4, 4.57] | [1.41, 4.56] | [1.43, 4.54] |
68$\%$ | [2.18, 3.3] | [2.4, 3.57] | [2.65, 3.76] | [3.21, 4] | [3.24, 4.04] | [3.38, 4.06] |
|  95$\%$ | [1.96, 4.2] | [2.18, 4.24] | [2.29, 4.27] | [2.94, 4.39] | [2.98, 4.38] | [3.19, 4.43] |
$R^{\rm heat}_{\rm peak}$ ($~{}h^{-1}\rm Mpc$) | Prior | [0, 198] | [0, 203] | [0, 201] | [0, 201] | [0, 216] | [0, 198] |
68$\%$ | [0, 4.] | [0, 3.3] | [0, 3] | [0,3.6] | [0,3.6] | [0,3.5] |
|  95$\%$ | [0, 12.] | [0, 10] | [0, 13] | [0,12] | [0,11.5] | [0, 11.5] |
$\Delta R^{\rm heat}_{\rm FWHM}$ ($~{}h^{-1}\rm Mpc$) | Prior | [0, 440] | [0, 445] | [0, 452] | [0, 452] | [0, 447] | [0, 452] |
68$\%$ | [0, 14] | [0, 9.1] | [0, 11.4] | [0, 10.] | [0, 11.2] | [0, 9.3] |
|  95$\%$ | [0, 29] | [0, 22] | [0, 26] | [0, 23] | [0, 25] | [0, 23] |

Figure 3 shows the posterior distributions of the IGM parameters at redshift
6.5. The constraints on such parameters are also listed in Table 4. We find
that an IGM with $\overline{x_{\rm HII}}\lesssim 0.56$, $\overline{T_{\rm
K}}\lesssim 250$ K, $f_{\rm heat}\lesssim 0.63$, $R^{\rm heat}_{\rm
peak}\lesssim 18~{}h^{-1}$ Mpc and $\Delta R^{\rm heat}_{\rm FWHM}\lesssim
50~{}h^{-1}$ Mpc falls within 95 per cent credible intervals of the ruled out
IGM scenarios. The same limit for the average brightness temperature is
$\overline{\delta T_{\rm b}}\lesssim-90$ mK. The PDFs of $\overline{x_{\rm
HII}}$, $f_{\rm heat}$ and $\overline{\delta T_{\rm b}}$ suggest that the set
of excluded models contains a significant number of models in which the IGM is
still neutral and cold at redshift $z\approx 6.5$. A similar conclusion was
reached by Greig et al. (2021a). However, it is hard to tell from the achieved
limits on the IGM parameters if these neutral models are completely cold or
contain patchy heated regions in an otherwise cold IGM.

The 68 and 95 per cent credible intervals limits of the excluded IGM
parameters for other redshifts are also listed in Table 4. As the upper limit
of the power spectrum from MWA is weaker at a higher redshift, the limits also
become weaker, and especially for $z\gtrsim 7.8$ their significance is very
small. In fact, we neither get any constraint on the source nor on the IGM
parameters at $z\approx 8.7$. For $z\approx 7.8$ and 8.2, the values of
log10($\mathcal{L}_{\rm ex}$) remain $\lesssim-5$, which corresponds to the
value of log10($\mathcal{L}_{\rm ex}$) when a modelled power spectrum is below
3$\sigma$ error of the observed upper limits at all three scales. Clearly, the
reason is that all the grizzly models have power spectra values that are well
below the upper limits (see Figure 1). As the significance of the limits at
redshifts $\gtrsim 7.8$ is poor, we do not show them in this work. Note that
the limits as stated in Table 4 only stand for the part of the parameter space
which has the largest probability to be ruled out.

###  3.2 Analysis with varying radio background

Now, we explore the five-dimensional simulation parameter space, i.e. $\zeta$,
$M_{\rm min}$, $M_{\rm min,X}$, $f_{X}$ and $A_{r}$ using MCMC at six
different redshifts. The explored ranges of the simulation parameters are
listed in Table 1.

Figure 4 shows the posterior distribution of the simulation parameters of the
models with high exclusion probability at $z\approx 6.5$. The 68 per cent
credible intervals limits of these model parameters are $\zeta\lesssim 1$,
$M_{\rm min}\gtrsim 2\times 10^{10}~{}{\rm M}_{\odot}$, $M_{\rm min,X}\gtrsim
4\times 10^{10}~{}{\rm M}_{\odot}$ and $f_{X}\lesssim 0.1$. As the amplitude
of the power spectrum increases with increasing $A_{r}$, the excluded models
prefer a higher value of $A_{r}$. From this five-dimensional MCMC analysis, in
fact, we find that $A_{r}\gtrsim 0.15$ ($\gtrsim 2.5$) is the most likely to
be excluded within the 95 (68) per cent credible intervals by the MWA upper
limit at redshift 6.5 alone.

Figure 4: Posterior distribution of the simulation parameters at $z\approx
6.5$ from the MCMC analysis. The color-bar shows the marginalized probability
of the models to be ruled out. The solid and dashed curves corresponds to the
68 and 95 per cent credible intervals of the models ruled out by the MWA upper
limit at that redshift. The diagonal panels refer to the marginalized
probability distributions for ruled out simulation parameters.  Figure 5:
Posterior distribution of the IGM parameters at $z\approx 6.5$ from the MCMC
analysis. The color-bar shows the marginalized probability of the models to be
ruled out. The solid and dashed curves corresponds to the 68 and 95 per cent
credible intervals of the models ruled out by the MWA upper limit at that
redshift. The diagonal panels refer to the marginalized probability
distributions for ruled out IGM parameters. Note that we have used Logarithm
of the absolute values of $\overline{\delta T_{\rm b}}$ as the range is large.
In reality, $\overline{\delta T_{\rm b}}$ values of these excluded models are
mostly negative.  Figure 6: limits on the IGM parameters as obtained from the
excluded models at six different redshifts using the recent MWA upper limits
within 68 and 95 per cent credible interval. These constraints have been
determined at each redshift independently. For $z\geq 7.8$, no constraints are
found in absence of an additional radio background ($A_{r}=0)$.

Figure 5 shows instead the posterior distributions of the IGM parameters at
the same redshift. These constraints are also listed in Table 5 and show that
an IGM with $\overline{x_{\rm HII}}\lesssim 0.6$, $\overline{T_{\rm
K}}\lesssim 10^{3}$ K, $f_{\rm heat}\lesssim 0.65$, $R^{\rm heat}_{\rm
peak}\lesssim 13~{}h^{-1}$ Mpc and $\Delta R^{\rm heat}_{\rm FWHM}\lesssim
30~{}h^{-1}$ Mpc falls within 95 per cent credible intervals of the ruled out
IGM scenarios. As these intervals indicate a larger value of $A_{r}$ (i.e.
$A_{r}\gtrsim 2.5$), the derived limit of the average brightness temperature,
$\overline{\delta T_{\rm b}}\gtrsim-10^{4}$ mK, is weaker compared to what we
had found for $A_{r}=0$.

The constraints on the IGM parameters for the other five redshifts are also
listed in Table 5 and presented in Fig. 6. We find that the 95 per cent
credible intervals for each of the quantities, $\overline{x_{\rm HII}}$,
$f_{\rm heat}$, $R^{\rm heat}_{\rm peak}$ and $\Delta R^{\rm heat}_{\rm
FWHM}$, are similar at all six redshifts. This shows that the ionization and
thermal states of the IGM which are given by the IGM parameter values within
these limits produce extreme amplitudes of the large-scale power spectra,
while the limits on the IGM parameters are independent of redshift. Note that
because the source populations differ between redshifts, it will require
different values of the source parameters to achieve similar thermal and
ionization state of the IGM at different redshifts. As the power spectrum
amplitude increases with increasing $A_{r}$ and the upper limits are weaker at
higher $z$, the 95 per cent credible interval limit for $A_{r}$ increases as
well, for example from $\approx 0.15$ at redshift 6.5 to $\approx 27$ at
redshift 8.7. These correspond to $T_{\gamma,\rm eff}\approx 26$ K and 140 K
at redshifts 6.5 and 8.7, respectively. The limit at redshift 6.5 is
equivalent to an excess radio background which is $0.008\%$ of the CMB at 1.42
GHz. Excess of $0.009\%,0.01\%,1.2\%,1.3\%$ and $1.5\%$ are found from the
limits at redshifts 6.8, 7.2, 7.8, 8.2 and 8.7, respectively.

##  4 Discussion

We highlight that the constraints on the IGM and the simulation parameters
obtained in the previous section do not correspond to models ruled out with a
high significance, but rather to models that have a high probability to be
excluded. Future stronger upper limits will improve the
significance/probability of the models to be ruled out.

It should be noted that the part of the parameter space which has high
probability to be excluded is associated with models that produce a high
amplitude of the large-scale power spectrum. Large ionized regions in a cold
IGM, and/or large heated regions in a non-uniform spin temperature IGM are the
obvious scenarios that come to mind. Alternatively, density fluctuations in a
cold IGM can also produce high values of the large-scale power spectrum,
especially at the lower redshifts where the density fluctuations are larger
and the IGM is expected to be colder in the absence of heating sources. The
same models also suggest a high amplitude of the excess radio background, as
this boosts the fluctuations at all scales. It is challenging to devise
scenarios alternative to those discussed here that have such a high chance of
exclusion without accounting for other non-standard physics, such as e.g. an
excess cooling mechanism that makes the gas in the IGM colder than in a
standard scenario.

As we consider each redshift independently, the models more likely to be
excluded show similar IGM conditions (e.g. $\overline{x_{\rm HII}}\lesssim
0.5$ and cold IGM) at all redshifts, although the constraints become obviously
stronger at redshifts with tighter upper limits. In our previous study (Ghara
et al., 2020), we obtained constraints on the IGM parameters at $z\approx 9.1$
using the upper limits from LOFAR (Mertens et al., 2020), finding that an IGM
with $f_{\rm heat}\lesssim 0.35$, $T_{\rm K}\lesssim 160$ K, $R^{\rm
heat}_{\rm peak}\lesssim 70~{}h^{-1}\rm Mpc$ and $\Delta R^{\rm heat}_{\rm
FWHM}\lesssim 110~{}h^{-1}\rm Mpc$ is the most likely to be ruled out within
95 per cent credible interval limits. This physical state of the IGM is
similar to those found to be more likely ruled out in the present study. This
is to be expected, as those states of the IGM which produce high amplitudes of
the large-scale power spectra are easy to exclude by the 21-cm signal upper
limits.

We note that, for some of the models with high exclusion probability, the
large-scale power spectrum is driven by $T_{\rm S}$ fluctuations (and/or
density fluctuations) only when the $x_{\rm HII}$ fluctuations are small.
These models have an ionization fraction $\overline{x_{\rm HII}}\approx 0$ at
$z\approx 6.5$, which is inconsistent with the combined constraints from
observations of the Thomson scattering optical depth and $z>6$ quasar spectra
(see e.g., Mitra et al., 2015), as well as $\rm{Ly{\alpha}}$ emitter
statistics (Ouchi et al., 2010), spectroscopy of gamma-ray bursts (Chornock et
al., 2013), dark pixel statistics of the spectroscopy of high-$z$ quasars
(McGreer et al., 2015), which all indicate that the IGM is not highly neutral
at redshift $\lesssim 7$.

The constraints presented in the previous section are obtained analysing only
one of the six upper limits/redshifts at a time. In the MCMC analysis the
upper limits are not all combined into a single likelihood, as otherwise this
would be mostly driven by the best upper limit, i.e. the one at $z\approx
6.5$, biasing the conclusions on the IGM properties at individual redshifts.

A recent study by Greig et al. (2021a) obtained constraints on source
parameters using the MWA upper limits by combining all redshifts into a single
likelihood for a set of redshift independent source parameters. The derived
constraints on the IGM parameters at different redshifts as obtained in fig 6
of Greig et al. (2021a) depend on the background source model. As a result,
the constraints on $\overline{x_{\rm HI}}$ follow a redshift evolution as
expected in such model. We on the other hand constrain the IGM state at each
redshift using the upper limits of the power spectra at that redshift only,
and thus the constraints do not depend on the redshift evolution in the
background model. The differences in the analysis methods make it hard to
directly compare our results with those in Greig et al. (2021a). However, as
the results of Greig et al. (2021a) are expected to be biased by the strongest
upper limits at redshift 6.5, one might expect that the IGM states of the
excluded models of both studies are qualitatively similar at redshift 6.5,
which we indeed have found to be the case in Section 3.1.

Using LOFAR upper limits at $z=9.1$ from Mertens et al. (2020), Mondal et al.
(2020) obtained constraints on the excess radio background as $0.1-9.6\%$ of
the CMB at 1.42 GHz. Note that these 95 per cent credible limits are obtained
from marginalizing models which are valid under the LOFAR upper limits. Our
limits are instead derived from the marginalization of the excluded models. It
is therefore not possible to compare the above constraint to the one we find,
namely that the 95 per cent credible limits of $A_{r}$ at redshift 8.7 (6.5)
correspond to an excess radio background which is $1.5\%$ ($0.008\%$) of the
CMB at 1.42 GHz.

##  5 Summary & Conclusions

Trott et al. (2020) has provided new upper limits on the dimensionless
spherically averaged power spectrum of the 21-cm signal from MWA at six
different redshifts, in the range $6.5-8.7$. The best 2$\sigma$ upper limit
from that paper, obtained from 110 hours of integration time, is
$\Delta^{2}(k=0.14~{}h~{}\rm Mpc^{-1})=(43)^{2}$ mK2 at $z\approx 6.5$. Here
we have used these limits to rule out possible reionization scenarios. The
main focus of this study is to constrain the ionization and thermal states of
the IGM.

We use the grizzly code to generate hundred of thousands of models for
different combinations of parameters, namely, ionization efficiency ($\zeta$),
minimum mass of the UV emitting halos ($M_{\rm min}$), minimum mass of X-ray
emitting halos ($M_{\rm min,X}$), X-ray heating efficiency ($f_{X}$) and
excess radio background efficiency ($A_{r}$). The outputs of these models are
power spectra and derived IGM parameters, such as volume averaged ionization
fraction ($\overline{x_{\rm HII}}$), volume averaged gas temperature of the
partially ionized IGM ($\overline{T_{\rm K}}$), mass averaged brightness
temperature ($\overline{\delta T_{\rm b}}$), volume fraction of the heated
region ($f_{\rm heat}$), size of the heated regions at which the PDF of the
sizes peaks ($R^{\rm heat}_{\rm peak})$ and the FWHM of the PDFs ($\Delta
R^{\rm heat}_{\rm FWHM})$. The outputs from the grizzly simulations are used
within an MCMC framework to constrain the simulation parameters using the
upper limits from MWA. The simulation parameters from the MCMC chains are
later used to generate lists of IGM parameters. The main findings of the paper
are the followings.

  * •

From the set of power spectra generated using grizzly, we do not find any
model which can be excluded in the absence of an extra radio background
component at redshifts 7.8, 8.2 and 8.7. On the other hand, we find a
significant number of models that can be ruled out at the other three
redshifts.

  * •

From the MCMC analysis, we find that the excluded models are more likely to
have a smaller number density of the UV as well as X-ray emitting sources. At
redshift $\approx 6.5$, for 68 per cent credible region of the excluded
models, the mass of the dark matter halo that hosts these sources is $\gtrsim
3\times 10^{10}~{}{\rm M}_{\odot}$, while $\zeta\lesssim 1.4$ and
$f_{X}\lesssim 0.1$.

  * •

The probability with which a reionization model can be excluded increases with
increasing $A_{r}$. We found $A_{r}\gtrsim 0.15$ for the 95 per cent
credibility interval of the excluded models at redshift 6.8. This limit
increases for higher redshifts as the upper limits are weaker there. At the
highest redshift of $z=8.7$ we find $A_{r}\gtrsim 27$ for the same credibility
interval. Within 95 per cent credible interval limits, the additional radio
background is at least $0.008\%,0.009\%,0.01\%,1.2\%,1.3\%$ and $1.5\%$ of the
CMB at 1.42 GHz at redshifts 6.5, 6.8, 7.2, 7.8, 8.2 and 8.7 respectively.
These values differ with redshifts because the upper limits have been
considered individually. In the presence of such radio background, the limits
on the IGM parameters for the 95 per cent credible interval of the excluded
models are $\overline{x_{\rm HII}}\lesssim 0.6$, $\overline{T_{\rm K}}\lesssim
10^{3}$ K, $f_{\rm heat}\lesssim 0.6$, $R^{\rm heat}_{\rm peak}\lesssim
10~{}h^{-1}~{}\rm Mpc$, $\Delta R^{\rm heat}_{\rm FWHM}\lesssim
30~{}h^{-1}~{}\rm Mpc$.

The limits on the IGM parameters are similar at all redshifts considered in
this study, indicating that certain IGM states correspond to large amplitudes
of the large-scale power spectra and thus are more easily ruled out by the
upper limits on the 21-cm signal power spectrum. More stringent upper limits
will provide stronger constraints on the IGM parameters. In addition,
combining these with other observational constraints such as the observation
of the global 21-cm signal and high-$z$ galaxies will strengthen the
understanding of the IGM during the EoR.

## Acknowledgements

We would like to thank the anonymous referee for insightful comments. We thank
Cathryn M. Trott for providing the lists of the MWA upper limits and the
corresponding errors, and Bradley Greig for useful discussion related to this
work. We acknowledge that the results in this paper have been achieved using
the PRACE Research Infrastructure resources Curie based at the Tr$\grave{\rm
e}$s Grand Centre de Calcul (TGCC) operated by CEA near Paris, France and
Marenostrum based in the Barcelona Supercomputing Center, Spain. Time on these
resources was awarded by PRACE under PRACE4LOFAR grants 2012061089 and
2014102339 as well as under the Multi-scale Reionization grants 2014102281 and
2015122822. The computations were partly enabled by resources provided by the
Swedish National Infrastructure for Computing (SNIC) at PDC partially funded
by the Swedish Research Council through grant agreement no. 2018-05973. RG and
SZ furthermore acknowledge support by the Israel Science Foundation (grant no.
255/18). GM is supported by Swedish Research Council grants 2016-03581 and
2020-04691.

## DATA AVAILABILITY

The derived data generated in this research will be shared on reasonable
request to the corresponding authors.

## References

  * Bañados et al. (2018) Bañados E., et al., 2018, Nature, 553, 473
  * Barkana (2018) Barkana R., 2018, Nature, 555, 71
  * Barkana & Loeb (2005) Barkana R., Loeb A., 2005, ApJ, 626, 1
  * Barry et al. (2016) Barry N., Hazelton B., Sullivan I., Morales M. F., Pober J. C., 2016, MNRAS, 461, 3135
  * Barry et al. (2019) Barry N., et al., 2019, ApJ, 884, 1
  * Behroozi & Silk (2015) Behroozi P. S., Silk J., 2015, ApJ, 799, 32
  * Bonaldi & Brown (2015) Bonaldi A., Brown M. L., 2015, MNRAS, 447, 1973
  * Bowman & Rogers (2010) Bowman J. D., Rogers A. E. E., 2010, Nature, 468, 796
  * Bowman et al. (2018) Bowman J. D., Rogers A. E. E., Monsalve R. A., Mozdzen T. J., Mahesh N., 2018, Nature, 555, 67
  * Bradley et al. (2019) Bradley R. F., Tauscher K., Rapetti D., Burns J. O., 2019, ApJ, 874, 153
  * Chapman et al. (2016) Chapman E., Zaroubi S., Abdalla F. B., Dulwich F., Jelić V., Mort B., 2016, MNRAS, 458, 2928
  * Cheng et al. (2018) Cheng C., et al., 2018, ApJ, 868, 26
  * Chornock et al. (2013) Chornock R., Berger E., Fox D. B., Lunnan R., Drout M. R., Fong W.-f., Laskar T., Roth K. C., 2013, ApJ, 774, 26
  * Cohen et al. (2020) Cohen A., Fialkov A., Barkana R., Monsalve R. A., 2020, MNRAS, 495, 4845
  * Datta et al. (2007) Datta K. K., Bharadwaj S., Choudhury T. R., 2007, MNRAS, 382, 809
  * Datta et al. (2010) Datta A., Bowman J. D., Carilli C. L., 2010, ApJ, 724, 526
  * Davies et al. (2018) Davies F. B., et al., 2018, ApJ, 864, 142
  * Dawoodbhoy et al. (2018) Dawoodbhoy T., et al., 2018, MNRAS, 480, 1740
  * DeBoer et al. (2017) DeBoer D. R., et al., 2017, Publications of the Astronomical Society of the Pacific, 129, 045001
  * Dixon et al. (2016) Dixon K. L., Iliev I. T., Mellema G., Ahn K., Shapiro P. R., 2016, Monthly Notices of the Royal Astronomical Society, 456, 3011
  * Dowell & Taylor (2018) Dowell J., Taylor G. B., 2018, ApJ, 858, L9
  * Draine & Miralda-Escudé (2018) Draine B. T., Miralda-Escudé J., 2018, ApJ, 858, L10
  * Eastwood et al. (2019) Eastwood M. W., et al., 2019, AJ, 158, 84
  * Ewall-Wice et al. (2018) Ewall-Wice A., Chang T. C., Lazio J., Doré O., Seiffert M., Monsalve R. A., 2018, ApJ, 868, 63
  * Fan et al. (2006) Fan X., et al., 2006, AJ, 132, 117
  * Feng & Holder (2018) Feng C., Holder G., 2018, ApJ, 858, L17
  * Fialkov & Barkana (2019) Fialkov A., Barkana R., 2019, MNRAS, 486, 1763
  * Fialkov et al. (2018) Fialkov A., Barkana R., Cohen A., 2018, Physical Review Letters, 121, 011101
  * Fioc & Rocca-Volmerange (1997) Fioc M., Rocca-Volmerange B., 1997, A&A, 326, 950
  * Fixsen et al. (2011) Fixsen D. J., et al., 2011, ApJ, 734, 5
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
  * Furlanetto et al. (2006) Furlanetto S. R., Oh S. P., Briggs F. H., 2006, Phys. Rep., 433, 181
  * Gehlot et al. (2019) Gehlot B. K., et al., 2019, MNRAS, 488, 4271
  * Gehlot et al. (2020) Gehlot B. K., et al., 2020, MNRAS, 499, 4158
  * Ghara & Choudhury (2020) Ghara R., Choudhury T. R., 2020, MNRAS, 496, 739
  * Ghara & Mellema (2020) Ghara R., Mellema G., 2020, MNRAS, 492, 634
  * Ghara et al. (2015a) Ghara R., Choudhury T. R., Datta K. K., 2015a, MNRAS, 447, 1806
  * Ghara et al. (2015b) Ghara R., Datta K. K., Choudhury T. R., 2015b, MNRAS, 453, 3143
  * Ghara et al. (2016) Ghara R., Choudhury T. R., Datta K. K., 2016, MNRAS, 460, 827
  * Ghara et al. (2017) Ghara R., Choudhury T. R., Datta K. K., Choudhuri S., 2017, MNRAS, 464, 2234
  * Ghara et al. (2018) Ghara R., Mellema G., Giri S. K., Choudhury T. R., Datta K. K., Majumdar S., 2018, MNRAS, 476, 1741
  * Ghara et al. (2020) Ghara R., et al., 2020, MNRAS, 493, 4728
  * (43) Giri S. K., , PhD thesis, Stockholm University, Department of Astronomy, 2019. , p. 68, ISBN 978-91-7797-611-0
  * Giri & Mellema (2020) Giri S. K., Mellema G., 2020, arXiv e-prints, p. arXiv:2012.12908
  * Giri et al. (2018a) Giri S. K., Mellema G., Dixon K. L., Iliev I. T., 2018a, MNRAS, 473, 2949
  * Giri et al. (2018b) Giri S. K., Mellema G., Ghara R., 2018b, MNRAS, 479, 5596
  * Giri et al. (2019a) Giri S. K., Mellema G., Aldheimer T., Dixon K. L., Iliev I. T., 2019a, MNRAS, 489, 1590
  * Giri et al. (2019b) Giri S. K., D’Aloisio A., Mellema G., Komatsu E., Ghara R., Majumdar S., 2019b, Journal of Cosmology and Astro-Particle Physics, 2019, 058
  * Greig & Mesinger (2015) Greig B., Mesinger A., 2015, Monthly Notices of the Royal Astronomical Society, 449, 4246
  * Greig & Mesinger (2017) Greig B., Mesinger A., 2017, MNRAS, 472, 2651
  * Greig et al. (2017) Greig B., Mesinger A., Haiman Z., Simcoe R. A., 2017, MNRAS, 466, 4239
  * Greig et al. (2019) Greig B., Mesinger A., Bañados E., 2019, MNRAS, 484, 5094
  * Greig et al. (2021a) Greig B., Trott C. M., Barry N., Mutch S. J., Pindor B., Webster R. L., Wyithe J. S. B., 2021a, MNRAS, 500, 5322
  * Greig et al. (2021b) Greig B., et al., 2021b, MNRAS, 501, 1
  * Gupta et al. (2017) Gupta Y., et al., 2017, Current Science, 113, 707
  * Harker et al. (2009) Harker G., et al., 2009, MNRAS, 397, 1138
  * Hasegawa & Semelin (2013) Hasegawa K., Semelin B., 2013, MNRAS, 428, 154
  * Hills et al. (2018) Hills R., Kulkarni G., Meerburg P. D., Puchwein E., 2018, Nature, 564, E32
  * Hinshaw et al. (2013) Hinshaw G., et al., 2013, ApJS, 208, 19
  * Hothi et al. (2021) Hothi I., et al., 2021, MNRAS, 500, 2264
  * Islam et al. (2019) Islam N., Ghara R., Paul B., Choudhury T. R., Nath B. B., 2019, MNRAS, 487, 2785
  * Jensen et al. (2013) Jensen H., et al., 2013, MNRAS, 435, 460
  * Kamran et al. (2021) Kamran M., Ghara R., Majumdar S., Mondal R., Mellema G., Bharadwaj S., Pritchard J. R., Iliev I. T., 2021, MNRAS, 502, 3800
  * Kapahtia et al. (2021) Kapahtia A., Chingangbam P., Ghara R., Appleby S., Choudhury T. R., 2021, arXiv e-prints, p. arXiv:2101.03962
  * Kern et al. (2019) Kern N. S., Parsons A. R., Dillon J. S., Lanman A. E., Fagnoni N., de Lera Acedo E., 2019, ApJ, 884, 105
  * Kolopanis et al. (2019) Kolopanis M., et al., 2019, ApJ, 883, 133
  * Krause et al. (2018) Krause F., Thomas R. M., Zaroubi S., Abdalla F. B., 2018, New Astron., 64, 9
  * Liu et al. (2014) Liu A., Parsons A. R., Trott C. M., 2014, Phys. Rev. D, 90, 023019
  * Madau et al. (1997) Madau P., Meiksin A., Rees M. J., 1997, ApJ, 475, 429
  * Majumdar et al. (2012) Majumdar S., Bharadwaj S., Choudhury T. R., 2012, MNRAS, 426, 3178
  * Majumdar et al. (2016) Majumdar S., Datta K. K., Ghara R., Mondal R., Choudhury T. R., Bharadwaj S., Ali S. S., Datta A., 2016, Journal of Astrophysics and Astronomy, 37, 32
  * McGreer et al. (2015) McGreer I. D., Mesinger A., D’Odorico V., 2015, MNRAS, 447, 499
  * Mebane et al. (2020) Mebane R. H., Mirocha J., Furlanetto S. R., 2020, Monthly Notices of the Royal Astronomical Society, 493, 1217
  * Mellema et al. (2015) Mellema G., Koopmans L., Shukla H., Datta K. K., Mesinger A., Majumdar S., 2015, Advancing Astrophysics with the Square Kilometre Array (AASKA14), p. 10
  * Mertens et al. (2018) Mertens F. G., Ghosh A., Koopmans L. V. E., 2018, MNRAS, 478, 3640
  * Mertens et al. (2020) Mertens F. G., et al., 2020, MNRAS, 493, 1662
  * Mesinger & Furlanetto (2007) Mesinger A., Furlanetto S., 2007, ApJ, 669, 663
  * Mevius et al. (2016) Mevius M., et al., 2016, Radio Science, 51, 927
  * Mineo et al. (2012) Mineo S., Gilfanov M., Sunyaev R., 2012, MNRAS, 419, 2095
  * Mitra et al. (2015) Mitra S., Choudhury T. R., Ferrara A., 2015, MNRAS, 454, L76
  * Mondal et al. (2020) Mondal R., et al., 2020, MNRAS, 498, 4178
  * Monsalve et al. (2017) Monsalve R. A., Rogers A. E. E., Bowman J. D., Mozdzen T. J., 2017, ApJ, 835, 49
  * Monsalve et al. (2018) Monsalve R. A., Greig B., Bowman J. D., Mesinger A., Rogers A. E. E., Mozdzen T. J., Kern N. S., Mahesh N., 2018, ApJ, 863, 11
  * Muñoz & Loeb (2018) Muñoz J. B., Loeb A., 2018, arXiv e-prints, p. arXiv:1802.10094
  * Nambissan et al. (2021) Nambissan J., et al., 2021, Experimental Astronomy, 1572-9508
  * Ouchi et al. (2010) Ouchi M., et al., 2010, ApJ, 723, 869
  * Paciga et al. (2013) Paciga G., et al., 2013, MNRAS, 433, 639
  * Park et al. (2019) Park J., Mesinger A., Greig B., Gillet N., 2019, Monthly Notices of the Royal Astronomical Society, 484, 933
  * Parsons et al. (2014) Parsons A. R., et al., 2014, ApJ, 788, 106
  * Patil et al. (2017) Patil A. H., et al., 2017, ApJ, 838, 65
  * Patra et al. (2015) Patra N., Subrahmanyan R., Sethi S., Udaya Shankar N., Raghunathan A., 2015, ApJ, 801, 138
  * Philip et al. (2019) Philip L., et al., 2019, Journal of Astronomical Instrumentation, 8, 1950004
  * Planck Collaboration et al. (2020) Planck Collaboration et al., 2020, A&A, 641, A6
  * Price et al. (2018) Price D. C., et al., 2018, MNRAS, 478, 4193
  * Pritchard & Loeb (2012) Pritchard J. R., Loeb A., 2012, Reports on Progress in Physics, 75, 086901
  * Rasmussen & Williams (2006) Rasmussen C. E., Williams C. K. I., 2006, Gaussian Processes for Machine Learning. The MIT Press, 2006. ISBN 0-262-18253-X
  * Reis et al. (2020a) Reis I., Fialkov A., Barkana R., 2020a, MNRAS, 499, 5993
  * Reis et al. (2020b) Reis I., Fialkov A., Barkana R., 2020b, MNRAS, 499, 5993
  * Ross et al. (2019) Ross H. E., Dixon K. L., Ghara R., Iliev I. T., Mellema G., 2019, MNRAS, p. 1183
  * Ross et al. (2020) Ross H. E., Giri S. K., Dixon K. L., Ghara R., Iliev I. T., Mellema G., 2020, arXiv e-prints, p. arXiv:2011.03558
  * Schneider et al. (2020) Schneider A., Giri S., Mirocha J., 2020, arXiv e-prints, p. arXiv:2011.12308
  * Singh & Subrahmanyan (2019) Singh S., Subrahmanyan R., 2019, ApJ, 880, 26
  * Singh et al. (2017) Singh S., et al., 2017, ApJ, 845, L12
  * Sokolowski et al. (2015) Sokolowski M., et al., 2015, Publ. Astron. Soc. Australia, 32, e004
  * Sun & Furlanetto (2016) Sun G., Furlanetto S. R., 2016, MNRAS, 460, 417
  * Tashiro et al. (2014) Tashiro H., Kadota K., Silk J., 2014, Phys. Rev. D, 90, 083522
  * Thomas & Zaroubi (2008) Thomas R. M., Zaroubi S., 2008, MNRAS, 384, 1080
  * Thomas & Zaroubi (2011) Thomas R. M., Zaroubi S., 2011, MNRAS, 410, 1377
  * Thomas et al. (2009) Thomas R. M., et al., 2009, MNRAS, 393, 32
  * Tingay et al. (2013) Tingay S. J., et al., 2013, Publications of the Astronomical Society of Australia (PASA), 30, 7
  * Trott et al. (2020) Trott C. M., et al., 2020, MNRAS, 493, 4711
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Voytek et al. (2014) Voytek T. C., Natarajan A., Jáuregui García J. M., Peterson J. B., López-Cruz O., 2014, ApJ, 782, L9
  * Wang et al. (2020) Wang F., et al., 2020, ApJ, 896, 23
  * Watson et al. (2013) Watson W. A., Iliev I. T., D’Aloisio A., Knebe A., Shapiro P. R., Yepes G., 2013, Monthly Notices of the Royal Astronomical Society, 433, 1230
  * Wayth et al. (2018) Wayth R. B., et al., 2018, Publ. Astron. Soc. Australia, 35, 33
  * Weiser & Zarantonello (1988) Weiser A., Zarantonello S. E., 1988, Mathematics of Computation, 50, 189
  * Wise et al. (2014) Wise J. H., Demchenko V. G., Halicek M. T., Norman M. L., Turk M. J., Abel T., Smith B. D., 2014, MNRAS, 442, 2560
  * Yang et al. (2020) Yang J., et al., 2020, ApJ, 904, 26
  * Zarka et al. (2012) Zarka P., Girard J. N., Tagger M., Denis L., 2012, in Boissier S., de Laverny P., Nardetto N., Samadi R., Valls-Gabaud D., Wozniak H., eds, SF2A-2012: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics. pp 687–694
  * Zaroubi (2013) Zaroubi S., 2013, in Wiklind T., Mobasher B., Bromm V., eds, Astrophysics and Space Science Library Vol. 396, The First Galaxies. p. 45 (arXiv:1206.0267), doi:10.1007/978-3-642-32362-1_2
  * van Haarlem et al. (2013) van Haarlem M. P., et al., 2013, Astronomy & Astrophysics, 556, A2