Cross-correlating Carbon Monoxide Line-intensity Maps with Spectroscopic and Photometric Galaxy Surveys

Table 1 describes the current anticipated parameters for the initial phase (or Phase I) of COMAP. The receiver is currently undergoing commissioning at the Owens Valley Radio Observatory (OVRO) in California, where we expect the Phase I instrument to undertake a two-year observing campaign.

While a wide range of predictions exist for the CO power spectrum at z ∼ 3 (Righi et al. 2008; Visbal & Loeb 2010; Pullen et al. 2013; Breysse et al. 2014; Li et al. 2016; Padmanabhan 2018), the sensitivity calculations made in Li et al. (2016)—given their fiducial model—place COMAP Phase I squarely in a regime where instrumental noise dominates over sample variance, which is still true after various changes to COMAP parameters made since the writing of Li et al. (2016). This dictates the optimal observing strategy to some extent, pushing COMAP toward surveying, at most, several small fields (as close as possible to the ∼1 deg² field of view) with maximum observing efficiency. In a 2D analysis assuming a total on-sky time of one year (∼9000 hr) split across four patches (for ∼2200 hr per patch), Breysse et al. (2014) found that a survey footprint of four patches with almost 4 deg² per patch would maximize the total signal-to-noise ratio (S/N). If the optimal area scales linearly with on-sky time, the fiducial area per patch of 2.5 deg² is close to ideal for a survey time of 1500 hr per patch as assumed in this work.

Conventional galaxy surveys are a natural target for cross-correlation with LIM, and the successful detection of 21 cm line emission from galaxies at z ∼ 1 comes from cross-correlation with spectroscopic galaxy surveys (Chang et al. 2010; Switzer et al. 2013). However, spectroscopic data are currently limited in depth and abundance at z ∼ 3, so we look to existing public photometric data sets.

The COSMOS2015 catalog contains half a million galaxies observed in $1\lt z\lt 6$ across 1.58 square degrees of sky near the celestial equator. The catalog is K_s-selected (the K_s band being at 2.2 mm), and as mentioned previously, the completeness limit is K_s = 24.0. The K_s magnitude correlates well with stellar mass up to z ∼ 4 (and magnitudes in longer-wavelength bands may be used at higher redshifts; see, e.g., Davidzon et al. 2017). The redshift distribution skews largely toward lower redshift, but the source abundances are still relatively high for the redshift range relevant to COMAP, within which we find just under 20,000 sources over 1.58 square degrees (of which 0.2 square degrees are masked due to saturated pixels) with K_s ≤ 24.

The critical limiting factor of the COSMOS2015 catalog for studies of 3D large-scale structure is the redshift accuracy. Laigle et al. (2016) quote photometric redshift errors at 3 < z < 6 to be ${\sigma }_{z}=0.021(1+z)$, with some fraction of catastrophic failures; certain subsets even reach ${\sigma }_{z}\,\lesssim 0.01(1+z)$. However, Davidzon et al. (2017) suggest that the error is higher for z ≳ 3 galaxies, and is closer to ${\sigma }_{z}=0.03(1+z)$.

Deep low- to medium-resolution spectroscopic follow-up exists in the COSMOS field, but the surveys either do not satisfactorily cover z > 2 or are limited in area. A recent catalog of ten thousand objects selected across the COSMOS field (Hasinger et al. 2018) only contains ∼10² objects in the redshift range of interest to COMAP, with the majority of spectroscopic redshifts well below (or above) that range. Meanwhile, the VIMOS Ultra-Deep Survey (VUDS) (Le Fèvre et al. 2015) reports one of the largest z > 2 emission-line galaxy samples, with ∼2800 spectroscopic redshifts at z = 2.5–3.5 down to i_AB ≃ 25 over a square degree of sky. However, the initial data release (the only public data release, at time of writing; see Tasca et al. 2017) covers less than 10% of this area, and even the full square degree is split across three patches covering ∼10³ square arcminutes each, including half a square degree of the COSMOS field, or around a third of the COSMOS2015 coverage (and a fifth of the expected area per COMAP patch). Such surveys are well-suited for measuring stellar mass functions and average spectral properties, but the areas covered are less than ideal for cross-correlation against line-intensity maps. Uncertainty in the resulting cross spectra is roughly proportional to the inverse square root of the survey volume, so factors of 3–5 in sky area can noticeably affect detection significance.

The central stated goal of the HETDEX survey is constraining the expansion history of the universe, and specifically detecting dark energy at 3σ significance, by identifying ∼10⁶ LAEs through a wide-field spectroscopic survey (Hill et al. 2008; Hill & HETDEX Consortium 2016). However, we note a few interesting differences between HETDEX and previous conventional spectroscopic galaxy surveys measuring dark energy, e.g., the Dark Energy Spectroscopic Instrument (DESI) (DESI Collaboration et al. 2016), the SDSS-IV Extended Baryon Oscillation Spectroscopic Survey (eBOSS) (Dawson et al. 2016), and WiggleZ (Drinkwater et al. 2018).

• 1.  
  HETDEX targets a redshift range of z = 1.9–3.5,¹⁰ well beyond the typical redshifts of z ∼ 1 of other dark-energy-centric optical and NIR surveys. (eBOSS and DESI will target quasars, and thus the Ly-α forest at z ≳ 2, but emission line galaxies only up to z ∼ 2.)
• 2.  
  HETDEX does not target specific points on the sky based on prior imaging, but rather samples its survey footprint with integral field spectroscopy, integrating for ∼20 minutes at each spot in the sky, and picks sources out from the noisy spectra.

The first point is of interest to us because, unlike many other dark-energy-centric surveys, future HETDEX detections will fall squarely in a redshift range relevant to COMAP Phase I. The second point is of interest because, in principle, the noisy spectra uniformly sampled across the survey footprint could be processed and analyzed as a line-intensity cube.

Due to the survey and instrument design of HETDEX, the survey footprint will not be completely filled in with data. Rather, the integral field unit (IFU) arrangement of the VIRUS instrument covers only 1/4.5 of the area of each 20'-diameter "shot," and while some dithering (in three exposures) will fill in areas between fibres, no attempt will be made to fill in the IFU spacing of 100'', given a greater need for survey volume than for capture of small scales. The patch used for the Spitzer/HETDEX Exploratory Large Area (SHELA) survey (Papovich et al. 2016) is an exception, and the area between IFUs will be filled in for this field only (Hill & HETDEX Consortium 2016). At 13 × 2 square degrees, the SHELA patch could entirely contain a single (appropriately oriented) COMAP patch.

HETDEX is subject to interloper emission, detecting ∼10⁶ galaxies from z < 0.5 emitting in the [O ii] doublet (Hill et al. 2008). When extracting individual emitters from LAE survey spectra, imposing a minimum equivalent width cutoff removes the low-z emitters (Cowie & Hu 1998; Adams et al. 2011), and more sophisticated classification using Bayesian methods may also recover more of the underlying LAE sample (Leung et al. 2017). While we find no literature discussing foreground removal strategies in the context of LIM with HETDEX, such literature does exist in the context of [C ii] observations (Cheng et al. 2016; Lidz & Taylor 2016; Sun et al. 2018) and some strategies may be applicable beyond their original context. Furthermore, the low-z [O ii] emission will have no corresponding component in COMAP data, potentially leading to its amelioration in COMAP–HETDEX line-intensity cross-correlation.

To conclude this section, Table 3 shows a summary of the coverage and redshift precision of all surveys discussed above. For simulation purposes, we will expand or truncate coverage as necessary to match the COMAP coverage, as explained in Section 3.3.

Table 3.  A Summary of Survey Coverage and Redshift Precision for all Experiments Considered Above

Experiment	Field size	Range of z	${\sigma }_{z}/(1+z)$	Source Selection	Source Count
or catalog	(deg2)				(per deg2 per Δ z = 1)
COMAP	2.5	2.4–3.4	∼1/2000	none (surveys aggregate CO emission)	...
COSMOS2015	1.58	1–6	∼0.02	Ks-band magnitude ≲24.0	∼13000
HETDEX	300 + 150	1.9–3.5	∼1/700	Lyα luminosity ≳3 × 1042 erg s−1	∼1400
(in SHELA)	(26)	...	...	...	∼6000

Note. The COMAP survey footprint will comprise two or more patches of 2.5 square degrees each; the HETDEX survey footprint includes a 300 deg² "Spring" field and a 150 deg² "Fall" field, with possible 50–60% extensions to each. The "HETDEX in SHELA" row describes HETDEX full-fill coverage of the 13 × 2 deg² SHELA field.

Download table as:  ASCIITypeset image

We use a cosmological N-body simulation as the basis for our simulations. In particular, we use the c400-2048 box, which is part of the Chinchilla suite of dark-matter-only simulations. Li et al. (2016), who used the same simulation, provide implementation details of the simulation and subsequent halo identification. The simulation spans 400h⁻¹ Mpc on each side, and has a dark matter particle mass of $5.9\times {10}^{8}{h}^{-1}\,{M}_{\odot };$ we include dark matter halos more massive than ${M}_{\mathrm{vir}}={10}^{10}\,{M}_{\odot }$ in our analysis, meaning that we assume halos with lower virial halo mass are not massive enough to host galaxies with substantial star formation activity. (Li et al. 2016 justify the same choice of cutoff mass for CO simulations in their Appendix A; we consider its effect on Lyα simulations in our Appendix A alongside other details of Lyα modeling.)

To simulate galaxies in our field of observation, we use dark matter halos identified in "lightcone" volumes, enclosing all halos within a given sky area and redshift range, with each lightcone based on arbitrary choices of observer origin and direction within the cosmological simulation. We use 100 lightcones spanning z = 1.5–3.5 and a flat-sky area of 100' × 100', each populated with ∼10⁶–10⁷ halos. These lightcones form the basis for our simulated observations in z = 2.4–3.4. Note that this redshift range spans approximately 1 Gpc, so some line-of-sight repetition of the N-body data will occur, as we exploit the periodic boundary conditions of the simulation box—which is only $400{h}^{-1}\approx 570$ Mpc along each side—to extend the lightcone beyond the actual simulated comoving volume. However, the lightcone extents are not so much greater than the simulation volume that we expect this periodicity to impact the results of our study.

We derive CO and Lyα line luminosities for each of the halos (above the 10¹⁰ M_⊙ cutoff mass) from their virial masses and redshifts, with Figure 1 showing the average halo mass–line luminosity relations at z = 2.8 (the COMAP mid-band redshift). In addition, we calculate stellar masses (for the galaxy survey selection) and SFRs (as an intermediate property for the line luminosities) for each halo. Below, we explain the derivation of each of these properties.

Zoom In Zoom Out Reset image size

Stellar mass—We assign a stellar mass M_* to each halo using the best-fit stellar mass–halo mass relation from Behroozi et al. (2013a, 2013b). We apply the mean relation and the redshift-dependent scatter from the model, which is ≈0.23 dex at the redshifts considered here. The stellar mass is a property in itself, and unlike the SFR, it does not influence other properties.

Star formation rate—We convert halo masses to SFRs for each halo via interpolation of data from Behroozi et al. (2013a, 2013b). The main focus of these papers is to constrain the stellar mass–halo mass relation and derived quantities by comparing simulation data with observational constraints, and the resulting data include the average SFR in a halo given its mass and redshift.

We approximate halo-to-halo scatter in SFR by adding 0.3 dex log-normal scatter to the SFR obtained above, preserving the linear mean. The assumption of 0.3 dex scatter is reasonable, given the 0.2–0.4 true or intrinsic scatter observed in the SFR–stellar mass relation (Speagle et al. 2014; Salmon et al. 2015), combined with the tight ∼0.2 dex scatter in the stellar mass–halo mass relation of Behroozi et al. (2013a).

We also re-express SFR as infrared (IR) luminosity, using a known tight correlation:

$$\begin{eqnarray}&&\displaystyle \frac{\mathrm{SFR}}{{M}_{\odot }\,{\mathrm{yr}}^{-1}}={10}^{-10}\left(\displaystyle \frac{{L}_{\mathrm{IR}}}{{L}_{\odot }}\right).\end{eqnarray} \tag{ 1 }$$

As in Behroozi et al. (2013a) and Li et al. (2016), we assume a Chabrier initial mass function (IMF) (Chabrier 2003).

CO luminosity—We convert between IR luminosity and observed CO luminosity through power-law fits to observed data, commonly given in the literature:

$$\begin{eqnarray}&&\mathrm{log}\left(\displaystyle \frac{{L}_{\mathrm{IR}}}{{L}_{\odot }}\right)=\alpha \mathrm{log}\left(\displaystyle \frac{{L}_{\mathrm{CO}}^{{\prime} }}{{\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{pc}}^{2}}\right)+\beta ,\end{eqnarray} \tag{ 2 }$$

where for our fiducial model, we take α = 1.37 and β = −1.74 from a fit to high-redshift galaxy data (z ≳ 1) given in Carilli & Walter (2013), following Li et al. (2016).

In this context, ${L}_{\mathrm{CO}}^{{\prime} }$ (or indeed, any ${L}_{\mathrm{line}}^{{\prime} }$) is the observed luminosity (or velocity- and area-integrated brightness temperature) of the halo, which we convert into an intrinsic luminosity for each halo, as in Li et al. (2016):

$$\begin{eqnarray}&&\displaystyle \frac{{L}_{\mathrm{line}}}{{L}_{\odot }}=4.9\times {10}^{-5}{\left(\displaystyle \frac{{\nu }_{\mathrm{line},\mathrm{rest}}}{115.27\mathrm{GHz}}\right)}^{3}\displaystyle \frac{{L}_{\mathrm{line}}^{{\prime} }}{{\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{pc}}^{2}}.\end{eqnarray} \tag{ 3 }$$

We also add 0.3 dex log-normal scatter in CO luminosity, again preserving the linear mean. We model this scatter as completely independent of the scatter in SFR. If we were simulating CO emission alone, as in Li et al. (2016), we could end up with the same log-normal distribution width with a total scatter of ${\sigma }_{\mathrm{tot}}={({\sigma }_{\mathrm{SFR}}^{2}/{\alpha }^{2}+{\sigma }_{{L}_{\mathrm{CO}}}^{2})}^{1/2}=0.37$ (with all σ in units of dex) on top of the mean SFR–L_CO relation. However, unlike in Li et al. (2016), the scattered SFR for a given halo informs both CO and Lyα line luminosities. Furthermore, our approach to scatter is subtly different from adding log-normal scatter to CO luminosity for a given halo mass and redshift; see Appendix B for details.

Note that, while the scatter exhibited is representative of the amount of scatter seen in the high-redshift galaxy data obtained via Carilli & Walter (2013), it may not be representative of the galaxy population at large that COMAP will study, and a larger scatter in CO luminosity than we have assumed here may reduce our ability to cross-correlate CO against galaxy catalogs. However, this means that, were COMAP to indeed confidently detect CO while finding little in the way of CO-galaxy cross-correlation, that in itself would lead to interesting insights about how stochastic CO is in these high-redshift galaxies (or at least, those galaxies with sufficiently high mass and luminosity to be cataloged in the cross-correlation sample).

Lyα luminosity—Before dust absorption and other attenuation mechanisms, the Lyα line is 8.7 times stronger (for case B recombination) than Hα, a frequently chosen emission-line tracer of star formation activity (Kennicutt & Evans 2012). We use an intrinsic Lyα–SFR calibration (via Hα) with an escape fraction encapsulating possible attenuation of the intrinsic Lyα luminosity. We outline the specifics behind our model in Appendix A, the end result of which is that, for a given halo, we calculate

$$\begin{eqnarray}&&{L}_{\mathrm{Ly}\alpha }=1.6\times {10}^{42}\left(\displaystyle \frac{\mathrm{SFR}}{{M}_{\odot }\,{\mathrm{yr}}^{-1}}\right){f}_{\mathrm{esc}}(\mathrm{SFR},z)\,\mathrm{erg}\,{{\rm{s}}}^{-1},\end{eqnarray} \tag{ 4 }$$

with

$$\begin{eqnarray}\begin{array}{rcl}{f}_{\mathrm{esc}}(\mathrm{SFR},z) & = & {\left(1+{e}^{-1.6z+5}\right)}^{-1/2}\\ & & \times {\left[0.18+\displaystyle \frac{0.82}{1+0.8{\left(\tfrac{\mathrm{SFR}}{{M}_{\odot }{\mathrm{yr}}^{-1}}\right)}^{0.875}}\right]}^{2}.\end{array}\end{eqnarray} \tag{ 5 }$$

The escape fraction f_esc increases with lower SFR and higher redshift, so as to allow this model to match observed LAE luminosity functions (LFs) (including those from Sobral et al. (2017) and Gronwall et al. (2007); see Appendix A for more references) in the redshift range of interest (z ∼ 2–4).

As with the CO luminosity, we add 0.3 dex log-normal scatter in Lyα luminosity. Due to the non–power-law nature of the SFR–L_Lyα relation, the net scatter in L_Lyα varies non-monotonically with halo mass (between 0.31 and 0.42 dex in the mean relation shown in Figure 1 at z = 2.8).

Note that our model is tuned to LAE observations and assigns luminosities only at the centers of dark matter halos; thus, it does not account for the finer details of Lyα radiative transfer beyond each LAE, per se—which would result, for example, in diffuse Lyα halos or blobs (Steidel et al. 2011). Incorporating such details would introduce additional components to the Lyα emission, which we discuss in more detail in Appendix C. The signal as simulated here could be modified by these components in ways that may be detectable in HETDEX data through LIM, but not necessarily through LAE identification. All of this suggests the need to study implications of diffuse Lyα emission for HETDEX and cross-correlation with COMAP, but we leave this for future work.

After the processing outlined in the last section, each halo has a sky position, redshift (excluding peculiar velocities, which have minimal effects on the results of this work¹¹ ), virial halo mass, stellar mass, CO luminosity, and Lyα luminosity. We now use these properties to simulate survey data for COMAP Phase I, a COSMOS2015-like mass-selected galaxy catalog, and a HETDEX-like Lyα survey.

For ease of analysis, all survey cubes are generated with the same grid of voxels, based on the COMAP observation. The angular extent of each voxel is ${\delta }_{x}={\delta }_{y}=0\buildrel{\,\prime}\over{.} 4$ or 1.16 × 10⁻⁴ rad in each direction (oversampling the COMAP beamwidth by a factor of 10 and the on-sky VIRUS IFU width by a factor of 2.1), and each voxel spans ${\delta }_{\nu }=15.625\,\mathrm{MHz}$ in COMAP frequency (equivalent to 335 GHz in HETDEX frequency) unless otherwise specified. The simulated cubes span 100' × 100' and a continuous 26–34 GHz band in COMAP frequency (557–729 THz in HETDEX frequency), enclosing a total comoving volume of 190 × 190 × 1000 = 3.6 × 10⁷ Mpc³.

3.3.1. CO Intensity Data

3.3.2. Galaxy Overdensity Field

We devise an ideal NIR-selected galaxy survey tracing galaxies down to a certain stellar mass limit. We claim that the galaxy–halo connection allows us to model this, starting with a catalog of halos and imposing stellar mass cuts corresponding to realistic magnitude limits.

Mass Completeness—To crudely simulate the K_s magnitude cut used by catalogs like COSMOS2015, we assume the K_s magnitude correlates reliably with the stellar mass in our redshift range, an assumption that Laigle et al. (2016) support—at least in relating completeness limits for the two quantities. Each step down in magnitude is a factor of 10^0.4 up in brightness, so a constant mass-to-light ratio would result in the same factor up in stellar mass. Laigle et al. (2016) find their ${K}_{s,\mathrm{lim}}=24.0$ limit to be equivalent to a stellar mass completeness limit of ${M}_{* ,\mathrm{lim}}={10}^{10}\,{M}_{\odot }$ in our redshift range. We extrapolate this to different completeness limits with the following relation:

$$\begin{eqnarray}&&\mathrm{log}{M}_{* ,\mathrm{lim}}=10.0-0.4({K}_{s,\mathrm{lim}}-24.0).\end{eqnarray} \tag{ 8 }$$

The limits of ${K}_{s,\mathrm{lim}}=(25.0,24.0,23.0,22.0)$ are then equal to $\mathrm{log}({M}_{* }/{M}_{\odot })=(9.6,10.0,10.4,10.8)$. We use these stellar mass cuts to select our mock galaxy survey sample in each lightcone. Of these, the $\mathrm{log}({M}_{* }/{M}_{\odot })\gt 10.0$ cut matches the COSMOS2015 source abundance within a factor of order unity, so we take this as our fiducial M_* cut.

Redshift accuracy—To simulate uncertainty in redshifts derived from imagery, we apply different levels of scatter in observed redshift relative to the true cosmological redshift. While we do simulate cross-correlations against a survey with perfect galaxy redshift knowledge, not even spectroscopic surveys have such information. Therefore, we simulate normal scatter of redshifts with ${\sigma }_{z}/(1+z)=0.0007$, 0.003, 0.01, 0.02, and 0.03. The first scenario meets the minimum redshift accuracy required for cosmological applications of the Subaru Prime Focus Spectrograph (PFS) (Takada et al. 2014). The second scenario corresponds to lower-resolution spectroscopy, as seen in HST grism surveys like 3D-HST (${\sigma }_{z}/(1+z)=0.003$) (Momcheva et al. 2016), prism surveys like PRIMUS (${\sigma }_{z}/(1+z)=0.005$) (Coil et al. 2011; Cool et al. 2013), or even narrow-band photometric surveys like the PAU Survey (${\sigma }_{z}/(1+z)=0.0037$) (Eriksen et al. 2019). The last three scenarios represent optimistic, fiducial, and pessimistic expectations for photometric redshift accuracy, based on the discussion in Section 2.2.

After applying the stellar mass cut and redshift scatter (if applicable), we calculate the galaxy overdensity across the voxel grid. We count the number of galaxies ${N}_{\mathrm{gal},\mathrm{vox}}({\boldsymbol{x}})$ in each voxel, divide by the comoving volume of the voxel to get the number density ${n}_{\mathrm{gal},\mathrm{vox}}={N}_{\mathrm{gal},\mathrm{vox}}/{V}_{\mathrm{vox}}$. The quantity we deal with then is normalized by the average number density ${\bar{n}}_{\mathrm{gal}}$ across all voxels observed:

$$\begin{eqnarray}&&{\delta }_{\mathrm{gal},\mathrm{vox}}=\displaystyle \frac{{n}_{\mathrm{gal},\mathrm{vox}}({\boldsymbol{x}})}{{\bar{n}}_{\mathrm{gal}}}-1.\end{eqnarray} \tag{ 9 }$$

3.3.3. Lyα Survey Simulation

We have now defined a grid of voxels and four quantities associated with each voxel: the CO brightness temperature ${T}_{\mathrm{CO}}$, the mass-selected galaxy overdensity ${\delta }_{\mathrm{gal}}$ (for four different mass cuts), the Lyα spectral intensity per log-frequency interval ${\nu }_{{\rm{L}}y\alpha }{I}_{\nu ,\mathrm{Ly}\alpha }$, and the luminosity-selected LAE overdensity ${\delta }_{\mathrm{LAE}}$ (for two different luminosity cuts). Following previous works (Visbal & Loeb 2010; Li et al. 2016), we use Fourier estimators of the auto and cross power spectra. If $\tilde{A}({\boldsymbol{k}})$ and $\tilde{B}({\boldsymbol{k}})$ are the Fourier transforms of the fields $A({\boldsymbol{x}})$ and $B({\boldsymbol{x}})$, then the full 3D auto spectra are

$$\begin{eqnarray}&&{P}_{{\rm{A}}}({\boldsymbol{k}})={V}_{\mathrm{surv}}^{-1}{\left|\tilde{A}({\boldsymbol{k}})\right|}^{2},\quad {P}_{{\rm{B}}}({\boldsymbol{k}})={V}_{\mathrm{surv}}^{-1}{\left|\tilde{B}({\boldsymbol{k}})\right|}^{2};\end{eqnarray} \tag{ 10 }$$

and the full 3D cross spectrum is

$$\begin{eqnarray}&&{P}_{{\rm{A}}\times {\rm{B}}}({\boldsymbol{k}})={V}_{\mathrm{surv}}^{-1}\mathrm{Re}[\tilde{A}({\boldsymbol{k}}){\tilde{B}}^{* }({\boldsymbol{k}})].\end{eqnarray} \tag{ 11 }$$

Because the fields are defined on a grid of voxels at discrete values of ${\boldsymbol{x}}$, the Fourier transforms and full 3D spectra are also defined at discrete ${\boldsymbol{k}}$.

With the assumption of isotropy, we then spherically average the power spectra in shells of $k=| {\boldsymbol{k}}|$ of width ${\rm{\Delta }}k=0.035$ Mpc⁻¹, each containing some number of modes N_modes(k) for which $k-{\rm{\Delta }}k/2\lt | {\boldsymbol{k}}| \lt k+{\rm{\Delta }}k/2$, to obtain the spherically averaged 3D power spectra P_A(k), P_B(k), and ${P}_{{\rm{A}}\times {\rm{B}}}(k)$. Here, we take $A({\boldsymbol{x}})={T}_{\mathrm{CO}}$, while $B({\boldsymbol{x}})$ can be any of the other three fields defined above.

Broadly speaking, we can consider each of the auto and cross P(k) to be the sum of a clustering component that dominates at low k and a constant shot-noise term that dominates at high k. The clustering component traces the matter power spectrum with some bias associated with the quantity being observed, while the shot-noise component arises from Poisson fluctuations. Cross shot noise can have interesting interpretations explored in previous work: Wolz et al. (2017) show that H i-galaxy cross shot noise may be used to infer the H i content of the cross-correlated galaxies, and Breysse & Rahman (2017) show that ¹²CO–¹³CO cross shot noise within COMAP may be used to learn about ¹²CO–¹³CO isotopologue ratios and ¹²CO saturation. One could extend the latter idea to cross-correlate between two separate line-intensity surveys, e.g., between COMAP and HETDEX, in order to learn about the molecular fraction of Lyα emitters (using CO as a proxy for molecular gas). However, as we will see, such astrophysical interpretation of the cross shot noise must account for attenuation both from sparse sampling, as seen in HETDEX, and from redshift errors in all galaxy samples. We leave a detailed investigation into effects of such attenuation on astrophysical inferences for future work.

For our purposes, we take the sources of uncertainty to be sample variance, thermal noise in the CO temperature or Lyα intensity field, and shot noise in the galaxy density field.¹² From Visbal & Loeb (2010),

$$\begin{eqnarray}&&{\sigma }_{{P}_{{\rm{A}}\times {\rm{B}}}}^{2}(k)=\displaystyle \frac{{P}_{{\rm{A}},\mathrm{total}}(k){P}_{{\rm{B}},\mathrm{total}}(k)+{P}_{{\rm{A}}\times {\rm{B}}}^{2}(k)}{2{N}_{\mathrm{modes}}(k)},\end{eqnarray} \tag{ 12 }$$

where the "total" power spectra include noise, interloper emission, and other components not necessarily correlated to the tracer. We ignore those components in this work, apart from instrumental noise:

$$\begin{eqnarray}&&{P}_{\mathrm{CO},\mathrm{total}}(k)={P}_{\mathrm{CO}}(k)+{P}_{n,\mathrm{COMAP}};\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{P}_{\mathrm{Ly}\alpha ,\mathrm{total}}(k)={P}_{\mathrm{Ly}\alpha }(k)+{P}_{n,\mathrm{HETDEX}}.\end{eqnarray} \tag{ 14 }$$

Therefore, it is best to treat the signal-to-noise estimates given in this work as upper bounds on what the actual surveys may ultimately achieve. In general, the instrumental (thermal) noise power spectrum is given by the root-mean-square temperature or intensity fluctuation per voxel ${\sigma }_{n}$ and the comoving voxel volume ${V}_{\mathrm{vox}}$; it is assumed to be pure white noise, and thus constant across all k (Lidz et al. 2011):

$$\begin{eqnarray}&&{P}_{n}={\sigma }_{n}^{2}{V}_{\mathrm{vox}}.\end{eqnarray} \tag{ 15 }$$

(This is analogous to the inverse of the weight per solid angle $w={({\sigma }_{\mathrm{pix}}^{2}{{\rm{\Omega }}}_{\mathrm{pix}})}^{-1}$ in the calculation of uncertainties in 2D ${C}_{{\ell }}$ analysis from Knox 1995.) We then need only calculate ${\sigma }_{n,\mathrm{COMAP}}$ and ${\sigma }_{n,\mathrm{HETDEX}}$ based on the expected instrumental and survey parameters.

Calculation of the COMAP instrumental noise follows the same procedure outlined in Appendix C of Li et al. (2016), using the parameters in Table 1. Specifically, ${\sigma }_{n,\mathrm{COMAP}}$ derives from the system temperature ${T}_{\mathrm{sys}}=40$ K, the number of feeds ${N}_{\mathrm{feeds}}=19$, the frequency resolution ${\delta }_{\nu }=15.625\,\mathrm{MHz}$, and the survey time per pixel ${\tau }_{\mathrm{pix}}=(1500\mathrm{hr})\cdot {\delta }_{x}{\delta }_{y}/(2.5{\deg }^{2})$ (being simply the total survey time per patch divided by the number of pixels per patch):

$$\begin{eqnarray}&&{\sigma }_{n,\mathrm{COMAP}}=\displaystyle \frac{{T}_{\mathrm{sys}}}{\sqrt{{N}_{\mathrm{feeds}}{\delta }_{\nu }{\tau }_{\mathrm{pix}}}}.\end{eqnarray} \tag{ 16 }$$

We simulate and assume only one patch of 2.5 deg², so the signal-to-noise estimates are also given per patch. If we fix the on-sky time per patch and the solid angle per patch, uncertainties from COMAP instrumental noise (which we expect to dominate total uncertainties) will decrease as the square root of the number of patches (from the linear increase in the number of modes averaged to obtain our best P(k) estimate).

For the HETDEX instrumental noise, we refer to the sensitivity metrics given in Hill et al. (2014). Each VIRUS fiber covers a solid angle of 1.8 square arcseconds at one time, and a dither pattern of three exposures allows the area within each IFU to be completely covered. The line sensitivity expected from each 20 minute shot is ≲4 × 10⁻¹⁷ erg s⁻¹ cm⁻², and dividing by the 5.4 square arcsecond solid angle per dithered fiber gives ${\sigma }_{n,\mathrm{HETDEX}}\,=3.15\times {10}^{-7}$ erg s⁻¹ cm⁻² sr⁻¹. Thus, we now have ${P}_{n,\mathrm{COMAP}}\,={\sigma }_{n,\mathrm{COMAP}}^{2}{V}_{\mathrm{vox}}$ and ${P}_{n,\mathrm{HETDEX}}={\sigma }_{n,\mathrm{HETDEX}}^{2}{V}_{\mathrm{vox}}$. For COMAP, the dependence on the voxel size (simulated or otherwise) cancels out because ${\sigma }_{n}^{2}$ is proportional to the inverse of the frequency bandwidth per channel as well as the inverse of the solid angle per pixel (via ${\tau }_{\mathrm{pix}}$). For HETDEX, we choose ${V}_{\mathrm{vox}}$ in the context of ${P}_{n,\mathrm{HETDEX}}$ to correspond to the 5.4 square arcsecond solid angle per dithered fiber (canceling the same factor we used to convert from line sensitivity to intensity fluctuation per voxel) and the spectral resolution of the instrument (or a redshift interval of $(1+z)/R\sim 0.005$), which comes out to ∼0.03 Mpc³. As with the cross power spectrum uncertainty in Equation (12), the errors on the individual auto power spectra are also given by dividing the "total" spectra by the square root of N_modes(k). In the case of galaxy or LAE overdensities, the "total" power spectrum is equal to the simulated power spectrum, as we never subtract the shot noise term of $1/\bar{n}$. In the other cases, we use the "total" spectra from above:

$$\begin{eqnarray}&&{\sigma }_{{P}_{\mathrm{CO}}}(k)=\displaystyle \frac{{P}_{\mathrm{CO}}(k)+{P}_{n,\mathrm{COMAP}}}{\sqrt{{N}_{\mathrm{modes}}(k)}};\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&{\sigma }_{{P}_{\mathrm{Ly}\alpha }}(k)=\displaystyle \frac{{P}_{\mathrm{Ly}\alpha }(k)+{P}_{n,\mathrm{HETDEX}}}{\sqrt{{N}_{\mathrm{modes}}(k)}}.\end{eqnarray} \tag{ 18 }$$

When calculating S/N, we also need to account for attenuation in the CO signal due to the beam. As discussed in Li et al. (2016), the beam resolution limit attenuates the Fourier-transformed CO temperature field ${\tilde{T}}_{\mathrm{CO}}({\boldsymbol{k}})$ at each ${\boldsymbol{k}}$ by a factor of $\exp (-{k}_{\perp }^{2}{\sigma }_{\mathrm{beam}}^{2}/2)$, where ${k}_{\perp }$ is the transverse component of ${\boldsymbol{k}}$ and ${\sigma }_{\mathrm{beam}}$ the width of the Gaussian profile of the beam (projected into the comoving survey volume). However, angular resolution limits differ significantly between the CO temperature field and the field being cross-correlated against, and the latter (which may be due to pixelization, galaxy survey limits, fiber diameters, and so on) will be much finer than the COMAP ${\sigma }_{\mathrm{beam}}$. Thus, while the full 3D CO auto spectrum ${P}_{\mathrm{CO}}{({\boldsymbol{k}})\propto {\tilde{T}}_{\mathrm{CO}}({\boldsymbol{k}})}^{2}$ is attenuated at each ${\boldsymbol{k}}$ by $\exp (-{k}_{\perp }^{2}{\sigma }_{\mathrm{beam}}^{2})$, the cross spectra only scale linearly with the attenuated CO signal and consequently are attenuated at each ${\boldsymbol{k}}$ by approximately $\exp (-{k}_{\perp }^{2}{\sigma }_{\mathrm{beam}}^{2}/2)$.

The spherically averaged auto and cross P(k) are correspondingly attenuated by

$$\begin{eqnarray}&&{W}^{2}(k)=\langle \exp (-{k}_{\perp }^{2}{\sigma }_{\perp }^{2}){\rangle }_{{\boldsymbol{k}}},\end{eqnarray} \tag{ 19 }$$

where the average is over all (discrete) ${\boldsymbol{k}}$ in the $P({\boldsymbol{k}})$ averaging that fall within the shell corresponding to k, and ${\sigma }_{\perp }$ is the applicable resolution limit for the P(k) being calculated.¹³ For the CO auto spectrum, ${\sigma }_{\perp }$ is simply ${\sigma }_{\mathrm{beam}}$, and given the large beam size of the COMAP telescope, effectively ${\sigma }_{\perp }\approx {\sigma }_{\mathrm{beam}}/\sqrt{2}$ for the cross spectra as discussed above.

Recall that, toward the end of Section 3.3, we also discussed attenuation in cross spectra between CO and Lyα intensity fluctuations due to the limited spectral resolution of HETDEX. This follows an average similar to that in Equation (19), but with ${k}_{\parallel }$, the line-of-sight component of ${\boldsymbol{k}}$:

$$\begin{eqnarray}&&{W}_{z}^{2}(k)=\langle \exp (-{k}_{\parallel }^{2}{\sigma }_{\parallel }^{2}){\rangle }_{{\boldsymbol{k}}}.\end{eqnarray} \tag{ 20 }$$

Were we dealing with the HETDEX auto spectrum, we would base ${\sigma }_{\parallel }={\sigma }_{\parallel ,\mathrm{HETDEX}}$ on the redshift-space resolution of $0.0015(1+z)$ at the average redshift of the survey volume, again based on the resolving power $R\sim 700$ of VIRUS (Hill et al. 2014) and observed Lyα component velocity offsets (Chonis et al. 2013) as discussed in Section 3.3. However, in cross-correlation with COMAP, which has much higher redshift precision, we can assume that ${\sigma }_{\parallel }\approx {\sigma }_{\parallel ,\mathrm{HETDEX}}/\sqrt{2}$ using arguments similar to those for ${\sigma }_{\perp }$.

In discussing our results, we will often quote the total S/N over "all" scales, meaning all scales $k\in (0.017,4.2)$ Mpc⁻¹ (in linear bins of width 0.034 Mpc⁻¹) that our simulations nominally access (with the minimum and maximum k respectively corresponding to the lightcone and voxel angular widths in comoving space). Following Li et al. (2016), we calculate this total ${\rm{S}}/{\rm{N}}$ as

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{S}}}{{\rm{N}}}={\left[\displaystyle \sum _{k}{\left(\displaystyle \frac{{\rm{S}}}{{\rm{N}}}(k)\right)}^{2}\right]}^{1/2}={\left[\displaystyle \sum _{k}{\left(\displaystyle \frac{{P}_{\mathrm{obs}}(k)}{{\sigma }_{P}(k)}\right)}^{2}\right]}^{1/2},\end{eqnarray} \tag{ 21 }$$

where ${P}_{\mathrm{obs}}(k)=P(k){W}^{2}(k)$ (or $P(k){W}^{2}(k){W}_{z}^{2}(k)$ in the case of the CO-Lyα intensity cross spectrum) and the sum is over all k-bins with central values in the previously specified range of $(0.017,4.2)$ Mpc⁻¹. Note that ${W}^{2}(k)$ and ${W}_{z}^{2}(k)$ also modify the auto and cross spectra (before thermal noise) in the expressions for ${\sigma }_{P}(k)$.

We explore two different variables in the galaxy survey: (1) the mass completeness of a perfect redshift survey (${\sigma }_{z}/(1+z)=0$), and (2) the redshift accuracy of a survey with the fiducial mass completeness cut ($\mathrm{log}({M}_{* }/{M}_{\odot })\gt 10.0$).

4.1.1. Mass Completeness

Figure 2 shows how r(k) varies based on the mass completeness of an ideal survey with perfectly determined redshifts. Note that because the CO emission traces faint galaxies well below halo masses of ${10}^{12}\,{M}_{\odot }$ in between the brighter galaxies with ${M}_{\mathrm{vir}}\gtrsim {10}^{12}\,{M}_{\odot }$, r(k) falls off with higher k as the CO and galaxy surveys begin to trace less similar fluctuations at smaller comoving scales. The falloff is greater with less complete surveys, and it impacts all scales significantly once we reach ∼4 × 10³ galaxies in our survey volume (or densities of ∼10³ galaxies per deg² per ${\rm{\Delta }}z=1$).

Zoom In Zoom Out Reset image size

We show signal-to-noise ratios for all simulated cross spectra in Table 4. All of these cross spectra—even the one with the lowest assumed density—might be detected with a signal-to-noise ratio above 20, far higher than the 4.6 expected from the CO auto spectrum alone, at least provided that the galaxy sample fully covers the COMAP volume.

Table 4.  Mean over All Simulated Observations of 100 Lightcones of Total ${\rm{S}}/{\rm{N}}$ for ${P}_{\mathrm{CO}\times \mathrm{gal}}(k)$ over All Modes

$\mathrm{log}({M}_{* ,\min }/{M}_{\odot })$	${\sigma }_{z}/(1+z)$	Median Galaxy Count	${\rm{S}}/{\rm{N}}$
9.6	0.	$5.4\times {10}^{4}$	33.2
10.0	0.	$2.9\times {10}^{4}$	32.4
10.4	0.	$1.3\times {10}^{4}$	29.6
10.8	0.	$3.5\times {10}^{3}$	22.7
9.6	0.0007	$5.4\times {10}^{4}$	30.1
10.0	0.0007	$2.9\times {10}^{4}$	29.1
10.4	0.0007	$1.3\times {10}^{4}$	26.6
10.8	0.0007	$3.5\times {10}^{3}$	20.4
10.0	0.003	$2.9\times {10}^{4}$	19.1
10.0	0.01	$2.8\times {10}^{4}$	10.6
10.0	0.02	$2.7\times {10}^{4}$	7.32
10.0	0.03	$2.7\times {10}^{4}$	5.93

Note. For comparison, the ${\rm{S}}/{\rm{N}}$ for P_CO(k) is 4.6. All signal-to-noise ratios are quoted for a single patch of 2.5 deg² observed for 1500 hr; we may expect an improvement of up to a factor of $\sqrt{2}$ if two equivalent patches are observed for 1500 hr each.

Download table as:  ASCIITypeset image

4.1.2. Redshift Accuracy

Now fixing the mass completeness cut at the fiducial value of $\mathrm{log}({M}_{* ,\min }/{M}_{\odot })=10.0$, we vary the redshift accuracy in the survey. Figure 3 shows the resulting r(k) values for each value of ${\sigma }_{z}/(1+z)$ assumed.

Zoom In Zoom Out Reset image size

Note that, even for ${\sigma }_{z}/(1+z)=0.0007$, we find significant attenuation of cross-correlation at large k, suggesting the effect is particularly great for the cross shot-noise component of the signal (which dominates over the clustering component for $k\gtrsim 1$ Mpc⁻¹).¹⁴ It is realistic to expect attenuation of this kind given the possible fitting errors in spectroscopic redshift determination, as well as mismatches in line-of-sight pixelization or resolution between the two survey data. This level of precision is also where redshift-space distortions (RSD) from peculiar velocities of galaxies begin to distort line-of-sight structure. As this precision is thus sufficient for galaxy redshift surveys looking for RSD, they have little motivation to pursue higher spectral resolution—e.g., Gaztañaga et al. (2012) and Eriksen & Gaztañaga (2015) demonstrate that ${\sigma }_{z}/(1+z)=0.003$ is adequate for the purposes of the PAU Survey (mentioned above in Section 3.3.2).

For both ${\sigma }_{z}/(1+z)=0.0007$ and ${\sigma }_{z}/(1+z)=0.003$ (the higher- and lower-resolution spectroscopic errors), we see significant attenuation of the cross-correlation at large k, i.e., the smallest scales simulated, but relatively little attenuation at the largest scales probed. This works to our advantage, as our single-dish line-intensity survey aims to detect CO fluctuations at these largest scales rather than the CO shot noise. The S/N reflects this, falling only from 32.4 to 29.1 if ${\sigma }_{z}/(1+z)\,=0.0007$, and then to 19.1 if we increase ${\sigma }_{z}/(1+z)$ to 0.003.

These numbers become more discouraging as we approach errors more typical of wide-band photometric surveys, settling in a range closer to 6–11. We show the effect graphically in Figure 4, and again summarize the signal-to-noise ratios calculated in Table 4. Even taking these ratios at face value, the high ${\sigma }_{z}$ significantly dulls the advantage of cross-correlation over auto-correlation in detection significance. We must also carefully consider the integration time of 1500 hr per patch assumed for all scenarios. In the specific case of COMAP, when observing from the site in California, this integration time takes 2–3 times longer to achieve on an equatorial field like COSMOS versus on a field at 50–70° decl. If we had fixed the "real" survey duration for all scenarios instead of the integration time, we would expect to see no advantage in detection significance from cross-correlating against a COSMOS2015-like galaxy catalog over CO autocorrelation in a field of our choice.

Zoom In Zoom Out Reset image size

One natural step we might take to ameliorate the problem of photometric redshift errors is to coarsen the line-of-sight resolution of the data, so as to make the redshift errors less relevant. While this will result in boosting the signal closer to its true value, by essentially removing attenuated line-of-sight modes from consideration, it will also result in increased uncertainties in the end result as the Fourier-space volume—and thus the number of modes decreases. The net result is largely a loss in total cross signal-to-noise, as we show in Figure 5, and only a slight gain for photometric scenarios of ${\sigma }_{z}/(1+z)\geqslant 0.01$ (which plateaus when ${\sigma }_{z}/(1+z)\approx {\delta }_{\nu }/\nu$).

Zoom In Zoom Out Reset image size

We consider two variations on cross-correlating against a HETDEX-like survey: one in which we cross-correlate against a Lyα intensity cube, and one in which we cross-correlate against the LAE overdensity field with different luminosity cuts, as discussed in Section 3.3. We fix ${\sigma }_{z}/(1+z)=0.0015$ in all cases, however.

We first consider the COMAP×LAE scenario, and show the simulated r(k) in Figure 6. For luminosity cuts of $(3\times {10}^{42},6\times {10}^{42})$ erg s⁻¹, we find an average of $(1.4\times {10}^{4},4.2\times {10}^{3})$ LAEs in the survey volume, with approximately one-fourth as many LAEs when the volume is sparsely sampled with a fill factor of 1/4. The number of LAEs with ${L}_{\mathrm{Ly}\alpha }\gt 3\times {10}^{42}$ erg s⁻¹ approximately matches the expected source abundance in Hill et al. (2008) of 8 × 10⁵ LAEs across ∼400 (sparsely sampled) square degrees in a redshift interval of ${\rm{\Delta }}z=1.6$.

Zoom In Zoom Out Reset image size

If the HETDEX data filled the survey volume completely, as in the SHELA field, the COMAP×LAE cross spectrum would be detectable with total S/N as high as 20.7 for ${L}_{\mathrm{Ly}\alpha }\gt 3\times {10}^{42}$ erg s⁻¹, even with the LAE redshifts scattered by ${\sigma }_{z}/(1+z)=0.0015$ (without which the ${\rm{S}}/{\rm{N}}$ might be higher by around 30%). However, with the quarter-fill factor, the ${\rm{S}}/{\rm{N}}$ does drop to 14.5. The numbers are lower by 20%–25% for the more stringent cut of ${L}_{\mathrm{Ly}\alpha }\gt 6\times {10}^{42}$ erg s⁻¹.

We now consider cross-correlating against a Lyα intensity cube generated from HETDEX, showing r(k) in Figure 7. Compared to cross-correlation against LAE overdensity, the roll-off of r(k) with greater k is slower, even with sparse sampling and limited spectral resolution. The signal is potentially detectable at an S/N of 29.1 without sparse sampling, and an S/N of 23.3 with sparse sampling. This results in a slight edge versus cross-correlating against individually identified LAEs, although the simulated advantage may change with the Lyα model; see Section 5.2 (and Appendix C) for further discussion.

Zoom In Zoom Out Reset image size

We summarize the signal-to-noise ratios from COMAP–HETDEX cross-correlation (and expected LAE counts for applicable scenarios) in Table 5, and compare the LAE cross-correlation S/N graphically against S/N from cross-correlation against a mass-selected galaxy sample in Figure 8. Note that, unlike the COSMOS field, the HETDEX survey footprint is partly well-matched with areas of relatively high observing efficiency for COMAP. Therefore, a CO observing campaign with sufficient data (in one patch of several) for a $\sim 5\sigma$ CO auto detection could readily overlap with HETDEX to generate a $\sim 15\sigma$ detection in cross-correlation, per the signal-to-noise ratios in Table 5.

Zoom In Zoom Out Reset image size

Table 5.  Mean over All Simulated Observations of 100 Lightcones of Total ${\rm{S}}/{\rm{N}}$ for ${P}_{\mathrm{CO}\times \mathrm{Ly}\alpha }(k)$ over All Modes

${L}_{\mathrm{Ly}{\alpha },\min }$ (erg s−1)	Median LAE count	${\rm{S}}/{\rm{N}}$
none (LIM)	⋯	29.1 (23.3)
3 × 1042	$1.4\times {10}^{4}$ ($3.5\times {10}^{3}$)	21.2 (14.5)
6 × 1042	$4.2\times {10}^{3}$ ($1.1\times {10}^{3}$)	16.7 (10.6)

Note. Counts and ${\rm{S}}/{\rm{N}}$ in parentheses are for quarter-fill sparse sampling; counts and ${\rm{S}}/{\rm{N}}$ not in parentheses are for full-fill sampling. We assume ${\sigma }_{z}/(1+z)=0.0015$ in all cases. For comparison, the ${\rm{S}}/{\rm{N}}$ for P_CO(k) is 4.6. All ${\rm{S}}/{\rm{N}}$ are quoted for a single patch of 2.5 deg² observed for 1500 hr; we may expect to see an increase of up to a factor of $\sqrt{2}$ if two equivalent patches are observed for 1500 hr each, and a further, roughly linear increase with more integration time.

Download table as:  ASCIITypeset image

To conclude this section, we show a plot of all auto and cross P(k) with sensitivities in Figure 9. Note the shape of the CO auto P(k), which flattens beyond $k\sim 1$ Mpc⁻¹ as the shot-noise component of the power spectrum begins to dominate over the clustering component following the underlying matter distribution. Any such shot-noise component in the cross spectra is far less apparent, as expected from the random redshift errors wiping out smaller-scale correlations. (The exception is the CO×HETDEX LIM cross P(k) plotted, which does not incorporate this effect, as it is folded into the accompanying sensitivity curve instead. Thus, a shot-noise component is visible for this set of P(k).) This matches what we also demonstrate in the r(k) plots of Figures 3 and 6.

Zoom In Zoom Out Reset image size

Note also that the shapes of the sensitivity curves clearly show the impact of cross-correlating against data with significantly higher angular resolution than the COMAP data (and thus reducing ${\sigma }_{\perp }$ in Equation (19) by a factor of $\sqrt{2}$). We also plot the signal-to-noise ratio at each k for all spectra considered, which shows the same.

Note finally that we find a different total signal-to-noise ratio if we consider uncertainties on the anisotropic power spectrum P(k, μ) (where $\mu ={k}_{\parallel }/k$ is the cosine of the k-space spherical polar angle, using ${k}_{\parallel }$ to describe the line-of-sight component of the vector ${\boldsymbol{k}}$) instead of the spherically averaged P(k), which may be higher for galaxy samples with $\sigma {z}/(1+z)\geqslant 0.01$. However, the enhancement is not enough to alter the fundamental conclusions of our work; see Appendix E for further discussion.

Our treatment of the effects of redshift errors is considerably simplified from what we may expect from real-world survey data. In particular, we simulate photometric redshifts that lack both bias and catastrophic outliers, either of which could potentially further affect the cross-correlation signal.

Furthermore, photometric redshift algorithms usually compute a redshift probability density function (PDF) for each galaxy, which can be reduced to the point estimate that we consider. Indeed, Asorey et al. (2016) propose using the full redshift PDF for each galaxy in studies of angular galaxy clustering within bins of roughly twice the typical PDF widths.

However, the technique is of limited applicability when using LIM to probe 3D clustering. In general, techniques using galaxy redshift PDFs will not (and are not meant to) recover the true redshift of the galaxy or the true galaxy density fluctuations in the survey volume. In fact, using the redshift PDFs for each galaxy instead of the estimated redshifts will merely convolve the PDF (around the estimated redshift) with the scatter distribution of the estimated redshift (versus the true redshift). This can only result in additional line-of-sight smearing and thus attenuation of the signal.

Furthermore, the detectability of the line-intensity signal relies in large part on the decrease in uncertainty due to being able to access the line-of-sight fluctuation modes. In Figure 5, where we actually explore the idea of wider frequency channels (albeit not the idea of using full redshift PDFs), the signal-to-noise remains low despite recovery of the signal due to the loss of line-of-sight information.

Ultimately, we must contend with some suppression of the cross spectrum signal from errors in galaxy redshifts, even for spectroscopic redshifts. In principle, we might consider calculating a transfer function to compensate for this suppression. However, the signal being suppressed in the first place diminishes our confidence in its detection, regardless of whether or not we can undo the suppression in analysis.

Furthermore, an accurate transfer function will require accurate and precise characterization of redshift errors—both of the mean offset or bias, as well as of the variance or scatter. For comparison, when binning populations of galaxies in a photometric survey to gauge cosmological parameters, both bias and scatter in redshift must be characterized within 0.003 to limit significant degradation in error, even using bins of ${\rm{\Delta }}z\sim 0.1$ (Huterer et al. 2006; Ma et al. 2006). In practice, LIM can probe fluctuations across bins of ${\rm{\Delta }}z$ as fine as ∼0.001, although the science is often more astrophysical than cosmological, especially in deep surveys with small sky fractions.

Such requirements in and of themselves should not preclude the possibility of studies of shot noise cross-correlation of the kind posited by Wolz et al. (2017) (which proposes cross-correlation against a spectroscopic galaxy survey, for which the errors may be more precisely characterized than for photometric redshifts) and Breysse & Rahman (2017) (if extended to cross-correlations between line-intensity surveys). However, future studies may wish to take careful inventory of expected redshift errors. A thorough study of the feasibility of such techniques and the impact of these effects on shot noise detection significance is beyond the scope of this particular work, but still highly desirable as these cross-correlations become increasingly viable.

Recall that we consider two possible outputs of HETDEX: the LAE overdensity ${\delta }_{\mathrm{LAE},\mathrm{vox}}$ for all emitters above a certain luminosity cut, and the total Lyα line-intensity quantified as ${\nu }_{{\rm{L}}y\alpha }{I}_{\nu ,\mathrm{Ly}\alpha }$. One might expect that, because the latter includes Lyα emission from emitters below a realistic luminosity cut, the line-intensity cube should trace structures absent in the ${\delta }_{\mathrm{LAE},\mathrm{vox}}$ cube—and potentially significantly improve detectability in cross-correlation against T_CO.

However, our simulated cross-correlation of ${T}_{\mathrm{CO}}$ against ${\nu }_{{\rm{L}}y\alpha }{I}_{\nu ,\mathrm{Ly}\alpha }$ only provides a slim advantage (40–60%, based on the figures in Table 5) in signal-to-noise ratio over cross-correlation against LAEs detected with ${L}_{\mathrm{Ly}\alpha }\gt 3\times {10}^{42}$ erg s⁻¹. This seems significant enough until one considers that the signal-to-noise ratio is still about the same as a conventional spectroscopic galaxy survey with ${M}_{* ,\min }\sim {10}^{10}\,{M}_{\odot }$. Furthermore, in our model, LAEs with ${L}_{\mathrm{Ly}\alpha }\gt 3\times {10}^{42}$ erg s⁻¹ are typically in halos with ${M}_{\mathrm{vir}}\gtrsim {10}^{12}\,{M}_{\odot }$ (looking at Figure 1), and the line-intensity data should trace twice as much (or more) signal by including emission from lower-mass halos (looking at Figure 13 in Appendix A.4).¹⁵ Note first that this ignores both the additional information that may be captured in LIM (but is beyond the scope of our LAE-based simulations) and the additional challenges that would be inherent to LIM, including contamination from zodiacal light and interloper emission in [O ii] (the latter of which would be rejected in COMAP–HETDEX cross-correlation but would nonetheless contribute additional uncertainty about that cross signal).¹⁶ With these caveats in mind, we should ask why overcoming such challenges and observing the full LAE population would appear to result in so little improvement in forecast ${\rm{S}}/{\rm{N}}$.

One possible explanation may be the halo mass–line luminosity relations in Figure 1. The large step down in Lyα escape fraction around $\mathrm{SFR}\sim 1\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ or halo mass ${M}_{\mathrm{vir}}\sim {10}^{11}\,{M}_{\odot }$ means that the slope of the halo mass–${L}_{\mathrm{Ly}\alpha }$ relation begins to decline above this mass, before the turnaround at ${10}^{12}\,{M}_{\odot }$ seen in both this relation and in the halo mass–${L}_{\mathrm{CO}}$ relation. This would lead to differences between the CO and Lyα signals in the relative contributions of halos with ${M}_{\mathrm{vir}}\lesssim {10}^{11}\,{M}_{\odot }$ and those with ${M}_{\mathrm{vir}}\gtrsim {10}^{12}\,{M}_{\odot }$, which may adversely affect the cross-correlation of the line-intensity signals. However, whether only the most massive halos behave similarly enough to co-correlate would certainly be a highly model-dependent effect. Furthermore, the Lyα relation is based solely on observed LAE densities, which may not result in a complete model of Lyα emission (as discussed briefly in Section 3.2).

Another (not mutually exclusive) possibility is that LAE detection, and not Lyα intensity mapping, may actually be the optimal observation for HETDEX, which has low instrumental noise and high angular resolution—and thus, low source confusion. Cheng et al. (2018) provide an overview of optimal observations in different noise and confusion regimes using an analytic source model. Their work suggests that detecting individual sources rather than aggregate LIM is the optimal observation for HETDEX, as well as for the higher-redshift Lyα observations that would be possible with the Cosmic Dawn Intensity Mapper (CDIM) (Cooray et al. 2016), although not for the SPHEREx concept (Doré et al. 2014).

However, Cheng et al. (2018) are also careful to note that the Lyα model used only incorporates point sources (just as our own model relies only on LAE LFs) and does not account for the expected extended emission from radiative transfer beyond galaxies. We have also discussed this as a limitation of our model, and repeat our caveat from Section 3.2 that a need exists for future work on implications of extended Lyα emission for HETDEX and COMAP–HETDEX cross-correlation.

The conversion is based on assuming a certain intrinsic Hα luminosity per SFR and a Lyα/Hα line ratio of 8.7. Stellar synthesis modeling done in Kennicutt et al. (1994) (via Kennicutt 1998) suggested that, for a Salpeter IMF,

$$\begin{eqnarray}&&\displaystyle \frac{{L}_{{\rm{H}}\alpha }}{\mathrm{SFR}}=\displaystyle \frac{1.26\times {10}^{41}\,\mathrm{erg}\,{{\rm{s}}}^{-1}}{{M}_{\odot }\,{\mathrm{yr}}^{-1}}.\end{eqnarray} \tag{ 23 }$$

The calibration arises from models of stellar evolutionary tracks, ionizing emissivity, and recombination rates for gas with $T={10}^{4}$ K; see Kennicutt et al. (1994) and Hummer & Storey (1987) for the last item.

Murphy et al. (2011) update this calibration with revised stellar synthesis models incorporating a Kroupa IMF, yielding a SFR per luminosity 0.68 times that of the old calibration, or a luminosity per SFR ${(0.68)}^{-1}=1.47$ times that of the old calibration. Any difference in this calibration due to using a Chabrier IMF is comparatively small, and we use the value from Murphy et al. (2011) without alteration.¹⁸

The convention in previous literature is to convert the above ${L}_{{\rm{H}}\alpha }/\mathrm{SFR}$ ratio into a ${L}_{\mathrm{Ly}\alpha }/\mathrm{SFR}$ ratio by assuming a Lyα/Hα ratio of 8.7. Common citations for this convention include

• 1.  
  Pengelly (1964) (communicated by Seaton), the first of three papers including Pengelly & Seaton (1964) and Seaton (1964);
• 2.  
  Brocklehurst (1971) (again communicated by Seaton, and in fact the content is very similar to Pengelly ( 1964));
• 3.  
  Hummer & Storey (1987);
• 4.  
  and Hu et al. (1998) (who cite Brocklehurst ( 1971), but are sometimes cited in isolation—for example, in Hayes (2015) and Bridge et al. (2018).

Of these, only the last citation is strictly appropriate, as it is the only one explicitly stating a Lyα/Hα ratio of 8.7. The first three deal with hydrogen and helium recombination rates and line ratios—including the Balmer series and specifically the Hα/Hβ ratio—for different possible gas densities and temperatures, and for different assumptions about whether the gas is optically thin (case A) or thick (case B) to the recombination lines. However, as Henry et al. (2015) note, there needs to be additional information to link the Lyα line to the Balmer series.

Osterbrock (1989) is a possible source for the Lyα/Hα ratio. In the low-density limit, 68% of recombinations lead to Lyα emission, and 45% lead to H-α emission; see Dijkstra et al. (2014) or Dijkstra (2017). Combining this with the ratio of photon energies, we obtain a Lyα/Hα flux ratio of 8.2.

However, this is in the low-density limit, and collisional excitations at higher densities result in an enhanced ratio. Typical assumptions for the electron density fall within the range of ${n}_{e}={10}^{2}$–10³ cm⁻³, and if we consult tables of line ratios as in Dopita & Sutherland (2003) (which Henry et al. (2015) consult for Hα/Hβ and Lyα/Hβ tables, and are synthesised from Storey & Hummer (1995)), we find that 8.2–9.1 is a reasonable range for the line ratio, given the density range at $T={10}^{4}$ K (and assuming case B recombination). The conventional value of 8.7 appears to have been chosen (or at least kept) as a happy intermediate.

The resulting intrinsic conversion between SFR and Lyα luminosity is

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{L}_{\mathrm{Ly}\alpha }}{\mathrm{SFR}} & = & \displaystyle \frac{{L}_{\mathrm{Ly}\alpha }}{{L}_{{\rm{H}}\alpha }}\displaystyle \frac{{\mathrm{SFR}}_{{\rm{K}}98}}{\mathrm{SFR}}\displaystyle \frac{{L}_{{\rm{H}}\alpha }}{{\mathrm{SFR}}_{{\rm{K}}98}}=\displaystyle \frac{8.7}{0.68}\displaystyle \frac{1.26\times {10}^{41}\,\mathrm{erg}\,{{\rm{s}}}^{-1}}{{M}_{\odot }\,{\mathrm{yr}}^{-1}}\\ & = & \displaystyle \frac{1.6\times {10}^{42}\,\mathrm{erg}\,{{\rm{s}}}^{-1}}{{M}_{\odot }\,{\mathrm{yr}}^{-1}}.\end{array}\end{eqnarray} \tag{ 24 }$$

This is the origin of our value for the numeric coefficient in Equation (4).

The above conversion operates under the assumption that recombination balances photoionization within H ii regions. However, this does not include the possibility of ionizing radiation being absorbed by dust before it is able to trigger a photoionization event, or the possibility of recombination line emission being absorbed by dust. There is also the possibility of ionizing photons escaping the galaxy without a photoionization event (in which case it will likely trigger an event in the intergalactic medium) or being absorbed in H i regions without triggering recombination line emission.

In this model, we ignore the last two possibilities for simplicity but model the first two, in an abbreviated version of the sort of model found in Cai et al. (2014):

$$\begin{eqnarray}&&{L}_{\mathrm{Ly}\alpha }(\mathrm{SFR},z)=C\mathrm{SFR}{f}_{\mathrm{esc}}^{\mathrm{ion}}{f}_{\mathrm{esc}}^{\mathrm{Ly}\alpha },\end{eqnarray} \tag{ 25 }$$

where C is the value obtained in Equation (24), and both escape fractions are functions of SFR and redshift. For this work, as shown in Equations (4) and (5), we lump the two escape fractions together into a single effective escape fraction relative to the intrinsic Lyα prediction.

We take ${f}_{\mathrm{esc}}^{\mathrm{ion}}\sim {f}_{\mathrm{esc}}^{\mathrm{Ly}\alpha }$, effectively squaring one escape fraction of UV photons against dust. Note that the escape/absorption mechanisms for ionizing photons (${\lambda }_{\mathrm{rest}}\leqslant 912$ Å) and Lyα photons (${\lambda }_{\mathrm{rest}}=1216$ Å) are in fact different, and any correlation between the two is subject to large amounts of scatter. The escape fractions are at least around the same order of magnitude, however—see the numbers obtained through simulations in Yajima et al. (2014) and the dust attenuation factors assumed in Cai et al. (2014) (the dust optical depth at 1216 Å is assumed to be 1.08 times the dust optical depth at 1350 Å, which itself is assumed to be $1/\gamma \simeq 1.18$ times the dust optical depth of ionizing photons).

The escape fraction is highly contrived, to two ends:

• 1.  
  It increases monotonically with redshift (converging to 1 as $z\to \infty$).
• 2.  
  It decreases with higher SFR.

The latter is easier to justify—it is natural to associate higher SFR with more dust, and there is observational evidence for this correlation (Santini et al. 2014). However, this in turn makes the former more difficult to justify: cosmic SFR density evolves nonmonotonically with redshift—increasing up to $z\sim 3$ before showing a clear decline after $z\sim 2$ (Madau & Dickinson 2014)—suggesting that the redshift evolution of the escape fraction cannot be monotonic. However, (a) the scope of our modeling is limited to $z\gtrsim 2$, where the behavior may as well be monotonic, and (b) it may be possible for factors other than SFR to influence dust content in older/late-type galaxies, although there is considerable uncertainty around the latter.

Assigning specific numbers requires either a sophisticated simulation incorporating radiative transfer (like Yajima et al. 2014) or observational constraints. Two simulation forecasts influence our choice of form and approximate parameter values, with subsequent fine-tuning based on observational constraints:

• 1.  
  Results in Yajima et al. (2014) show a median ${f}_{\mathrm{esc}}^{\mathrm{ion}}\sim 0.2$ evolving very weakly with redshift and a nonmonotonic evolution of ${f}_{\mathrm{esc}}^{\mathrm{Ly}\alpha }$, with median values ranging between ∼0.3 and ∼0.9. Scatter around the median is quite large, however.
• 2.  
  Results in Garel et al. (2012) show a simulated escape fraction of near-unity for most simulated galaxies with SFR less than 1 ${M}_{\odot }$ yr⁻¹. The distribution of ${f}_{\mathrm{esc}}(\mathrm{SFR})$ evolves strongly into a flat one with higher SFR, with an average of 21% for simulated galaxies with SFR greater than 10 ${M}_{\odot }$ yr⁻¹.

We model the total escape fraction ${f}_{\mathrm{esc}}\equiv {f}_{\mathrm{esc}}^{\mathrm{ion}}{f}_{\mathrm{esc}}^{\mathrm{Ly}\alpha }$ as the squared product of a generalized logistic function in redshift and an algebraic function in SFR:

$$\begin{eqnarray}\begin{array}{rcl}{f}_{\mathrm{esc}}(\mathrm{SFR},z) & = & \left[{\left(1+{e}^{-\xi (z-{z}_{0})}\right)}^{-\zeta }\right.\\ & & \times {\left.\left({f}_{0}+\displaystyle \frac{1-{f}_{0}}{1+{\left(\mathrm{SFR}/{\mathrm{SFR}}_{0}\right)}^{\varsigma }}\right)\right]}^{2}.\end{array}\end{eqnarray} \tag{ 26 }$$

This form combines S-shaped curves in each variable, and shows the desired asymptotic behavior discussed above. The function in redshift is an overall normalization between 0 and 1, with a characteristic redshift z₀ acting as an inflection point of the redshift evolution, and ξ and ζ controlling the shape. Meanwhile, the function in SFR changes between f₀ and 1 around a characteristic ${\mathrm{SFR}}_{0}$, with ς again controlling the shape.

With the model above (including 0.3 dex log-normal scatter in SFR and in Lyα luminosity), we are able to translate the analytic halo mass function fit from Behroozi et al. (2013a, 2013b) into simulated LFs at different redshifts. We use a brute-force technique, randomly drawing from the halo mass function and applying the Lyα model to the masses drawn, then binning the resulting luminosities to obtain an LF.

We tune the escape fraction parameter values based on comparing these simulated data to observed LFs at four different redshifts: $z\sim 0.3$ from Cowie et al. (2010), $z\sim 0.92$ from Barger et al. (2012), $z\sim 2.23$ from Sobral et al. (2017), and $z\sim 3.1$ from Gronwall et al. (2007). The resulting parameters are $\xi =1.6$, z₀ = 3.125, $\zeta =1/4$, f₀ = 0.18, ${\mathrm{SFR}}_{0}\approx 1.29\,{M}_{\odot }$ yr⁻¹, and $\varsigma =0.875$. We match the higher-redshift data better than the lower-redshift data, suggesting that our model cannot completely describe the strong evolution of the LAE LF from $z\sim 0.3$ to $z\sim 2$ (which is expected, given the monotonic redshift evolution of the escape fraction, as discussed in the previous section). However, note also that there is support for a composite Schechter/power-law LF for low-to-intermediate redshift LAEs, and by and large we are trying to match only the Schechter part of this. We show a comparison of the simulated LFs to these observed data in Figure 10.

Zoom In Zoom Out Reset image size

Without further tuning, we also compare against LFs derived in Sobral et al. (2018) from a compilation of deep and wide LAE surveys (dubbed S-SC4K), and the plots in Figure 11 show a reasonable match up to $z\sim 5$.

Zoom In Zoom Out Reset image size

Because the Chinchilla lightcones used in this work span z = 1.5–3.5, we can use these to simulate Lyα fluctuations and power spectra at $z\sim 2$ through the same methods used in the main work. We compare these in Figure 12 to power spectra in previous work in Pullen et al. (2014) (using only the halo contribution) and Fonseca et al. (2017), and find our model yields predicted power spectra squarely in between the two previous works.

Zoom In Zoom Out Reset image size

As was the case for Li et al. (2016), our choice to assign no line luminosities to halos below ${10}^{10}\,{M}_{\odot }$ in virial halo mass is partly pragmatic. From the point of view of simulation constraints, because we use a cosmological N-body box whose dark matter particle mass is only $5.9\times {10}^{8}{h}^{-1}\,{M}_{\odot }=8.4\times {10}^{8}\,{M}_{\odot }$, the halo population is severely incomplete for ${M}_{\mathrm{vir}}\lesssim {10}^{10}\,{M}_{\odot }$. Unlike CO, however, Lyα emission does not require a particularly dusty or high-metallicity environment, so any physical mass cutoff for Lyα emission would likely be much lower than for CO emission.

Note, however, that this cutoff mainly affects our simulations of HETDEX line-intensity cubes, because the Lyα luminosity cutoffs used for our mock LAE catalogs correspond typically to halo masses well above a ${10}^{10}\,{M}_{\odot }$ cutoff. Furthermore, even considering line-intensity cubes, our model ${L}_{{\rm{L}}y\alpha }({M}_{\mathrm{vir}})$ relation falls off quite sharply for ${M}_{\mathrm{vir}}\lesssim {10}^{11}\,{M}_{\odot }$. Therefore, even if we assigned model luminosities to a well-represented population of halos with ${M}_{\mathrm{vir}}\lesssim {10}^{10}\,{M}_{\odot }$, they would likely not contribute significantly to the signal.

We show this in Figure 13 via analytic calculations of contributions to average line intensity from different ranges of halo masses. We calculate the line luminosity per volume ${{dL}}_{\mathrm{line}}/{dV}=\int {L}_{\mathrm{line}}(M)({dn}/{dM}){dM}$, where dn/dM is the halo mass function fit in Behroozi et al. (2013b) at the appropriate redshift, and convert this into observer quantities of CO brightness temperature and Lyα $\nu {I}_{\nu }$. For both CO and Lyα emission, the contribution to mean line intensity falls off below ${10}^{11}\,{M}_{\odot }$ in halo mass for both lines, but especially rapidly for Lyα emission; the analytic results suggest that we have captured a great majority of any expected signal using our cutoff halo mass of ${10}^{10}\,{M}_{\odot }$, at least for our assumed model.

Zoom In Zoom Out Reset image size

Note in particular that, even with the same 3D $P({\boldsymbol{k}})$, there is no reason to expect the total signal-to-noise to be the same for P(k) and P(k, μ). For simplicity, take the CO auto spectrum (ignoring or absorbing ${W}^{2}(k)$) as an example. Suppose that we have n_μ μ-bins in each k-shell, such that

$$\begin{eqnarray}&&{N}_{\mathrm{modes}}(k)=\displaystyle \sum _{i}^{{n}_{\mu }}{N}_{m,i}(k)\end{eqnarray} \tag{ 40 }$$

and ${N}_{m,i}(k)$ represents the number of modes in the μ-bin centered at ${\mu }_{i}$ falling within the k-shell centered at k. (From here on, summing over i always implies a sum over all n_μ applicable values.)

In the main text, we calculate the spherically averaged spectrum P(k) and then the signal-to-noise ratio at each k (which is then added in quadrature over all k to obtain total S/N across all modes):

$$\begin{eqnarray}&&{\left[\displaystyle \frac{{\rm{S}}}{{\rm{N}}}(k)\right]}_{\mathrm{sph}}^{2}=\displaystyle \frac{P{\left(k\right)}^{2}{N}_{\mathrm{modes}}(k)}{{\left[{P}_{n}+P(k)\right]}^{2}},\end{eqnarray} \tag{ 41 }$$

where $P(k)=[{\sum }_{i}P(k,{\mu }_{i}){N}_{m,i}(k)]/{N}_{\mathrm{modes}}(k)$. We substitute and simplify to obtain

$$\begin{eqnarray}&&{\left[\displaystyle \frac{{\rm{S}}}{{\rm{N}}}(k)\right]}_{\mathrm{sph}}^{2}=\displaystyle \frac{{\left[{\displaystyle \sum }_{i}P(k,{\mu }_{i}){N}_{m,i}(k)\right]}^{2}/{N}_{\mathrm{modes}}(k)}{{\left[{P}_{n}+{\displaystyle \sum }_{i}P(k,{\mu }_{i}){N}_{m,i}(k)/{N}_{\mathrm{modes}}(k)\right]}^{2}}.\end{eqnarray} \tag{ 42 }$$

However, if we took the S/N for each ${\mu }_{i}$ and then averaged, we would have

$$\begin{eqnarray}&&{\left[\displaystyle \frac{{\rm{S}}}{{\rm{N}}}(k)\right]}_{\mathrm{aniso}}^{2}=\displaystyle \sum _{i}{\left[\displaystyle \frac{{\rm{S}}}{{\rm{N}}}(k,{\mu }_{i})\right]}^{2}=\displaystyle \sum _{i}\displaystyle \frac{P{\left(k,{\mu }_{i}\right)}^{2}{N}_{m,i}(k)}{{\left[{P}_{n}+P(k,{\mu }_{i})\right]}^{2}}.\end{eqnarray} \tag{ 43 }$$

It is difficult to see a way that the two can be generally equal. It is reasonable to approximate ${N}_{m,i}(k)\approx {N}_{\mathrm{modes}}(k)/{n}_{\mu }$ (intervals in μ are roughly linear with intervals in the polar angle up to $\mu \sim 0.5$, so the number of modes in each μ-bin should be similar), in which case

$$\begin{eqnarray}&&{\left[\displaystyle \frac{{\rm{S}}}{{\rm{N}}}(k)\right]}_{\mathrm{sph}}^{2}=\displaystyle \frac{{\left[{\displaystyle \sum }_{i}P(k,{\mu }_{i})\right]}^{2}{N}_{\mathrm{modes}}(k)/{n}_{\mu }^{2}}{{\left[{P}_{n}+{\displaystyle \sum }_{i}P(k,{\mu }_{i})/{n}_{\mu }\right]}^{2}},\end{eqnarray} \tag{ 44 }$$

and

$$\begin{eqnarray}&&{\left[\displaystyle \frac{{\rm{S}}}{{\rm{N}}}(k)\right]}_{\mathrm{aniso}}^{2}=\displaystyle \sum _{i}\displaystyle \frac{P{\left(k,{\mu }_{i}\right)}^{2}{N}_{\mathrm{modes}}(k)/{n}_{\mu }}{{\left[{P}_{n}+P(k,{\mu }_{i})\right]}^{2}}.\end{eqnarray} \tag{ 45 }$$

In the case that sample variance dominates our uncertainties, i.e., ${P}_{n}\ll P(k,{\mu }_{i})$,

$$\begin{eqnarray}&&{\left[\displaystyle \frac{{\rm{S}}}{{\rm{N}}}(k)\right]}_{\mathrm{aniso}}^{2}=\displaystyle \sum _{i}\displaystyle \frac{{N}_{\mathrm{modes}}(k)}{{n}_{\mu }}={N}_{\mathrm{modes}}(k)={\left[\displaystyle \frac{{\rm{S}}}{{\rm{N}}}(k)\right]}_{\mathrm{sph}}^{2}.\end{eqnarray} \tag{ 46 }$$

However, we now demonstrate that, in the extreme case where instrumental noise dominates the uncertainties and only one of n_μ bins has a nonzero value of $P(k,{\mu }_{i})$, the signal-to-noise ratio is significantly higher for P(k,μ) than for P(k). (Note that it would be misleading to assume ${P}_{n}\ll P(k,{\mu }_{i})$ for all i when the right-hand side is zero for most i.) Take ${\mu }_{0}$ to be the bin with nonzero P(k,μ):

$$\begin{eqnarray}&&{\left[\displaystyle \frac{{\rm{S}}}{{\rm{N}}}(k)\right]}_{\mathrm{sph}}^{2}=\displaystyle \frac{P{\left(k,{\mu }_{0}\right)}^{2}{N}_{\mathrm{modes}}(k)/{n}_{\mu }^{2}}{{\left[{P}_{n}+P(k,{\mu }_{0})/{n}_{\mu }\right]}^{2}},\end{eqnarray} \tag{ 47 }$$

and

$$\begin{eqnarray}&&{\left[\displaystyle \frac{{\rm{S}}}{{\rm{N}}}(k)\right]}_{\mathrm{aniso}}^{2}=\displaystyle \frac{P{\left(k,{\mu }_{0}\right)}^{2}{N}_{\mathrm{modes}}(k)/{n}_{\mu }}{{\left[{P}_{n}+P(k,{\mu }_{0})\right]}^{2}}.\end{eqnarray} \tag{ 48 }$$

If ${P}_{n}\gg P(k,{\mu }_{0})$, as is typical in our simulated surveys,

$$\begin{eqnarray}&&{\left[\displaystyle \frac{{\rm{S}}}{{\rm{N}}}(k)\right]}_{\mathrm{sph}}^{2}=\displaystyle \frac{P{\left(k,{\mu }_{0}]\right)}^{2}{N}_{\mathrm{modes}}(k)}{{n}_{\mu }^{2}{P}_{n}}=\displaystyle \frac{1}{{n}_{\mu }}{\left[\displaystyle \frac{{\rm{S}}}{{\rm{N}}}(k)\right]}_{\mathrm{aniso}}^{2},\end{eqnarray} \tag{ 49 }$$

resulting in a difference of a factor of $\sqrt{{n}_{\mu }}$ in signal-to-noise at this k.

For one lightcone out of our 100, we simulate cross-correlations between a CO cube and galaxy sample with ${M}_{* }\gt {10}^{10}\,{M}_{\odot }$ for all ${\sigma }_{z}/(1+z)$ values considered in the main text. However, we now obtain the anisotropic P(k,μ), which we show in Figure 15. As anticipated, the signal loss with increasing ${\sigma }_{z}/(1+z)$ is highly anisotropic, and is greater for higher μ at any given k.

Zoom In Zoom Out Reset image size

Table 6 also shows the total signal-to-noise ratios across all modes for this realization, comparing between P(k) and P(k,μ). The difference is notable for high values of ${\sigma }_{z}/(1+z)$, but not as great as the $\sqrt{{n}_{\mu }}$ in the simple calculation above (which would have been a factor of almost 8). Furthermore, as we show in Figure 16, the change in signal-to-noise with frequency resolution differs from what we show in Figure 5 for P(k). Working with P(k,μ), by definition, separates the line-of-sight modes from the transverse modes, and thus the different degrees of attenuation experienced due to redshift error. Thus, the only effect of decreasing the number of frequency channels in the survey volume is to decrease the number of modes averaged and consequently to increase uncertainties; the slight gain in P(k) signal-to-noise shown in Figure 5 for ${\sigma }_{z}/(1+z)\geqslant 0.01$ is absent in the P(k,μ) signal-to-noise curves in Figure 16. We thus consider our basic conclusion—that photometric errors significantly reduce any advantage in detection significance from cross-correlation—to be unchanged.

Zoom In Zoom Out Reset image size

Table 6.  Total Signal-to-noise Ratio over All Modes for Spherically Averaged Power Spectra (${\rm{S}}/{{\rm{N}}}_{\mathrm{sph}}$) and Anisotropic Power Spectra (${\rm{S}}/{{\rm{N}}}_{\mathrm{aniso}}$) in One Realization

Power Spectrum	${\sigma }_{z}/(1+z)$ for galaxies	${\rm{S}}/{{\rm{N}}}_{\mathrm{sph}}$	${\rm{S}}/{{\rm{N}}}_{\mathrm{aniso}}$
CO-galaxy cross	$0.$	33.2	37.7
CO-galaxy cross	0.0007	29.9	28.6
CO-galaxy cross	0.003	19.9	19.9
CO-galaxy cross	0.01	11.3	13.3
CO-galaxy cross	0.02	7.7	10.6
CO-galaxy cross	0.03	6.6	9.5
CO auto	...	4.7	5.0

Note. For the galaxy sample, $\mathrm{log}({M}_{* ,\min }/{M}_{\odot })=10$, with a galaxy count of $2.9\times {10}^{4}$ without redshift errors. All signal-to-noise ratios are still quoted for a single patch of 2.5 deg² observed for 1500 hr; we may expect to see this increase by up to a factor of $\sqrt{2}$ if two equivalent patches are observed for 1500 hr each, and a further, roughly linear increase with more integration time.

Download table as:  ASCIITypeset image