Joint analysis constraints on the physics of the first galaxies with low-frequency radio astronomy data

ABSTRACT

The first billion years of cosmic history remains largely unobserved. We demonstrate, using a novel machine learning technique, how combining upper limits on the spatial fluctuations in the 21-cm signal with observations of the sky-averaged 21-cm signal from neutral hydrogen can improve our understanding of this epoch. By jointly analysing data from the Shaped Antenna measurement of the background RAdio Spectrum (SARAS3, redshift z ≈ 15−25) and limits from the Hydrogen Epoch of Reionization Array (HERA, z ≈ 8 and 10), we show that such a synergetic analysis provides tighter constraints on the astrophysics of galaxies 200 million years after the big bang than can be achieved with the individual data sets. Although our constraints are weak, this is the first time data from a sky-averaged 21-cm experiment and power spectrum experiment have been analysed together. In synergy, the two experiments leave only |$64.9^{+0.3}_{-0.1}$| per cent of the explored broad theoretical parameter space to be consistent with the joint data set, in comparison to |$92.3^{+0.3}_{-0.1}$| per cent for SARAS3 and |$79.0^{+0.5}_{-0.2}$| per cent for HERA alone. We use the joint analysis to constrain star formation efficiency, minimum halo mass for star formation, X-ray luminosity of early emitters, and the radio luminosity of early galaxies. The joint analysis disfavours at 68 per cent confidence a combination of galaxies with X-ray emission that is ≲33 and radio emission that is ≳32 times as efficient as present day galaxies. We disfavour at 95 per cent confidence scenarios in which power spectra are ≥126 mK² at z = 25 and the sky-averaged signals are ≤−277 mK.

1 INTRODUCTION

The period of cosmic history covering the birth of primordial stars, formation of the very first black holes, and assembly of the earliest galaxies is largely unconstrained by observations. Covering the first billion years after the formation of the cosmic microwave background (CMB), the period is currently out of reach of the modern telescopes, but theoretical models suggest the first star formed between z ∼ 20 and 60, i.e. around 35−200 million years after the big bang (Bond 1981; Bromm & Larson 2004; Klessen 2019).

Observational prospects of this epoch are improving with the exciting new observations from the JWST (Windhorst et al. 2006; Park et al. 2020; Robertson 2022), which is expected to be able to observe bright galaxies out to redshift z ≈ 20 (e.g. Windhorst et al. 2006; Naidu et al. 2022; Hainline et al. 2023) and is much more sensitive than the previous generation of telescopes such as the Hubble Space Telescope. The highest spectroscopically confirmed galaxy observed with JWST at the time of writing is at z = 13.20 (Robertson et al. 2023). This picture of the infant Universe will potentially be supplemented by proposed X-ray observatories such as ATHENA (The Athena Collaboration 2013), AXIS (Mushotzky 2018), and LYNX (The Lynx Team 2018) that will probe X-ray sources and the hot diffused gas in the infant Universe, improving our understanding of some of the first X-ray emitting objects. In addition, radio telescopes such as the Square Kilometre Array (SKA; Mellema et al. 2013; Koopmans et al. 2015; Greig, Mesinger & Koopmans 2020) will provide information about the state of the Intergalactic Medium through tomographic observations of the 21-cm line (Furlanetto, Oh & Briggs 2006; Barkana 2016; Mesinger 2019) during the Cosmic Dawn and Epoch of Reionization (EoR). In the coming decades, through a joined effort over a range of different wavelengths, researchers will reveal the first billion years of cosmic history between the formation of the CMB and the current epoch of the galaxies.

In 2018, the EDGES collaboration reported a tentative detection of the sky-averaged (or global) 21-cm signal (Bowman et al. 2018) which led to an increased interest in the field. The reported absorption feature is much deeper than conventional theoretical models (Mesinger, Furlanetto & Cen 2011; Cohen et al. 2017; Mirocha, Furlanetto & Sun 2017; Reis, Fialkov & Barkana 2021; Muñoz et al. 2022) and implies a rapid period of efficient star formation followed by a delayed onset of strong heating of the neutral gas. Exotic processes such as cooling of the gas through interactions with charged cold dark matter (Barkana 2018; Barkana et al. 2018; Berlin et al. 2018; Fialkov, Barkana & Cohen 2018; Kovetz et al. 2018; Muñoz & Loeb 2018; Slatyer & Wu 2018; Liu et al. 2019; Li et al. 2021; Wu, Xu & Zheng 2022; Caputo et al. 2023) and models with excess radio backgrounds above the CMB (Ewall-Wice et al. 2018; Fialkov & Barkana 2019; Jana, Nath & Biermann 2019; Reis, Fialkov & Barkana 2020; Acharya, Dhandha & Chluba 2022; Mittal & Kulkarni 2022b; Acharya, Cyr & Chluba 2023) have been proposed to explain the depth of the absorption feature. However, explaining the implied rapid onset of star formation, delayed and then rapid heating is much harder and a number of works have suggested that there are unaccounted for systematics in the data (Hills et al. 2018; Bradley et al. 2019; Singh & Subrahmanyan 2019; Sims & Pober 2020; Bevins et al. 2021a).

A number of other experimental efforts are underway to try to detect the high-redshift 21-cm signal. Upper limits on the 21-cm power spectrum, which measures fluctuations in the 21-cm field of time and space, have been reported by PAPER (Jacobs et al. 2015), MWA (Trott et al. 2020; Kolopanis et al. 2023), LOFAR (Patil et al. 2017; Mertens et al. 2020), AARTFAAC (Gehlot et al. 2020), HERA (The HERA Collaboration 2022; Abdurashidova et al. 2022a, b), LWA (Eastwood et al. 2019), and LEDA (Garsden et al. 2021). Similarly, limits on the magnitude of the global 21-cm signal have been reported by several experiments including SARAS (Singh et al. 2017, 2018; Bevins et al. 2022c, d), EDGES High Band (Monsalve et al. 2017, 2018, 2019), and LEDA (Bernardi et al. 2016). The limits have been used to place constraints on the properties of the first galaxies (Monsalve et al. 2017; Singh et al. 2017; Monsalve et al. 2018; Singh et al. 2018; Monsalve et al. 2019; Mondal et al. 2020; The HERA Collaboration 2022; Abdurashidova et al. 2022a; Bevins et al. 2022c, d) and are improving significantly over time. Future space-based missions such as DAPPER and FARSIDE (Burns et al. 2021) as well as ongoing and upcoming experiments from Earth such as SARAS (Singh et al. 2022), MIST (Monsalve et al. 2023), REACH (de Lera Acedo et al. 2022), LOFAR (Koopmans 2017), NenuFAR (Zarka et al. 2018), HERA (DeBoer et al. 2017), the SKA all hope to further improve our understanding of the early Universe through the 21-cm line.

In this paper, we take upper limits from 18 nights of observations with the HERA interferometer on the magnitude of the 21-cm power spectrum at redshifts z ≈ 8 and ≈10. HERA provides a window to the EoR when the neutral hydrogen gas was efficiently ionized by ultraviolet emission from the first massive galaxies (Abdurashidova et al. 2022a), and the SARAS3 radiometer (Singh et al. 2022) that probes the sky-averaged 21-cm emission at redshifts between z ≈ 15−25. SARAS3 observes the Cosmic Dawn when the first small galaxies, containing the first stars and X-ray emitting objects, are expected to have formed. For the first time, we jointly analyse the upper limits from HERA and SARAS3 to improve our understanding about the properties of the Intergalactic Medium during the Cosmic Dawn and the properties of the first galaxies. While in the HERA analysis the authors explore a number of different models of the 21-cm signal, we focus here on the models that include an excess radio background above the CMB from high-redshift radio galaxies, Ly-α heating, CMB heating, multiple scattering of Ly-α photons and X-ray heating (see Reis, Fialkov & Barkana 2020, 2021). Since these models were also used in the SARAS3 analysis, it is convenient to use them here to showcase the power of the joint analysis. We refer to the HERA constraints on this specific 21-cm model as ‘the HERA results’ throughout the work but stress that the constraints are model dependent, and we make efforts to discuss constraints on derived parameters such as the radio background temperature and spin temperature allowing for comparisons across different models of the 21-cm signal where possible.

We introduce a machine learning enhanced Bayesian workflow to efficiently perform the joint analysis. When considered together, we find that the results from HERA and SARAS3 provide a better leverage on theoretical scenarios across a wider redshift range than the constraints from each of the individual experiment.

In Section 2, we review the data and the specifics of the experiments incorporated in our analysis. This is followed by a discussion of our machine learning enhanced Bayesian workflow and about the synergies between the power spectrum and sky-averaged 21-cm experiments. We present the implications of our work for the astrophysical constraints in Section 3. This is followed by conclusions in Section 4.

2 METHODOLOGY

2.1 Data

Observing from a lake in Southern India, the global 21-cm experiment SARAS3 took 15 h of usable data over the frequency band 55−85 MHz. After correcting for the impact of the environment, man-made radio signals and the efficiency of the antenna, a measurement of the average sky temperature was attained. The data are expected to include foregrounds from the Galaxy and other galaxies as well as the 21-cm signal. The data were used to constrain the properties of the first galaxies (Bevins et al. 2022c) and the collaboration reported a null detection of the EDGES absorption feature with 95.3 per cent confidence (Singh et al. 2022).

A previous iteration of the SARAS experiment, SARAS2, observed at higher frequencies, 110−200 MHz. The data contained systematic structures, thought to be from the environment, but was used to place constraints on the magnitude of the 21-cm signal and properties of early galaxies (Singh et al. 2017, 2018; Bevins et al. 2022c). In this work, we combine constraints from Bevins et al. (2022c) with constraints from SARAS3 and HERA.

The best upper limits observed by interferometric experiments to date have come from HERA interferometer (Abdurashidova et al. 2022a), followed by MWA in the redshift range z = 6.5−8.7 (Trott et al. 2020) and LOFAR, which provide the tightest upper limits at redshifts z ∼ 9.1 (Patil et al. 2017) and z ∼ 9.3−10.6 (Mertens et al. 2020). HERA is located in South Africa in the Karoo Desert (DeBoer et al. 2017) and in this work we use spherically averaged power spectra upper limits from Internal Data Release 2 (Abdurashidova et al. 2022b), which include observations over 18 nights with 39 antennas. The data cover the wavenumber range |$k=0.128\, h\mathrm{Mpc}^{-1}$| to |$0.960\, h\mathrm{Mpc}^{-1}$| and from the two bands focusing on redshifts z = 7.9 and 10.3. We leave the inclusion of constraints from the analysis of Internal Data Release 3 that was recently published (The HERA Collaboration et al. 2022) to future work, noting that there is an overlap between the data sets. In Section 3.7, we consider the implications of including the upper limits on the power spectrum from MWA and LOFAR in our joint analysis with SARAS3 and HERA.

2.2 Machine learning enhanced Bayesian workflow

In 21-cm cosmology, data analysis efforts are increasingly employing Bayes theorem:

to derive constraints on the astrophysical scenarios of the early Universe. Here the likelihood, |$\mathcal {L}(\theta)$|⁠, represents the probability that we observe the data, D, from SARAS3 or HERA, given a particular model, |$\mathcal {M}$|⁠. The prior, π(θ), represents our assumed knowledge before we consider any data and the evidence, |$\mathcal {Z}$|⁠, is a normalization constant that can be used for model comparisons. The posterior, |$\mathcal {P}(\theta |D, \mathcal {M})$|⁠, tells us which parts of the parameter space, θ, given the data and chosen model, are more probable than others. In many 21-cm experiments, θ is composed of parameters that describe instrumental effects, θ_I, foregrounds, θ_fg, and the astrophysical processes that influence the 21-cm signal, θ₂₁. We typically refer to the set of θ_I and θ_fg as the nuisance parameters. Since we are only interested in the 21-cm signal, we work with the marginal or nuisance-free posteriors and nuisance-free likelihoods, |$\mathcal {L}(\theta _{21})$|⁠, which can be estimated using normalizing flows and the marginal Bayesian statistics code margarine (Bevins et al. 2022a,b).

margarine uses neural networks known as masked autoregressive flows (MAFs) to model the probability distribution, |$\mathcal {P}(\theta |D, \mathcal {M})$|⁠, of a set of samples, here θ₂₁, marginalizing over nuisance parameters describing the instrumental systematics, θ_I, and foregrounds, θ_fg in the process. It does this by shifting and scaling a base standard normal distribution, making the probability of the distribution easily tractable, to replicate the target samples, where the shifting and scaling are determined by the outputs of the MAFs.

This can subsequently be used to calculate a nuisance-free likelihood given a prior distribution and the Bayesian evidence using

where α = {θ_I, θ_fg} (Bevins et al. 2022b). In instances when the prior is flat then the following is true up to a normalization constant

With margarine we can then evaluate |$\log (\mathcal {L}(\theta _{21}))$| for any set of θ₂₁ for any existing posterior distribution, |$\mathcal {P}(\theta |D, \mathcal {M})$|⁠, from previous analysis of a data set like HERA or SARAS3

where the log  is base 10. The likelihood evaluations can then be combined,

to be sampled over using Markov chain Monte Carlo (MCMC) methods or in our case Nested Sampling implemented with polychord (Handley, Hobson & Lasenby 2015a,b).

MCMC sampling methods approximate the unnormalized posterior distribution by directly sampling the product |$\mathcal {P}(\theta) \approx \mathcal {L}(\theta)\pi (\theta)$| and do not provide estimates of the evidence, |$\mathcal {Z}$|⁠. The family of algorithms typically use random walkers to traverse the parameter space, with points accepted and rejected based on some probabilistic criteria in an effort to explore the space fully. The HERA analysis (Abdurashidova et al. 2022a) of the excess radio background models explored here used the emcee (Foreman-Mackey et al. 2013) implementation of MCMC sampling (Abdurashidova et al. 2022a).

Nested sampling (Skilling 2004) numerically approximates the integral

which can be derived from equation (1) and the requirement that the posterior must integrate to 1, by evolving a series of live points to higher and higher likelihood values. Through approximating the evidence, the algorithm produces samples on the normalized posterior distribution, which we subsequently use to determine preferred and disfavoured regions of the parameter space. A comprehensive review of Nested Sampling can be found in (Ashton et al. 2022).

2.3 Signal Modelling

In order to realistically model the range of time covered by the Cosmic Dawn and Epoch of Reionization, we need a consistent modelling of the cosmological and astrophysical processes from redshift 60, when star formation might have begun, all the way to redshift 5 at the end of the EoR. To that end we employ seminumerical simulations (Visbal et al. 2012; Fialkov et al. 2013; Fialkov & Barkana 2014; Fialkov, Barkana & Visbal 2014; Reis, Fialkov & Barkana 2021) that have previously been used in the HERA (Abdurashidova et al. 2022a), SARAS2 (Bevins et al. 2022d), and SARAS3 (Bevins et al. 2022c) analyses as well as similar analysis of the LOFAR (Mondal et al. 2020) and EDGES high-band (Monsalve et al. 2019) data. We model the three-dimensional 21-cm field as a function of cosmic time during the infant Universe taking into account important astrophysical process including the Wouthuysen–Field (WF) effect (Wouthuysen 1952; Field 1959; Fialkov & Barkana 2014), heating of the intergalactic medium by X-ray (Fialkov, Barkana & Visbal 2014), Ly-α (Reis, Fialkov & Barkana 2021) and CMB photons (Fialkov & Barkana 2019), multiple scattering of Ly-α (Reis, Fialkov & Barkana 2021), relative velocity between dark matter and gas (Visbal et al. 2012), feedback of Lyman-Werner radiation on star formation (Fialkov et al. 2013), and radio emission from galaxies (Reis, Fialkov & Barkana 2020).

At high redshifts around z ≈ 20−30, the first stars begin to form and produce Ly-α photons that interact with the baryonic matter, predominantly composed of neutral hydrogen, in the Universe. Neutral hydrogen atoms absorb and re-emit ambient Ly-α photons in a process known as the Wouthuysen–Field (WF) effect (Wouthuysen 1952; Field 1959), introduced above, that drives the relative number of atoms with aligned and anti-aligned proton and electron spins. This process couples the spin temperature of the neutral hydrogen, T_s, which describes the distribution of hydrogen atom spins, to the gas temperature, T_k, which is cooling at a faster rate than the radio background as the Universe expands. Further, interactions between the neutral hydrogen and Ly-α emission result in the transfer of kinetic energy that raises the gas temperature, and coupled spin temperature, in a process known as Ly-α heating (Madau, Meiksin & Rees 1997; Chuzhoy & Shapiro 2007) preventing the gas from cooling adiabatically. However, despite the heating, the gas temperature remains cooler than the radio background at high redshifts, leading to an absorption feature in the sky-averaged 21-cm signal. Further, it leads to a peak in the power spectrum at high z ≈ 25 on angular scales corresponding to the effective horizon of the Ly-α emission and the distribution of galaxies, which disappears when the coupling becomes saturated (Cohen, Fialkov & Barkana 2018). The intensity and spatial fluctuations of the Ly-α emission evolve with the population of early galaxies and consequently it is dependent on their star formation rate and the minimum halo mass for star formation. We parametrize these quantities with the star formation efficiency, f_*, which quantifies the percentage of the baryonic mass in the star forming haloes that is converted into stars, and the minimum circular velocity, V_c, which is proportional to the cube root of the minimum halo mass for star formation, M.

At intermediate redshifts of z ≈ 10−20 the gas is further heated by X-ray binaries (Fragos et al. 2013; Fialkov, Barkana & Visbal 2014), continuing Ly-α heating, CMB heating (Venumadhav et al. 2018), and heating through structure formation. X-ray heating is dependent on processes such as black hole binary formation and X-ray production in high redshift galaxies. This means that it has to be separately parametrized and in our seminumerical simulations, we model the X-ray production efficiency, f_X, which is directly proportional to the X-ray luminosity per star formation rate, L_X/SFR measured in erg s⁻¹ M|$_\odot ^{-1}$| yr, between 0.2 and 95 keV. The heating affects the redshift of the minimum and depth of the sky-averaged 21-cm signal. If sufficiently efficient it can raise the temperature of the gas, and, consequently, the coupled 21-cm brightness temperature, above the radio background resulting in an emission feature in the sky-averaged signal at low redshifts. The heating rate has a direct impact on how fast the brightness temperature transforms from absorption to 0 K or emission. In the power spectrum, due to the non-uniformity of heating, the various mechanisms can produce a peak in the signal at around z ≈ 15. Although this redshift is model dependent (see Cohen, Fialkov & Barkana 2018) and in some cases when heating is done by hard X-rays (energies >1 keV with long mean free paths) X-ray heating is smooth and no peak is imprinted in the power spectrum (Fialkov & Barkana 2014; Fialkov, Barkana & Visbal 2014).

Finally, at more recent times, z ≈ 5−15, ultraviolet emission from the first massive galaxies begins to ionize the neutral hydrogen, stripping the abundant gas of its electrons. This reduces the sky-averaged 21-cm signal and when reionization is complete, the signal disappears. The process produces a peak in the power spectrum at the scales corresponding to the typical size of ionized bubbles, but again the signal is destroyed once reionization is complete. The process of reionization is highly dependent on the ionizing efficiency of sources, ζ, which in the models explored here is normalized by the CMB optical depth, τ. While τ has been weakly constrained by cosmological experiments such as Planck (Planck Collaboration VI 2020), we treat it as a free parameter and note that 21-cm cosmology offers a means by which to break degeneracies between τ and other cosmological parameters.

Throughout our analysis, we explore a broad range of radio luminosities that produce radio excesses above the CMB of between ≈0.5 and 270 times the CMB temperature at z = 20 and ≈1−32000 times the CMB at z = 10. While moderate radio emission might be expected (Mirocha & Furlanetto 2019), the extreme values of radio production efficiencies are not expected in nature and are usually invoked to explain the anomalously deep EDGES signal (Bowman et al. 2018; Ewall-Wice et al. 2018; Feng & Holder 2018; Fialkov & Barkana 2019; Jana, Nath & Biermann 2019; Mirocha & Furlanetto 2019; Reis, Fialkov & Barkana 2020) but struggle to explain the rapid star formation and rapid heating of the gas (Mittal & Kulkarni 2022a) that is implied by the shape of the EDGES signal.

To summarize, the key parameters in the astrophysical model are the star formation efficiency, f_*, which quantifies the percentage of the baryonic mass in the star-forming haloes that is converted into stars; the minimum virial circular velocity, V_c, which is proportional to the cube root of the halo mass, M; the X-ray production efficiency, f_X, which is directly proportional to the X-ray luminosity per star formation rate, L_X/SFR, measured in erg s⁻¹ M|$_\odot ^{-1}$| yr, between 0.2 and 95 keV; the CMB optical depth, τ; finally, the radio production efficiency, f_radio, of high-redshift radio galaxies which is proportional to the radio luminosity per star formation rate, L_r, measured in W Hz⁻¹ M|$_\odot ^{-1}$| yr at 150 MHz,

and normalized such that it has a value of one for the present day population of radio galaxies (Reis, Fialkov & Barkana 2020). In our Bayesian analysis, θ₂₁, therefore, comprises the set of parameters {f_*, V_c, f_X, τ, f_radio} or equivalently {|$f_*, M, L_X/\rm {SFR}, \tau , L_\mathrm{r}/\rm {SFR}$|}.

We explore wide prior ranges on all the parameters in an attempt to let the data inform us about the high-redshift astrophysical processes. Specifically, we constrain f_* in the range of 0.001−0.5, V_c in the range of 4.2−100, f_X in the range of 10⁻⁴−1000, τ in the range of 0.035−0.077 and finally f_radio in the range of 1−99750. Typically, the radio background is assumed to be equal to the CMB and contributions from radio-luminous galaxies are not conventionally considered. Here, we expect that early galaxies will contribute to the radio background, thus increasing the amplitude of both the sky-averaged 21-cm signal (Feng & Holder 2018) and the power spectrum (Fialkov & Barkana 2019; Reis, Fialkov & Barkana 2020).

Both the sky-averaged 21-cm signal and the power spectrum rely on the same underlying physics, and constraints from experiments targeting the different probes can be effectively combined to improve our knowledge of the infant Universe.

2.4 Emulators

The seminumerical simulations used to model the signal take of order a few hours per parameter set, which is impractical for Nested Sampling or MCMC runs which need of order millions of signal evaluations (the SARAS3 constraints used in this work required ≈15 million signal evaluations to derive). However, as both the sky-averaged signal and the power spectrum change smoothly when we vary parameters, we can interpolate values at intermediate parameters from existing simulations. An increasingly common practice to achieve this are the fast and precise neural network based emulators. Typically, these networks take of order a few tens of milliseconds to evaluate meaning they are much more well suited for computationally intensive fitting algorithms. The specific emulators used in this analysis are detailed below.

In practice, it is important to include a theory error associated with the accuracy of our emulators in the likelihood functions and the uncertainty in the theoretical modelling (Greig & Mesinger 2017). An emulator error was included in the HERA likelihood in the original analysis (Abdurashidova et al. 2022a) and is done when analysing data from MWA and LOFAR in Section 3.7 but not done in the original SARAS2 and SARAS3 analysis. Since we are using the posterior distributions from the original HERA, SARAS2, and SARAS3 analysis and density estimators to produce our marginal or nuisance-free likelihood functions, we are unable to include the emulator error in the SARAS2 and three likelihoods here. However, the error in the global signal emulator is lower than the noise in the SARAS3 and SARAS2 data and the noise is expected to dominate over the both the uncertainty in the theoretical modelling and the emulator error. Theory errors should be included in future analysis, particularly as limits on the magnitude of the power spectrum and global signal improve, before performing joint analysis.

2.4.1 The sky-averaged 21-cm signal

To emulate the sky-averaged 21-cm signal, we use the publicly available emulator globalemu (Bevins et al. 2021b) trained on sets of astrophysical simulations.

The emulator has previously been used in the individual analysis of data from SARAS2 (Bevins et al. 2022d) and SARAS3 (Bevins et al. 2022c) and detailed in the corresponding papers. We therefore only briefly summarize the accuracy of the neural network here. The original set of simulations comprise approximately 10 000 models however the explored range of τ, the CMB optical depth, is large and when training our neural network we filter out models that have a value of τ in the range given by Planck ±3σ. This results in a training set of approximately 4300 models and the network is tested on approximately 500 models. For the test data set, a root mean squared error (RMSE) of 5.11 mK is found and a 95 percentile RMSE of 20.53 mK indicating a high level of accuracy (see table M.2 in Bevins et al. 2022c).

2.4.2 The 21-cm power spectrum

For comparing models with the HERA measurements we use the same 21-cm power spectrum emulator used in the HERA analysis (Abdurashidova et al. 2022a), trained on the same suit of simulations that give rise to the sky-averaged signal used in the SARAS2 and SARAS3 analysis. Based on an input of the five model parameters, this emulator returns the 21-cm power spectrum with a relative accuracy of 20  per cent at the wave numbers and redshifts observed by HERA. For HERA the impact of τ is more significant because the data corresponds to lower redshifts and so the full prior range on the parameter is used for training. This results in a training set of ∼8000 simulations and another 2000 independent samples for testing. Full details and tests of this emulator can be found in the HERA analysis paper (Abdurashidova et al. 2022a, appendix B).

We note that when performing our joint analysis we use the narrow prior on τ defined by Planck, since the SARAS3 posterior is not defined in the original analysis for values outside this range.

2.4.3 Temperatures

To emulate physical properties such as the spin temperature and temperature of the radio background (as seen in Fig. 1) we use a globalemu-style emulator. In this framework, the emulator takes in the five astrophysical parameters and a single redshift and returns a single corresponding temperature. In practice, this means that vectorized calls have to be made to emulate the spin temperature, T_s(z), and the radio background temperature, T_r(z), as a function of redshift, but the method is found to be more accurate, quicker and allows for interpolation at a range of different redshifts. We use the same training and test sets as used for the 21-cm power spectrum emulator, and a similar architecture of four layers, with a reduced size of 100, 30, 10, and 5 nodes per layer. This allows us to emulate the spin temperature T_s within |$\pm 6~{{\ \rm per\ cent}}$| accuracy (95 per cent confidence interval), and the radio background temperature T_r within |$\pm 4~{{\ \rm per\ cent}}$| accuracy (95 per cent confidence interval).

3 RESULTS

Although the constraints presented here are weak and model dependent, through the novelty of combining the previously reported HERA and SARAS3 constraints we produce the tightest constraints to date on the properties of the infant Universe, as detailed below. This is the first time a joint analysis between a global signal data and interferometric limits has been attempted.

3.1 The ratio of T_s and T_r

To visualize the importance of combining the constraining power of HERA and SARAS3 we show, in the top panel of Fig. 1, constraints on the ratio of the spin temperature of neutral hydrogen, T_s, and background radiation temperature, T_r. The background radiation temperature is the sum of the CMB temperature and radio background from galaxies. The ratio determines the maximum absorption of the sky-averaged 21-cm signal, the smaller the ratio the larger the signal can be. The SARAS3 limits on the 21-cm signal (grey markers) correspond to lower limits on T_s/T_r, with the corresponding 1 and 2σ contours shown as lines and extrapolated out of the SARAS3 band. Our simulations provide a natural link between the power spectra and the global quantities, e.g. T_s/T_r, meaning that we can use the limits on the power spectrum from HERA to derive an equivalent constraint on T_s/T_r. These constraints are shown in Fig. 1 (blue markers and lines). The joint constraint, as shown by the green-shaded regions, provides the strongest constraints on the ratio, and in particular at redshifts z = 15−20, gives better constraints than either experiment alone. Further, the combination of the two experimental data sets improves the constraints at intermediate redshifts over a pure extrapolation of each one of the sets of constraints. As a guideline, the dashed black line in the figure shows the ratio for T_r = T_CMB, i.e. in the absence of radio emission from galaxies, and assuming adiabatically cooled gas in an expanding Universe in the absence of any astrophysical heating sources but with saturated coupling between the 21-cm spin temperature and the gas kinetic temperature. This limit is often used in the literature to give context to the constraints (e.g. The HERA Collaboration et al. 2022; Abdurashidova et al. 2022a). However, we note that for the models tested here the gas does not cool adiabatically because of the CMB and Ly-α heating at early times, and we have an excess radio background above the CMB.

3.2 Functional constraints on Δ² and T₂₁

By taking the samples of θ₂₁ output from our fits and using the neural network emulators to produce corresponding global signals and power spectra we can look at constraints on the temperature of the 21-cm signal and on the power spectrum. We explore the functional constraints in the T₂₁ − z and Δ² − z planes as shown in Fig. 2. We see that for both the power spectra and the sky-averaged signals, although more clearly for the latter, the range of plausible models is reduced by our joint analysis. The signals that are inconsistent with the data typically have strong power spectra and a corresponding deep absorption trough in the global signal, and belong to scenarios with weak X-ray heating and strong radio luminosity. Further, we see that the 2σ region for the functional constraints corresponds to signals with the power spectra ≲10^2.1 mK² at z = 25, which via our modelling maps to global 21-cm signals shallower than −277 mK. We note that the 3σ limit on the power spectrum at z = 25 of ≲10^3.2 mK² is approximately equivalent to the expected sensitivity of NenuFAR from 1000 h of observations (Mertens, Semelin & Koopmans 2021). For the global 21-cm signal 3σ limit on the magnitude reduces from ≈−2630 mK from HERA to ≈−1770 mK from the joint analysis. Remarkably, we find that for the sky-averaged signal, the 2σ limit is very close to the minimum depth of the ‘standard astrophysical’ models, ≈−165 mK (Reis, Fialkov & Barkana 2021), where the radio background is equated to the CMB, the contributions from radio galaxies and X-ray heating sources are assumed to be negligible, while the CMB and Ly-α heating are present. It is clear from our analysis, that the joint constraint improves limits on the range of plausible global signals and power spectra.

3.3 Parameter constraints

We now interpret the constraints on the temperatures in terms of constraints on the properties of the first galaxies. The key constraints from the joint analysis are shown in the bottom panels of Fig. 1 including the limits on the high-redshift star formation that drives the high-redshift portion of the 21-cm signal via the process of Ly-α coupling (constraints in the plane V_c − f_*) and the constraints on the luminosity of X-ray and radio sources (L_r − L_X) that primarily regulates the depth of the absorption trough. The full marginalized 1D and 2D posteriors corresponding to the joint analysis are shown in Fig. 3 and the key numerical results are summarized in Table 1 alongside the individual constraints from SARAS3 and HERA. In the Appendix, we show the individual constraints from each experiment and the joint constraint on L_r − L_X and V_c − f_* for comparison. We find that the combination of the two experiments leads to stronger constraints in the two-dimensional probability distribution of L_r − L_X than either of the two experiments individually. Where HERA constrains the population of high-redshift radio-luminous galaxies to be ≲400 times brighter in the radio band than the current population, the combination of the data sets constrains the galaxies to be ≲ 300 times brighter when marginalizing over the other parameters (V_c, f_*, f_X, and τ). Similarly, HERA disfavours at 68 per cent confidence galaxies with an X-ray luminosity ≲0.25 times the present-day value in combination with the radio luminosity of galaxies in the early universe that is ≳400 times the present-day value. The joint analysis provides a stronger constraint, ruling out scenarios where the X-ray luminosity is ≲33 times the present-day value and the radio luminosity of the first galaxies is ≳32 times the present-day value at 68  per cent confidence.

We find comparable constraints on f_* and minimum mass of star forming haloes, M, as was found with SARAS3 alone (Bevins et al. 2022c) when combining the data sets. Since these parameters predominantly affect the high-redshift signal during the Cosmic Dawn they are better constrained by SARAS3, and we do not gain any information about star formation from HERA, which observes at z ≈ 8 and 10 during the Epoch of Reionization. Marginalizing over the radio and X-ray luminosities, we disfavour at 68 per cent confidence galaxies in which |$\gtrsim 2~{{\ \rm per\ cent}}$| of gas is converted into stars and the minimum mass for star forming haloes is ≲ 45 million solar masses.

3.4 Constraints on the radio and X-ray backgrounds

The product f_*f_radio is proportional to the total radio background created by radio-luminous galaxies, and, equivalently, f_*f_X is a proxy for the total X-ray background created by the early population of X-ray sources. This is demonstrably true for our model parametrization as the star formation rate is proportional to f_* and the radio luminosity and X-ray luminosities per star formation rate are proportional to f_radio and f_X, respectively. Both X-ray and radio backgrounds are responsible for regulating the depth of the absorption feature, and they can also be observed independently by other telescopes (e.g. observations of the unresolved X-ray background by Chandra (Lehmer et al. 2012) and of the low-frequency radio background by ARCADE2/LWA (Fixsen et al. 2011; Dowell & Taylor 2018). In Fig. 4, we show constraints on the values of f_* f_X and f_* f_radio achieved by SARAS3, HERA and the joint analysis. Since these combinations of the parameters regulate the absorption depth of the global 21-cm signal, we can also condition our prior on the astrophysical parameters to produce signals with the same central frequency and depth as the absorption feature found in the EDGES data (Fialkov & Barkana 2019; Reis, Fialkov & Barkana 2021; Bevins et al. 2022c). In each panel of Fig. 4, we show black contours corresponding to these EDGES-like physical signals. We see that while HERA and SARAS3 allow for combinations of f_X f_* and f_radio f_* that could partially explain EDGES, the combination of the two experiments, which produces a tighter constraint on the X-ray and radio luminosities of early galaxies, disfavours a large portion of the EDGES-like parameter space (i.e. most of the EDGES-like parameter space is beyond the 95  per cent contours of the joint constraints, while it is well within the 95  per cent contours for SARAS3 and HERA individually). This demonstrates further the power of combining different data sets. However, we note that the explored theoretical signals do not fit EDGES data well as none of them closely reproduces the flattened Gaussian-like feature found in the data (Bowman et al. 2018).

3.5 Constraints on phemenological parameters

We explore the structure of the global 21-cm signals which are consistent with the data from SARAS3, HERA and the two sets together. In order to do this, we look at the minimum temperature, T_min, the corresponding central redshift, z₀, and an approximate full width at half-maximum, Δz, of each absorption trough as defined in the top right of Fig. 5. More specifically, for each parameter set θ₂₁ (either from the prior parameter range, or from the posterior ranges consistent with SARAS3, HERA, and the joint analysis) we generate a global signal using the neural network emulator globalemu (Bevins et al. 2021b). We then measure the values of T_min, z₀, and Δz for each signal producing probability distributions for each parameter corresponding to the constraints from each data set and showing the extent of the prior. We compare these distributions with the values used to parametrize the phenomenological EDGES flattened Gaussian signal with the EDGES 99 per cent confidence ranges shown by the black crosses in the corner plot in Fig. 5.

When considered individually, the SARAS3 and HERA experiments allow for astrophysically motivated sky-averaged 21-cm signals that have a minimum temperature, location and width in agreement with the EDGES detection, while the joint analysis rules out models with a depth that is consistent with EDGES at greater than a 2σ significance. The joint analysis has a preference for shallower (lower values of |T_min|) and narrower signals with higher central frequencies, as can be seen in the corresponding 1D PDFs, which is driven largely by the HERA constraints.

3.6 The impact of SARAS2

Previous analysis of the SARAS2 data revealed some weak constraints, most notably in the plane of L_X − L_r in agreement with HERA and SARAS3, on the properties of galaxies in the infant Universe (Bevins et al. 2022d). SARAS2 is at much lower redshifts than SARAS3 but overlaps with the redshifts probed by HERA having recorded observations in the band z ≈ 7−12. The data are contaminated by a sinusoidal systematic and a number of different models were fitted to this feature. The systematic models correspond to a signal introduced prior to the antenna possibly from ground emission or some unknown component of the foreground, and separately a signal introduced in the system electronics potentially from cable reflections. The sinusoidal systematic was fitted alongside a signal model generated with globalemu and a foreground model that is conditioned to be smooth preventing it fitting out non-smooth systematics or signals in the data (Bevins et al. 2021a). Here, we take the best-fitting model, with a systematic from ground emission or a non-smooth component from the foreground, with the highest evidence from the original analysis (Bevins et al. 2022d) and combine the corresponding constraints on the astrophysical parameters V_c, f_*, f_X, f_radio, and τ with the joint constraints from HERA and SARAS3 to assess the impact.

As can be seen in Fig. 6, the addition of SARAS2 to our analysis washes out the constraint in the plane f_* − V_c. One possible explanation for this is that the addition of SARAS2, while constraining the properties that affect the signal at low redshifts, increases the envelope of possible models at higher redshifts, where star formation is more important, that are plausible even given the constraints from the SARAS3 data. Despite this, we note that we maintain the constraint in the plane L_X − L_r when we add SARAS2 into our analysis.

We can quantify the impact of SARAS2 on our analysis by looking at the percentage of the astrophysical prior volume which is consistent with the different combinations of the three different data sets. To calculate this percentage, we use margarine to calculate the marginal Kullback–Leibler (KL) divergence (Kullback & Leibler 1951), |$\mathcal {D}$|⁠, between the flat prior on the five astrophysical parameters in the set θ₂₁ and the corresponding posteriors. The KL divergence gives a measure of the information gained when moving from the prior to the posterior and is related to the percentage via

where V_π is the prior volume and |$V_\mathcal {P}$| is the posterior volume. This quantity is useful as it quantifies the constraining power of the different data sets in all five dimensions, including correlations that may not be visible in the one and two-dimensional projections used to produce the corner plots in this paper and in the literature.

We show in Fig. 7 the percentage of the astrophysical parameter prior volume that is consistent with different combinations of the data sets discussed in this work (including additional interferometric measurements of the power spectrum discussed in Section 3.7). We see that the combination of either or both of the SARAS data sets with HERA lead to a percentage consistency with the data of |$\approx 63 - 65~{{\ \rm per\ cent}}$| and this is likely dominated by HERA. Individually, HERA allows for |$\approx 80~{{\ \rm per\ cent}}$| of the astrophysical parameter space, SARAS2 for |$\approx 90~{{\ \rm per\ cent}}$| and SARAS3 for |$\approx 92~{{\ \rm per\ cent}}$|⁠. Due to the uncertainty in the modelling of the systematics in the SARAS2 analysis, we leave SARAS2 out of the main results.

3.7 The impact of MWA and LOFAR

In Fig. 8, we show the projected posteriors derived using HERA data alone (left panel) and HERA, MWA and LOFAR together (right panel). We note that the constraints from the different interferometers are all at low redshifts between z ≈ 6 and 10 and varying wavenumbers or angular scales. These are detailed in Table 2 along with the constraints from the individual experiments on key parameters.

We derived parameter constraints from the MWA and LOFAR data using the approach taken in the orginal HERA analysis. Specifically, we take the measured upper limits, the mean power spectrum and uncertainty, and treat it as a measurement of cosmological signal plus systematics. As in HERA (Abdurashidova et al. 2022a) we take this uncertainty to be Gaussian and marginalize over a uniform prior on the systematics, yielding the likelihood

where N_d represents the number of data points, d_i and σ_i correspond to the mean and an error term which includes a contribution from the standard deviation in each data point and an emulator error, and m_i(θ₂₁) is the model prediction for that redshift and wave number. Thus, a model prediction m ≫ d gives |$\mathcal {L}\approx 0$| while m ≪ d gives a constant. This likelihood is effectively a step function that disfavours models above a given amplitude. A full discussion of its derivation can be found in Abdurashidova et al. (2022a).

We performed this joint analysis using a full analytic likelihood approach (independent of margarine) since there are no nuisance parameters describing the systematics. We find that each of the experiments disfavours individually similar regions of the L_r−L_X plane. However, the joint analysis does not improve the results derived from HERA data alone (as we summarize in Table 2 and is further illustrated in Fig. 7) that motivates our decision to only use HERA in our main results with SARAS3.

4 CONCLUSION

Through a combination of constraints on fluctuations and the sky-averaged 21-cm signal of neutral hydrogen, we have improved our understanding of the first galaxies that formed in the infant Universe between 200 and 700 million years after the big bang. This is the first time the data from the two different 21-cm probes have been combined to derive constraints on the astrophysical properties of the early galaxies. Even though the existing constraints are weak and model dependent, we develop novel methodology and outline the approach, which will become increasingly more useful as the next-generation experiments deliver stronger observational constraints.

Considering a wide space of plausible astrophysical models including high-redshift sources of ultraviolet, X-ray and radio photons which affect the 21-cm signal, we calculate corresponding sky-averaged spectra as well as the power spectra of fluctuations. Using an upper limit on the fluctuations from the HERA interferometer and non-detection of the global 21-cm signal by the SARAS3 radiometer, we find that only |$64.9^{+0.3}_{-0.1}$| per cent of the explored theoretical parameter space is consistent with the joint SARAS3 and HERA constraint, which is a significant improvement over the individual values of |$92.3^{+0.3}_{-0.1}$| per cent and 78.7 ± 0.2 per cent, respectively.

Using the newly developed methodology, we place the tightest current constraints on the properties of cosmic gas, such as the spin temperature of the 21-cm hydrogen line (closely related, but not equal, to the gas kinetic temperature) and the radio background temperature, T_r, as well as on the radio and X-ray luminosities of the first galaxies disfavouring at 68 per cent confidence galaxies that are approximately 32 times more efficient radio emitters than present galaxies and simultaneously are less than 33 times bright in the X-ray band. This work reports an increased degree of confidence over a wider range of redshifts than previous works which typically extrapolate outside the redshifts targeted by individual experiments, while here we interpolated between the observations of SARAS3 at z = 15−25 and the HERA limits at lower redshifts z ∼ 8 and 10.

In this work, we also considered the addition of interferometric data sets from MWA and LOFAR to our analysis, which, owing to the current weaker limits reported by these experiments, led to a negligible improvement in the results. Similarly to HERA, these experiments probe the physics of the EoR covering a similar redshift range to HERA (between z ≈ 6 and 10), with the current HERA data providing the tightest constraints on our models. Of the current global or sky-averaged 21-cm experiments only SARAS2, SARAS3, and EDGES were able to place limits on the astrophysics of the infant Universe. Our main focus is on the SARAS3 limits (although we also consider SARAS2), as there are concerns surrounding the cosmological nature of the signal reported in EDGES and a degree of uncertainty in the modelling of systematics in the SARAS2 data. We note that the SARAS2 data covers a similar redshift range as the HERA data and based on our analysis does not lead to a significant improvement in our constraints.

The new methodology developed in this paper will allow for synergies between the upcoming observations, e.g. of the power spectrum from cosmic dawn measured by the NenuFAR (Zarka et al. 2018), HERA, LOFAR, or the SKA (Mellema et al. 2013), as well as the measurements of the global signal by the wide-band REACH experiment covering the redshift range z ≈ 7−28 (de Lera Acedo et al. 2022), PRIZM (Philip et al. 2019), MIST (Monsalve et al. 2023), and missions to the moon (Burns et al. 2021) among others.

ACKNOWLEDGEMENTS

HTJB acknowledges the support of the Science and Technology Facilities Council (STFC) through grant number ST/T505997/1. IAC acknowledges support from the Institute of Astronomy Summer Internship Program at the University of Cambridge. WJH and AF were supported by Royal Society University Research Fellowships. EdLA was supported by the STFC through the Ernest Rutherford Fellowship. RB acknowledges the support of the Israel Science Foundation (grant no. 2359/20), The Ambrose Monell Foundation and the Institute for Advanced Study.

DATA AVAILABILITY

The SARAS3 and SARAS2 data are available upon reasonable request to SS. The HERA data is publicly available at https://reionization.org. The MWA and LOFAR data were provided to IAC through private correspondence with the relevant authors.

margarine is available at https://github.com/htjb/margarine. The semi-analytic simulations of the 21-cm signal discussed in Section 2 are not publicly available. globalemu is available at https://github.com/htjb/globalemu and the power spectrum emulators used are not publicly available but are built with sklearn (https://scikit-learn.org/stable/install.html). polychord is available at https://github.com/PolyChord/PolyChordLite. The graphs in the paper were made with anesthetic that is available at https://github.com/handley-lab/anesthetic and fgivenx which is available at https://github.com/handley-lab/fgivenx.

References

APPENDIX: INDIVIDUAL CONSTRAINTS FROM HERA AND SARAS3 ON KEY PARAMETERS

In Fig. A1, we show the individual constraints from HERA and SARAS3 along with the joint constraints on the two constrained parts of the parameter space L_X versus L_r and V_c versus f_*. All other parameters have been marginalized out in these plots.

From the figure, we see that HERA dominates the constraints on the strength of the radio and X-ray backgrounds. However, SARAS3 provides disfavours at 68 per cent confidence higher L_X in combination with high L_r than HERA leading to a tighter constraint on the parameter space in the joint analysis.

For the properties governing star formation, namely f_* and V_c, we see that SARAS3 dominates the joint analysis. The SARAS3 observations are at higher redshifts than HERA and cover a period of Cosmic history where the 21-cm signal is more sensitive to the processes of star formation. The apparent constraint from HERA on the star formation properties is washed out by SARAS3 since SARAS3 is more constraining.

The figure shows how the two probes are complimentary, providing additional information because they cover different redshifts dominated by different physical processes.