HERA Phase I Limits on the Cosmic 21-cm Signal: Constraints on Astrophysics and Cosmology During the Epoch of Reionization

Recently, the Hydrogen Epoch of Reionization Array (HERA) collaboration has produced the experiment's first upper limits on the power spectrum of 21-cm fluctuations at z~8 and 10. Here, we use several independent theoretical models to infer constraints on the intergalactic medium (IGM) and galaxies during the epoch of reionization (EoR) from these limits. We find that the IGM must have been heated above the adiabatic cooling threshold by z~8, independent of uncertainties about the IGM ionization state and the nature of the radio background. Combining HERA limits with galaxy and EoR observations constrains the spin temperature of the z~8 neutral IGM to 27 K<T_S<630 K (2.3 K<T_S<640 K) at 68% (95%) confidence. They therefore also place a lower bound on X-ray heating, a previously unconstrained aspects of early galaxies. For example, if the CMB dominates the z~8 radio background, the new HERA limits imply that the first galaxies produced X-rays more efficiently than local ones (with soft band X-ray luminosities per star formation rate constrained to L_X/SFR = { 10^40.2, 10^41.9 } erg/s/(M_sun/yr) at 68% confidence), consistent with expectations of X-ray binaries in low-metallicity environments. The z~10 limits require even earlier heating if dark-matter interactions (e.g., through millicharges) cool down the hydrogen gas. Using a model in which an extra radio background is produced by galaxies, we rule out (at 95% confidence) the combination of high radio and low X-ray luminosities of L_{r,\nu}/SFR>3.9 x 10^24 W/Hz/(M_sun/yr) and L_X/SFR<10^40 erg/s/(M_sun/yr). The new HERA upper limits neither support nor disfavor a cosmological interpretation of the recent EDGES detection. The analysis framework described here provides a foundation for the interpretation of future HERA results.


INTRODUCTION
One of the final frontiers of observational cosmology is the Cosmic Dawn, during which the first luminous sources formed and grew into galaxies. This era ended with the reionization of the intergalactic medium (IGM), when ultraviolet photons from these sources ionized virtually all the neutral hydrogen -and hence when stars and black holes affected every baryon in the Universe. This constitutes the last baryonic phase transition in the Universe's history and has important implications for later generations of galaxies.
Observations are now beginning to probe this era. Measurements of the large-scale polarization of the cosmic microwave background (CMB) imply that reionization reached its midpoint at z ∼ 7-8 (Planck Collaboration 2018;de Belsunce et al. 2021;Heinrich & Hu 2021). Models of Lyman-α emission lines of galaxies (Stark et al. 2010;Schenker et al. 2012;Jensen et al. 2013a;Caruana et al. 2014;Pentericci et al. 2014;Mesinger et al. 2015;Mason et al. 2018Mason et al. , 2019 and quasars (Mesinger & Haiman 2004;Bolton et al. 2011;Greig et al. 2017;Davies et al. 2018;Yang et al. 2020;Wang et al. 2020;Greig et al. 2019) also suggest a relatively large neutral fraction at z ∼ 7. While the conventional wisdom has long held that the reionization process ends at z ∼ 6 (e.g. McGreer et al. 2015a; though see Lidz et al. 2006;Mesinger 2010), recent measurements of the Lyman-α forest suggest that it may continue to somewhat later times (Becker et al. 2015; Choudhury et al. 2021b;Qin et al. 2021a). However, our understanding of this era is still incomplete: models and empirical extrapolations suggest that even the deepest Hubble Space Telescope (and upcoming JWST) observations probe only a fraction of the total star formation in the early Universe (Robertson et al. 2015;Behroozi & Silk 2015;Mason et al. 2015;Furlanetto et al. 2017;Gillet et al. 2020). This could mean that the galaxies providing most of the reionizing photons will remain unseen. Moreover, while reionization is the most dramatic effect of the first galaxies, their Xray and ultraviolet radiation fields can affect the IGM even while it remains neutral -a phase that cannot be observed directly by many cosmological probes.
A complete understanding of the Cosmic Dawn therefore requires complementary measurements of the IGM gas. The most powerful potential probe is the 21-cm spin-flip line of neutral hydrogen (Field 1959;Madau et al. 1997). The 21-cm line is particularly sensitive to (Furlanetto 2006;Morales et al. 2012;Pritchard & Loeb 2012): (1) structure formation in the Universe, which can be observed through density fluctuations; (2) the reionization process, which eliminates the 21-cm signal inside the large ionized bubbles that grow throughout that era; (3) the X-ray background (or other exotic heating or cooling mechanisms), which likely sets the IGM temperature before reionization and hence determines whether the 21-cm line is seen in absorption or emission; (4) the non-ionizing ultraviolet background, as photons that redshift into the hydrogen Lyman-α transition mix the hyperfine level populations; and (5) the radio background at high redshifts, including the CMB but also potential contributions from astrophysical sources or exotic processes.
Because the spin-flip cosmological signal is very weak compared to other astrophysical radio backgrounds, mapping these IGM fluctuations is extremely challenging, and early efforts to observe it have focused on two complementary directions. One is the "global" all-sky signal, measuring the sky averaged spectral signature of the line, covering Gpc 3 -sized comoving volumes (Shaver et al. 1999;Muñoz & Cyr-Racine 2021). Several such experiments are underway (Price et al. 2017;Singh et al. 2017;Philip et al. 2019;Voytek et al. 2014;DiLullo et al. 2020). Of these, only the EDGES collaboration has made a tentative detection (Bowman et al. 2018), although the cosmological interpretation of the measurement is subject to significant instrumental and systematic uncertainties (e.g. Hills et al. 2018;Sims & Pober 2019;Bradley et al. 2019;Singh & Subrahmanyan 2019;Tauscher et al. 2020). Interestingly, the claimed signal is much stronger than expected, requiring either that the IGM temperature is smaller than allowed by adiabatic cooling (from, e.g., energy exchange with dark matter; Barkana 2018;Berlin et al. 2018;Slatyer & Wu 2018;Kovetz et al. 2018), or that an additional radio background (beyond the CMB) is present in the early Universe (e.g., Feng & Holder 2018;Pospelov et al. 2018;Ewall-Wice et al. 2018a;Fialkov & Barkana 2019;Mebane et al. 2020).
A number of other experiments hope to use interferometers to measure statistical fluctuations in the 21-cm background, most often quantified through the power spectrum, which measures the variance in the field as a function of smoothing scale. Several experiments have now published upper limits from z ∼ 6-10, though these limits so far probe only a small fraction of the parameter space spanned by "standard" models of early galaxies (Mondal et al. 2020;Ghara et al. 2020;Greig et al. 2021b,a;Ghara et al. 2021).
Recently in HERA Collaboration (2021) (hereafter H21), we presented the first upper limits on the 21-cm power spectrum from a new experiment, the Hydrogen Epoch of Reionization Array (HERA). HERA is now under construction in the Karoo Desert of South Africa (DeBoer et al. 2017). Its phased construction allowed an initial observing campaign in 2017-18, whose results are considered here. Note that we base these results on data from just 39 antennas; HERA is now expanding to ∼ 350 antennas, so the interpretation here provides a framework for improved analyses in the future. This paper is organized as follows. We introduce the physics of the 21-cm signal in section 2. Then, in section 3, we describe HERA's limits and our inference tools. In section 4, we use a very simple model to motivate the most important implications of HERA's upper limit. In the following four sections, we present several complementary interpretations to elucidate these results: we use the 21cmMC code to infer constraints on early-galaxy populations and the IGM (section 5) and a phenomenological model that directly parametrizes IGM properties to better understand the IGM constraints (section 6). Then, we examine the implications of the HERA limits for exotic dark-matter models (section 7), and finally we consider constraints derived from models with an enhanced radio background (section 8). In section 9, we summarize these results and their implications for the epoch of reionization.
Throughout this work, we assume a standard flat ΛCDM cosmology, consistent with the latest CMB measurements (Planck Collaboration 2018). The separate analyses use slightly different cosmological parameters, but these have little effect on our constraints. We denote comoving Mpc with cMpc.

THE 21-CM SIGNAL
HERA and other low-frequency instruments aim to observe emission or absorption of the neutral hydrogen hyperfine transition at an observed wavelength of λ obs = 21(1 + z) cm. The intensity of this line is conventionally expressed as the differential brightness temperature, δT 21 , relative to the low-frequency radio background, which we assume has a brightness temperature T rad at the relevant frequency. Then the brightness of a patch of the IGM can be expressed approximately as (Madau et al. 1997;Furlanetto 2006) δT 21 (ν) = T S − T rad 1 + z 1 − e −τν 0 where T 0 = 27[(1 + z)/10] 1/2 mK is the overall normalization, H is the Hubble parameter at the appropriate redshift, and we have assumed Ω m h 2 = 0.15 and Ω b h 2 = 0.023 (with H 0 = 100h km/s/Mpc). Here, x HI is the neutral fraction of the patch, δ = (ρ −ρ)/ρ is its fractional overdensity, T S is the spin temperature (or the excitation temperature of the 21-cm transition), and dv r /dr is the gradient of the proper velocity along the line of sight.
The spin temperature is determined by (e.g., Madau et al. 1997;Furlanetto 2006;Pritchard & Loeb 2012;Venumadhav et al. 2018a) where x rad , x α , and x c are coupling constants describing the strength of the relevant interactions. This equation 4 The HERA Collaboration reflects the competition between several processes: (1) interactions with radio photons tend to drive T S to T rad , with a coupling constant x rad ; (2) collisions drive T S toward the kinetic temperature of the gas, T K with a coupling constant x c ; and (3) absorption and re-emission of Lyman-α photons mixes the hyperfine states and also drives T S toward T K with a coupling x α , through a process known as the Wouthuysen-Field effect (Wouthuysen 1952;Field 1958Field , 1959Hirata 2006). Meanwhile, the kinetic temperature is affected by the expansion cooling of the IGM and interactions with several radiation backgrounds -most importantly, prior to reionization, any X-ray background generated by early sources. A proper accounting of the temperature requires tracking both the IGM properties and the radiation backgrounds generated by galaxy formation or exotic processes in the early Universe.
It is important to note that, in the standard picture, reionization by UV photons is an inhomogeneous process -(nearly) completely ionized regions around the first galaxies expand into (nearly) completely neutral IGM patches as the source population grows. The values of x Hi we quote below can therefore be considered as approximately corresponding to the volume filling factor of the remaining neutral IGM patches during the EoR. However, X-rays have much longer mean free paths than UV photons and can deposit their energy in the neutral IGM, partially ionizing and heating that phase, so the relation between the true neutral fraction and the filling factor of the ionized bubbles is not exact.
Given the sensitivity of current experiments, the focus of interferometric observations to date has been on measuring the spatial power spectrum of the 21-cm signal, δ T 21 (k 1 )δT 21 (k 2 ) = (2π) 3 δ D (k 1 + k 2 )P 21 (k 1 ), (3) where tildes denote Fourier transforms, angular brackets denote ensemble averages, and δ D is the Dirac delta function.
The velocity term in equation (1) accounts for the mapping between redshift and real space, which is complicated by redshift-space distortions (RSDs ;Kaiser 1987;Bharadwaj & Ali 2004;Barkana & Loeb 2005a). Crudely, overdense regions expand more slowly than the average Universe, so they appear compressed along the radial direction, while underdense regions appear larger in that direction. Because these distortions occur only along the line of sight, they make the power spectrum anisotropic. The modes used in the HERA analysis are mostly aligned along the line of sight and care must be Note that in the present analysis we only include every other k bin from the limits quoted in H21 in order to mitigate the effect of non-zero covariance between neighboring k modes. We start this decimation at k = 0.179/cMpc and k = 0.134/cMpc for Bands 1 and 2, respectively. Filled points represent positive measurements, and errorbars without points represent negative measurements. The errorbars show ±1σ uncertainty.
taken when comparing to the theoretical models, as we discuss further below.

HERA PHASE I POWER SPECTRUM LIMITS
Next we establish the formalism that will be used to interpret our observables. Section 3.1 describes HERA's data products that are used in this paper, and section 3.2 defines the likelihood that links these observable quantities to theoretical models. The goal is therefore to provide the necessary machinery to interpret our measurements in a model-agnostic way before introducing our theoretical models in subsequent sections.

Observational Campaign
The power spectrum upper limits analyzed in this paper have been published in H21. Here, we describe some of the essential features of the data for convenience, but we refer the reader to H21 for more details.
The upper limits relevant to the paper are reproduced in Figure 1. These were based on 18 nights of data (Julian Dates 2458098 to 2458116) taken as part of an observing campaign from October 2017 to April 2018 when HERA was in its Phase I observing configuration. In Phase I, HERA observed with "hybrid" antenna elements which consisted of HERA's 14-m parabolic antennae with modified cross-dipole feeds and a front-end from the PAPER experiment (Parsons et al. 2010;De-Boer et al. 2017). HERA Phase I also inherited PA-PER's back-end system, which processed 100 MHz of bandwidth from 100-200 MHz. For these observations, HERA consisted of 52 operating antennas, 39 of which were deemed science-ready after passing our data quality metrics (H21). Note that these 52 antennas make up a small fraction of the experiment at full capacity of ∼350 antennas, which will observe from 50 -225 MHz (Dillon & Parsons 2016;DeBoer et al. 2017).
The analysis and reduction of these data are discussed in H21 and in several supporting papers in more detail (Kern et al. 2020a,b;Dillon et al. 2020;Tan et al. 2021;Aguirre et al. 2021). For the purposes of this work, the important takeaway is that, while nearly the full band is processed in the data reduction pipeline, only two portions of the band are largely free of radio frequency interference (RFI), which sets the redshift ranges studied in this work (Band 2, centered at z = 7.9, and Band 1, centered at z = 10.4). Additionally, the power spectra studied in this work come from only one of the fields reported in H21. Because HERA observes in a drift scan mode, it surveys a ∼ 10 • -wide stripe centered on declination −30.7 • . However, to avoid the brightest portions of the sky (including foregrounds from our Galaxy as well as bright sources such as Fornax A), H21 made further cuts to the data in local sidereal time (LST). This yields three fields (with LST ranges from 1.25 to 2.7 hours, 4.5 to 6.5 hours, and 8.5 to 10.75 hours) worth of data that were propagated through to the power spectrum pipeline. The parameter inference discussed in this work comes solely from the limits presented from the first cut (Field 1; see Fig. 1 in H21), as these showed the least amount of foreground contamination and therefore produced the most stringent limits.
For the z ∼ 8 band, the data presented in H21 provide the most sensitive upper limits on the 21 cm power spectrum to date, improving upon previous limits at that redshift by roughly one order of magnitude. Another important feature of the H21 analysis is that they report measurements consistent with the thermal noise floor at intermediate and high Fourier k wavevectors. The dynamic range between that noise floor and the peak measured foreground signal is ∼ 10 9 in power, in spite of the fact that they perform no explicit foreground subtraction in their analysis. Upper limits on the 21 cm power spectrum are also currently best constrained by the MWA at lower redshifts (Trott et al. 2020) and by LOFAR at higher redshifts (Gehlot et al. 2019;Mertens et al. 2020).

Data Likelihood
To relate our power spectrum measurements to theoretical models, we first group our data at all k bins and redshifts into a column vector, i.e., which has a length of N d = N k × N z . In this work, we use the power spectrum data tabulated in H21 Tables  3 and 4 for Field 1 only, spanning a k range of 0.13-0.64 cMpc −1 and the two redshift bins z = 10.4 and z = 7.9. Furthermore, we also make use of the associated window function and covariance matrices, which are included with the data and will be publicly accessible. In this work, we assume the thermal noise on the data to be Gaussian distributed and thus adopt a Gaussian likelihood. This is a fair approximation as the large amounts of averaging performed in the analysis Gaussianizes the data due to the central limit theorem. Having adopted a model M for the cosmic 21 cm signal (e.g. one of the simulations described in later sections), m, and a model for any extant systematics, u, we can write the probability distribution for the data given the parameters (i.e., the likelihood function ) as where r(θ, u) = d − u − Wm(θ), θ are the parameters of M, m is the simulation's deterministic prediction of the data vector mean given θ, W is the N d × N d window function matrix of the data 1 , and Γ = Σ −1 is the N d × N d precision matrix, which is the inverse of the covariance matrix of the data. The window functions account for the corrections to the predicted mean vector due to the telescope measurement and data reduction process (c.f. Tegmark 1997;Liu & Tegmark 2011;Dillon et al. 2013;Liu & Shaw 2020;; in other 6 The HERA Collaboration words, it is the point spread function of the power spectrum measurement in Fourier k space. The covariance matrix accounts for the variance of the measured power spectrum and the correlation of that uncertainty between band powers, irrespective of non-thermal systematics. This covariance is assumed to be diagonal given the analysis methods in H21 (see Section 3.2.1 for details). The on-diagonal elements (i.e., the variances) are estimated using antenna auto-correlation data to model the instrument noise. Because the power spectrum is a quadratic statistic, the sky signal enters in various signal-noise cross terms even if our variance model is due entirely to instrumental noise. For this contribution it is the total sky signal (including foregrounds) that matters, and we model this using the empirically measured power spectrum, as detailed in Tan et al. (2021). Note that while we write Σ as model independent, there are some terms that can be model dependent, and thus it can take on an explicit dependence on θ (cosmic variance, for example, is dependent on the amplitude of the predicted mean signal). For the current limits, we do not expect cosmic variance to be important (H21). Ultimately, one is interested in the probability distribution of the parameters θ given the data d, i.e. the posterior probability distribution p(θ, u|d, M). This is related to the likelihood via Bayes' theorem: where p(θ|M) is our prior distribution on the parameters and p(u) is our prior on the systematics (assumed to be independent of the physical parameter prior).

Marginalizing over systematics
The likelihood as expressed in equation (5) (and thus the posterior in eq. 6) has a dependence on the systematics, u. In this paper, we have no explicit way of modelling u, so we desire a likelihood which is dependent only on the astrophysical parameters. This suggests marginalizing over the prior range of the unknown systematics. In principle, we would express u = u(φ), i.e. we would have some physically motivated set of parameters φ that produce a set of systematics in ∆ 2 21 (k, z), and we would marginalize the posterior over these parameters. In the absence of such a physically motivated model, we marginalize directly over the binned values u: In particular, taking a multivariate uniform prior on u gives Note the assumptions that have been made in obtaining this expression. Here, our multivariate uniform prior on u allows each k and z bin to vary independently, thus allowing random fluctuations of arbitrary form. Although this is the form we employ for this paper, future analyses would be considerably improved with detailed physical models for systematics that might be present, for example by imposing smoothness priors (in k and/or z) when appropriate.
If we also assume that Γ is diagonal, and writing t = d − Wm(θ), then this equation reduces to where the second line follows due to the separability of the factors in u i , and σ i = (Γ −1/2 ) ii is the standard deviation of d i . In this paper, we utilize bandpowers that are widely separated in wavenumber (see below) so that a diagonal covariance matrix is a good approximation.
In order to provide a systematics-marginalized likelihood, we must choose prior ranges for the systematics on each (k, z)-bin. Allowing for unbounded (possibly negative) systematics would not allow us to constrain the cosmological signal as the systematic would be completely degenerate with the model. Thus, we should look to our understanding of the data analysis process to set this prior. Calibration errors causing residual phase differences or chromatic effects can lead to biases that are positive or negative. Though negative biases or systematics in the power spectrum have been observed in previous experiments (see, e.g., Kolopanis et al. 2019), for the present H21 dataset and analysis pipeline the most likely causes of such issues have been mitigated by the application of absolute calibration and other improvements; large negative detections are not observed in null tests or validation simulations. The most likely remaining source of systematic bias is unmodeled signal chain chromaticity common to all elements, and this would couple positive foreground power beyond the wedge. Given this expectation of positive-only systematics, we set the prior constraint that u ≥ 0, yielding where erf is the error function. It is worth making clear that this form of the posterior is relatively flat once t i σ i . Since t i represents the data minus the theory model, in effect this means that our posterior produces close to equal probability for any scenario in which the model is less than or equal to measured values (within error bars). Our treatment of systematics therefore leads to a well-defined posterior that naturally treats data points as "upper limits". This result is the same form as that derived in Appendix B of Ghara et al. (2020) in the interpretation of LOFAR data. A similar derivation of the marginal upper-limit likelihood can also be found in Appendix A of Li et al. (2019). If the off-diagonal components of Γ are not zero, the integral of Equation 7 is not tractable in closed form. For the HERA data used in this work, we specifically use band powers that are widely separated in k such that their error correlations are negligibly small; concretely we use only every second k-bin ("decimation"). Quantitatively, Figure 20 of HERA Collaboration (2021) shows an example of the normalized covariance between k bins, demonstrating that after decimation the remaining modes have negligible covariance, on the order of 1-2%. In decimating, we could in principle choose either the even or odd k-bins from each band, and as each of the four choices would be a slight underestimation of the constraints, we choose the combination providing the strongest limits. This includes k = 0.17 cMpc −1 in Band 1 and k = 0.13 cMpc −1 in Band 2; as a matter of convention we refer to the former as even and the latter as odd k-bins. As we will show later (see Fig. 3 and 19), the constraints on realistic models are primarily driven by the two most stringent limits, so we can expect the decimation to have a negligible effect.

"Inverse" Likelihood
In practice, given that the upper limits presented in H21 are still roughly two orders of magnitude above fiducial 21 cm models, the majority of the parameter space for standard models is left unconstrained. One way to illustrate how the new limits help is by combining them with existing constraints. Alternatively, to provide a clearer picture of the model parameter choices that exceed the HERA limits we also consider an "inverse likelihood" defined as where L 0 is the maximum of L m . 2 With the inverse likelihood, the resulting marginalized distributions identify the parameter combinations that can be ruled out by the HERA limits alone (see Fig. 4 for an example). However, these distributions must be treated with caution: models that lie inside of the projections of the full 2 This typically is the likelihood of a model power spectrum equal to zero. Including L 0 here makes sure L inv is independent of the normalization of Lm, e.g. of the number of data points used. . Lower limits on T S from HERA data (95% C.L., green and purple arrows) compared to the upper limit from EDGES (blue arrow). The green HERA limits have been obtained by assuming that the IGM is fully neutral and at a constant temperature (aside from small fluctuations due to adiabatic expansion). The purple HERA limit is from full galaxy models at z = 7.9 with 21cmMC (see section 5); note that by construction the 21cmMC models cannot cool below the adiabatic prediction, so we do not show a 21cmMC limit at z = 10.4. Both assume spherically symmetric RSDs. The black line shows the standard-model prediction assuming full Wouthuysen-Field coupling (TS = TK ) without any heating. The HERA Band 2 data can rule out adiabatic cooling in both approaches, requiring some heating to take place before z = 7.9.
distribution are not necessarily excluded. The inverse likelihood should only be used to gain intuition about the utility of the HERA limits and parameters that are necessary (but not sufficient) to drive a power spectrum beyond the HERA limits.

BUILDING PHYSICAL INTUITION: A DENSITY-DRIVEN BIAS APPROACH
Before studying galaxy-driven models, we begin with a simple bias analysis, which will allow us to build intuition about the implications of the HERA measurements. As these limits are well above predictions of "vanilla" models of the reionization era, the most important parameter that can be constrained is the IGM temperature, as the temperature ratio term in equation (1) can become arbitrarily large for gas that is very cold. In the spirit of simplicity, throughout this section we will assume efficient Wouthuysen-Field coupling (T K = T S ), sourced by a non-ionizing ultraviolet background from early star formation. We then infer constraints from the most stringent HERA k-bin at each redshift, which we will interpret in terms of changes to the gas kinetic temperature (i.e., we will set T rad = T CMB ).

8
The HERA Collaboration Let us begin by supposing that the 21-cm power spectrum traces the matter power spectrum, which is appropriate when the fluctuations are sourced by the matter fluctuations. The key assumption here is that the bias parameter is scale-independent -which is exact for ionization and temperature that vary linearly with density, and can be extended beyond this approximation (McQuinn et al. 2005). We then use the HERA measurements to constrain the bias parameter b m . We compute the linear matter power spectrum from CAMB 3 (Lewis & Bridle 2002) at each z, and we find the 95% CL limits of b m < {156, 529} mK for z = {7.9, 10.4} (an analogous analysis for the relativevelocity power spectrum can be found in Appendix A). These can be translated into lower limits on the ratio of the spin-to-radio temperatures through the relation derived from equation (1) (see e.g., Pritchard & Loeb 2008), where C T is the adiabatic index, which accounts for the preferential cooling of underdense regions; and µ is the line-of-sight cosine of the wavenumbers observed, which accounts for RSDs. Here, and throughout this text, an overbar represents average over volume. We obtain the adiabatic index C T = δT K /(T K δ) ≈ 0.6 as a function of z following Muñoz et al. (2015a), which we correct for kinetic temperatures above the adiabatic threshold (T K > T ad K ) by writing C T (T K ) → C T × min(1, T ad K /T K ), assuming homogeneous heating. We further assume negligible ionizations, setting x HI = 1, and for RSDs we take spherically averaged modes (µ = 0.6) through this section to match the common procedure done in simulations. We will show how the constraints shift when altering these two assumptions later in section 7.
Under our assumptions, the b m upper limits translate into lower limits for the spin temperature of at 95% confidence, where we re-emphasize we have assumed x HI = 1. We show these limits in Fig. 2, along with the adiabatic-cooling prediction in the standard CDM model. These T S values have interesting implications for the thermal state of the IGM at high redshifts. As is clear from Figure 2, the HERA Band 2 (z = 7.9) 95% confidence limit is above the adiabatic-cooling prediction, which demands that some heating must have occurred before z = 7.9. Moreover, the HERA limits for Band 1 (z = 10.4), while below the adiabatic limit at that z, can be used to clarify the state of the IGM in comparison with the claimed EDGES detection (also shown in Fig. 2), which we will explore in section 7. We emphasize that these limits rest on three strong assumptions, which we will highlight here and will, in the following sections, explore with more physics-rich models. First, the limits assume full Wouthuysen-Field coupling (T S = T K ), which is all but guaranteed by the redshifts we consider (z ≤ 10.4, see section 5). Second, they assume a value of x HI = 1, which can be varied at each z in equation (12), though only homogeneously. Lastly, in this analysis we have performed a spherical average of RSDs, whereas HERA data mostly contains modes along the line of sight (µ ≈ 1, see section 2). Properly accounting for RSDs can result in stronger limits, as we will show in section 7.
In summary, the bias approach here outlined is useful for building intuition, although reionization models and observations suggest that the spatial fluctuations in the ionization field (rather than the matter field) should drive the 21-cm signal at z ∼ 8. We will explore such models in detail in the following sections, but for now we show the limit from 21cmMC in Fig. 2. We describe below how this limit was obtained, but we see already that there is general agreement (∼ factor of few) with the density-driven bias limit at z ≈ 8 (indeed our density-driven bias limit is very close to the analogous density-driven 21cmMC limit, denoted by the red contours in Figure 5). The reason these two approaches yield similar results -despite their vastly different assumptions about the EoR -is that the density and ionization power spectra are of the same magnitude at z ∼ 8 and k ∼ 0.1 Mpc −1 (e.g. Furlanetto et al. 2004a). Astrophysical models can only modify the peak power during the EoR by a factor of a few (e.g. Greig & Mesinger 2015). The only way to reach the power spectrum amplitudes probed by the HERA limits is by having a large ∼ (1 − T rad /T S ) 2 pre-factor, i.e., requiring T S T rad . In this regime, model differences can be easily compensated by relatively small changes in T S . Therefore, in the regime of current HERA limits, constraints on T S are of the same magnitude whether the 21-cm power spectrum tracks density or ionization fluctuations.

GALAXY AND IGM PROPERTIES INFERRED FROM HERA OBSERVATIONS
We next consider the HERA limits in light of "standard" galaxy formation models using data-constrained 21cmFAST semi-numerical simulations.

Galaxy-driven models of the cosmic 21-cm signal
Here we briefly summarize how the 21-cm signal is computed using the galaxy-driven models of 21cmFAST 4 (Mesinger & Furlanetto 2007;Mesinger et al. 2011;Murray et al. 2020). The main ansatz of these 21-cm models is that cosmic radiation fields are sourced by galaxies, hosted by dark matter halos (whose relation to the large-scale matter field is comparably well-understood). We generate Eulerian density and velocity fields with second-order Lagrangian perturbation theory (2LPT; e.g. Scoccimarro 1998). Galaxy properties are then assigned to dark-matter halos via scaling relations with halo mass. In section 6, we explore toy models in which radiation fields are not directly associated with galaxies in order to study the robustness of our inferences and the "value-added" by explicit models of structure formation.
Specifically, we use the empirical galaxy relations of Park et al. (2019), capable of reproducing the observed UV luminosity functions of galaxies during the EoR (z = 6 -10), as well as the spatial distribution of IGM opacities seen in Lyα forest spectra at z = 5 -6 (Qin et al. 2021a). Consistent with semi-analytic models and hydrodynamic simulations of high-z galaxies (e.g. Moster et al. 2013;Xu et al. 2016;Sun & Furlanetto 2016;Mutch et al. 2016;Tacchella et al. 2018;Behroozi et al. 2019;Yung et al. 2019;Ma et al. 2020), we describe the mean stellar to halo mass relation, M * /M h , with a power law: where (Ω b /Ω m ) is the mean baryon fraction, and the stellar fraction, f * = f * ,10 (M h /10 10 M ) α * , is restricted to be between 0 and 1. The corresponding starformation rate assumes a characteristic star-formation time-scale that scales with the Hubble time, H −1 (which during matter domination is equivalent to scaling with the halo free-fall time): We then compute the galactic emissivities (soft UV, ionizing UV and X-ray), assuming they scale with the star formation rates. We identify ionized regions with an excursion set approach (Furlanetto et al. 2004b), comparing the cumulative (local) numbers of emitted photons and recombinations. We slightly adjust the number of emitted photons to correct for the non-conservation of ionizing photons in excursion set algorithms (e.g. Zahn et al. 2007;Paranjape & Choudhury 2014a; for details see Park et al. in prep). Sub-grid IGM recombinations are tracked according to Sobacchi & Mesinger (2014). We assume PopII stellar SEDs for the ionizing and soft UV emission, corresponding to ∼5000 ionizing photons produced per stellar baryon (e.g. Leitherer et al. 1999;Barkana & Loeb 2005b) 5 . A fraction 1 − f esc of these photons is absorbed within the galaxy itself, and does not reach the IGM. We allow the ionizing escape fraction to also scale with the halo mass: where f esc,10 is the normalization and α esc is a powerlaw index. The ionizing escape fraction is also restricted to values between 0 and 1. Although there is currently no consensus on the ionizing escape fraction or its dependence on galaxy properties, simulations suggest that such a generic power law is an acceptable characterization of the population-averaged values (e.g. Paardekooper et al. 2015;Kimm et al. 2017;Lewis et al. 2020).
In contrast to ionizing UV photons, the soft UV and X-ray photons responsible for coupling the gas and spin temperatures and heating the gas, can have long mean free paths through even the neutral IGM. We follow the corresponding ionization and heating rates for each simulation cell, by integrating the specific emissivities back along the lightcone, attenuated by the corresponding opacities. Our simulations track the spatial fluctuations in the X-ray and Lyman series backgrounds, with the IGM opacity computed assuming a standard "picketfence" absorption for Lyman series photons and absorption from partially-ionized hydrogen and helium in a two-phased IGM for X-ray photons (e.g. Mesinger et al. 2011Mesinger et al. , 2013Qin et al. 2020a). The X-ray SED emerging from galaxies is approximated as a power-law whose luminosity scales with the SFR. This is consistent with theoretical models and observations of local star form- The HERA Collaboration ing galaxies, for which X-ray emission is dominated by high-mass X-ray binaries (HMXBs) and/or the hot ISM (e.g. Fragos et al. 2013;Mineo et al. 2012;Pacucci et al. 2014;Brorby et al. 2014;Lehmer et al. 2016). Specifically, we parametrize the typical emerging X-ray SED of high-z galaxies via their integrated soft-band (< 2 keV) luminosity per SFR (in units of erg s −1 M −1 yr), where L X /SFR is the specific X-ray luminosity per unit star formation escaping the host galaxies in units of erg s −1 keV −1 M −1 yr, taken here to be a power law with energy index α X and E 0 is the minimum energy for Xrays to be able to emerge from the galaxy and not be absorbed locally in the ISM. For reference, the typical value of E 0 ∼ 0.5 keV found in the simulations of Das et al. (2017) corresponds to an HI column density of ∼ 10 21.4 cm −2 , assuming zero metallicity. In summary, our 21cmFAST galaxy models have nine free parameters: 1. f * ,10 , the normalization of the stellar mass-halo mass relation, evaluated at M h = 10 10 M 2. α , the power law index of the stellar mass-halo mass relation 3. f esc,10 , the normalization of the ionizing escape fraction-halo mass relation, evaluated at M h = 10 10 M 4. α esc , the power law index of the ionizing escape fraction -halo mass relation 5. M turn , the characteristic halo mass scale below which the abundance of active galaxies is exponentially suppressed 6. t * , the characteristic star formation time scale, expressed in units of the Hubble time 7. L X<2 keV /SFR, the soft-band X-ray luminosity per unit SFR 8. E 0 , the minimum X-ray energy of photons capable of escaping their host galaxies 9. α X , the energy power law index of the X-ray SED We emphasize that this flexible galaxy parametrization used in 21cmFAST enables us to set physically meaningful priors over the free parameters and use high-z galaxy observations in our inference. For instance, the common simplification of a constant stellar to halo mass relation is inconsistent with galaxy SFR and LF observations, and can thus bias parameter inference (c.f. Mirocha et al. 2017; Fig. 1 in Park et al. 2019). Our galaxy model therefore allows us to use existing highz observations, in addition to HERA, when computing the model likelihood (see section 5.4). This quantifies the "added value" of HERA, given that existing observations already exclude a significant prior volume (e.g. Park et al. 2019). Without them our posterior would strongly depend on our priors.

Inference
To perform Bayesian inference, we use 21cmMC 6 (Greig & Mesinger 2015, 2017a with the recently implemented Multinest-based (Feroz et al. 2009) sampler (Qin et al. 2021b; see also Binnie & Pritchard 2019). For a given sample of astrophysical parameters, we compute 4D realizations of the 21-cm signal in a cubic volume with a periodic boundary condition and a length of 250 cMpc. The initial conditions and 2LPT are calculated on a 512 3 grid, while the final radiation fields are computed on a 128 3 grid. Choosing a line-of-sight axis, we account for nonlinear RSDs via the real-to-redshift space sub-grid transformation described in Greig & Mesinger (2018), and first introduced in Mao et al. (2012); Jensen et al. (2013b). 7 When evaluating the likelihood according to equation (5), we add in quadrature a conservative 20% modeling error (c.f. Zahn et al. 2011) as well as the sample variance from our simulation. In contrast to other simulation-based inference codes, 21cmMC forward models 4D realizations of the 21cm signal. We compute the power spectra (PS) on-the-fly from these 4D realizations without emulators, over our 9-dimensional parameter space; we therefore do not include emulator error/bias in our likelihood.
When we include other observational constraints in our inference procedure (see Section 5.4), we calculate the total likelihood with L total = L m ×L LFs ×L DF ×L τe , where the last three terms reflect the comparison be- with RSD without RSD Figure 3. Two examples of galaxy-driven models that are ruled out by HERA 2021 limits. The rows correspond to Band 1 (z = 10.4; top) and Band 2 (z = 7.9; bottom). The columns correspond to (from left to right): (i) the IGM density; (ii) the brightness temperature of a "reionization driven" model (xHi=0.73, T S=1.42 at z=7.9; black patches are cosmic HII regions); (iii) brightness temperature of a "density driven" model (xHi=0.98, T S=1.99 at z=7.9); and (iv) the corresponding power spectra together with the H21 limits. Note that the power spectra of the models are much flatter (with k) than the observational limits; thus the constraining power is entirely provided by the two filled squares at low k. In the bottom right panel, we also show power spectra ignoring redshift space distortions; RSDs are important for the density driven models but much less so for the reionization driven models. The slices are 1 cGpc on a side and 2 cMpc thick and were generated with 21cmFAST v3.  (Qin et al. 2021a), seem to favor a slightly later end to reionization (a delay of ∆z ∼ 0.5).
Since reionization-driven PS amplitudes are maximized around the midpoint of the EoR, which for current observations occurs right around HERA's Band 2 at z ∼ 8, we expect that shifting the EoR towards later times could slightly weaken the HERA constraints we derive below.
In this section we highlight the astrophysical models disfavored by the current HERA limits. To do this, we use the inverse likelihood from equation (10). Because the inverse likelihood is only illustrative, we also confine the analysis to the two most stringent limits at z = 7.9 and z = 10.3. In any case, these two data points provide all the constraining power because the observed limits rise much more steeply with k than the model predictions. This allows us to compare to similar analysis of recent LOFAR and MWA data, which also used an inverse likelihood and the same galaxy models (Greig et al. 2021a,b).
Before showing the full distribution of models, in Fig.  3 we show examples of two classes of models capable of exceeding the HERA upper limits. The top row corresponds to Band 1 (z = 10.4) and the bottom to Band 2 (z = 7.9). Slices through the density field and 21-cm brightness temperature fields are shown on the left, with the 21-cm power spectra shown together with the data in the rightmost panels. For visualization purposes, the maps are generated from larger boxes than used in the 12 The HERA Collaboration . Distribution of models disfavored by H21, calculated using the inverse likelihood (eq. 10) and using only the two most constraining data points in Bands 1 and 2 (c.f. filled squares in Fig. 3). The 1D and 2D marginalized distributions were generated by assuming flat priors over the ranges shown by the figure; we caution that these marginalized inverse likelihood results should not be interpreted strictly as a "posterior", but instead serve to illustrate where the models disfavored by H21 sit in astrophysical parameter space. In the bottom left we show the 2D and 1D distributions, while in the top right we show the EoR history, global signal evolution and power spectrum evolution at k = 0.13 cMpc −1 . Red / blue curves denote "density driven" / "reionization driven" models, classified according to the value of the neutral fraction at z = 6; shaded regions enclose 68% of the disfavored models for each mode. In the power spectrum evolution plot we also show the two H21 data points used to compute these distributions (note that the Band 1 data point is at a slightly higher wavenumber of k = 0.17 cMpc −1 ); this highlights that the Band 2 (z = 8) data point has all of the constraining power. In the EoR history panel, we also include the QSO dark fraction upper limits from McGreer et al. (2015b) (empty square). In the bottom right, we also include the corresponding PDFs of the CMB optical depth, τe, from both modes; the gray region spans 68% C.L. of the observed value, implied by the galaxy-model recovery of Planck Collaboration (2018) EE power spectra described in Qin et al. (2020b). These two EoR observations were not used in the inverted likelihood; unlike the "density driven" modes, the "reionization driven" modes are largely consistent with these limits.
inference, corresponding to 1 cGpc on a side, but with the same 2 cMpc resolution. 9 The 21-cm power spectra in the two classes of models exceeding these upper limits are driven by spatial fluctuations in either: (i) the IGM ionized fraction, which we will refer to as "reionization driven" (also referred to as "cold reionization" in the literature; e.g. Mesinger et al. 2014) ; or (ii) the gas density, which we will refer to as "density driven" (see section 4, and Greig et al. 2021a for the same qualitative result using recent MWA limits). 10 Both scenarios require a cold IGM, which sets a lower limit on the heating rate (and hence on the X-ray emissivity within these models).
In the right panels, we confirm that the 21-cm power spectra of both scenarios are flatter than the observational limits. Thus when the observational limits are consistent with thermal noise, the constraining power comes entirely from the deepest limits (c.f. Mertens et al. 2020;Trott et al. 2020;Ghara et al. 2021) (in our case, primarily the deepest data point at z = 7.9).
In the bottom right panel, we also show how the power spectra depend on RSDs. For the "reionization driven" model, RSDs are not important since the first HII regions are highly biased, zeroing out the signal from the densest regions with the strongest RSDs Jensen et al. 2013a;Ghara et al. 2015;Ross et al. 2020). However, the "density driven" models have a negligible contribution from ionization and heating, with the 21-cm power spectrum driven entirely by the nonlinear matter field. By comparing the solid and dashed red curves, we see that non-linear RSDs can boost the spherically averaged power by factors of ∼ 2-3, in excess of the linear prediction of 1.87 (e.g. Bharadwaj & Ali 2004;Barkana & Loeb 2005a). Indeed without RSDs, this density driven model is consistent with the data at ∼ 1 σ. We explore the effect of different RSD assumptions for density-driven models in section 7.1.
In Figure 4 we show a corner plot corresponding to the inverted likelihood from equation (10). We caution that our parameter ranges in this figure/subsection 9 We confirm that the power spectra in the 1 cGpc and 250 cMpc runs are converged to the percent level or better for the relevant wave-numbers, k > 0.1 cMpc −1 . This level of convergence is consistent with the results of Kaur et al. (2020), who quantified the bias and scatter in the 21-cm signal resulting from missing large-scale modes (see also Iliev et al. 2014), and is orders of magnitude smaller than the observational uncertainties. 10 Ghara et al. (2020) also consider a model in which highly biased AGN with luminous, soft X-ray SEDs but negligible UV emission dominate the radiation background. Such extreme scenarios might also produce very strong temperature fluctuations, capable of exceeding the HERA limits; however, such models are not inside our prior volume.
do not correspond to a "prior" belief of the distribution of disfavored models, and marginalizing over an inverse likelihood is different from an inversion of the 2D marginalized Bayesian posteriors. Therefore Figure 4 should not be interpreted as a Bayesian posterior of disfavored models, and it is difficult to formally relate it to the normal likelihood results in the next section. However, the figure illustrates where the models that exceed HERA reside in our parameter space. In the top right, we draw from these distributions the redshift evolution of the mean neutral fraction, the mean 21-cm signal, and 21-cm power spectrum at k = 0.13 cMpc −1 .
Here we highlight the two modes discussed above: red and blue curves denote the "density driven" and "reionization driven" models, classified on the basis of whether the Universe is mostly neutral or mostly ionized at z = 6. 11 The shaded regions enclose 68% of the distributions. Astrophysically, the two modes are most easily distinguished by the ionizing escape fraction parameter, f esc,10 , and to a lesser degree by their star formation efficiencies, here parametrized by the ratio f * ,10 /t * . All of the models require that the IGM was not heated significantly, as seen by the upper limits on the X-ray luminosity per SFR, L X,<2kev /SFR.
The upper right panels show that the "density driven" models are already ruled out by other observations, since they fail to reionize the Universe early enough. In particular, we show the observed upper limit on the neutral fraction from the dark pixels in the Lyman forests (Mc-Greer et al. 2015b), as well as the Compton scattering optical depth from Planck 2018 (Qin et al. 2020b). Note that these observations were not used in computing the inverse likelihood. However, some of the "reionizationdriven" models are consistent with current observations. We return to this in the next subsection.
In Figure 5, we show where these HERA-disfavored models sit in the marginalized 2D space of x Hi vs T S . 12 11 There is a clear bimodality in the z ∼ 6 neutral fraction of disfavored models (c.f. top right panel of Fig. 4), allowing us to easily distinguish the "density driven" and "reionization driven" modes.
In this early EoR regime, the negative contribution of the ionization-density cross-correlation can result in a decrease of large-scale 21-cm power (e.g. Lidz et al. 2008;Zahn et al. 2011), making it difficult for those models to exceed the HERA limits. Thus the highest power is achieved when only one variable is dominating the fluctuations and the cross terms can be ignored. 12 Note that the values of T S we quote throughout section 5 are averaged only over the neutral IGM component (T S is undefined for ionized gas). Because in the standard picture reionization is approximately "inside out" on large scales, averaging over the neutral IGM means that the T S limits are slightly biased towards underdense volumes.
14 The HERA Collaboration The left and right panels correspond to z = 7.9 and 10.4. The two modes discussed above are clearly seen to emerge by z = 8; at present, the lower-redshift data provide most of the constraining power. At z ≈ 8 we see that HERA-disfavored models have low spin temperatures: T S ∼ < 3 K (or more generically for any uniform radio background T S /T radio ∼ < 0.1 for 0.1 ∼ < x Hi ∼ < 0.9. These constraints are somewhat tighter than analogous ones based on recent LOFAR (Mertens et al. 2020) and MWA (Trott et al. 2020) upper limits: T S ∼ < 2-2.5 K, over narrower ranges in x Hi (c.f. Fig. 4 in Greig et al. 2021b and Fig. 6 in Greig et al. 2021a). Thus, as expected from the stronger PS upper limits, the H21 limits rule out more models than previous power spectrum limits. Furthermore, the density-driven modes were not ruled out by the previous LOFAR limits, which had a larger amplitude and were performed at a higher redshift (z = 9.1; at which the adiabatic-cooling temperature is larger by a factor of (1 + 9.1) 2 /(1 + 7.9) 2 ). At x Hi ∼ > 0.9, the range of temperatures for the disfavored models broadens. This is due to the negative contribution of the ionization-density cross power term, that dominates the large-scale 21-cm power in this regime (Lidz et al. 2008;Zahn et al. 2011). The first galaxies drive HII regions that are very biased in the early stages of the EoR. These quickly cover up the largest matter overdensities, which had earlier dominated the 21-cm power spectrum. Thus for models with negligible temperature fluctuations, the large-scale power drops in the early stages of the EoR before rising again as it transitions from being sourced by the matter fluctuations to ionization fluctuations.

How do the HERA limits improve upon previous complementary data?
As already mentioned, many of the models that are disfavored by the current HERA limits are already inconsistent with existing observations of the z > 6 Universe. Here we put the HERA constraints in context with these other observations by computing the Bayesian posterior over our parameter space with and without the new HERA limits. In particular, we run two inferences (see also section 5.2): • without HERA: This run corresponds roughly to our current state of knowledge, without including 21-cm observations. 13 As detailed in section 5.2, the likelihood incorporates observations of (i) the 13 Here we restrict ourselves to arguably the most modelindependent EoR constraints. In the future, as the 21-cm data improves, we will fold in additional constraints from Lyman-α emitting galaxies, QSO damping wing analysis, opacity fluctuations in the Lyman forests, the patchy kinetic Sunyaev-Zeldovich effect (e.g. Stark et al. 2010;Schenker et al. 2012;Pentericci et al. 2014;Mason et al. 2017;Becker et al. 2015;Bosman et al. 2018;Bañados et al. 2018;Wang et al. 2021;Reichardt et al. 2021). These require more subtle modeling of associated systematics, but could have a non-negligible impact on the recovered EoR history (e.g. Greig & Mesinger 2017b;Dai et al. 2019;Qin et al. 2021a;Choudhury et al. 2021a). We also do not include previous 21-cm upper limits from MWA and LOFAR since these weaker PS limits would not change our with HERA posterior (see the PS evolution inset in Fig. 6). Thus by comparing without HERA and with HERA, we highlight the impact of 21-cm measurements.
galaxy UV luminosity functions at z=6−10; (ii) the upper limit on x Hi at z ∼ 5.9, inferred from the Lyman forest dark fraction; and (iii) the CMB optical depth τ e .
• with HERA: Here the likelihood is computed using both the complementary observations in without HERA above, as well as the HERA limits from Bands 1 and 2. Specifically, we use the (regular) HERA likelihood as defined in equation (5).
5.4.1. Galaxy properties: disfavoring X-ray faint galaxies with HERA The corner plot of these two posteriors is shown in Figure 6, with tan (purple) denoting without HERA (with HERA). As discussed in detail in Park et al. (2019), we see that current observations (without HERA) already rule out a significant fraction of our prior volume, which highlights the power of our 21cmMC approach's inclusion of complementary galaxy observations. Observations of high-z UV luminosity functions shown in the topmiddle sub-panels of Figure 6 constrain the stellar-tohalo mass relation and its scaling with halo mass (f * ,10 and α * ), as well as place an upper limit on the characteristic turnover scale (M turn ). On the other hand, observations of the EoR timing through the CMB optical depth (c.f. bottom-right sub-panel) and the Lyman forest dark fraction (c.f. upper limit in the EoR history sub-panel) constrain the ionizing escape fraction normalization (f esc,10 ) to within 1 dex and place very weak constraints on its evolution with halo mass (α esc ). Using such complementary observations in the likelihood is especially important when sampling from a high dimensional parameter space with flat priors, for which most of the prior volume is sourced by "extreme" corners of parameter space that are already ruled out by existing observations (as is immediately evident from Fig. 6).
Comparing the without HERA and with HERA posteriors, we see that the H21 limits do not have a notable impact over most of the astrophysical parameter space. The new models that HERA rules out, discussed in the previous section, occupy a modest prior volume. 14 However, note that the three X-ray parameters (L X<2keV /SFR, E 0 , α X ) are largely unconstrained by the complementary observations over our prior ranges, because none of the without HERA observations are 14 We use a narrower prior range on L X<2keV /SFR and Mturn in Fig. 6 compared to Fig. 4. This is because Fig. 6 is a true posterior requiring physically reasonable prior ranges, which we discuss further below when presenting galaxy inference. In contrast, Fig. 4 is only meant to illustrate where HERA-disfavored models are expected to reside in our parameter space.
sensitive to the IGM temperature, the observable most strongly affected by the X-ray emissivity. In this part of parameter space, HERA does have a notable impact by ruling out models with weak X-ray heating, which in our parametrization is predominately determined by the integrated soft-band X-ray luminosity to SFR, L X<2keV /SFR. The exclusion of these models is also evident in the 21-cm panels at the upper right, where the recovered signal ranges decrease significantly when including HERA data.
We show a zoom-in of the marginalized 1D PDFs of L X<2keV /SFR in Figure 7. The marginalized without HERA posterior is consistent with the flat prior over the range shown. Current observations do not constrain this quantity aside from disfavoring extreme values of L X<2keV /SFR ∼ > 10 42 erg s −1 M −1 yr, which is so large that X-rays can significantly contribute to reionization (e.g. Mesinger et al. 2013), making it too early in many models. However, the with HERA posterior is able to rule out the lower end of this range, resulting in a 68% highest posterior density (HPD) confidence interval of L X<2keV /SFR = {10 40.2 , 10 41.9 } erg s −1 M −1 yr. H21 is the first observation to place constraints over this range; the analogous analysis of MWA and LOFAR observations (c.f. Fig. 1 in Greig et al. 2021b and Fig. 2 in Greig et al. 2021a) disfavored models with lower luminosities. 15 In Figure 7 we also compare the with HERA limits with estimates based on high-mass X-ray binaries (HMXBs), thought to be the dominant X-ray sources in high-z galaxies (e.g. Fragos et al. 2013). The left vertical line denotes the average value observed from HMXBs in local, metal-enriched, star-forming galaxies (Mineo et al. 2012; see also e.g. Lehmer et al. 2010). Because the HMXB luminosity increases with decreasing metallicity (e.g. Basu-Zych et al. 2013;Douna et al. 2015;Brorby et al. 2016), we do not expect the first, metal-poor galaxies to sit on the left side of this line. And indeed, this local scaling relation is outside of the with HERA 68% confidence interval; thus HERA data already suggests that the first galaxies were more X-ray luminous than their local counterparts. In contrast, the right vertical line in Figure 7 corresponds to the theoretical result from Fragos et al. (2013) for a metal-free HMXB population, expected to be more representative of the first galaxies. Our recovered 1D posterior of L X,<2keV /SFR supports theoretical predictions (e.g. Fragos et al. 2013) and the observed evolution with metalicity and redshift  2013) for a metal-free HMXB population, expected to be more representative of the first galaxies. HERA is the first observation to constrain the X-ray luminosities of Cosmic Dawn galaxies over this range, disfavoring the values seen in local, metalenriched galaxies at > 1σ. Douna et al. 2015;Brorby et al. 2016;Lehmer et al. 2016), that this quantity increases towards high redshifts. Finally, we caution that our limits on L X<2keV /SFR could weaken if alternate heating mechanisms play a significant role. Although we include adiabatic, ionization, X-ray and Compton heating/cooling, in some extreme models alternate heating sources could dominate. These could include shock heating (e.g.  Chuzhoy & Shapiro 2007). However, the amount of heating required by the HERA limits at z ≈ 8 is generally beyond what most of these alternate sources can achieve without violating constraints from other highz observations in our model likelihood. For example, Lyman-α heating only dominates for a relatively large, slowly evolving star-formation density coupled with a low X-ray efficiency. This region of parameter space is ruled out by the combination of complementary observations and HERA limits (e.g. compare the narrower range of our with HERA posterior in the top right pan-els of Fig. 6 to the range of blue curves in Fig. 10 of Reis et al. 2021). Thus it is unimportant for the dataconstrained with HERA posterior in this section, though it can be important in ruling out extreme models when not considering complementary observational data (c.f. Reis et al. 2021 and §8).

IGM properties: disfavoring a cold IGM with HERA
In Figure 8 we show the marginalized without HERA and with HERA posteriors in the space of (x Hi , T S ) (tan and purple regions, respectively). In gray we also show the prior distribution over this space. Comparing the tan to the gray regions, we see that previous observations disfavor a notable prior volume also in the space of IGM properties. 16 Most notably, current observations shift the posterior so that the midpoint of the EoR is occurring around z ∼ 8 to match EoR constraints from Planck and QSO spectra. Now introducing the H21 limits with the purple curves, we see that the HERA disfavors this region of low temperatures for 0.4 ∼ < x Hi ∼ < 0.8 at z = 8. These are the previously-mentioned "reionization-driven" models: having large fluctuations in the ionization field combined with a cold IGM. The impact of HERA is most strongly seen in the marginalized temperature PDFs in the right side panel: with HERA and without HERA exhibit qualitatively different distributions, with the HERA limits strongly disfavoring the low T S peak seen in the posterior without HERA. This demonstrates that the HERA limits are ruling out otherwise viable models.
In Figure 9, we further investigate the physical origins of the temperature PDFs, plotting the spin temperature distributions in the bottom panel and the corresponding kinetic temperature distributions in the top panel. Both are averaged only over the neutral IGM, specifically those cells with x Hi > 0.95. We see that the kinetic temperature of the neutral IGM smoothly extends to T K ∼ 10 4 K, without the bimodality seen in the spin temperature distributions for without HERA. This is because the spin temperature is inversely weighted between the kinetic (T K ) and radio background (T rad ) temperatures (c.f. eq. 2). As T K → ∞, the spin temperature asymptotes to T S → (1 + x α )T rad ∼ (1 + x α ) × 24 K for the standard assumption of a CMB-dominated radio background at z = 8. Although x α scales with the Lyman-α background, it cannot exceed values of 16 We note that our priors over galaxy parameters do not translate into flat priors over (x Hi , T S ). It is easier to theoretically and empirically motivate priors on (fundamental) galaxy properties than on (derived) IGM properties. Thus choosing flat priors directly over mean IGM properties could result in biased posteriors when using weakly constraining data (e.g. Ghara et al. 2020Ghara et al. , 2021.

18
The HERA Collaboration x α ∼ 300 in our data-constrained models of f esc without the gas in the simulation cell becoming ionized. This results in the sharp upper limit of T S ∼ < 600-10 3 K for the neutral IGM seen in Figure 9. 17 Comparing the purple and the tan curves in Fig. 9, we reach the main conclusions of this subsection. H21 observations substantially improve our understanding of the z = 8 neutral IGM temperatures, 18 allowing us to place 68% (95%) high posterior density confidence intervals on the spin temperature of 27 K < T S < 630 K (2.3 K< T S < 640 K) and the kinetic temperature of 8.9 K <T K < 1.3 × 10 3 K (1.5 K <T K < 3.3 × 10 3 K). Other observations of the early Universe and high-z galaxies are unable to constrain these temperatures on the low end.
Indeed because these temperature constraints of the neutral IGM come almost exclusively from the 21-cm signal (where they depend only on the ratio T S /T rad ; 17 Indeed the marginalized prior on T S (shown with the gray curve in the bottom panel of Figure 9) extends out to T S ∼ 10 4 K as the prior volume includes low values of fesc that do not reionize the Universe. Observations exclude these models from the posterior. 18 We want to re-emphasise that our temperatures are averaged over c.f. eq. 1), we can generalize our temperature limits for any homogeneous radio background even if the standard assumption of T rad = T CMB is incorrect. In the regime of T rad > T CMB our with HERA limits can thus be generalized as 1.1 (0.095) < T S /T rad < 26 (26) and 0.37 (0.062) < T K /T rad < 54 (140) at 68% (95%) confidence.

CONSTRAINTS ON IGM PROPERTIES USING A REIONIZATION-DRIVEN PHENOMENOLOGICAL MODEL
Here we introduce simple, phenomenological models for reionization-driven 21-cm power spectra and compare the resulting constraints on IGM properties to those obtained with 21cmFAST in the previous section. Although very simple, these phenomenological models help build physical intuition for the most important effects to consider when interpreting upper limits on the 21-cm power spectrum. We summarize the functionality of this model briefly below, and we defer a more complete description to Mirocha et al., in preparation. Our principal goal is to examine a model built directly from IGM structures rather than galaxy models, so that we do not make any explicit assumptions about the heating and ionizing sources during reionization. To that end, we parameterize the process not with galaxy properties but with the IGM temperature and with the ionized bubble size distribution (BSD) directly. Note that this approach does make implicit assumptions about the sources of reionization, e.g., through the assumed BSD parameterization; it is just non-trivial to determine what these assumptions are. However, they are certainly different from physical models like 21cmFAST, and as a re- sult, help to determine how robust IGM constraints are to modeling assumptions. For an idealized two-phase IGM in which the BSD is known, the two-point statistics of the ionization field can be worked out analytically following Furlanetto et al. (2004b). In 21cmFAST and similar models, the excursion set approach is used to forward model the BSD, but we parameterize it more flexibly here with a log-normal distribution and allow the characteristic bubble size, R b , and the width of the distribution, σ b , to vary as free parameters. Note that BSDs derived from excursion set or semi-numerical models generally have broader tails to low R b than even a log-normal (Furlanetto et al. 2004a;Paranjape & Choudhury 2014b;Ghara et al. 2020), but for fits to a single k and a wide prior on σ b , as we perform here, we do not expect the detailed shape of the BSD to be important. We further assume that the "bulk IGM" outside of bubbles is fully neutral and is of uniform spin temperature, T S . The fourth and final free parameter is the volume-filling fraction of ionized gas, Q HII ≡ 1−x Hi , which normalizes the BSD.
To model the 21-cm power spectrum within this simplified framework, one must model the ionization field and its correlation (or anti-correlation) with the density field. Because we abstract away assumptions about astrophysical sources completely, and instead work in terms of the BSD and mean IGM properties Q HII and T S , it is not immediately obvious how to do this. While it is possible to estimate the behaviour of cross-terms using the halo model (Furlanetto et al. 2004b) or perturbation theory ), here, we take a simpler approach that avoids explicit assumptions about astrophysical sources. If we assume for simplicity that the structure of the density field mirrors that of the ionization field, i.e., it is a binary field, cross-terms involving ionization and density can be re-written in terms of the ionization power spectrum given the typical density of ionized regions, δ i . To estimate δ i , we assume that reionization is "inside-out," or in other words that the ionized volume fraction Q HII is made up of the densest fraction Q HII of the volume. Then, to complete this "volume matching" procedure, we assume the density PDF is log-normal (Coles & Jones 1991) with a variance given by the density field smoothed on the scale at which the BSD peaks. This naturally leads to a model in which the typical bubble density declines with time, so the importance of cross-terms is greatest in the early stages of reionization. Finally, as in §4, we assume µ = 0.6 to match the spherical averaging done in 21cmMC simulations.
We perform two MCMC fits using emcee (Foreman-Mackey et al. 2013) -one using the inverse likelihood (eq. 10) and one with the regular likelihood (eq. 9) -to the k = 0.134 cMpc −1 limit from Band 2 at z = 7.9 using 192 walkers for a total of ∼500,000 steps. We adopt flat priors on each model parameter: 0 ≤ T S /K ≤ 10 3 , 0 ≤ Q HII ≤ 1, 0 ≤ R b /cMpc ≤ 30, and 0.5 ≤ σ ≤ 2. Note that while the 21-cm signal is insensitive to T S once T S T γ , our lower limits on T S are sensitive to the prior range. Our choice of T S ≤ 10 3 K is motivated by the maximum allowed spin temperature in standard scenarios (see section 5.4 and Fig. 9), though we broaden the lower bound from the expected adiabatic cooling limit of T S 1.7 K to zero so that more exotic scenarios may be considered.
In the top panel of Figure 10, we show constraints on the mean spin temperature and ionized fraction of the IGM obtained from this model after marginalizing over the parameters of the bubble size distribution (R b and σ b ). We obtain 95% (68%) lower limits on the spin temperature of the z = 7.9 IGM of T S ≥ 5.3 K (T S ≥ 25 20 The HERA Collaboration K). Qualitatively, these results are in good agreement with those derived using 21cmMC (the with HERA posterior is shown with purple contours; see also Fig. 8). As discussed in the previous section, the data-constrained 21cmMC posterior is dominated by "reionization driven" fluctuations, since the "density driven" models have a neutral fraction at z ∼ 8 that is too high and are disfavored by EoR observations. It is therefore encouraging that our "reionization driven" phenomenological model is broadly consistent with the with HERA posterior from 21cmMC. This implies that our claims of HERA's upper limits disfavouring models in which the IGM has not been heated at z ∼ 8 are not sensitive to the nature of the EoR fluctuations.
In the bottom panel, we show the results obtained via the inverse likelihood. Note again that we require only that the mean temperature of the IGM be positive, which is why the disfavoured region in this panel extends all the way to T S = 0. This is one of the advantages of the phenomenological approach: it can constrain more exotic scenarios without invoking a particular physical model (c.f. section 7, where we introduce some such physically motivated models). Here again, our results are broadly consistent with the analogous ones from 21cmMC (the "reionization driven" modes are shown with the blue curves; note the red contours are "density driven" modes that are not considered by our phenomenological BSD model).
To further explore this agreement, in the bottom panel of Figure 10 we show iso-power contours for several different bubble sizes, holding fixed the width of the BSD at σ b = 0.8. The rationale here is simple: isolikelihood contours should trace iso-power contours for inference based on a single k mode. From this plot we see that if bubbles are generally small, R b 2 cMpc, the phenomenological model predicts that warmer temperatures are needed to preserve the large-scale power as x Hi → 1. However, if bubbles are generally larger, with R b 2 cMpc, this trend is reversed. These results suggest that physical models like 21cmFAST effectively have a low prior probability assigned to models with large bubbles at early times. Indeed, excursion set calculations suggest that typical bubbles sizes R b 2 cMpc generally do not emerge until reionization is underway at the ∼ 20 − 30% level (Furlanetto et al. 2004b). However, because the phenomenological model can have arbitrarily large bubbles at any time, the density-driven mode is washed out when marginalizing over R b and σ b . Though the "density-driven" models are ultimately disfavoured given that they do not complete reionization by z ∼ 6 (see Fig. 4), they serve as interesting test case nonetheless (see section 4). Figure 10. Constraints on the mean properties of the z ∼ 8 IGM using phenomenological models (see §6) compared to the 21cmFAST results. Top: Filled cyan contours are 68% (dark) and 95% (light) confidence levels obtained with the phenomenological model, while purple contours are those from 21cmMC (as in Fig. 8). Bottom: As in Fig. 6, the blue and red open contours in each panel correspond to reionization-and density-driven scenarios, respectively, while the filled contours show the disfavoured region determined by the phenomenological model. Additional black contours in this panel trace the phenomenological model's predictions at fixed k = 0.134 power, ∆ 2 21 = 24.07 2 mK 2 , corresponding to the HERA measurement +1σ, for three different bubble sizes assuming a fixed σ b = 0.6. The cross-hatching along the bottom of each panel indicates regions with temperature below that of a homogeneous and unheated high-z IGM.

CONSTRAINTS ON DARK MATTER AND ADIABATIC COOLING USING DENSITY-DRIVEN MODELS
In previous sections we have obtained limits on the IGM spin temperature T S using different approaches. Here we study how these limits compare to predictions in the standard CDM cosmological model, as well as models of millicharged DM (mQDM). We will also briefly explore how our assumptions about RSDs affect the limits imposed on the IGM. Throughout this section we will use our density-driven phenomenological model, equation (12) (in all cases assuming T S = T K ). While this approach has limited validity, it provides a useful test bed of our assumptions, as it allows us to obtain analytic limits under different RSD assumptions, as well as (f dm = 0.5%) Figure 11. Analytically derived lower limits on T S/Trad from HERA data (95% confidence, green arrows) compared to the upper limit from EDGES (blue arrow, which also implies the two lower limits at z = 15 and 21 given their profile shape). The H21 limits have been obtained by assuming density-driven fluctuations (see section 4 for details) and two different constant values of the neutral-hydrogen fraction xHi = 1 (dark green) and 0.5 (light green). The black line shows the standard CDM prediction assuming full coupling (TS = TK ), with the dashed line corresponding to zero heating and the solid line to a toy model of X-ray heating. The HERA Band 2 data can rule out adiabatic cooling under our assumptions, requiring some X-ray heating to take place before z = 7.9. The red line includes a fraction f dm = 0.5% of millicharged DM (mQDM), so as to explain the EDGES depth. Without any heating (dashed), this curve is ruled out by HERA Band 1, showing that there must be heating between z = 17 and z = 10.4 if EDGES is explained by mQDM. The conclusion is similar for a model with excess radio emission (with a radio fraction fr = 9000, see section 8), shown as the tan line. The empty symbols represent spherically averaged RSDs (purple, which were shown in Fig. 2). extend the temperature range studied below the adiabatic cooling threshold to probe mQDM models, neither of which are currently included in the usual 21cmFAST simulation-based approach of section 5.

The Impact of RSDs on the T S Limits
First we study how our analytic limits change under different RSD assumptions. Within our bias approach this can be readily implemented by varying µ in equation (12). For simulations, on the other hand, it is challenging to study the µ → 1 limit, given the geometry of the Fourier modes populating a square box. The analytic limits obtained in section 4 (T S ≥ {7.8, 1.9} K at z = {7.9, 10.4}), assumed µ = 0.6 to match the spherical averaging done in 21cmMC simulations. Under the assumption that modes lie predominantly along the line of sight (µ ≈ 1), as actually observed by HERA (see section 3.1), these limits strengthen to T S ≥ {11, 2.6} K (for z = {7.9, 10.4}) at 95% confidence, which are ∼ 50% stronger, as shown in Figure 11. If we had ignored RSDs (µ = 0), but kept the same assumptions otherwise, the 95% CL limits would weaken to T S > {3.1, 0.74} K at z = {7.9, 10.4}, a factor of ∼ 3 smaller. The difference between these three assumptions highlights the importance of properly modelling RSDs in 21-cm powerspectrum analyses. We note, however, that these results assume the density field drives the 21-cm fluctuations, in which case RSDs always increase the 21-cm power spectrum. This trend can be reversed if radiation fields are the main source of anisotropy (e.g. in "reionization driven" scenarios as in Fig. 3), though it is not expected to change our conclusions, see Sec. 5.2.
We also show the impact of varying the neutralhydrogen fraction x Hi on our analytic results. Unlike the galaxy-driven models of previous sections -in which patchy reionization enhances the 21-cm power spectrum because of the bubble structure -here we assume uniform reionization (which could result from exotic processes; e.g. Evoli et al. 2014;, in which case x Hi < 1 suppresses the 21-cm power spectrum, as clear from equation (12). Had we assumed x Hi = 0.5 (instead of fixing x Hi = 1), we would arrive at the 95% confidence limits T S > {4.1, 1.2} K at z = {7.9, 10.4} (both with µ ≈ 1). While it is unlikely that x Hi deviates significantly from unity at z = 10.4, a global value of x Hi = 0.5 is in line with our expectations for z = 7.9.
These limits have interesting implications for the thermal state of the IGM at high redshifts, as well as for the first EDGES claimed detection (Bowman et al. 2018). We compare all the T S limits (divided by T CMB ) in Figure 11 against the T S /T rad prediction for the standard CDM model, both in the absence of heating and with a fiducial X-ray heating model, akin to the ones implemented within 21cmFAST in previous sections. The HERA Band 2 95% confidence limit is above the adiabatic cooling prediction at z = 7.9, both for x Hi = 0.5 and 1 (and in fact for any x Hi ≥ 0.3 in this bias approach). Thus, HERA requires some heating by z = 7.9 given our assumptions. Moreover, the HERA limits for Band 1 (z = 10.4), while below the adiabatic limit at that z, can be used to set constraints on dark-matter induced cooling of the gas, which we now explore.

Dark matter-baryon interactions
The first claimed 21-cm detection from the Cosmic Dawn in Bowman et al. (2018) shows a surprisingly deep absorption feature at z ∼ 17. The depth of this absorption, if interpreted to be cosmological (see, however, e.g. Hills et al. 2018;Sims & Pober 2019;Tauscher 22 The HERA Collaboration Hatched regions are ruled out by different experiments, in brown we show limits from the SLAC experiment (Prinz et al. 1998), in orange the most conservative BBN constraints (Jaeckel & Ringwald 2010), in red the constraint from cooling of SN1987A from Chang et al. (2018) (see, however, Bar et al. 2020 for criticisms), and in cyan the region disfavored if there is a dark photon A mediating the millicharge ) (see also Vogel & Redondo 2014). Blue shows the EDGES-preferred region (z ≈ 17, following , and the green region is ruled out by HERA band 1 (z = 10.4) at 95% confidence, assuming density-driven fluctuations. An EDGES detection of millicharged DM is only compatible with HERA if heating takes place between z ≈ 17 and z = 10.4. We have taken a fraction f dm = 0.5% of DM to be millicharged, but our conclusions extend to all relevant fractions. The red star is the point that gives rise to the red line in Fig. 11. We remind the reader that the HERA Band 2 data (z = 7.9) already rules out adiabatic cooling at 95% confidence, so by construction it also rules out any DM model that produces additional cooling.
et al. 2020), can be translated into a requirement that T S /T rad ≤ 0.08 at z = 17, a factor of two smaller than allowed by the standard cosmological model. Reducing the Wouthuysen-Field coupling in this case only exacerbates the tension, as it would bring the spin temperature closer to that of the CMB, producing shallower absorption. A possible explanation for this anomalous depth consists of lowering the temperature of baryons in the IGM by allowing them to interact with the cosmologically abundant-and kinetically cold-DM. Elastic scattering between these two fluids would bring them closer to thermal equilibrium, cooling down the baryons and heating up the DM. These interactions could take the form of a new fundamental force (Tashiro et al. 2014;Muñoz et al. 2015b;Barkana 2018;, which however would be in conflict with fifth-force constraints and stellar-cooling bounds. Alternatively, part of the DM can be electrically charged, for instance through a dark-photon portal (Holdom 1986), a scenario dubbed millicharged DM (mQDM). In this case there are no new charges for baryons, and therefore fifth-force and stellarcooling bounds are naturally evaded Barkana et al. 2018;Slatyer & Wu 2018;Berlin et al. 2018;Kovetz et al. 2018;Liu et al. 2019). Here we briefly study how well DM-baryon interactions, in the form of mQDM 19 , can be constrained by the H21 limits.
To illustrate the effect of mQDM, we show in Figure 11 the prediction for an example model using the software developed in . We solve the coupled differential equations for the mQDM and hydrogen-gas temperatures starting at recombination. The interactions due to millicharges produce thermalization of the (initially cold) DM and the hydrogen gas, therefore cooling the latter. For this figure we have chosen mQDM with a charge Q χ = 1.3 × 10 −5 e, where e is the electron charge, and mass m χ = 10 MeV, composing a fraction 0.5% of the total DM. These parameters are chosen to (barely) explain the EDGES depth and, as is clear from the figure, the cooling induced at later times lowers T S below the HERA limit both at z = 10.4 and z = 7.9, ruling out this model in the absence of heating.
We generalize this result by performing a 2D scan of mQDM charges Q χ and masses m χ , assuming that mQDM particles compose a 0.5% fraction of the DM, which is at the edge of the 95% confidence interval region allowed by CMB constraints (Boddy et al. 2018), with the remainder being neutral and non-interacting CDM. We show the results in Figure 12, where we also show the region that produces enough cooling to explain the EDGES depth (defined to be T S < 4 K as in . This region is entirely contained by the HERA Band 1 constraint (T S > 2.6 K at 95% confidence), which shows that all the mQDM models that explain the EDGES depth also require heating before z = 10.4 in order to avoid conflict with HERA. Our conclusions hold for all other mQDM fractions in the relevant range f dm = 0.05-5%.

ASTROPHYSICAL CONSTRAINTS IN MODELS
WITH AN EXTRA RADIO BACKGROUND

The 21-cm Signal in the Presence of Radio Sources
In this section we use HERA data to constrain models in which either astrophysical or exotic high-redshift radio sources contribute to the total radio background, in addition to the CMB. Such an excess radio background above the CMB level has been observed at z = 0, with the data consistent with a synchrotron radio background of a spectral index −2.58 and a brightness temperature ∼ 603 K at the rest-frame 21-cm frequency Seiffert et al. 2011;Dowell & Taylor 2018). The nature of this excess is still undetermined (e.g. Subrahmanyan & Cowsik 2013), and it could partially be accounted for by a population of unresolved high-redshift sources of either astrophysical or exotic origin (Ewall-Wice et al. 2018b;Jana et al. 2019;Fraser et al. 2018;Pospelov et al. 2018;Brandenberger et al. 2019;Thériault et al. 2021).
An excess high-z radio background would have important implications for the 21-cm signal, because the stronger background amplifies the absorption (via the temperature term in eq. 1 including an effect on coupling coefficients in eq. 2; see complete discussion in Fialkov & Barkana 2019;Reis et al. 2020). Such models have been presented as potential explanations of the anomalously strong EDGES Low Band detection (Bowman et al. 2018); for example, Fialkov & Barkana (2019) found that the EDGES signal can be explained if the cosmological (high redshift) contribution of such a background is between 0.1% and 22% of the CMB at 1.42 GHz (see also Mirocha & Furlanetto 2019;Jana et al. 2019;Ewall-Wice et al. 2018bMebane et al. 2020;Reis et al. 2020;Thériault et al. 2021). These explanations are challenging, however, requiring either unconfirmed exotic sources or astrophysical sources that are far stronger than expected based on local observations (Ewall-Wice et al. 2018b;Mirocha & Furlanetto 2019;Mebane et al. 2020;Ewall-Wice et al. 2020) and necessitating rapid X-ray heating to match the steep recovery in the EDGES signal (Reis et al. 2020).
More interestingly for our purposes, the presence of a radio background can also enhance fluctuations in the 21-cm signal (Fialkov & Barkana 2019;Reis et al. 2020), so that H21 can place limits on such a background (whether generated by discrete sources or more exotic processes) at z ∼ 8 and 10. In this section we will consider both such scenarios, including the resulting limits in the context of other observations of the low-frequency radio background Dowell & Taylor 2018) and X-ray background (Lehmer et al. 2012).

Modelling
We generate spherically averaged 21-cm power spectra as a function of astrophysical parameters using the semi-numerical simulation method described in Visbal et al. (2012);Fialkov et al. (2014); Fialkov & Barkana (2019); Cohen et al. (2020); Reis et al. (2020Reis et al. ( , 2021. Our simulations are 384 cMpc on a side and have a resolution of 3 cMpc. Initial large-scale perturbations in density and relative velocity between dark matter and gas ) are linearly evolved from the Dark Ages (z ∼ 60) to z = 5. Using the modified Press-Schechter mass function which takes into account the effect of large-scale overdensity and velocity fields (Barkana & Loeb 2004;Fialkov et al. 2012a), we calculate the halo abundance in each voxel of the simulation. Each halo is then populated by stars, and emissivities in different bands are calculated (see, e.g. Cohen et al. 2020, for details). RSDs are computed by multiplying the real space isotropic 21-cm signal by (dv r /dr) −1 which is the radial component of the velocity gradient created by structure formation . Using coeval simulation cubes we calculate the spherically-averaged power spectrum at every redshift.
The key radiation backgrounds are all driven by the cosmic star formation rate, which in the simulations depends on two parameters. First, we choose a minimum circular velocity of star forming halos V c , which determines the halo population that can form stars. We then choose a star formation efficiency f * , which measures the amount of collapsed gas that turns into stars for halos above the atomic cooling limit, imposing an extra suppression in smaller halos (e.g. Cohen et al. 2020). The code includes the suppression of star formation by Lyman-Werner feedback , relative velocities between dark matter and gas (Fialkov et al. 2012a) and photoheating feedback (Cohen et al. 2016). To calculate the Ly-α background, we assume Pop II star formation following Barkana & Loeb (2005b). For completeness we note that here we include multiple scattering of Ly-α photons and Poisson fluctuations in the number of first star forming halos (however, these effects are not significant at the redshift range observed by HERA, Reis et al. 2021).
X-ray heating of the IGM is powered by a population of X-ray binaries with the ratio of bolometric luminosity to SFR of calculated between 0.2 keV and 95 keV assuming a hard X-ray SED of X-ray binaries . The free parameter f X normalizes the X-ray efficiency relative to a population of present day binaries (but including an order-of-magnitude increase in this ratio at the 24 The HERA Collaboration low metallicity expected for high-redshift galaxies, Fragos et al. 2013). The assumed SED is relatively hard, peaking at ∼ 1 keV Fialkov et al. 2014). Note that assuming a different SED could affect our final results (Fialkov et al. 2014;Pacucci et al. 2014;Das et al. 2017;Reis et al. 2020). The unresolved X-ray background observed by Chandra imposes an upper limit on this contribution , as we will explore later, but we allow a broad range of f X between 10 −4 and 10 3 in the estimation framework. This model also includes heating by scattering of Lyα photons (Chen & Miralda-Escudé 2004;Chuzhoy & Shapiro 2007;Ciardi et al. 2010;Mittal & Kulkarni 2021;Reis et al. 2021) and the CMB (Hirata & Sigurdson 2007;Venumadhav et al. 2018b;Fialkov & Barkana 2019;Reis et al. 2021, though see Meiksin 2021. With the onset of the first stellar population, the extra heating processes raise the IGM temperature above the adiabatic limit even in the absence of X-ray heating, reducing the 21-cm background at the relevant redshifts in some scenarios by a factor of a few (see Reis et al. 2021, for more details).
The process of reionization is implemented using the excursion set formalism (Furlanetto et al. 2004a) and is described by two parameters: the ionizing efficiency of sources ζ, which is normalized via the total CMB optical depth τ , and the horizon of ionizing photons, R mfp . Although the latter parameter does affect the intensity of the 21-cm fluctuations at the end of the EoR, we fix it at 40 Mpc here as it plays a secondary role in our constraints.
Finally, we explore two types of radio backgrounds (beyond the CMB): 1. A fluctuating, time-variable radio background generated by galaxies, parameterized by f r . We assume that the galaxy radio luminosity per unit frequency in units of W Hz −1 is proportional to the SFR (following Mirocha & Furlanetto 2019) where α r is the spectral index in the radio band, which we set to the typical value of 0.7 (Mirocha & Furlanetto 2019), which is compatible with observations (Hardcastle et al. 2016;Gürkan et al. 2018). We calculate T rad at redshift z by summing up over the past light-cone contribution of all the radio galaxies (see Reis et al. 2020, for more details).

A smooth synchrotron background that decays with time, for which we replace T CMB by
where ν obs is the observed frequency, A r is defined relative to the CMB temperature and β = −2.6 is the spectral index in agreement with the AR-CADE2 ) and LWA1 (Dowell & Taylor 2018) observations. Here we treat this background as phenomenological, but it could have been produced by exotic radio sources, e.g. radiative decay of relic neutrinos into sterile neutrinos (Chianese et al. 2018), light dark matter decays (Fraser et al. 2018;Pospelov et al. 2018) and superconducting cosmic strings (Brandenberger et al. 2019).
In this work, we allow a broad range of f r and A r parameters. However, as we discuss later (and as was shown by Fialkov & Barkana 2019;Reis et al. 2020) models with strong radio backgrounds, e.g. f r × f * > 10 3 − 10 4 for the radio from galaxies, are constrained by AR-CADE2/LWA1 data.
To summarize, the models considered here include an extra radio background in addition to the CMB either produced by radio galaxies or emitted by exotic sources. Our models build on the following parameters: V c varied between 4.2 and 100 km s −1 , f * between 0.001 and 0.5, f X in the range between 10 −4 and 10 3 , τ between 0.035 and 0.088, f r from 1 to 10 5 , and A r between 10 −2 and 10 5 . Owing to the large dynamic ranges, we assume uniform priors on the parameters log 10 f * , log 10 V c , log 10 f X , τ , and log 10 f r or log 10 A r . For completeness we also include our constraints on the standard models (i.e. with no extra radio background above the CMB).
Although similar in spirit to the 21cmFAST simulation suite described in section 5, there are differences between the two sets of simulations. We refer the reader to the relevant papers for details on the physics and implementation differences. Broadly, the simulations described in this section include some additional heating processes, such as Lyman-α heating (which can affect the IGM temperature when it is very cold in some models) and (of course) radio emission, but they have a more prescriptive source model with fewer free parameters (and they are not constrained by ancillary observations such as the galaxy luminosity functions). A detailed code comparison is beyond the scope of this paper; instead, we focus on how these distinct codes can address the issues to which they are each best suited.

Parameter estimation
We explore the parameter space of models compatible with the HERA upper limits based on the likelihood L m defined in equation (9). We also decimate the power spectra as described in Section 3.2.1, using the even wave numbers (k = 0.086, 0.17, ... cMpc −1 ) of Band 1 and odd wave numbers (k = 0.13, 0.21, ... cMpc −1 ) of Band 2. Note that we only compute the power spectrum up to k = 1.1 cMpc −1 , limited by the simulation resolution. Larger wavenumbers (smaller scales), however, are irrelevant as the HERA limits rise much steeper at small scales than realistic models so those scales do not contribute towards the constraints.
Because individual simulations take a few hours to complete, we instead use an emulator to interpolate the power spectra from ∼ 10 4 existing simulation runs (for each of the two types of radio background that we investigate here). We implement the emulator using neural networks: taking the astrophysical parameters as an input, a network predicts the logarithm of the power spectra for the HERA redshifts (z = 7.93 and 10.37) and wave numbers (from k = 0.086 to 1.1 cMpc −1 ). The architecture of the emulator includes a multi-layer perceptron with 4 hidden layers of 100 nodes each, implemented using scikit-learn (Pedregosa et al. 2011). The power spectra are predicted with a relative error of 20%; we take this uncertainty into account by adding it in quadrature to the observational error σ i in equation (9). Although this is an approximation, the associated error is negligible in the context of current analysis. A detailed discussion of the emulator and its accuracy can be found in Appendix B.
We explore the parameter space using the MCMC Ensemble sampler emcee (Foreman-Mackey et al. 2013), and we visualize and analyze the results using anesthetic (Handley 2019) and GetDist (Lewis 2019). 8.4. Results: A radio background generated by galaxies As was alluded to above, models with an additional radio background can easily, unlike most standard scenarios, exceed the HERA upper limits. To illustrate this point, we show a random subset of the simulated power spectra for the case of a radio background from galaxies in Figure 13. The power spectra are shown at z = 7.93 and colored with respect to their compatibility with HERA constraints, indicating the difference in loglikelihood ∆ log L m compared to the best fit (which is ∆ 2 21 ≈ 0 mK 2 ). For comparison, we also plot the current HERA limits marked by data points with error bars. Clearly, a substantial fraction of the models (shades of orange) are excluded by the current HERA limits with high significance. In comparison, corresponding standard models (no additional radio background) have Figure 13. Power spectra of 1000 randomly selected models with an extra radio background created by galaxies, at z = 7.93 (HERA Band 2). A number of these models can be ruled out by the HERA data as shown by the color of the lines indicating the likelihood of each model. We use the decimated data points as described in section 3.2.1, taking only the "even" points of Band 1 and "odd" points of Band 2. We show the Band 2 points with black circles for the points we take into account and white circles for the unused data. The error bars show 1σ errors and the crosses indicate the 2σ upper limits. For comparison, the thick dashed line shows the maximal possible amplitude of the ensemble of standard models at each wave number k (i.e. the envelope). much lower amplitudes. We show the envelope of these models (i.e. the maximal possible amplitude of the ensemble of standard models at each k) with the thick dashed line in Figure 13 and discuss astrophysical implications of HERA for these cases in section 8.6.
Using the HERA likelihood alone (L m ) we show the marginalized constraints on the parameters in Figure  14. The diagonal panels show the 1D marginalized posterior PDFs while the others show the 2D marginalized PDFs with the dotted and dashed lines indicating the 68 and 95% confidence contours (containing 68 and 95% of the 2D posterior probability, respectively). The 2D marginalized posteriors involving f * , V c , and τ are relatively flat with the ratio between minimum and maximum posterior probability 20 being between 0.3 and 0.6. This results in confidence contours which could be easily affected by fluctuations due to the random sampling and 20 These are measured from the bins in Figure 14 using the minimum/maximum bin sample count (∝ posterior value). A value close to one shows that the PDF is largely flat and therefore does not provide a strong constraint.

26
The HERA Collaboration are strongly dependent on the prior. On the contrary, we find the 2D posterior in the f X -f r plane to show a strong contrast between minimum and maximum regions (with minimum/maximum posterior ratio of 0.02, i.e. dropping by more than three e-folds). There is a vanishing probability for models with both a strong radio background (large f r ) and weak X-ray heating (low f X ). The large contrast in the probability across this sub-space indicates that the constraints on the combination of f X and f r are expected to be robust, i.e. even for different priors (which would cause a small shift in the contour lines) most of the high-f r low-f X region will still be excluded.
Marginalizing over the rest of the model parameters, we constrain f X to be greater than 0.25 and f r less than 397 at 68% confidence level individually, which maps to the excluded (at 68%) L X, 0.2−95 keV /SFR < 7.6 × 10 39 erg s −1 M −1 yr and L r,ν /SFR > 4×10 24 W Hz −1 M −1 yr (calculated at reference frequency ν = 150 MHz). The 2D region where both f X < 0.25 and f r > 397 (marked by orange solid lines on the corresponding 2D PDF in Figure 14) approximately corresponds to the region excluded at 95% confidence level (exact contours shown as dashed black line).
To relate the constraints on model parameters to the physical state of the IGM at the observed redshifts, we calculate limits on the radiation temperature T rad and gas temperature T K at redshift z = 8 (note that these are derived, i.e. not sampled, parameters in the context of this model). The temperatures are averaged over all the regions that are not fully ionized (i.e. over the volume outside of HII regions). 21 These temperatures are stored for every simulation and emulated similarly to the power spectra (see Appendix B for details). We obtain 2D posterior PDFs for these temperatures, marginalized over all other parameters, and we show the results in the top right panel of Figure 14 (blue contours). Because the sampling is uniform in terms of the simulation parameters (or their log for f * , V c , f X , τ , f r ), the priors on the derived parameters (temperatures) are not uniform and shown as orange contours in the background.
The H21 limits significantly affect the temperature distribution, excluding high values of T rad when T K 1000 K. This is expected, as, typically, a combination 21 For completeness, here the model suite includes scenarios that have extremely strong heating corresponding to f * × f X 1 which are ruled out by Chandra observations. In such scenarios gas outside of the HII regions could be partially ionized up to ∼ 50% by X-ray photons (e.g. Fialkov et al. 2017) which leads to a suppressed 21-cm signal. The maximal possible temperature of such partially ionized gas is set to be the ionizing temperature of hydrogen, 1.57 × 10 5 K, equivalent to the energy of 13.6 eV. of high T rad and low T K leads to a strong 21-cm signal exceeding the observed HERA limits. Exceptions are scenarios in which gas is either significantly ionized already by z = 8 or in which the radio background has a non-negligible effect on the spin temperature T S , driving it towards the background radiation temperature T rad , and, as a result, suppressing the total power spectrum. In such models the 21-cm power spectrum is below the HERA limits even when gas is predominantly neutral and is colder than the radio background. These models populate the lower-left corner of the temperature plot and are not ruled out by HERA even when T K T rad . Overall we derive a constraint on the ratio log 10 (T rad /T K ) < 1.1 at z = 8 (at 95% confidence), compared to the prior constraint of log 10 (T rad /T K ) < 2.6 (95% confidence), noting that these limits can be exceeded at low T K due to the effect described above. The bounds are shown as dashed lines in the temperature plot in Figure 14. Reis et al. (2020) showed that models with a strong radio background, corresponding to large values of the product f * × f r , violate the ARCADE2/LWA1 measurement with only a weak dependence on the values of other astrophysical parameters. Additionally, models with large values of f * × f X are ruled out by the Chandra X-ray background , though in detail the constraints are sensitive to the assumed SED at high energies. To compare HERA constraints to the ARCADE2/LWA1 and Chandra limits we calculate the contributions to the observed present-day radio and Xray backgrounds created by a population of sources at z > 8 and flag as rejected models in which these highredshift contributions exceed the observations. In Figure 15 we show that the H21 limits are complementary to both the ARCADE2/LWA1 and Chandra X-ray background limits in that they rule out a different section of the parameter space, which is characterized by a high radio background and low X-ray heating (orange filled circles).

Synchrotron radio background
Next, we calculate limits on the model parameters for the case of a smooth phenomenological synchrotron radio background (given by eq. 19 with ν obs = 1420/(1+z) MHz) that is stronger at higher redshifts. Here we build on the formalism developed by Fialkov & Barkana (2019) of the excess radio background, taking into account its effect on the background intensity, spin temperature and coupling coefficients, as well as heating by radio/CMB photons. We expand over the model used in Fialkov & Barkana (2019)  HERA constraints on models with extra radio background from galaxies. The triangle plot shows the marginalized 1D and 2D posteriors for all model parameters, with 68% and 95% confidence contours indicated by black dotted and dashed lines, respectively. The backgrounds show the 2D histograms of posterior values to illustrate how strongly a region is excluded, white background indicates a posterior close to zero while dark blue corresponds to the maximum posterior value of a panel.
The grey areas in diagonal panels, and the solid orange lines show the individual 68% confidence limits on fX and fr (numbers quoted in text). Combined, these approximately correspond to our main result, the 95% confidence excluded region in the fX-fr space, as shown in the center bottom panel. The inset in the top right corner shows the corresponding constraints on the derived parameters, T rad and T K , at z = 8 for both prior and posterior. The prior is uniform in the model parameters (log 10 f * , log 10 Vc, log 10 fX, τ , log 10 fr) but therefore non-uniform in the derived parameters, T rad and T K . The blue and orange dashed lines show upper bounds (at 95% confidence) on log 10 (T rad /T K ) from prior and posterior, respectively.

28
The HERA Collaboration 10 4 10 3 10 2 10 1 10 0 10 1 10 2 X-ray heating f X f * . HERA constraints compared to the radio background constraints from LWA1 and the X-ray background constraints from Chandra. The plot color-codes the models by their HERA likelihood Lm (orange corresponds to the models excluded at ∆ log Lm ≥ 20 relative to the best fit and purple are the allowed models), and models excluded by LWA and Chandra are shown as empty circles and crosses, respectively. The latter constraints are computed from the level of X-ray and radio background produced population of sources at z > 8 and approximately correspond to fX ·f * 1 and fr · f * 10 3 , respectively. We see that HERA can exclude a substantial fraction of the otherwise unconstrained parameter space.
fluctuations as described in Reis et al. (2021). We use the HERA likelihood alone and obtain qualitatively similar results (Figure 16) to the ones found in models where the radio background is generated by galaxies. Differences arise because the synchrotron background decays with time, while the radio background from galaxies grows. In addition, the former background is smooth, while the latter traces galaxies and, thus, includes spatial fluctuations.
We find that models with a high excess radio background of A r > 31 (calculated relative to the CMB at the reference frequency of 78 MHz, corresponding to 1.6% of the CMB at 1.42 GHz) are excluded at 68% confidence; and so are models with a low X-ray efficiency of f X < 0.20 (corresponding to the X-ray luminosity per SFR of L X, 0.2−95 keV /SFR < 5.9 × 10 39 erg s −1 M −1 yr). Models fulfilling both criteria are clearly excluded by the data as seen in the f X -A r panel in Figure 16, they mostly lie beyond the 95% 2D confidence contour (dashed line).
We also constrain the ratio log 10 (T rad /T K ) < 1.7 at z = 8 at 95% confidence, compared to the prior constraint of log 10 (T rad /T K ) < 3.5, again noting that these limits are less strict due to the shape of the posterior distribution as shown in the temperature plot in Figure  16.
A similar analysis with these models was done by Mondal et al. (2020) using the upper limit at z = 9.1 derived from 141 hours of LOFAR HB data (Mertens et al. 2020). These data rule out (at 68% confidence) a strong contribution above 0.83% of the CMB at 1.42 GHz (corresponding to A r > 15.9) of the high-redshift Universe to the ARCADE2 and LWA1 measurements as well as scenarios with weak X-ray heating f X < 0.01. Because the assumed synchrotron radio background decays with time, tighter constraints at higher redshifts are expected. However, the theoretical framework used here has been upgraded since the analysis of the LO-FAR data (specifically, RSDs and extra heating by Lyα and radio were added). This difference impedes a quantitative comparison between the two works. Moreover, compared to Mondal et al. (2020), here we used a broader prior range on A r , which could partially explain our higher limit.

Implications for standard models
For completeness we note that only a small fraction of the prior volume used for this simulation suite is ruled out if the only source of the background radio photons is the CMB. We find that every excluded model (at the level of ∆ log L m > 1) belongs to the region of parameter space with f X 1 (corresponding to X-ray luminosity per unit SFR of 3 × 10 40 erg s −1 M −1 yr), V c 50 km s −1 , τ 0.064 and f * 0.07. Note that varying R mfp , which is fixed here to 40 Mpc, could slightly modify the standard model's constraints (see Mondal et al. 2020). The excluded models are characterized by a low X-ray heating efficiency and star formation rate, resulting in a partially reionized IGM (close to 50%) at z = 8 and creating the strongest reionization-driven power spectra at the frequencies observed by HERA. We do not rule out density-driven models with HERA observations, likely due to Ly-α heating which lowers the power spectrum amplitude in those models.

EDGES-motivated scenarios
As mentioned in section 8.1, models with a strong extra radio background were proposed as a potential explanation of the anomalously deep absorption trough detected by EDGES. In such models, an accompanying enhancement of the power spectra (compared to the corresponding cases with the CMB as the background radiation) is expected (Fialkov & Barkana 2019;Reis et al. 2020, see also   Posterior PDF Posterior PDF Excluded region Figure 16. HERA constraints on models with a smooth synchrotron extra radio background, with analogous notation to Figure 14. Dotted and dashed lines indicate 68% and 95% confidence contours, the background shows the histogram of the 2D posterior. The grey areas and orange solid lines show the 1D 68% percent limits on fX and Ar individually, the combination of the latter approximately corresponds to the 95% combined exclusion region, our main result. The upper right inset shows the prior (orange) and posterior (blue) distributions of the derived model parameters T K and T rad at z = 8, with the dashed lines indicating the respective 95% confidence limits on log 10 (T rad /T K ). HERA at k = 0.13 cMpc 1 HERA at k = 0.17 cMpc 1 Figure 17. Models compatible with a cosmological interpretation of the EDGES detection in the context of HERA limits. Top: Global signals of randomly selected models (grey) and those that are compatible with EDGES (purple), selected by their absorption feature close to 78 MHz and the narrowness of the profile. The models here include an extra radio background from radio galaxies. Bottom: Corresponding power spectra at k = 0.13 cMpc −1 . HERA measurements at k = 0.13 and 0.17 cMpc −1 with 1σ error bars are shown for comparison. These wavenumbers are the most constrained by HERA but still cannot constrain the EDGES compatible models as their power spectra peak at a much higher redshift.
in the case of the EDGES-compatible scenarios, the expected enhancement of the power spectrum mostly falls around the central frequency of the detected global signal (i.e. 78 MHz, corresponding to z ∼ 17), and, thus, far outside the HERA band. This is because, in order to create a deep and narrow signal at z ∼ 17 to explain the EDGES detection, a combination of a strong radio background and strong X-ray heating is required (Fialkov & Barkana 2019;Reis et al. 2020), which is a different corner of the astrophysical parameter space than that excluded by HERA. Models with strong X-ray heating tend to be suppressed at low redshifts, and, therefore, the EDGES-compatible models are also compatible with the limits set by HERA.
To illustrate this behavior, we check which of our simulated global signals are broadly compatible with EDGES 22 and mark them as purple lines in Figure 17. The rest of the models are painted in grey. In the bottom panel we plot corresponding power spectra. As we can see from the figure, indeed, the power spectra of the EDGES-compatible models (purple) peak at much higher redshifts than those probed by HERA and have low power in the HERA band. None of the EDGES compatible models (neither with the radio background from galaxies as shown here, nor with the phenomenological synchrotron radio background) are excluded by the HERA data.

DISCUSSION
The HERA Collaboration has recently presented our first power spectrum limits at z = 7.9 and z = 10.4 in H21. In this paper, we have used a suite of theoretical models to quantify the implications of these upper limits for the landscape of IGM and galaxy formation models of the Cosmic Dawn. The most fundamental of our conclusions, about the IGM temperature at z ∼ 8, are illustrated in Figure 18. We have shown that: • The HERA limits at z = 7.9 strongly disfavor otherwise viable models with weak or no IGM heating by astrophysical sources. In the fiducial scenario where the CMB dominates the high-z radio background, the HERA limits combined with galaxy and EoR observations constrain the spin temperature of the neutral IGM to 27 (2.3) K < T S < 630 (640) K at 68% (95%) highest posterior density confidence interval. This corresponds to kinetic temperatures of the neutral gas of 8.9 K < T K < 1.3 × 10 3 K at 68% confidence. We have demonstrated this with flexible galaxy evolution models (using 21cmFAST; Fig. 9), a phenomenological model that directly parameterizes the IGM properties (Fig. 10), and a simple bias model (Fig. 11).
• Using 21cmFAST and 21cmMC in which astrophysical heating pre-reionization is sourced by galactic X-ray emission, we found that the combination of other high-z observations and the recent limits from HERA allows us to constrain the X-ray luminosity per unit star formation rate of Cos- Figure 18. Summary of lower limits on the IGM spin temperature at z = 7.9 enabled by HERA H21 limits. From top to bottom, we show the 95% lower limits derived via density-driven models ( §4), phenomenological reionizationdriven models ( §6), 21cmMC ( §5), and radio background models ( §8). The red region shows temperatures below the adiabatic cooling limit. Black points show our final constraints, while gray points indicate priors, or constraints derived with simplifying assumptions or without HERA, e.g., the top panel adopts two contrasting scenarios, a fully neutral case and a case with uniform 50% ionization. The bottom three panels all indicate the prior range, as well as the limits obtained via independent measurements alone in the 21cmMC case. The bottom panel includes both homogeneous ( §8.5) and inhomogeneous ( §8.4) radio backgrounds -in the latter case, the correspondence between top and bottom x axes assumes a uniform T rad = Tγ(z = 7.9), but we caution that TS/T rad constraints cannot simply be converted to T S constraints.
• Using a different set of simulations (section 8), we showed that HERA's limit (on its own) places a constraint on an early radio background for models in which galaxies have substantial radio emission but weak X-ray heating. We find that the models with L r,ν /SFR > 4 × 10 24 W Hz −1 M −1 yr (calculated at reference frequency ν = 150 MHz) and L X, 0.2−95 keV /SFR < 7.6 × 10 39 erg s −1 M −1 yr are excluded at 95% confidence. The ratio of radio to gas temperatures must satisfy log 10 (T rad /T K ) < 1.1 at z = 8 to be consistent with HERA at 95% confidence.
Considering a phenomenological model where the radio background (in addition to the CMB) has a synchrotron spectrum, we find that an extra radio background of 1.6% of the CMB at 1.42 GHz or stronger and weak heating by X-ray binaries with luminosity L X, 0.2−95 keV SFR < 5.9 × 10 39 erg s −1 M −1 yr, are robustly excluded by the data (at 95 % confidence). For this model we find that the ratio log 10 (T rad /T K ) < 1.7 at z = 8 is required for consistency with HERA at 95% confidence.
The HERA constraints on X-ray and radio efficiencies are complementary to the constraints from the unresolved X-ray background (measured by Chandra) and the low-frequency radio background (detected by ARCADE2/LWA), ruling out otherwise unconstrained parameter space.
• There has been considerable recent interest in exotic cosmological models that might explain the EDGES measurement at z ∼ 18 (Bowman et al. 2018), either through cooling the IGM below the adiabatic limit or imposing an additional radio background. The H21 limits cannot confirm or disprove the EDGES signal, simply because it observes at a lower redshifts. However, the limits do independently rule out a very cold IGM at z = 10.4 (as well as an IGM at or below the adiabatic limit at z = 7.9), so that heating must occur before z ∼ 10 if millicharged dark matter is behind the anomalous EDGES depth. The EDGES signal itself, if taken at face value, already requires such heating at z 15, so cosmological explanations of that signal are compatible with the H21 limits. Future measurements from HERA at lower frequencies will be important for further testing such scenarios. These conclusions, validated by multiple independent approaches, demonstrate that the H21 upper limits provide important astrophysical constraints on early galaxies. This is the first time IGM heating has been required by observations, independent of any assumptions about galaxy formation during this era. Nevertheless, the upper limits are sufficient to show how state-of-the-art galaxy models can offer new insights, as demonstrated by the difference in temperature constraints in "densitydriven" models.
Moreover, our results emphasize how 21-cm measurements are highly complementary to other probes of the reionization era. No other existing probe can constrain the X-ray or radio emissivity of early-galaxy populations as effectively. The requirement of at least minimal heating (or substantially more if radio emission is strong or if dark-matter interactions are at play) translates to a lower limit on radiative heating by star-forming galaxies (conventionally interpreted as X-ray heating, though at very low temperatures Lyman-α heating also plays a role). While the current upper limits do not strongly constrain the expected parameter space of early galaxies, they show the promise of 21-cm interferometers. We also note that the combination of HERA data with complementary data is most powerful when they probe similar redshifts, so that assumptions about redshift evolution can be avoided.
Of course, there are several caveats to this interpretation and ways in which our methods can be improved: (i) they are driven largely by a very small number of measurements (single bandpowers at z = 7.9 and z = 10.4) so are more subject to systematic concerns; (ii) our constraints do not fully capture the correlations between bandpowers; and (iii) our simulationbased analyses also do not treat the angular dependence of the signal properly, as current observations are sensitive only to modes nearly along the line-of-sight, where the power spectrum is most strongly affected by redshiftspace distortions (see Fig. 11). Our assumptions are strictly conservative for density-driven scenarios, but the effect of redshift-space distortions is more complicated during reionization. We expect to improve the match between simulations and the data for future observing campaigns.
The H21 constraints are based on a single observing season with a small number of antennae and the first iteration of HERA's infrastructure. HERA has contin-ued to expand and improve, and later campaigns will offer much higher sensitivity. This analysis provides a foundation for the interpretation of these forthcoming datasets. this research were supported by the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013. MGS acknowledges support from the South African Radio Astronomy Observatory(SARAO) and the South Africa National Research Foundation (Grant No. 84156).

CONTRIBUTIONS OF THE AUTHORS
S.R.F. leads the HERA Theory team and organized this effort. Y.Q. led the 21cmFAST and 21cmMC simulation and analyses, with particular support from B.G. and A.M. ( §5). J.M. led the analysis using phenomenological models ( §6), with particular support from J.B.M. and S.R.F. J.B.M. led the bias analysis ( §4) and the dark-matter interpretation ( §7). S.H. led the analysis described in §8 using models created by A.F., and developed the emulators based on the architecture by S.S. S.G.M. and N.K. developed the derivation of the likelihood ( §3.2). R.B. is in the Builder's list of the 21-cm code used in section 8, I.R. led the development of the radio-background from galaxies model used in Section 8. The HERA Builder's list is included as authors because of the dependence of this work on the combined efforts of collaborators on the various software repositories as well as the necessity of using HERA data to constrain the models used in this paper.

A. FLUCTUATIONS TRACING THE VELOCITY FIELD
In section 4, we assumed that the 21-cm signal traced the matter power spectrum for simplicity. However, in some regimes the 21-cm power spectrum may more closely trace the velocity-induced acoustic oscillations (VAOs) from the dark matter-baryon relative velocities v cb . This regime can be triggered by a strong dependence of galaxy formation on these streaming velocities at small scales Dalal et al. 2010;Naoz et al. 2012), which however are not expected to play a key role at low redshifts (Fialkov et al. 2012b;Muñoz 2019a) (see however Park et al. 2021;Cain et al. 2020 for the effect during the epoch of reionization). Exotic interactions between dark matter and baryons can modulate the gas temperature with the same velocity feature (Muñoz et al. 2015b;. Regardless of whether its origin is astrophysical or exotic, a detection of VAOs would provide us with a standard ruler (Muñoz 2019b).
To model such a scenario, we replace the matter power spectrum in equation (11) with ∆ 2 v cb (k) and constrain the corresponding b v . The quantity ∆ 2 v cb (k) is the power spectrum of δ v cb = 3/2 (v 2 cb /v 2 rms − 1), which is dimensionless, z-independent, and has unit variance (Dalal et al. 2010;Ali-Haïmoud et al. 2014). Using the H21 data we can set the 95% confidence limits b v cb < {49, 175} mK at z = {7.9, 10.4}, respectively. For comparison, a typical VAO contribution to the 21-cm power spectrum found in the simulations of Muñoz (2019a) and  is roughly ∆ 2 21 ≈ 10 mK 2 (which peaks at z ≈ 15 during their X-ray heating epoch). This would convert into a bias b v cb ≈ 5 mK, much below the limits above.
B. EMULATOR FOR EXTRA RADIO BACKGROUND SIMULATIONS Emulators are widely used to rapidly evaluate 21-cm power spectra across the broad astrophysical parameter space (Kern et al. 2017;Schmit & Pritchard 2018;Jennings et al. 2019;Ghara et al. 2020;Mondal et al. 2020). This technique allows us to interpolate existing simulations instead of running a new simulation at every point of the parameter space.
Here we choose to emulate the power spectrum using a neural network (multi layer perceptron) regression, implemented in scikit-learn (Pedregosa et al. 2011).
We train emulators for the two cases of radio backgrounds separately: one data set corresponding to the models with radio background created by galaxies (10700 samples), the other corresponding to the smooth synchrotron radio background (10300 samples). These sets were created by randomly sampling the parameters V c , f * , f X , τ , and f r or A r within the bounds described in section 8.2. To calculate the likelihood of each parameter set we go through a two-step process: First, we emulate the power spectrum assuming no RSDs (emulators based on the corresponding ∼ 10000 models mentioned above), then we use another emulator to account for the effect of RSDs (using another ∼ 2000 simulations calculated separately for each type of radio background), boosting or suppressing the power spectra by a factor between 0.8 to 4 depending on the values of the astrophysical parameters. The final power spectra are used in the likelihood calculation. We also emulate the values of spin, gas, and radiation temperatures (training data ∼ 10000 cases for each type of radio background) to infer constraints on the physical properties of the IGM (e.g. Figure 11 and log 10 f r All even (band 2) and odd (band 1) data points Only lowest data point per band Figure 19. Top: Accuracy of the neural network predictions. The triangle plots show the relative accuracy of predicting the power spectrum at z = 7.93, k = 0.13 cMpc −1 , for all parameters in the test set. The relative emulator error is shown by the color of the points, from yellow (no error) to red (50% error). For comparison we add the black contour lines from Figure 14 corresponding to 68% and 95% parameter confidence intervals, respectively. The histogram insert shows the overall distribution of the error σ rel = ∆ 2 predict − ∆ 2 test / ∆ 2 predict + 1 mK 2 , for z = 7.93, k = 0.13 cMpc −1 (orange histogram) and for all points combined (blue). The orange histogram can be approximated by a normal distribution with scale σ = 0.2 (pink curve). Botton: The left panel shows the covariance matrix of the relative emulator error (σ rel as used above, variance of 0.03 corresponds to 17% error) over both bands and the first few wavenumbers. The emulator error correlation between the different bands is negligible, but the errors within a band can be strongly correlated with a high correlation coefficient r ≈ 0.9. The middle panel shows the overall impact of taking the emulator uncertainty into account by looking at the marginalized fr, fX and τ constraints. We see that assuming no emulator error (orange), 20% (black, as used in results) or 40% (green) does not noticeably affect the parameter constraints as the observational uncertainty dominates the error budget in the relevant region. The right panel shows parameter constraints using a restricted likelihood, based on just a single data point per band (red contours, using k = 0.17 and 0.13 cMpc −1 for band 1 and 2, respectively) compared to the full likelihood (black contours, as used for our scientific results). The constraints with the single data points are just slightly weaker than using the full likelihood which shows that the parameter constraints are mainly driven by the lowest point in each band.