Improved Constraints on the 21 cm EoR Power Spectrum and the X-Ray Heating of the IGM with HERA Phase I Observations

Of the 182 nights during this season of simultaneous construction, commissioning, and observing, a significant fraction of nights was discarded for a variety of reasons. Most of these were hardware failures, including network outages, power outages, too many low- and/or high-power antennas, a briefly broadcasting antenna, broken receivers, and broken X-engines. Some were due to site issues, including high winds, a lightning storm, and excess radio frequency interference (RFI). While all nights have significant narrowband RFI contamination from FM radio, TV broadcasts, and ORBCOMM satellites (see Section 3.2.3), most nights that were completely flagged for excess RFI showed consistent broadband emission contaminating channels typically free of RFI. These cuts first reduced the 182 nights to 104 nights using contemporaneous observing logs and real-time analysis. After inspecting hundreds of jupyter notebooks ³³ summarizing the nightly results of the data analysis pipeline after each key stage (see Section 3), this was reduced to 94 nights. For more details on the precise selection of nights, see Dillon (2021a).

In Figure 2, we show how the nights passing our various data quality checks span the observing season. The breaks in good data (due to a network outage, a correlator malfunction, and a broadcasting antenna) naturally divided the season into four epochs. The data used in H22a were a subset of Epoch 1. While in theory each night could be analyzed independently before binning all of them together at constant local sidereal time (LST), we found it useful to analyze epochs individually, both for systematics mitigation (see Section 3.2.4) and for statistical tests on subsets of the data (see Section 6.3).

Zoom In Zoom Out Reset image size

Because we observed during the Southern summer, most 12 hr "nights" included some data taken before sunset or after sunrise. These were flagged in our analysis. Further, a number of partial nights were flagged, usually due to excess RFI. The majority of these were due to broadband RFI during the first few hours of the night, possibly attributable to construction activity on site. A few other partial nights were flagged due to suspicious nightly calibration solutions, especially in Epoch 3 when the galaxy was rising at the end of each night. More detail on the precise subset of nights flagged is given in Dillon (2021a).

The end result of our data cuts is a set of observations that, when LST-binned together, are significantly deeper than those in H22a, and cover over 21 hr of LST. As we can see in Figure 3, this data set peaks at ∼70 nights' observing around 7 hr of LST, nearly four times deeper than the observations used in H22a. Roughly speaking, this sets an upper bound on the factor by which our limits might improve due purely to the increase in sensitivity, since the noise on P(k) scales inversely with observing time.

Zoom In Zoom Out Reset image size

Antenna selection began with the nightly data quality monitoring system described in H22a. It identified malfunctioning antennas by looking for antennas participating in baselines that were either outliers in total visibility power, or had visibility amplitude spectra significantly different from other baselines measuring the same physical separation on the ground. This procedure ultimately informed the most rigorous identification of malfunctioning antennas described in Storer et al. (2022), which was applied to HERA Phase II data. The results of nightly analysis were synthesized into a set of per-antenna, per-night flags by the HERA commissioning team as part of an internal data release.

Because our metrics for antenna quality are relative ones, we decided to expand and harmonize this list of flagged antenna-nights under the conservative presumption that antennas that misbehave consistently enough are probably also anomalous at some lower level on other nights. If, in any given epoch, an antenna was flagged more than 10% of the days, we flagged it for the whole epoch. If an antenna was completely flagged for more than 60% of the epochs that it appeared in (i.e., three of four epochs for antennas that were observing for the whole season), then it was completely removed.

Antennas passing this first series of cuts were then used for an initial round of redundant-baseline calibration where per-antenna gains and per-unique-baseline visibilities are solved for simultaneously as part of a large overdetermined system of equations (Liu et al. 2010). Antennas outside the southwest sector of HERA's split-core configuration (155, 156, 180, 181, 182, and 183—see Figure 4—as well as two outriggers not pictured) were excluded as well because they would introduce extra tip-tilt degeneracies (Liu et al. 2010; Zheng et al. 2014, 2016; Dillon et al. 2018) and thus complicate a subsequent sky-based absolute calibration (Li et al. 2018; Kern et al. 2020a).

Zoom In Zoom Out Reset image size

As Dillon et al. (2020) described, redundant-baseline calibration can be cast as a χ²-minimization problem where χ² quantifies how consistent deviations from redundancy are with thermal noise (see Section 3.1 and Equation (2)). If one attributes each baseline's contribution to χ² to both of its constituent antennas equally, one can form a χ² per antenna statistic that is sensitive to particularly nonredundant antennas. In our first round of redundant-baseline calibration, antennas that made significant excess contributions to χ² are flagged in 60-integration (i.e., 10.7 minute) chunks, and then calibration is performed again, iteratively, until no outliers remain. These were then likewise harmonized; antennas flagged for nonredundancy more than 20% of a night or 5% of an epoch were flagged for the whole night/epoch.

In Figure 4, we show the per-antenna flagging fraction after all of these cuts. The overall impact is that while the array was growing, the number of antennas included in our analysis remained largely static (see also Figure 2). Likely this set of nightly flags and antenna flags is overly conservative and some good data were thrown out. Due to the extreme dynamic range challenge faced by all of 21 cm cosmology, we adopted the stance that it was far safer to throw out possibly good data than to risk including bad data. Importantly, all of these decisions about data selection were performed without reference to final power spectra and are therefore less likely to introduce experimenter bias. Once the set of good antennas and nights was selected, it was not subsequently changed.

H22a gives a detailed description of the analysis steps that take raw visibilities all the way to power spectra. Here we provide a high-level overview of each step in order to give context for the changes and updates in this analysis. We refer the reader to H22a and its supporting papers for more detail. The steps in our data reduction and systematics pipeline are as follows:

• 1.  
  Redundant-baseline calibration: We begin with direction-independent calibration, wherein our observed visibilities, ${V}_{{ij}}^{\mathrm{obs}}$, are modeled as

  $$\begin{eqnarray}&&{V}_{{ij}}^{\mathrm{obs}}={g}_{i}{g}_{j}^{* }{V}_{{ij}}^{\mathrm{true}}+{n}_{{ij}}.\end{eqnarray} \tag{ 1 }$$

  Here g_i is a complex, time- and frequency-dependent gain associated with the ith antenna, and n_ij is the noise on that visibility. Redundant-baseline calibration assumes that all baselines b _ij with the same physical separation and orientation should observe the same true visibility, and thus solves for both gains and unique visibilities as a large χ²-minimization problem, where

  $$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \sum _{i\lt j}\displaystyle \frac{{\left|{V}_{{ij}}^{\mathrm{obs}}-{g}_{i}{g}_{j}^{* }{V}_{i-j}^{\mathrm{sol}}\right|}^{2}}{{\sigma }_{{ij}}^{2}}.\end{eqnarray} \tag{ 2 }$$

  ${V}_{i-j}^{\mathrm{sol}}$ here is the visibility solution for all redundant baselines with the same physical separation as V_ij . We calibrate by minimizing χ² for every time and frequency independently, using only internal degrees of freedom and without reference to any model of the instrument or the sky. For an exploration of redundant-baseline calibration in the context of HERA, see Dillon et al.(2020).
• 2.  
  Absolute calibration: The internal consistency of redundant baselines alone cannot solve for three important degrees of freedom, namely the overall gain amplitude and two-phase tip/tilt terms. These degeneracies must be solved, per time and frequency, by reference to externally calibrated (or simulated) visibilities (Kern et al. 2020a). In H22a and in this work, we performed absolute calibration using a set of visibilities synthesized from three nights with LSTs spanning this data set (JDs 2458042, 2458116, and 2458207). These were calibrated on three separate fields using the MWA GLEAM catalog (Hurley-Walker et al. 2017) and CASA (McMullin et al. 2007), as described in Section 3.3 of H22a.
• 3.  
  RFI flagging: RFI is identified and flagged using an iterative outlier detection algorithm described in H22a. Essentially, several sets of waterfalls (visibilities, gains solutions, etc.) are converted into time- and frequency-dependent measures of "outlierness" expressed as a z-score or modified z-score. This is done by looking at a 17 × 17 pixel region centered on each pixel of the waterfall and measuring how consistent each is with its neighbors in time and frequency. After averaging together z-scores of baselines or antennas to improve the signal-to-noise ratio (S/N), 5σ outliers and 2σ outliers neighboring 5σ outliers are flagged. This is done independently for each set of waterfalls, and the flags are then combined. We start RFI flagging with median filters and modified z-scores to reduce the effect of really bright RFI, then use those flags as a prior on a second round using mean filters and standard z-scores. Finally, we examine these statistics over the whole day, looking for whole channels or whole integrations that are 7σ outliers and their 3σ outlier neighbors. The same flags are applied to all baselines.
• 4.  
  Gain smoothing: After flagging, all gains solutions are smoothed to mitigate the effect of noise and calibration errors, taking as a prior that the gains should be stable and relatively smooth in frequency. This smoothing is performed with a CLEAN-like deconvolution algorithm (Högbom 1974; Parsons & Backer 2009), filtering gains in 2D Fourier space on a 6 hr timescale and a 10 MHz frequency scale (or equivalently, 100 ns delay scale). These are the scales on which we have evidence for intrinsic gain variation in time (Dillon et al. 2020) and frequency (Kern et al. 2020a). For more implementation details, see H22a.
• 5.  
  LST-binning: Having calibrated and flagged each night's data, we then coherently average nights together on a fixed LST grid. This 21.4 s cadence grid—double the integration time of raw visibilities—is created by assigning each observation to the nearest gridded time and then rephasing that visibility to account for sky rotation due to the difference in LST. An additional round of flagging is performed on a per-LST, per-frequency, and per-baseline basis, looking for 5σ outliers ³⁴ in modified z-score among the list of rephased visibilities to be averaged together. This cut is designed to identify low-level residual RFI and calibration failures; it is highly unlikely to flag outliers due to noise.
• 6.  
  Hand-flagging: After LST-binning, a final set of additional flags are added by manually examining high-pass delay-filtered residuals. This filtering was performed on using an iterative delay CLEAN to remove power below the 2000 ns scale in order to highlight spectrally compact features. Clear outlier channels and/or times that are consistent across baselines are flagged upon visual inspection. The same additional flags are applied to all baselines.
• 7.  
  Inpainting flagged channels: When using fast Fourier transforms (FFTs) to form power spectra, flagged channels introduce discontinuities that leak foreground power to high delays. To avoid this, we use the same delay-based CLEAN algorithm to low-pass filter the data on 2000 ns scales and inpaint the flagged channels with the filtered result. Entirely flagged times are not inpainted. For more details and a demonstration of this procedure, see H22a.
• 8.  
  Cable reflection calibration: When a signal bounces off of both ends of a cable before being transmitted, the result is a copy of the signal at a fixed delay associated with twice the light travel time along the cable. Kern et al. (2020b) showed how the 20 m and 150 m cables in the signal chain produce these reflections, which can be represented as complex gain terms. While these gains are in principle calibratable with redundant-baseline calibration, our gain smoothing procedure completely removes them. Thus, following the procedure outlined in Kern et al. (2019), we iteratively model and calibrate out reflection and subreflection terms using autocorrelations, which have much higher S/Ns than cross-correlations. Since cable reflections are stable over many nights, this is done after LST-binning.
• 9.  
  Crosstalk subtraction: Kern et al. (2020b) also demonstrated the pernicious impact of crosstalk systematics across a large range of delay modes in HERA Phase I data. It was hypothesized that this was due to over-the-air crosstalk that led to autocorrelations leaking to cross-correlations at high delays. In the Appendix, we show how this effect's delay and amplitude structure can be explained by an emitter in the signal chain. Because autocorrelations are nonfringing and thus quite stable in time, the effect can be mitigated by modeling each baseline's excess power near zero fringe-rate. Kern et al. (2019) did this using singular value decomposition (SVD) to find the delay and time modes affected and used Gaussian process regression to limit the range of fringe-rates modeled and subtracted so as to avoid EoR signal loss. That fringe-rate range is symmetric about zero and limited by the east–west projected baseline length ∣b_E−W∣ such that

  $$\begin{eqnarray}&&\left|{f}_{\max }\right|=0.024\,\mathrm{mHz}\left(\displaystyle \frac{| {b}_{{\rm{E}}-{\rm{W}}}| }{1\,{\rm{m}}}\right)-0.28\,\mathrm{mHz}.\end{eqnarray} \tag{ 3 }$$

  The signal loss due to crosstalk subtraction was calculated in Aguirre et al. (2022) and corrected for in H22a; we repeat the calculation in Section 4.3 and find similar results. Crosstalk subtraction is not applied to baselines with projected east–west distances less than 14 m, since the zero fringe-rate mode overlaps with the fringe-rates associated with the main lobe of the primary beam (Parsons et al. 2016).
• 10.  
  Coherent time averaging: Following H22a, we next coherently average visibilities from the 21.4 s cadence after LST-binning down to a 214 s cadence, using the same rephasing procedure described above to account for sky rotation. This timescale was chosen to keep signal loss at the ∼1% level (Aguirre et al. 2022); we repeat that calculation for this data set in Section 4.3 and find consistent results.
• 11.  
  Forming pseudo-Stokes I visibilities: Before forming power spectra, we construct pseudo-Stokes visibilities. As H22a showed, this limits the leakage of Faraday-rotated foregrounds into high delays, though Q → I leakage from primary beam asymmetry is still possible (Moore et al. 2013; Asad et al. 2016; Kohn et al. 2016; Nunhokee et al. 2017; Asad et al. 2018). For the pseudo-I channel, this step consists of simply averaging together the EE- and NN-polarized visibilities for the same baseline, where E and N denote the east- and north-aligned linearly polarized feeds, respectively (Kohn et al. 2019).
• 12.  
  Power spectrum estimation: We compute power spectra using the delay approximation, in which we substitute a Fourier transform along the frequency axis of a visibility (i.e., a delay transform) for a line-of-sight Fourier transform (Parsons et al. 2012a). This strategy avoids mapmaking entirely (Dillon et al. 2013, 2015b; Xu et al. 2022). We thus approximate τ and u (the magnitude of the baseline in units of wavelengths) as mapping linearly to line-of-sight Fourier modes, k_∥, and transverse Fourier modes, k_⊥, respectively. The power spectra are estimated by taking the real part of

  $$\begin{eqnarray}&&\hat{P}({k}_{\perp },{k}_{\parallel })=\displaystyle \frac{{X}^{2}Y}{{{\rm{\Omega }}}_{{pp}}B}{\widetilde{V}}_{1}(u,\tau ){\widetilde{V}}_{2}^{* }(u,\tau ),\end{eqnarray} \tag{ 4 }$$

  where $\widetilde{V}$ is a Fourier transformed visibility in frequency. Ω_pp is the full-sky integral of the squared primary beam response—we use the beam simulated in Fagnoni et al. (2021)—and B is the bandwidth. X and Y are scalars mapping angles and frequency to cosmological distances, defined via k_∥ = 2π τ/X, k_⊥ = 2π∣ b ∣/(Y λ), with X = D(z) and Y = c(1 + z)²[H(z)]⁻¹/(1420 MHz), where H(z) is the Hubble parameter, and D(z) is the transverse comoving distance (Hogg 1999). ³⁵ $\hat{P}({k}_{\perp },{k}_{\parallel })$ here is in units of mK² h⁻³ Mpc³, though it conventional to report the "dimensionless" power spectrum (which has units of mK²):

  $$\begin{eqnarray}&&{{\rm{\Delta }}}^{2}(k)\equiv P(k)\displaystyle \frac{{k}^{3}}{2{\pi }^{2}}.\end{eqnarray} \tag{ 5 }$$

  H22a shows how this can be recast into the language of quadratic estimators (Tegmark 1997; Liu & Tegmark 2011; Dillon et al. 2013; Trott et al. 2016). In that formalism, we use a diagonal normalization matrix M (i.e., no decorrelation of bandpower uncertainties). In lieu of any inverse covariance weighting, we use a Blackman-Harris taper to prevent foreground leakage into the EoR window. When computing Equation (4), we cross-multiply Fourier transformed visibilities from alternating 214 s blocks of time (i.e., ${\widetilde{V}}_{1}$ and ${\widetilde{V}}_{2}$), using all pairs of baselines in a redundant-baseline group. In the (quite accurate) approximation that visibilities interleaved in this way have uncorrelated noise, this produces an estimate of the power spectrum free of noise bias.
• 13.  
  Error estimation: As in H22a, we employ the noise estimation formalism of Tan et al. (2021). The noise power spectrum is given by

  $$\begin{eqnarray}&&{P}_{N}=\displaystyle \frac{{X}^{2}Y{{\rm{\Omega }}}_{\mathrm{eff}}{T}_{\mathrm{sys}}^{2}}{{t}_{\mathrm{int}}{N}_{\mathrm{coherent}}\sqrt{2{N}_{\mathrm{incoherent}}}},\end{eqnarray} \tag{ 6 }$$

  where T_sys is the system temperature, t_int is the integration time, and N_coherent and N_incoherent are the numbers of integrations averaged together coherently or incoherently—i.e., averaged as visibilities with phase information or averaged as power spectra. Ω_eff is the effective beam area, defined as ${{\rm{\Omega }}}_{\mathrm{eff}}\equiv {{\rm{\Omega }}}_{p}^{2}/{{\rm{\Omega }}}_{{pp}}$ in Appendix B of Parsons et al. (2014). We use the inverse square of the noise power spectrum to perform inverse-variance-weighted averaging of the power spectra. For reporting final error bars on power spectra, we use an unbiased estimator of the noise and signal-noise cross terms developed by Tan et al. (2021),

  $$\begin{eqnarray}\begin{array}{rcl}{\hat{P}}_{\mathrm{SN}} & = & {P}_{N}\left(\sqrt{\sqrt{2}(\hat{P}/{P}_{N})+1}\right.\\ & & \left.\,\left[\sqrt{1/\sqrt{\pi }+1}-1\right]\right).\end{array}\end{eqnarray} \tag{ 7 }$$
• 14.  
  Incoherent power spectrum averaging: Here and in H22a, power spectra are averaged incoherently over several axes to produce the final limits. First, all baseline-pairs within a redundant-baseline group are averaged, ignoring auto-baseline pairs (power spectra formed from the same pair of antennas at interleaved times). This preserves most of the sensitivity of coherently averaging visibilities within a redundant group before forming power spectra, while excluding the pairs most likely to exhibit correlated systematics. Then power spectra are averaged incoherently in several disjoint LST ranges, which we call "fields" since they correspond to different parts of the sky at zenith. Power spectra are estimated independently for the LST ranges in the separate fields, and we perform no further averaging in H22a or this analysis when reporting power spectrum upper limits. Finally, we perform a binning to spherical $k=\sqrt{{k}_{\parallel }^{2}+{k}_{\perp }^{2}}$, excising baselines based on their length and certain sets of delay modes based on their proximity to the "horizon wedge," which is set by the light travel time along the baseline

  $$\begin{eqnarray}&&{\tau }_{\mathrm{wedge}}=| {\boldsymbol{b}}| /c.\end{eqnarray} \tag{ 8 }$$

  P_N and ${\hat{P}}_{\mathrm{SN}}$ are propagated through each of these averaging steps.

For more details on the implementation of these techniques, we refer the interested reader to H22a and its supporting papers.

With the full context of the analysis pipeline established above, we now detail the ways it has changed since H22a. While most are relatively minor tweaks (Dillon 2021a), we detail them here for completeness and reproducibility.

3.2.1. Updates to Redundant-baseline Calibration

3.2.2. Updates to Absolute Calibration

3.2.3. Updates to RFI Excision

While the fundamental algorithm for RFI excision remains unchanged, we made two updates to how it was performed on a nightly basis. The first is related to how the analysis dealt with data file boundaries. Previously, every data file was analyzed in parallel for RFI. Since the outlier identification algorithm relies on neighboring times and frequencies, this became less reliable near file boundaries where there are roughly half as many neighbors to compare to, and likely led to some of the ∼10 minutes periodicity we saw in the flags in H22a. In this analysis, we parallelized the pipeline in overlapping time chunks, so that every point was compared to exactly the same number of neighbors—except at the beginning and end of the night and at the edges of the band.

Second, we modified the set of data products used in searching for outliers. In both analyses we used raw visibilities (albeit only in the mean filter round); gains from both redundant-baseline calibration and absolute calibration; and χ² from both calibration steps. Based on experiments we performed on which data products were providing unique information and not leading to over-flagging, we removed a global cut on outliers in χ²—which likely led to the over-flagging around Fornax A in H22a—and added uncalibrated autocorrelations for their high S/N and computational tractability compared to the full set of visibilities. The result is still a very expansive set of flags and likely contains a significant number of false-positive identifications of RFI, especially around ORBCOMM at 137 MHz and the clock line at 150 MHz (see, e.g., Figures 5 and 8). Given the extreme dynamic range requirements of 21 cm cosmology, it is far safer to over-flag than under-flag (Kerrigan et al. 2019).

Zoom In Zoom Out Reset image size

In that spirit, we also revisited how a final set of by-hand flags were identified. In H22a, this flagging was performed by examining a handful of key baselines after high-pass filtering to remove structure below the 2000 ns scale. In this work, the four epochs were first combined together without any inpainting, reflection calibration, or crosstalk subtraction. After performing the same per-baseline high-pass filter on every baseline, their amplitudes were averaged together, inverse-variance weighted by noise, and then corrected for the noise bias.

The result, shown in the left-hand panel of Figure 5, highlights residual frequency structure. Much of this has low N_samples and/or borders on previously identified RFI, indicating its origin as low-level, inconsistently flagged RFI. As the right-hand panel shows, outlier channels and integrations were first identified by averaging along each axis. Next, individual areas of concern were flagged by hand, by converting the waterfall to a bitmap image and individually marking times and frequencies in Adobe Photoshop. An effort was made to flag coherent rectangular regions near previously identified flags to avoid cherry-picking, though fundamentally this step involved a series of subjective judgment calls. Once the final flagging waterfall was developed, it was not revisited after estimating power spectra in order to avoid experimenter bias.

3.2.4. Updates to Reflection Systematics and Crosstalk Mitigation

In H22a, all of the steps in Section 3.1 after LST-binning were performed on the full-sensitivity, 18 nights data set. At first, we attempted performing the same analysis with all 94 nights binned together but found that the level of residual crosstalk had substantially higher S/N than anything seen in H22a. Some baselines were particularly bad, exhibiting a delay- and time-averaged S/N greater than 10 in the affected delay range. Upon examining the pre-subtraction waterfalls of the baselines where crosstalk subtraction performed the worst, we found clear evidence for temporal structure in the delays contaminated by crosstalk. In Figure 6, we show one such baseline. Plotting the real part of its waterfall in delay space shows clear temporal structure at certain delays, including some where it flips from positive to negative and vice versa. This is potentially disastrous for the crosstalk subtraction technique of Kern et al. (2019), which relies on stability in delay space.

Zoom In Zoom Out Reset image size

Furthermore, there appears to be a correlation between discontinuities in N_samples (right-hand panel of Figure 6) and discontinuities in the delay structure of ${\tilde{V}}_{{ij}}$. Since the former are largely attributable to epoch boundaries, which affect how often each LST is observed, we hypothesized that the changing and growing array was affecting the precise structure of the crosstalk. This ultimately pointed us toward a new understanding of the physical origin of the effect, namely that all antennas' signals were being broadcast from one point on the west side of the array, likely the refrigerated enclosures that contained the analog receivers. We discuss this model and the evidence for it in detail in the Appendix.

The upshot of this result is twofold. First, it confirms that the model of Kern et al. (2019), of autocorrelations leaking into cross-correlations, is correct. Second, it implies that as long as the array is stable, the effect should be stable in LST-binned data as well. We therefore decided to perform inpainting, cable reflection calibration, and crosstalk subtraction on a per-epoch basis before binning together the four epochs. This proved a substantially better approach—as Figure 7 shows—and got us much closer to consistency with noise after crosstalk subtraction.

Zoom In Zoom Out Reset image size

A few minor tweaks were also made to cable reflection calibration and crosstalk subtraction. The total number of terms used to fit the cable reflections was increased from 28 terms between 75 and 1500 ns to 35 terms between 75 and 2500 ns, to both better model the 20 m cables and to be able to model a few extra long cables whose reflection timescale was larger than 1500 ns but which did not appear in the H22a data set. The SVD used in crosstalk subtraction was previously limited to 30 delay modes; we increased it to 50 based on experiments where it made the crosstalk residuals a bit more consistent with noise. Finally, we revised how weights were applied before computing the SVD. First, we weighted each time by the frequency-averaged number of samples. H22a used an unweighted SVD, but as Figure 3 shows, the approximation of weights that are constant in time breaks down when considering such a large and discontinuous data set. Second, we also set the weight in the SVD to zero from 15.3 to 21 hr to prevent bright galactic emission at the edge of the data set from introducing temporal structure that caused the crosstalk subtraction procedure to perform worse for the most sensitively measured LSTs.

Just as with the data reduction pipeline, we sought to apply the same power spectrum analysis procedures and choices as were used in H22a. However, differences in LSTs observed, flagging, and systematics removal motivated slightly different approaches to how the data should be reduced and cut. These decisions were made without reference to the final power spectra in order to minimize experimenter bias.

3.3.1. Picking Bands and Fields

The two frequency bands analyzed in H22a were Blackman-Harris tapered ranges from 117.09–132.62 MHz and 150.29–167.77 MHz. These were motivated by the two contiguous regions of low flag occupancy (see Figure 12 in H22a). We reproduce that same plot in Figure 8 and, following the same logic (albeit with somewhat different flags), pick Band 1 and Band 2 to range from 117.19–133.11 MHz and 152.25–167.97 MHz, respectively. The bands still center on approximately the same redshifts: z = 10.4 and z = 7.9.

Zoom In Zoom Out Reset image size

Because we observed a larger range of LSTs than in H22a, we need to define new fields (i.e., new LST ranges) in which to independently estimate power spectra. In order to motivate the choice of fields without reference the final power spectra, we looked at two other statistics, which we plot in Figure 9. These two statistics are the inverse-variance-weighted, baseline-averaged delay-zero power, P(τ = 0)—a proxy for foreground power—and P_N , which is flat in delay and tells us about both foreground power and total observation time. We used these to define a total of five fields: A (21.5–0.0 hr), B (0.75–2.75 hr), C (4.0–6.25 hr), D (6.25–9.25 hr), and E (9.25–14.75 hr). ³⁶ We restricted our field boundaries to quarter-hour increments to avoid cherry-picking integrations.

Zoom In Zoom Out Reset image size

The rationale for defining the fields is as follows. We wanted fields B, C, and D, to correspond reasonably well to fields 1, 2, and 3 in H22a, so to cover the new LST ranges, we added fields A and E. Field A was set by the flagging gaps at either end, intentionally avoiding the last integration before the flagging gap between 0 and 1 hr due to the potential for signal loss from crosstalk subtraction (Aguirre et al. 2022). Likewise, field B was defined to exclude the first integration after the gap and keep Fornax A in the main lobe no brighter than its brightest point in the first sidelobe around 2 hr. Field C was defined to start at a roughly symmetrical place to where Field B ended with respect to Fornax A and to include the range of maximum sensitivity from roughly 4–6 hr. Thus, the upper field boundary was set by the sidelobe of the Galactic plane at 6.25 hr. The boundary between fields D and E was set to include the roughly symmetrical sidelobe at ∼9 hr within Field D, keeping the Galactic plane contained to a single field. Field E ends a bit before where the crosstalk subtraction gets zero weight in the SVD. Once these field definitions were established, they were not allowed to change.

3.3.2. Excluding Baseline Pairs with Substantial Residual Crosstalk

Despite subtracting crosstalk on a per-epoch basis, we still found strong evidence for residual crosstalk on certain baselines. In the baseline examined in Figure 7, we can see clear residual power as an excess S/N in the delay range of 800 < ∣τ∣ < 1500 ns (right-hand panel). This is most prominent near the Galactic center, which got zero weight in the SVD, and near Fornax A, but there appears to be a slight excess at other LSTs as well. To quantify this, we averaged ∣S/N∣ over that delay range and over the three most sensitive fields, B, C, and D. This average was performed separately for positive and negative delays, since we now know that those two signals have independent origins (see the Appendix).

We computed this averaged ∣S/N∣ for every auto-baseline pair (i.e., power spectra formed from the same baseline at interleaved times, rather than power spectra formed from different but redundant baselines) and plotted a histogram of it for each band in Figure 10, treating positive and negative delays as independent samples. Compared to an equivalently sized delay range from 2800 < ∣τ∣ < 3500 ns (solid black), we see evidence for a mild excess on most baselines. This is perhaps not too surprising; the crosstalk subtraction algorithm of Kern et al. (2019) attempts to model and subtract the crosstalk down to the noise—in our case, the noise in a single epoch. To the extent that crosstalk remains correlated from epoch to epoch, integrating down should reveal more crosstalk. That said, there is a tail of outliers in Figure 10 that motivated us to perform a cut at ∣S/N∣ > 1.5. The cut was performed separately for positive and negative delays, so some baseline pairs are "half-flagged." The vast majority of auto-baseline pairs (>95%) were kept. More baseline-pairs were cut at negative delays than positive delays because the antenna ordering means that negative delays were more often associated with antennas nearer the crosstalk source (see the Appendix).

Zoom In Zoom Out Reset image size

We also computed the ∣S/N∣ for cross-baseline pairs as well and found that they were highly correlated with the ∣S/N∣ of the two corresponding auto-baseline pairs. However, we decided to more conservatively use only the auto-baseline pairs—which are not included in the final power spectra—for our cut. Any cross-baseline pair with one baseline participating in a flagged auto-baseline pair (and delay sign) was flagged. This cut is the most surgical and perhaps more worrisome analysis change from H22a, in the terms of removing individual power spectra before averaging. However, by only looking at high ∣τ∣—well above the corresponding k values where we set our tightest upper limits—and by using only the auto-baseline-pairs, we insulate ourselves from the risk of cherry-picking and signal loss.

One other key change from H22a is the shortest east–west projected baseline length allowed to be included in our final spherical power spectra. Even after averaging cross-baseline pairs within redundant groups, we still see a substantial uptick in ∣S/N∣ in the crosstalk delay range for baselines with 14.6 m projected east–west baselines, the shortest baselines used in H22a (see Figure 11). This makes sense physically—these baselines have their main lobe closest to zero fringe-rate, where the crosstalk subtraction algorithm only operates on a very narrow range of fringe-rates (Equation (3)) for fear of signal loss. While the crosstalk is centered at 0 mHz, it has some width in fringe-rate space. We expect, therefore, that these baselines should be the first to show residual crosstalk as we integrate down. It was also likely true in H22a; the lower noise level in these data simply makes the systematics clearer.

Zoom In Zoom Out Reset image size

We decided to conservatively increase that cut to 15 m, throwing out several redundant-baseline groups, including the most sensitive single baseline group used in H22a: the single-unit 14.6 m east–west baseline. To keep this baseline without simply accepting excess systematics, we would have had to find a way to more aggressively filter crosstalk, which would then have necessitated a more thorough and precise quantification of baseline-dependent signal loss using end-to-end simulations. Since our aim in this work is to apply the analysis of H22a as directly as possible, we defer such an investigation to future work.

3.3.3. Changes to k Cuts and Bins

3.3.4. More Precisely Calculating Power Spectrum Window Functions

The expectation value of the estimated power spectrum for a given baseline and delay, $\hat{P}({\boldsymbol{u}},\tau )$, is actually a weighted sum of the neighboring true bandpowers P( k ). These weights are usually referred to as the window functions W, defined through

$$\begin{eqnarray}&&\hat{P}(u,\tau )\propto \int {{\rm{d}}}^{3}{\boldsymbol{k}}\,P({\boldsymbol{k}})W({\boldsymbol{k}};{\boldsymbol{u}},\tau ).\end{eqnarray} \tag{ 9 }$$

In H22a, the horizontal error bars on the spherical power spectrum were evaluated with the same assumptions made to analyze the data, i.e., the delay approximation, in which the delays and line-of-sight Fourier modes are treated interchangeably. This leads to underestimating the tails of the window functions, and in particular the foreground power leaking into the EoR window at low k_∥ (Liu et al. 2014a). In this work, we estimate the exact window functions by lifting this approximation in the derivative of the covariance with respect to each bandpower, and hence obtain an accurate description of the mapping between instrumental space ( u , τ) and cosmological space (k_⊥, k_∥). In doing so, we can account for the delay approximation when comparing theory to data because we now know exactly which k_⊥ and k_∥ modes contribute to a given bandpower and in what proportion. Note that, in order to account for the frequency dependence of the HERA primary beam, we have used the simulations introduced in Fagnoni et al. (2021). For details on the derivation of the window functions, as well as a complete illustration of their importance in the analysis of low-frequency radio data in general, and of the HERA data in particular, we refer the interested reader to Gorce et al. (2023). We discuss the impact of the improved window function calculation in Section 5.2.

3.3.5. Quantifying Decoherence due to Nonredundancy

Because we average together power spectra of pairs of different baselines within a redundant group, we must quantify the effect of nonredundancy on the power spectrum. We know HERA's putatively redundant baselines are not quite redundant (Carilli et al. 2020; Dillon et al. 2020), so we should expect some level of decoherence when cross-multiplying baselines that see a slightly different beam-weighted sky. Following H22a, we compare incoherent power spectrum averages—which are decoherence-free by construction because they only use auto-baseline pairs—to forming power spectra from visibilities coherently averaged within a redundant-baseline group. This yields a metric for decoherence given by

$$\begin{eqnarray}&&{\rm{\Delta }}\kappa (t,\tau )\equiv \displaystyle \frac{{P}_{\mathrm{coherent}}(t,\tau )-{P}_{\mathrm{incoherent}}(t,\tau )}{\langle {P}_{\mathrm{incoherent}}(t,\tau )\rangle },\end{eqnarray} \tag{ 10 }$$

where the angle brackets indicate a rolling time average over 1 hr timescales to ameliorate the effects of nulls in power for certain baselines (as was done in H22a).

In particular, we examine Δκ(t, τ = 0), which we take as a foreground-dominated (and thus high S/N) metric of decoherence of sky signal in the primary beam. In Figure 12, we show the histogram of decoherence levels at zero delay, using different LSTs and unique baselines as samples of decoherence. H22a performed this analysis over a short range of LSTs around the galactic plane crossing (7.2–8.3 hr) and found a ∼1% signal loss due to decoherence. We measure Δκ by taking a median over a wider range of LSTs—all five fields—in part to account for LST- or JD-dependent gain errors and the possibility that the array became more or less redundant over the season. Assuming that the signal loss due to nonredundancy should relatively be stable in LST, and with the knowledge that nulls in power on certain baselines can create spurious temporal outliers in Δκ that always create extra apparent decoherence, even with the 1 hr smoothing, we take the median of this histogram as our signal loss estimate. This yields a slightly larger estimate of the loss: 1.9% for Band 1 and 2.4% for Band 2. This will be accounted for in our final power spectrum upper limits, as we will discuss presently.

Zoom In Zoom Out Reset image size

3.3.6. Correcting for Potential Biases and Signal Loss

Our primary method for validating the analysis pipeline presented in Section 3 is via an end-to-end simulation, wherein we generate realistic visibility simulations of foregrounds, noise, and a boosted EoR-like signal that should be easily detectable given our sensitivity. This allows us to holistically evaluate the performance of the analysis and identify any unknown sources of bias. One could also approach the same problem by injecting a signal of known amplitude (larger than the real EoR) into the visibilities and analyzing that data in the same way. However, this technique requires high confidence in one's calibration solutions (injected visibilities must be "uncalibrated" before injection), and it is difficult to disentangle residual systematics form power spectrum biases.

Our sky simulations consist of three unpolarized components: diffuse Galactic emission, a point-source catalog, and an EoR analog created as a Gaussian random field drawn from a known power spectrum. We present each of these here and discuss how they differ from the corresponding components in Aguirre et al. (2022). For the diffuse Galactic emission, we use the Global Sky Model (GSM; de Oliveira-Costa et al. 2008), computed at every frequency we measure using pygdsm. ³⁷ The result is an HEALPix map (Górski et al. 2005) with resolution N_side = 512. We smooth the map with a 1° Gaussian kernel and degrade it to N_side = 128 to speed up our simulations, since that provides sufficient resolution to cover the angular scales HERA is sensitive to. Each HEALPix pixel is treated as a single point source at its center using the pixel area for normalization.

Our model of point-source emission differs significantly from that of Aguirre et al. (2022) in our handling of spatial gaps in the MWA GLEAM catalog (Hurley-Walker et al. 2017) due to the Galactic plane. Since H22a used observations spanning a much smaller range of LSTs than those used in this work, it was decided to slightly further restrict the range of LSTs simulated to 1.5–7 hr in order to avoid those gaps. In this work, we wanted to be able to simulate the full 24 hr of LST, so we developed a GLEAM-analog with simulated point sources. We created 14,073,688 random synthetic sources with uniformly random positions across the sky. Their fluxes were drawn to match the GLEAM source count distribution in Franzen et al. (2019) from 0.001–87 Jy. Their spectral indices were drawn from a Gaussian distribution with mean −0.79 and standard deviation 0.05, again following Franzen et al. (2019; see also, e.g., Carroll et al. 2016; Offringa et al. 2016). The 200,000 brightest sources were simulated at their random positions, the remainder (those below 0.1975 Jy) were treated as confusion noise. For the sake of computational expediency, these were added to the GSM at the nearest HEALPix grid point. Above 87 Jy, we added several real sources in their true positions. This includes GLEAM J215706-694117, GLEAM J043704+294009, GLEAM J122906+020251, GLEAM J172031-005845, and the "A-Team" sources reported in Table 2 of Hurley-Walker et al. (2017). We also include Fornax A, which we model as a single 750 Jy source with a spectral index of −0.825 at an R.A. of 3^h22^m42^s and a decl. of $-37^\circ 12^{\prime} 2^{\prime\prime}$ (McKinley et al. 2015).

Finally, our boosted EoR analog is created with the same simulator ³⁸ used in Aguirre et al. (2022). We used a higher amplitude mock EoR ³⁹ and a slightly different power spectrum slope: P(k) ∝ k^−2.7. That EoR simulation was binned to the same HEALPix grid as the diffuse foregrounds, where again each pixel center is treated as a point source.

To actually simulate visibilities, we use vis_cpu (Bull 2021), ⁴⁰ a fast visibility simulator validated against pyuvsim (Lanman et al. 2019), ⁴¹ a reference simulator designed for accuracy (R. Pascua et al. 2023, in preparation). The simulator used in Aguirre et al. (2022), RIMEz, ⁴² calculates visibilities in spherical harmonic space. vis_cpu takes a much simpler approach. It calculates a per-antenna visibility factor for each pointlike sky component—essentially the square root of the source flux with a phase factor that depends on frequency and antenna position—and multiplies them by a Jones matrix. These are then cross-multiplied to form visibilities and summed over sources. We use the primary beam calculated in Fagnoni et al. (2021) and interpolated in azimuth and zenith angle. We simulate the full 24 hr of LST at a 5 s cadence for each unique baseline and frequency observed by HERA.

Our end-to-end validation began with producing two sets of 94 nights of simulated data, one with the mock boosted EoR and one without. Both included foregrounds, noise, and systematics. For each night, we interpolated the original set of simulated visibilities using a cubic spline onto each night's 10.7 s cadence, since the LST grid of the data varied from night to night. Unlike in Aguirre et al. (2022), where we used only a subset of the antennas, in this work we inflated the redundant baselines to produce data files with all baselines not flagged in the real data. The final simulated data set completely matched the real data in terms of baselines, times, and frequencies. No nonredundancy due to antenna-to-antenna variation of beams or positions offsets (Orosz et al. 2019; Dillon et al. 2020; Choudhuri et al. 2021) was simulated.

We then "uncalibrate" the data by applying per-antenna complex bandpasses and cable reflection terms and then add noise to the visibilities, both steps performed exactly the same way as in Aguirre et al. (2022). Finally, our simulations have per-baseline crosstalk added to them. Again, this is done with nearly the same procedure as in the prior validation work. For each baseline, we model the crosstalk as a series of copies of the autocorrelation of each antenna in the baseline, each multiplied by a complex delay term and an amplitude that decreases exponentially with delay by a factor of 100 from the peak delay out to 2000 ns. In this work, the amplitudes and delays of the crosstalk peak, and how they depend on position in the array, are motivated by the new physical model of the crosstalk (see the Appendix) in which we attribute crosstalk to an emitter on the west side of the array. (In our model we use α = −2.29 and τ_offset = 0 ns; see Equation (A4) for definitions.) The crosstalk structure is allowed to change per epoch, but not per night. The amplitudes are on the low end of what is observed in real data in order to avoid crosstalk effects on certain baselines becoming crosstalk-dominated, since the visibility simulation itself is also somewhat underpowered on long baselines relative to the real sky. We do not explicitly model any of the multipath effects hypothesized to be responsible for the breadth of the crosstalk spectrum, relying instead on that series of peaks as in Aguirre et al. (2022). All of the techniques for visibility corruption, along with an interface to vis_cpu, are packaged into hera_sim. ⁴³

With our procedure for turning visibility simulations into full nights of "uncalibrated," systematics-corrupted data matching the real data, we then apply our analysis pipeline almost exactly as described in Section 3. We perform redundant-baseline calibration and absolute calibration, then flag, then smooth our gain solutions. We next bin together individual epochs and perform inpainting, cable reflection calibration, and crosstalk subtraction. After binning together the four epochs, we form pseudo-Stokes visibilities, time average, and estimate power spectra. We run the entire end-to-end pipeline twice—once without the mock EoR and once with it—including the various power spectrum cuts described in Section 3.3. In order to faithfully validate the analysis pipeline, the same software packages—pyuvdata ⁴⁴ (Hazelton et al. 2017), hera_cal, ⁴⁵ hera_qm, ⁴⁶ and hera_pspec ⁴⁷ —with the same git hashes were used to run the end-to-end validation. ⁴⁸

The one major step that was performed differently was RFI flagging. Following Aguirre et al. (2022), we do not inject RFI and then attempt to detect it. Rather, we simply cross-apply the real data's flags to the simulated data at the same step in the pipeline. Since the times, frequencies, and baselines match perfectly, this step is very straightforward. Likewise, we take the final flagging mask generated from a manual inspection of the data (see Section 3.2.3) and apply it to each epoch of simulated data at precisely the same point in the real pipeline that those flags are applied to real data.

Before we discuss the final power spectra, we can now evaluate whether the fix offered in Section 3.2.2 eliminates the bias in absolute calibration discovered in Aguirre et al. (2022) and corrected for in H22a. In Figure 13, we compare the calibration solutions produced by the pipeline for Epoch 1 (the longest epoch) with the known simulated gains, averaging over antennas, times, and nights. By comparison to Figure 11 of Aguirre et al. (2022), we can see that the absolute calibration bias is largely eliminated and that the phases are still well-recovered across the band. The gain errors are largest at the band edges, due to edge effects of the low-pass filter used in gain smoothing.

Zoom In Zoom Out Reset image size

Interestingly, we actually now see a slight negative bias for most frequencies (0.6% for Band 1, 0.9% for Band 2), which could lead to overestimating power spectra and error bars by a few percent, since gain errors impact power spectra quartically. One possible origin of the effect is a rare failure mode of the absolute calibration bias fix described in Section 3.2.2. When fitting for a single overall amplitude, solutions are biased but quite stable in time. When fitting for a complex number and then taking the absolute value of that, we have seen rare instances where the data and reference visibilities are so far apart that the gains are driven to 0 in order to minimize χ². In real data, these sorts of collapses are easily identified as discontinuities and flagged as RFI. This justifies our conservative choice in Section 3.3.6 to not correct for any remaining absolute calibration bias since we do not expect this effect to be as large on our final upper limits. However, in our end-to-end pipeline validation, we flag precisely where we did in the real data, ignoring any potential new discontinuities. These artificially low gains can thus be the consequence of a rare calibration failure getting spread out and diluted by gain smoothing.

We turn now to the final result of our end-to-end test: spherically averaged power spectra for both bands and all five fields, with and without our mock boosted EoR. In Figure 14, we show one field for each band (corresponding to our lowest limits; see Section 5). The power spectra include signal loss corrections for crosstalk subtraction and coherent time averaging, as described in Section 3.3.6, but no correction for nonredundancy as none is included in the simulation. For comparison, we also show the results of a simulation run with only the mock EoR with no foregrounds, noise, systematics, or flags (black line). It is interpolated directly onto the grid of the final LST-binned data set, time-averaged, converted to pseudo-Stokes I, and formed into power spectra.

Zoom In Zoom Out Reset image size

In general, we find that the results of our end-to-end test are in good agreement with the mock EoR-only power spectrum, as the bottom panels of Figure 14 show. Though the error bars on the power spectrum with EoR signal may be underestimates (since sample variance is not accounted for in our real pipeline), there is no evidence for an additional, unaccounted-for contribution to signal loss that might lead us to report artificially low upper limits. If anything, we see some evidence that our results are biased slightly high. This is likely due to a number of factors: the bias high due to rare absolute calibration failures, a possible overestimation of the impact of signal loss, and the effect of flagging and inpainting. When the same pipeline is run without the EoR, the results are consistent with noise at almost all k—a result of the somewhat aggressive cuts we discuss in Sections 3.2.4 and 3.3.3—though we unsurprisingly do see some residual foregrounds at very low k, especially in Band 1. This is consistent with what we see in the real data in most fields (see Section 5.2).

Separately from our end-to-end simulations, we also repeat two tests of signal loss performed in Aguirre et al. (2022) that were used in H22a to correct the final power spectrum measurements (see Section 3.3.6). The first is the effect of coherently time averaging LST-binned visibilities from 21.4–214 s before cross-multiplying interleaved visibilities to form power spectra. This is done by interpolating mock EoR-only visibilities onto the 21.4 s grid of the final LST-binned data set, forming pseudo-Stokes I, and then either forming power spectra directly or forming power spectra after coherently averaging to a 214 s cadence with rephasing as described in Section 3.1. If we take the result with 21.4 s integrations to be loss-free, which it should be to a good approximation, then we find a 1.2% signal loss for Band 1 and a 1.5% signal loss for Band 2 (see Figure 15).

Zoom In Zoom Out Reset image size

It is not surprising that the signal loss should be slightly higher for Band 2, since the primary beam is smaller at higher frequencies and the rephasing of visibilities before averaging cannot account for the changing beam-weighted sky. Likewise, if we were to open up additional degrees of freedom in our signal loss correction, we would likely find that this decoherence depends on baseline length and orientation—baselines that fringe more quickly should see more loss here. Since we are using the same weighting of baselines here as in analysis of the real data, the overall loss estimate should be correct to first order. Therefore, because the overall signal loss is quite small, it is not necessary to further complicate the signal loss correction by making it baseline dependent. Regardless, these values are consistent with the 1% loss figure used in H22a and calculated in Aguirre et al. (2022), which did not separate out the two bands.

The second specialized test is to examine the impact of crosstalk subtraction (described above in Sections 3.1 and 3.2.4). The crosstalk subtraction algorithm devised in Kern et al. (2019), demonstrated in Kern et al. (2020b), and employed for H22a removes power near fringe-rate zero. The maximum extent of this removal in fringe-rate space (Equation (3)) is designed to keep signal loss at the ∼1% level. To measure this effect, we interpolate to get one data set per epoch with foregrounds, mock EoR, and crosstalk injected, but no noise or calibration errors. To this data set we apply our final flagging mask, inpaint, subtract crosstalk, LST-bin the four epochs together, form pseudo-Stokes visibilities, coherently time average, and then form power spectra. We spherically average those power spectra over baselines using the same noise-based weights as were applied to the data. In Figure 16, we compare those per-time-step power spectra to mock EoR-only power spectra with the same averaging performed. We take a median over delay up to 4000 ns, the highest delay in the crosstalk subtraction SVD, to produce a single bias estimate per LST and per band. As H22a argues, we expect the crosstalk subtraction bias to be scale-independent.

Zoom In Zoom Out Reset image size

The final result is a bit difficult to interpret. Clearly the result is LST-dependent; just as in Aguirre et al. (2022), we see higher levels of signal loss near gaps in the data—this was one of the reasons we avoided the range from 0–0.75 hr of LST. We also see some evidence under-subtraction around Fornax A (in Band 1) and the Galactic center, though the latter is expected due the weighting in the SVD. To avoid the effect of outliers in both directions, we estimate our per-band signal loss by taking the median over LSTs in the five fields. This yields 2.4% signal loss in Band 1 and 3.3% loss in Band 2, which is basically consistent with the 1% and 3% used in H22a.

We begin with cylindrically averaged power spectra. These are first averaged over all cross-baseline pairs within each redundant-baseline group (excluding those with high residual crosstalk; see Section 3.2.4). Next they are averaged incoherently in time over each field. Finally, they are incoherently averaged in ∣ u ∣, combining together baselines with the same length. All averages use P_N (Equation (6)) to perform inverse-variance weighting; P_N and ${\hat{P}}_{{SN}}$ (Equation (7)) errors are propagated during each step. This produces power spectra in τ and u, which is equivalent in the delay approximation to k_∥ and k_⊥ (Parsons et al. 2012a). These cylindrical power spectra are the most sensitive data products that still keep k_∥ and k_⊥ separate, which is useful because different scales along those two axes are measured in a fundamentally different way instrumentally and correspond to different ways of measuring distance cosmologically. While we are not yet sensitive enough to constrain the cosmological signal's dependence on line-of-sight velocity effects, cylindrical power spectra are still useful for evaluating how foregrounds and other systematics appear upon deep integration.

We show our cylindrically averaged power spectra for each of the five fields in Figure 17 for Band 1 and Figure 18 for Band 2. In the top row of each figure, we show the real part of the power spectrum. ⁴⁹ We also show the power spectrum S/N, P/P_N in the bottom row. In regions that are noise-dominated, we expect roughly half of the power spectrum bins to be negative (white), and we expect the S/N to have mean 0 and standard deviation 1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

What we see is fairly consistent with that expectation across most of the EoR window. However, the foreground-dominated region clearly extends well beyond the horizon wedge (gray dashed line). This justifies the need for a buffer beyond the horizon wedge and our choice specifically to expand it from 200 ns as in H22a to 300 ns (black solid line). As we discussed in Section 3.3.3, our choice of the buffer was set by examining these S/N plots and trying to balance the amount of foreground leakage into the window at low k_⊥ with the amount of noise-dominated modes excised at high k_⊥.

Note that a number of data points that we excise from our spherically averaged power spectra are still shown in Figures 17 and 18. We have not yet removed the modes under the horizon wedge and buffer, nor have we removed baselines with projected east–west distance less than 15 m (though baselines with a projected east–west distance less than 14 m have already been removed because crosstalk subtraction could not be performed on them without substantial signal loss; see Section 3.1).

We now turn to our spherically averaged power spectrum measurements and upper limits. These are produced from averaging together all baselines incoherently after excising modes below the horizon wedge (Equation (8)) plus a 300 ns buffer and all baselines with a projected east–west length less than 15 m (which includes the single-unit east–west spacing; see Figure 4). The result, for all bands and fields, is shown as a "dimensionless" power spectrum (Equation (5)) in Figure 19. ⁵⁰

Zoom In Zoom Out Reset image size

We estimate vertical errors by propagating ${\hat{P}}_{{SN}}(k)$ (Equation (7)) and show the measurements ±2σ after converting to Δ² using Equation (5). When measurements are negative, we adopt the strong prior that the true power spectrum must be zero or positive and plot upper limits at precisely 2σ. The 1σ error level derived from ${\hat{P}}_{{SN}}(k)$ is shown as a black dashed line. Our absolute lowest 2σ upper limits on the power spectrum are Δ²(k = 0.36 h Mpc⁻¹) ≤ 3496 mK² at z = 10.4 (Band 1, Field D) and Δ²(k = 0.34 h Mpc⁻¹) ≤ 457 mK² at z = 7.9 (Band 2, Field C). These limits are 2.6 and 2.1 times deeper, respectively, than those presented in H22a. We report the measured power spectra, 1σ error bars, and 2σ upper limits results at most values of k—including low k, which is where most of our astrophysical constraining power comes from—in Table 3 for Band 1 and Table 4 for Band 2. ⁵¹

Table 3. Band 1 (z = 10.4) Power Spectra, Errors, and Upper Limits from Figure 19

	k	Δ2	1σ	${{\rm{\Delta }}}_{\mathrm{UL}}^{2}$	k	Δ2	1σ	${{\rm{\Delta }}}_{\mathrm{UL}}^{2}$	k	Δ2	1σ	${{\rm{\Delta }}}_{\mathrm{UL}}^{2}$
	(h Mpc−1)	(mK)2	(mK)2	(mK)2	(h Mpc−1)	(mK)2	(mK)2	(mK)2	(h Mpc−1)	(mK)2	(mK)2	(mK)2
Field A	0.23	153,691	2933	159,556	0.73	−32,771	58,377	116,753	1.22	588,465	274,340	1,137,146
	0.29	107,667	4271	116,208	0.79	61,900	74,547	210,994	1.28	406,694	318,088	1,042,869
	0.36	45,688	7105	59,898	0.85	43,307	93,159	229,625	1.35	669,812	366,722	1,403,255
	0.42	1182	11,336	23,854	0.91	158,855	115,453	389,760	1.41	−22,002	418,809	837,617
	0.48	−6756	17,231	34,462	0.97	−2569	140,221	280,442	1.47	721,387	476,106	1,673,599
	0.54	28,245	24,666	77,577	1.04	−25,200	167,921	335,842	1.53	190,383	537,247	1,264,877
	0.60	115,543	33,997	183,537	1.10	−206,058	199,136	398,273	1.59	−63,312	604,314	1,208,628
	0.66	101,535	45,017	191,569	1.16	412,511	235,335	883,180	1.65	440,694	676,991	1,794,676
Field B	0.23	6414	878	8169	0.73	−16,887	19,741	39,482	1.22	−3,270	90,946	181,893
	0.29	3352	1343	6038	0.79	38,162	25,246	88,654	1.28	−23,054	105,349	210,698
	0.36	4815	2346	9506	0.85	14,007	31,608	77,223	1.35	208,370	121,329	451,027
	0.42	18,593	3852	26,297	0.91	103,161	38,899	180,959	1.41	312,314	138,131	588,576
	0.48	26,651	5896	38,443	0.97	117,071	46,962	210,995	1.47	52,833	157,295	367,424
	0.54	54,508	8475	71,457	1.04	5188	55,544	116,277	1.53	28,767	178,728	386,224
	0.60	30,557	11,578	53,714	1.10	−21,580	66,145	132,291	1.59	202,144	200,619	603,381
	0.66	6089	15,308	36,705	1.16	−38,769	77,809	155,618	1.65	46,875	225,195	497,265
Field C	0.23	5565	472	6509	0.73	−5911	10,003	20,006	1.22	101,849	46,514	194,876
	0.29	3763	687	5137	0.79	−1040	12,731	25,461	1.28	−42,117	53,624	107,247
	0.36	1525	1185	3895	0.85	−555	15,963	31,926	1.35	109,485	62,018	233,520
	0.42	6806	1993	10,793	0.91	−35,281	19,656	39,312	1.41	−59,321	70,760	141,521
	0.48	−6007	3067	6134	0.97	−51,204	23,825	47,651	1.47	−42,848	80,462	160,924
	0.54	−7550	4320	8641	1.04	−27,920	28,398	56,795	1.53	148,253	90,978	330,208
	0.60	19,718	5898	31,514	1.10	43,689	33,761	111,212	1.59	227,669	102,654	432,978
	0.66	18,353	7755	33,864	1.16	43,764	39,661	123,085	1.65	130,546	114,700	359,945
Field D	0.23	13,006	677	14,360	0.73	−11,090	14,727	29,454	1.22	197,836	68,914	335,664
	0.29	2325	1011	4347	0.79	−5471	18,791	37,582	1.28	−45,122	79,627	159,254
	0.36	−207	1748	3496	0.85	−29,815	23,439	46,878	1.35	62,242	91,821	245,885
	0.42	6866	2888	12,643	0.91	−21,743	29,005	58,010	1.41	−95,510	105,107	210,215
	0.48	−1255	4423	8846	0.97	26,219	35,219	96,657	1.47	158,372	119,620	397,612
	0.54	25,094	6301	37,696	1.04	18,461	42,035	102,530	1.53	69,220	134,761	338,742
	0.60	17,496	8,595	34,687	1.10	112,195	50,065	212,326	1.59	−3969	151,666	303,332
	0.66	14,751	11,370	37,490	1.16	−46,370	58,914	117,827	1.65	−256,395	169,976	339,953
Field E	0.23	28,531	957	30,444	0.73	−29,956	20,066	40,133	1.22	−4695	94,165	188,330
	0.29	14,846	1394	17,634	0.79	27,200	25,575	78,350	1.28	70,547	109,087	288,722
	0.36	6793	2382	11,557	0.85	45,734	32,039	109,812	1.35	32,410	125,598	283,606
	0.42	2510	3946	10,401	0.91	2109	39,559	81,228	1.41	220,926	143,720	508,367
	0.48	−2082	6016	12,031	0.97	48,558	48,118	144,795	1.47	129,091	163,068	455,227
	0.54	8265	8530	25,324	1.04	−99,003	57,725	115,450	1.53	370,465	185,109	740,683
	0.60	−21,997	11,634	23,267	1.10	−58,513	68,608	137,217	1.59	153,201	207,953	569,107
	0.66	14,858	15,432	45,722	1.16	−78,339	80,512	161,025	1.65	175,314	232,387	640,089

Note. The upper limit is taken to be 2σ above the power spectrum measurement or above 0, whichever is greater. The lowest upper limit is in bold.

Download table as:  ASCIITypeset image

Table 4. Band 2 (z = 7.9) Power Spectra, Errors, and Upper Limits from Figure 19

	k	Δ2	1σ	${{\rm{\Delta }}}_{\mathrm{UL}}^{2}$	k	Δ2	1σ	${{\rm{\Delta }}}_{\mathrm{UL}}^{2}$	k	Δ2	1σ	${{\rm{\Delta }}}_{\mathrm{UL}}^{2}$
	(h Mpc−1)	(mK)2	(mK)2	(mK)2	(h Mpc−1)	(mK)2	(mK)2	(mK)2	(h Mpc−1)	(mK)2	(mK)2	(mK)2
Field A	0.27	9331	649	10,628	0.83	−14,403	14,954	29,908	1.40	63,705	68,760	201,225
	0.34	6479	1027	8532	0.90	12,644	18,997	50,637	1.47	104,906	79,266	263,438
	0.41	4731	1774	8279	0.97	−7283	23,849	47,699	1.54	−115,042	91,161	182,322
	0.48	−2629	2857	5713	1.05	76,346	29,519	135,384	1.61	−93,039	104,349	208,699
	0.55	1803	4376	10,556	1.12	−4423	35,449	70,899	1.68	−80,538	118,562	237,124
	0.62	−12,812	6296	12,592	1.19	−95,045	41,893	83,785	1.75	227,088	134,070	495,227
	0.69	−2143	8621	17,243	1.26	−19,638	49,987	99,974	1.83	157,690	151,649	460,988
	0.76	−4635	11,486	22,973	1.33	−44,250	58,807	117,614	1.90	576,790	169,926	916,643
Field B	0.27	747	229	1204	0.83	10,939	5601	22,140	1.40	61,788	25,729	113,245
	0.34	1149	377	1902	0.90	4164	7130	18,424	1.47	48,349	29,834	108,017
	0.41	−770	657	1314	0.97	24,889	8943	42,774	1.54	23,599	34,335	92,268
	0.48	1254	1080	3415	1.05	10,417	11,027	32,470	1.61	15,477	39,195	93,868
	0.55	−1770	1655	3310	1.12	−33,006	13,209	26,418	1.68	52,495	44,690	141,876
	0.62	1916	2369	6654	1.19	−11,140	15,659	31,319	1.75	125,943	50,570	227,083
	0.69	9099	3231	15,562	1.26	33,775	18,651	71,077	1.83	−41,091	56,640	113,280
	0.76	11,389	4315	20,019	1.33	2963	22,054	47,071	1.90	−14,827	63,497	126,993
Field C	0.27	330	129	587	0.83	8578	3048	14,674	1.40	19,670	13,891	47,452
	0.34	44	206	457	0.90	10,570	3868	18,306	1.47	−14,754	16,094	32,188
	0.41	−258	360	720	0.97	1360	4824	11,009	1.54	−7811	18,511	37,022
	0.48	736	595	1925	1.05	12,516	5958	24,431	1.61	−40,971	21,178	42,355
	0.55	1174	921	3016	1.12	−13,271	7158	14,315	1.68	−41,603	24,006	48,012
	0.62	−1277	1305	2609	1.19	−18,287	8473	16,946	1.75	−21,962	27,134	54,269
	0.69	−379	1778	3556	1.26	−23,016	10,058	20,116	1.83	25,762	30,636	87,034
	0.76	1636	2341	6317	1.33	4083	11,855	27,794	1.90	−47,835	34,297	68,593
Field D	0.27	413	181	775	0.83	3556	4295	12,146	1.40	−5572	19,750	39,500
	0.34	1912	294	2500	0.90	−609	5482	10,963	1.47	38,934	22,954	84,841
	0.41	2898	511	3921	0.97	−12,554	6839	13,678	1.54	13,747	26,422	66,590
	0.48	1079	836	2751	1.05	−3677	8465	16,931	1.61	15,991	30,140	76,271
	0.55	1161	1282	3725	1.12	3285	10,238	23,761	1.68	−28,833	34,159	68,319
	0.62	4365	1824	8013	1.19	29,719	12,123	53,965	1.75	14,543	38,550	91,642
	0.69	6178	2487	11,151	1.26	15,449	14,403	44,255	1.83	95,139	43,587	182,313
	0.76	9354	3307	15,967	1.33	10,098	16,922	43,942	1.90	68,406	48,901	166,208
Field E	0.27	1437	240	1918	0.83	23,107	5493	34,092	1.40	−616	25,170	50,340
	0.34	984	375	1733	0.90	12,550	6995	26,541	1.47	−40,653	29,194	58,387
	0.41	507	649	1806	0.97	−6756	8755	17,510	1.54	40,604	33,594	107,793
	0.48	−1879	1070	2141	1.05	1256	10,779	22,815	1.61	6847	38,395	83,636
	0.55	−3158	1632	3264	1.12	−7231	13,009	26,017	1.68	−2468	43,729	87,457
	0.62	−5496	2316	4631	1.19	27,219	15,431	58,081	1.75	15,566	49,408	114,382
	0.69	−542	3164	6328	1.26	11,508	18,321	48,149	1.83	26,283	55,545	137,374
	0.76	5617	4209	14,036	1.33	10,852	21,545	53,943	1.90	18,011	62,323	142,657

Note. The upper limit is taken to be 2σ above the power spectrum measurement or above 0, whichever is greater. The lowest upper limit is in bold.

Download table as:  ASCIITypeset image

In general, most of our measurements are consistent with the expected noise level. At low k, especially in Band 1, we see evidence for residual foregrounds beyond our horizon plus buffer cut. This makes sense; foregrounds are brighter at low frequency. For a few fields and bands, there is some evidence for residual crosstalk, which should appear between roughly 0.4–1.0 h Mpc⁻¹. Close inspection of Tables 3 and 4 also reveals a handful of large negative power spectra at the roughly −2σ level, especially in Band 2, Field E. There are more such points than might be expected from random noise, even after accounting for that fact that, as H22a showed, errors are correlated between neighboring k bins at the ∼25% level.

This is potentially concerning, since ignoring auto-baseline-pairs as we do here and as was done in H22a has the potential to introduce negative power spectrum biases (Morales et al. 2022), which would then lead to artificially low upper limits. Tracing back these points to their corresponding cylindrical bandpowers in Figures 17 and 18 reveals only a handful of negative power spectra in excess of what is expected from thermal noise, the largest of which are all near ∣τ∣ = 900 ns. Because these excesses appear at higher k than the best limits and because we conservatively chose to use 2σ as the upper limit wherever Δ² is negative, it is unlikely that our deepest power spectra or the likelihoods used in Section 7 are significantly affected. That the systematics in this range remain mostly marginal speaks to the efficacy of our crosstalk subtraction and baseline selection techniques (see Section 3.2.4).

We use horizontal error bars in Figure 19 to represent the range of k values whose contribution sums to 68% of each measurement. To examine both the window functions and our deepest upper limit more closely, we reproduce our Band 2, Field C results from Figure 19 in the upper panel of Figure 20, along with the window functions themselves in the lower panel. The window function matrix W , defined in Equation (9) and propagated to the spherical power spectra, tells us the extent to which different bandpower measurements are expected to be the linear combination of the true bandpowers across k modes. Each row—shown as a different colored line in the bottom row of Figure 20—sums to 1 by definition. The fact that each row is peaked along the diagonal of the matrix means dominant contribution to each measurement is in fact the k mode at which the measurement is reported.

Zoom In Zoom Out Reset image size

As discussed in Section 3.3.4, these window functions more precisely take into account the delay approximation than the simpler calculation used in H22a. While they may look similar to the window functions presented in H22a because of the linear scale used in the figure, there are as many as 3 orders of magnitude of difference in the amplitude of the tails, with the delay approximation leading to a large underestimate (Gorce et al. 2023). This is particularly important if the true 21 cm power spectrum has features in k space or deviates substantially from a roughly flat power law in Δ²(k); although, that is not the case in most "vanilla" models of reionization.

As in H22a, we performed a series of noise integration tests to determine whether bandpowers outside the wedge region are consistent with being noise-dominated. By cumulatively averaging the samples that went into the fully averaged power spectra shown in Figure 19, we are able to test whether the samples average together as would be expected for bandpowers formed from uncorrelated (white) random noise, or whether an additional correlated or nonrandom signal, such as a source of systematic contamination, may be present.

Our first test is to simply integrate all bandpowers above a fixed k. Following H22a, we randomly draw 300,000 pure-noise realizations of the bandpowers using the error covariance matrix propagated through our various averaging steps and used to set the error bars on our power spectra. We then compute the fraction of noise realizations that are larger than the data, which is by definition the p-value of the measurement if the null hypothesis states that measured values are drawn from the noise covariance and have mean zero. In Table 5, we show the results of this test for both the real and imaginary parts of both bands' and all fields' power spectra.

In general, we find our imaginary power spectra consistent with noise across k, as expected. The one exception is in Field A for Band 1, where we are likely seeing noise-foreground cross terms dominating the imaginary power spectrum. Our real power spectra are another matter. We see strong evidence for inconsistency with noise at low k, which is expected given the foreground leakage just outside the wedge buffer in both cylindrical power spectra (Figures 17 and 18) and at low k in spherical power spectra (Figure 19). In Band 1, we also see more marginal evidence for inconsistency with noise even after k < 0.5 hMpc⁻¹ is excluded. This is likely attributable to residual crosstalk, which appears in our final power spectra at the few-σ level.

While the p-values are a useful test of the null hypothesis on the final power, a lot can be obscured by collapsing an entire power spectrum to just a few numbers. To further probe the consistency of our results with noise, especially at high delay, we perform two additional tests. The first assesses how baselines integrate down within a redundant group. In Figure 21 we examine how a single delay mode at τ ≈ ±3000 ns integrates down as more baselines are added to the cumulative average. We use our most sensitive redundant group: the 29 m east–west baseline. We normalize the cumulative average by the expected variance of that cumulative average derived from P_N . We show 1σ and 2σ regions where we expect 68% and 95% of pure-noise realizations to fall in gray. We have validated this calculation with a set of 100 white noise-only simulations that have been passed through the same flagging, weighting, and power spectrum estimation steps as the real data. The simulations are constructed to have noise variance that depends on time, frequency, and baseline in the same way as the data, according to our empirical estimates of the noise power spectrum P_N .

Zoom In Zoom Out Reset image size

The overall conclusion from Figure 21 is that there are no strong deviations from the expected noise-like behavior for any band or field in our data set at ∣τ∣ ≈ 3000 ns, according to this statistic. Band 1, Field D has the most substantial deviations, but these remain largely within the 95% region. Band 2, Field A also approaches the edge of the 95% region. Bands 1 and 2 of Field C have averages that are more consistently close to zero than for any of the others, though this is not necessarily a statistically significant anomaly.

Another similar test is shown in Figure 22 where we examine the same baseline group, but average over baselines in the group and over a wider range of delays ( ± 2000–4000 ns) to look for evidence of LST-dependent systematics. We compare each average to the distribution of 100 noise realizations that are drawn from P_N and reflect the real data's sampling and cuts. In the fields of interest, we see little evidence for significant outliers, indicating that there is no strongly LST-dependent high-delay systematic affecting our results for this baseline group. Outside our fields, we do see some large outliers from Fornax A and the Galactic center, which might be concerning in the future if they appear elsewhere in HERA data as we integrate deeper.

Zoom In Zoom Out Reset image size

Though we have taken several steps to mitigate the effect of crosstalk—including subtracting it from our visibilities and cutting baselines that exhibit substantial residuals (see Section 3.2.4)—it is useful to understand how those choices may have affected the distribution of bandpowers. In H22a we performed a one-sample Kolmogorov–Smirnov (K-S) test comparing the cumulative distribution function (CDF) of power spectra over 0.4 hr of LST to a validated analytic noise model.

Here we iterate slightly on that test with a two-sample K-S test comparing the CDF of measured bandpowers for the redundantly averaged 29 m east–west baseline group to that of the same 100 noise realizations as in our z-score test above. Each CDF is computed over LSTs in Field C and over a series of 200 ns-wide delay bins (except the first bin, which spans 0–500 ns). In Figure 23, we show the results of that work for both bands and for data with (dashed lines) and without crosstalk subtraction (solid lines).

Zoom In Zoom Out Reset image size

To infer the expected range of distances between CDFs, we also compute a distribution of two-sample K-S statistics for each delay ranges by comparing many pairs of independent noise realizations. We show the 2σ range of that distribution in the gray band in Figure 23. We can see that while the distribution of bandpowers was highly inconsistent with noise before crosstalk subtraction, afterwards they are consistent beyond 500 ns for this baseline group.

In order to check the consistency of the epochs (defined in Section 2.1) at the bandpower level and thus justify (post hoc) the decision to combine them to arrive at our final limits, it is useful to begin by looking at single-epoch power spectra. These are formed by simply going straight from crosstalk subtraction to forming pseudo-Stokes I and coherently averaging. In Figure 24, we show one such power spectrum from Band 2, Field C—our deepest limit—and compare measurements and vertical error bars from individual epochs to the full data set. The full data set's power spectrum (which is identical to that in Figure 20), is generally lower than any individual epoch. Since most individual epochs are consistent with noise, this is not too surprising. Since the epochs are binned together with weights proportional to the number of visibility samples, it is not trivial to look at this figure and read off the impact of each epoch on any given k mode. Therefore, while there are no obvious inconsistencies visible here, it behooves us to approach the question more quantitatively.

Zoom In Zoom Out Reset image size

We answer that question by performing a Bayesian jackknife test across epochs. This test is described in detail in Wilensky et al. (2022); we summarize our specific implementation of it here. For each band, field, and spherically averaged k-mode, we consider 2^N hypothetical models for the N ≤ 4 epochs that contribute to the particular band and field's power spectrum at the specified k-mode (Epoch 0, for example, does not contain any observations in Field E). In each hypothetical model, we propose that a bandpower is either unbiased—i.e., consistent with zero-mean noise—or strongly biased—i.e., consistent with measuring a signal—by an unspecified amount within a relatively constrained prior range. To compute the posterior for each model, we marginalize over a multivariate Gaussian bias prior that is centered at 6σ_i with a width of σ_i in each direction, where σ_i is the square root of the estimated variance for the ith bandpower. This is intended to represent a strong but not excessive bias with an a priori S/N in the 4–8 range. We then use a series of Bayesian model comparisons to identify which epoch(s) is most likely to be biased for any given bandpower.

Figure 25 shows the results of this test for each band, field, and k mode. To calculate the inferred odds of there being at least one significantly biased epoch—albeit only when the large biases reflected by the alternate hypotheses are considered—we have evaluated the posterior probability of each bias configuration, summed those in which at least one bias is present, and divided by the posterior probability of the null hypothesis that no biases were present. We use a flat prior for the bias configurations. We generally find that for the majority of bands and fields, most spherically averaged Fourier modes are more consistent with pure noise than a scenario with an epoch-dependent bias, particularly for higher-k modes (see the caption of Figure 25 for a discussion of fields/epochs that potentially deviate from the null hypothesis). Since crosstalk subtraction is performed on a per-epoch basis, it is likely that residual crosstalk is a major source of bias that varies from epoch to epoch. This majority-null result justifies the decision to average the epochs together into one final power spectrum for each band and field.

Zoom In Zoom Out Reset image size

This is not to say that we strongly suspect there is no underlying signal due to residual systematics or one of cosmological origin. Since only large biases are considered in the alternate hypotheses, our conclusion is that if biases are present, then they are sufficiently small to be difficult to distinguish from the expected statistical fluctuations in the bandpowers. Indeed, if we marginalize over bias priors centered at smaller biases (e.g., less than σ_i ), the posterior probability of the bias configurations is more evenly diffused over the hypotheses, demonstrating a lack of certainty that arises from an inability to finely distinguish weak bias configurations with so few data.

To interpret the power spectra and upper limits reported in H22a, all four techniques examined in H22b employ the same statistical approach to evaluating the posterior probability of model parameters in a Bayesian framework. Each theoretical model ${ \mathcal M }$ with a set of parameters θ compares the data d to modeled power spectra m ( θ ) convolved with the window function W . We write this difference as t ≡ d − W m ( θ ). We next assume that each measurement is due to some unknown combination of 21 cm signal, noise, and systematics. Further we assume that systematics (typically residual foregrounds and crosstalk) can only add power. These assumptions are only appropriate when claiming upper limits on the 21 cm power spectrum, as we do in this work. To claim a detection, one would have to impose additional constraints on the relative contribution of systematics to the measurement. To do that credibly, one likely needs high-S/N detections at multiple k modes, redshifts, and fields, which can be subjected to rigorous jackknives and other statistical tests for internal consistency, along with a theoretical framework able to match the data well.

Marginalizing over systematics, as we show in H22b, yields a posterior probability of the form

$$\begin{eqnarray}&&p({\boldsymbol{\theta }}| {\boldsymbol{d}},{ \mathcal M })\propto p({\boldsymbol{\theta }}| { \mathcal M })\displaystyle \prod _{i}^{{N}_{d}}\displaystyle \frac{1}{2}\left(1+\mathrm{erf}\left[\displaystyle \frac{{t}_{i}}{\sqrt{2}{\sigma }_{i}}\right]\right),\end{eqnarray} \tag{ 11 }$$

where the product on the right-hand side is the likelihood. Here N_d is the number of data points considered and σ_i is the standard deviation of d_i . This result is consistent with similar approaches in the literature (Li et al. 2019; Ghara et al. 2020). In this derivation, we also assumed that measurements have uncorrelated errors. While neighboring measurements in k are in fact quite correlated, we follow H22b and throw out every other measurement in k—keeping the absolute lowest limits—to largely eliminate that correlation.

To understand the effect of the form of the likelihood in Equation (11), we calculate it using a power law Δ² ∝ k^α over several orders of magnitude in Δ² evaluated at k = 0.35 h Mpc⁻¹ in Figure 26. At z = 10.4, we use α = 1.3; at z = 7.9 we use α = 0.0. These values of α were derived by interpolating the 21cmMC power spectrum realizations in our posterior without HERA (see Section 7.4.3) at the k values of our two deepest measurements, fitting a power law, and then taking the median α over models that had not yet completely reionized. The flattening of Δ² with decreasing redshift is a generic consequence of inside-out reionization, which suppresses small-scale power.

Zoom In Zoom Out Reset image size

The first thing to note is that our likelihoods are essentially flat for Δ² values much less than the upper limit; the signal could fall anywhere in that range because the systematics could fall anywhere in that range as well. Our results can therefore be heavily impacted by the prior, $p({\boldsymbol{\theta }}| { \mathcal M })$, and what it considers the smallest viable power spectrum.

Figure 26 also highlights the fact that, in contrast to H22b, our result features fairly comparable limits from multiple fields and k modes. Because our likelihoods are not Gaussian, the results do not scale like $1/\sqrt{N}$ as one might expect. Limits with larger measurements and smaller error bars look different in Equation (11) than measurements with the same 2σ upper limits but with smaller measured Δ² and larger error bars. So, while it is not surprising that most of the information in Band 2 comes from Field C, it is interesting and not necessarily intuitive that Fields C and D contribute equal amounts of information in Band 1. Likewise, we can see that just using a single k mode (dotted black lines)—even though it is the single best upper limit and combines together all fields—yields somewhat worse constraints than combining multiple ks. The vast majority of the constraining power comes from k < 1 h Mpc⁻¹ (solid black line). Of course, combining k modes together in this way is only technically correct if we are evaluating the likelihood in Equation (11) for a model with Δ² ∝ k^α . In general, our model posteriors use all fields and all k-values, but compare them to the data at the proper k values.

In H22b, we interpreted the results of H22a using a suite of four different theoretical models to infer constraints on the IGM and high-z galaxies:

• 1.  
  a simple "density-driven" linear bias model, in which the 21 cm fluctuations are assumed to trace the density fluctuations, multiplied by a bias factor that depends on the average thermal and ionization state of the IGM;
• 2.  
  a phenomenological "reionization-driven" model, which parameterizes the ionized bubble distribution and IGM spin temperature directly without making any explicit assumptions about galaxies (Mirocha et al. 2022);
• 3.  
  a seminumeric model of the 21 cm signal, 21cmFAST (Mesinger et al. 2011), along with the inference engine 21cmMC (Greig & Mesinger 2017) that uses it to explore and constrain a range of parameterized galaxy properties (e.g., Park et al. 2019) with multiwavelength probes; and
• 4.  
  an independent seminumeric model that models the galaxies differently and allows for a radio background in excess of the CMB at high redshifts (Reis et al. 2020).

Our models are constructed in two qualitatively different ways. The density-driven and reionization-driven models interpret the 21 cm signal directly as a function of IGM properties. They avoid making explicit assumptions about the sources generating the radiation fields that determine the 21 cm signal (though their priors on IGM properties carry some implicit assumptions), and as such have a lighter implementation. The seminumeric models, on the other hand, start with a model of galaxy evolution and simulate a cosmological volume in order to predict the 21 cm signal. The latter are significantly more complicated but have the distinct advantage of a parameter set that is rooted in our current understanding of galaxy properties, which allows us to combine 21 cm measurements with other constraints on the galaxy population (such as UV luminosity functions and the X-ray background). The disadvantage of the seminumeric approach is that the resulting constraints on galaxy physics are only sensible if their parameterization is flexible enough to include all of the relevant physics. In H22b we considered two independent seminumeric models to mitigate this problem.

The key result of H22b, across all of these models, was that the IGM must have been heated above the expectation for an adiabatically cooling IGM from recombination to z = 7.9 at >95% credibility. ⁵² If this heating is interpreted in the conventional way—as the result of X-rays generated by accretion onto black hole remnants of early star formation and then depositing that thermal energy in the IGM—this level of heating suggested that early galaxies were more efficient X-ray emitters than their local counterparts. If an excess radio background were allowed, the observations jointly constrained the efficiencies of radio and X-ray emission from galaxies, also restricting the parameter space of otherwise viable models of these mechanisms during the EoR.

In this work, we apply those same four techniques—with some minor adjustments, as we will discuss below—to our data from both bands. As H22b argued, the easiest point of comparison for the models is the inferred average spin temperature of the neutral IGM, ${\overline{T}}_{S}$. In Figure 27 we compare the 68% and 95% highest posterior density (HPD) credible interval on ${\overline{T}}_{S}$ in all four models and show how those constraints have improved since H22b. ⁵³ The improvements are modest but consequential, especially at z = 10.4: all four models now independently require the IGM to be heated above the adiabatic cooling limit at >95% credibility. The precise values of lower limits vary from model to model, reflecting their different assumptions and approximations. In particular, 21cmMC shows the most dramatic increase in the spin temperature constraint at z = 7.9, a fact that deserves further examination. ⁵⁴ As we discuss the specific models and their detailed results throughout this section, we will look to understand what precisely drives these differences. In Section 7.6, we will summarize our takeaways and put these spin temperature constraints in the broader cosmological context.

Zoom In Zoom Out Reset image size

Before we discuss the detailed results of our more complex techniques for astrophysical inference based on seminumerical simulations, here we briefly detail the two simpler methods used in H22b to infer constraints on the spin temperature and other IGM properties.

In the first, which we refer to as the "density-driven" model, the 21 cm power spectrum is assumed to follow that of matter, multiplied by a bias parameter squared (as is standard in perturbation theory; see, e.g., McQuinn & D'Aloisio 2018; Georgiev et al. 2022; Qin et al. 2022; Sailer et al. 2022 for more complex approaches). Limits on this bias can then be translated into lower bounds on ${\overline{T}}_{S}/{T}_{\mathrm{CMB}}$ under some assumptions, including redshift-space distortions as a function of the line-of-sight cosine μ (see H22b) and that the IGM properties are roughly homogeneous (as the model ignores ionization and temperature fluctuations beyond adiabatic; see H22b for details). Our limits on the bias b_m using 94 nights (and all fields and every other k, as discussed above) are ∣b_m ∣ < 60 mK at z = 7.9 and ∣b_m ∣ < 160 mK at z = 10.4 at 95% credibility, which assuming x_HI = 1 translate into the IGM limits shown in Figure 27 and in Figure 34. Note that these limits are on the absolute value ∣b_m ∣. This is unimportant for cases where ${\overline{T}}_{S}\ll {T}_{\mathrm{radio}}$, but as the limits tighten, we will obtain contours around ${\overline{T}}_{S}={T}_{\mathrm{radio}}$ (which yields no 21 cm signal and thus b_m = 0). In fact, the 68% credible limit on the z = 7.9 bias is ∣b_m ∣ < 30 mK, which translates into the double-sided region from ${\overline{T}}_{S}=15$ K to ${\overline{T}}_{S}=60$ K (assuming T_radio = T_CMB, as is standard).

In the second, which we refer to as the "reionization-driven" model, the IGM is modeled as a two-phase medium, with fully ionized bubbles drawn from a log-normal size distribution and the "bulk" IGM outside bubbles assumed to be of uniform temperature Mirocha et al. 2022. The advantage of this approach ⁵⁵ is that it works directly with IGM quantities, which makes it easy to interpret, and avoids making explicit assumptions about galaxies. ⁵⁶ In its simplest form, it requires four parameters: the volume filling fraction of ionized gas, Q ≡ (1 − x_{H I}), the IGM spin temperature, ${\overline{T}}_{S}$, the characteristic bubble size, R_b , and log-normal dispersion, σ_b . In this work, when jointly fitting both HERA bands, we also use a seven-parameter version of the model in which the ionized fraction and characteristic bubble size are allowed to evolve with redshift as power laws. We require only that Q and R_b increase with time, that reionization completes at z ≥ 5.3, and that all parameters are positive. We perform our Markov Chain Monte Carlo (MCMC) inference using emcee (Foreman-Mackey et al. 2013).

Now, with 94 nights of data, we infer spin temperatures in excess of the adiabatic cooling limit at z = 7.9 and 10.4, both for fits that consider each band separately and the joint fit to both bands using a seven-parameter version of the model (see Figure 27). At 95% (68%) credibility, we obtain ${\overline{T}}_{S}\gt 11.0$ (35.2) K at z = 7.9, and ${\overline{T}}_{S}\gt 6.2$ (26.4) K at z = 10.4, when fitting both bands simultaneously. Note that, in this joint fit, lower limits on ${\overline{T}}_{S}$ grow at z = 7.9, as expected, but actually slightly decrease at z = 10.4 (at 95% credibility) relative to the results of single-band fits. This is a result of using an HPD estimate of the credible interval, combined with change in the shape of the ${\overline{T}}_{S}$ posterior, which goes from being roughly flat as ${\overline{T}}_{S}\to {10}^{3}\ {\rm{K}}$ to a more "peaked" distribution, slightly away from the maximal value of ${\overline{T}}_{S}$. This is consistent with the results of the density-driven model, suggesting that our limits are beginning to disfavor scenarios with saturated ${\overline{T}}_{S}$ at z = 7.9.

The differences between these approaches explains some of the differing conclusions in Figure 27. For example, the density-driven model yields higher ${\overline{T}}_{S}$ limits than the reionization-driven model, because the latter assumes ionization and density are positively correlated. As a result, the IGM must be made colder to compensate for the loss of the densest regions to ionization. We see also in the reionization-driven model the power of jointly fitting multiple bands; the z = 7.9 limit increases by roughly a factor of 2 in the joint fit because a given temperature at z = 7.9 is only viable if the temperature evolution is consistent with the z = 10.4 data.

21cmFAST (Mesinger et al. 2011; Murray et al. 2020) is a seminumerical simulator for computing the evolution of the 21 cm signal across cosmic time by assuming that star-forming galaxies, hosted by dark matter halos, drive the cosmic radiation fields that heat and ionize the IGM. It uses empirical scaling relations to assign galaxy properties to dark matter halos based on their halo masses. These include the stellar mass to halo mass ratio, the UV ionizing escape fraction, the star formation rate, and the X-ray luminosity. 21cmFAST simulates reionization using the excursion-set formalism (Furlanetto et al.2004). It accounts for inhomogeneous recombinations (Sobacchi & Mesinger 2014) and includes a numerical correction for photon conservation (Park et al. 2021).

In Section 5.1 of H22b, we discussed the nine physically meaningful free parameters that go into the scaling relations describing galaxy properties (also see Park et al. 2019). We adopted either flat linear or logarithmic priors—depending on the parameter's dynamic range and how it enters into the model—within our 9D hypercubic parameter space. The ranges on these priors were picked to allow a broad range of physically plausible values while not strongly constraining the parameters most constrained by other high-redshift probes. These include UV luminosity functions (Bouwens et al. 2016, 2015; Oesch et al. 2017) and measurements of the IGM opacities during the EoR such as the Lyα forest (Bosman & Becker 2015; McGreer et al. 2015; Qin et al. 2021a) and CMB polarization and optical depth (Planck Collaboration 2020; Qin et al. 2020). The Lyα forest and the CMB's τ provide important clues as to the timing and duration of reionization, which then constrain the ionizing escape fraction, given the observed UV luminosity function. We use 21cmMC (Greig & Mesinger 2015, 2017) and its MultiNest sampler (Feroz et al. 2009; Qin et al. 2021b) to perform Bayesian inference in this work. The posterior probability distribution without HERA serves as the starting point for comparing models against HERA data.

One key result from these other high-redshift probes is the strong constraint on the star formation efficiency of dark matter halos at high redshift, which is determined by the UV luminosity functions at z ∼ 6–8; existing observations constrain both the peak efficiency and show that it declines toward small halo masses (e.g., Tacchella et al. 2013; Mason et al. 2015; Mirocha et al. 2017; Park et al. 2019; Sabti et al. 2022). The version of 21cmFAST used here assumes that this behavior can be extrapolated to higher redshifts and smaller halos. As a result, the range of galaxy formation models that are allowed by 21cmMC is relatively restricted, and all display a rapid increase in the stellar mass density between the two redshifts measured by HERA. The resulting constraints must be interpreted with this in mind, as more complex galaxy evolution histories (which break the assumed extrapolation by, for example, introducing a new population of sources, see, e.g., Muñoz 2019; Qin et al. 2020; Muñoz et al. 2022 for implementations in 21cmFAST) are not included in the prior and are left for future work.

7.4.1.  21cmMC Constraints on X-Ray Luminosity

In H22b we explored how adding HERA affected the full posterior parameter covariance. In this work, we focus on the parameter most constrained by HERA, the ratio of the integrated soft-band (<2 keV) X-ray luminosity to the star formation rate. Since 21cmFAST assumes that X-ray photons govern the thermal history of the neutral IGM, this L_{X<2 keV}/SFR parameter essentially describes the heating power of EoR galaxies per unit of star formation. In H22b, we obtained the first observational evidence for an enhanced X-ray luminosity of high-redshift (z > 6) galaxies, with a 68% HPD credible interval of L_{X<2 keV}/SFR ∼ 10^39.9–10^41.6 erg s⁻¹ ${M}_{\odot }^{-1}$ yr. This disfavors a relationship between star formation and soft X-ray luminosity at or below the one seen in local, metal-enriched galaxies at >68% credibility.

As Figure 28 shows, we find that the full season of HERA observing constrains the 95% credible interval on L_{X<2 keV}/SFR to the range 10^40.4–10^41.8 erg s⁻¹ ${M}_{\odot }^{-1}$ yr. This result assumes as a prior that L_{X<2 keV}/SFR < 10⁴² erg s⁻¹ ${M}_{\odot }^{-1}$ yr, beyond which X-rays reionize the universe too quickly (Mesinger et al. 2013). More than 99% of the 21cmMC posterior volume excludes the possibility of the local relation for high-mass X-ray binaries (HMXBs; Mineo et al. 2012) continuing to hold at high redshift. It is consistent, however, with models of extremely low-metallicity galaxies, where high-mass stars have less mass loss from line-driven winds than their solar-metallicity counterparts (Fragos et al. 2013).

Zoom In Zoom Out Reset image size

7.4.2.  21cmMC Constraints on the IGM's Thermal History

Our constraints on the soft X-ray efficiency are themselves a consequence of our ability to use our upper limits to exclude a range of scenarios with low levels of IGM heating. In Figure 29 we show our updated marginalized posteriors on the predicted average IGM temperatures—both the spin temperature, ${\overline{T}}_{S}$, and the kinetic temperature, ${\overline{T}}_{K}$—along with results from H22b for comparison. As demonstrated in H22b, current EoR constraints from Planck and quasar spectra already disfavor a large number of models in the prior volume that predict either highly ionized IGM at z ≥ 10.4 or a completely neutral one at z ≤ 10.4. These constraints also have a slight impact on the average IGM temperature, excluding models with high ${\overline{T}}_{{\rm{K}}}$ or ${\overline{T}}_{{\rm{S}}}$ at these redshifts. However, because a decently sized fraction of parameter space with an unheated IGM at these redshifts is not ruled by these probes, and since 21cmMC cannot produce spin temperatures below the adiabatic limit, our posterior without HERA shows a pileup of probability right around that limit.

Zoom In Zoom Out Reset image size

When we incorporate the HERA limits, a significant range of models with low IGM temperatures can be further excluded. We showed in H22b how HERA observations substantially improve our understanding of the neutral IGM at z = 7.9. However, there was still a small fraction of the total posterior with low values of ${\overline{T}}_{{\rm{S}}}$, so H22b could not completely rule out an unheated IGM the observed redshifts. With the improved limits presented in this work, we now find that an unheated IGM is disfavored at greater than 99% credibility at both z = 10.4 and 7.9. The new HPD 95% credible intervals on the spin and kinetic temperatures are 4.7 ${\rm{K}}\lt {\overline{T}}_{S}\lt 171.2$ K and 3.2 ${\rm{K}}\lt {\overline{T}}_{K}\lt 313.2$ K at z = 10.4 and 15.6 ${\rm{K}}\lt {\overline{T}}_{S}\lt 656.7$ K and 13.0 ${\rm{K}}\lt {\overline{T}}_{K}\lt 4768$ K at z = 7.9. ⁵⁷

Since the 21 cm brightness temperature is proportional to (1 − T_CMB/T_S), any model where T_S < 0.5T_CMB will show an enhanced power spectrum amplitude relative to one where T_S ≫ T_CMB. Since the CMB is 24.3 K at z = 7.9, our results rule out the range of cold reionization scenarios with enhanced power spectra at z = 7.9 with greater than 95% credibility. While we cannot definitively rule out any power spectrum enhancement due to a cold IGM z = 10.4, we have ruled out the large swath of parameter space that produces the largest power spectra.

Figure 29 brings into sharper relief the question raised by Figure 27 of why the 21cmMC results show such a large change in ${\overline{T}}_{S}$. Once we reach the regime where ${\overline{T}}_{S}\gt {T}_{\mathrm{CMB}}$, large increases in spin temperature have only a modest impact on the power spectrum. It follows then that the updated power spectrum limits at z = 7.9 are not primarily driving the spin temperature constraint at z = 7.9. Just as in the reionization-driven phenomenological model where jointly fitting the two redshifts produced tighter ${\overline{T}}_{S}$ constraints, the combination of our two redshift bands is crucial. In particular, the star formation histories inferred by 21cmMC from UV luminosity functions (in which the stellar mass density increases significantly from z = 10.4 to z = 7.9) causes substantial IGM heating between the two bands. Thus the modest heating demanded by HERA at z = 10.4 translates to a strong inference about the temperature at the later times. We will see quite a similar effect when we examine the inferred power spectrum posteriors in Figure 31.

7.4.3.  21cmMC Derived Constraints on the 21 cm Power Spectrum

Because we have ruled out very large negative 21 cm brightness temperatures arising from cold reionization, we can significantly constrain the range of possible values of the 21 cm power spectrum in our model. Figure 30 shows the inferred 68% and 95% posterior equal-tailed credible intervals on the power spectrum at redshifts 12, 10, 8, and 6 after incorporating results of a full season of HERA observation. One can think of these contours as the testable—and therefore falsifiable—predictions of 21cmMC given the HERA data, the other astrophysical constraints, and the assumptions of the model. In Figure 30, we also show the 95% equal-tailed credible interval from H22b and from our posterior without HERA. Upper limits from H22a and this work as well as a number of other previously published measurements from GMRT (Paciga et al. 2013), PAPER (Kolopanis et al. 2019), MWA (Li et al. 2019; Trott et al. 2020), and LOFAR (Mertens et al. 2020) are also presented in Figure 30.

Zoom In Zoom Out Reset image size

It is evident that the HERA limits have been significantly improved since H22a and that this work represents our best constraint so far on the neutral IGM during the EoR. As a result, the posterior distribution of 21 cm power spectra in our model has also become much tighter, which is consistent with our exclusion of an unheated IGM at these redshifts. That said, one might wonder why the 95% posteriors are well below our 2σ power spectrum upper limits in Figure 30. The reason is threefold. First, recall that our likelihood model is not Gaussian; it is an error function that asymptotes to a flat probability at small model Δ² (Equation (11)). As Figure 26 shows, power spectra just below our limits are far less likely than power spectra well below them. This concentrates posterior probability at lower power spectrum values. The effect is especially important when multiple fields and k modes contribute significantly to the likelihood, instead of the two measurements from a single field that dominated the results in H22b.

The second reason is that, as already discussed, our inference is heavily informed by other high-redshift observations and the galaxy model assumed by 21cmFAST. This is particularly relevant in the context of the third reason, which is the influence of constraints from the two bands on each other. To better understand the impact of the model, we show in the top row of Figure 31 the prior and full posterior probability distributions from the three inferences, i.e., without HERA and with HERA after 18 (H22b) and 94 nights (this work). We only show the distributions at the k values of our deepest limits, roughly 0.35 h Mpc⁻¹. First consider the PDF without incorporating HERA. It has two clear peaks: the one at smaller Δ² corresponds to models with ${\overline{T}}_{S}\gg {T}_{\mathrm{CMB}}$ (with abundant X-ray heating) while the other is mostly "cold reionization" with little heating. In between, there is a valley in the distribution because one must fine-tune the heating and ionization to get a signal between these two extremes, which is comparatively unlikely given our priors and the other high-redshift constraints.

Zoom In Zoom Out Reset image size

The H22b results ruled out the extreme end of the cold reionization peak z = 7.9, but at z = 10.4 much of that stronger peak was still viable. With our new limits, an IGM near the adiabatic limit is essentially excluded at z = 10.4, which clearly favors models with considerable X-ray heating. Because that heating is assumed by the models to continue through z = 7.9 (when the stellar mass density, and hence density of X-ray sources, has increased by a factor of a few), the higher-redshift measurement helps to constrain the lower-redshift one as well. In particular, it eliminates the last little tail of large amplitude power spectra in the posterior, which is why the 95% contour moves so much between H22b and this work in Figure 30, though the 68% contour (not shown here for H22b) does not change very much.

As further evidence for the impact of the z = 10.4 measurements on the z = 7.9 posterior, we show two additional sets of curves in the bottom row of Figure 31. The first are the ratios of the posterior after 94 nights with HERA to the pre-HERA posterior. The second are the likelihoods that we get from Equation (11) by combining together measurements from all fields and k modes from 0.2 < k < 1.0 h Mpc⁻¹ (dropping every other k, as discussed in Section 7.1). The likelihoods are evaluated assuming power law Δ² ∝ k^1.3 at z = 10.4 and Δ² ∝ k^0.0 at z = 7.9 (see Section 7.1 for details). Both curves are normalized to plateau at 1. Since the pre-HERA posterior is treated as the prior for the post-HERA inference, we should expect this ratio to match the likelihood by Bayes' theorem. It should be noted that using a power-law power spectrum is only an approximation. To understand the precise disagreements between the likelihood and the ratio of the posteriors, we would have to account for the detailed dependence of Δ² on k, the ways in which the z = 7.9 measurement constrains the z = 10.4 posterior (and not just vice versa), and sampling noise.

With those caveats, the match looks reasonable at z = 10.4, but not nearly as good at z = 7.9. It follows then that the inferred constraints on both ${\overline{T}}_{S}$ and Δ² at z = 7.9 are driven externally to the data at z = 7.9, which is all that goes into the black likelihood curve. Therefore, the z = 7.9 results must be a consequence of the z = 10.4 limits and the way 21cmFAST models the evolution the X-ray luminosity by tying it to star formation and then constraining that star formation rate with other probes, most notably the UV luminosity function. In Section 7.6, we will return to the question of the impact of the specific modeling choices of 21cmMC and how they compare to the other three techniques.

In H22b, we reported parameter constraints on an alternate seminumerical model that allows for a significant excess radio background. EoR scenarios with high levels of radio background at rest-frame 21 cm wavelengths can potentially produce strong 21 cm absorption signals, since the brightness temperature is proportional to (1 − T_radio/T_S) (Feng & Holder 2018). We know there exists today a radio background well in excess of the CMB from observations with ARCADE 2 (Fixsen et al. 2011; Seiffert et al. 2011) and the LWA (Dowell & Taylor 2018). It has been theorized that if this excess is sourced by a population of unresolved high-redshift, potentially exotic sources (Ewall-Wice et al. 2018; Fraser et al. 2018; Pospelov et al. 2018; Brandenberger et al. 2019; Jana et al. 2019; Thériault et al. 2021), it could explain the anomalously strong absorption signal seen by EDGES (Bowman et al. 2018). Such explanations are not without difficulty; they would have to feature sources far stronger than what would be expected from local observations (Ewall-Wice et al. 2018, 2020; Mirocha & Furlanetto 2019; Mebane et al. 2020) and would need to include rapid X-ray heating at 16 ≳ z ≳ 14 to explain the high-frequency side of the EDGES trough (e.g., Mittal & Kulkarni 2022).

Since a radio background can also increase the amplitude of 21 cm fluctuations (Ewall-Wice et al. 2014; Fialkov & Barkana 2019; Reis et al. 2020), limits from HERA can constrain astrophysical parameters describing models with excess radio background. In H22b, we used a seminumerical simulation (Visbal et al. 2012; Fialkov et al. 2014; Fialkov & Barkana 2019; Cohen et al. 2020; Reis et al. 2020, 2021) in which the key radiation fields are all driven by the cosmic star formation rate. This is set (in part) by the star formation efficiency f_* with which collapsed gas in halos turns into stars, and the circular velocity V_c , which determines the minimum mass for star-forming halos. Just as in H22b, the efficiencies of X-ray and radio background luminosity relative to the star formation rate are parameterized by f_X and f_r , and reionization is parameterized by the CMB optical depth τ. For more details on the implementation of the model and its parameterization, see Section 8 of H22b. ⁵⁸

There are a few differences between the precise inference procedure used here and the one used in H22b. More parameters are now allowed to vary rather than fixing them to specific values. Most significantly, we no longer limit ourselves to the X-ray spectral energy distribution (SED) of X-ray binaries (Fragos et al. 2013; Fialkov et al. 2014) and instead parameterize the SED by a power-law spectral index α_X (either 1.0, 1.3, or 1.5) and a minimum frequency cutoff ${\nu }_{\min }$, which we vary from 0.1–3 keV with log-uniform prior. Additionally, we no longer fix the mean free path to R_mfp = 40 cMpc but vary it between 10 and 70 cMpc. ⁵⁹ The new parameterization of the X-Ray SED implies that more varieties of X-ray heating are captured. For the same value of f_X , a softer SED (with a larger number of low-energy photons, e.g., lower ${\nu }_{\min }$) leads to a more efficient heating with shorter X-ray mean free path, while a harder SED (lower number of low-energy photons, e.g., higher values of ${\nu }_{\min }$) results in less heating with a larger fraction of photons remaining unabsorbed.

Another improvement over the analysis performed in H22b is the application of more accurate emulators. Our new emulator, based on the GLOBALEMU formalism (Bevins et al.2021) and implemented using scikit-learn (Pedregosa et al. 2011), uses a more general neural network reaching an accuracy of ${}_{-7 \% }^{+11 \% }$ (at 68% confidence) at z = 7.9 and k = 0.34 h Mpc⁻¹. Additionally, instead of an MCMC, we now use nested sampling with PolyChord (Handley et al. 2015a, 2015b).

The updated analysis and the improved HERA upper limit together result in tightened constraints, two of which are of particular interest. The biggest impact is on the X-ray efficiency f_X , which is consistent with the 21cmMC results on L_X /SFR (see Section 7.4.1). In general, HERA excludes models with high f_r and low f_X , since they would produce the brightest amplitude of 21 cm fluctuations. The lower bound on the 68% HPD credible interval of f_X increases from f_X > 0.03, using the data from H22a, to f_X > 0.18 using the limits presented here. However, the upper bound on the 68% HPD credible interval of f_r decreases only slightly, from f_r < 586 to f_r < 575. Note that these numbers differ from the values quoted in H22b due to the changed X-ray SED and priors.

In Figure 32, we show the region of the parameter space disfavored by HERA in this model, where we include the impact of f_* on both X-ray heating and the total radio background. This space was already constrained by other measurements of astrophysical backgrounds. Reis et al. (2020) showed that models with strong f_* · f_r are inconsistent with the radio background from LWA (Dowell & Taylor 2018) and ARCADE 2 (Fixsen et al. 2011; Seiffert et al. 2011), where the lower-frequency measurements from LWA are the more constraining of the two (Reis et al. 2020). Likewise, Fialkov et al. (2017) showed that large values of f_* · f_X are generally ruled out by Chandra X-ray background measurements in the 0.5–2 keV band (Lehmer et al. 2012). We show those limits in Figure 32 in yellow and blue, respectively. Compared to the results in H22b, our new 21 cm constraints leave little room for the LWA radio background to be explained entirely as a cosmological signal originating from z > 8. Such models are not entirely excluded, but they are mostly disfavored at 68% credibility or greater.

Zoom In Zoom Out Reset image size

The second result, which is also an effect of the higher X-ray efficiency, can be seen via our derived constraints on the spin temperature of the IGM. In Figure 33, we show how HERA data updates our model's prior on ${\overline{T}}_{S}$ and ${\overline{T}}_{\mathrm{radio}}$. Both of these quantities are averaged over the neutral IGM at z = 7.9.

Zoom In Zoom Out Reset image size

In general, our results exclude models with large radio background temperatures without also featuring large spin temperatures. Specifically, our 95% credible interval requires that ${\overline{T}}_{S}/{\overline{T}}_{\mathrm{radio}}\gt 0.21$ at z = 7.9. By performing the same calculation for Band 1, we can also show that at z = 10.4, ${\overline{T}}_{S}/{\overline{T}}_{\mathrm{radio}}\gt 0.18$ with 95% credibility. While the posterior largely prefers relatively low radio background temperatures, there remain a set of models where a large radio background becomes coupled to the spin temperature, yielding ${\overline{T}}_{S}/{\overline{T}}_{\mathrm{radio}}\approx 1$ and thus small power spectra.

Finally, we can now more fully compare the results of our four models and their implications for the thermal history of the IGM. As we saw in Figure 27, the results from 21cmMC exhibit a more dramatic improvement at z = 7.9 in this work relative to H22b than either of the phenomenological models (see Section 7.3) or the radio background model (see Section 7.5). In the discussion around Figure 31, we saw clear evidence that 21cmMC's change in ${\overline{T}}_{S}$ and thus in the inferred posterior probability distribution of Δ² was driven by the power spectrum constraints at z = 10.4. This tight connection between redshifts requires rather rapid thermal evolution of the IGM, which is generally the case for the range of 21cmMC models we consider.

In the other models, the spin temperature history can evolve quite gradually, and thus yield similarly gradual 21 cm signal evolution that will be more difficult to disfavor with improved z = 10.4 limits. For the reionization-driven model, this is because ${\overline{T}}_{S}$ is parameterized as a power law in redshift with a uniform prior on the power-law index. In the radio background model, there are two reasons for gradual evolution. First, the stellar-mass–halo-mass relation is constant in the radio background model, which gives rise to a much more gradual cosmic star formation history than, e.g., models in which the relation is itself a function of halo mass. Second, because the radio and X-ray backgrounds are both sourced by galaxies, both ${\overline{T}}_{\mathrm{radio}}$ and ${\overline{T}}_{S}$ will grow similarly as the cosmic star formation rate density rises, resulting in relatively gradual evolution in the 21 cm signal. By contrast, in our 21cmMC simulations, T_radio = T_CMB, which of course declines as (1 + z). This means ${T}_{\mathrm{radio}}/{\overline{T}}_{S}$ will evolve more rapidly than in than in the radio background models, even if both share an identical ${\overline{T}}_{S}$ history.

Another reason for the difference between models is the way 21cmMC directly incorporates data from other wavelengths (especially the galaxy luminosity function)—a major advantage of the approach. Models that fit the luminosity function at all redshifts require the star formation efficiency to increase with galaxy mass (at least up to galaxy masses comparable to the Milky Way, see, e.g., Mirocha et al. 2017; Park et al. 2019). The 21cmMC models therefore also favor massive galaxies, whose abundance evolves quite rapidly at high redshifts. Because these galaxies are also the sources of IGM heating in this model, their rapidly evolving abundance also leads to rapid IGM heating. Because the luminosity function fits "bake in" this kind of behavior, 21cmMC models that feature modest heating at z = 10.4 require stronger heating at z = 7.9. The expected increase is roughly the ratio of the stellar mass density at the two redshifts, which is a factor of ∼3. This explains why the strong ${\overline{T}}_{S}$ limit in 21cmMC at z = 7.9 is actually informed in large part by the ostensibly weaker temperature limits at z = 10.4.

Of course, this conclusion is reliant on the parameterization of the galaxy luminosity function within 21cmMC. While a wide range of models agree that a physically motivated extrapolation of the observed luminosity function to higher redshifts behaves similarly, it is of course reasonable to suppose that the extrapolation breaks down—perhaps because of processes unique to high redshifts, or perhaps simply because young dwarf galaxies behave differently than their larger cousins (see, e.g., Qin et al. 2020). One could always broaden the allowed set of galaxy models to accommodate such exotic astrophysics. We have, of course, explored one class of such models by examining the effect of an excess radio background in Section 7.5. One could also imagine adding another source of early heating like Population III remnants that could inject sufficient heat into the IGM at z = 10.4 without necessarily affecting the constraints at z = 7.9, though this might require some degree of fine tuning—a question left for future work. Regardless, it is useful to understand how we should interpret the HERA limits in the context of "vanilla" galaxy models. Here, the 21cmMC results demonstrate that the "known" galaxy population requires substantial IGM heating by z ∼ 8.

Regardless of the differences between our models, we should not lose sight of the fact that they are all pointing in fundamentally the same direction: a heated IGM during the EoR. To put them in the broader cosmological context of other 21 cm experiments, we compare them to a few illustrative scenarios for the thermal history of the IGM in Figure 34. The arrows show the lower bounds of the 95% credible regions of ${\overline{T}}_{S}/{T}_{\mathrm{radio}}$ at z = 7.9 and z = 10.4 inferred by each of the four models. We also show the expectation for an IGM temperature set by adiabatic cooling since recombination with no additional heating from astrophysical sources. Just as Figure 27 showed, the lower limits at both redshifts are above this curve, requiring substantial X-ray heating (as illustrated by the solid black line where L_{X<2 keV}/SFR = 10⁴¹ erg s⁻¹ ${M}_{\odot }^{-1}$ yr). This is a significant improvement from H22b, which could not demonstrate any heating beyond the adiabatic cooling curve at z = 10.4. More interestingly, they are all inconsistent with an IGM thermal history that one would predict by extrapolating the local relation between X-ray luminosity and star formation rate to high redshift (black dotted–dashed line; see Figure 28). The accordance of multiple models here strengthens the argument in Section 7.4.1 that our results broadly favor an IGM heated by low-metallicity HMXBs.

Zoom In Zoom Out Reset image size

Figure 34 also shows measurements at higher redshifts. The EDGES result (Bowman et al. 2018) found a surprisingly large absorption signal at z ∼ 17, requiring either substantial cooling of the gas below the adiabatic limit or a larger radio background than provided by the CMB. The implied range of ${\overline{T}}_{{\rm{S}}}/{T}_{\mathrm{radio}}$ is shown by the red error bars, but note that EDGES found this absorption trough to be very narrow, implying that the universe was heated shortly afterward. However, the recent SARAS 3 measurement (Singh et al. 2022) did not detect such a signal; the upper limit is subtle to express quantitatively, but here we show the limit their measurement places on a feature with the shape of the best-fit EDGES signal with an unknown amplitude. Though we show the SARAS 3 limit offset in redshift for clarity, the measurement spanned a wide frequency range corresponding to 15.8 < z < 24.6.

Because HERA observed at much smaller redshifts than EDGES and SARAS 3, it is difficult to compare directly—even if the deep EDGES trough is real, any heating between z ∼ 15 and z ∼ 10 (as indeed the EDGES best-fit model requires!) would make the two measurements consistent. We show an example of such a model (green solid line) where the IGM is cooled by interactions with a fraction (0.5%) of millicharged dark matter with mass 10 MeV and charge 10⁻⁵ times that of the electron, following Muñoz & Loeb (2018). The gas is subsequently heated by HMXBs. This model is designed to explain EDGES and still be consistent with lower-redshift observations like HERA's, though it is in mild tension with SARAS 3. More independent measurements like SARAS 3, as well as low-frequency power spectra from HERA Phase II and other interferometers, will be extremely valuable in understanding this early era and any new physics that may have impacted it.