A Space-based Observational Strategy for Characterizing the First Stars and Galaxies Using the Redshifted 21 cm Global Spectrum

One of the last frontiers of observational cosmology is the time period stretching from the end of the dark ages through cosmic dawn (≈80−500 million years after the big bang). This is a virtually unobserved yet key epoch in the early universe. During this interval, the first luminous objects including stars, galaxies, and accreting black holes "turned on" (e.g., Loeb & Furlanetto 2013). Furthermore, this time period saw the birth of structural complexity in the universe. At the beginning of the dark ages, corresponding to the epoch of recombination, the universe was smooth to 1 part in 10⁵ as evidenced by the cosmic microwave background (CMB; e.g., Mather et al. 2013). However, less than a billion years later, the universe was teeming with complex structure. Thus, this transition time in the universe is vital to understanding how the core components and structures of today's universe came to be.

The highly redshifted 21 cm spectral line of neutral hydrogen, produced by a spin-flip hyperfine transition (Ewen & Purcell 1951; Field 1958), provides an observable window into the early universe's intergalactic medium (IGM) before the epoch of reionization (EoR) was complete (Madau et al. 1997). The heating and ionization caused by the "first objects to light up the universe"¹⁴ serve as indirect probes of the nature of the first stars and galaxies. With an effective optical depth of ≈1% and sensitivity to low temperatures, the resulting signal measured against the CMB permits us to investigate a large evolutionary range from the dark ages through the end of the EoR (e.g., Furlanetto et al. 2006; Morales & Wyithe 2010; Pritchard & Loeb 2012).

The 21 cm all-sky or global signal (Shaver et al. 1999; Pritchard & Loeb 2010; Figure 1) is an attractive observational target for either a single antenna (e.g., Bowman & Rogers 2010; Burns et al. 2012; Patra et al. 2013; Voytek et al. 2014; Bernardi et al. 2015; Sokolowski et al. 2015a) or a small, compact array of antennas (e.g., Presley et al. 2015; Singh et al. 2015; Vedantham et al. 2015). Features in the spectrum may provide the first constraints on the birth and nature of the first luminous objects (e.g., Furlanetto 2006). Such an experiment for 21 cm cosmology is analogous to the COBE measurement of the CMB blackbody spectrum, which set the stage for more detailed studies of spatial fluctuations by WMAP and Planck.

Zoom In Zoom Out Reset image size

In this paper, we describe a space-based strategy for observations of the 21 cm global signal that probes the time of formation and the characteristics of the first stars and galaxies. We demonstrate how signal extraction using a realizable radiometer system and Bayesian statistical techniques, in the presence of strong galactic and extragalactic foregrounds, can measure spectral features and constrain the physical properties of the first luminous objects. We use the new detailed design of the Dark Ages Radio Explorer (DARE) to illustrate how the 21 cm spectrum can be extracted from the foreground using a feasible observational strategy.

DARE is proposed to conduct observations between 40 and 120 MHz in an orbit around the Moon with data taken only above the lunar farside. On Earth, the ionosphere corrupts low-frequency observations (see, e.g., Vedantham et al. 2014; Rogers et al. 2015; Sokolowski et al. 2015b; Vedantham & Koopmans 2015; Datta et al. 2016, and references therein) due to refraction, absorption, and emission driven by solar emissions and the solar wind (Liu et al. 2011; Polygiannakis et al. 2003). At 50 km above the lunar farside, >90 dB of radio frequency interference (RFI) attenuation produces an environment quiet to <1 mK (e.g., McKinley et al. 2013). In addition, the Moon shields the instrument (about half the time) from variable solar emission caused by flares and coronal activity (e.g., Mercier & Trottet 1997). Therefore, observations above the night-time, pristine, radio-quiet lunar farside (as verified by RAE-2, Alexander & Kaiser 1976) bypass the challenges presented by the Earth and the Sun and provide an optimal site for measurements of the global 21 cm signal.

The key insight permitting the cosmic dawn signal to be detected in the presence of bright foregrounds is that once the Moon blocks solar effects and terrestrial RFI, the foregrounds are significantly different in their characteristics from the expected 21 cm spectral signal. The 21 cm monopole strength is about four orders of magnitude weaker than the Galactic foreground. However, the 21 cm signal is separable from the foreground because it is spatially uniform at angular scales ≳10° (e.g., Bittner & Loeb 2011; Liu et al. 2013), unpolarized, and has distinct spectral features, whereas the observed foreground varies spatially, exhibits polarized emission, and is spectrally featureless. The 21 cm cosmological signal can then be extracted using algorithms similar to those employed for CMB observations implemented via a Markov Chain Monte Carlo framework (Harker et al. 2012, 2016).

The paper is organized as follows. In Section 2, we introduce and summarize the space-based observational strategy. In Section 3, an overview of the stellar models for the sky-averaged 21 cm signal used to develop the observational strategy is presented. Section 4 describes the nature and brightness of astronomical foregrounds, which must be considered in efforts to measure the much weaker Cosmic Dawn signal. Section 5 provides a synopsis of the new design for DARE. Section 6 describes our software pipeline for signal extraction. In Section 7, we discuss the physical parameters (and their uncertainties) associated with the first stars, black holes, and galaxies that are expected to be measured using the 21 cm all-sky spectrum. Section 8 presents a summary of the potential use of the 21 cm background to detect the first luminous objects in the early universe.

Here we briefly describe key aspects of our observational strategy. The following sections will provide details about each item, as well as their relevance to the overall tactic. The core components of the strategy are as follows.

• 1.  
  We incorporate a wide range of theoretical models ($\gt 1.5\times {10}^{4}$) from two different classes of possible signals, differing by the generation of stars whose contribution dominates the behavior of the signal (Pop II or Pop III). We show that the DARE instrument in orbit of the Moon can effectively differentiate between these models using our Bayesian inference pipeline.
• 2.  
  We realistically model the diffuse foregrounds accounting for spatial variations of their spectral index, which is estimated from all-sky, publicly available maps at two frequencies (45 and 408 MHz).
• 3.  
  The new DARE reference instrument design incorporates (1) an optimized antenna with on-orbit beam calibration, (2) the replacement of Dicke switches for bandpass calibration with a pilot frequency tone system capable of high dynamic range monitoring of gain variations and measurements of the system reflection coefficients, and (3) polarimetric observations to provide a model-independent measure of the beam-averaged foregrounds. The observations, performed from the radio-quiet zone above the Moon's farside, will be enabled through a unique "frozen" 50 × 125 km lunar equatorial orbit (Plice et al. 2017). The nominal observation time corresponds to 800 hr, which results in radiometric noise integration to the 1.7 mK level at 60 MHz. The instrument provides data with the frequency range (40–120 MHz), spectral resolution (50 kHz), beam characteristics ($\approx 60^\circ$ FWHM at 60 MHz), and polarization required to measure the spectral features expected from the wide range of theoretical models considered. DARE's present-day observing strategy utilizes four quiet-sky pointing directions away from the galactic center.
• 4.  
  Our signal extraction pipeline is centered around a singular value decomposition (SVD) approach, which allows us to robustly separate the 21 cm signal from the additional contributions to the measurement by using orthogonal modes of variation of each component. These modes are determined from well-characterized training sets constructed from either theory or measurements.

We constructed a detailed end-to-end observation model that generates simulated antenna temperatures using our models for the diffuse foregrounds and the 21 cm spectrum, and the predicted telescope pointing, Moon location, and instrument characteristics. Unlike previous papers (e.g., Mirocha et al. 2015; Harker et al. 2016), which assumed perfect knowledge of the instrument, this new process accounts for and propagates the uncertainty in the instrumental parameters to the signal extraction pipeline. Our instrument sensitivity metric is defined as the rms uncertainty of the extracted 21 cm spectrum, averaged over the observation band. Our requirement for this metric is to keep it below 20 mK for all models processed with our pipeline.

In this section, we discuss the global 21 cm signal and describe the broad set of physical models that are incorporated into our analysis strategy.

The 21 cm global signal arises from the radiation effects produced by the first stars, accreting black holes and galaxies on the surrounding IGM. X-ray and UV emission from these objects and their descendants heated and ionized the tenuous gas that lies between galaxies, culminating in the EoR several hundred million later. The 21 cm background can be used to measure these radiation effects with the hyperfine line of the neutral hydrogen (H i) gas pervading the universe. The expansion of the universe redshifts these photons from earlier epochs to lower observed frequencies, $\nu =1420/(1+z)$ MHz (e.g., at z = 30, ν = 45 MHz). Importantly, this frequency–redshift relation enables a direct reconstruction of the history of the universe as a function of time from the 21 cm spectrum.

Figure 1 shows some example predictions (among those currently allowed) for the 21 cm spectrum during the dark ages and cosmic dawn. The brightness temperature of this 21 cm signal is given by (e.g., Madau et al. 1997; Shaver et al. 1999; Furlanetto et al. 2006)

$$\begin{eqnarray}\delta {T}_{b} & \simeq & 27{x}_{{\rm{H}}{\rm{I}}}\left(\displaystyle \frac{{T}_{s}-{T}_{\mathrm{CMB}}}{{T}_{s}}\right){\left(\displaystyle \frac{1+z}{10}\right)}^{1/2}\\ & & \times (1+{\delta }_{B}){\left[\displaystyle \frac{{\partial }_{r}{v}_{r}}{(1+z)H(z)}\right]}^{-1}\mathrm{mK},\end{eqnarray} \tag{ 1 }$$

where ${x}_{{\rm{H}}{\rm{I}}}$ is the fraction of neutral gas, T_s is the 21 cm spin temperature, ${T}_{\mathrm{CMB}}$ is the CMB temperature, ${\delta }_{B}$ is the baryon overdensity (taken here to be ${\delta }_{B}\sim 0$), and H(z) is the Hubble parameter. The last term in this equation includes the effect of the peculiar velocities with line-of-sight velocity derivative ${\partial }_{r}{v}_{r}$. Since we will measure the spatially averaged $\delta {T}_{b}$, the effects of the last term in Equation (1) are negligible for observations of the 21 cm global signal (e.g., Bharadwaj & Ali 2004; Barkana & Loeb 2005).

Several important physical processes drive the evolution of $\delta {T}_{b}$ with redshift. These include (1) UV radiation from the first stars, which "activates" the spin-flip signal through the Wouthuysen–Field mechanism (Wouthuysen 1952; Field 1958); (2) X-ray heating, likely generated by gas accretion onto the first black holes; and (3) ionizing photons from the first galaxies (which destroy the neutral hydrogen).

The relevant radiation backgrounds grow at different times, so their interplay creates distinct features in the spectrum (Furlanetto et al. 2006; Pritchard & Loeb 2010; Mesinger et al. 2011). When the first stars appear, their UV radiation drives T_s toward the cold temperatures that are characteristic of IGM gas ($z\sim 35\mbox{--}22$ across our range of models; Region B in Figure 1), triggering a deep absorption trough (Madau et al. 1997). Shortly after, black holes likely formed, e.g., as remnants of the first stars ($z\sim 25\mbox{--}12$ across our range of models; Region C). The energetic X-ray photons from these accreting black holes travel great distances, eventually ionize H and He atoms, and produce photo-electrons that deposit some of their energy as heat in the IGM (Shull & van Steenberg 1985; Furlanetto & Johnson Stoever 2010), transforming the 21 cm signal from absorption into emission as the gas becomes hotter than the CMB (Region D). The emission peaks as photons from these stars and black holes ionize the IGM gas ($z\lt 12$), eventually eliminating the spin-flip signal.

The dashed red curve in Figure 1 assumes that the efficiency and properties of star formation in early galaxy populations (Sun & Furlanetto 2016; Mirocha et al. 2017) and the relationship between X-ray luminosity and star formation rate are the same as at later times (Mineo et al. 2012b). There are several reasons to expect that this Pop II model is conservative, i.e., that it underestimates the total production rate of UV and X-ray photons. For example, it assumes solar metallicity, though stars in high-z galaxies are likely forming in metal-poor environments, which can boost their UV (Eldridge & Stanway 2009) and X-ray outputs (Brorby et al. 2016). The solid red curve in Figure 1 assumes that galaxies have metallicities (Z) 5% of solar, which results in a shallower absorption feature due to enhanced X-ray emission (assuming the Brorby et al. 2016 L_X-SFR-Z relation). Alternatively, the black curves include a simple model for Pop III stars, in which low-mass halos (below atomic cooling threshold) can produce UV and X-ray photons (neglected by red curves). Boosts of 100 in the efficiency of the UV (dashed black) and also the X-ray luminosity (solid black) of Pop III stars relative to Pop II result in a variety of qualitatively different predictions for the global 21 cm signal. Pop III models that resemble our black curves should be relatively straightforward to distinguish from Pop II-only models for an experiment like DARE (see Section 6 and Figure 7). At this stage, our ability to label each set of curves as being Pop II- or Pop III-dominated assumes that the current Pop II models are well calibrated (Mirocha et al. 2017). New measurements (by, e.g., JWST) in the coming years can be immediately incorporated into the model and will act to mitigate degeneracies between Pop II and Pop III sources. More subtle features of the signal, such as its asymmetry, may also reveal the presence of Pop III despite uncertainties in the calibration of Pop II models (J. Mirocha et al. 2017, in preparation).

It is also worth noting that the 21 cm global signal traces the collective effects of all sources in the redshift ranges illustrated in Figure 1, which form a mostly unresolved sea of fainter objects that likely dominate the total emissivity of the early universe. The red curves in Figure 1 are calibrated to match the latest luminosity function measurements from HST (which probe relatively bright galaxies that can be resolved) and CMB optical depth (${\tau }_{e}$) measurements from Planck (Mirocha et al. 2017), with variations arising solely due to differences in the adopted properties of galaxies beyond the current detection threshold. JWST and future CMB missions will further constrain the bright-end of the galaxy luminosity function and ${\tau }_{e}$, respectively, and will thus enhance the sensitivity of the 21 cm global signal to Pop III stars and their remnants in faint galaxies.

The signal models described in this section are used to create the signal training set, an essential component of our observational strategy from which the signal extraction pipeline calculates the main modes of variation of the signal. The new data expected from JWST and from CMB missions will constrain parameter space, which will allow us to restrict the training set and reduce parameter degeneracies and covariances. See Section 6 for more details on the training set and its effect on the uncertainty of our signal estimate.

Here, we discuss the origin and properties of the foregrounds expected in the 21 cm measurement from lunar orbit, which are modeled and accounted for in our signal extraction pipeline.

4.1. Galaxy/Extragalactic Foregrounds

Beam-averaged diffuse sky foregrounds represent the strongest contributors to any highly redshifted 21 cm measurement for a space-based experiment. The most important arises from our Galaxy (Shaver et al. 1999). In addition, a "sea" of extragalactic sources appear as another diffuse, spectrally featureless power-law foreground (at DARE's resolution) and contributes ∼10% of the total sky brightness temperature (Figure 2). The emission from these foregrounds is produced by synchrotron radiation that intrinsically has a smooth frequency spectrum (e.g., Petrovic & Oh 2011; Bernardi et al. 2015). On top of the spectral smoothness, the foregrounds are spatially variable (inset in Figure 2). Their featureless spectrum and spatial variability contrast with the spectral features and spatial uniformity of the 21 cm spectrum, making them separable (Liu et al. 2013; Switzer & Liu 2014).

Zoom In Zoom Out Reset image size

Theoretical models predict that the foreground is well approximated by a third-order polynomial to levels below the amplitude of the 21 cm global signal, especially over low-foreground regions (Bernardi et al. 2015). Smoothness over a frequency range much broader than ∼40–120 MHz is supported by sky models produced from measurements that cover the range of 10 MHz–5 THz (de Oliveira-Costa et al. 2008; Sathyanarayana Rao et al. 2017; Zheng et al. 2017). These models rely on, at most, five components to describe the spectral content of the foreground over several decades in frequency. Global measurements from the Experiment to Detect the Global EoR Signature (EDGES), Sonda Cosmológica de las Islas para la Detección de Hidrógeno Neutro (SCI-HI), Shaped Antenna measurement of the background RAdio Spectrum (SARAS), and Large-Aperture Experiment to Detect the Dark Ages (LEDA) provide further validation of the intrinsic foreground smoothness (Rogers & Bowman 2008; Voytek et al. 2014; Patra et al. 2015; Bernardi et al. 2016; Mozdzen et al. 2017).

In our strategy, we use a diffuse foreground model produced from all-sky observations taken at two frequencies, 45 and 408 MHz (Haslam et al. 1982; Guzmán et al. 2011), in order to account for spatial variations in the spectral index.

The spectrally smooth foreground is altered via the frequency-dependent antenna response (Vedantham et al. 2014; Bernardi et al. 2015; Mozdzen et al. 2016). The beam directivity of finite-sized, wideband antennas does not remain constant across frequency (Rumsey 1966). This beam "chromaticity" impacts the observed spectrum of the spatially dependent foregrounds. The variation with frequency of the beam shape and directivity imprints spectral structure into the beam-averaged response that is not intrinsic to the foregrounds.

As part of our strategy, chromaticity is addressed by minimizing instrumental design effects and making precise beam measurements on the ground and on-orbit. We also estimate the beam chromaticity by modulating the beam-averaged foregrounds through rotation of the antenna about the boresight axis. This technique is discussed in Section 5. This represents a significant advancement over previous simulations. For instance, in Harker et al. (2016), the beam was taken to be Gaussian and the integrated foreground was assumed to perfectly follow a polynomial of the form $\mathrm{ln}{(T)={\sum }_{i=0}^{5}{a}_{i}\mathrm{ln}(\nu )}^{i}$.

Finally, we note that the low-foreground areas of the sky are polarized to a few percent ($\lesssim 5 \%$) (Jelić et al. 2014, 2015; Lenc et al. 2016). Our dual polarization instrument directly measures this intrinsic sky polarization. At the same time, this polarization is minimized through dilution produced by our wide antenna beam, and also averaged down by our scanning strategy, which includes antenna rotation.

4.2. Other Foregrounds

To illustrate how the cosmological 21 cm spectrum can be extracted from the foregrounds, we use the new science instrument proposed for DARE (Figure 3). In Burns et al. (2012), we outlined a very basic approach for idealized, lunar-based 21 cm cosmology observations. We have now advanced the fidelity of the instrument model to evaluate measurements of the spectrum and constrain parameters for the first luminous objects at the significance level presented in Section 2, in the presence of realistic and well modeled uncertainties.

Zoom In Zoom Out Reset image size

Figure 4 shows a block diagram of the current DARE instrument design, which consists of four subsystems: (1) an antenna composed of a pair of crossed biconical dipoles above a ground plane that provides dual polarization with low reflection coefficient (−12 dB average across the band) and beam chromaticity (the beam directivity spectral knowledge goal is ∼20 ppm, see below); (2) a thermally controlled receiver with a calibration architecture that utilizes precise continuous-wave frequency tones optimized to yield a frequency response that meets DARE's rms criterion; (3) a spectrometer with a wide bandwidth and digital receiver that provides the spectral resolution and Stokes processing of the V and H channels; and (4) an instrument electronics subsystem to interface with the spacecraft. The expectation for the instrument envelopes the hardware performance of systems on the ground (e.g., EDGES, Bowman & Rogers 2010; Cosmic Twilight Polarimeter, CTP, Nhan et al. 2017) and in space (Global Precipitation Measurement Microwave Imager, GMI¹⁵ ).

Zoom In Zoom Out Reset image size

We model the forward instrument response following that used by EDGES (Monsalve et al. 2017) as:

$$\begin{eqnarray}&&P=g[| F{| }^{2}({\eta }_{l}{T}_{{\rm{A}}}+(1-{\eta }_{l}){T}_{\mathrm{Ap}})(1-| {{\rm{\Gamma }}}_{{\rm{A}}}{| }^{2})+{T}_{\mathrm{off}}],\end{eqnarray} \tag{ 2 }$$

where P is the raw power measured by the instrument, g and ${T}_{\mathrm{off}}$ represent the system gain and radiometric offset, respectively, ${\eta }_{l}$ accounts for the antenna and balun losses at physical temperature T_Ap, $1-| {{\rm{\Gamma }}}_{{\rm{A}}}{| }^{2}$ accounts for the reflection coefficient of the antenna, and $| F{| }^{2}$ is the throughput of the receiver front end accounting for multiple reflections between the receiver and antenna. The instrument calibration activities consist of using ground, on-board, and on-orbit calibration to invert the forward instrument response model and provide an estimate of the antenna temperature T_A.

During science observations, the receiver is calibrated continuously using the pilot tone injection receiver architecture. The calibration system generates tones at $\approx 5$ frequencies simultaneously across the band to adequately sample the frequency range. The tones are each within a single 50 kHz spectrometer bin, and thus produce negligible degradation in spectral performance. The nominal calibration cycle consists of a sequence of four states, which are enabled for 10 s each: (1) high-level tones directed toward the receiver, (2) low-level tones directed toward the receiver, (3) high-level tones directed toward the antenna, and (4) low-level tones directed toward the antenna. The gain of the receiver is computed by differencing instrument-measured power from the high- and low-level injected tones toward the receiver divided by the difference in effective input brightness temperature of the tones characterized during pre-flight calibration. Likewise, the tones injected toward the antenna afford an on-board measurement of the antenna reflection computed in a similar fashion.

The terms ${T}_{\mathrm{off}}$, ${\eta }_{l}$, and $| F{| }^{2}$ in Equation (2) are computed based on ground measurements and the on-orbit trending of the receiver gain and reflection coefficient.

The DARE antenna and receiver are designed to minimize temperature variations by limiting exposure to the solar flux and lunar albedo. For the antenna, thermal baffles, as shown in Figure 3, result in a predicted physical temperature change over the lunar orbit of $10^\circ {\rm{C}}$ allowing DARE to maintain a nearly constant beam directivity. The front-end receiver includes a proportional-integral-derivative temperature control to provide predicted thermal stability of $0.1^\circ {\rm{C}}$, thus reducing receiver systematics to meet DARE's calibration and stability requirements.

A novel feature of the current design is on-orbit calibration of the antenna directivity. Measurements of the beam will be obtained by receiving $\geqslant 3$ narrow-band, circularly polarized signals spaced across the band, transmitted from a large antenna on the Earth as DARE orbits the Moon above the nearside. The spacecraft (and antenna) is slowly rotated while it continuously receives these signals. The received signal power at each frequency as a function of antenna pointing will produce a slice through the beam power pattern. The transmitted signals will also reflect off the lunar regolith and return to the same antenna on Earth to correct for ionospheric effects. The in situ beam measurement system is currently being baselined to use the 140-foot radio telescope at the Green Bank Observatory, operating with 50% aperture efficiency, transmitting 10 kW of power, and using 10 s averaging.

Another innovation in the current design of DARE is polarization measurements to constrain and distinguish the beam-averaged foregrounds from the unpolarized H i signal (Nhan et al. 2017). Our observation strategy incorporates rotation of the antenna about the boresight axis to modulate the signals captured by the two polarization arms of the antenna. This modulation results in induced polarization that tracks and measures the beam-averaged foreground spectrum, without relying on, e.g., polynomial model fits, and is insensitive to the spatially uniform 21 cm signal. CMB polarization measurements use analogous modulation approaches, achieving stability and systematic control required for $\mu K$ polarimetric sensitivity (e.g., Bennett et al. 2003; Bischoff et al. 2013).

The on-orbit measurements of the antenna directivity and the induced polarization technique enable us to anticipate a knowledge of the beam-averaged foregrounds at a level of ∼20 ppm. This represents an important advancement that allows us to achieve our goal of <20 mK rms spectral uncertainty on the extracted 21 cm models.

In this section, we demonstrate an essential aspect of our observational strategy: how our data analysis pipeline is able to separate the 21 cm signal from foregrounds measured through a realizable instrument. We model each of these components, signal along with the foreground and instrument systematics, as described below.

In our previous work (Harker et al. 2012, 2016, hereafter H12 and H16), we developed a foundation for a 21 cm signal extraction pipeline using a Markov Chain Monte Carlo (MCMC) framework. However, we assumed an idealized instrument with exact knowledge of most instrument systematics and the form of the beam-averaged foreground. We have now significantly expanded the initial analyses of H12 and H16 by implementing a robust SVD modeling scheme based upon a pragmatic end-to-end instrument model (Section 5). Specifically, the current pipeline (which will be released to the community in a later publication) incorporates the following aspects for the first time.

• 1.  
  Full simulations of the antenna beam-weighted foreground. These simulations are based upon beam patterns calculated by the CST electromagnetic simulation package¹⁶ and our diffuse foreground model, described in Section 4.
• 2.  
  A calibration model, based upon expected lab measurements and uncertainties, that includes all parameters in Equation (2). The instrument model described in H12 included only the antenna reflection coefficient.
• 3.  
  A modeling scheme, detailed below, based upon the implementation of SVD on well-characterized training sets for both the signal and a complete set of instrument and foreground systematics. The SVD technique independently determines the main modes of variation in the signal and systematics. The MCMC algorithm then simultaneously fits all the coefficients associated with the SVD modes to extract the signal. This is a major improvement over our previous use of polynomials (Fourier series) to fit the foreground (reflection coefficient).

The MCMC algorithm in the pipeline samples the likelihood function

$$\begin{eqnarray}&&\mathrm{ln}L(\gamma )=-\displaystyle \frac{1}{2}\displaystyle \sum _{r=1}^{{N}_{r}}\displaystyle \sum _{i=1}^{{N}_{\nu }}{\left[\displaystyle \frac{{T}_{{\rm{A}},{\rm{D}}}^{(r)}({\nu }_{i})-{T}_{{\rm{A}},{\rm{M}}}^{(r)}({\nu }_{i},\gamma )}{{\sigma }_{r}({\nu }_{i})}\right]}^{2},\end{eqnarray} \tag{ 3 }$$

where ${T}_{{\rm{A}},{\rm{D}}}^{(r)}(\nu )$ and ${T}_{{\rm{A}},{\rm{M}}}^{(r)}(\nu )$ are the antenna temperature spectra for the data (D) and the model (M), r and ν indicate the sky direction and frequency channel, respectively, and ${\sigma }_{r}(\nu )\,={T}_{{\rm{A}},{\rm{D}}}^{(r)}(\nu )/\sqrt{{\rm{\Delta }}\nu {\rm{\Delta }}t}$ is the thermal noise level in the data for a given frequency bin of width ${\rm{\Delta }}\nu$ centered on ν integrated over a time ${\rm{\Delta }}t$.

We model ${T}_{{\rm{A}},{\rm{M}}}^{(r)}(\nu )$ as a linear combination of (dimensionless) principal modes derived from SVD analyses (Paciga et al. 2013; Switzer & Liu 2014; Vedantham et al. 2014),

$$\begin{eqnarray}&&{T}_{{\rm{A}},{\rm{M}}}^{(r)}(\nu ,\gamma )=\displaystyle \sum _{i=1}^{n}{({\gamma }_{21})}_{i}{f}_{i}(\nu )+\displaystyle \sum _{j=1}^{m}{({\gamma }_{\mathrm{sys}})}_{j}^{(r)}{g}_{j}(\nu ),\end{eqnarray} \tag{ 4 }$$

where ${f}_{i}(\nu )$ and ${g}_{j}(\nu )$ are the SVD signal and systematic modes, respectively, and ${({\gamma }_{21})}_{i}$ and ${({\gamma }_{\mathrm{sys}})}_{j}$ (both with units K) are the coefficients associated with each of them. We fit the entire parameter space, $\gamma =[{\gamma }_{21},{\gamma }_{\mathrm{sys}}]$, simultaneously (using the emcee code; Foreman-Mackey et al. 2013) in order to account for the covariance between all parameters and ensure self-consistency. This MCMC calculation efficiently and robustly obtains the full posterior distribution.

In this work, we utilize n = 6 (signal) and m = 7 (systematic) SVD modes because they are able to fit our fiducial models to within the thermal noise level achieved through 800 hr of integration. For future analyses, however, we are developing a novel technique that will choose the optimal number of modes to use in the pipeline. The details of this key advancement will be described in forthcoming works (D. Tauscher et al. 2017, in preparation; K. Rapetti et al. 2017, in preparation).

The systematic modes ${g}_{j}(\nu )$ are derived from 10,000 simulated data sets, which vary the foreground and instrument within expected uncertainties. This process utilizes Equation (2), its inverse, and the fiducial values of the calibration and beam-weighted foreground parameters (D. Tauscher et al. 2017, in preparation). Currently, the signal modes ${f}_{i}(\nu )$ are derived from input training sets of 21 cm spectrum simulations (15,000 and 960, respectively) based on two well-motivated ranges of physical models (primordial Pop II and Pop III stars; see Section 3).¹⁷ In future work, the signal models will be combined into a single training set.

The Bayesian nature of the MCMC permits the incorporation of key prior knowledge on the instrument calibration and foregrounds when retrieving the posterior probability distribution of the model parameters. In the instrument simulations, we account for all the identified uncertainties and priors, including a 50 mK prior on the beam-averaged foregrounds from measurements of the induced polarization. Even though, at this stage, the induced polarization is used only as a prior on the antenna temperature, T_A. In future work, all four Stokes parameters will be included in the likelihood function.

Figure 5 shows the SVD modes used in this work. The left and middle panels contain the signal modes for the models of primordial Pop II and Pop III stars, respectively. For the purpose of reducing the covariance between our parameters, γ, it is important that the SVD systematic modes (right panel) are as orthogonal as possible (i.e., have a minimal dot product) with the SVD signal modes (left or middle panel). Note that the order of the modes is color-coded in the same way for each panel to ease the comparison of modes. The extracted signals in Figure 7 show that the shapes of the modes are sufficiently different to enable a clean extraction of the signal.

Zoom In Zoom Out Reset image size

Figure 6 shows the covariance matrix of the six signal parameters, ${({\gamma }_{21})}_{i}$ with $i\in \{1,2,\ldots ,\,6\}$, and the seven systematic parameters, ${({\gamma }_{\mathrm{sys}})}_{j}$ with $j\in \{1,2,\ldots ,\,7\}$. The top left corners within each of the four regions in both panels of Figure 6 demonstrate that the lowest order signal modes have enough similarities in shape with the first three to four systematic modes to generate only modest covariances.

Zoom In Zoom Out Reset image size

By simultaneously fitting all parameters, γ, and marginalizing over the systematic parameters, ${\gamma }_{\mathrm{sys}}$, we are able to clearly separate the signal from the systematics despite those covariances, as shown in Figure 7. The widths of the uncertainty bands result from the covariances between the signal parameters. These covariances depend upon the overlap between signal modes and all other modes. By construction, the signal modes are orthogonal among themselves, but there is no enforced orthogonality between them and the systematic modes. This leads to covariances between pairs of signal modes, which overlap with the same systematic mode.¹⁸ These details will be explored in an upcoming work (D. Tauscher et al. 2017, in preparation).

Zoom In Zoom Out Reset image size

In summary, for the DARE instrument parameters discussed in Section 5 and 800 hr of total integration above the lunar farside, our signal extraction pipeline recovers the spectra and uncertainties for two representative models (Pop II and Pop III star models) shown in Figure 7. In addition to the 21 cm signal, the pipeline simultaneously fits the SVD modes of the receiver, beam, and foreground utilizing prior information and on-orbit measurements. With an average rms of $\approx 17$ mK, we recover the major features in the spectra and differentiate between different stellar population models.

With the calibrated spectra and uncertainties in Figure 7, it is straightforward to estimate when the first luminous objects ignited and began reionizing the universe. Since redshift maps directly to frequency, measurements of the extrema frequencies from the 21 cm spectrum determine when major events occurred. The frequency of the Region B extremum (${\nu }_{B}$) determines the z at which the UV background activates the 21 cm transition (i.e., first star's ignition). This clean and accurate measurement delineates the nature of the first stars, especially considering that no observational bounds currently exist. Using a Pop III model as a working example, DARE will extract ${\nu }_{B}$ with a 1% (0.4 MHz) uncertainty (68% confidence). Similarly, the redshift when the first black holes began accretion is measured from the Region C extremum frequency (${\nu }_{C}$) with 1% (0.6 MHz) uncertainty. The redshift of the beginning of reionization is measured from the extremum ${\nu }_{D}$ with 2% (2 MHz) uncertainty. Different models (Figure 1) yield similar uncertainties for the extrema frequencies.¹⁹

As long as signal extraction is performed with a sufficiently flexible parametric function, the extrema frequencies—and thus the timing of major events—can be recovered in a fairly model-independent fashion (Harker et al. 2016). However, extracting information about the characteristics of the first stars and galaxies, along with the history of the early universe, require a physical model of the calibrated spectrum. We use the models of Mirocha et al. (2017) to determine the ionization history, thermal history, and properties of galaxies required to match the calibrated spectrum.

The history of reionization in the early universe is characterized by the H i fraction (${x}_{{\rm{H}}{\rm{I}}}$) and the IGM kinetic temperature (T_K) at $z\sim 11$. Our modeling of features in Region D using DARE's sensitivity yields uncertainties of 5% and 10% for ${x}_{{\rm{H}}{\rm{I}}}$ and T_K, respectively.

Next, the features in the 21 cm spectrum at the lowest frequencies depend upon the stellar populations that dominate the UV background; if, for example, Pop III star formation is efficient, we should expect features of the signal to occur at lower frequencies (higher redshift) than if Pop II stars dominate the background because Pop II stars form in more massive halos, which do not become abundant until relatively late times (low redshift). DARE's sensitivity can separate the effects of the broad classes of Pop II and Pop III stellar models considered in this work (see Figures 1 and 7), subject to the assumed calibration of the Pop II contribution (see Section 3; Mirocha et al. 2017) and the model for Pop III stars. A useful metric for gauging the influence of Pop III is the ratio of UV production efficiencies for Pop III compared to Pop II stars, ${\xi }_{\alpha }\equiv {\xi }_{\alpha ,\mathrm{III}}/{\xi }_{\alpha ,\mathrm{II}}$. The value of ${\xi }_{\alpha ,\mathrm{II}}$ is drawn from the BPASS models (Eldridge & Stanway 2009), assuming solar metallicity, while ${\xi }_{\alpha ,\mathrm{III}}$ is allowed to vary freely. DARE constrains ${\xi }_{\alpha ,\mathrm{III}}$ with 25% uncertainty in Figure 8.

Zoom In Zoom Out Reset image size

The characteristics of the first X-ray sources (Region C in Figure 1) are inferred from the ratio of X-ray heating efficiencies between Pop III and Pop II stars. Analogous to the UV constraints, the Pop III X-ray efficiency, ${\xi }_{{\rm{X}},\mathrm{III}}$, is allowed to vary freely, while ${\xi }_{{\rm{X}},\mathrm{II}}$ is anchored to the local relation between X-ray luminosity and SFR (Mineo et al. 2012a) assuming high-mass X-ray binaries are the dominant source. DARE can measure ${\xi }_{{\rm{X}},\mathrm{III}}$ with 15% uncertainty (see Figure 8). Further modeling plus multi-wavelength observations (e.g., the cosmic X-ray background; Fialkov et al. 2017) may help to better constrain the identity of the universe's first X-ray sources, whether they be black hole X-ray binaries, hot gas in star-forming galaxies, or proto-quasars.

Finally, before concluding, we emphasize that these Pop III models are quite simple, as, for example, they neglect an explicit treatment of feedback. As a result, the interpretation of the precise value of ${\xi }_{\alpha ,\mathrm{III}}/{\xi }_{\alpha ,\mathrm{II}}$ may be considerably more complex in practice, but the finding that both values are non-zero is robust.

To achieve the science potential of 21 cm global spectral observations, we proposed an observational strategy that carefully considers the local environment, the instrument, and the methods for signal extraction. A lunar-orbiting experiment above the Moon's farside has the best probability of measuring the 21 cm spectrum since this environ is free of ionospheric effects and human-generated radio frequency interference.

Signal extraction in the presence of bright foregrounds is the greatest challenge for all observations of the 21 cm cosmological spectrum. Utilizing singular value decomposition to model the foreground and instrument along with a Markov Chain Monte Carlo numerical inference technique to survey parameter space, we showed that it is possible to accurately recover the expected features in the spectrum in the presence of bright foregrounds with the instrument characteristics of the Dark Ages Radio Explorer (DARE) for $\approx 800\,\mathrm{hr}$ of integration. To separate the signal from the foreground, the antenna system must be well-characterized requiring temperature control and precise beam directivity measurements on the ground and in space. In addition, a model-independent constraint on the foreground from polarimetric observations is an important element in the signal extraction.

From the extracted 21 cm spectrum (including confidence intervals), we showed that meaningful constraints can be placed upon the physical parameters of primordial radiating objects. The redshift for the commencement of first star formation and X-ray emission from the first accreting black holes along with the redshift of the beginning of reionization can be inferred to within ≈a few percent. The 21 cm signal is also uniquely sensitive to the different radiation effects produced by Pop II and Pop III stellar models. Specifically, the UV production and X-ray heating efficiencies can be constrained, thus determining which stellar population was dominant within the first galaxies. Finally, the history of reionization in the early universe can be characterized by the redshift evolution of the H i ionization fraction (x_{H i}) inferred from the 21 cm spectrum.

Accurate parameter estimation is a core capability required for 21 cm global signal observations and interpretations. Bayesian methods have significant potential for 21 cm observations (Bernardi et al. 2016; Greig et al. 2016). They have proven to be successful for other experiments targeting weak signals, including CMB observations (Planck Collaboration et al. 2016a, 2016b) and the LIGO gravitational wave detections (Veitch et al. 2015; Abbott et al. 2016). The next step in the analyses of the global 21 cm spectrum is to construct a likelihood function allowing differentiation between differing physical models for the first halos. Similarly, modeling different levels of structure in the beam-convolved foregrounds needs a refined Bayesian approach. This is a highly computationally intensive process. We are refining and extending our SVD modeling approach toward these goals. In addition, recent developments of nested sampling algorithms for high dimensional parameter spaces, which operate in massively parallel computer architectures (Handley et al. 2015) have great potential for 21 cm cosmology applications.

In conclusion, measurements of spectral features in the 21 cm spectrum will answer key science questions from the NRC Astrophysics Decadal Survey: "What were the first objects to light up the universe and when did they do it?" With a clean measurement within the radio-quiet environs of the Moon's farside and proven technology, the 21 cm global signal will open a new window of discovery into the early universe.

This research was supported by the NASA Ames Research Center via grants NNA09DB30A, NNX15AD20A, NNX16AF59G, as well as by a NASA SSERVI Cooperative Agreement to J.B. The theoretical work was partially funded by NASA ATP grant NNX15AK80G to S.F. and J.B. D.R. is supported by a NASA Postdoctoral Program Senior Fellowship at the NASA Ames Research Center, administered by the Universities Space Research Association under contract with NASA. We thank the following colleagues for their important contributions including Jill Bauman, Jay Bookbinder, Matthew D'Ortenzio, Robert Hanel, Butler Hine, Joseph Lazio, Stephanie Morse, Ken Galal, Laura Plice, Gary Rait, Jeremy Stober, and Eugene Tu. We also benefited from helpful input from James Aguirre, Mina Cappuccio, James Condon, Steven Ellingson, Bryna Hazelton, Danny Jacobs, Mark LaPole, Adam Lidz, John Mather, Harry Partridge, Aaron Parsons, Jonathan Pober, Greg Schmidt, Greg Taylor, and Sandy Weinreb.