Tight Constraints on the Excess Radio Background at $z = 9.1$ from LOFAR

The ARCADE2 and LWA1 experiments have claimed an excess over the Cosmic Microwave Background (CMB) at low radio frequencies. If the cosmological high-redshift contribution to this radio background is between 0.1% and 22% of the CMB at $1.42\,$GHz, it could explain the tentative EDGES Low-Band detection of the anomalously deep absorption in the 21-cm signal of neutral hydrogen (Fialkov&Barkana 2019). We use the upper limit on the 21-cm signal from the Epoch of Reionization ($z=9.1$) based on $141\,$hours of observations with LOFAR (Mertens et al. 2020) to evaluate the contribution of the high redshift Universe to the detected radio background. Marginalizing over astrophysical properties of star-forming halos, we find (at 68% C.L.) that the cosmological radio background can be at most 0.7% of the CMB at $1.42\,$GHz. This limit rules out strong contribution of the high-redshift Universe to the ARCADE2 and LWA1 measurements. Even though LOFAR places limit on the extra radio background, excess of $0.1-0.7$% over the CMB (at $1.42\,$GHz) is still allowed and could explain the EDGES Low-Band detection. We also constrain the thermal and ionization state of the gas at $z = 9.1$, and put limits on the properties of the first star-forming objects. We find that, in agreement with the limits from EDGES High-Band data, LOFAR data disfavour scenarios with inefficient X-ray sources and cases where the Universe was ionized by massive halos only.

do not constrain properties of the first population of starforming objects such as their star formation efficiency, feedback mechanisms that regulated primordial star formation, and the properties of the first sources of heat (e.g., Xray binaries). These properties can be probed using lowfrequency radio observations of the redshifted 21-cm signal of neutral hydrogen (e.g., Pober et al. 2014;Greig et al. 2016;Singh et al. 2017;Monsalve et al. 2018Monsalve et al. , 2019. The 21-cm signal is produced by atomic hydrogen in the IGM. The hyper-fine splitting of the lowest energy level of a hydrogen atom gives rise to the rest-frame ν21 = 1.42 GHz radio signal with the equivalent wavelength of about 21 cm (see Barkana 2016, for a recent review). Owing to its dependence on the underlying astrophysics and cosmology, this signal is a powerful tool to characterise the formation and the evolution of the first populations of astrophysical sources and, potentially, properties of dark matter, across cosmic time. Because the 21-cm signal is measured against the diffused radio background, usually assumed to be only the Cosmic Microwave Background (CMB), this signal can also be used to constrain properties of any excess background radiation at low radio frequencies.
Recently, a detection of the global 21-cm signal from z ∼ 17 was reported by the EDGES collaboration . The reported signal significantly deviates from standard astrophysical models (e.g., Cohen et al. 2017, show a large set of viable 21-cm global signals varying astrophysical parameters in the broadest possible range) and concerns about the signal being of cosmological origin have therefore been expressed (Hills et al. 2018;Sims & Pober 2019;Singh & Subrahmanyan 2019;Spinelli et al. 2019;Bradley et al. 2019). Despite these concerns, several theories have been proposed to explain the stronger than expected absorption. Over-cooling of hydrogen gas by dark matter (Barkana 2018) has been proposed as a possible solution. Alternatively, the existence of a new component of radio background at low radio frequencies in addition to the CMB could also lead to a deeper 21-cm absorption feature due to the stronger contrast between the temperatures of the background and the gas (e.g., Bowman et al. 2018;Feng & Holder 2018). Astrophysical sources such as supernovae or accreting supermassive black holes (Biermann et al. 2014;Ewall-Wice et al. 2018Jana & Nath 2018;Mirocha & Furlanetto 2019) could produce such an extra radio background. However, these sources would need to be several orders of magnitude more efficient in producing synchrotron radiation than corresponding sources at low redshifts, which is not very likely. An extra radio background can also be created by more excessbackground agents such as active neutrinos (Chianese et al. 2018), dark matter (Fraser et al. 2018;Pospelov et al. 2018) or superconducting cosmic strings (Brandenberger et al. 2019). Interestingly, excess radio background at low radio frequencies was claimed by the ARCADE2 collaboration at 3 − 90 GHz (Fixsen et al. 2011) as well as by LWA1 at 40 − 80 MHz (Dowell & Taylor 2018). Specifically, the latter measurement shows that the excess can be fitted by a power law with a spectral index of −2.58±0.05 and a brightness temperature of 603 +102 −92 mK at the reference frequency 1.42 GHz. However, the nature of this excess is still debated (Subrahmanyan & Cowsik 2013).
Apart from the EDGES Low-Band, several other global signal experiments report upper limits: the Large-Aperture Experiment to Detect the Dark Ages (LEDA, Price et al. 2018) reported an upper limit of 890 mK on the amplitude of the 21-cm signal at z ∼ 20 (Bernardi et al. 2016); at lower redshifts both the EDGES High-Band collaboration (z ∼ 6.5 − 14.8, Monsalve et al. 2017) and the Shaped Antenna measurement of the background RAdio Spectrum (SARAS, z ∼ 6.1 − 11.9, Singh et al. 2017) reported non-detections. The data collected by the latter two experiments rule out scenarios with negligible X-ray heating, placing limits on the properties of first star forming objects and X-ray sources (Monsalve et al. , 2019Singh et al. 2017Singh et al. , 2018. In parallel, interferometric radio arrays are placing upper limits on the fluctuations of the 21-cm signal, including the Low-Frequency Array (LOFAR 1 , Patil et al. 2017;Gehlot et al. 2019;Mertens et al. 2020 The recently reported LOFAR measurements (Mertens et al. 2020) are based on 141 hours of observations and are currently the tightest upper limits on the 21-cm power spectrum from z = 9.1, allowing us to rule out scenarios with cold IGM for models with the CMB as the background radiation (Ghara et al. 2020). In this paper we use these upper limits to constrain any excess radio background. We also derive limits on astrophysical parameters and the properties of the IGM with and without the excess radio contribution. This paper is structured as follows: In Section 2, we describe the simulations used to generate the mock data sets of the 21-cm signals. In Section 3 we describe the mock data set and the ranges of parameters probed. In Section 4, we discuss the statistical analysis employed to constrain the model parameters. In Section 5, we report our constraints on the amplitude of the excess radio background and compare it to the values that could explain the EDGES Low Band detection. We also place limits on the thermal and ionization state of the gas at z = 9.1, and on the properties of the first star-forming objects. We provide a qualitative comparison with the results of Ghara et al. (2020) in Section 6. Finally, we conclude in Section 7.

SIMULATED 21-CM SIGNAL
The 21-cm signal of neutral hydrogen observed against a background radiation of the brightness temperature T rad (at 1.42 GHz at redshift z), depends on the processes of cosmic heating and ionization. The brightness temperature of the 21-cm signal is given by: where TS is the spin temperature of the transition which at Cosmic-Dawn redshifts is coupled to the kinetic temperature of the gas, Tgas, through Ly-α photons produced by stellar sources (Wouthuysen 1952;Field 1958). The value of τ21 is the optical depth at redshift z given by where dv/dr = H(z)/(1 + z) is the gradient of the line of sight component of the comoving velocity field, H(z) is the Hubble parameter at z, and nH is the neutral hydrogen number density at z which depends on the ionization fraction and is driven by both ultraviolet and X-ray photons. The spin temperature encodes complex astrophysical dependencies and can be written as where xC is the collisional coupling coefficient and xα is the Wouthousen-Field coupling coefficient (Wouthuysen 1952;Field 1958). Both xC and xα depend on the value of T rad : with Pα being the total rate (per atom) at which Ly-α photons are scattered within the gas and T * is the effective temperature of the 21-cm transition (0.068 K). The collisional coupling coefficient is where κ i 10 is the rate coefficient for spin de-excitation in collisions with the species of type i of density ni, where we sum over species i (see e.g., Barkana 2016, for a recent review).

Radio Background
Usually the CMB is assumed to be sole contributor to the background radiation and T rad = TCMB(1 + z), where TCMB is the present day value of the CMB temperature, 2.725 K. However, as was mentioned in Section 1, the anomalously strong EDGES Low-Band signal has encouraged the development of alternative models in which the radio background is enhanced (e.g., Bowman et al. 2018;Feng & Holder 2018). Here we adopt a phenomenological global extra radio background with a synchrotron spectrum in agreement with observations by LWA1. The total radio background at redshift z is then given by where ν obs is the observed frequency, Ar is the amplitude defined relative to the CMB temperature and calculated at the reference frequency of 78 MHz (which is the centre of the absorption trough reported by the EDGES collaboration) and β = −2.6 is the spectral index (in agreement with the LWA1 observation). We vary the value of Ar between 0 and 400 at 78 MHz with the upper limit being close to the LWA1 limit and corresponds to 21% of the CMB at 1.42 GHz. All values of Ar between 1.9 (equivalent to 0.1% of the CMB at 1.42 GHz) and 400 were shown to explain the EDGES Low detection (for a tuned set of astrophysical parameters; see more details of the modelling in Fialkov & Barkana 2019).

Astrophysical Parameters
Astrophysical processes affect the 21-cm signal by regulating the thermal and ionization states of the gas. In our modelling, we account for the effect of radiation (Ly-α, Lyman-Werner, X-ray and ionizing radiation) produced by stars and stellar remnants on the 21-cm signal. The process of star formation is parameterized by two parameters. The first parameter is the star formation efficiency, f * , defined as the amount of gas in halos that is converted into stars, which we vary in the range f * = 0.1% to 50%. The second parameter is the value of circular velocity of dark matter halos, Vc, which is varied between 4.2 km s −1 (molecular hydrogen cooling limit, corresponding to the dark matter mass of M h = 1.5 × 10 6 M ⊙ at z = 9.1) and 100 km s −1 (M h = 2×10 10 M ⊙ at z = 9.1). The high values of Vc implicitly take into account various chemical and mechanical feedback effects (e.g., the supernovae feedback which is expected to expel gas from small halos thus rising the threshold mass for star formation) which we do not include explicitly. Cooling of gas via molecular hydrogen cooling channel, and subsequent star formation, happens in small halos of circular velocity 4.2 km s −1 < Vc < 16.5 km s −1 (M h ∼ 10 5 − 10 7 M ⊙ ). Abundance of molecular hydrogen is suppressed by Lyman-Werner (LW) radiation (Haiman et al. 1997;Fialkov et al. 2013). Additional inhomogeneous suppression is introduced by the relative velocity between dark matter and baryons, v bc (Tseliakhovich & Hirata 2010), which imprints the pattern of Baryon Acoustic Oscillations (BAO) in the 21-cm signal (Dalal et al. 2010;Maio et al. 2011;Visbal et al. 2012). Higher mass halos (Vc > 16.5 km s −1 ) form stars owing to atomic hydrogen cooling and are sensitive to neither the LW feedback nor to the effect of v bc , but are affected by photoheating feedback (Sobacchi & Mesinger 2013;Cohen et al. 2016;Sullivan et al. 2018). X-ray sources re-heat and mildly re-ionise the gas after the period of adiabatic cooling. We assume hard X-ray spectral energy distribution (SED) typical for a population of high redshift X-ray binaries (a complex function of X-ray energy with a peak at ∼ 2 keV adopted from Fragos et al. 2013). We vary the normalisation constant fX defined as the ratio of the total X-ray luminosity to the star formation rate, A value fX = 1 yields LX normalized to observations of X-ray binaries found in low metallicity regions today (see Fragos et al. 2013). Here we explore the wide range fX = 10 −6 − 100. Values of fX 100 are unlikely as such a population would saturate the unresolved X-ray background observed by the Chandra X-ray Observatory , while values fX 10 −6 contribute negligible X-ray heating.
In our simulations the effects of ionizing radiation (ultraviolet radiation from stars) are defined by two parameters: the mean free path of ionizing photons, R mfp = 10 − 70 comoving Mpc, and the ionizing efficiency of sources, ζ, which is tuned to yield the CMB optical depth τ in the range between 0.045 and 0.1. For a fixed set of astrophysical parameters, either ζ or τ can be used as parameter (for more details on the relation between ζ and τ see Cohen et al. 2019a). Here we choose to use the latter as it is directly probed by the CMB experiments. The latest values of τ measured by the Planck satellite are τ = 0.054 ± 0.007 (e.g., Planck Collaboration et al. 2018). However, because in this paper we focus on the constraints imposed by the LOFAR upper limits, we explore a broader range of values (0.045 − 0.1) including higher values of τ which can be constrained by the LOFAR data. We

MOCK DATA SETS AND PARAMETER SETS
We use a hybrid computational framework (similar to, but independent of the publicly available 21cmFAST code of Mesinger et al. 2011) to estimate the evolution of the large scale 21-cm signal. The code takes into account the effects of inhomogeneous star formation, thermal and ionization histories. Processes on scales below the resolution scale of 3 comoving Mpc (such as star formation, LW and photoheating feedback effects, effects of v bc ) are implemented using subgrid prescriptions. Radiation produced by stars and stellar remnants is propagated accounting for the effects of redshift on the energy of the photons and absorption in the IGM. Reionization is implemented using an excursion set formalism (Furlanetto et al. 2004). Astrophysical parameters (f * , Vc, fX, τ, R mfp , Ar, and the spectral energy distribution of X-ray photons taken from Fragos et al. 2013) are received as an input. The code generates cubes of the 21-cm signal at every redshift along with the temperature of the neutral IGM, ionization state, intensity of the Ly-α and LW backgrounds (see more details in Visbal et al. 2012;Cohen et al. 2017;. The comoving volume of each simulation box is 384 3 Mpc 3 . The simulation is run from z = 60 to z = 6. Using the same set of initial conditions for the distribution and velocities of dark matter and baryons (the fiducial IC), and then varying the astrophysical and background parameters, we construct a total of 20762 models, 6515 of which have a boosted radio background with respect to the CMB and are referred to as the excess-background models, while the remainder are reference standard models with Ar = 0 (used as a separate data set).
For each simulation, we calculate the values of the spherically averaged binned 21-cm power spectra P (kc), where kc is the centre of a wave-number bin chosen by Mertens et al. (2020). The power spectrum is averaged over redshifts z = 8.7 − 9.6 (to account for the LOFAR bandwidth), and binned over wave-numbers in agreement with the LOFAR observational setup (see Table 1 for the details of the wave-number binning). From each simulation, we also extract: the mean temperature of the gas in neutral regions Table 1. Summary of LOFAR measurements directly taken from Table 4 of Mertens et al. (2020). From left to right: central mode of each bin in units of h Mpc −1 , the extent of each k bin, spherically averaged power spectrum in each bin, 1 − σ error in the binned power spectrum. at z = 9.1, Tgas, the mean ionization fraction at z = 9.1, xH ii, the redshift at which the ionization fraction (of volume) is 50%, zre, and the duration of reionization ∆z, defined as the redshift range between the epoch when the mean ionization fraction was 90% and 10%.

Sample Variance
The lowest wave-number observed by Mertens et al. (2020) with LOFAR is kc = 0.075 h Mpc −1 , which corresponds to the scale of ∼ 125 comoving Mpc and is a significant fraction of the size of our simulation box (384 Mpc). Therefore, power spectrum in the lowest k-bin is subject to strong statistical fluctuations due to sample variance, as is shown in Figure 1. For the set of initial conditions that we used to generate the entire data-set (our fiducial IC) the bin-averaged power spectrum in the lowest k-bin is 1.1σ away from the mean calculated over 18 realizations. We correct for this systematic offset by introducing a bias factor. We perform a new suite of simulations to systematically estimate the effect of sample variance. For each set of astrophysical parameters out of 360 selected combinations 6 , 10 simulations with different initial conditions, including the fiducial set, were performed. The bias in the binned power spectrum was subsequently calculated at every k-bins as the ratio of the binned power spectrum averaged over 10 realizations to the one derived from the fiducial set: We find that at z = 9.1 (close to the mid-point of reionization for the models that can be constrained by LOFAR in the standard case) the bias varies as a function of the reionization parameters τ and R mfp , while it has a very weak dependence on the rest of the parameters (Vc, f * , fX and Ar). We jointly fit the bias as a second order polynomial in τ times a linear function of R mfp . Because the entire data set described in Section 3 was created using the fiducial IC set, we apply the corresponding parameter-dependent and kc-dependent bias factor to all the simulated results to compensate for the effect of sample variance. Multiplying by the bias factor is essentially equivalent to averaging over 10 simulations. Furthermore, we fit the variation in the simulated power in each bin (σSV,sim(kc), blue error bars in Figure 1), as a function of astrophysical parameters. We find that the fractional standard deviation, σSV,sim(kc)/P fiducial (kc), can be fitted with a quadratic function of τ times a linear function of R mfp , similarly to bSV(kc). The variation due to sample variance has a very weak dependence on Vc, f * , fX and Ar. The error in the power spectrum (after it has been corrected by the bias factor) is then given by σSV,sim(kc)/ √ 10. Finally, in order to account for theoretical uncertainty in modelling 7 , we impose a lower limit of 10% on the relative error of the power spectrum of each individual simulation (Ghara et al. 2020). In the total error budget of the corrected power spectrum this error should also be divided by √ 10.
The total theoretical parameter-dependent error in the 7 The values of the 21-cm signal generated by the numerical simulation are subject to uncertainty. This is because some of the effects of order ∼ (1 + δ), where δ is the stochastic dimensionless perturbation of the density field, have not been taken into account. For example, at the moment we assume linear growth of structure on large scales (> 3 Mpc).
binned power spectrum is, thus, given by where ∆ 2 th (kc) = P fiducial (kc)k 3 c / 2π 2 is the calculated power spectrum in mK 2 units.

Binning over the Model Parameters
Our goal is to derive constraints on the excess radio background and also explore implications for the rest of the model parameters, as well as for the thermal and ionization states of the IGM. Based on the value of the power spectrum for each set of model parameters, we evaluate the likelihood of each point in the parameter space θ as described in the next Section. We, therefore, need to bin the parameter space θ and calculate the binned power spectra ∆ 2 th (kc, θ) and the corresponding theoretical error σ th (kc, θ). To remind the reader, ∆ 2 th (kc, θ) and σ th (kc, θ) are binned in redshift, wave-number and θ.
We explore two discrete sets of the parameter spaces with θ defined as either the model parameters θ = [f * , Vc, fX , τ, R mfp , Ar] or the derived quantities θ = [Tgas,xH ii, zre, ∆z]. The range of each parameter is divided into 10 equally spaced bins, and each bin is tagged by the bin-averaged value of relevant parameters. Due to the large ranges, the binning is logarithmic for f * , Vc, fX , Ar and Tgas, and linear for τ , R mfp ,xH ii, zre and ∆z. We assume flat priors on each of the parameters across the entire allowed range (see Section 2): 0.001 f * 0.5, 4.2 km s −1 Vc 100 km s −1 , 10 −6 fX 100, 0.045 τ 0.1, 10 R mfp 70 comoving Mpc and zero outside these ranges. In the standard case Ar = 0 and in the excess background case, we vary 0.2 Ar 400 (thus covering the range 0.01%-21% of the CMB at 1.42 GHz). The priors on [Tgas,xH ii, zre, ∆z] are defined based on the ranges of these parameters found in our simulations (Section 3): 2.2 K Tgas 400 K (the lower limit is the temperature of the gas in an adiabatically expanding universe at z = 9), 0.02 xH ii 1.00, 6 zre 10 and 2 ∆z 5, and zero outside these ranges.
For θ = [f * , Vc, fX , τ, R mfp , Ar] this binning gives rise to 10 5 bins in the standard case and 10 6 bins in the excessbackground case; for θ = [Tgas,xH ii, zre, ∆z] there are 10 4 bins in each case. Due to the relatively small number of models, not all bins are populated. To solve this issue, we use the model sets to train Artificial Neural Networks (ANN, see Appendix A for details) and use that to construct an emulator (similar approach has been taken by Cohen et al. 2019b;Kern et al. 2017), which we then use to interpolate the empty bins.

STATISTICAL ANALYSIS METHODOLOGY
In general, the 21-cm signal is expected to be a non-Gaussian field (Bharadwaj & Pandey 2005;Mellema et al. 2006;Mondal et al. 2015) and the non-Gaussian effects will play a significant role in the error estimates of 21-cm power spectrum (Mondal et al. 2016(Mondal et al. , 2017. In addition, the data in LOFAR bins are slightly correlated due to the finite station size. Therefore, the power spectrum error-covariance matrix is expected to be non-diagonal. However, in reality bins show very weak correlation because the bins are chosen relatively wide compared to the footprint of a LOFAR station which acts as a spatial convolution kernel. With minimal error, we can therefore assume that the bins are uncorrelated and the covariance matrix is diagonal. The probability of a model (tagged by θ) given data can then be written as a product of the probabilities in each individual wave-number bin kc ∈ ki. In addition, because of the bin-averaging and large scales considered, we can assume that the signal is close to a Gaussian random field.
The LOFAR measurements reported by Mertens et al. (2020) are upper limits. Therefore, following Ghara et al. (2020), we can represent the probability of a model given the observed power spectrum values using the error function: where ∆ 2 21 (ki) is the measured power spectrum in the i-th kc bin with uncertainty ∆ 2 21, err (ki) listed in Table 1. The total variance in the bin is given by According to this definition, the probability of a model is close to unity when its power spectrum at z = 9.1 is less than [∆ 2 21 (ki)−σ(ki, θ)] for all ki, and the probability is close to zero when ∆ 2 th (ki, θ) is greater than [∆ 2 21 (ki) + σ(ki, θ)] for any ki.
As an illustration, in Figure 2 we show the complete set of excess-background power spectra (6515 models in total) colour-coded by the probability that the data is consistent with the model. For comparison, we also show the maximum power of the models in the standard case (white line). The upper limits from Mertens et al. (2020) are plotted for reference. As we see from the figure, the current observational limits from LOFAR are strong enough to rule out a significant fraction of the explored excess-background scenarios (all corresponding to a cold IGM with ∼ 50% ionization at z = 9.1, as we will see later). However, for the standard astrophysical scenarios where the values of the power spectra are lower, only the most extreme models can be ruled out, and only in the lowest k-bin. A set of corresponding thermal histories is plotted in the right panel of Figure 2. The LOFAR upper limits by Mertens et al. (2020) disfavour a late X-ray heating which leaves the IGM cold for most of the EoR. Scenarios with early X-ray heating cannot be ruled out by the data as, typically, the corresponding power spectrum values are low.

RESULTS
Using the predicted values of the spherically averaged binned power spectrum in all seven k-bins we can rule out scenarios which yield strong fluctuations at z = 9.1. A few factors need to come together to ensure maximum power. First, the spin temperature has to be fully coupled to the gas temperature, which, for realistic star formation scenarios, is guaranteed to be the case at z = 9.1 (e.g., Cohen et al. 2018). Second, the larger the contrast between Tgas and T rad , the stronger the signal. Because we assume the evolution of T rad to be independent of star formation, the strongest contrast between Tgas and T rad is reached in cases of cold IGM (for both excess-background and standard models). In addition, owing to the higher radio background temperature, the signals are enhanced in our excess-background case compared to the standard case. Finally, fluctuations in the gas temperature and the neutral fraction play a role. Because here we have chosen a hard X-ray spectrum (Fragos et al. 2013;, heating is nearly homogeneous, and the dominant source of fluctuations at z = 9.1 is the nonuniform process of reionization with peak power at ∼ 50% ionization fraction. For a fixed thermal history, nearly homogeneous reionization would result in a smoother signal and, thus, lower power of the 21-cm fluctuations, compared to a patchy reionization scenario.

Limits on the Excess-Radio Background
Using L( θ), we calculate the normalized probability for each of the parameters, θ = [f * , Vc, fX , τ, R mfp , Ar], and parameter pairs, marginalising over the rest of the parameter space. The resulting probability distributions are normalized using the criterion that the total probability (area under the curve) is 1 within the considered prior ranges. The resulting two-dimensional and one-dimensional probabilities of all the model parameters are shown in Figure 3, where we divided each probability function by its peak value to show the marginalised likelihood of all possible combinations uniformly. Using one-dimensional probabilities we find the 68% confidence intervals for each one of the parameters (see Table 2). We calculate the 68% confidence levels (C.L.) by selecting parameter-bins with the highest probability up to a cumulative probability of 0.68. We also find the regions where the one-dimensional probabilities are below exp(−1/2) of the peak.
Marginalising over the residual model parameters (f * , Vc, fX , τ, R mfp ), we derive constraints on Ar finding that LOFAR upper limit rules out Ar > 13.2 at 68%, equivalent to 0.7% of the CMB at 1.42 GHz. The likelihood, which peaks at low values of Ar, drops by a factor of exp(−1/2) by Ar = 60.9, which corresponds to 3.2% of the CMB at 1.42 GHz. In our analysis we have fixed the value of the spectral index of the radio background to β = −2.6. We have checked that the uncertainty in the spectral index ∆β = 0.05, reported by LWA1 (Dowell & Taylor 2018), would lead to only up to ∼ 3 per cent variation in the intensity of the excess radio background at the frequency of 140 MHz corresponding to z = 9.1.  showed that the global signal reported by EDGES Low-Band can be produced by adding an extra radio background with 1.9 < Ar < 418 relative to the CMB at the 78 MHz reference frequency (corresponding to 0.1 − 22 per cent of the CMB at 1.42 GHz). Even though part of this range is now ruled out by the new LOFAR limits, models with values of Ar between 0.1% and 0.7% of the CMB at 1.42 GHz are still allowed and could fit the EDGES Low-Band detection. Such a small contribution is within the measurement error of LWA1 (Dowell & Taylor 2018, report excess background of 603 +102 −92 mK at the 21-cm rest frame frequency of 1.42 GHz) and would remain a plausible ex- Figure 2. We show the excess-background models colour-coded with respect to the probability that the data is consistent with the model (Eq. 10) as is indicated on the colour bar. Left: Binned power spectra vs wavenumber (in units of Mpc −1 , where we have assumed h = 0.6704 for conversion from Table 1) at z = 9.1. The white dashed line shows the maximum power of the models in the standard case (L = 0.4898). Magenta data points correspond to the LOFAR data from Table 1 (two-sided error bars). Right: corresponding thermal histories, i.e., evolution of the mean temperature of neutral intergalactic gas with redshift. Each curve is shown down to the (model-dependent) redshift of end of reionization. Table 2. Limits on astrophysical parameters and the derived IGM parameters. From left to right: type of model; mean temperature of neutral gas at z = 9.1 in K; ionization fraction of the IGM; duration of reionization defined as the redshift interval between 90% neutral IGM and 10% neutral; redshift of the mid-point of reionization, zre (defined as the redshift at which neutral fraction is 50%); star formation efficiency; minimum circular velocity of star forming halos in km s −1 ; X-ray heating efficiency; CMB optical depth; mean free path of ionizing photons in comoving Mpc; amplitude of the excess radio background compared to the CMB at the reference frequency of 78 MHz (as defined in Eq. 6). For the case of excess radio background (Ex. bck. in the table) we show both 68% limits (top row) and we find the parameter values at which the likelihood drops to exp(−1/2) of the peak value (middle row). In the standard case we can only show the 68% limits, as the 1D PDFs are rather flat.

Model
Tgasx planation for the detected EDGES signal even if the excess measured by ARCADE2 and LWA1 is due to an erroneous Galactic modelling (Subrahmanyan & Cowsik 2013). In Figure 4, as an illustration, we show global 21-cm signals for those excess-background models from our data set that are broadly consistent with the tentative EDGES Low-Band detection. In order to define this consistency we follow the simple approach taken by   −min [T21(68 < ν < 88)] < 1 K. (13) The signals in Figure 4 are colour-coded with respect to the LOFAR likelihood (same as in Figure 2). All the signals consistent with EDGES Low-Band have relatively high LOFAR likelihood, L 0.31. This is because the EDGES detection implies an early onset of the Ly-α cou-pling (Schauer et al. 2019) due to efficient star formation (f * > 2.8%) in lower-mass halos with circular velocity below Vc = 45 km s −1 (corresponding to M h < 7.8×10 8 M ⊙ at z = 17 . In such models the IGM is heated and partially ionized by z = 9.1, resulting in relatively low-intensity 21-cm signals in the LOFAR band.

Astrophysical Limits
Next, we explore the implications of the LOFAR upper limits for the rest of the model parameters (f * , Vc, fX , τ, R mfp ). In this work we assume hierarchical structure formation with a simple prescription for the formation of stars and X-ray binaries. Therefore, LOFAR limits at z = 9.1 can be used to constrain properties of the first star forming halos (appearing at z ∼ 30 − 60 in our simulations) and first sources of light at Cosmic Dawn. The resulting two-dimensional and one-dimensional probabilities are shown in Figure 3 and the limits are summarised in Table 2. In the limiting case of the negligible radio background our results converge to the standard case with the CMB as the background radiation. This trend is demonstrated in Figure 3 where the two-dimensional probabilities of models with Ar = 0 are appended below the  Figure 3. 1D and 2D marginalised likelihood of the astrophysical parameters (f * , Vc, f X , τ, R mfp , Ar) obtained using excessbackground models. In addition, we append the normalized likelihood values for our standard models below the white band of the bottom row to highlight the consistency with the excess-background case. The standard-case normalized likelihood was calculated by using a joined set of the excess-background models and standard models. Regions of 2D marginalised likelihoods which are on the darker (red, purple and black) side of the solid lines are disfavoured with more than 68% confidence, and the regions which are on the darker side of the dashed lines are below exp(−1) of the peak probability (similar to the 2D Gaussian 1-sigma definition). The grey regions in the 1D likelihood distribution are also disfavoured at the 68% confidence level, and the black regions are below exp(−1/2) of the peak probability. Threshold parameter values at which likelihood drops by a factor of exp(−1/2) and the 68% limits are listed in Table 2. white band. For completeness we also explore the set of standard models separately, showing their two-dimensional and one-dimensional probabilities in Figure 5 and listing the corresponding 68% limits in Table 2. All disfavoured models feature efficient star formation with f * 5% at 68% C.L. (Table 2). However, the corresponding 1D marginalised likelihood is rather flat and never drops below a factor of exp(−1/2) relatively to its peak value. Higher values of f * result in stronger fluctuations which are easier to rule out. Higher values of f * also imply stronger Ly-α background and, thus, an earlier onset of Ly-α coupling which yields signals with larger amplitudes (e.g., Cohen et al. 2019a).
Another model parameter related to star formation in first halos is Vc. Higher Vc is equivalent to larger minimum mass of star forming halos which are more strongly clustered, thus yielding stronger fluctuations. In the hierarchical model of star formation that we adopt here, higher Vc also implies later onset of star formation and X-ray heating. In such models chances are that fluctuations (e.g., heating and Ly-α) are not yet saturated by z = 9.1 resulting in stronger 21-cm signals that can be ruled out by LOFAR. We find that values of Vc above 28 km s −1 (corresponding to 4.5×10 8 M ⊙ at z = 9.1) are disfavoured by the data at 68% (the corresponding 1D marginalised likelihood is rather flat and never drops below the threshold value of exp(−1/2) relatively to its peak value). The standard-physics limit is 36 km s −1 , or 9.5 × 10 8 M ⊙ at z = 9.1. Even though the limits on Vc are weak at the moment, the LOFAR data favour the existence of low-mass halos (in agreement with EDGES High Band results, Monsalve et al. 2019). In our models gas temperature is regulated by the interplay between several cooling and heating mechanisms with the major roles played by adiabatic cooling due to the expansion of the universe and X-ray heating by X-ray binaries, although the latter is partially degenerate with f * and Vc which regulate the number of X-ray binaries 8 . Therefore, the X-ray efficiency of the first X-ray binaries is directly constrained by LOFAR with a value fX < 1 × 10 −2 disfavoured at 68% C.L., implying a lower limit on the total X-ray luminosity per star formation rate (Eq. 7) of 3 × 10 38 erg s −1 M −1 ⊙ yr. The 1D likelihood, which peaks at high fX values, is steep enough and drops below the threshold exp(−1/2) of its peak value at fX = 8 × 10 −4 (corresponding to 2.4 × 10 37 erg s −1 M −1 ⊙ yr). In the standard case only the 68% limit can be calculated and is fX < 5 × 10 −3 (1.5 × 10 38 erg s −1 M −1 ⊙ yr respectively). The current LOFAR data also disfavour models with mid-point of reionization at z ∼ 9. In such models the peakpower from ionizing fluctuations falls within the currentlyanalysed LOFAR band, and, consequently, such models are relatively easy to exclude. This constraint can be mapped on to limits on τ : scenarios with τ > 0.076 (excess background) or τ > 0.080 (standard models) are disfavoured at 68%. In both theories, the 1D likelihood curves of τ peak at low values of τ but do not drop below the threshold value of exp(−1/2) within the prior ranges. Finally, we find that the constraints on the model parameter R mfp are very weak, with the 1D marginalised likelihood being very flat. This means that our model power spectrum is not sensitive to the changes in R mfp value at z ∼ 9.

Comparison with EDGES
Focusing on the standard models we can compare the LO-FAR limits reported above to the limits extracted from the data of the global 21-cm instrument EDGES High-Band (90 − 190 MHz). Using a similar set of standard models and similar prior ranges of parameters as we explore here, Monsalve et al. (2019) found that the EDGES High-Band data favour (at 68% confidence) the following parameter ranges assuming a fixed X-ray SED (softer than what we use here; however, the global signal constraints prove to be nearly insensitive to the X-ray SED, Monsalve et al. 2019): R mfp < 36.1 Mpc, Vc < 21.5 km s −1 (equivalent to 2 × 10 8 M ⊙ at z = 9.1), fX > 2.5 × 10 −3 , f * < 0.4% or f * > 3.6% (signals with both lower and higher values of f * are likely to be outside of the band of EDGES High), τ < 0.072 or 0.074 < τ < 0.079 (where the second band is most likely due to the instrumental systematic). Overall LOFAR and the EDGES High-Band are in agreement ruling out scenarios with inefficient X-ray heating and models in which the Universe was ionized by massive halos only (of mass few ×10 8 M ⊙ or higher, at z ∼ 9.1). Similar trends were found with the SARAS2 data (although only ∼ 300 models were examined in that case, Singh et al. 2017).

Limits on the Thermal and Reionization Histories
We use the LOFAR upper limits on the 21-cm power spectrum to put limits on the thermal and ionization state of the IGM at z = 9.1. We repeat the likelihood calculation applying it to the IGM parameters θ = [Tgas,xH ii, zre, ∆z]. The resulting two-dimensional and one-dimensional probabilities of Tgas,xH ii, zre and ∆z are shown in Figure 6 (left panel shows the case of the extra radio background, while the standard case is shown on the right for comparison). Our results are also summarised in Table 2 where the 68% limits are listed. We have also tabulated the limits obtained from the regions where the one-dimensional probabilities are below exp(−1/2) of the peak (similar to the Gaussian 1-sigma definition).
As we see from the figure and the table, the LOFAR data indeed disfavour scenarios with cold IGM. The lower limit on the temperature of neutral gas at z = 9.1 is 19.2 K at 68% C. L. (while it is only 10.1 K in the standard case). The likelihood, which peaks at high values of Tgas, drops by a factor of exp(−1/2) at Tgas = 6 K in the excess-background Ar Figure 5. 1D and 2D marginalised likelihood of the astrophysical parameters (f * , Vc, f X , τ, R mfp ) obtained using standard models (Ar = 0). The regions of 2D marginalised likelihoods which are on the darker side of the solid lines are disfavoured with more than 68% confidence. The grey regions in the 1D likelihood distribution are also disfavoured at the 68% confidence level (also listed in Table 2).
case. As expected, there is some degree of degeneracy between the constraints on the thermal and reionization histories with the strongest limits on temperature coming from the cases with mid-point of reionization occurring at z ∼ 9.
Through marginalising over the thermal histories we can put limits on the process of reionization ( Figure 6 and Table 2). We find that the LOFAR limits disfavour fast reionization scenarios (with ∆z 3) with ionized fractions between ∼ 38% and ∼ 72% at z = 9.1. The high end of the allowedxH ii values (xH ii > 72% at z = 9.1) is inconsistent with other probes of reionization and would be ruled out if joined constraints were considered: e.g., Ly-α damping wing absorption in the spectrum of the quasar at z = 7.54 suggests that the Universe is ∼ 60% neutral at that redshift (ionization fraction less than 40% Bañados et al. 2018;Davies et al. 2018). The quantitative joint analysis, however, is beyond the scope of this paper. The grey regions in the 1D likelihood distribution are also disfavoured at the 68% confidence level, and (for excess-background models) the black regions are below exp(−1/2) of the peak probability. Threshold parameter values at which likelihood drops by a factor of exp(−1/2) and the 68% limits are listed in Table 2.

Ghara
cal constraints on the IGM properties assuming standardphysics models (with the CMB as the background radiation). We verify the consistency of our conclusions with Ghara et al. (2020) by qualitatively comparing our standard case results for the thermal and ionization states of the IGM. A quantitative comparison between the two works is not possible at this stage because of the different choices of modelling, parameterization and priors. Moreover, because the 21-cm signal is sensitive to the thermal and ionization histories the values of the gas temperature and ionization fraction can be directly constrained using the data. However, the mapping between these quantities and the astrophysical properties of the UV and X-ray sources (in our case f * , Vc, fX , τ and R mfp ) is model-dependent. Therefore, in this paper we refrain from comparing the astrophysical constraints leaving it for future work.
In their work Ghara et al. (2020) explored two scenarios: (1) homogeneous spin temperature, which implies saturated Ly-α background and homogeneous X-ray heating. The parameters that are varied in this case include gas temperature (or, equivalently, spin temperature), minimum halo mass and ionizing efficiency. (2) Inhomogeneous heating by soft X-ray sources with power-law SED where Mmin and the spectral index of X-ray sources were kept fixed; the parameters that were varied are ionizing efficiency, X-ray efficiency (defined differently than in our work) and minimum mass of X-ray emitting halos. In all cases the value of star formation efficiency was kept fixed at f * = 2%. In comparison, we explore the popular case of heating by a realistic population of X-ray binaries with hard SED. In this case heating is inefficient and fluctuations are smoothed out ). Therefore, we expect our results to be closer to case (1) of Ghara et al. (2020). Moreover, in our work all the parameters (except for X-ray SED) are allowed to vary over a wide range, e.g., f * is varied between 0.1% and 50%.
Despite these differences in our modelling, qualitatively our work is consistent with Ghara et al. (2020). Both works rule out a cold IGM with an ionization fraction close to 50%. Namely, in their case (1)xH ii ∼ 0.2 − 0.75 and Tgas 3 K are disfavoured; while we find thatxH ii ∼ 0.38 − 0.72 and Tgas 10.1 K are disfavoured (at 68%).

CONCLUSIONS
In this paper we have used the upper limit on the 21-cm signal from z = 9.1 based on 141 hours of observations with LOFAR (Mertens et al. 2020) to evaluate the contribution of the high redshift Universe to the excess radio background over the CMB detected by ARCADE2 (Fixsen et al. 2011) and LWA1 (Dowell & Taylor 2018). Assuming synchrotron spectrum of the radio background with spectral index β = −2.6 and marginalising over the astrophysical properties of first star-forming sources, we find (at 68% C.L.) the contribution above the CMB level to be less than a factor of 13.2 at the reference frequency of 78 MHz, equivalent to 0.7% of the CMB at 1.42 GHz. This limit, for the first time, rules out strong contribution of the high-redshift Universe to the excess detected by ARCADE2 and LWA1. At the level below 0.7% of the CMB, the extra radio background could, on one hand, be strong enough to explain the tentative EDGES Low-Band detection which requires an excess of at least 0.1% of the CMB ). On the other hand, such a small contribution would be within the measurement error of the radio telescopes. Hence, it would remain a plausible explanation for the detected EDGES signal, even if the excess radio background measured by AR-CADE2 and LWA1 is due to an erroneous Galactic modelling (Subrahmanyan & Cowsik 2013).
We also use LOFAR data to constrain thermal and ionization state of the IGM at z = 9.1 in models with and without the extra radio background over the CMB. If such an extra radio background is present at z = 9.1, the fluctuations in the 21-cm signal are boosted compared to the standard case, which gives LOFAR a larger lever to reject models. Therefore, for the models with excess radio background, constraints on the astrophysical properties and the properties of the IGM are tighter than in the standard case. In particular, warmer IGM scenarios with mean neutral gas temperature of up to 19.2 K are disfavoured (at 68% C.L.) compared to only up to 10.1 K in the standard case.
Using the LOFAR data we have also derived 68% C.L. limits on the astrophysical parameters of Cosmic Dawn and EoR. The data disfavour very efficient star formation above 5%, imply the existence of small halos at early times (of masses below few×10 8 M ⊙ at z = 9.1), require the presence of X-ray sources, and disfavour a CMB optical depth above τ ∼ 0.076. For the suite of standard models, we point out that the LOFAR data rule out similar type of models as those rejected by the global signal experiments, namely the EDGES High-Band (Monsalve et al. 2019) and SARAS2 (Singh et al. 2018). Finally, we note that our constraints of the standard-physics parameters are in a qualitative agreement with the results reported by Ghara et al. (2020). A detailed comparison between these two works is beyond the scope of this paper.
Although other high redshift probes (e.g., the Planck measurement of the CMB optical, high redshift quasars and galaxies) allow to put tighter constraints on the ionization history, the 21-cm observations provide a unique way to probe the thermal history of the Universe, constrain properties of first star forming halos and test the nature of the radio background.  Figure 6 (right)