The Redshifted 21 cm Signal in the EDGES Low-band Spectrum

During the dark ages and the subsequent cosmic dawn (CD), the differential cooling of matter against the cosmic microwave background (CMB) is expected to give rise to redshifted global 21 cm absorption features from neutral hydrogen at those times. The precise profile of the absorption is strongly tied to a number of astrophysical parameters (Cohen et al. 2017; Mirocha et al. 2017). The first ever detection of such an absorption feature was recently reported by the EDGES experiment (Bowman et al. 2018a, hereafter BRM18). The profile of the absorption feature was centered at 78 MHz and appeared flattened. The FWHM was 19 MHz, and the depth of the absorption appeared to be 0.5 K.

The reported detection was a complete surprise: both the amplitude and the flattening of the profile were unexpected and inconsistent with previous theoretical work based on astrophysical cosmology constrained by other data (Cohen et al. 2017). Therefore, the reported signal, if confirmed and genuine, has wide-ranging implications for galaxy formation models (Fialkov & Barkana 2019; Mirocha & Furlanetto 2019), dark matter models (Lidz & Hui 2018; Liu & Slatyer 2018; Mitridate & Podo 2018; Safarzadeh et al. 2018; Schneider 2018), and also for expectations for the 21 cm power spectrum from the CD (Fialkov et al. 2018; Kaurov et al. 2018; Muñoz et al. 2018).

The amplitude is a factor of 2 larger than the maximum amplitude considered possible based on standard theoretical predictions and cannot be explained by astrophysical models in the ΛCDM framework (Barkana 2018; Witte et al. 2018) or without finely tuned modifications to the background cosmology (Hill & Baxter 2018). Physical explanations for the depth of the absorption require invoking mechanisms that might cool the gas to temperatures below that achievable by adiabatic cooling (Dvorkin et al. 2014; Tashiro et al. 2014; Muñoz et al. 2015; Barkana 2018), or else by invoking an additional radio background apart from the CMB at the CD (Slatyer & Wu 2018; Feng & Holder 2018; BRM18).

In recent literature motivated by the reported detection, the interaction between baryonic and dark matter has been considered as a new mechanism that might supplement adiabatic expansion and enhance the cooling of matter. Nonstandard Coulomb-like interactions between dark matter particles and baryons can lead to energy and momentum exchange between them and can cool the gas substantially (Barkana 2018; Muñoz & Loeb 2018; Slatyer & Wu 2018). However, such interactions are severely constrained by other astronomical observations and laboratory experiments, which indicate that the dominant dark matter component is unlikely to be able to cool the gas to the extent indicated by the EDGES absorption (Barkana et al. 2018; Berlin et al. 2018; Kovetz et al. 2018; Muñoz & Loeb 2018).

The second option, which postulates an excess radio background in the high-redshift universe at the CD, has also been explored (Dowell & Taylor 2018; Feng & Holder 2018). Accretion onto the first black holes has been suggested as a plausible source for the additional background. Such a scenario would also contribute to X-ray and UV backgrounds and thus require the sources to be sited in significantly obscured environments so that these radiations do not disallow the hyperfine populations to be consistent with the EDGES signal (Ewall-Wice et al. 2018); additionally, their energetics are severely constrained by the substantial inverse-Compton losses at high redshifts (Sharma 2018). Mechanisms involving particle interactions, often invoking exotic physics, have also been proposed to produce an excess radiation background that might explain the EDGES absorption (Aristizabal Sierra & Fong 2018; Lawson & Zhitnitsky 2019; Moroi et al. 2018).

Detecting a redshifted 21 cm signal from the CD and the Epoch of Reionization (EoR) is an extremely challenging experiment in long wavelength astronomy, requiring unprecedented measurement accuracy of the radiometer and control of systematics. Owing to the difficulty of such measurements, there have been concerns over the modeling of the EDGES data (Hills et al. 2018, hereafter HKM18) suggesting that resonance modes in the ground beneath the antenna might lead to spurious absorption signatures (Bradley et al. 2019), and concerns regarding contamination from spinning dust (Draine & Miralda-Escudé 2018). Motivated by the unexpected nature of the EDGES claim and the difficulty with finding astrophysical explanations without introducing either new physics or new source populations as extensions to standard models for cosmology and particle physics, we have explored an alternate approach in the analysis of the EDGES spectrum.

In this paper, we use publicly released EDGES data⁴ and employ a different foreground model compared to that adopted by BRM18. Our approach is based on our earlier study (Sathyanarayana Rao et al. 2017b) on modeling the foregrounds for such global CD/EoR measurements. In Section 2, we briefly describe this maximally smooth function approach and present results from its application to the data. In Section 3, we carry out joint modeling of the data using maximally smooth functions to represent the foreground together with flattened absorption with the form of the signal that is reported by EDGES. We compare it with the possibility that the calibrated data have residual systematics of sinusoidal form. Adopting the hypothesis that there is such a residual in the data, in Section 4 we perform joint fitting of a maximally smooth function to represent the foreground, plus a sinusoidal systematic and examine it for an additional 21 cm signal. Section 5 presents a discussion on the physical parameter space populated by favored signals. We present a summary and conclusions in Section 6.

We first estimate the instrument/measurement noise in the EDGES spectrum. The spectrum in the log₁₀(temperature) versus log₁₀(frequency) domain is fitted with polynomials of increasing order, the least-squares fit is subtracted from the data, and the standard deviation in the residual is computed. Additionally, for each residual spectrum, differences between every spectral value and the value in the adjacent spectral bin is computed, and the standard deviation in these differences is computed. For a zero-mean Gaussian random signal, this standard deviation in the differences is $\sqrt{2}$ times the standard deviation in the signal, so we reduce these computed standard deviations by $\sqrt{2}$. Figure 1 shows the run of the standard deviation in the residual on subtracting out best-fit polynomials of progressively greater order. In computing these standard deviations, it may be noted that we have taken into account the progressive reduction in the number of degrees of freedom in the residual data when fitting with polynomials of progressively increasing order. As shown by the figure, the standard deviation of residuals decreases with increasing order, indicating that the EDGES spectrum does have complex structure that requires a very high-order polynomial to model to the sensitivity of the observations. The standard deviations estimated from the differences between adjacent channels vary little with increasing order of the fit polynomial, as expected. The analysis suggests that the instrument measurement noise in the spectrum corresponds to rms spectral fluctuations of about 17 mK.

Zoom In Zoom Out Reset image size

Radiometers, like EDGES, have wide beams that sample a large sky solid angle and hence their spectra represent the average over a large number of individual sources—compact and diffuse—which may have different emission mechanisms and spectral indices and shapes. An average spectrum may no longer be accurately represented with a physical model that is appropriate for a single source. For example, if the spectrum is the sum of the spectra of two sources, and the pair has straight synchrotron spectra but with different spectral indices, the average will not fit a model of a synchrotron source with a single power-law spectrum. In this case, the spectrum at low frequencies would be dominated by the source with the steeper spectrum, and at high frequencies, the source with a flat spectrum would dominate the total. In the log temperature versus log frequency domain, in which the individual sources have linear spectra, the total spectrum would have a curvature that would require a polynomial with increasingly higher degrees to fit with increasing accuracy and lower residuals. For this reason, an attempt to use a physical model with multiple parameters to fit the total spectrum as measured by a radiometer may be limited in accuracy.

There have been a number of parametric and nonparametric approaches to foreground modeling (Harker et al. 2009; Bernardi et al. 2015; Chapman et al. 2015), most of which have no physical basis. BRM18 have adopted a five-term parametric model to represent the foreground in the EDGES spectrum. The five-term model is motivated by the physics of radiation mechanisms relevant to the frequency range of the spectrum, including synchrotron radiation with a curved spectral index, ionosphere emission, and absorption. HKM18 point out that the best-fit parameters have unphysical values; as discussed in BRM18 and Bowman et al. (2018b), there are residual instrumental effects and calibration errors that are expected to be subsumed by their five-term foreground model, and therefore, the fit parameters may not be physical even if the foreground were completely described by their physical model.

Sathyanarayana Rao et al. (2017b) created mock spectra that might be representative of sky spectra observed by radiometers with wide beams. Using all-sky maps that are available at discrete frequencies, the spectra toward individual sky pixels were modeled by fitting a physical model (GMOSS; Sathyanarayana Rao et al. 2017a) to the discrete measurements at that pixel. The physical model for every pixel consisted of a synchrotron power law with a spectral break or the sum of two power laws, depending on whether the pixel spectrum was convex or concave, respectively. To this was added an optically thin thermal emission component, and the pixel spectrum was also allowed to have a low-frequency thermal absorption. Mock spectra were generated, assuming observations with wide beams, by a beam-weighted averaging over the physical sky spectra of the individual pixels. The work showed that the radiometer observations could be fit to mK accuracy using what was called "maximally smooth" functions, and that the residuals to such fits preserved the turning points inherent in any embedded 21 cm signals. Because fitting with maximally smooth functions reveals in the residuals the turning points in any embedded signal, it can also be used as a diagnostic to inspect for nonsmooth frequency structures in the data that may result from instrument-generated additive signals or spectral distortions introduced in the intrinsic foreground spectrum by nonideal, uncorrected properties of the antenna.

As described in Sathyanarayana Rao et al. (2015, 2017b), maximally smooth functions are a special class of polynomials, where the second- and higher-order derivatives of the polynomial are not allowed to have zero crossings across the band of interest. As a result, the constrained polynomial fits only the smooth component of the spectrum and does not subsume any embedded wiggle or higher-order variation present in the data. If, for example, the foreground were a summation of synchrotron spectra with a distribution in spectral indices, and had a 21 cm signal with multiple turning points as an additive, then fitting with a maximally smooth function would fit out the foreground and the residual would reveal the turning points. As the order of the constrained polynomial-form maximally smooth function is increased, the smooth component of the spectrum is fitted with greater accuracy without also fitting out the 21 cm signal. In contrast, fitting with unconstrained polynomials of increasing order would fit the foreground with increasing accuracy; however, the 21 cm signal would also be fitted out with greater accuracy, and the residual would not obviously reveal the turning points.

Shown in the top panel of Figure 2 is the best-fit maximally smooth function. The fit in Figure 2, as well as the fits discussed below for different models, are least-squares fits assuming uniform errors across the spectrum. Ideally, the fits need to account for the variation in errors across the spectrum that arise from varying sky brightness and also the systematic errors, which are difficult to estimate and may dominate the thermal noise from the system temperature. For comparison and clarity in displaying the characteristics of this fit, the best-fit power-law spectrum with a constant spectral index is shown for reference in red, and the best-fit maximally smooth function is shown in blue with the deviations from the best-fit power law enhanced by factor of 30. Clearly, the EDGES spectrum is convex and steepens toward higher frequencies. The run of the spectral index is shown in the middle panel, where it is seen that the spectral index changes from about −2.545 at 50 MHz to about −2.585 at 100 MHz. It is evident from the middle panel of Figure 2 that the variation of the spectral index with frequency has an inflection point. The prominent inflection in the spectrum in the log₁₀(frequency)–log₁₀(brightness temperature) domain is avoided by representing the spectrum in the log₁₀(frequency)–brightness temperature domain; therefore, in this work we fit the EDGES spectrum with maximally smooth functions in the latter domain. It may also be noted here that throughout the analysis presented in this work, we use the EDGES spectrum over their full band of 50–100 MHz, which is the band adopted for analysis in BRM18.

Zoom In Zoom Out Reset image size

We find that we require a maximally smooth function represented by a sixth-order constrained polynomial, with seven polynomial terms, to model the smooth component in this spectrum to within the rms noise in the data. Unlike unconstrained polynomials, the effective number of degrees of freedom in the case of maximally smooth functions would almost always be less than the number of polynomial terms, provided that a sufficiently high-order polynomial form is used. This is because of the constraints on higher-order derivatives, and hence on the coefficients, imposed by the condition that the polynomial be maximally smooth. We thus find that increasing the number of terms of the maximally smooth polynomial beyond seven terms does not significantly improve the fit or change the best-fit parameters.

The residual (data minus fit) is shown in the bottom panel of Figure 2 and has a standard deviation of 44 mK. The standard deviation in the residual is significantly greater than the measurement noise in the data, which is about 17 mK, as discussed above. The large residual clearly shows that the sky spectrum is inconsistent with expectations based exclusively on large-angle averages over physically motivated foreground models. The residual to the maximally smooth fit shows that the EDGES data almost certainly does have an embedded component with multiple turning points, apart from the expected dominant foreground.

We have demonstrated and discussed above that there is significant spectral structure in the EDGES data apart from a maximally smooth foreground. We know of no algorithm that can invert the measured spectrum to decompose it into a maximally smooth foreground plus a unique additive signal. The turning points suggested by the fit presented in the previous section provide a clue, which can serve as a guess of the parametric forms for the additive embedded signal. The data may then be jointly fitted with such parametric forms along with a maximally smooth polynomial that represents the foreground. Any additive embedded signal, which when fitted to the spectrum jointly with a maximally smooth polynomial yields a residual that is consistent with thermal noise, is a plausible signal. Below, we examine the residuals from fits that adopt signal profiles corresponding to the flattened absorption suggested by BRM18 and compare with residuals from adopting a sinusoidal function for the embedded signal.

We begin this section by first assuming that the data has a flattened absorption profile with the form

$$\begin{eqnarray}&&{T}_{21}(\nu )=-A\left(\displaystyle \frac{1-{e}^{-\tau {e}^{B}}}{1-{e}^{-\tau }}\right),\end{eqnarray} \tag{ 1 }$$

with

$$\begin{eqnarray}&&B=\displaystyle \frac{4{(\nu -{\nu }_{0})}^{2}}{{w}^{2}}\mathrm{log}\left[-\displaystyle \frac{1}{\tau }\mathrm{log}\left(\displaystyle \frac{1+{e}^{-\tau }}{2}\right)\right],\end{eqnarray} \tag{ 2 }$$

and parameters suggested by BRM18: A = 0.50, ν₀ = 78.0, w = 19.0, and τ = 7.0. Having first subtracted this profile from the EDGES spectrum, we then fit the remainder with a maximally smooth function. The final residual after subtracting both the flattened absorption with the BRM18 parameters plus this best-fit maximally smooth function is shown in Figure 3. The residual has standard deviation 23.1 mK.

Zoom In Zoom Out Reset image size

Smoothing the residuals enhances the signal-to-noise ratio for broader-scale frequency structures that are coherent across frequency channels (Press et al. 2007). Therefore, the residual, Hanning-smoothed (Oppenheim 1999) to successively lower frequency resolutions, is also shown in Figure 3. For comparison and as a reference, we show in Figure 4 sample residuals containing thermal noise realizations drawn from a Gaussian distribution with zero mean and 17 mK standard deviation, along with smoothed versions of this synthetic spectrum. The standard deviation in this case reduces to 4.8 mK at the maximum smoothing scale of 3.125 MHz. Comparing the residuals in Figure 3 with the expectation in Figure 1, there is clearly unmodeled structure present in the residual, and the standard deviation does not attain that expected from measurement noise. Additionally, comparing the smoothed versions of the residual in Figure 3 with the reference set in Figure 4, it is clear that smoothing the residual does not decrease the standard deviation of the smoothed spectrum as expected given the width of the smoothing functions: Hanning-smoothed with kernel of width 3.125 MHz, the standard deviation in the residual reduces to only 10.5 mK.

Zoom In Zoom Out Reset image size

We next perform a joint fit of a maximally smooth function and the flattened absorption profile, this time allowing the parameters of the absorption to be optimized. The flattened profile is described by its amplitude A, center frequency ν₀, full width at half-power w, and flattening factor τ, and we keep all four parameters unconstrained in the modeling. The outcome of the fit is shown in Figure 5, and the residuals, smoothed to successively lower spectral resolutions, are shown in Figure 6. The best-fit amplitude is 921 ± 35 mK centered at frequency 78.2 ± 0.1 MHz. The best-fit FWHM is 19.2 ± 0.2 MHz, and the flattening factor is 4 ± 0.28. The resulting standard deviation in the residual structures is 22.5 mK, which reduces to 9.0 mK when Hanning-smoothed with a kernel of 3.125 MHz. It may be noted here that the best-fit depth of the absorption we have obtained differs significantly from that reported by BRM18, lying at the boundary of their 99% confidence interval. If the foreground were modeled as a maximally smooth function instead of the five-term physical model proposed by BRM18, the absorption is inferred to have an amplitude of about 0.9 K instead of the 0.5 K inferred by BRM18. As we have argued above, the foreground component in wide-beam radiometer measurements is not expected to follow physically motivated forms appropriate for individual sources, but is expected to be fit with maximally smooth functions. Therefore, if the 21 cm signal from the CD has the form of a flattened absorption profile, the depth of the absorption is 0.9 K, and the tension with previous theoretical expectations is all the more worsened.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Motivated by the above finding, and also the periodic sinusoidal nature of the residual structure in Figure 2, we have attempted to model the spectrum jointly with a maximally smooth function and a sinusoid with amplitude, frequency, and phase as free parameters. In Figure 7, we show the best-fit results, and in Figure 8, we show the residuals Hanning-smoothed to lower resolutions. The best-fit amplitude of the sine function is 60 ± 10 mK, with a period of 12.3 ± 0.1 MHz. This is similar to the result of HKM18, where the foreground model was the same as that adopted by BRM18. The standard deviation in the residuals after subtracting the best-fit maximally smooth foreground and best-fit sinusoid is 22.9 mK, employing three free parameters to model the nonsmooth component. The resulting standard deviation is only marginally higher than that obtained assuming a flattened absorption profile for the nonsmooth component, which uses four free parameters. Smoothed with a Hanning kernel of width 3.125 MHz, the residuals in the case assuming a sinusoidal form for the signal have a standard deviation of 9.4 mK, which is again marginally higher than that for the case assuming a flattened Gaussian-form signal.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Although the residuals in the case assuming a sinusoidal form for the signal have a slightly higher standard deviation, sinusoidal signals have a smaller number of free parameters compared to a flattened-absorption-type signal. Therefore, to compare the goodness of fits and usefulness of the modeling in the two cases, we compute the Bayesian information criterion (BIC) for the two cases.

Assuming a normal distribution of the model errors, the BIC can be defined as

$$\begin{eqnarray}&&\mathrm{BIC}=n\,\mathrm{ln}\left(\displaystyle \frac{1}{n}\displaystyle \sum _{i=1}^{n}{({x}_{i}-\hat{{x}_{i}})}^{2}\right)+k\,\mathrm{ln}(n),\end{eqnarray} \tag{ 3 }$$

where n is the number of data samples, k is the number of parameters, and ${\sum }_{i=1}^{n}{({x}_{i}-\hat{{x}_{i}})}^{2}$ is the reduced sum of squares (Schwarz 1978; Priestley 1981). The model with the lower BIC is preferred. We take the number of free parameters for the maximally smooth function to be the number of terms in its polynomial expansion, though the effective number of free parameters would be fewer as discussed above in Section 2. We thus present the worst-case scenario for the BIC by making this assumption. Because we use the full band data set, the number of data samples remains the same, n = 123, in all cases. We adopt the scale in Kass & Raftery (1995) for interpreting BIC values, in which it is argued that if the difference in BIC values for two models, ΔBIC, exceeds 6, then there is strong evidence against the model with the higher BIC. Unsurprisingly, both the sinusoidal and flattened absorption profiles result in very similar BIC values of −876.3 and −876.1, respectively (ΔBIC ∼ 0.2), which confirm the finding that neither of these models ought to be preferred over the other, at least in so far as this formal statistical comparison is concerned (Kass & Raftery 1995). It is therefore equally likely that the data may well have an unwanted sinusoidal systematic of hitherto unknown origin embedded in them, instead of an unexpectedly deep flattened 21 cm absorption.

Sinusoidal systematics may arise due to impedance mismatch between the antenna and subsequent electronics in the receiver, or along the analog signal path. As a result, all receiver noise in the signal path within the receiver electronics will suffer multipath propagation owing to reflection at the antenna and also at impedance mismatches in the signal path. Multipath propagation within the signal path leads to sinusoidal spectral structure in the measured spectrum (Meys 1978). If the foreground sky signal suffers from multipath propagation, the resulting sinusoid amplitude would directly scale with the sky brightness temperature. However, when the receiver noise suffers from multipath propagation, the amplitude of the sinusoid would be constant in time to the extent that the receiver noise is constant. Because Bowman et al. (2018b) do not observe a change in the amplitude of the flattened Gaussian model with a change in sky brightness temperature, we infer that the sinusoid systematic might be an additive component arising due to the multipath propagation of receiver noise.

It may be noted here that the best-fit period of the sinusoid is close to the spectral periodicity expected in the antenna gains, as shown in Figures 4(b) and (c) in BRM18. The period for the sinusoid that we infer to be likely present as a spurious component in the EDGES low-band data is close to half of the period of the spectral ripple that appears in the residuals of the published EDGES high-band results (see Figure 4(b) in Monsalve et al. 2017). The amplitude of the best-fit sinusoid in the low-band case is about 60 mK; this is very similar to the amplitude of the sinusoidal feature in the residual spectrum of Monsalve et al. (2017). The low-band antenna is a factor of 2 scaled-up version of the high-band antenna. Such a scaling would have the same impedance properties in corresponding bands, and hence the amplitude of the sinusoid is expected to be the same in the two bands. Given that this is consistent with what we observe, it is possible that the sinusoidal ripple is an inherent characteristic related to the antenna structure or its interaction with the environment.

In summary, (i) on fitting with just a maximally smooth function (Figure 2), the residual was highly suggestive of an additional sinusoidal component, (ii) multipath propagation and reflections at impedance mismatches within the signal path both lead to a spurious sinusoidal component, (iii) interpreting the EDGES data as implying an expectedly deep absorption causes substantial tension with prevalent theoretical models based on previous observations, and (iv) the BIC computed in the two cases, taking into account the standard deviations in the residuals in the two cases along with the numbers of free parameters in describing the nonsmooth flattened absorption profile and sinusoidal profile, suggest that neither of these models is preferred over the other. For all these reasons taken together, we model the data as containing a foreground of a maximally smooth functional form plus a spurious sinusoidal component, and test in Section 4 below for the presence of any additional cosmological 21 cm signals.

As discussed above, we adopt a combination of a maximally smooth function along with a three-parameter sinusoid to represent, respectively, the foregrounds and residual systematic structure in the EDGES spectrum. With this model, we inspect the data for the presence of any additional physically motivated 21 cm signals. Because most accepted theoretical models predict an absorption feature in the CD as a result of, first, a Wouthuysen–Field coupling of the spin temperature to the low kinetic temperature, followed by X-ray heating, we examine the data for an additional Gaussian-form absorption signal. We parameterize the absorption signal to be of Gaussian form with the amplitude, frequency, and width as free parameters (Bernardi et al. 2016). A joint fit of a maximally smooth function, plus a sinusoid, plus a Gaussian-form absorption yields a best-fit amplitude of 133 ± 60 mK for the Gaussian, with width at half-power of 9 ± 3 MHz and centered at 72.5 ± 0.8 MHz. This best-fit model is shown in Figure 9, with smoothed versions of the residuals shown in Figure 10. We find that including the Gaussian appreciably reduces further the standard deviation of the residuals to 20.1 mK; smoothing the residuals with a Hanning kernel of width 3.125 MHz yields a residual with standard deviation 4.6 mK. In comparison, the joint fit of the five-term physical model describing the foreground together with a four-term flattened absorption profile, adopted by BRM18, results in a residual with a standard deviation of 24.6 mK. Though the number of free parameters in the model adopted here is larger, information metrics like the BIC that take into consideration the reduction in standard deviation along with the number of free parameters strongly favor the model suggested here, which includes a sine-form spurious additive, over the model adopted in BRM18 (ΔBIC ∼ 26). Comparing different models for the spectrum based on the BIC, the model adopted in BRM18 is least preferred, followed by a maximally smooth function plus a flattened Gaussian profile (ΔBIC ∼ 8). The most preferred model is a combination of a maximally smooth function, sinusoid, and Gaussian: when compared to maximally smooth function plus a flattened Gaussian profile, this model is preferred by the Bayesian information metric with a ΔBIC ∼ 18. Table 1 summarizes the models explored, along with the standard deviation of the residuals and the BIC.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Cohen et al. (2017) provided an atlas of 265 theoretically plausible global 21 cm signal templates that were generated by seminumerical simulations consistent with previous observations that constrain the physics and sources through the CD and reionization. The atlas sparsely covers the signal space and spans the range of absorption profiles expected in the CD. We examined whether these templates might be allowed by the EDGES data, apart from a maximally smooth foreground and spurious sinusoidal component. The EDGES spectrum is examined for each of the templates in the atlas, one by one, by performing a joint fit of the EDGES data to a maximally smooth function, sinusoid, plus the 21 cm signal template S₀(ν) multiplied by a scale factor c, which is a free parameter in the joint fit. The order of the maximally smooth function is chosen to be large enough so that the resulting residuals do not change with further increase in the number of terms of the polynomial. This approach is similar to the method adopted in Monsalve et al. (2017) and Singh et al. (2018). We select those 21 cm signal templates of the atlas for which the scale factor c is consistent with unity but inconsistent with zero. The selection criterion is

$$\begin{eqnarray}&&0\leqslant (c-\delta c)\leqslant 1\leqslant (c+\delta c),\end{eqnarray} \tag{ 4 }$$

where c is the best-fit scale factor and δc is the associated error derived from the uncertainty in the least-squares fit. We find that 7 of the 265 signals in the atlas satisfy the above criterion, with varying levels of significance.

We use the standard score, $c/\delta c$, to represent the significance. It measures the power of rejection of the null hypothesis, i.e., the power with which the hypothesis that the scale factor is consistent with zero is rejected. The standard score has units of standard deviation. We accept only those signals as being favored by the data in which the standard score is more than 2, which implies that these signals are consistent with the data and preferred by the data with greater than a 2 standard-deviation confidence compared to the null hypothesis. With this additional criterion, we find that three out of the aforementioned seven signals are accepted, and they have standard scores of 2.9, 3.3, and 3.4 respectively. The standard deviation of the residuals remains at ∼22.8 mK in all three cases.

In Figure 11, we show the three 21 cm templates that are favored by the EDGES data, out of the 265 atlas templates considered. All three yield best-fit scale factors close to unity in the joint fit of a maximally smooth function, sinusoid, and the template multiplied by the scale factor, which is a free parameter in the fit. Figure 12 shows the best-fit scale factors for the three signals selected by the data along with the extents of their errors. We list the astrophysical parameters for these signals in Table 2 and discuss the allowed parameter space in Section 5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The central frequency, the depth, and the width of the absorption trough are regulated by high-redshift stars and black hole formation. The low-frequency end of the absorption trough is related to the onset of star formation and the build-up of the Lyα background radiation, while the birth of the first X-ray sources triggers cosmic heating and marks the point of the deepest absorption. The relative timing of these events, as well as the rates of star formation and heating, determines the width and depth of the absorption feature. Therefore, astrophysical properties such as star formation efficiency (f_*), minimum mass of the first star-forming halos (M_min or, equivalently, the circular velocity of these halos, V_c), and properties of the first X-ray sources such as the X-ray heating efficiency (f_X) and spectral energy distribution (SED) can be constrained with the high-redshift 21 cm signal. Because the process of reionization is inefficient at high redshifts in the absence of massive galaxies, parameters such as the total CMB optical depth (τ) and the mean free path of the ionizing photons (R_mfp) do not play a significant role at high redshifts. For each high-redshift star formation scenario, the reionization parameters (τ and R_mfp) can be tuned to render the reionization history in agreement with available observations (e.g., Greig & Mesinger 2017; Planck Collaboration et al. 2018; Greig et al. 2019).

We examined, above, 265 standard astrophysical models presented by Cohen et al. (2017) and found that three of them are consistent with the EDGES low-band data. To produce this set of models, ΛCDM cosmology was assumed, and the following astrophysical parameters were varied in the broadest range allowed by observations: f_*, V_c, f_X, τ, and R_mfp; in addition, hard versus soft X-ray SEDs were explored. Earlier, 25 models out of this set, sharing low X-ray efficiency (f_X = 0 or f_X = 0.1) and rapid reionization, were found inconsistent with SARAS 2 data in the 110–200 MHz frequency range (Singh et al. 2017, 2018). Extending the data set of astrophysical models and varying, in addition to the above-mentioned parameters, the slope of the X-ray spectrum and the minimum frequency of X-ray photons, a larger set of models was produced and examined with the EDGES high-band data (A. Cohen et al. 2019, in preparation; Monsalve et al. 2019). Marginalizing over the rest of the parameters, Monsalve et al. (2019) found that the high-band data alone disfavor low X-ray heating efficiency (f_X < 0.0042) and high minimum halo mass (V_c > 19.3 km s⁻¹, equivalent to ∼1.3 × 10⁸ M_⊙ at z = 10).

The location of the absorption feature within the EDGES low-band frequency range implies early onset of star formation and efficient X-ray heating, which is consistent with the nondetection of the 21 cm signal by SARAS 2 and the EDGES high band. Interpreting the EDGES detection, Mirocha & Furlanetto (2019) indicated that in order to ensure the absorption feature is within the band, high star formation efficiency of a few percent in faint galaxies below the current detection threshold (halo masses of 10⁸–10¹⁰ M_⊙) is required in addition to a luminous X-ray population, while Schauer et al. (2019) found that inefficient star formation in mini halos is required. In agreement with these studies, we find that the three astrophysical models selected by our analysis above feature star formation efficiency of f_* = 5% (case A in Table 2) and f_* = 15.6% (cases B and C), minimum mass of star-forming halos of ∼4 × 10⁸ M_⊙ at z = 17 (corresponding to V_c = 35.5 km s⁻¹), and moderately high X-ray heating rate of f_X = 14.7 (case A) and f_X = 1.58 (cases B and C). Out of these three cases, the τ of A and B are consistent with the latest Planck data (within 2σ of the best-fit τ = 0.054 ± 0.007; Planck Collaboration et al. 2018), while the τ of C is within 5σ.

It may be noted here that the above exploration of the astrophysical parameter space is far from being exhaustive and a more rigorous analysis is required to establish the ranges of astrophysical parameters favored by the EDGES low-band data. In particular, recognizing that the standard deviation of the residuals in the case of the three models selected is higher than that for a Gaussian-form absorption, and recognizing that the location of the peak absorption and also the FWHM in the case of the three selected 21 cm profiles are somewhat different from those for the Gaussian-form signal, further exploration of the physical parameter space is needed to generate models that might provide residuals with standard deviation closer to that for the case of a Gaussian absorption signal or, better still, closer to the expected measurement noise in the EDGES low-band data.

The EDGES claim of a detection of a redshifted 21 cm signal from the CD as a flattened absorption profile centered at 78 MHz, with amplitude of 0.5 K, has evoked a wide range of nonstandard theoretical explanations for the signal's origin. However, it has also triggered discussions on the modeling of the data, and questions whether the understanding of known systematics was adequate to account for, model, and subtract their contributions. Systematic errors might be caused by an inadequate foreground model, uncalibrated/unmodeled responses of the instrument to the sky, ground, or receiver noise, or a combination of these. Accurate calibration, understanding, and thereby modeling of systematics, and analysis that separates systematics from any cosmological 21 cm signal, are challenges in making an unambiguous detection.

We use the publicly available calibrated and time-averaged spectrum from EDGES low band in the frequency range of 50–100 MHz. The wide beams of radiometers such as EDGES average sky spectra over large angles and hence over a large number of sources with a distribution of radiation mechanisms and parameters; therefore, the average spectra cannot necessarily be fit with physical models with physically acceptable parameters. Previous work by Sathyanarayana Rao et al. (2017b) suggests that such average spectra may be modeled using maximally smooth functions to the required accuracy. Therefore, we adopt the maximally smooth function as the representation for the foreground, in which the inherent constraints on higher-order derivatives, and thereby the coefficients, make it impossible to overfit, or fit out embedded high-order spectral structures. Adopting such a foreground model, we are led to the conclusion that a flattened absorption signal in the data would have a best-fit amplitude of 0.9 K, substantially deeper than the 0.5 K depth inferred by BRM18.

Because the absorption depth inferred from our analysis, and also the depth inferred by BRM18, is unexpected and require nonstandard explanations, we have explored whether there may be uncalibrated systematics in the data that are the cause of the unexpected result. The residuals to the fit with a maximally smooth foreground model indicate that there may be a sinusoidal form systematic in the data; this has also been suggested previously by HKM18. Sinusoidal residuals are indeed commonly seen in spectral radiometers and often arise from multipath propagation of wide-band noise within the signal path. We find that the EDGES data are equally well fit with a foreground model plus sinusoidal systematic, compared to a foreground model plus flattened absorption.

Assuming that such a sinusoidal systematic does exist in the data, along with the foreground, we find that the data then favor the presence of an additional Gaussian absorption feature with amplitude 133 mK: a depth that is well within the expectations from standard astrophysical cosmology. Further, such a model is favored by information metrics over the model adopted by BRM18, due to the 34% lower variance in the residuals. We then examined whether the data might allow for any of the theoretically computed 21 cm profiles, using the atlas of templates provided by Cohen et al. (2017) that are based on standard cosmology, constrained by previous observations, and represent a sparse sampling of the signal space. We performed joint fits to the data using a maximally smooth function for the foreground and allowing for a sinusoidal systematic, and find that the data indeed favor a set of three templates with close to three-standard-deviation significance standard scores. Standard astrophysical models with the background ΛCDM cosmology, therefore, are found to fit the data well. The central frequency, the width, and the depth of the best-fit absorption profile imply early onset of star formation in small dark matter halos of M_h ≲ 10⁹ M_⊙ at few percent efficiency. In addition, moderately high X-ray heating is required with f_X ≳ 1, which is consistent with standard assumptions for a high-redshift population of X-ray binaries.

We conclude that the EDGES data are indeed likely to have detected a 21 cm signal from the CD, although with an absorption depth that is consistent with standard models.

We thank Anastasia Fialkov for useful discussions and sharing the set of 21 cm signal templates used in exploring the consistency of the data with plausible models. We also thank Richard Hills for his comments and suggestions on an earlier version of this manuscript. We thank the EDGES collaboration for making the data publicly available. S.S. acknowledges support from McGill Astrophysics postdoctoral fellowship.