Results from EDGES High-band. II. Constraints on Parameters of Early Galaxies

Here, we describe the computation of the EDGES data likelihood for the 21 cm signals associated to each astrophysical model. We start by modeling the uncertainty in the EDGES integrated spectrum as Gaussian and, thus, write the likelihood of the data as

$$\begin{eqnarray}{ \mathcal L }(d| {\boldsymbol{\theta }}) & = & \displaystyle \frac{1}{\sqrt{{(2\pi )}^{{N}_{\nu }}| {\rm{\Sigma }}| }}\\ & & \times \exp \left\{-\displaystyle \frac{1}{2}{[d-m({\boldsymbol{\theta }})]}^{T}{{\rm{\Sigma }}}^{-1}[d-m({\boldsymbol{\theta }})]\right\},\end{eqnarray} \tag{ 1 }$$

where d is our spectrum with ${N}_{\nu }$ frequency channels, m is the model for the spectrum with parameters ${\boldsymbol{\theta }}$, and Σ is the covariance matrix of the data, of size ${N}_{\nu }\times {N}_{\nu }$. The model of the spectrum is the sum of the model for the 21 cm signal and the model for the foreground, which are functions of different parameters:

$$\begin{eqnarray}&&m({\boldsymbol{\theta }})={m}_{21}({{\boldsymbol{\theta }}}_{21})+{m}_{\mathrm{fg}}({{\boldsymbol{\theta }}}_{\mathrm{fg}}).\end{eqnarray} \tag{ 2 }$$

Here, ${{\boldsymbol{\theta }}}_{21}$ represents the four-element vector of astrophysical parameters associated with the 21 cm signal, and ${{\boldsymbol{\theta }}}_{\mathrm{fg}}$ is the vector of foreground parameters.

We model the foreground with the polynomial expression (Mozdzen et al. 2016; Monsalve et al. 2017a, 2017b)

$$\begin{eqnarray}&&{m}_{\mathrm{fg}}({{\boldsymbol{\theta }}}_{\mathrm{fg}})=\displaystyle \sum _{i=0}^{{N}_{\mathrm{fg}}-1}{a}_{i}{\nu }^{-2.5+i}=A{{\boldsymbol{\theta }}}_{\mathrm{fg}}.\end{eqnarray} \tag{ 3 }$$

The A matrix has a size ${N}_{\nu }\times {N}_{\mathrm{fg}}$, where the ${N}_{\mathrm{fg}}$ columns correspond to the ${\nu }^{-2.5+i}$ basis functions. The first function, with i = 0, primarily accounts for the contribution to the brightness temperature by Galactic synchrotron radiation. Other terms fit additional spectral structure in the sky or contributions from residual calibration effects. To determine the optimum number of terms in the foreground model over $90\mbox{--}190\,\mathrm{MHz}$, we compute the weighted residuals rms for a range of ${N}_{\mathrm{fg}}$ values. We see the rms going through a knee at ${N}_{\mathrm{fg}}=5$. Specifically, for ${N}_{\mathrm{fg}}=1,2,3,4$, the rms are 7100, 758, 200, and 57 mK, respectively. At ${N}_{\mathrm{fg}}=5$ the rms is 17 mK, and adding more terms decreases the rms by $\lt 1$ mK. Thus, in this analysis we use ${N}_{\mathrm{fg}}=5$.

Other models for the diffuse foreground have been proposed. The most popular is the $\mathrm{log}{(T)=\sum {a}_{i}(\mathrm{log}\nu )}^{i}$ expression, which also tries to capture the temperature, spectral index, and departures of the spectrum from a power law (Pritchard & Loeb 2010; Harker et al. 2012, 2016; Voytek et al. 2014; Bernardi et al. 2015, 2016; Harker 2015). Non-parametric alternatives have also been suggested (Switzer & Liu 2014; Vedantham et al. 2014), some of which rely on training sets from simulated observations using instrument models and all-sky radio maps (Tauscher et al. 2018). Unfortunately, the absolute accuracy of the models and maps is not yet high enough to reach agreement between simulations and measurements at the desired level of a few tens of mK or better (i.e., Haslam et al. 1982; de Oliveira-Costa et al. 2008; Guzmán et al. 2011; Sathyanarayana Rao et al. 2017a; Zheng et al. 2017). Although we considered using other models, including the linearized physical foreground model introduced in Bowman et al. (2018), whose basis functions approximate the expected contributions from Galactic synchrotron radiation and the ionosphere, we ultimately chose to use Equation (3). This is a generic expression that with few terms can model our spectrum, which is potentially subject to low-level calibration residuals (Mozdzen et al. 2016; Monsalve et al. 2017a). We note that although Singh et al. (2017, 2018) use "maximally smooth" functions (Sathyanarayana Rao et al. 2015, 2017b) in laboratory tests and simulations, during the analysis of sky data they also rely on generic polynomials to model their instrument-filtered foreground signal.

We can evaluate Equation (3) analytically at any foreground parameter values. However, the 21 cm signals in our library are fixed. For this reason, we cannot compute Equation (1) for random values of the full parameter vector (astrophysical + foreground). Instead, to explore the astrophysical parameter space, we compute the likelihood for each 21 cm signal in our library as follows. (i) First, we subtract the 21 cm signal from the spectrum. (ii) Then, we subtract from (i) the foreground model that best fits the difference in (i). (iii) Finally, we compute the likelihood for the 21 cm signal from the residuals in (ii) after marginalizing over the uncertainty in the foreground parameters. This approach represents an alternative to the simultaneous exploration of the full parameter space with techniques such as Markov Chain Monte Carlo (MCMC). We leave the MCMC analysis for future work, as we are currently testing tools for efficient evaluation of 21 cm signals for random combinations of astrophysical parameters.

Now, we provide details of the process outlined in the previous paragraph. Taking advantage of the linearity of the foreground parameters, we compute their maximum likelihood estimator (MLE) as follows:

$$\begin{eqnarray}&&{\hat{{\boldsymbol{\theta }}}}_{\mathrm{fg}}={({A}^{T}{{\rm{\Sigma }}}^{-1}A)}^{-1}{A}^{T}{{\rm{\Sigma }}}^{-1}(d-{m}_{21}({{\boldsymbol{\theta }}}_{21})).\end{eqnarray} \tag{ 4 }$$

This MLE maximizes the likelihood along the foreground parameter axes having conditioned on a selection of astrophysical parameters, ${{\boldsymbol{\theta }}}_{21}$. The foreground parameters can now be written in terms of the MLE and a deviation term,

$$\begin{eqnarray}&&{{\boldsymbol{\theta }}}_{\mathrm{fg}}={\hat{{\boldsymbol{\theta }}}}_{\mathrm{fg}}+{{\boldsymbol{\delta }}}_{\mathrm{fg}}.\end{eqnarray} \tag{ 5 }$$

Correspondingly, the foreground model can be expanded as follows:

$$\begin{eqnarray}&&{m}_{\mathrm{fg}}({{\boldsymbol{\theta }}}_{\mathrm{fg}})={m}_{\mathrm{fg}}({\hat{{\boldsymbol{\theta }}}}_{\mathrm{fg}})+{m}_{\mathrm{fg}}({{\boldsymbol{\delta }}}_{\mathrm{fg}})=A{\hat{{\boldsymbol{\theta }}}}_{\mathrm{fg}}+A{{\boldsymbol{\delta }}}_{\mathrm{fg}}.\end{eqnarray} \tag{ 6 }$$

With this foreground model, whose only free parameters are ${{\boldsymbol{\delta }}}_{\mathrm{fg}}$, Equation (1) becomes a function of ${{\boldsymbol{\theta }}}_{21}$ and ${{\boldsymbol{\delta }}}_{\mathrm{fg}}$, and can be written as

$$\begin{eqnarray}{ \mathcal L }(d| {{\boldsymbol{\theta }}}_{21},{{\boldsymbol{\delta }}}_{\mathrm{fg}}) & = & \displaystyle \frac{1}{\sqrt{{(2\pi )}^{{N}_{\nu }}| {\rm{\Sigma }}| }}\\ & & \times \exp \left\{-\displaystyle \frac{1}{2}{({d}_{\star }-A{{\boldsymbol{\delta }}}_{\mathrm{fg}})}^{T}{{\rm{\Sigma }}}^{-1}({d}_{\star }-A{{\boldsymbol{\delta }}}_{\mathrm{fg}})\right\},\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{d}_{\star }=d-{m}_{21}({{\boldsymbol{\theta }}}_{21})-{m}_{\mathrm{fg}}({\hat{{\boldsymbol{\theta }}}}_{\mathrm{fg}}).\end{eqnarray} \tag{ 8 }$$

We show in Appendix A that the marginalization of this likelihood over the foreground parameters, ${{\boldsymbol{\delta }}}_{\mathrm{fg}}$, gives the following marginal likelihood for the 21 cm astrophysical model,

$$\begin{eqnarray}{ \mathcal L }(d| {{\boldsymbol{\theta }}}_{21}) & = & \int { \mathcal L }(d| {{\boldsymbol{\theta }}}_{21},{{\boldsymbol{\delta }}}_{\mathrm{fg}})d{{\boldsymbol{\delta }}}_{\mathrm{fg}}\\ & = & \sqrt{\displaystyle \frac{{(2\pi )}^{{N}_{\mathrm{fg}}-{N}_{\nu }}}{| {\rm{\Sigma }}| | {C}^{-1}| }}\exp \left\{-\displaystyle \frac{1}{2}{d}_{\star }^{T}{({\rm{\Sigma }}+V)}^{-1}{d}_{\star }\right\},\end{eqnarray} \tag{ 9 }$$

where $C={({A}^{T}{{\rm{\Sigma }}}^{-1}A)}^{-1}$ is the 5 × 5 covariance matrix of the foreground parameters, $V={({{\rm{\Sigma }}}_{\mathrm{fg}}^{-1}-{{\rm{\Sigma }}}^{-1})}^{-1}$, and ${{\rm{\Sigma }}}_{\mathrm{fg}}={{ACA}}^{T}$ is the ${N}_{\nu }\times {N}_{\nu }$ covariance matrix of the foreground model in the frequency domain. The dependence of ${ \mathcal L }$ on the astrophysical parameters is only through d_⋆. C and V have to be calculated only once, as they depend only on (1) the covariance matrix of the data and (2) the foreground model. This enables us to process many 21 cm signals very efficiently.

The covariance matrix of the data, Σ, is the sum of the covariance matrix of the noise and a matrix that accounts for systematic uncertainty. The noise covariance matrix is diagonal and describes the channel variance resulting from the effective system temperature and integration time. For reference, the channel standard deviation at 140 MHz with a channel width of 390.6 kHz is 6 mK.

Our systematics covariance matrix represents the variations in the spectral residuals to a five-term foreground model, for different values of the calibration parameters within their uncertainty range. To estimate this matrix, we propagate the uncertainty in the calibration parameters to the spectrum in a Monte Carlo (MC) fashion. In each MC realization, we calibrate the data by randomly perturbing all the calibration parameters simultaneously in the ranges listed in Table 1 of Monsalve et al. (2017b). Then, we compute the rms residuals across the spectrum to our five-term foreground model. Finally, we identify the 68% upper limit for the MC rms distribution and subtract the minimum rms of the distribution, which occurs for our nominal calibration. From this computation, we obtain a 68% limit of 35 mK. Due to the wide range of perturbation possibilities resulting from all uncertainty sources affecting simultaneously, we categorize the channel-to-channel systematic uncertainty as uncorrelated. Thus, we arrive at the diagonal systematics covariance matrix ${{\rm{\Sigma }}}_{\mathrm{syst}}={(35\mathrm{mK})}^{2}{{\boldsymbol{I}}}_{{N}_{\nu }}$. We are leaving for future work the estimation of a more refined covariance matrix. See Section 5.4.

Using a more restrictive foreground model that only accounts for the expected physical foreground contributions could reduce the covariance with the 21 cm signals. However, it could also lead to (1) higher residuals from the fit to the nominal spectrum if low-level calibration errors are present and (2) a higher systematic uncertainty estimate, derived as described in the previous paragraph. As an example, fitting the High-band spectrum with the linearized physically motivated foreground model of Bowman et al. (2018), Equation (1), produces higher residuals rms than Equation (3) with five terms. Thus, using such a model might ultimately result in weaker overall constraints and larger effects from fit residuals. For our model of Equation (3), these effects are discussed in Sections 5.1 and 5.2.

In our second analysis scenario, we apply three observational constraints in addition to the EDGES data. They correspond to

• 1.  
  The constraint on the electron-scattering optical depth from Planck, ${\tau }_{e}=0.058\pm 0.012$ (Planck Collaboration XLVII 2016). We model this constraint as Gaussian with center and standard deviation as provided by that estimate.
• 2.  
  A constraint on ${\bar{x}}_{{\rm{H}}{\rm{I}}}$ from the "dark fraction," i.e., the fraction of pixels dark in both Lyα and Lyβ, in the spectra of high-redshift quasars. This constraint is independent from reionization modeling. McGreer et al. (2015) recently obtained a 1σ upper limit of ${\bar{x}}_{{\rm{H}}{\rm{I}}}\lesssim 0.06+0.05$ at z = 5.9. We model this limit as a uniform probability distribution over ${\bar{x}}_{{\rm{H}}{\rm{I}}}\lesssim 0.06$ and a one-sided Gaussian, with σ = 0.05 for ${\bar{x}}_{{\rm{H}}{\rm{I}}}\gt 0.06$.
• 3.  
  The first¹⁰ detection (2σ) of an ongoing reionization, obtained from a Lyα damping wing: the ${\bar{x}}_{{\rm{H}}{\rm{I}}}$ PDF from the analysis of the z = 7.1 quasar ULASJ1120+0641 (Mortlock et al. 2011; Greig et al. 2017) with the "Small H ii" model¹¹ of reionization.

In this combined analysis scenario, the likelihood for each astrophysical model is obtained by multiplying the likelihood of the EDGES data for the 21 cm signal (Section 4.1), by the likelihood of τ_e, ${\bar{x}}_{{\rm{H}}{\rm{I}}}$ at z = 5.9, and ${\bar{x}}_{{\rm{H}}{\rm{I}}}$ at z = 7.1, for the same model, computed from the constraints in points 1–3 above.

We interpolate the likelihood computed for our sample of astrophysical parameter values onto a grid of 40 equally spaced bins along each of the four dimensions. The interpolation is performed using an inverse distance weighting scheme that uses the simple Shepard's method, where we weigh the 4D distance between the interpolated position and the values originally sampled, r, as $w={r}^{-12}$. The exponent of −12 was selected to minimize numerical artifacts in our interpolated probability space. We tested values −4, −8, and −12, and verified that the choice does not alter the shapes of the PDFs or the marginalized 68% and 95% probability limits.

We then assume uniform priors for the astrophysical parameters (uniform in log₁₀ scale for ${T}_{\mathrm{vir}}^{\min }$ and ${L}_{{\rm{X}}\lt 2\mathrm{keV}}/\mathrm{SFR}$) and numerically integrate the gridded 4D likelihood to obtain marginalized 1D and 2D posterior PDFs for each parameter and parameter pair. We smooth these PDFs (both 1D and 2D) with a Gaussian kernel of scale 1.5 in order to minimize small-scale numerical artifacts that resulted from the interpolation procedure.

In Figure 2, we see that the EDGES High-band spectrum provides significant discrimination across the parameter space. The most interesting constraints are those for ${T}_{\mathrm{vir}}^{\min }$ and ${L}_{{\rm{X}}\lt 2\mathrm{keV}}/\mathrm{SFR}$ (second and third parameters), including their joint PDF, because changes in the value of these parameters have the largest impact on the global signal.

Zoom In Zoom Out Reset image size

${T}_{\mathrm{vir}}^{\min }$ regulates the minimum mass of halos that produce star-forming galaxies, which strongly affects the timing of the IGM heating and reionization. Lower values of ${T}_{\mathrm{vir}}^{\min }$ result in spectral features produced earlier, with a significant fraction of the absorption trough occurring at frequencies below the EDGES High-band range. Higher values produce troughs to which we have higher sensitivity. This leads to the expectation of lower probabilities as ${T}_{\mathrm{vir}}^{\min }$ increases.

${L}_{{\rm{X}}\lt 2\mathrm{keV}}/\mathrm{SFR}$ controls the IGM heating, therefore it has a significant impact on the timing, depth, and width of the absorption trough. For this parameter, we expect our highest sensitivity to be below the middle of the range, where the trough is relatively deep and narrow. For high values, the troughs are narrow but shallow, while low values result in late, deep, and wide troughs, as the heating is inefficient and the IGM temperature continues its descending trend before it slowly reacts to the heating.

We see in the ${T}_{\mathrm{vir}}^{\min }$ and ${L}_{{\rm{X}}\lt 2\mathrm{keV}}/\mathrm{SFR}$ panels of Figure 2 that the PDFs are broadly consistent with expectations. The data disfavor high values of ${T}_{\mathrm{vir}}^{\min }$ and intermediate values of ${L}_{{\rm{X}}\lt 2\mathrm{keV}}/\mathrm{SFR}$. The high-probability region of low ${T}_{\mathrm{vir}}^{\min }$ and high ${L}_{{\rm{X}}\lt 2\mathrm{keV}}/\mathrm{SFR}$, condensed in the upper-left corner of the panel ${T}_{\mathrm{vir}}^{\min }\mbox{--}{L}_{{\rm{X}}\lt 2\mathrm{keV}}/\mathrm{SFR}$, corresponds to signals with most of the absorption trough below the EDGES High-band range. The models in the lower-right corner of the same panel (${\mathrm{log}}_{10}({T}_{\mathrm{vir}}^{\min })\gtrsim 5.7$ and ${\mathrm{log}}_{10}({L}_{{\rm{X}}\lt 2\mathrm{keV}}/\mathrm{SFR})\lesssim 39$) have high probabilities because although they represent "cold" EoR scenarios with deep absorption troughs, the troughs fall close to the high-frequency limit of the EDGES band, where we have low sensitivity.

To develop a better intuition for the structure in the PDFs of the astrophysical parameters, especially ${T}_{\mathrm{vir}}^{\min }$ and ${L}_{{\rm{X}}\lt 2\mathrm{keV}}/\mathrm{SFR}$, we (1) correlated the probability of the 21 cm signals with the phenomenological parameters of their absorption troughs (center, depth, and width) and (2) computed the PDFs for a variety of simulated EDGES High-band spectra. Specifically, through point (2), we studied the impact of the residuals of the measured spectrum to a five-term foreground model fit. These residuals are shown in Figure 4(b) of Monsalve et al. (2017b), and can be described as ripples across the spectrum with a period of $\approx 20\,\mathrm{MHz}$, an amplitude that decreases with frequency, and a weighted rms over frequency of 17 mK. In these tests, the simulated spectra consisted of the best-fit foreground model plus noise and ripples that match the real data. In some tests, we varied the amplitude of the ripples, and in others we added ripples only within specific sub-bands of the spectrum. The simulations also included the systematic uncertainty assigned to the real data. For reference, in Appendix B, Figure 5, we show the PDFs resulting from a simulated spectrum that only includes the foreground model and noise. Comparing that figure with Figure 2 illustrates the impact of the ripples on our results. From our simulations, and referring to Figures 2 and 5, we derive the following main conclusions:

• 1.  
  The high-probability region of low ${T}_{\mathrm{vir}}^{\min }$ and high ${L}_{{\rm{X}}\lt 2\mathrm{keV}}/\mathrm{SFR}$ is not sensitive to the structure in the spectrum above the smooth foreground model.
• 2.  
  The high-probability blob in the upper-right corner of the panel ${T}_{\mathrm{vir}}^{\min }-{L}_{{\rm{X}}\lt 2\mathrm{keV}}/\mathrm{SFR}$ (including its lower-probability tail along the diagonal), is produced by ripples in the measured spectrum below ∼110 MHz. Most signals corresponding to this blob have shallow ($\lesssim 100$ mK) absorption troughs centered around that frequency.
• 3.  
  Ripples in the range ∼110–140 MHz are responsible for the high-probability intersection of $5.0\lesssim {\mathrm{log}}_{10}({T}_{\mathrm{vir}}^{\min })\lesssim 5.5$, ${\mathrm{log}}_{10}({L}_{{\rm{X}}\lt 2\mathrm{keV}}/\mathrm{SFR})\leqslant 39.5$, and ${E}_{0}\gtrsim 0.6$. The corresponding signals have deep ($\gtrsim 200$ mK), but wide ($\mathrm{FWHM}\,\gtrsim 50$ MHz) troughs.

Parameters ζ and E₀ have a lower impact on the global 21 cm signal, but the EDGES data can still discriminate across the value ranges explored. ζ represents the ionizing efficiency and thus has control over the timing and the sharpness of the signal throughout reionization, both while in absorption (reionization happening in parallel to X-ray heating) and emission. As the value of this parameter increases, the reionization transition is sharper, therefore higher values are disfavored by the data. E₀ represents the low-energy limit of the heating X-ray spectrum. It affects the signal in a way related to, but subtler than, ${L}_{{\rm{X}}\lt 2\mathrm{keV}}/\mathrm{SFR}$ (Greig & Mesinger 2017b). Lower values of this parameter tend to produce sharper signals and are disfavored by the data.

The broad high-probability patterns in Figure 2 are not likely to change significantly with comparable data sets that cover the same frequency range unless there is a significant reduction in noise level and uncertainty assumed. On the other hand, the high-probability blobs described above for ${T}_{\mathrm{vir}}^{\min }$ and ${L}_{{\rm{X}}\lt 2\mathrm{keV}}/\mathrm{SFR}$ could be reduced in the future with a smoother spectrum, assuming that the ripples in the current EDGES High-band data are due to instrumental artifacts or an incomplete foreground model.

In the upper half of Table 1, we present the 68% and 95% marginalized limits for the four astrophysical parameters from EDGES data alone. As a reference result, we reject $39.4\,\lt {\mathrm{log}}_{10}({L}_{{\rm{X}}\lt 2\mathrm{keV}}/\mathrm{SFR})\lt 39.8$ at 95% confidence.

Figure 3 builds on the findings of Greig & Mesinger (2017a), who derived constraints on the parameters that are most sensitive to the reionization history, ζ, ${R}_{\mathrm{mfp}}$, and ${T}_{\mathrm{vir}}^{\min }$, using constraints on τ_e from Planck and ${\bar{x}}_{{\rm{H}}{\rm{I}}}$ from high-z quasars. We point the reader to that work for details and context, in particular their Figures 5 and 8, which are the most relevant for comparison with our results. In our Figure 3, we show results after incorporating the same τ_e and ${\bar{x}}_{{\rm{H}}{\rm{I}}}$ constraints to our EDGES analysis. As mentioned in Section 3, due to the low sensitivity of the global signal to changes in ${R}_{\mathrm{mfp}}$, we do not include this parameter in our analysis, but instead we explore ${L}_{{\rm{X}}\lt 2\mathrm{keV}}/\mathrm{SFR}$ and E₀, which characterize the IGM X-ray heating.

Zoom In Zoom Out Reset image size

Out of the four parameters, our interest here focuses on ζ and ${T}_{\mathrm{vir}}^{\min }$. Both regulate the timing of reionization, thus they suffer from a degeneracy when constrained using just τ_e and ${\bar{x}}_{{\rm{H}}{\rm{I}}}$. As shown in Greig & Mesinger (2017a), these constraints reduce the probability of high ζ and, even more strongly, the probability of low ${T}_{\mathrm{vir}}^{\min }$. We see in Figures 2 and 3 that the EDGES measurement contributes to the combined result by disfavoring high values of both, ζ and ${T}_{\mathrm{vir}}^{\min }$. As a result, the joint PDF for these parameters in Figure 3 (second row, first column) has a narrower high-probability region than in Greig & Mesinger (2017a). In particular, we estimate $4.5\,\leqslant {\mathrm{log}}_{10}({T}_{\mathrm{vir}}^{\min }/{\rm{K}})\leqslant 5.7$ at 95% confidence.

The combined analysis disfavors low values of E₀ more strongly than EDGES alone. For ${L}_{{\rm{X}}\lt 2\mathrm{keV}}/\mathrm{SFR}$, this analysis still disfavors values in the middle of the range, but the probability of the lowest and especially the highest values of this parameter is reduced relative to Figure 2.

From the simulations described in Section 5.1, we also identify the effects on Figure 3 of the ripples in the spectrum above the foreground model. Specifically, reducing the ripple amplitude at ∼110 MHz shifts the high-probability peak at ${\mathrm{log}}_{10}({L}_{{\rm{X}}\lt 2\mathrm{keV}}/\mathrm{SFR})\approx 41.5$ to slightly higher values (favoring signals with the absorption trough at lower frequencies), while reducing the ripples at ∼130 MHz reduces the probability at ${\mathrm{log}}_{10}({L}_{{\rm{X}}\lt 2\mathrm{keV}}/\mathrm{SFR})\lesssim 39$. All of this also results in lower probability for high values of ζ and ${T}_{\mathrm{vir}}^{\min }$, reducing their covariance.

In the lower half of Table 1, we present the 68% and 95% marginalized limits from the combined analysis.

The combined analysis assigns low probability to a large fraction of global 21 cm signals. Most of the signals that remain with high probability can be identified with one of the two high-probability regions in the PDF of ${L}_{{\rm{X}}\lt 2\mathrm{keV}}/\mathrm{SFR}$. This results in two main "bands" of signals that remain favored. We see hints of this in Figure 1(b), where the majority of the 50 signals have very low probability (yellow) and the favored signals (green and blue) either have most features at low frequencies, corresponding to high ${L}_{{\rm{X}}\lt 2\mathrm{keV}}/\mathrm{SFR}$, or have troughs that are deep and wide within the EDGES High-band range, corresponding to low ${L}_{{\rm{X}}\lt 2\mathrm{keV}}/\mathrm{SFR}$.

We present a more complete picture of this in Figure 4, where we plot the 10,000 signals used in this analysis and highlight those that have the highest 5% probability (500 signals). These high-probability signals have a brightness temperature within ±25 mK at $z\approx 6$ and, although they span a wide range of features, with absorption peaks between ∼70 MHz and ∼130 MHz, we can identify two prominent high-probability bands: signals with troughs centered below and above ∼110 MHz. Most of these signals are smooth over $90\mbox{--}190\,\mathrm{MHz}$ and can be well fitted by our foreground model to within the uncertainties. However, some signals with peaks at ∼105 and ∼130 MHz are in the top 5% because they match ripples in the spectrum above the foreground model close to those frequencies. As discussed in Sections 5.1 and 5.2, we identified the effect of this spectral structure through simulations. When this structure is removed in simulated data, the signals with peaks at ∼105 and ∼130 MHz are replaced in the top 5% by signals with peaks at $\lesssim 90\,\mathrm{MHz}$. Despite this change in some specific signals, the two-band pattern is also present in these simulations of spectrally smooth data, which indicates that this is a robust result from the combined analysis.

Zoom In Zoom Out Reset image size

In light of the claim by Bowman et al. (2018) of an absorption feature centered at 78 MHz, it is interesting to note that several high-probability signals in Figure 4 have troughs at $\sim 70\mbox{--}90\,\mathrm{MHz}$. Regardless of the non-standard depth of the feature detected, having an absorption trough in that range would increase the probability of those signals. However, those signals have ${T}_{\mathrm{vir}}^{\min }\sim {10}^{4}$, which is at the corner of our astrophysical parameter space. In our combined analysis that only uses the High-band spectrum, those values for ${T}_{\mathrm{vir}}^{\min }$ are disfavored at $\gt 95 \%$ confidence after marginalization of the posterior over the parameter space (see Figure 3 and Table 1).

In this paper, we have improved on Monsalve et al. (2017b) by incorporating the systematic uncertainty in the spectrum into the likelihood computation. However, although the adopted uncertainty level is realistic, its covariance matrix is still generic in that it treats all spectral variations possible within that envelope as equally likely. To refine the systematics covariance matrix, we plan to conduct a more rigorous characterization of possible spectral perturbation modes of the instrument. This characterization could be joined with analysis techniques that model the instrument passband in terms of orthogonal modes derived from perturbation training sets (Burns et al. 2017; Tauscher et al. 2018).

We also leave for future work the evaluation of 21cmFAST and other 21 cm models—especially those proposed recently to account for the large absorption feature at 78 MHz—simultaneously using data from the EDGES Low- and High-band instruments, which would provide a total frequency coverage of $50\mbox{--}190\,\mathrm{MHz}$.