Foreground and Sensitivity Analysis for Broadband (2D) 21 cm–Lyα and 21 cm–Hα Correlation Experiments Probing the Epoch of Reionization

We define the 3D power spectrum $P({\boldsymbol{k}})$ of the image cube $I({\boldsymbol{x}})$ as

$$\begin{eqnarray}&&P({\boldsymbol{k}})=\displaystyle \frac{\langle | \tilde{I}({\boldsymbol{k}}){| }^{2}\rangle }{V},\end{eqnarray} \tag{ 1 }$$

where $\tilde{I}({\boldsymbol{k}})$ is given by

$$\begin{eqnarray}&&\tilde{I}({\boldsymbol{k}})={dV}\displaystyle \sum _{{\boldsymbol{x}}}I({\boldsymbol{x}}){e}^{i{\boldsymbol{k}}\cdot {\boldsymbol{x}}},\end{eqnarray} \tag{ 2 }$$

Here V is the survey volume, and dV is the voxel size. Note that $P({\boldsymbol{k}})$ has units of ${[I]}^{2}\cdot {{\rm{Mpc}}}^{3}$ and that we often plot instead the more intuitive quantity ${\rm{\Delta }}(k)={[{k}^{3}P(k)/2{\pi }^{2}]}^{1/2}$, where P(k) is the average of $P({\boldsymbol{k}})$ within the 1D power spectrum bin k.

Similarly, over narrow fields of view, the angular power spectrum $C({\boldsymbol{\ell }})$ of a 2D (e.g., broadband) image $I({\boldsymbol{\theta }})$ can be shown to be approximately

$$\begin{eqnarray}&&C({\boldsymbol{\ell }})=\displaystyle \frac{\langle | \tilde{I}({\boldsymbol{\ell }}){| }^{2}\rangle }{{\rm{\Omega }}},\end{eqnarray} \tag{ 3 }$$

where $\tilde{I}({\boldsymbol{\ell }})$ is given by

$$\begin{eqnarray}&&\tilde{I}({\boldsymbol{\ell }})=d{\rm{\Omega }}\displaystyle \sum _{{\boldsymbol{\theta }}}I({\boldsymbol{\theta }}){e}^{i{\boldsymbol{\ell }}\cdot {\boldsymbol{\theta }}},\end{eqnarray} \tag{ 4 }$$

where Ω is the survey solid angle and $d{\rm{\Omega }}$ is the pixel size. Thus, over a narrow field of view, we need only evaluate a Fourier transform to estimate the angular power spectrum. Writing this out in detail gives⁵

$$\begin{eqnarray}&&C({\boldsymbol{\ell }}(a,b))={\left|\displaystyle \sum _{m,n}I(m,n)\exp \left(-\displaystyle \frac{2\pi i}{N}({am}+{bn})\right)\right|}^{2}\displaystyle \frac{d{\theta }^{2}}{{N}^{2}},\end{eqnarray} \tag{ 5 }$$

where $d\theta =d{\theta }_{x}=d{\theta }_{y}$ is the pixel size, $N\equiv {N}_{x}={N}_{y}$ is the number of pixels on a side of a square image, and $-N/2\leqslant a,b\leqslant N/2$ are integers. Further, ${{\ell }}^{2}={{\ell }}_{x}^{2}+{{\ell }}_{y}^{2}$, and

$$\begin{eqnarray}&&{{\ell }}_{x}=2\pi a/{Nd}\theta ,\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{{\ell }}_{y}=2\pi b/{Nd}\theta .\end{eqnarray} \tag{ 7 }$$

Note that $C({\boldsymbol{\ell }})$ has the units of ${[I]}^{2}\cdot d{\theta }^{2}$, and we often work with ${\rm{\Delta }}({\ell })={[{\ell }({\ell }+1)C({\ell })/2\pi ]}^{1/2}$, which has the same units⁶ as I. Here $C({\ell })$ is the average of $C({\boldsymbol{\ell }})$ within the 1D power spectrum bin ℓ.

The 3D 21 cm power spectrum is often cylindrically binned from 3D ${\boldsymbol{k}}$ space to 2D $({k}_{\perp },{k}_{\parallel })$ space, where ${k}_{\perp }^{2}\equiv {k}_{x}^{2}+{k}_{y}^{2}$ represents modes perpendicular to the line of sight and ${k}_{\parallel }={k}_{z}$ represents modes along the line of sight. We show in Appendix B that this cylindrically binned power spectrum is related to the angular power spectrum of a broadband image (over a narrow field of view) as

$$\begin{eqnarray}&&P({k}_{\perp },{k}_{\parallel }=0)={D}_{c}^{2}{\rm{\Delta }}{D}_{c}C({\ell }({k}_{\perp })).\end{eqnarray} \tag{ 8 }$$

Here ${\ell }={D}_{c}{k}_{\perp }$, where D_c is the comoving line-of-sight distance to the center of the cube and ${\rm{\Delta }}{D}_{c}$ is the comoving depth of the cube.

The 3D and 2D cross spectra are defined, extending Equations (1) and (3) to the cross spectrum, as

$$\begin{eqnarray}&&{P}_{12}({\boldsymbol{k}})=\displaystyle \frac{\langle {\tilde{I}}_{1}^{* }({\boldsymbol{k}}){\tilde{I}}_{2}({\boldsymbol{k}})\rangle }{V},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{C}_{12}({\boldsymbol{\ell }})=\displaystyle \frac{\langle {\tilde{I}}_{1}^{* }({\boldsymbol{\ell }}){\tilde{I}}_{2}({\boldsymbol{\ell }})\rangle }{{\rm{\Omega }}},\end{eqnarray} \tag{ 10 }$$

where 1 and 2 denote the 21 cm and the IR fields, respectively. The cross spectrum is a quantity that ranges between $\pm {({C}_{1}({\boldsymbol{\ell }}){C}_{2}({\boldsymbol{\ell }}))}^{1/2}$ in the 2D case, depending on how correlated, uncorrelated, or anticorrelated the two fields are. It is thus often renormalized as

$$\begin{eqnarray}&&{c}_{12}({\boldsymbol{\ell }})\equiv \displaystyle \frac{{C}_{12}({\boldsymbol{\ell }})}{\sqrt{{C}_{1}({\boldsymbol{\ell }}){C}_{2}({\boldsymbol{\ell }})}},\end{eqnarray} \tag{ 11 }$$

where c is known as the coherence and is insensitive to a simple rescaling of either field. However, uncorrelated foreground residuals in either field will substantially bias the coherence toward zero (Furlanetto Lidz 2007; Lidz et al. 2009), whereas they merely contribute a zero mean noise to the cross spectrum. Slight foreground correlations, of course, will bias the cross spectrum as well, as we explore later.

The MWA is a low-frequency radio interferometer in Western Australia consisting of 128 phased array tiles, each with $\sim 30^\circ \times (150\text{MHz}/f)$ beams (FWHM) and steerable in few-degree increments with a delay line beamformer. We use low-frequency observations of a quiet field centered at (R.A., decl.) = (0°, $-27^\circ$) J2000 recorded over 30.72 MHz bandwidth centered at 186 MHz, corresponding to $z=6.0\mbox{--}7.3$ for rest-frame 21 cm.

We use MWA image products produced by Beardsley et al. (2016). The MWA observations are recorded as 2-minute "snapshots" that are flagged for radio frequency interference using COTTER (Offringa et al. 2015) and then calibrated and imaged using Fast Holographic Deconvolution (FHD).⁷ Model visibilities are simulated from a foreground model of diffuse and point-source (Carroll et al. 2016) emission in the field and used for both calibration and foreground subtraction. For each snapshot, FHD produces naturally weighted data and model image cubes, as well as primary and synthesized beam cubes. FHD outputs these "cubes" in HEALPix format per frequency. Note that this processing is performed in parallel on "odd" and "even" data cubes whose data are interleaved at a 2 s cadence for the purpose of avoiding a thermal noise bias in autospectrum analyses. In cross-spectrum analyses, the noise between the radio and IR images is independent, so in principle we should average the odd and even cubes together to achieve the lowest noise. However, for simplicity, we use only the even cube in this work given that thermal noise is much smaller than residual foregrounds.

Following Dillon et al. (2015), we rotate these HEALPix maps so that the MWA field center lies at the north pole, and then we project the pixels onto the xy plane to obtain naturally weighted image space cubes of the raw data (${I}_{{\rm{nat}}}$), the model data (${I}_{{\rm{nat,mod}}}$), the synthesized beam (I_w) (i.e., the Fourier transform of the uv weights), and the primary beam, all in orthographic projection. We flag the upper and lower 80 kHz channels in each of 24 coarse channels across the band to mitigate aliasing, average each cube over frequency to make broadband images, and then apply uniform weighting using

$$\begin{eqnarray}&&{I}_{{\rm{uni}}}({\boldsymbol{\theta }})=\displaystyle \frac{{10}^{-26}{\lambda }^{2}}{{k}_{B}}\displaystyle \sum _{{\boldsymbol{u}}}\displaystyle \frac{{\tilde{I}}_{{\rm{nat}}}({\boldsymbol{u}})}{{\tilde{I}}_{w}({\boldsymbol{u}})}{e}^{-2\pi i{\boldsymbol{\theta }}\cdot {\boldsymbol{u}}}{d}^{2}u,\end{eqnarray} \tag{ 12 }$$

where ${\tilde{I}}_{i}({\boldsymbol{u}})={\sum }_{{\boldsymbol{\theta }}}{I}_{i}({\boldsymbol{\theta }}){e}^{2\pi i{\boldsymbol{\theta }}\cdot {\boldsymbol{u}}}d{\rm{\Omega }}$ for $i={\rm{nat}},w$. The units of ${I}_{{\rm{nat}}}$ and I_w are Jy sr⁻¹ and sr⁻¹, respectively, as seen from their approximate definitions of ${I}_{{\rm{nat}}}({\boldsymbol{\theta }})\approx {\sum }_{j}V({{\boldsymbol{u}}}_{j}){e}^{-2\pi i{{\boldsymbol{u}}}_{j}\cdot {\boldsymbol{\theta }}}{{du}}^{2}$ and ${I}_{w}({\boldsymbol{\theta }})\approx {\sum }_{j}{e}^{-2\pi i{{\boldsymbol{u}}}_{j}\cdot {\boldsymbol{\theta }}}{{du}}^{2}$, where the sums are over all measured visibilities $V({{\boldsymbol{u}}}_{j})$ in units of Jy. These definitions are only approximate because FHD performs corrections to account for wide-field effects.

ATLAS is a system of multiple 0.5 m f/2 wide-field telescopes (Tonry 2011) designed for near-Earth asteroid detection and tracking. Two telescopes are currently in operation located on the islands of Maui and Hawaii. Each telescope has a single $\mathrm{10,560}\times \mathrm{10,560}$ STA1600 CCD sensor, with a pixel scale of $1\buildrel{\prime\prime}\over{.} 86$, giving a field of view of $5\buildrel{\circ}\over{.} 5$ on a side. We observe in the Johnson I band (810 nm center with 150 nm FWHM) with the Maui telescope, corresponding to $z=5.1\mbox{--}6.3$ for Lyα. While this redshift range does not exactly match that of our radio observations, it overlaps sufficiently for our purpose of characterizing the noise and foregrounds in 21 cm–Lyα cross-correlation experiments. The automatic processing pipeline provides single sky-flattened (using a median of the nightly science exposures) images registered to nominal Two Micron All Sky Survey astrometry.

We perform two separate observing campaigns, which we illustrate in Figure 1. We first perform a wide survey to best characterize bright foregrounds. We raster-scan a roughly $20^\circ \times 20^\circ$ grid with 5° spacing over the MWA field (dashed black circle), integrating for 2.5 minutes at each pointing (blue filled squares). The observations were conducted between 2016 September 07 22:00 and 2016 September 08 00:50 Hawaii–Aleutian Standard Time, when the Moon was 36% illuminated. We use swarp⁸ (Bertin et al. 2002) to stack all these frames over a 20° orthographic field centered on (R.A., decl.) = $(0^\circ ,-30^\circ )$ (large blue square) with $1\buildrel{\prime\prime}\over{.} 86$ resolution, using the default background subtraction settings to mitigate temporal and spatial background variation.

Zoom In Zoom Out Reset image size

Our second campaign is a slightly deeper survey designed to better mitigate airglow fluctuations and CCD systematics for the purpose of studying faint foregrounds below the detection limit. We select four 5 deg fields positioned around the MWA beam peak for best cross-correlation precision: (R.A., decl.) = $(-2\buildrel{\circ}\over{.} 5,-24\buildrel{\circ}\over{.} 5)$, $(2\buildrel{\circ}\over{.} 5,-24\buildrel{\circ}\over{.} 5)$, $(-2\buildrel{\circ}\over{.} 5,-29\buildrel{\circ}\over{.} 5)$, $(2\buildrel{\circ}\over{.} 5,-29\buildrel{\circ}\over{.} 5)$ (J2000). We raster-scan a 3 × 3 grid of 30 s observations within each field (red filled circles) intended to mitigate slight amplifier nonuniformities across the CCD array. The observations were conducted on 2016 November 02 between 21:47 and 23:11 Hawaii–Aleutian Standard Time, when the Moon was 5% illuminated.

We stack the frames in each of the four deep fields using swarp over only the central $4^\circ \times 4^\circ$ region, over which all nine 30 s frames overlap (red squares). Otherwise, slight background discontinuities would be introduced by the different temporal coverage of different regions of the stack. In this stacking we disable background subtraction for the purpose of studying the effects of airglow-induced diffuse backgrounds.

To characterize the bright sources relevant to broadband 21 cm–Lyα and 21 cm–Hα intensity mapping correlation measurements, we calculate the correlations between catalogs at 185 MHz, 850 nm, and 4.5 μm as a function of mask depth. These bands correspond roughly to 21 cm, Lyα, and Hα, respectively, from $z\sim 6\mbox{--}7$.

We use the 185 MHz catalog reduced from deep observations of the MWA field depicted in Figure 1 by Carroll et al. (2016). Figure 2 (left panel) shows a histogram of source fluxes in this field. The survey depth varies somewhat over the MWA field owing to the varying primary beam, resulting in a catalog that is 50% complete down to 150 mJy and 95% complete down to 250 mJy. These completeness levels are shallower than the 70 mJy completeness quoted by Carroll et al. (2016) owing to our large rectangular field. The MWA's intrinsic astrometry is at the $2^{\prime} \mbox{--}3^{\prime}$ level, though Carroll et al. (2016) cross-match with higher-frequency catalogs to achieve order $10^{\prime\prime}$ astrometry.

Zoom In Zoom Out Reset image size

We use the W2 band of ALLWISE (Wright et al. 2010; Cutri et al. 2013) as our 4.5 μm catalog. We download the list of sources within the MWA field using the All Sky Search on the NASA/IPAC Infrared Science Archive⁹ and plot the histogram of source fluxes in Figure 2 (middle panel). This ALLWISE band is specified to be 95% complete down to 88 μJy (15.7 AB mag), though it has slight sky coverage nonuniformities due to satellite coverage.

Lastly, we run SExtractor¹⁰ (Bertin & Arnouts 1996) on our wide 20° ATLAS composite image to generate an 850 nm catalog. We allow local background bias and noise estimation, set pixel saturation at 20,000 counts to avoid artifacts, and use the AUTO aperture profile. We extract sources down to $3\sigma$ above the background, in order to achieve the most complete point-source mask. Given that our ATLAS observations have been calibrated and imaged through a preliminary pipeline, we cross-match these sources with sources closer than $1^{\prime\prime}$ in the AAVSO¹¹ Photometric All Sky Survey¹² (Henden et al. 2016), which is complete to ∼3 mJy. We find matches for $\sim 20 \%$ of ATLAS detections, not unreasonable given the deeper flux limit of ATLAS compared to APASS. The top panel of Figure 3 shows a 2D histogram of APASS versus ATLAS magnitude as a function of ATLAS magnitude. We fit a Gaussian to the relative magnitude for sources brighter than 13 mag and find that our roughly calibrated ATLAS sources are too bright by 0.279 ± 0.003 mag. Applying this correction, we plot a histogram of ATLAS source fluxes in the right panel of Figure 2, finding that our survey is complete to roughly 0.5 mJy, a factor of 6 deeper than the APASS survey.

Zoom In Zoom Out Reset image size

Having prepared catalogs of point-source foregrounds in our three bands, we proceed to study how they manifest in intensity mapping correlation experiments. Traditionally, radio/infrared correlations have been studied by cross-matching high-frequency radio detections with infrared sources coincident within a few arcseconds and then plotting radio versus infrared luminosity. Such studies have revealed the well-known radio–far-infrared correlation thought to be due to massive star formation (e.g., de Jong et al. 1985; Helou et al. 1985; Xu et al. 1994; Yun et al. 2001; Willott et al. 2003; Mauch & Sadler 2007). Massive stars blow out ionized bubbles, generating radio free–free emission correlated with the ionizing flux. Some fraction of these ionizing photons are absorbed by dust clouds and reprocessed into far-infrared emission (Xu et al. 1994). At radio frequencies lower than ∼10 GHz, synchrotron dominates over free–free emission, and the correlation is thought to arise from the acceleration of cosmic-ray electrons in these stars' supernovae.

Our approach is different. For all the advantages of broadband intensity mapping, sources cannot be localized to specific redshifts, meaning that it is foreground fluxes, not luminosities, whose correlations are of interest. Of course, compact foregrounds may be masked or subtracted to some residual level, but any correlation of these residual foreground fluxes could bury the EOR correlation. We begin in this section by analyzing foreground fluxes as a function of masking depth, and in the next section we turn to the foregrounds in residual images below the detection limit of these catalogs.

We begin by gridding all three catalog fluxes in Jy to the $20^\circ \times 20^\circ$ grid centered at (R.A., decl.) = $(0,-30^\circ )$ depicted in Figure 1 at $5^{\prime}$ resolution and calculating the zero-lag correlations between the images as

$$\begin{eqnarray}&&c=\displaystyle \frac{\langle {I}_{{\rm{rad}}}{I}_{{\rm{IR}}}\rangle -\langle {I}_{{\rm{rad}}}\rangle \langle {I}_{{\rm{IR}}}\rangle }{\sqrt{(\langle {I}_{{\rm{rad}}}^{2}\rangle -\langle {I}_{{\rm{rad}}}{\rangle }^{2})(\langle {I}_{{\rm{IR}}}^{2}\rangle -\langle {I}_{{\rm{IR}}}{\rangle }^{2})}},\end{eqnarray} \tag{ 13 }$$

where $\langle \rangle$ denotes an average over pixels, and the uncertainty due to sample variance is ${\rm{\Delta }}c\approx {N}_{{\rm{pix}}}^{-1/2}$, where ${N}_{{\rm{pix}}}$ is the total number of pixels in the image. Between MWA and WISE catalogs we find $c=-0.003\pm 0.005$, and between MWA and ATLAS catalogs we find $c=0.001\pm 0.005$. Both are consistent with zero, as expected, as the brightest sources in both infrared catalogs are likely stars, whose radio emission is expected to be negligible. As a first experiment, we recalculate these correlations after excluding the brightest 10% of sources in all three catalogs, effectively masking down to ${10}^{-3.75}$ Jy at 4.5 μm and 10⁻² Jy at 850 nm, and find ${c}_{{\rm{MWA}}-{WISE}}=0.031\pm 0.005$ and ${c}_{{\rm{MWA-ATLAS}}}=0.0086\,\pm 0.005$. The former is a $6\sigma$ detection and merits some investigation. How does this apparent correlation depend on the flux cut? What is it due to? And what does it mean for broadband correlation experiments? Further, does the MWA–ATLAS correlation remain consistent with zero at stricter flux cuts?

To begin to answer these questions, we plot in Figure 4 the 185 MHz–4.5 μm correlation (top left) and 185 MHz–850 nm correlation (top right) as a function of the masking depth (i.e., the maximum flux of remaining sources). We plot the S/Ns of these correlation measurements, taking the noise to be ${N}_{{\rm{pix}}}^{-1/2}$ as described above, in the next row. Note that adjacent cells in the correlation matrix plots are somewhat correlated, so a consistent positive sign is not in and of itself evidence of significance. We assess significance by comparing each correlation measurement individually with the expected noise (the S/N), as well as by checking that the correlation vanishes when the 185 MHz image is flipped (bottom two rows).

Zoom In Zoom Out Reset image size

The 185 MHz and 4.5 μm catalogs exhibit a positive correlation peaking at 0.0332 ± 0.005 after masking infrared sources down to 10⁻⁴ Jy (18.9 mag) and radio sources down to 1 Jy, and remain significant down to the completeness limits of these catalogs. There is no significant correlation detection after flipping the 185 MHz image, indicating that this detection is not an artifact of the analysis or of primary beam or vignetting effects. The 185 MHz and 850 nm catalogs exhibit a marginal $3\sigma$ correlation after masking infrared sources down to 10⁻³ Jy (16.4 mag) and radio sources down to 0.3 Jy, though it does not appear significant in comparison to the level of correlation noise in the flipped image.

To understand these findings, we begin by investigating which 4.5 μm sources are responsible for this correlation. We select the subset of sources detected in the WISE 3.4 μm, 4.5 μm, and 12 μm bands and plot them (Figure 5, left panel) in the ${W}_{23}\,\equiv$ [4.6 μm]–[12 μm] versus ${W}_{12}\,\equiv$ [3.4 μm]–[4.6 μm] color–color space used by Wright et al. (2010) to illustrate the separation between different types of sources. Nikutta et al. (2014) study more quantitatively how sources separate in this space, finding that stars are isolated in the region ${W}_{12}=-0.04\pm 0.03$, ${W}_{23}=0.05\pm 0.04$ (1σ). In the right panel, we plot the faintest 90% of sources (fainter than 18.25 mag at 4.6 μm) in the same color–color space and observe that this cut effectively cleanly excludes nearly all the stars. This explains the detection of a 185 MHz–4.5 μm correlation only after masking the brightest 10% of sources.

Zoom In Zoom Out Reset image size

To further probe which mid-infrared sources are responsible for this correlation, we make a rough cut to separate quasars and active galactic nuclei (AGNs; ${W}_{12}\gt 0.6$) from star-forming (SF) galaxies (${W}_{12}\lt 0.6$; Mingo et al. 2016). In Figure 6 we plot the power spectrum of 185 MHz sources (left panel), 4.5 μm AGNs (middle panel, blue points), and 4.5 μm SF galaxies (middle panel, red points). In the right panel we plot the coherence (i.e., normalized cross spectrum) of the 185 MHz catalog with AGNs (blue) and with SF galaxies (red). The AGN cut exhibits no significant correlation with the 185 MHz sources, while the SF galaxy cut exhibits a significant correlation rising from a few percent at ${\ell }\sim 7000$ to 16% at ${\ell }\sim 300$.

Zoom In Zoom Out Reset image size

The fall of the correlation toward high ℓ  is likely due to the MWA's $3^{\prime}$ resolution at 185 MHz, corresponding to a maximum ℓ  of roughly 4000. Both the falling 4.5 μm catalog power spectrum and the relatively flat 185 MHz power spectrum are functions of the detailed properties of these surveys. Tegmark et al. (2002) and Dodelson et al. (2002) show that the galaxy angular power spectrum $C({\ell })$ is approximately equal to the 3D matter power spectrum $P(k({\ell }))$ convolved with a window function that depends on the redshift coverage and flux limit of the sample. The matter power spectrum is known to rise as k¹ for $k\lesssim 0.02\,{\rm{h}}$ Mpc⁻¹ before falling as ${k}^{-3}$. Galaxy surveys typically probe the regime just after the turnover where the slope is transitioning from 0 to −3 (Tegmark & Zaldarriaga 2002). In order to maximize its sensitivity to low surface EOR 21 cm emission, the MWA was designed as a relatively compact array in comparison to higher-resolution radio interferometers such as the Very Large Array. This low resolution makes the MWA catalog severely flux limited (Carroll et al. 2016), which in turn effectively masks many galaxies that would otherwise be seen. This large masked volume translates into a wide Fourier convolution kernel, explaining why the MWA catalog power spectrum is so flat.

Lastly, we hypothesize that the absence of an observed 185 MHz–850 nm correlation is due to the larger fraction of stars in 850 nm images than in 4.5 μm images. To check this, we study the fraction of stars and galaxies in a similar Galactic field observed by the Sloan Digital Sky Survey (SDSS; Eisenstein et al. 2011) DR13 (SDSS Collaboration et al. 2016). SDSS catalog¹³ objects are classified as either stars or galaxies using ugriz photometry, though unfortunately SDSS does not reach the declination of our MWA field at (R.A., decl.) = ($0,-27^\circ$). We thus use a 5° field centered at (R.A., decl.) = (205°, 22°), which has the same Galactic longitude but is flipped to the other side of the Galactic plane, giving a similar line of sight through the Galaxy.

In Figure 7 we plot a histogram of the fluxes of the objects marked as stars (red) and galaxies (blue) in this field at the SDSS i band $(762\pm 65)$ nm. We observe that stars dominate the field above 10⁻⁴ Jy; this is a factor of a few below the survey depth of our wide ATLAS 850 nm image, and thus none of the flux cuts explored above were deep enough to reveal the extragalactic sources. This is consistent with the fact that the 185 MHz–4.6 μm correlation only appeared after the stars were removed, a procedure easier in the mid-infrared than the near-infrared.

Zoom In Zoom Out Reset image size

Let us now consider why the 4.5 μm SF sample is 5%–15% correlated with the 185 MHz catalog, while the AGN sample is not. Of course, some slight correlation is expected simply because brighter AGNs typically reside in more massive galaxies, which are typically brighter in stars (see, e.g., Figure 1 in Seymour et al. 2007, or Figure 4 in Willott et al. 2003), but Mauch & Sadler (2007) find no strong correlation between radio and near-infrared luminosities. As discussed above, though, broadband correlation intensity mapping experiments are affected not only by luminosity correlations but by flux correlations as well. We show in this section that fluxes in two different bands may appear correlated owing to geometric effects even when their intrinsic luminosities are completely independent of each other. By geometric effects we refer to to the fact that more distant objects are generally weaker in all bands than nearer objects.

We first make a few approximations to build intuition and then simulate the effect as a function of source masking depth. Consider a sky survey with fixed field of view of a set of objects with uncorrelated infrared and radio emission. By uncorrelated we mean that the infrared and radio luminosities are independent random variables determined by the infrared and radio luminosity functions, respectively. Assume that the objects are uniformly distributed in space out to $z\sim 0.5$, and work in Cartesian space for simplicity. We are interested in the effective correlation between radio and infrared fluxes in the same sky pixels, but let us approximate this by calculating the correlation between source fluxes in the two bands. Starting from Equation (13), we have

$$\begin{eqnarray}&&c=\displaystyle \frac{\langle {F}_{{\rm{rad}}}{F}_{{\rm{IR}}}\rangle -\langle {F}_{{\rm{rad}}}\rangle \langle {F}_{{\rm{IR}}}\rangle }{\sqrt{(\langle {F}_{{\rm{rad}}}^{2}\rangle -\langle {F}_{{\rm{rad}}}{\rangle }^{2})(\langle {F}_{{\rm{IR}}}^{2}\rangle -\langle {F}_{{\rm{IR}}}{\rangle }^{2})}},\end{eqnarray} \tag{ 14 }$$

where $\langle \rangle$ denote an average over sources in the catalog. Then making the approximation that the radio and infrared luminosity functions are independent of line-of-sight distance, we can rewrite this equation in terms of moments of these luminosity functions and the distribution of object distances using ${F}_{i}={L}_{i}/4\pi {D}^{2}$ for $i\,=$ radio, IR

$$\begin{eqnarray}&&c=\displaystyle \frac{\beta -1}{\sqrt{(\beta {\alpha }_{{\rm{rad}}}-1)(\beta {\alpha }_{{\rm{IR}}}-1)}}\approx \displaystyle \frac{1}{\sqrt{{\alpha }_{{\rm{rad}}}{\alpha }_{{\rm{IR}}}}},\end{eqnarray} \tag{ 15 }$$

where $\beta \equiv \langle {D}^{-4}\rangle /\langle {D}^{-2}{\rangle }^{2}$ and ${\alpha }_{i}=\langle {L}_{i}^{2}\rangle /\langle {L}_{i}{\rangle }^{2}$, and again $\langle \rangle$ denote an average over sources in the catalog. The approximation in the last part of this equation results from the fact that for typical parameters and luminosity functions (see below) $\beta \gt \gt 1$ and $\alpha \gt 1$.

For a survey of a fixed angular field of view, uniform spatial distribution of objects, and Cartesian spacetime, the distribution of object distances is $\rho (D)={\rho }_{0}{D}^{2}$. Assuming that all the objects are between distances ${D}_{{\rm{\min }}}$ and ${D}_{{\rm{\max }}}$, the normalization constant is ${\rho }_{0}=3/({D}_{{\rm{\max }}}^{3}-{D}_{{\rm{\min }}}^{3})$, and

$$\begin{eqnarray}&&\beta =\displaystyle \frac{{\int }_{{D}_{{\rm{\min }}}}^{{D}_{{\rm{\max }}}}{D}^{-4}\rho (D){dD}}{{\left({\int }_{{D}_{{\rm{\min }}}}^{{D}_{{\rm{\max }}}}{D}^{-2}\rho (D){dD}\right)}^{2}}\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&=\,\displaystyle \frac{{D}_{{\rm{\max }}}^{2}+{D}_{{\rm{\max }}}{D}_{{\rm{\min }}}+{D}_{{\rm{\min }}}^{2}}{3{D}_{{\rm{\max }}}{D}_{{\rm{\min }}}}.\end{eqnarray} \tag{ 17 }$$

We observe that the radio–infrared flux correlation c for some type of objects is a function of their radio and infrared luminosity functions. In fact, we can see immediately that if the luminosity distributions are wide, their α values are large, and c is small. Conversely, if the luminosity functions are narrow, then the distance to the sources plays a more significant role in determining their fluxes, and so c is larger.

To quantify whether this effect can explain our measured radio–infrared correlation in SF galaxies and the lack of one in AGNs, we use AGN and SF luminosity functions at 1.4 GHz from Mauch & Sadler (2007) (Figure 8, right panel) and at 8 μm from Fu et al. (2010) (left panel). The former describe galaxies at $z\lt 0.3$, while the latter describe galaxies at $z\sim 0.6$. In principle, we should use luminosity functions at our actual radio and infrared bands, 185 MHz and 4.5 μm, and of course this analysis could be extended using proper redshift-dependent luminosity functions, though we find that our simplified analysis suffices to explain our earlier correlation measurements. We leave a more detailed study for future work. Indeed, Prescott et al. (2016) find that the AGN and SF galaxy radio luminosity functions at lower frequencies, specifically at 325 MHz, closely follow those at 1.4 GHz up to an overall scaling that cancels out of our correlation coefficient. We use approximately the same range of luminosities used by Mauch & Sadler (2007) and Fu et al. (2010) and adjust the minimum luminosities slightly to achieve the same number density of each type of object in both radio and infrared surveys. Above these limiting magnitudes (see Figure 8), the number density of AGNs is ∼0.0020 Mpc⁻³ and that of SF is ∼0.00011 Mpc⁻³. In the end we find that our results are only weakly sensitive to these luminosity minima, as their faint ends become less and less significant in real, flux-limited samples.

Zoom In Zoom Out Reset image size

We pick fiducial survey parameters of ${d}_{{\rm{\min }}}=20$ Mpc and ${d}_{{\rm{\max }}}=3000\,\mathrm{Mpc}$ (${z}_{{\rm{\max }}}=0.75$), giving $\beta \approx 50$. Using the above luminosity functions, we find ${\alpha }_{{\rm{SF,IR}}}=1.474$, ${\alpha }_{{\rm{SF,rad}}}\,=14.56$, ${\alpha }_{{\rm{AGN,IR}}}=22.97$, and ${\alpha }_{{\rm{AGN,rad}}}=257.5$. These values agree with qualitative observation that the AGN luminosity function is wider than the SF luminosity function in both radio and infrared bands (Figure 8). These values give a predicted radio–infrared correlation of 0.21 for SF galaxies and 0.01 for AGNs, agreeing with our finding of a significant radio–infrared correlation for SF and near-zero correlation for AGNs. The exact values deviate from our measurements for a number of reasons. The MWA and WISE catalogs are not matched in depth or redshift coverage and thus do not survey an exactly overlapping set of radio and infrared sources. Further, these calculations assume a volume-limited survey, in contrast to our flux-limited radio and IR surveys. Additionally, real-world luminosity functions can exhibit redshift evolution. Of course, our measurement in Section 4.2 did not even split up the MWA catalog into separate AGN and SF subsets, as such detailed characterization of low-frequency radio foregrounds remains an active area of research. Lastly, the limited MWA resolution pushes the observed correlation with infrared images to zero at high ℓ, suppressing the overall correlation computed in image space. Future work will be needed to quantify these effects in greater detail and assess their significance in EOR cross-spectrum measurements.

Can this unwanted radio–infrared foreground correlation be mitigated by masking the brightest sources? Using the luminosity functions presented above, we simulate radio and infrared surveys for each of AGNs and SF galaxies. We begin by generating the mock radio catalogs of AGNs and SF galaxies, choosing a Poisson random number of each in each of 400 logarithmic luminosity bins. Using logarithmic bins best samples the large dynamic range of the luminosity functions. We distribute the objects uniformly over a volume ${D}_{{\rm{\max }}}={{cz}}_{{\rm{\max }}}/{H}_{0}=3212$ Mpc deep and ${\theta }_{{\rm{FOV}}}{D}_{{\rm{\max }}}=1121$ Mpc wide and then pick a random infrared luminosity for each radio object from the appropriate infrared luminosity function. Finally, we plot in Figure 9, along the lines of Figure 4, the predicted 1.4 GHz–8 μm correlation of our mock AGN and SF catalogs after masking down to a maximum radio and infrared flux. As we saw above, without any flux cut we find a roughly 20% radio–infrared correlation for SF galaxies and negligible correlation for AGNs. As we mask fainter and fainter sources, the AGN correlation generally increases to the 5%–10% level, while the SF correlation first increases and then decreases after masking down to 10⁻⁴ Jy. With increasing mask depth, these correlations do not necessarily approach zero monotonically, and more detailed modeling of effective foreground flux correlations will be necessary in real-world intensity mapping correlation experiments probing the EOR. In the next section we move beyond the bright sources and study the magnitudes and correlation properties of the residual radio and infrared foregrounds in our MWA and ATLAS observations.

Zoom In Zoom Out Reset image size

We begin by quantifying the 185 MHz foreground residuals in angular power spectrum measurements. In broadband (i.e., multifrequency synthesis) images, thermal noise quickly integrates below foreground residuals. Reaching the cosmological signal should therefore be a matter of foreground mitigation and not the long time averages needed to measure the 3D power spectrum (e.g., Beardsley et al. 2013; Pober et al. 2014). We check this hypothesis, asking how much observation time is required to achieve the best foreground subtraction.

Working with the MWA image products presented in Section 3.1, we compute the angular power spectrum as

$$\begin{eqnarray}&&C({\ell })=\displaystyle \frac{{\displaystyle \sum }_{{\boldsymbol{\ell }}}| {\tilde{I}}_{{\rm{uni}}}({\boldsymbol{\ell }}){\tilde{I}}_{w}({\boldsymbol{\ell }}){| }^{2}}{{\displaystyle \sum }_{{\boldsymbol{\ell }}}| {\tilde{I}}_{w}({\boldsymbol{\ell }}){| }^{2}}\displaystyle \frac{d{\theta }^{2}}{{N}^{2}}\end{eqnarray} \tag{ 18 }$$

summing over all the ${\boldsymbol{\ell }}$ values in each ℓ bin, where N is the number of pixels of size $d\theta$ on each side of the square image. Note that ${\ell }=2\pi u$, where u is the Fourier dual to angle from the field center assuming the small-angle approximation. We estimate the thermal noise power spectrum by computing the power spectrum of the difference between the interleaved odd and even cubes discussed earlier, which contains only thermal noise.

We plot in Figure 10 the power spectra of 185 MHz broadband images from 3 hr data selections, spread over 1 (red solid), 5 (green solid), and 19 (blue solid) nights. We also make a broadband image of the same 32 hr data set used by Beardsley et al. (2016) (black dashed). We plot the power spectra of the raw (pre-foreground subtraction), residual (post-foreground subtraction), and noise (difference between successive integrations) images out to ${\ell }=2600$, corresponding to a maximum baseline length of ∼700 m. Beyond this, the uv coverage becomes sparse, introducing artifacts in the application of gridding and uniform weighting, though a more sophisticated analysis could likely use these longer baselines. We cross-check our imaging and power spectrum analysis by comparing the ${k}_{\parallel }=0$ bin of the 3D power spectrum of Beardsley et al. (2016) converted to an angular power spectrum using Equation (8) (magenta dashed).

Zoom In Zoom Out Reset image size

We find that the raw power spectra of the different 3 hr data sets agree with each other and with that of the 32 hr data set, as expected given that foregrounds overwhelm thermal noise in a broadband image. Interestingly, the residual power spectrum decreases as the 3 hr are spread over more and more nights until it reaches the level of the deep 32 hr integration, a factor of ∼4 lower in power than in the single-night analysis. These findings could be explained by slight ionosphere-related errors that limit the accuracy of each night's calibration and thus of its foreground subtraction. Further work would be needed to understand this effect in detail, but for now our conclusion is that the 32 hr integration has the best foreground subtraction, and we use this deep cube in our later analyses. Note that, as expected, the thermal noise of the deep integration is 10 times lower in power than the 3 hr integrations. Because these are band-averaged images, even the 3 hr thermal noise is at least 100 times lower than its foreground residuals.

We proceed to generate foreground-masked 850 nm images of each of the four deep ATLAS integration fields shown in Figure 1. Each of these fields is a stack of nine 30 s exposures with $5^\circ$ field of view, dithered such that the overlap is a $\sim 4^\circ$ region. By confining ourselves to this overlap region, we avoid the background discontinuities that affect many nominally wide-field infrared image data sets whose mosaicking introduces significant background patchiness. Mitchell-Wynne et al. (2015) have demonstrated that complex fitting (Fixsen et al. 2000) can help reduce such patchiness in small mosaics (∼10 arcminutes), but further work is required to study whether such techniques can be applied to much larger fields.

Our approach is to mask each image at ATLAS's native $1\buildrel{\prime\prime}\over{.} 86$ resolution and then coarse grid down to $3\buildrel{\,\prime}\over{.} 5$, approximately the resolution of the MWA, taking each coarse pixel's value as the average of all unmasked fine pixels within it. If fewer than 10% of its fine pixels remain unmasked, we consider the whole coarse pixel masked to avoid introducing too much noise variation between different coarse pixels. For illustrative purposes, we proceed in the following four stages, which we illustrate in Figure 11 for the ATLAS field centered at (R.A., decl.) $=\,(2\buildrel{\circ}\over{.} 74,-24\buildrel{\circ}\over{.} 79$) (the top right red box in Figure 1). Each row shows the result of an additional masking stage, as outlined below. The left column shows a typical $9^{\prime}$ field to illustrate the masking up close, the middle column shows the resulting coarse binned image with $3\buildrel{\,\prime}\over{.} 5$ resolution, and the right column shows the fast Fourier transform (FFT) of the center image (plotted as $\mathrm{log}[{\rm{\Delta }}({\ell })/({\text{kJy sr}}^{-1})]$) in order to identify detector systematics.

• 1.  
  Mask saturated regions (Figure 11, row 1). Nearly saturated pixels are associated with nearby bright stars that would dominate fluctuation measurements; thus, we mask $4^{\prime}$ around all pixels within 30% of saturation (white regions in left and middle columns). This wide mask removes the broad wings of the point-spread function (PSF) revealed by these extremely bright sources. Roughly 92% of fine pixels remain unmasked after this step, and 96% of coarse pixels remain.
• 2.  
  Mask sources to 5σ (Figure 11, row 2). We mask circular regions with radius equal to five times the profile rms along the minor axis of each source as measured by SExtractor (see Section 4.1). The reason we use the minor-axis rms is that vertical charge leakage in the CCD array results in unrealistically large major-axis rms measurements for bright sources. After this stage ∼87% of fine pixels remain unmasked, and the fraction of coarse pixels remaining is unaffected.
• 3.  
  Mask sources to 12σ and mask other emission above 70 kJy sr⁻¹ over the background (Figure 11, row 3). We use a larger masking radius to remove the PSF wings around the 5σ source masks. We also determine that the vertical CCD charge leakage can be isolated by looking for emission brighter than 70 kJy sr⁻¹ over the background, and that it can be flagged without cutting into the shot noise between sources. A total of 58% of fine pixels remain unmasked after this step, and again the fraction of coarse pixels remaining is unaffected.
• 4.  
  Mask horizontal and vertical Fourier modes (Figure 11, row 4). The previous two stages revealed detector artifacts within ${\rm{\Delta }}{\ell }\sim 100$ of ${{\ell }}_{x}=0$ and ${{\ell }}_{y}=0$. These compact Fourier systematics correspond to slight horizontal and vertical discontinuities in the center image due to imperfect gain matching between 16 different amplifiers that process different rectangular regions of the CCD array. We conservatively mask Fourier modes with $| {{\ell }}_{x}| \lt 200$ or $| {{\ell }}_{y}| \lt 200$ to eliminate this effect.

Zoom In Zoom Out Reset image size

Note that these 850 nm deep observations were recorded during near new moon conditions, and we find that the mean airglow in source-free regions is ∼3 × 10³ kJy sr⁻¹, of order 19 AB mag arcsec⁻². For comparison, Sullivan & Simcoe (2012) measure a 1020 nm continuum airglow brightness (i.e., after spectrally masking the OH lines) of $20\pm 0.5$ AB mag arcsec⁻² far away from the Moon.

We proceed to characterize the residual infrared fluctuations in power spectrum space, using the optimal quadratic estimator to properly account for the masking of coarse pixels. This estimator was introduced to astronomy by Tegmark (1997) to recover CMB power spectra from maps with arbitrary survey geometries and noise properties and has recently been revived for 3D power spectrum analysis of 21 cm data by Dillon et al. (2014, 2015), Liu & Tegmark (2011), Dillon et al. (2013), and Ali et al. (2015). We employ this estimator in a manner more similar to the original CMB case, using it to account for pixel masking in broad bandpower spectrum estimation. We briefly summarize the estimator here and refer to Dillon et al. (2014) for a more detailed description.

We label the normalized estimator of the power in bin α as p_α, related to the unnormalized estimator q_α as ${\boldsymbol{p}}={\boldsymbol{Mq}}$. Bold lowercase letters are vectors, while bold uppercase letters are matrices. The unnormalized estimator is given by

$$\begin{eqnarray}&&{q}_{\alpha }=\displaystyle \frac{1}{2}{({\boldsymbol{x}}-\langle {\boldsymbol{x}}\rangle )}^{t}{{\boldsymbol{C}}}^{-1}{{\boldsymbol{C}}}_{,\alpha }{{\boldsymbol{C}}}^{-1}({\boldsymbol{x}}-\langle {\boldsymbol{x}}\rangle ),\end{eqnarray} \tag{ 19 }$$

where ${\boldsymbol{x}}$ is a column vector containing all ${N}_{{\rm{pix}}}$ pixel measurements, ${\boldsymbol{C}}$ is the pixel–pixel covariance matrix, and ${{\boldsymbol{C}}}_{,\alpha }\equiv d{\boldsymbol{C}}/{{dp}}_{\alpha }$ is the derivative of the covariance with respect to the power in bin α. Note that ${{\boldsymbol{C}}}_{,\alpha }={{\boldsymbol{A}}}^{\dagger }{\boldsymbol{A}}$, where ${\boldsymbol{A}}$ is an ${N}_{\alpha }\times {N}_{{\rm{pix}}}$ with elements ${A}_{{ij}}=\exp (i{{\boldsymbol{\theta }}}_{j}\cdot {{\boldsymbol{\ell }}}_{i})$. Here ${{\boldsymbol{\ell }}}_{i}$ refers to the ith out of N_α ${\boldsymbol{\ell }}$ modes in bin α. Note that ^t denotes a transpose and † donates a conjugate transpose.

The matrix of window functions (i.e., horizontal error bars) of the bandpowers p_α, defined such that ${{\boldsymbol{p}}}_{{\rm{estimated}}}={{\boldsymbol{Wp}}}_{{\rm{true}}}$, is given by ${\boldsymbol{W}}={\boldsymbol{MF}}$, where ${\boldsymbol{F}}$ is the Fisher matrix and ${\boldsymbol{M}}$ is an arbitrary invertible normalization function encoding the compromise between horizontal and vertical error bars. The covariance between the measured p_α values is given by ${\boldsymbol{\Sigma }}={{\boldsymbol{MFM}}}^{t}$. Dillon et al. (2014) argue that taking ${\boldsymbol{M}}\propto {{\boldsymbol{F}}}^{-1/2}$ is a good compromise between small horizontal error bars and small vertical error bars. For simplicity, we take ${\boldsymbol{M}}={{\boldsymbol{F}}}^{-1/2}$ and correct the normalization at the end by dividing each element of ${\boldsymbol{p}}$ by the peak of the appropriate row of ${\boldsymbol{W}}$. In this case, the bandpower variances are given by the reciprocals of the peaks of the window functions used to normalize the bandpowers. We find that using the sum instead of the peak significantly biases the recovered bandpowers downward when the power spectrum is nonflat, as in our case. Lastly, the elements of the Fisher matrix are given by

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{\alpha \beta }=\displaystyle \frac{1}{2}{\rm{tr}}({{\boldsymbol{C}}}^{-1}{{\boldsymbol{C}}}_{,\alpha }{{\boldsymbol{C}}}^{-1}{{\boldsymbol{C}}}_{,\beta }).\end{eqnarray} \tag{ 20 }$$

Now we turn to application of this formalism to power spectrum estimation from our masked IR images. Later we will adapt it to estimation of the 21 cm–IR cross spectrum. We take ${\boldsymbol{x}}$ to be a vector of all IR coarse ($3\buildrel{\,\prime}\over{.} 5$) pixel values, with masked pixels set to zero. After gridding the high-resolution images to reach this resolution, photon shot noise is negligible, and the image space covariance is the sum of the sample variance ${{\boldsymbol{C}}}_{{\rm{signal}}}$ and the masking covariance ${{\boldsymbol{C}}}_{{\rm{mask}}}$. ${{\boldsymbol{C}}}_{{\rm{mask}}}$ is a diagonal matrix with $\infty$ for masked pixels and 0 otherwise. In practice, we replace $\infty$ with a number 10⁷ times larger than the largest eigenvalue of ${{\boldsymbol{C}}}_{{\rm{signal}}}$, finding that the results are not sensitive to this parameter. The sample variance is easily obtained by writing it in Fourier space, ${{\boldsymbol{C}}}_{{\rm{ft}}}$, where it is a diagonal matrix with a guess of the true power spectrum on the diagonal, and then Fourier transforming it into image space with Fourier transform matrix ${ \mathcal F }$. Putting these together gives

$$\begin{eqnarray}&&{\boldsymbol{C}}={{ \mathcal F }}^{\dagger }{{\boldsymbol{C}}}_{{\rm{ft}}}{ \mathcal F }+{{\boldsymbol{C}}}_{{\rm{mask}}}.\end{eqnarray} \tag{ 21 }$$

Note that ${ \mathcal F }$ is an ${N}_{{\rm{pix}}}\times {N}_{{\rm{pix}}}$ matrix with elements ${{ \mathcal F }}_{{ij}}=\exp (-i{{\boldsymbol{\theta }}}_{i}\cdot {{\boldsymbol{\ell }}}_{j})$, where i runs over all pixels and j runs over all Fourier cells. Said differently, a guess of the power spectrum is necessary to optimally downweight the sample variance noise on the estimated power spectrum. If the accuracy of the guess were in question, one could always iterate by feeding the estimated power spectrum back into the quadratic estimator, though in practice we find that this is not necessary in this work.

Using this estimator, we calculate the power spectrum after each stage of masking and plot the mean spectrum over all four deep ATLAS fields in Figure 12 (left panel). Instead of predicting the bandpower errors from the input covariances, we conservatively bootstrap the error bars by computing the standard deviation of each bandpower over the four fields. The power spectrum of all 850 nm sources, masking only saturated regions, rises proportionally to ℓ (red dots), as expected for Poisson source counts when the power spectrum is plotted as ${\rm{\Delta }}=\sqrt{{{\ell }}^{2}{C}_{{\ell }}/2\pi }$. Masking sources out to $5\sigma$ removes two orders of magnitude in power (green dots), and masking out to $12\sigma$ and above 70 kJy sr⁻¹ over the background removes a factor of a few more in power (blue dots). Lastly, excluding modes with $| {{\ell }}_{x}| \lt 200$ or $| {{\ell }}_{y}| \lt 200$ gives our final 850 nm anisotropy spectrum (black dots). We note that more stringent masking does not significantly alter this final result. After the first masking stage, we use a flat power spectrum in $C({\ell })$ as our guess in the quadratic estimator formalism, and after the other three masking steps we use the 1.1 μm power spectrum from Zemcov et al. (2014).

Zoom In Zoom Out Reset image size

In Figure 12 (right panel) we compare our final residual power spectrum measurement with other measurements¹⁴ in the literature. Our 850 nm ATLAS spectrum (black dots) agrees very well with the 1.1 ± 0.25 μm CIBER spectrum (Zemcov et al. 2014) (blue dots), with much smaller error bars at the ℓ modes we can access (${\ell }\lesssim 4000$) due to ATLAS's larger field of view.

Zemcov et al. (2014) argue that their spectrum is limited by foregrounds at all ℓ: by diffuse Galactic light (DGL) (i.e., dust) at ${\ell }\lesssim 1000$, by intrahalo light (IHL; cyan dashed line) at ${\ell }\sim 1000$, and by galaxy number counts below the flux limit at larger ℓ. ATLAS's $2^{\prime\prime}$ resolution is only a factor of 3 finer than CIBER's $7^{\prime\prime}$ resolution, so perhaps it is not surprising that the level of galaxy number counts is not appreciably different. The ATLAS points at ${\ell }\gt 1000$ are in fact slightly above CIBER's, though the measurements were conducted over slightly different bands. Mitchell-Wynne et al. (2015) demonstrate a dramatic improvement in foreground reduction ℓ ≳ 10³ in $10^{\prime}$ HST mosaics with $0\buildrel{\prime\prime}\over{.} 1$ resolution. Their power spectrum (red points) is a factor of ∼10 in amplitude below the previous IHL model, and their new fitting suggests that Zemcov et al. (2014) were in fact limited almost exclusively by DGL (red solid line). For comparison, we plot their new IHL model (red dashed) and EOR model (gray area). These results show that finer angular resolution can help reduce the foregrounds in the ${\ell }\sim {10}^{2}\mbox{--}{10}^{4}$ modes that can be cross-correlated with 21 cm EOR intensity mapping experiments.

Before turning to 21 cm–Lyα cross-spectrum measurements with our 185 MHz and 850 nm images, we generate optimistic and pessimistic theoretical cross spectra for comparison. We simulate 21 cm and Lyα cubes using 21 cm power spectra from Pober et al. (2014), Lyα power spectra from Gong et al. (2014), and the coherence between the two fields from Heneka et al. (2016). Combining simulations from all these sources allows us to better estimate the modeling uncertainty. Future work is needed to more self-consistently model these fields and their correlation over a range of possible reionization scenarios and to infer astrophysical parameters from observed cross spectra.

Gong et al. (2014) model the Lyα cross spectrum and plot an uncertainty region over a range of likely values of escape fraction of ionizing photons, fraction of radiation emitted at Lyα, SF rate, and intergalactic medium (IGM) clumping factor (their Figure 1). We take the upper and lower edges of their uncertainty region as our optimistic and pessimistic power spectra, respectively.

Similarly, Pober et al. (2014) simulate 21 cm power spectra over a range of reionization scenarios using Mesinger et al. (2011), with various values of ionizing efficiency (ζ) (including the escape fraction), the minimum halo virial temperature (${T}_{{\rm{vir}}}$), and the ionizing photon mean free path in the IGM (${R}_{{\rm{mfp}}}$). Whether the signal at our redshift of interest (z = 7) is large or not depends mostly on whether reionization is already largely finished by that time or not. So we take as our pessimistic model the z = 8 power spectrum of the ($\zeta =31.5$, ${T}_{{\rm{vir}}}=1.5\times {10}^{4}$ K, ${R}_{{\rm{mfp}}}=30$ Mpc) scenario, whose reionization midpoint is z = 9.5. We take as our optimistic model a late reionization scenario with ${T}_{{\rm{vir}}}=3\times {10}^{5}$ K (other parameters unchanged), whose reionization midpoint is z = 5.5.

From these power spectra and coherence functions, we generate approximate 21 cm and Lyα cubes assuming Gaussian statistics. This is an approximation, given that the 21 cm field is expected to become less and less Gaussian as reionization proceeds (Wyithe & Morales 2007), and more work is needed to understand this effect. The simulated cubes have (1 Mpc)³ resolution over a (218 Mpc)³ volume at z = 7, corresponding to ${\rm{\Delta }}z=0.6$ and $0\buildrel{\,\prime}\over{.} 4$ angular resolution. Our 21 cm and Lyα cubes have units of mK and kJy sr⁻¹, respectively, and we average them in the line-of-sight direction to produce broadband images.

We plot the original 3D spherically averaged power spectra and coherence function as solid lines in Figure 13, and we plot the computed 2D power spectra from the line-of-sight averaged cubes, as well as their coherence, as dashed lines. The 3D and 2D power spectra are plotted as ${\rm{\Delta }}(k)$ and ${\rm{\Delta }}({\ell })$, respectively, as defined in Section 2.1. As expected, line-of-sight averaging tends to remove power on short spatial scales (compared to the line-of-sight depth of the cube), but it acts similarly on both cubes, so the coherence function is preserved. Note that the cross spectrum ${{\rm{\Delta }}}_{12}$ is defined in terms of the coherence as ${{\rm{\Delta }}}_{12}^{2}=c{[{{\rm{\Delta }}}_{1}^{2}{{\rm{\Delta }}}_{2}^{2}]}^{1/2}$ for both the 3D and 2D cases, and we use the same coherence function from Heneka et al. (2016) in the optimistic and pessimistic scenarios. These results show that 2D experiments can detect much of the 3D fluctuation power at $k\lt 0.2$ Mpc⁻¹, and motivate more complete and self-consistent simulations of the 21 cm and Lyα fields throughout the EOR over a range of optimistic and pessimistic scenarios.

Zoom In Zoom Out Reset image size

Having characterized the residual sky power spectrum at 185 MHz and 850 nm after applying the best foreground masking and subtraction permitted by our data sets, we now search for a correlation between these radio and foreground residuals. To account for the nonuniform uv sampling of the MWA image, as well as for the image space masking of the ATLAS image, we again use the optimal quadratic estimator. In Appendix C we show that with the approximation that the correlation between the two images is small, the proper extension of the optimal quadratic estimator from the autospectrum case to the cross-spectrum case is given by

$$\begin{eqnarray}&&{q}_{\alpha }\approx {({{\boldsymbol{x}}}_{21}-\langle {{\boldsymbol{x}}}_{21}\rangle )}^{t}{{\boldsymbol{C}}}_{21}^{-1}{{\boldsymbol{C}}}_{,\alpha }{{\boldsymbol{C}}}_{{\rm{IR}}}^{-1}({{\boldsymbol{x}}}_{{\rm{IR}}}-\langle {{\boldsymbol{x}}}_{{\rm{IR}}}\rangle ),\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{\alpha \beta }\approx {\rm{tr}}({{\boldsymbol{C}}}_{21}^{-1}{{\boldsymbol{C}}}_{,\alpha }{{\boldsymbol{C}}}_{{\rm{IR}}}^{-1}{{\boldsymbol{C}}}_{,\beta }),\end{eqnarray} \tag{ 23 }$$

where ${{\boldsymbol{C}}}_{,\alpha }$ is the same matrix used in Section 5.2, and ${{\boldsymbol{x}}}_{21}$ and ${{\boldsymbol{x}}}_{{\rm{IR}}}$ are column vectors of 21 cm and IR pixel values. Observe that in this approximation we model the radio covariance separately from the IR covariance. This separation is also ideal for our current study, where we seek to characterize any cross spectrum between the two, rather than model a cross spectrum and downweight by it.

As before, we use normalization ${\boldsymbol{M}}={{\boldsymbol{F}}}^{-1/2}$ in ${\boldsymbol{p}}={\boldsymbol{Mq}}$ and perform a final normalization using the peaks of the window functions ${\boldsymbol{W}}={\boldsymbol{MF}}$. We use the same IR covariance matrix given in Equation (21) and model the radio covariance matrix as solely due to sample variance. We construct this covariance using the first term on the right-hand side of Equation (21), taking the power spectrum of the residual MWA broadband images as our guess. We emphasize that in both cases random noise is subdominant to the foreground residuals, so we do not include photon shot noise or thermal noise in these covariances.

We apply this formalism to each of the four 4° deep ATLAS fields shown in Figure 1, pairing each IR image with the overlapping region of the MWA image cropped from the naturally weighted image. We crop the image space synthesized beam over the same field of view and then use it to apply uniform weighting to the cropped MWA image using Equation (12). In Figure 14 we plot the resulting cross spectrum (red markers) averaged over all four fields with $2\sigma$ error bars. Most are consistent with zero, and roughly half the estimated bandpowers are negative (triangles) and positive (circles), as expected if there is no measured correlation. Our bins are evenly spaced in $\mathrm{log}{\ell }$, and it is beneficial to overlap them slightly since our normalization matrix $M\sim {F}^{-1/2}$ ensures that bin errors are always uncorrelated. We take the tops of the error bars as 95% upper limits on the cross spectrum of residual foregrounds of 21 cm and Lyα emission from the EOR, achieving a tightest upper limit of ${{\rm{\Delta }}}^{2}\lt 181$ $({\text{kJy sr}}^{-1}\,{\rm{mK}})$ (95%) at ${\ell }\sim 800$. For comparison, we plot optimistic (black solid) and pessimistic (black dashed) model cross spectra as discussed in the previous section.

Zoom In Zoom Out Reset image size

We check these quadratic estimator results (red points) against FFT-based power spectrum results on the same images (gray points). In fact, the two spectra are quite similar. In most bins the quadratic estimator limits are lower but are not significantly closer to the fiducial 1% radio/IR foreground correlation. To understand why the FFT results are only slightly worse, consider that after our IR foreground masking only ∼5% of coarse pixels end up masked. With different masking parameters more or fewer coarse pixels would get masked, in some cases necessitating the quadratic estimator to deal with large masked regions. But with the best masking parameters we settled on, the improvement is modest.

Also in Figure 14 we plot the radio–infrared foreground cross spectrum that a 1% flux correlation would produce (green line), as motivated by our results in Sections 4.3 and 4.2. Such a slight geometric flux correlation is thus still predicted to be slightly below the sensitivity of our MWA–ATLAS experiment, but is it significant for other classes of experiments? To study this, we predict for a range of current and future radio–IR surveys both the foreground sample variance noise and the foreground cross spectrum for a slight geometric correlation.

We describe the radio and IR surveys we study here and summarize them side by side in Table 1.

Table 1.  Current and Future Radio and IR Surveys

Radio Survey	Resolution	${S}_{{\rm{\max }}}$
MWA	$6^{\prime}$	150 mJy
MWA (GMRT Catalog)	$6^{\prime}$	25 mJy
HERA (GMRT Catalog)	$12^{\prime}$	25 mJy
GMRT	$40^{\prime\prime}$	25 mJy
LOFAR	$20^{\prime\prime}$	3 mJy
SKA	$10^{\prime\prime}$	0.1 mJy
IR Survey	Field of View	${f}_{{\rm{IR}}}$ a
ATLAS	8°	1.0
WideATLAS	40°	1.0
WFIRST Mosaic	4°	0.01
DES Mosaic	40°	0.01

Note.

^aFraction of IR foreground power remaining relative to our ATLAS analysis.

Download table as:  ASCIITypeset image

Radio surveys

• 1.  
  MWA. We use the $6^{\prime}$ resolution permitted by our MWA analysis, corresponding to a maximum baseline length of 700 m and a survey depth at 50% completeness of 150 mJy (see Section 4.1).
• 2.  
  HERA. Primarily designed for 3D power spectrum measurements that are severely thermal noise limited, HERA is a compact redundant array, and even the outriggers only improve the resolution to $12^{\prime}$ (DeBoer et al. 2017). We assume that the GMRT catalog is used for source subtraction (see below).
• 3.  
  GMRT. We assume a maximum baseline length of 10 km, within which the uv plane is filled relatively uniformly; this corresponds to a $40^{\prime\prime}$ resolution. Intema et al. (2017) demonstrate a catalog depth of 25 mJy at 50% completeness.
• 4.  
  LOFAR. Similarly, we use a maximum baseline length of 20 km, to maintain a relatively filled uv plane, corresponding to a $20^{\prime\prime}$ resolution. Yatawatta et al. (2013) report a catalog depth of 3 mJy.
• 5.  
  SKA. Prandoni & Seymour (2015) report a nominal SKA-Low confusion-limited survey depth of 0.1 mJy at $10^{\prime\prime}$ resolution.

IR surveys

• 1.  
  ATLAS. Our four ATLAS stacking fields give a total field of view of 8°, though we are planning a future wider survey with a 40° field of view.
• 2.  
  WFIRST mosaic. Mitchell-Wynne et al. (2015) have demonstrated high-fidelity mosaicking using the technique of Fixsen et al. (2000), producing $10^{\prime}$ wide mosaics from the $\sim 2^{\prime}$ HST field of view. Along these lines, we assume that a 4° WFIRST field can be used through some combination of mosaicking (making larger IR images) and tiling (averaging the cross spectrum over many small fields). The instantaneous field of view of WFIRST is roughly 100 times that of HST with comparable resolution. Motivated by the results of Mitchell-Wynne et al. (2015), we assume a 100×  improvement in IR foreground removal over our ATLAS analysis.
• 3.  
  DES mosaic. Along these lines, we assume that a similar analysis can be conducted using the Dark Energy Survey (Dark Energy Survey Collaboration et al. 2016) over a field of view that is 10×  larger on a side, due to the instrument's correspondingly larger field of view. We assume the same improvement in foreground removal as in the HST mosaic discussed previously.

Neglecting the small geometric foreground correlation, the foreground sample variance noise due to uncorrelated foregrounds is

$$\begin{eqnarray}{{\rm{\Delta }}}_{{\rm{noise}}}^{2}({\ell }) & = & \sqrt{\displaystyle \frac{{{\rm{\Delta }}}_{{\rm{IR}}}^{2}({\ell }){{\rm{\Delta }}}_{{\rm{rad}}}^{2}({\ell })}{2{N}_{{\ell }}}}\\ & = & \sqrt{\displaystyle \frac{{{\rm{\Delta }}}_{{\rm{ATLAS}}}^{2}({\ell }){{\rm{\Delta }}}_{{\rm{MWA}}}^{2}({\ell }){f}_{{\rm{rad}}}{f}_{{\rm{IR}}}}{2{N}_{{\ell }}}}\\ & & \times \left(\displaystyle \frac{{\theta }_{{\rm{FOV}},{\rm{ATLAS}}}}{{\theta }_{{\rm{FOV}}}}\right)\left(\displaystyle \frac{d\theta }{d{\theta }_{{\rm{MWA}}}}\right),\end{eqnarray} \tag{ 24 }$$

where ${f}_{{\rm{IR}}}$ and ${f}_{{\rm{rad}}}$ are the fractions of IR and radio foreground power remaining relative to those in our ATLAS and MWA analyses. Using the empirical formula for radio source counts given in Di Matteo et al. (2002), ${f}_{{\rm{rad}}}$ scales as ${f}_{{\rm{rad}}}\propto {S}_{{\rm{\max }}}^{1.25}$, where ${S}_{{\rm{\max }}}$ is the flux limit of the subtraction. Scaling this relative to our MWA analysis gives ${f}_{{\rm{rad}}}={({S}_{{\rm{\max }}}/{S}_{{\rm{\max }},{\rm{MWA}}})}^{1.25}$.

Similarly, the foreground cross spectrum for a correlation c is given by

$$\begin{eqnarray}{{\rm{\Delta }}}_{{\rm{FG}}}^{2}({\ell }) & = & c\sqrt{{{\rm{\Delta }}}_{{\rm{IR}}}^{2}({\ell }){{\rm{\Delta }}}_{{\rm{rad}}}^{2}({\ell })}\\ & = & c\sqrt{{{\rm{\Delta }}}_{{\rm{ATLAS}}}^{2}({\ell }){{\rm{\Delta }}}_{{\rm{MWA}}}^{2}({\ell }){f}_{{\rm{rad}}}{f}_{{\rm{IR}}}}.\end{eqnarray} \tag{ 25 }$$

In Figure 15 we plot the foreground sample variance noise versus foreground cross spectrum for a range of current and future radio and infrared survey pairs. We calculate these for the ${\ell }=400$ bin and assume a 1% geometric foreground correlation. The diagonal black line separates the regions within which each of these effects dominates. We observe that our MWA/ATLAS analysis (black circle) is the only experiment limited by sample variance noise of uncorrelated foregrounds. An improved analysis using the GMRT catalog for radio source subtraction from MWA images and widening the ATLAS survey to a wide 40° (black star) falls barely in the foreground correlation-limited regime. A similar analysis using HERA instead of the MWA (yellow star) produces similar results, as HERA is not optimized for the wide-field, high-resolution radio surveys that broadband correlation experiments require.

Zoom In Zoom Out Reset image size

LOFAR reaches within a factor of a few in power of the optimistic EOR cross-spectrum model when correlated against the HST (red triangle) or DES (red square) mosaics, and the SKA coupled with the same IR surveys achieves a near detection of the optimistic model (white square and white triangle). None of the experiments reach the pessimistic model, but beginning to probe optimistic models already starts to exclude late reionization scenarios, and performing the same correlation experiment at higher redshifts can likely give more leverage on earlier ones.

Interestingly, we note that virtually all the intermediate experiments between this work and the ultimate detection will be limited by the geometric foreground correlations. Motivated by our results in Section 4.2, we have assumed here that the coherence remains at the percent level irrespective of mask depth, though more work is needed to model this effect more completely and self-consistently. In any case, this geometric foreground correlation appears to be strictly positive, while the predicted correlation of 21 cm–Lyα from the EOR is negative, giving us a test of whether foregrounds remain in the data.