Formulating and critically examining the assumptions of global 21-cm signal analyses: How to avoid the false troughs that can appear in single spectrum fits

The EDGES team has worked extensively on the calibration of their receiver (Monsalve et al., 2017a, b). Assuming solar, weather, RFI, and other transient events are removed, due to the rigor of their lab measurements and calibration strategy, it is reasonable to assume that, to the mK level, their data consists solely of beam-weighted foreground and global 21-cm signal.

Although the covariance distribution of the noise is not known, the residuals presented in Bowman et al. (2018) seem to have a relatively flat magnitude of 20 mK when many smooth modes are removed. So, for our simulations we use ${\boldsymbol{C}}=\sigma^{2}{\boldsymbol{I}}$ where $\sigma=20$ mK and ${\boldsymbol{I}}$ is the identity matrix.

The EDGES analysis involves multiple linear foreground models, but a common model in the literature is a polynomial in $\ln{(\nu/\nu_{0})}$ multiplied by $(\nu/\nu_{0})^{-2.5}$ where $\nu$ is the observed frequency and $\nu_{0}$ is a reference frequency,

This is equivalent to assuming that ${\boldsymbol{F}}$ is a matrix with the terms $(\nu/\nu_{0})^{-2.5}[\ln{(\nu/\nu_{0})}]^{k-1}$ as its columns and the foreground coefficient vector ${\boldsymbol{a}}$ is given by ${\boldsymbol{a}}^{T}=\begin{bmatrix}a_{1}&a_{2}&\cdots&a_{N_{t}}\end{bmatrix}$. This model is based on the Taylor series approximation of $T_{0}(\nu/\nu_{0})^{\beta}$ around $\beta=-2.5$, which should adequately describe the intrinsic synchrotron foreground in each sky direction (see also Hibbard et al., in preparation).

In Bowman et al. (2018), the EDGES team uses a flattened Gaussian signal model of the form

where

In this section, we adopt the same model.

The antenna beam used in our simulations is a Gaussian beam with a frequency dependent Full Width at Half Maximum (FWHM), $\alpha(\nu)$, i.e.

The FWHMs we use are parabolas with fixed endpoints at $(50\text{ MHz},83^{\circ})$ and $(100\text{ MHz},104^{\circ})$⁴⁴4These endpoints are similar to those of the beam of the blade antenna used by EDGES (see extended data figure 4a of Bowman et al., 2018). and varying curvature, $c$. This means, they have the form

where $\alpha_{l}(\nu)$ is the line connecting the two endpoints and $\alpha_{c}(\nu)$ is the parabola with unit curvature and zeros at the endpoint frequencies. $\alpha_{l}(\nu)$ and $\alpha_{c}(\nu)$ are given by

For illustration purposes, we use $c/(1^{\circ}/\text{MHz}^{2})$ values of $8.0\times 10^{-3}$, $1.6\times 10^{-2}$, $2.4\times 10^{-2}$, and $3.4\times 10^{-2}$. The FWHM curves of the four beams used in our simulations are shown in the upper left panel of Figure 2. To ensure that results are not tainted by any fraction of the beam below the horizon, the beams are normalized such that

where $\mathcal{H}$ represents the half of the sphere above the horizon, i.e. $0\leq\theta\leq\pi/2$ and $0\leq\phi<2\pi$.

For the simulations used in this paper, we use the 408 MHz map from Haslam et al. (1982), scaled to lower frequencies via multiplication by a power law with spectral index -2.5, i.e.

The time dependence indicated in both $T$ and $T_{\text{Haslam}}$ is determined by the rotation of the Earth from the EDGES latitude as a function of sidereal time. The Haslam maps, rotated to match the zenith direction at 0, 3, and 6 hr LST (local sidereal time) at $26.7^{\circ}$ S latitude, are shown in the top left, top right, and bottom left panels of Figure 1.

The simulations were performed by averaging together equally all beam-weighted foregrounds between LST hours 0-6 from the EDGES latitude. As the central bulge of the Milky Way is closest to the zenith between LST hours 17 and 18, this strategy amounts to averaging during low foreground times. This observation strategy leads to an effective foreground given by

and shown at 80 MHz in the bottom right panel of Figure 1.

No 21-cm global signal was included in the simulations. All fits shown in Figure 2 are performed on the beam-weighted foreground alone with noise.

To each simulated measurement of the beam-weighted foreground we add the same realization of noise at the 20 mK level (see Section 3.1.2) to control for the effect of noise without ignoring it.

Using identical and independent basis sets for each spectrum leads to a degeneracy between the foreground and signal models, causing infinite uncertainties. To see this, suppose that each spectrum has its own polynomial basis, represented by the matrix ${\boldsymbol{F}}_{0}$. In this case, the full foreground basis matrix is

Multiplying this basis by a coefficient vector of the form ${\boldsymbol{x}}^{({\text{fg}})}=\begin{bmatrix}{\boldsymbol{\zeta}}^{T}&{\boldsymbol{\zeta}}^{T}&\cdots&{\boldsymbol{\zeta}}^{T}\end{bmatrix}^{T}$, leads to a foreground vector of the form ${\boldsymbol{F}}{\boldsymbol{x}}^{({\text{fg}})}={\boldsymbol{\Psi}}{\boldsymbol{F}}_{0}{\boldsymbol{\zeta}}$. Since, in this case, the foreground basis can generate vectors in the column space of ${\boldsymbol{\Psi}}$, the foreground and signal bases are degenerate and the uncertainties diverge to infinity. This is true even if the ${\boldsymbol{F}}_{0}$ basis can exactly fit the foregrounds in each spectrum. Note that this explanation holds no matter what the true beam-weighted foreground is because it only depends on the form of the beam-weighted foreground model.

We wish to find a basis that satisfies ${\boldsymbol{F}}^{T}{\boldsymbol{C}}^{-1}{\boldsymbol{F}}={\boldsymbol{I}}$ and can accurately fit the beam-weighted foreground expected for each spectrum. We suggest generating a simulated training set of foregrounds of the form

where ${\boldsymbol{b}}^{(m)}_{n}$ is the $n^{\text{th}}$ spectrum of the $m^{\text{th}}$ simulation and there are $N_{t}$ total simulations. For a given number of basis vectors $N_{F}$, weighted Singular Value Decomposition (SVD) of ${\boldsymbol{B}}$ as in Tauscher et al. (2018) yields a basis ${\boldsymbol{F}}$ satisfying the normalization conditions.⁶⁶6Weighted SVD refers to a decomposition of the form ${\boldsymbol{B}}={\boldsymbol{U}}{\boldsymbol{\Sigma}}{\boldsymbol{V}}^{T}$ where ${\boldsymbol{U}}^{T}{\boldsymbol{C}}^{-1}{\boldsymbol{U}}={\boldsymbol{I}}$, ${\boldsymbol{V}}^{T}{\boldsymbol{V}}={\boldsymbol{I}}$, and ${\boldsymbol{\Sigma}}$ is a rectangular diagonal matrix with decreasing non-negative elements on the diagonal. The basis matrix ${\boldsymbol{F}}$ is made from the first $N_{F}$ columns of ${\boldsymbol{U}}$. This basis will encode expected correlations between the different spectra in the foreground model directly, which is a crucial aspect of the analysis necessary to obtain finite errors (see also Tauscher et al., 2020) when allowing for any signal to be fit as discussed in Section 4.1.

The fits in Section 3 (Figure 2) show that false troughs can appear when fitting a single spectrum. Some may argue that troughs cannot be created by the foreground because the spectral dependence of the foreground does not allow it. However, one needs to keep in mind that what is measured is not the intrinsic spectral structure of the foreground radiation, but rather the spectral structure of the beam-weighted foreground. Since the angular structure of the beam weighting changes from frequency to frequency the foreground component of the measured spectrum is composed of different ratios of each direction’s radiation at each frequency.

This means that a linear model that fits the intrinsic foreground spectrum of each pixel, like the model from Equation 6, will not necessarily fit the beam-weighted foreground. To illustrate this, suppose that every sky direction’s intrinsic foreground spectrum can be fit by a model using a basis matrix ${\boldsymbol{F}}_{0}$, i.e. ${\boldsymbol{T}}(\theta,\phi)={\boldsymbol{F}}_{0}{\boldsymbol{x}}(\theta,\phi)$ where the $k^{\text{th}}$ element of ${\boldsymbol{T}}(\theta,\phi)$ is the foreground spectrum at the $k^{\text{th}}$ frequency. If the beam is achromatic, then the beam-weighted foreground is given by

where $B(\theta,\phi)$ is the beam pattern at every frequency that satisfies $\int B(\theta,\phi)\ d\Omega=1$. This means that a linear model that fits ${\boldsymbol{T}}(\theta,\phi)$ for every angle also fits the beam-weighted foreground ${\boldsymbol{y}}^{({\text{fg}})}_{\text{achromatic}}$, i.e. the foreground model residual is ${\boldsymbol{\delta}}_{\text{achromatic}}=({\boldsymbol{I}}-{\boldsymbol{\Phi}}){\boldsymbol{y}}^{({\text{fg}})}_{\text{achromatic}}={\boldsymbol{0}}$. In the case of beam chromaticity, the beam is a diagonal matrix with the diagonal entries being the beam patterns at each frequency, satisfying $\int{\boldsymbol{B}}(\theta,\phi)\ d\Omega={\boldsymbol{I}}$. In this case, the beam-weighted foreground is

and the residual

is not necessarily zero,⁹⁹9The only situation in which ${\boldsymbol{\delta}}^{\text{chromatic}}$ is necessarily zero for all beams, ${\boldsymbol{B}}(\theta,\phi)$, occurs when ${\boldsymbol{I}}={\boldsymbol{\Phi}}$. But, this implies that projecting into the span of the foreground basis has no effect, which can only happen if the foreground basis is complete. In this case, the signal model will be degenerate with the foreground model and therefore cannot be extracted. leading to a possible failure of the assumption that the foreground model is adequate (Assumption 3). Therefore, while troughs are unlikely to appear in intrinsic foreground spectra, they may well appear when beam-weighted foregrounds are fit with a model of a trough and an intrinsic foreground spectrum model simultaneously.

In several of their works (see, e.g., Bowman et al., 2018; Monsalve et al., 2018; Mozdzen et al., 2019), the EDGES team uses what they term the “beam chromaticity factor” (see Equation 14 of Monsalve et al., 2017a), defined at each LST through

where the LST dependence comes from how the temperature map is rotated with respect to the beam and $B_{\text{ref}}(\nu,\theta,\phi)$ and $T_{\text{ref}}(\nu,\theta,\phi)$ are the assumed forms of the beam and foreground. The goal of this factor is to make the spectrum appear as if it was measured with an achromatic beam (specifically, the beam at $\nu=\nu_{\text{ref}}$) so that, as discussed in Section 5.1.1, basis vectors that fit the intrinsic foreground spectra can be used to fit the beam-weighted foreground spectrum. If we denote a spectrum measured at a single LST by $y(\nu)$, then the EDGES beam calibration forms a corrected spectrum, $y^{\prime}(\nu)=y(\nu)/C(\nu)$. Neglecting noise and possible receiver errors for the sake of clarity, we can express the measured beam-weighted foreground spectrum through

where $B(\nu,\theta,\phi)$ and $T(\nu,\theta,\phi)$ are the unknown forms of the true beam and foreground. With these definitions, the corrected spectrum is given by

The intention is for the fraction on the second line of Equation 26 to be 1 so that $y^{\prime}(\nu)$ is simply the reference foreground weighted by the reference beam at the reference frequency. However, in order to assume that the fraction is 1, one must assume that $B_{\text{ref}}(\nu,\theta,\phi)=B(\nu,\theta,\phi)$ and $T_{\text{ref}}(\nu,\theta,\phi)=T(\nu,\theta,\phi)$; but, this is not a reasonable assumption as foreground and beam models are likely imperfect.¹⁰¹⁰10If the reference beam and foreground were equal to the true beam and foreground, then one would know the exact beam-weighted foreground and should simply subtract it from the measured spectrum, leaving only noise and the 21-cm signal. The beam chromaticity factor calibration method is therefore prone to introduce significant biases. On the contrary, our method of deriving modes describing the beam-weighted foreground from a large sample of simulations and observations allows uncertainties in both the beam and foreground to be included in the analysis. However, as will be elaborated on in Section 5.2, it is important to note that our method is contingent on the ability of the simulations and observations to accurately describe the data being analyzed. SVD can only provide eigenmodes in the space spanned by the training set curves.

In Section 3, in addition to the Haslam 408 MHz map, we also tested a toy model of galactic emission described by

This is a model that treats the Galactic plane as a $8^{\circ}$ FWHM Gaussian in colatitude, $\theta$, with a maximum given by 300 K at the Galactic center and 200 K at the anticenter and asymptotes to 25 K away from the plane. This map is then scaled by $[\nu/(408\text{ MHz})]^{-2.5}$. Due to the fact that the false troughs remain after this change (see e.g. the dashed green line in the lower left panel of Figure 2), we infer that the chromaticity of the quadratic beams used in Section 3 interacts with the galactic plane as seen from the EDGES latitude to generate residuals that are better fit with a flattened Gaussian plus foreground model than the foreground model alone, which could lead to analyses falsely concluding that there are troughs in the sky-averaged radio spectrum.

Since the false troughs found in Section 3 are fits to aspects of the beam-weighted foreground, we urge EDGES to perform the same experiment and analysis from the northern hemisphere where the galaxy behaves very differently in the sky than from the southern hemisphere. Because of the different orientation of the galactic plane, the distortions caused by chromatic beams in our simulations are more easily fit by foreground models like that in Equation 6 at northern latitudes than southern latitudes.

The natural endpoint of the EDGES-like analysis is not a maximum likelihood fit like the one shown in Figure 1 of Bowman et al. (2018), on which Figure 2 was based. Instead, one wishes to explore the parameter distribution using a sampling method like Markov Chain Monte Carlo (MCMC) exploration. However, the resulting distribution is only meaningful if the pieces of the model (i.e. foreground and signal models) accurately describe their respective components. In fact, parameter distributions from global 21-cm signal fits to observations are generated by a complex combination of the following four factors: 1) the foreground model parameterization; 2) the bias of the foreground model; 3) the signal model parameterization; 4) the bias of the signal model.

Ideally, factors 2 and 4 are unimportant because the foreground and signal models can fit the true foreground and signal to within the noise level. However, it is impossible to fully verify that this is the case in practice; so, their effects must be considered. Section 3 is not meant to explain the exact form and parameter distribution of the trough observed by EDGES, but instead to show that generic polynomial foreground models like those used by EDGES do not account for beam chromaticity and may lead to untrustable results. Essentially, we are pointing out that factor 2 is likely affecting the EDGES analysis.

In addition to generating bias in fits, residuals from fits with inaccurate foreground models may have led EDGES to choose the unjustified flattened Gaussian signal model, possibly furthering the signal bias of factor 4 (with respect to a more conventional signal model) silently while significantly reducing overall residuals.¹¹¹¹11If true, this would be a case of factor 2 affecting factor 3. For robustness reasons, it is important to justify the signal model in a way external to the observations themselves.

Due to the complex interactions between the four factors in generating a posterior distribution, MCMC exploration of distributions is outside the scope of this paper and is left for an extended examination in future work.

While moving the experiment to a different location could provide evidence that the reported trough is a product of an inadequate foreground model, this does not solve the underlying problem—that the foreground model does not account for beam chromaticity.

The MAA proposed in Section 4 avoids this problem by making a basis specific to the given antenna by performing SVD on a training set of simulations made with that antenna’s beam (see Section 4.2.2).

In addition, the MAA allows any possible global signal, avoiding unwanted bias from generating a signal model that interacts with the foreground model bias to produce false results. In fact, under the assumption that the foreground model generated by SVD is adequate, the MAA signal mean ${\boldsymbol{\gamma}}^{(21)}$ and covariance ${\boldsymbol{\Delta}}^{(21)}$ constitute a direct measurement not of a signal model, but of the signal itself. In fact, since all degrees of freedom are retained, a physical model can be fit directly to the constraints implied by ${\boldsymbol{\gamma}}^{(21)}$ and ${\boldsymbol{\Delta}}^{(21)}$ (such as those shown in the bottom panel of Figure 3).

Note that, for simplicity, the cosmic microwave background (CMB) is not included in the training set used in Section 4, which merely uses the Haslam map scaled with a spectral index of -2.5. If it was included (as it should be when analyzing data from a real experiment),¹²¹²12More precisely, the CMB should be subtracted from the foreground model used to generate the training set so that it can be found by the signal fitting and subtracted after the fact. then the signal mean ${\boldsymbol{\gamma}}^{(21)}$ would include a 2.725 K flat spectrum component which would need to be subtracted out.

Another important note about the MAA is that it is predicated upon the beam-weighted foreground training set provided to it. The results may change if the training set is changed, even if the data being analyzed remain the same. The training set used in Section 4 contained a wide variety of beam FWHMs but only one intrinsic foreground map. A different training set built using many intrinsic foreground maps and one beam might produce errors that differ significantly from the ones presented in Figure 3. However, even as different training sets would lead to different levels of precision, the accuracy of the MAA should remain steady as the training set changes as long as the true beam-weighted foreground is encompassed by the SVD eigenmodes of the training sets.¹³¹³13Here, we use precision to refer to the size of uncertainties, which is known even when the true 21-cm signal is unknown and accuracy to refer to the size of the difference between the signal reconstruction and the true signal with respect to the size of the uncertainties. One aspect of this effect is that some training sets will require more modes to be included in the foreground model than others, forcing one to degrade precision in order to retain accuracy.

Slosar (2017) pointed out that, much like with the Cosmic Microwave Background (CMB), Doppler shifting of the 21-cm background induces a dipole component that is related to the monopole. In particular, the dipole component, $\delta T_{b}^{(\text{dip})}$, of the 21-cm signal is related to the monopole spectrum, $\delta T_{b}^{(\text{mon})}$, through¹⁴¹⁴14Note that this assumes that there is no intrinsic (i.e. non-motion-induced) dipole of the 21-cm signal.

where $\beta$ is the magnitude of our velocity with respect to the CMB as a fraction of the speed of light and $\psi$ is the angle between that velocity and the line of sight. Based on CMB observations, $\beta\approx 1.2\times 10^{-3}$. Deshpande (2018) suggested that Equation 28 should be considered an essential qualifying test of any measurement of the global signal. Due to the fact that (by Equation 28), $\delta T_{b}^{(\text{dip})}$ is linear in $\delta T_{b}^{(\text{mon})}$, it could conceivably be included in the MAA in order for any fit to pass this test by default, essentially extending the minimal signal assumption introduced in Section 4.1. To do so, we note that the sum of the monopole and dipole components is

where $s(\nu)=\delta T_{b}^{(\text{mon})}(\nu)$. To determine how this would appear to a single antenna experiment, we must multiply by the beam and integrate over all angles. Defining $t_{k}(\nu)\equiv\int B_{k}(\nu,\theta,\phi)\ \delta T_{b}^{(\text{mon}+\text{dip})}(\nu,\psi)\ d\Omega$ as the signature of the 21-cm signal in the spectrum measured at the $k^{\text{th}}$ time period (or, equivalently, pointing angle), it is clear that

where $a_{k}(\nu)=\int B_{k}(\nu,\theta,\phi)\ \cos{\psi}\ d\Omega$.¹⁵¹⁵15Note that $a_{k}(\nu)$ will not be known exactly because the beam will not be known exactly. However, since the dipole component will be on the order of 10 mK (see Slosar, 2017) and $a_{k}(\nu)$ will be of order unity, the expected imprecision levels in the beam, which should be at the sub-percent level, should not impact the results significantly. This can be written in matrix notation where frequency channels correspond to the elements of vectors. In this way, the signature of the 21-cm signal in the $k^{\text{th}}$ spectrum is

where ${\boldsymbol{A}}_{k}$ is the diagonal matrix with nonzero elements given by ${\boldsymbol{a}}_{k}$, ${\boldsymbol{N}}$ is the diagonal matrix with the frequencies of the data channels on the diagonal, and ${\boldsymbol{D}}$ is a derivative matrix (see Appendix E for subtleties involved in taking derivatives with a matrix). Now, combining all $N_{s}$ spectra into one vector leaves ${\boldsymbol{t}}={\boldsymbol{\Psi}}{\boldsymbol{s}}$ where

To extend the MAA to include the motion-induced dipole, one must merely plug this ${\boldsymbol{\Psi}}$ into Equations 19a and 19b. The effect of including the motion-induced dipole on the results of the MAA is left for future work.