COMAP Early Science. VII. Prospects for CO Intensity Mapping at Reionization

The Epoch of Reionization (EoR) remains one of the least-explored periods of cosmic history. During this final cosmic phase transition, photons from the first luminous sources ionized the intergalactic medium (IGM) for the first time since the emission of the cosmic microwave background (CMB; Loeb & Barkana 2001; McQuinn 2016). Observations of the optical depth to the CMB have placed integrated limits on the redshift and duration of reionization (Planck Collaboration et al. 2020), but the details of the process are still largely unobserved.

The EoR has long been a popular target for line intensity mapping (LIM) survey proposals (Kovetz et al. 2017). Because they are sensitive to the aggregate emission from all emitting objects at a given redshift, intensity maps are less limited by the extreme faintness of the individual redshift z ≳ 6 sources that reionized the universe. By mapping line emission at different observing frequencies and therefore different redshifts, it is in principle possible to map the three-dimensional structure of the universe as reionization proceeds.

The first target for LIM surveys at the EoR was the 21 cm hyperfine transition from neutral hydrogen (Pritchard & Loeb 2012). Several experiments have sought or are seeking to use the 21 cm line to map the gradual disappearance of the neutral IGM across the EoR (Ali et al. 2015; Beardsley et al. 2016; DeBoer et al. 2017). The 21 cm line, however, carries little sensitivity to the actual ionizing sources themselves. Models of star formation and the IGM during the EoR are highly sensitive to assumptions about the interstellar media of star-forming galaxies and extrapolations of the luminosity function of faint, undetectable sources (McQuinn 2016). These properties could be probed by an intensity mapping survey specifically focused on the ionizing sources rather than the IGM.

Intensity mapping using rotational transitions of carbon monoxide (CO), first discussed in Righi et al. (2008) as a possible CMB foreground, traces aggregate emission from the dense molecular gas in which most star formation occurs. Lidz et al. (2011) demonstrated that such an observation would be particularly powerful during reionization as a complement to 21 cm surveys by probing the formation of new stars thought to provide the bulk of ionizing photons. A sufficiently high-redshift CO intensity mapping survey would be uniquely able to measure the total abundance of molecular gas during reionization. Moreover, because LIM measurements are sensitive to the faint end of the galaxy luminosity function, they can determine which galaxies contribute most to that measurement, whether reionization is dominated by rare bright objects or numerous fainter ones.

Realizing this potential for EoR intensity mapping of CO is a major goal of the Carbon Monoxide Mapping Array Project (COMAP). As described in prior papers in this series (Cleary et al. 2022), the currently observing COMAP Pathfinder is pursuing CO intensity mapping over three 4 deg² fields in a frequency band centered at 30 GHz using a 10.4 m antenna at the Owens Valley Radio Observatory (OVRO). The CO signal in this band is expected to be dominated by CO(1–0) emission from redshifts z = 2.4–3.4, the exploration of which makes up the primary science goals of the Pathfinder (Li et al. 2016; Chung et al. 2022). These observing frequencies also contain subdominant emission from the CO(2–1) line emitted at z ∼ 6–8, spanning a large portion of the EoR. We introduce here an experimental concept designed to isolate the CO intensity mapping signal from reionization and produce high-quality maps of star-forming molecular gas during that epoch.

Our planned extension of the existing COMAP Pathfinder, which we term COMAP-EoR, will accomplish this in two ways. First, more 30 GHz detectors will increase sensitivity to CO(1–0) from z ∼ 3 and CO(2–1) from z ∼ 6–8 (during the galaxy assembly and reionization epochs, respectively). Second, we will add a second frequency band centered at 16 GHz which will have access to the CO(1–0) transition at the EoR. We can autocorrelate the 16 GHz maps to measure the EoR signal directly. We can also cross correlate the 16 and 30 GHz bands to isolate the EoR CO(2–1) signal from the dominant lower-redshift line. This will also have the benefit of minimizing any other foreground or systematic effects present in either band.

In Chung et al. (2022), we provided a state-of-the art phenomenological model of CO emission during the z ∼ 3 galaxy assembly epoch, which is the focus of the COMAP Pathfinder. As we will see in this paper, our understanding of reionization remains far too limited to make a similar attempt at z ∼ 7. The paucity of directly detected CO emitters at these redshifts means that we cannot fit a useful empirical luminosity function, so we are left with scaling relations and ISM models that are even less certain than they are at z ∼ 3. Other effects like CMB heating and metallicity evolution further complicate a first-principles modeling effort. We will instead adopt a method similar to that of Breysse et al. (2014), and examine the space of existing literature models for the CO signal, using the range of signal amplitudes as a proxy for the detailed uncertainty calculations presented earlier in this series. We model the auto-power spectra of CO(1–0) and CO(2–1) at reionization, the galaxy assembly era CO(1–0) signal that dominates the 30 GHz band, and the cross spectrum between the two bands. For all but one of the literature models we predict a highly significant detection of EoR CO emission using the COMAP-EoR design. We go on to show that this measurement can place tight constraints on the abundance of high-redshift molecular gas, and the population of galaxies below the detection threshold of conventional surveys. We also show a dramatic improvement of the z ∼ 3 CO measurements compared to the Pathfinder forecasts, enabling extremely precise study of molecular gas and star formation in this later epoch.

We also include forecasts for a hypothetical third-generation Expanded Reionization Array stage of COMAP, termed COMAP-ERA, intended to follow the Pathfinder and COMAP-EoR surveys. This survey further increases the sensitivity at both the 16 and 30 GHz bands. We discuss how this extra depth can allow us fulfill the promise of tomographic 21 cm and CO cross correlation, tracing the coevolution of the interstellar and intergalactic media during the EoR.

For one of our literature models, based on semianalytic simulations presented in Yang et al. (2021), we find a CO signal that is considerably fainter than all of the others, is effectively undetected in COMAP-EoR, and is only seen at the lowest redshifts by COMAP-ERA. We use this model to demonstrate that the CO LIM observations proposed here remain scientifically useful even if only upper limits are obtained. At these sensitivities, our z ∼ 3 CO measurements are quite sensitive to contamination from EoR emission, so directly placing an EoR limit is necessary to reach the potential of the z ∼ 3 measurements, even if no EoR signal is detected. In addition, a LIM upper limit at z ∼ 6–8 in combination with a direct-detection CO survey would serve as proof that there is no significant unseen reservoir of molecular gas during reionization beyond that which is directly imaged.

Section 2 outlines the experimental design of COMAP-EoR, and we describe our power spectrum formalism in Section 3. Section 4 summarizes the literature models we use for our forecasts, the results of which appear in Section 5. We explore the scientific implications of a COMAP-EoR detection in Section 6. Further discussion appears in Section 7, and we conclude in Section 8. Throughout this work, we assume a flat Λ cold dark matter cosmology (CDM) consistent with the Planck 2018 results (Planck Collaboration et al. 2020). More detail on the COMAP Pathfinder can be found in the other papers in this series, including discussions of the instrumental hardware (Lamb et al. 2022), the data reduction pipeline (Foss et al. 2022), the power spectrum analysis (Ihle et al. 2022), the science and modeling implications (Chung et al. 2022), and the auxiliary Galactic plane observations (Rennie et al. 2022).

As with other papers in this series, we assume a flat ΛCDM cosmology with parameters Ω_m = 0.286, Ω_Λ = 0.714, Ω_b = 0.047, H₀ = 100h km s⁻¹ Mpc⁻¹ with h = 0.7, σ₈ =0.82, and n_s = 0.96, to maintain consistency with previous COMAP simulations (Li et al. 2016; Ihle et al. 2019). The cosmology is also broadly consistent with 9 yr Wilkinson Microwave Anisotropy Probe results (Hinshaw et al. 2013). Distances carry an implicit h⁻¹ dependence throughout, which propagates through masses (all based on virial halo masses, proportional to h⁻¹) and volume densities (∝h³).

Here we will outline the design of the COMAP-EoR instrument and survey. Continued observations with the current 30 GHz Pathfinder instrument will be supplemented with two additional 30 GHz receivers mounted on existing 10.4 m antennas at OVRO along with a new 16 GHz receiver mounted on a new 18 m antenna designed as a prototype for the next-generation Very Large Array (ngVLA). The basic parameters of these instruments can be found in Table 1. The predicted system temperature for the 30 GHz observations is based on the existing Pathfinder, while that for the 16 GHz instrument is based on expectations for ngVLA given in Selina et al. (2018).

Table 1. Parameters of the COMAP-EoR Instruments, Including the Central Frequency ν_obs and Width Δν of Each Band, the Redshift Ranges z₍₁₋₀₎ and z₍₂₋₁₎ of Each Line, the System Temperature T_sys, the Number of Feeds N_feeds per Dish, the Beam FWHM θ_FWHM, and the Channel Width δ ν

Band	Frequency	Bandwidth	CO(1–0) Range	CO(2–1) Range	System Temp.	Feeds/Dish	Beam FWHM	Channel Width
	νobs	Δν	z(1−0)	z(2−1)	Tsys	Nfeeds	θFWHM	δ ν
1	12.5 GHz	1 GHz	⋯	7.8–8.6	20 K	38 a	42	2 MHz
2	14 GHz	2 GHz	⋯	6.7–7.8	20 K	38 a	40	2 MHz
3	16 GHz	2 GHz	⋯	5.8–6.7	22 K	38 a	37	2 MHz
4	18.5 GHz	3 GHz	⋯	4.8–5.8	27 K	38 a	33	2 MHz
2	28 GHz	4 GHz	6.7–7.8	2.8–3.4	44 K b	19	45	2 MHz
3	32 GHz	4 GHz	5.8–6.7	2.4–2.8	44 K b	19	39	2 MHz

Notes.

^a 19 dual-polarization feeds. ^b The current Pathfinder has a system temperature of 44 K; we expect additional instruments to have an improved value of 34 K. This will be accounted for below in our effective observing time.

Download table as:  ASCIITypeset image

For the early Pathfinder observations reported in other papers in this series, we have assumed a measurement averaged over the entire 8 GHz wide band of the 30 GHz instrument. Because the Pathfinder is focused on galaxy assembly era measurements near the peak of cosmic star formation, we do not expect the z ∼ 3 CO signal to evolve dramatically over this frequency range, at least at the relatively low sensitivity of the Pathfinder. For COMAP-EoR, however, we have much more sensitivity to work with, and we have a reionization-era signal that may evolve quite dramatically over our frequency range. Thus, in Table 1 and throughout this paper we have divided our observations into several redshift bins. It is challenging at this point to determine an optimal redshift binning for the COMAP-EoR measurement, as both the overall amplitude and the redshift dependence of the CO signal are highly uncertain. We will make the somewhat arbitrary choice here to use a total of four frequency bands: two centered at 28 and 32 GHz in the high-frequency instrument and 14 and 16 GHz in the low-frequency that span the overlapping volume between the CO(1–0) and (2–1) lines at EoR (see Figure 1), and two other bands centered at 12.5 and 18.5 GHz that account for the additional frequency coverage of the low-frequency instrument. We leave for future work a more detailed calculation of the optimal redshift binning.

Zoom In Zoom Out Reset image size

Given the availability of improved low-noise amplifiers since the deployment of the Pathfinder, we expect the new 30 GHz instruments to have a somewhat lower system temperature, T_sys = 34 K compared to the Pathfinder's 44 K. Below, for ease of forecasting, we will assign T_sys = 44 K to all three 30 GHz receivers, and we will scale the effective observing time to account for the improved sensitivity. The noise level σ_N in a map scales as

$$\begin{eqnarray}&&{\sigma }_{N}\propto \displaystyle \frac{{T}_{\mathrm{sys}}}{\sqrt{{t}_{\mathrm{obs}}}},\end{eqnarray} \tag{ 1 }$$

for integration time t_obs, so we can write our total effective observing time as

$$\begin{eqnarray}&&{t}_{\mathrm{obs}}^{\mathrm{eff}}={t}_{\mathrm{obs}}^{\mathrm{PF}}+2{\left(\displaystyle \frac{44\ {\rm{K}}}{34\ {\rm{K}}}\right)}^{2}{t}_{\mathrm{obs}}^{\mathrm{new}},\end{eqnarray} \tag{ 2 }$$

where ${t}_{\mathrm{obs}}^{\mathrm{PF}}$ is the total observing time with the 44 K Pathfinder instrument, and ${t}_{\mathrm{obs}}^{\mathrm{new}}$ is the observing time on each of the new dishes.

For our forecasts here, we will assume that the 16 GHz instrument comes online immediately following the end of the current 5 yr Pathfinder campaign, with the two new 30 GHz receivers following an additional 2 yr after that. The nominal COMAP-EoR campaign will then consist of five more years of operation with the full four instruments. Under this time line, at the end of those 5 yr we will have accumulated 12 yr of time on the Pathfinder, 7 yr on the 16 GHz dish, and 5 yr on each of the new 30 GHz dishes. For the COMAP-EoR survey, we plan to continue to target the same three fields as the current Pathfinder observations. Assuming 1000 hr per year per field of available time, this gives us a total of 29,000 dish hours per field at 30 GHz, accounting for the T_sys adjustment, and 7000 dish hours per field at 16 GHz.

We will also provide forecasts for a hypothetical third-generation COMAP-ERA to provide an idea of what could be accomplished with even further increases in sensitivity. As this concept is relatively far in the future, we will model it fairly simply here. We will assume that, at the end of the above COMAP-EoR survey, we increase to 10 dishes at each frequency and observe for an additional 5 yr with all 20. Assuming these new instruments are identical to the COMAP-EoR equivalents, this would give us 110,000 dish hours at 30 GHz and 57,000 dish hours at 16 GHz. For simplicity, here we will assume all of this time is spent on continued observations of the same three Pathfinder fields. In practice, we may consider other fields in order to cross correlate with other EoR data, a possibility we will briefly discuss in Section 7.

Because COMAP-EoR will observe at 16 and 30 GHz simultaneously, our task in modeling the power spectrum is necessarily more complicated than for the Pathfinder. We will need to model the reionization-era auto spectra of the CO(1–0) and (2–1) lines as well as their cross correlation. In addition, we require a model of the galaxy assembly era CO(1–0) auto spectrum. Though this lower-redshift line is the primary signal for the current Pathfinder observations, it serves as an important foreground to the COMAP-EoR CO(2–1) measurement. Interloper power spectra in intensity mapping surveys become distorted anisotropically when projected into the frame of a higher-redshift signal (Visbal & Loeb 2010; Cheng et al. 2016; Lidz & Taylor 2016), so we will need to model the full angular behavior of the power spectra. Our formalism here is primarily based on Bernal et al. (2019), reproduced here for the convenience of the reader.

Throughout this paper, we will make the simplifying assumption that the CO emission in a given dark matter halo comes from a single point source at its center, i.e., we neglect one-halo contributions. We will also assume that an intensity mapping signal does not evolve across a given frequency band, and will calculate all of our power spectra at a redshift corresponding to the band center. We also neglect for now the line-broadening effects discussed in Chung et al. (2021). As we will see below, we expect any inaccuracies arising due to these assumptions to be small compared to the overall modeling uncertainty.

3.1. Auto Spectra

The analytic form for an intensity mapping auto spectrum P(k, μ, z) can be written as

$$\begin{eqnarray}&&P(k,\mu ,z)={P}^{\mathrm{clust}}(k,\mu ,z)+{P}^{\mathrm{shot}}(z),\end{eqnarray} \tag{ 3 }$$

where k is the magnitude of the wavevector of a given Fourier mode and μ is the cosine of the angle between that mode and the line of sight.

On large scales, the power spectrum is dominated by clustered emission, which traces the large-scale structure that takes the form

$$\begin{eqnarray}&&{P}^{\mathrm{clust}}={\left\langle {Tb}\right\rangle }^{2}(z){F}_{\mathrm{RSD}}^{2}(k,\mu ,z){P}_{m}(k,z).\end{eqnarray} \tag{ 4 }$$

The shape of the power spectrum is set by the dark matter power spectrum P_m , computed here using CAMB (Lewis & Bridle 2002). The overall amplitude of the intensity mapping power is set by the luminosity-weighted bias of the target line emitters, given by

$$\begin{eqnarray}&&\left\langle {Tb}\right\rangle ={C}_{\mathrm{LT}}(z){\int }_{M\min }^{\infty }L(M,z)b(M,z){f}_{\mathrm{duty}}(M,z)\displaystyle \frac{{dn}}{{dM}}{dM},\end{eqnarray} \tag{ 5 }$$

where L(M, z) is the mean CO luminosity of a halo with mass M, f_duty(M, z) is the fraction of halos with mass M, which emit CO at any given time, b(M, z) is the bias for a halo of mass M (Tinker et al. 2010), and dn/dM is the halo mass function, for which we assume the form of Tinker et al. (2008). The clustering amplitude $\left\langle {Tb}\right\rangle$ is often expressed as a product of the mean line intensity $\left\langle T\right\rangle$, here expressed in brightness temperature units, and the bias b of the emitting galaxies, which is defined below. In order to ensure that the mass function integral converges, we assume that only halos with masses above ${M}_{\min }$ emit CO, with the value of ${M}_{\min }$ set by the model under consideration. The factor

$$\begin{eqnarray}&&{C}_{\mathrm{LT}}(z)=\displaystyle \frac{{c}^{3}(1+z)}{8\pi {k}_{{\rm{B}}}{\nu }^{3}H(z)}\end{eqnarray} \tag{ 6 }$$

is the conversion factor between luminosity density and brightness temperature, where c is the speed of light, k_B is Boltzmann's constant, H(z) is the Hubble parameter, and ν is the rest frequency of the target line. We can also separate out an average bias factor

$$\begin{eqnarray}&&b(z)=\displaystyle \frac{\left\langle {Tb}\right\rangle }{\left\langle T\right\rangle }=\displaystyle \frac{{\int }_{M\min }^{\infty }L(M,z)b(M,z){dn}/{dMdM}}{{\int }_{M\min }L(M,z){dn}/{dMdM}},\end{eqnarray} \tag{ 7 }$$

which gives the degree to which the galaxies are more strongly clustered than the underlying dark matter.

Because intensity maps are made in redshift space, the observed intensity field is distorted anisotropically compared to the true field (see Hamilton 1998, for a review). The effect of these distortions on the power spectrum is encoded in F_RSD, which is given by

$$\begin{eqnarray}&&{F}_{\mathrm{RSD}}(k,\mu ,z)=\left(1+\displaystyle \frac{f(z)}{b(z)}{\mu }^{2}\right)\displaystyle \frac{1}{1+{\left(k\mu {\sigma }_{\mathrm{FoG}}\right)}^{2}/2}.\end{eqnarray} \tag{ 8 }$$

The first term gives the linear Kaiser effect (Kaiser 1987), which dominates on large scales with amplitude set by the logarithmic derivative of the growth factor f(z). The second term describes the fingers-of-God effect due to small-scale peculiar velocities, for which we assume a Lorentzian form with σ_FoG = 7 Mpc (Bernal et al. 2019).

Because our CO emission is sourced by a population of discrete galaxies rather than a continuous background, our power spectrum also includes a scale-independent shot noise component given by

$$\begin{eqnarray}\begin{array}{rcl}{P}^{\mathrm{shot}}(z) & = & {C}_{\mathrm{LT}}^{2}(z){\displaystyle \int }_{M\min }^{\infty }{L}^{2}(M,z){f}_{\mathrm{duty}}(M,z)\\ & & \times \,\displaystyle \frac{{dn}}{{dM}}{e}^{{\left[{\sigma }_{\mathrm{sc}}(M,z)\mathrm{ln}(10)\right]}^{2}}{dM}.\end{array}\end{eqnarray} \tag{ 9 }$$

The exponential factor comes from the fact that some models account for a scatter in line luminosity among halos with a given mass, rather than assigning all halos the mean L(M, z) (Li et al. 2016). We make the common assumption that this scatter has a lognormal form with standard deviation σ_sc(M, z) in units of dex (Li et al. 2016), which yields the expression above (see the discussion in Breysse et al. 2022).

3.2. Cross Spectrum

For COMAP-EoR, we also need to model the cross correlation between reionization-era CO(1–0) and (2–1). We can again write the cross-power spectrum P_× between two intensity mapping lines in the form

$$\begin{eqnarray}&&{P}_{\times }(k,\mu ,z)={P}_{\times }^{\mathrm{clust}}(k,\mu ,z)+{P}_{\times }^{\mathrm{shot}}(z).\end{eqnarray} \tag{ 10 }$$

The basic forms of the clustering and shot noise terms are derived in the appendix of Liu & Breysse (2021). The cross-clustering term is

$$\begin{eqnarray}\begin{array}{rcl}{P}_{\times }^{\mathrm{clust}} & = & {\left\langle {Tb}\right\rangle }_{1}(z){\left\langle {Tb}\right\rangle }_{2}(z)\\ & & \times \,{F}_{\mathrm{RSD},1}(k,\mu ,z){F}_{\mathrm{RSD},2}(k,\mu ,z){P}_{m}(k,z),\end{array}\end{eqnarray} \tag{ 11 }$$

where subscripts 1 and 2 indicate values computed for CO(1–0) and CO(2–1), respectively. Since both lines are sourced by the same discrete sources, we also have a cross-shot noise term

$$\begin{eqnarray}\begin{array}{rcl}{P}_{\times }^{\mathrm{shot}} & = & {C}_{\mathrm{LT},1}(z){C}_{\mathrm{LT},2}(z){e}^{{r}_{\mathrm{sc}}{\sigma }_{\mathrm{sc},1}{\sigma }_{\mathrm{sc},2}{\left[\mathrm{ln}(10)\right]}^{2}}\\ & & \times {\displaystyle \int }_{M\min }^{\infty }{L}_{1}(M,z){L}_{2}(M,z){f}_{\mathrm{duty}}(M,z)\displaystyle \frac{{dn}}{{dM}}{dM}.\end{array}\end{eqnarray} \tag{ 12 }$$

The exponential factor again comes from the increased shot noise due to the scatter about the mean L(M) relations. For the cross-shot case, we introduce the correlation coefficient r_sc between the scatter in the two lines, where r_sc = 1 corresponds to the perfectly correlated case, r_sc = 0 corresponds to the case where the scatter in the two lines is completely independent, and r_sc = −1 would correspond to the case where the scatter in the two lines is anticorrelated. For a derivation of this cross-scatter effect, see Appendix A. For the remainder of this work, we assume r_sc = 1 following Schaan & White (2021) and Yang et al. (2022).

3.3. CO(1–0) Interloper

The 30 GHz component of the COMAP-EoR survey will include a substantial contribution from galaxy assembly era CO(1–0) emitters. Formally, this z ∼ 3 emission constitutes a form of foreground contamination to the EoR signal, which will complicate attempts to measure EoR-era CO(2–1). The galaxy assembly contribution to the overall 30 GHz power spectrum can easily be an order of magnitude or more greater than the EoR contribution (see Figure 4 below). However, as seen in other papers in this series, this interloper provides significant scientific value in its own right. For our forecasts here, we will therefore simultaneously model the contributions from z ∼ 7 and z ∼ 3 so that our final molecular gas forecasts cover the full accessible redshift range.

In our 30 GHz data, low-redshift CO(1–0) and high-redshift CO(2–1) will be mapped into a common coordinate system. Thus, we must account for projection effects when dealing with these two lines. Because this paper is primarily (though not entirely) focused on reionization, here we will project the galaxy assembly line into the higher-redshift coordinate system. Thus, our actual observed power spectrum at 30 GHz will take the form

$$\begin{eqnarray}&&P(k,\mu )={P}_{2}^{\mathrm{EoR}}(k,\mu )+{P}_{1,\mathrm{proj}}^{\mathrm{GA}}(k,\mu ),\end{eqnarray} \tag{ 13 }$$

where the notations "EoR" and "GA" (for "galaxy assembly") denote z ∼ 7 and z ∼ 3 quantities, respectively. The apparent projected CO(1–0) power spectrum can be written as

$$\begin{eqnarray}&&{P}_{1,\mathrm{proj}}^{\mathrm{GA}}({k}_{\parallel },{k}_{\perp })=\displaystyle \frac{1}{{\alpha }_{\parallel }{\alpha }_{\perp }^{2}}{P}_{1}^{\mathrm{GA}}\left(\displaystyle \frac{{k}_{\parallel }}{{\alpha }_{\parallel }},\displaystyle \frac{{k}_{\perp }}{{\alpha }_{\perp }}\right),\end{eqnarray} \tag{ 14 }$$

where

$$\begin{eqnarray}&&{\alpha }_{\parallel }=\displaystyle \frac{H({z}_{\mathrm{EoR}})}{H({z}_{\mathrm{GA}})}\displaystyle \frac{1+{z}_{\mathrm{GA}}}{1+{z}_{\mathrm{EoR}}},\end{eqnarray} \tag{ 15 }$$

and

$$\begin{eqnarray}&&{\alpha }_{\perp }=\displaystyle \frac{{D}_{A}({z}_{\mathrm{GA}})}{{D}_{A}({z}_{\mathrm{EoR}})}\end{eqnarray} \tag{ 16 }$$

are the scalings for modes oriented parallel and perpendicular to the line of sight assuming the Hubble parameter H(z) and the comoving angular diameter distance D_A (z). The original, unprojected ${P}_{1}^{\mathrm{GA}}$ auto spectrum can be computed using an assumed L(M) model in the same manner as the high-redshift lines.

These projection effects mean that, even though we expect the high-redshift CO(2–1) line to be subdominant to the lower-redshift CO(1–0), we can still hope to obtain some information about the CO(2–1) auto spectrum. The projection adds extra anisotropy to the projected power spectrum, which can be used to separate the two signals (Cheng et al. 2016; Lidz & Taylor 2016). We leave for future work a discussion of other methods which could further improve the accuracy of the CO(2–1) auto spectrum, though we note the extensive literature on the matter particularly in the context of [C ii] intensity mapping (see, e.g., Gong et al. 2014; Breysse et al. 2015; Silva et al. 2015; Yue et al. 2015; Sun et al. 2018; Cheng et al. 2020). In particular, COMAP may be able to make use of cross correlations with Lyα emitters from the HETDEX survey (Hill et al. 2008, 2021; Gebhardt et al. 2021) to isolate the two signals (Chung et al. 2019, 2022; Silva et al. 2021). We also neglect for now any other possible contaminating lines. Chung et al. (2017) demonstrated that at 30 GHz, CO(1–0) at z ∼ 3 should dominate any other foreground lines. We do not expect this fact to change for the EoR CO(1–0) line.

3.4. Survey Sensitivity

To finalize our forecasts, we need to model the expected uncertainty on the above power spectra. Assuming pure white noise, the noise level σ_N in a given map voxel takes the form

$$\begin{eqnarray}&&{\sigma }_{N}=\displaystyle \frac{{T}_{\mathrm{sys}}}{\sqrt{{N}_{\mathrm{feeds}}\,\delta \nu {t}_{\mathrm{pix}}}},\end{eqnarray} \tag{ 17 }$$

where ${t}_{\mathrm{pix}}={t}_{\mathrm{obs}}^{\mathrm{eff}}({\sigma }_{\mathrm{beam}}^{2}/{{\rm{\Omega }}}_{\mathrm{field}})$ is the effective observing time per sky pixel. Note that we are using the effective observing time from Equation (2), so we are continuing to assume a single effective system temperature. N_feeds is the number of feeds in a single dish. Our noise field then has power spectrum

$$\begin{eqnarray}&&{P}_{N}={\sigma }_{N}^{2}{V}_{\mathrm{vox}},\end{eqnarray} \tag{ 18 }$$

where V_vox is the volume of a voxel, which is assumed to be a square ${\sigma }_{\mathrm{beam}}={\theta }_{\mathrm{FWHM}}/\sqrt{8\mathrm{ln}2}$ on a side and a single frequency-channel deep.

We also need to account for the high- and low-k cutoffs in our sensitivity induced by the limited spatial resolution and survey volume. For a given theoretical power spectrum P(k, μ), we can construct an observer-space spectrum

$$\begin{eqnarray}&&\tilde{P}(k,\mu )={W}_{\mathrm{vol}}(k,\mu ){W}_{\mathrm{res}}(k,\mu )P(k,\mu ),\end{eqnarray} \tag{ 19 }$$

where

$$\begin{eqnarray}\begin{array}{rcl}{W}_{\mathrm{vol}} & = & \left(1-\exp \left[-{\left(\displaystyle \frac{k}{{k}_{\perp }^{\min }}\right)}^{2}(1-{\mu }^{2})\right]\right)\\ & & \times \,\left(1-\exp \left[-{\left(\displaystyle \frac{k}{{k}_{\parallel }^{\min }}\right)}^{2}{\mu }^{2}\right]\right),\end{array}\end{eqnarray} \tag{ 20 }$$

and

$$\begin{eqnarray}&&{W}_{\mathrm{res}}={\eta }^{2}\exp \left\{-{k}^{2}\left[{\sigma }_{\perp }^{2}(1-{\mu }^{2})+{\sigma }_{\parallel }^{2}{\mu }^{2}\right]\right\}\end{eqnarray} \tag{ 21 }$$

cut off the measured power at low- and high-k, respectively. We use the forms of ${k}_{\perp }^{\min }$, ${k}_{\parallel }^{\min }$, σ_⊥, and σ_∥ from Bernal et al. (2019). In order to match effects seen in the Pathfinder survey, we have added an additional main beam efficiency factor η to account for the loss of power in the main beam. We adopt the Pathfinder value of η = 0.72 (Rennie et al. 2022) for the 30 GHz instruments, and assume the new 16 GHz instrument will have η ≈ 1. The W_vol and ${W}_{\mathrm{res}}$ values we assume at 30 GHz produce a sensitivity curve in good alignment with the optimistic 5 yr Pathfinder forecasts (Foss et al. 2022).

We have neglected in this section any explicit discussion of continuum foreground emission from sources like Milky Way dust and synchrotron or extragalactic radio point sources. We do not however expect continuum contamination to significantly affect a COMAP-EoR measurement. Any source of continuum emission will by definition appear as a signal strongly correlated between different frequency channels. As discussed in Section 3.3.3 of Foss et al. (2022), such correlated signals are effectively removed by the COMAP pipeline. An analogous argument was made in Keating et al. (2015) for the CO maps in a similar frequency range. This cleaning involves some loss of large-scale information, which is captured in our forecasts through the form of W_vol.

With our resolution effects in mind, we can now write the errors on our 16 and 30 GHz auto spectra as

$$\begin{eqnarray}&&{\sigma }_{16\,\mathrm{GHz}}^{2}(k,\mu )=\displaystyle \frac{1}{{N}_{\mathrm{modes}}(k,\mu )}{\left({\tilde{P}}_{1}^{\mathrm{EoR}}(k,\mu )+{P}_{N}^{16\,\mathrm{GHz}}\right)}^{2},\end{eqnarray} \tag{ 22 }$$

and

$$\begin{eqnarray}&&\begin{array}{l}{\sigma }_{30\,\mathrm{GHz}}^{2}(k,\mu )=\displaystyle \frac{1}{{N}_{\mathrm{modes}}(k,\mu )}\\ \ \times {\left({\tilde{P}}_{2}^{\mathrm{EoR}}(k,\mu )+{\tilde{P}}_{1,\mathrm{proj}}^{\mathrm{GA}}(k,\mu )+{P}_{N}^{30\,\mathrm{GHz}}\right)}^{2}.\end{array}\end{eqnarray} \tag{ 23 }$$

The number of modes N_modes available in a bin centered at k, μ is given by

$$\begin{eqnarray}&&{N}_{\mathrm{modes}}(k,\mu )={N}_{\mathrm{field}}\displaystyle \frac{{k}^{2}{\rm{\Delta }}k{\rm{\Delta }}\mu }{8{\pi }^{2}}{V}_{\mathrm{field}},\end{eqnarray} \tag{ 24 }$$

where V_field is the comoving volume of a single field and N_field = 3 accounts for the information from the three COMAP fields. We can similarly write the error on the cross spectrum as

$$\begin{eqnarray}&&{\sigma }_{\times }^{2}=\displaystyle \frac{1}{2}\left(\displaystyle \frac{{\tilde{P}}_{\times }^{2}}{{N}_{\mathrm{modes}}}+{\sigma }_{16\,\mathrm{GHz}}{\sigma }_{30\,\mathrm{GHz}}\right).\end{eqnarray} \tag{ 25 }$$

Following similar arguments to those in Bernal et al. (2019), we can after some algebra write the covariances between the auto and cross spectra

$$\begin{eqnarray}\begin{array}{rcl}{\sigma }_{\times -16\,\mathrm{GHz}}^{2}(k,\mu ) & = & \displaystyle \frac{1}{{N}_{\mathrm{modes}}}{\tilde{P}}_{\times }(k,\mu )\\ & & \times \,\left({\tilde{P}}_{1}^{\mathrm{EoR}}(k,\mu )+{P}_{N}^{16\,\mathrm{GHz}}\right),\end{array}\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&\begin{array}{l}{\sigma }_{\times -30\,\mathrm{GHz}}^{2}(k,\mu )=\displaystyle \frac{1}{{N}_{\mathrm{modes}}}{\tilde{P}}_{\times }(k,\mu )\\ \ \ \times \left({\tilde{P}}_{2}^{\mathrm{EoR}}(k,\mu )+{\tilde{P}}_{1,\mathrm{proj}}^{\mathrm{GA}}(k,\mu ),+{P}_{N}^{16\,\mathrm{GHz}}\right),\end{array}\end{eqnarray} \tag{ 27 }$$

and between the two auto spectra

$$\begin{eqnarray}&&{\sigma }_{16\,\mathrm{GHz}-30\,\mathrm{GHz}}^{2}=\displaystyle \frac{1}{{N}_{\mathrm{modes}}}{\tilde{P}}_{\times }^{2}.\end{eqnarray} \tag{ 28 }$$

Though we have worked in (k, μ) coordinates to this point, it is typically not possible to construct a well-defined μ from a curved-sky observation where the line-of-sight direction changes with telescope pointing. We will therefore express our sensitivity forecasts in terms of the multipoles of the power spectrum (Yamamoto et al. 2006; Bernal et al. 2019; Chung 2019). We can write the multipole ℓ of a given power spectrum as

$$\begin{eqnarray}&&{P}_{{\ell }}(k)=\displaystyle \frac{2{\ell }+1}{2}{\int }_{-1}^{1}P(k,\mu ){{ \mathcal L }}_{{\ell }}(\mu )d\mu ,\end{eqnarray} \tag{ 29 }$$

where ${{ \mathcal L }}_{{\ell }}(\mu )$ is the Legendre polynomial of degree ℓ. We will consider here the first three multipoles of each power spectrum, which we expect to contain the vast majority of the usable information. Our data thus consist of the monopoles, quadrupoles, and hexadecapoles of each of the 16 and 30 GHz auto spectra as well as the cross spectrum. We can write the covariance between each of these components as

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal C }}_{{\ell }{\ell }^{\prime} }^{{XY}}(k) & = & \displaystyle \frac{(2{\ell }+1)(2{\ell }^{\prime} +1)}{2}\\ & & \times \,{\displaystyle \int }_{-1}^{1}{\sigma }_{X}(k,\mu ){\sigma }_{Y}(k,\mu ){{ \mathcal L }}_{{\ell }}(\mu ){{ \mathcal L }}_{{\ell }^{\prime} }(\mu )d\mu .\end{array}\end{eqnarray} \tag{ 30 }$$

As an example to clarify this somewhat cumbersome notation, the covariance between the monopole (ℓ = 0) of the 30 GHz spectrum and the hexadecapole (${\ell }^{\prime} =4$) of the cross spectrum is given by

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal C }}_{04}^{30\,\mathrm{GHz}-\times }(k) & = & \displaystyle \frac{9}{2}{\displaystyle \int }_{-1}^{1}{\sigma }_{30\,\mathrm{GHz}}(k,\mu ){\sigma }_{\times }(k,\mu )\\ & & \times \,{{ \mathcal L }}_{0}(\mu ){{ \mathcal L }}_{4}(\mu )d\mu .\end{array}\end{eqnarray} \tag{ 31 }$$

The signal-to-noise ratio (S/N) for a given measurement is then

$$\begin{eqnarray}&&{\rm{S}}/{\rm{N}}={\left[\sum _{{ij}}{{\boldsymbol{d}}}_{i}^{T}{{ \mathcal C }}_{{ij}}{{\boldsymbol{d}}}_{j}\right]}^{1/2},\end{eqnarray} \tag{ 32 }$$

where d _i is the data vector constructed from the multipoles of the three power spectra in bins centered at k_i and ${{ \mathcal C }}_{{ij}}$ is the covariance matrix constructed from the components of Equation (30).

Here we will briefly summarize the models from the literature that we use to predict the range of possible CO signals at the EoR. For the full details behind each of these models, please see the relevant references. When a model does not contain all of the information we need for a full COMAP-EoR forecast (for example, only predicting CO(1–0) not CO(2–1)), we will make modest extensions to produce the necessary power spectra. As these results are intended simply to illustrate the range of possible S/Ns, the exact choice of these extensions should not significantly affect the overall picture. For our more detailed forecasts in later sections we will confine ourselves to models that inherently provide all of the relevant quantities. Figure 2 compiles the L(M) relations assumed for the CO(1–0) and CO(2–1) lines at reionization.

Zoom In Zoom Out Reset image size

4.1. Lidz et al. (2011)

Most of the models presented here use a method originally put forth by Lidz et al. (2011) to predict CO emission. Starting from a relationship between star formation rate (SFR) and halo mass, they use empirical scalings between SFR and infrared luminosity L_IR to get an estimate of L_IR(M), then use a final empirical scaling to go from L_IR to L_CO resulting in the following relation between CO luminosity and halo mass:

$$\begin{eqnarray}&&\displaystyle \frac{{L}_{\mathrm{CO}}(M)}{{L}_{\odot }}=2.8\times {10}^{3}\left(\displaystyle \frac{M}{{10}^{8}\ {M}_{\odot }}\right).\end{eqnarray} \tag{ 33 }$$

They assume a mass-independent duty cycle f_duty = 0.1, do not include any scatter about the mean relation, and assign the same mass–luminosity relation to both CO(1–0) and (2–1). This latter assumption means that the mean intensity of the CO(2–1) line will be eight times lower than the EoR CO(1–0) due to the factor of ν⁻³ from Equation (6). They explore a number of different ${M}_{\min }$ values and we assume their lowest value of 10⁸ M_☉ here. We also assume the same L(M) models for the galaxy assembly and EoR forecasts, though the signal will still evolve through the mass function. As discussed in Chung et al. (2022), this is likely a conservative estimate for this model, at least with regard to the EoR predictions, as the models from Pullen et al. (2013) predict a lower galaxy assembly era signal using effectively the same method.

4.2. Gong et al. (2011)

Gong et al. (2011) make use of a simulated CO catalog (Obreschkow et al. 2009). They fit the mean CO mass–luminosity relation with a double power law

$$\begin{eqnarray}&&L(M)={L}_{0}\left(\displaystyle \frac{{M}^{b}}{{M}_{c}}\right){\left(1+\displaystyle \frac{M}{{M}_{c}}\right)}^{-d},\end{eqnarray} \tag{ 34 }$$

where L₀, b, M_c , and d are fit parameters given in their Section 2. Fit values are given at z = 6, 7, and 8, we interpolate between them to estimate the fit at other redshifts. We find that this interpolation gives results in good agreement with their $\left\langle {Tb}\right\rangle$ calculations. For our higher-redshift CO(1–0) bin which is centered just above z = 8, we use their highest redshift values. Gong et al. (2011) do not assume a duty cycle or a scatter, and they set ${M}_{\min }={10}^{8}\ {M}_{\odot }$. Since their forecasts are redshift dependent, but do not extend to z ∼ 3, we assume the below model derived from Li et al. (2016) and Keating et al. (2020) (hereafter, Li16/Keating20) for ${P}_{1}^{\mathrm{GA}}$. As Li16/Keating20 is one of our brighter models, this will be a conservative choice as far as EoR sensitivity is concerned.

As with several of these models, Gong et al. (2011) only make predictions for the CO(1–0) line, so we will need to extrapolate their predictions to CO(2–1). On the low end, we could follow Lidz et al. (2011) and conservatively assign the same L(M) model to both lines. On the other hand, Pullen et al. (2013) argue that in the limit of high temperature and optical depth the CO(2–1) mean intensity would be eight times higher than that of CO(1–0). Given one of these two limits results in ${\left\langle {Tb}\right\rangle }_{2}=8{\left\langle {Tb}\right\rangle }_{1}$, and the other gives ${\left\langle {Tb}\right\rangle }_{2}={\left\langle {Tb}\right\rangle }_{1}/8$, we will split the difference here and assume the same mean brightness temperature for the two lines. This corresponds to increasing the L(M) fit from Equation (34) by a factor of 8. This will be our assumption for all models that do not include an explicit CO(2–1) prediction.

4.3. Mashian et al. (2015)

4.4. Sun et al. (2019)

4.5. Li et al. (2016)/Keating et al. (2020)

4.6. Yang et al. (2022)

Our final CO model, from Yang et al. (2022; hereafter Yang22), provides fitting functions optimized for intensity mapping based on semianalytic models (SAMs) from Yang et al. (2021). Unlike most of the above models, which rely heavily on empirical scalings, it attempts to self-consistently model the underlying physics that gives rise to CO emission. The SAMs are calibrated to a wide variety of galaxy observations, including lower-redshift CO lines. By providing fitting functions, this model enables easy application of the SAM results to intensity mapping forecasts like our work here. Mass-luminosity functions here take the form

$$\begin{eqnarray}&&\displaystyle \frac{L}{{L}_{\odot }}=2N\displaystyle \frac{M}{{M}_{\odot }}{\left[{\left(\displaystyle \frac{M}{{M}_{1}}\right)}^{-\alpha }+{\left(\displaystyle \frac{M}{{M}_{1}}\right)}^{\beta }\right]}^{-1}.\end{eqnarray} \tag{ 35 }$$

The double-power-law shape of Equation (35) is common to many data-driven treatments of star formation tracers (see, e.g., Moster et al. 2010; Padmanabhan 2018). Values of N, M₁, α, and β are provided for both CO lines, as well as separate fitting functions for σ_sc and f_duty. We assume ${M}_{\min }={10}^{10}\ {M}_{\odot }$, which here is set by the resolution limits of the semianalytic simulations.

We can now apply our power spectrum formalism to each of the above models to see how well the signals they predict can be detected by COMAP-EoR and COMAP-ERA. Figure 3 compiles all of the power spectrum forecasts from the above models, in both theoretical and observational form (i.e., with and without resolution effects) compared to sensitivity curves for our two survey concepts. We can clearly see the primary result of this exercise, that the available models produce an enormous range of possible signals, spanning over four orders of magnitude for CO(1–0). This highlights the extreme paucity of information about this period of cosmic history. It also justifies many of the simplifying assumptions made above, as their effect is almost certainly smaller than the range of signals seen here.

Zoom In Zoom Out Reset image size

Beyond the overall uncertainty, we can see several interesting features in the different models. The Lidz11 and Yang22 models bracket the range, particularly for CO(1–0) where the extra factor of 8 in the Lidz11 model significantly increases CO(1–0) relative to CO(2–1). For the other models, even where we have not forced the mean intensities to be identical they still predict similar levels. Compared to the other models, Sun19 has by far the most shot noise compared to its clustering amplitude due to its quite steep L(M) model, which peaks at higher halo masses than the others. Finally, the last three models are reasonably consistent with one another, which is perhaps surprising given the overall uncertainty.

Figure 3 only shows the EoR power spectra. As an example of the CO(1–0) interloper effect we compare the CO(2–1) auto spectrum from the Li16/Keating20 model to the projected galaxy assembly CO(1–0) in Figure 4. The effect is relatively small, but we can see by eye that the shapes of the two quadrupole spectra are slightly different, which lets us achieve some modest separation of the two signals. As shown in the bottom row of Figure 3, COMAP-EoR and particularly COMAP-ERA should have some quadrupole sensitivity on most of these models, so we should be able to take advantage of this effect.

Zoom In Zoom Out Reset image size

Table 2 provides the S/Ns obtained for COMAP-EoR and COMAP-ERA for the three power spectra in the z = 6.2 band. For the CO(2–1) auto S/N, we have effectively treated the interloper as an extra noise component. For example, this is why the Lidz11 model produces such a low S/N for the CO(2–1) auto spectrum, as the factor of 8 difference between CO(1–0) and (2–1) means that the low-redshift interloper is much brighter than the EoR line. This is not precisely correct, as we can perform internal cross correlations within our raw data to remove any overall bias from instrument noise (Ihle et al. 2022), while we cannot do the same with a signal on the sky. However, as mentioned above in a real analysis we could likely do more to mask out the lower-redshift line, which we have not accounted for here. We thus believe that treating the interloper as noise provides a sufficient approximation for our current level of detail.

From Figure 3 and Table 2, we can see that, in this redshift bin, COMAP-EoR performs quite well for all of the models except that of Yang22. The Lidz11 model is an outlier in the other direction with very high S/N, while the other models cluster in the S/N = 10–20 range. Table 3 gives the combined auto+cross S/Ns for each of the four frequency bands defined in Table 1 (where Band 3 corresponds to the redshift range for Table 2). Assuming there is negligible covariance between the four bands, we then sum the results in quadrature to obtain a total EoR detection significance for each model.

5.1. Parameter Constraints

Having studied the overall sensitivity of our planned survey, let us now examine how well we can measure the different components of the power spectrum individually. We will continue to follow a procedure analogous to Bernal et al. (2019), but given the very large uncertainty in the galaxy-evolution modeling part of this exercise we will hold all fundamental cosmological parameters fixed. We have access to three observables (CO(1–0) at galaxy assembly and EoR and CO(2–1) at EoR), and we have a clustering amplitude $\left\langle {Tb}\right\rangle$, average bias b, and shot noise amplitude P^shot for each, in addition to the cross-shot power ${P}_{\times }^{\mathrm{shot}}$. Each of these factors provide unique information about the luminosity distribution of the emitting halos.

We note, we believe for the first time in the literature, that two of these parameters are exactly degenerate even with the benefit of cross correlation. Specifically, the CO(2–1) shot power cannot be separated at our current level of detail from the lower-redshift CO(1–0) shot power. Because Poisson power is scale independent, we can only measure a single, overall shot amplitude in the 30 GHz power spectrum. There are not scale-dependent features to shift in the redshift projection, and the shot power only appears in the monopole, so they cannot be separated through the power spectrum anisotropy. We therefore introduce a new composite quantity

$$\begin{eqnarray}&&{S}_{\mathrm{tot}}={P}_{2}^{\mathrm{shot},\mathrm{EoR}}+\displaystyle \frac{1}{{\alpha }_{\parallel }{\alpha }_{\perp }^{2}}{P}_{1}^{\mathrm{shot},\mathrm{GA}},\end{eqnarray} \tag{ 36 }$$

which is the actual value measurable from the set of power spectra we consider here.

We thus have a total of nine free parameters:

$$\begin{eqnarray}&&\begin{array}{l}{\theta }_{i}=\left[{\left\langle {Tb}\right\rangle }_{1}^{\mathrm{EoR}},{\left\langle {Tb}\right\rangle }_{2}^{\mathrm{EoR}},{\left\langle {Tb}\right\rangle }_{1}^{\mathrm{GA}},{P}_{1}^{\mathrm{shot},\mathrm{EoR}},{S}_{\mathrm{tot}},\right.\\ \qquad \times \,\left.{P}_{\times }^{\mathrm{shot}},{b}_{1}^{\mathrm{EoR}},{b}_{2}^{\mathrm{EoR}},{b}_{1}^{\mathrm{GA}}\right].\end{array}\end{eqnarray} \tag{ 37 }$$

We will primarily use the Li16/Keating20 model to demonstrate the kinds of parameter constraints and science results we could obtain from COMAP-EoR. This is a relatively bright (though not the brightest) model, but more importantly it was shown to be broadly consistent with the best existing LIM data for CO(2–1) and above from mmIME (Keating et al. 2020), and perhaps slightly underestimates CO(1–0) LIM observations from COPSS (Keating et al. 2016). For some of our examples, we will also include predictions from the Yang22 semianalytic model. As we will see, this model appears to be quite pessimistic even compared to the handful of existing high-redshift CO direct observations, so it provides something of a worst-case scenario for CO at EoR. These are also the most recently published models from our set, and they both provide all of the information we need for our multi-line, multi-redshift forecasts without any additional assumptions.

We can obtain quantitative parameter forecasts by computing the Fisher matrix,

$$\begin{eqnarray}&&{F}_{{ij}}=\displaystyle \frac{\partial {{\boldsymbol{d}}}^{T}}{\partial {\theta }_{i}}{{ \mathcal C }}^{-1}\displaystyle \frac{\partial {\boldsymbol{d}}}{\partial {\theta }_{j}},\end{eqnarray} \tag{ 38 }$$

where d and ${ \mathcal C }$ are the complete data vector and covariance matrix from Section 3. The Fisher matrix is then the inverse covariance matrix for the chosen parameters θ_i . Given the huge modeling uncertainty above, we assume completely flat priors on all parameters.

The full nine-parameter forecast for the Li16/Keating20 model can be found in Appendix B; for the sake of readability we will highlight some important aspects of it here separately. As can be seen in Figure 3, the COMAP sensitivity is highest in the lower-k clustering regime of the power spectrum (see also Chung et al. 2022). Figure 5 shows the Fisher constraints on the three clustering amplitudes $\left\langle {Tb}\right\rangle$ that appear in our measurement. As expected from the strong overall detection of this model, COMAP-EoR obtains a quite strong detection of all three clustering amplitudes. We have also highlighted the constraints we would obtain in this space if we only had access to the two auto spectra at 16 and 30 GHz, i.e., we neglected P_× in our forecast. As a result, the EoR CO(2–1) and galaxy assembly CO(1–0) become quite strongly correlated (though not perfectly correlated, due to the anisotropy). This illustrates the unique benefit of the cross-correlation ability of the COMAP-EoR plan.

Zoom In Zoom Out Reset image size

Though our focus in this work is primarily on reionization, this cross correlation has crucial benefits for the lower-redshift science as well, even in the case of the Yang22 model where the EoR signal is undetected. Figure 6 shows the constraints on the clustering amplitude and bias of CO(1–0) at galaxy assembly from COMAP-EoR for our two demonstration models. The lower-redshift signal only appears at 30 GHz, so what we have here is effectively the difference between adding and not adding in the 16 GHz instrument. Since the high-redshift signal is so uncertain, it adds quite a bit of extra error to our attempts to make a precise z ∼ 3 observation. We have neglected this effect in the other papers in this series, as the current sensitivity is still relatively low, but at COMAP-EoR sensitivity this could be a quite substantial effect. For the brighter Li16/Keating20 model, it makes the difference between separating out the mean intensity and bias and not. For the Yang22 model, this demonstrates that there are significant science benefits to the 16 GHz observation even if it does not make a strong detection itself.

Zoom In Zoom Out Reset image size

Though COMAP is in general a clustering-focused measurement, we do see some constraints on the shot powers as well. The S_tot combined 30 GHz shot amplitude is quite strongly constrained, though as stated above it cannot be separated into low- and high-redshift components. We see a ∼ 2σ detection of the CO(1–0) and cross-shot powers at EoR under the Li16/Keating20 model. COMAP-EoR does not have the sensitivity to separate out the mean intensity and bias factors at EoR, but COMAP-ERA does, at least under this fairly optimistic model.

Up to this point we have focused on measuring the CO power spectra for our various models. Now we will move on to examine what these power spectrum measurements will tell us about the nature of high-redshift galaxies. We explore two critically important questions: what is the total abundance of star-forming molecular gas during reionization, and what type of galaxies contribute most to that measurement?

6.1. Cosmic Molecular Gas Abundance

One of the primary uses of any CO observation, including this one, is to use the molecule as a proxy for molecular gas, which itself is highly correlated with star formation activity. We can thus convert our CO power spectrum constraints into a measurement of the cosmic molecular gas history.

Our space of possible CO models is quite broad, and has neither the consistent parameterization nor the solid empirical basis we would need to carry out a full Bayesian calculation along the lines of the early Pathfinder results (Chung et al. 2022). Instead we will follow an analogous procedure to that used in the mmIME observations (Keating et al. 2020), wherein we assume for now that we precisely know the relationship between CO emission and molecular gas abundance. We will assume we know the shape of L(M) up to an overall amplitude, and we assume a known constant ratio α_CO between CO luminosity and molecular gas mass. This will result in a somewhat optimistic molecular gas forecast, but one that we argue is most consistent with previously published results. By using the mmIME prescription our forecasts will be directly comparable to their results. Measurements based on direct galaxy surveys also typically (Aravena et al. 2019; Riechers et al. 2019) assume a constant α_CO, and the assumption of known L(M) is roughly analogous to the assumption that the properties of low-luminosity CO sources can be extrapolated from bright detected objects. We explore some of the impact of these assumptions in Appendix C, and more detailed discussions of this type of calculation can be found in Breysse et al. (2022).

Quantitatively, we assume that the observed mass–luminosity relation for a given model is given by

$$\begin{eqnarray}&&L(M,z)=A(z){L}_{0}(M,z),\end{eqnarray} \tag{ 39 }$$

where L₀(M, z) is the default model relation used above, and any difference between the predicted and observed CO signal enters through the new parameter A(z). The cosmic molecular gas abundance is then given by

$$\begin{eqnarray}&&{\rho }_{{\rm{H}}2}(z)=A(z){\alpha }_{\mathrm{CO}}{\int }_{M\min }^{\infty }L{{\prime} }_{\mathrm{CO}(1-0)}(M,z)\displaystyle \frac{{dn}}{{dM}}{dM},\end{eqnarray} \tag{ 40 }$$

where

$$\begin{eqnarray}&&\displaystyle \frac{{L}_{\mathrm{CO}(1-0)}}{{L}_{\odot }}=4.9\times {10}^{-5}\displaystyle \frac{L{{\prime} }_{\mathrm{CO}(1-0)}}{{\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1}\ {\mathrm{pc}}^{2}}\end{eqnarray} \tag{ 41 }$$

relates the CO luminosity in physical units to that in the observer units commonly used for α_CO. We can thus perform a new Fisher forecast for A(z) and easily convert to ρ_H2. Here we will make the same assumption as mmIME that ${\alpha }_{\mathrm{CO}}\,=3.6\ {M}_{\odot }{({\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1}\ {\mathrm{pc}}^{2})}^{-1}$, based on measurements of individual high-redshift star-forming galaxies (Daddi et al. 2010). Appendix C discusses some of the uncertainty in this quantity, as does Breysse et al. (2022).

This is of course a very simplistic way to deal with a highly complex galaxy evolution problem. Breysse et al. (2022) showed that for mmIME, when carrying out this type of procedure, reasonable changes to the underlying CO model altered the final ρ_H2 result by considerably more than the statistical errors. The α_CO scaling alone is far more than a single redshift- and mass-independent constant (Bolatto et al. 2013). However, as we have repeated many times, much of this uncertainty will be swept up in the multiple-order-of-magnitude difference in signal between different models. In addition, this procedure is qualitatively similar to the assumptions made in common direct-detection CO analyses, where, for example, the choice of α_CO contributes significant systematic uncertainty (see, e.g., Boogaard et al. 2021)

Figure 7 shows forecasts for ρ_H2 measurements from COMAP-EoR and COMAP-ERA assuming the Li16/Keating20 and Yang22 models as a function of redshift. Forecasts are compared to both existing galaxy surveys using CO and dust observations and to the COPSS and mmIME CO LIM results. Broadly speaking, the Li16/Keating20 model is most consistent with the z ∼ 3 LIM data, while the Yang22 model is closer to the z ∼ 3 direct data. It should be noted that the Yang et al. (2021) semianalytic model makes its own prediction for ρ_H2(z), which may differ from what is plotted here, as the semianalytic models provide their own α_CO values, which are allowed to scale with mass and redshift. Plotted here is what we would infer from the Yang22 CO abundance assuming a constant α_CO.

Zoom In Zoom Out Reset image size

Comparing the two EoR predictions, we see that the Li16/Keating20 model evolves quite shallowly with redshift, while the Yang22 model falls off even more steeply than the direct observation points. The baseline COMAP-EoR survey provides an extremely tight measurement for the brighter model, but only a weak constraint on the lowest-redshift bin in our worst-case model. It takes COMAP-ERA to trace the full evolution of the extremely faint Yang22 case. Both models at z ∼ 3 perform dramatically better than any current measurement, though it is important to remember that the previously discussed model dependence is not included in either the current or our forecasted error bars.

6.2. Faint Star-forming Galaxies

In Figure 7, the Li16/Keating20 model with constant α_CO implies significantly more molecular gas than has been directly detected to date, in particular compared to the COLDz survey (Riechers et al. 2019), which provides the highest redshift data to date. A similarly flat evolution of the cosmic star formation rate density has been postulated using gamma-ray burst counts (Kistler et al. 2009). In order for that to be the case, there would need to be a significant reservoir of molecular gas present in galaxies too faint to appear in COLDz. This would in turn have important implications for the nature of reionization, as we would expect quite a bit more ionization-producing star formation activity with corresponding requirements on escape fraction (McQuinn 2016), and that activity would be concentrated in smaller but more numerous sources. This possibility gets to one of the key motivations for the concept of LIM in general, which is that it can constrain galaxy populations too faint to observe directly. In fact, there is already weak evidence at z ∼ 3 for excess star formation appearing in CO and [C ii] LIM data compared to what we would predict based on galaxy surveys alone (Breysse et al. 2022 and Yang et al. 2021 also compare the COPSS and mmIME ρ_H2 values in Figure 7 to direct measurements at the same redshifts). In this section we will quantify the search for excess faint emission.

Figure 8 shows the luminosity functions of our two demonstration models in the z = 6.2 redshift bin. We can compare these models to the detection limit of a hypothetical CO(1–0) deep field observed with the ngVLA, as described in Decarli et al. (2018). Then we can ask what is the total CO intensity we would obtain from only those galaxies brighter than the ngVLA limit. In other words, by comparing the luminosity functions measured from the LIM survey and the direct-detection survey we can see how much CO emission (and therefore star formation) is missed in the direct observations. Figure 9 shows what happens when we do this. Since this is a LIM-focused paper, we compute the clustering amplitude $\left\langle {Tb}\right\rangle$ for both the entire galaxy population and for only those galaxies brighter than the ngVLA limit. We then compare the difference between the two to the COMAP errors. We also make a rough approximation of the error on the ngVLA-determined amplitude assuming Poisson statistics and the 2 deg² COSMOS-spanning survey described in Decarli et al. (2018). For this survey area and the Decarli et al. (2018) detection limits, the Li16/Keating20 model predicts a few thousand CO(1–0) detections, while the Yang22 model predicts a few tens of detections, again highlighting the huge uncertainty in this model space.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For the Li16/Keating20 model, we see that there is a substantial amount of CO being missed by the direct ngVLA survey. One could always attempt to extrapolate the ngVLA luminosity function to lower values, but only a LIM survey like COMAP-EoR could see these faint sources directly. In the extremely faint Yang22 model, COMAP-EoR instead provides a definitive upper limit on how many faint sources there could possibly be. Even in the worst-case scenario for EoR detection, this upper limit is still extremely scientifically valuable in ruling out the population seen in Li16/Keating20.

Prospects for CO intensity mapping during the EoR are clearly highly dependent on the highly uncertain signal amplitude. The current COMAP observing strategy, targeting a handful of relatively small fields, is designed to optimize for the detection of faint signal. If the true signal is as faint as the Yang22 model predicts, EoR observations will likely remain in this regime up to at least the COMAP-ERA timescale. For the rest of the models we consider here, however, there may be motivation in the long term to expand the selected fields to enable more cross-correlation opportunities.

COMAP-EoR is primarily designed to map the cross correlation between CO(1–0) and CO(2–1), but this is far from the only interesting cross correlation in this redshift range. As discussed in Silva et al. (2021), a cross correlation with z ∼ 3 Lyα emitters from HETDEX can expand COMAP measurements of both lower-redshift and higher-redshift CO. Lidz et al. (2011) originally proposed a CO LIM project with the goal of cross correlating with 21 cm experiments like the Hydrogen Epoch of Reionization Array (DeBoer et al. 2017). These two surveys in combination could uniquely measure the typical size of ionized bubbles during the EoR. Though there are practical difficulties with the differing resolutions between CO and 21 cm observations, the high sensitivity of COMAP-ERA may make such a measurement possible.

LIM surveys of [C ii] like the Tomographic Intensity Mapping Experiment (Crites et al. 2014), the Carbon [C ii] line in the Post-reionization and Reionization epoch survey (Lagache 2018), and the Fred Young Submillimeter Telescope (CCAT-Prime Collaboration et al. 2021) at z ∼ 7 and the Experiment for Cryogenic Large-Aperture Intensity Mapping (Cataldo et al. 2021) at z ∼ 3 will also map unresolved emission from star-forming galaxies, providing a complement to CO. Recent literature has proposed intensity mapping of other species at reionization as well, including fine-structure lines of O iii (Padmanabhan et al. 2021) and rotational transitions of hydrogen deuteride (Breysse et al. 2021). While true cross correlations between all of these surveys may be practically difficult, even combining all of these measurements in same model of galaxy evolution will provide powerful insight into the high-redshift ISM.

From our Fisher forecasts in Figures 5, 6, and 10, we can clearly see the benefits of the multi-line cross-correlation approach of COMAP-EoR, advantages that are not so readily available to other LIM targets. The allowed volume of parameter space reduces dramatically when cross correlating the two CO lines over the case in which the two auto spectra are measured separately. Even when the EoR signal is too faint to detect, it is necessary to actually carry out the observation at both frequencies to make a high-precision measurement of the z ∼ 3 galaxy assembly era.

As mentioned above, however, we have also identified a new limitation of this approach, and indeed all similar LIM cross correlations. Because the shot noise in the cross spectrum has a unique value which cannot be directly determined from the two auto-shot amplitudes, there is no way to separate the shot noise levels between a target line and an interloping foreground line. This fact may be particularly relevant for the above [C ii] surveys, as they contain several different CO rotational transitions as interlopers to their EoR signals. It also may have implications for mmIME-like small-area surveys, which are only sensitive to shot power, as this limits their ability to separate out power spectra of different lines. In both of these cases, individual cross-shot powers will be accessible by cross-correlating pairs of tracers, but at least with cross correlation and anisotropy alone it will not be possible to isolate any individual auto-shot noise amplitude. We leave for future work an examination of other interloper-cleaning methods such as voxel masking.

Given the huge range of models shown above and the comparative faintness of the Yang22 model, it is clearly possible that EoR CO emitters will be too faint and rare to be detectable by COMAP-EoR. Beyond the simple prescriptions discussed here, effects like metallicity evolution and CMB backlighting (da Cunha et al. 2013) may act to push the signal in the fainter direction, particularly for low-mass galaxies. Despite this possibility, we argue here that even an upper limit set by a COMAP-EoR-type measurement is still extremely valuable. All but one of the existing LIM models we consider are bright enough to detect, so we would need a LIM observation to rule them out if nothing else. As shown in Figure 6, a high-redshift upper limit is also critically important to maximizing the scientific gain from the z ∼ 3 CO maps discussed in the rest of this series. Finally, even once we have access to high-quality direct observations from ngVLA (which are currently scheduled to appear on the timescale of COMAP-ERA, after the nominal COMAP-EoR campaign), a LIM upper limit provides confidence that the direct survey has indeed detected the bulk of cosmic star formation during reionization.

The CO power spectra we discuss here do not constitute the entirety of the information available to us in the COMAP-EoR data. Power spectra are two-point statistics, but as shown in Ihle et al. (2019) we can significantly improve a LIM measurement by including the one-point statistics as well, using the voxel intensity distribution (VID) formalism developed in Breysse et al. (2017). The power spectra above constrain only the first two moments of the CO luminosity functions, whereas the VID is sensitive to the full distribution. A joint one- and two-point analysis would improve our ability to constrain more sophisticated models of galaxy evolution, and provide a more detailed description of any faint populations like those shown in Section 6.2. More recently, Breysse et al. (2019) proposed a conditional one-point formalism that acts as a VID equivalent for the cross correlation. Such a conditional VID estimator or a continuous analog thereof would allow further non-Gaussian probes of the relationship between the two lines we study here.

We have presented an overview of the COMAP-EoR experiment, the next phase of the COMAP effort. By adding additional sensitivity at 30 GHz and a new observing band at 16 GHz, we can leverage the ladder of CO rotational transitions to map the CO(1–0) and (2–1) lines in overlapping volumes during the EoR. By examining CO emission models from the literature, we have shown that the strength of CO emission at z ∼ 5–8 is highly uncertain, a consequence of our poor understanding of the interstellar medium in these early galaxies.

In order to predict how well COMAP-EoR, and its hypothetical successor COMAP-ERA, could detect these models, we carried out the most detailed forecast to date of a cross correlation between two intensity maps. For all but one of the models we consider, we predict a highly significant detection of the EoR signal, with S/N ≳ 20. With such a strong detection, we can isolate the different components of the LIM signal, measuring multiple moments of the two CO luminosity functions. Using this measurement, we can make a uniquely complete measurement of the total reservoir of star-forming molecular gas during the EoR, and dramatically improve our measurement of the same quantity at z ∼ 3 over what is possible with the COMAP Pathfinder (Chung et al. 2022). For the faintest model we consider, COMAP-EoR obtains only upper limits on the reionization signal, but we show that even in this worst-case scenario it is still critically important to carry out the measurement in order to remove degeneracies on the z ∼ 3 measurement and to conclusively rule out the possibility of CO emission from galaxies too faint to observe directly.

Under our most pessimistic assumptions, a survey with the sensitivity of COMAP-ERA will be required to approach a CO detection. For all of the other above models, however, a ∼100σ measurement with COMAP-ERA would open up numerous opportunities for detailed study of the high-z ISM. At this level, the redshift evolution of the EoR could be followed extremely precisely, and a Lidz et al. (2011)-style cross correlation between CO and 21 cm may become possible. More detailed work will be necessary to catalog the possibilities of such a deep LIM observation.

Our extreme ignorance about the CO LIM signal at reionization should serve in its own right as motivation for a survey like COMAP-EoR. The brightest and faintest models we consider here represent wildly different predictions for the nature of star formation during reionization, and therefore the nature of the EoR itself. As we have shown, COMAP-EoR will radically reduce this uncertainty and open a key new window into this latest cosmic phase transition.

This material is based upon work supported by the National Science Foundation under grant Nos. 1517108, 1517288, 1517598, 1518282, and 1910999, and by the Keck Institute for Space Studies under "The First Billion Years: A Technical Development Program for Spectral Line Observations." Parts of the work were carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration, and funded through the internal Research and Technology Development program.

P.C.B. is supported by the James Arthur Postdoctoral Fellowship. D.T.C. is supported by a CITA/Dunlap Institute postdoctoral fellowship. The Dunlap Institute is funded through an endowment established by the David Dunlap family and the University of Toronto. H.I. acknowledges support from the Research Council of Norway through grant 251328. H.P. acknowledges support from the Swiss National Science Foundation through Ambizione Grant PZ00P2_179934. Work at the University of Oslo is supported by the Research Council of Norway through grants 251328 and 274990, and from the European Research Council (ERC) under the Horizon 2020 Research and Innovation Program (grant agreement No. 819478, Cosmoglobe). J.G. acknowledges support from the Keck Institute for Space Science, NSF AST-1517108 and University of Miami and Hugh Medrano for assistance with the cryostat design. J.L. acknowledges support from NSF Awards 1518282 and 1910999. L.K. was supported by the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No. 885990. We thank Isu Ravi for contributions to the warm electronics and antenna drive characterization.

The authors would like to thank Shengqi Yang, Anthony Pullen, Rachel Somerville, and Abhishek Maniyar for useful discussions.

Finally, we would like to thank an anonymous referee whose comments and suggestions significantly improved this manuscript.

Here we derive the form of the scattering effect on the cross-shot power seen in Equation (12). We follow similar derivations of the effect on the auto spectrum from Li et al. (2016), Sun et al. (2019), and Breysse et al. (2022).

Assume the observed CO line luminosities of a galaxy with halo mass M are ${L}_{1}={x}_{1}{L}_{1}^{0}(M)$ and ${L}_{2}={x}_{2}{L}_{2}^{0}(M)$ for CO(1–0) and CO(2–1), respectively, where L⁰(M) is the mean luminosity for a galaxy in a halo with mass M and x₁ and x₂ account for all of the factors besides halo mass that determine a galaxy's luminosity. Following Li et al. (2016) and our assumptions for the auto spectrum, we assume x₁ and x₂ are distributed lognormally. Formally, we model $\mathrm{log}({x}_{1})$ and $\mathrm{log}({x}_{2})$ as draws from a bivariate normal distribution with covariance matrix

$$\begin{eqnarray}{{ \mathcal C }}_{12}^{\mathrm{sc}}=\left(\begin{array}{cc}{\sigma }_{\mathrm{sc},1}^{2} & {r}_{\mathrm{sc}}{\sigma }_{\mathrm{sc},1}{\sigma }_{\mathrm{sc},2}\\ {r}_{\mathrm{sc}}{\sigma }_{\mathrm{sc},1}{\sigma }_{\mathrm{sc},2} & {\sigma }_{\mathrm{sc},2}^{2}\end{array}\right),\end{eqnarray} \tag{ A1 }$$

where r_sc is the correlation coefficient between the two distributions.

We can then write the scattered cross-shot term as

$$\begin{eqnarray}&&\begin{array}{l}{P}_{\times }^{\mathrm{shot}}={C}_{\mathrm{LT},1}(z){C}_{\mathrm{LT},2}(z){\displaystyle \int }_{M\min }^{\infty }{\displaystyle \int }_{0}^{\infty }{\displaystyle \int }_{0}^{\infty }\\ \ \ \times \,{x}_{1}{L}_{1}^{0}(M,z){x}_{2}{L}_{2}^{0}(M,z){f}_{\mathrm{duty}}(M,z)\\ \ \ \times \,\displaystyle \frac{{dn}}{{dM}}{{ \mathcal P }}_{\mathrm{sc}}({x}_{1},{x}_{2}){{dx}}_{1}{{dx}}_{2}{dM}.\end{array}\end{eqnarray} \tag{ A2 }$$

The joint lognormal probability distribution over x₁ and x₂ is

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal P }}_{\mathrm{sc}}({x}_{1},{x}_{2}) & = & \exp \left[\displaystyle \frac{1}{2}{\left({\boldsymbol{x}}-{\boldsymbol{\mu }}\right)}^{T}{{ \mathcal C }}_{12}^{\mathrm{sc}}\left({\boldsymbol{x}}-{\boldsymbol{\mu }}\right)\right]\\ & & \times \,\left(2\pi {x}_{1}{x}_{2}\sqrt{\det {{ \mathcal C }}_{12}^{\mathrm{sc}}}\right),\end{array}\end{eqnarray} \tag{ A3 }$$

where x = (x₁, x₂) and ${\boldsymbol{\mu }}=({\sigma }_{\mathrm{sc},1}^{2}\mathrm{log}(10)/2,{\sigma }_{\mathrm{sc},2}^{2}\mathrm{log}(10)/2)$ is the mean value of x set to enforce $\left\langle {{xL}}^{0}\right\rangle =\left\langle L\right\rangle$. Carrying out the integrals over x₁ and x₂ in Equation (A2) yields

$$\begin{eqnarray}&&\begin{array}{l}{P}_{\times }^{\mathrm{shot}}={C}_{\mathrm{LT},1}(z){C}_{\mathrm{LT},2}(z){e}^{{r}_{\mathrm{sc}}{\sigma }_{\mathrm{sc},1}{\sigma }_{\mathrm{sc},2}{\left[\mathrm{ln}(10)\right]}^{2}}\\ \ \ \times \,{\displaystyle \int }_{M\min }^{\infty }{L}_{1}(M,z){L}_{2}(M,z){f}_{\mathrm{duty}}(M,z)\displaystyle \frac{{dn}}{{dM}}{dM},\end{array}\end{eqnarray} \tag{ A4 }$$

which is the form from Equation (12).

Figure 10 shows the full nine-parameter Fisher matrix output for the Li16/Keating20 model at z = 6.2. Results for the worst-case Yang22 model are qualitatively similar, but with proportionally lower significance.

Zoom In Zoom Out Reset image size

The constraints on the cosmic molecular gas history from Section 6.1 make fairly strong assumptions about our knowledge of the relation between CO emission and ρ_H2(z). Though our choice of assumptions is generally in line with other forecasts and measurements in the literature (most notably Keating et al. 2020), it is nevertheless worth examining the model dependence of our conclusions.

Any conversion between CO luminosity and molecular gas mass will depend strongly on the assumed value of α_CO. In Figure 7, the existing constraints from intensity mapping (Keating et al. 2020) and many of those from direct surveys (e.g., Riechers et al. 2019) assume a constant α_CO across all halos and redshifts. In reality, α_CO is known to have quite strong variations not captured by this simple assumption (Bolatto et al. 2013).

To get an idea of how much this oversimplification might affect a COMAP-EoR measurement, the left panel of Figure 11 shows our Li16/Keating20 forecast for different values of α_CO, including a Milky Way-like value of ${\alpha }_{\mathrm{CO}}=4.3\ {M}_{\odot }\ {({\rm{K}}\ \mathrm{km}\,{{\rm{s}}}^{-1}\ {\mathrm{pc}}^{2})}^{-1}$ (Frerking et al. 1982; Dame et al. 2001; Bolatto et al. 2013) and a ULIRG-like ${\alpha }_{\mathrm{CO}}=0.8\ {M}_{\odot }\ {({\rm{K}}\ \mathrm{km}\,{{\rm{s}}}^{-1}\ {\mathrm{pc}}^{2})}^{-1}$ (Downes & Solomon 1998). These constant assumptions are still simpler than reality (see, e.g., Figure 7 of Breysse et al. 2022), but they provide an idea of the scale of the possible effect.

Zoom In Zoom Out Reset image size

On a similar note, in Equation (39), we assumed for the purposes of molecular gas forecasting that the shape of L(M) was known up to an overall amplitude. Though this is the same assumption used in the mmIME results in Keating et al. (2020), this is obviously not the case as the different models we consider here have quite different L(M) shapes. Qualitatively, this is similar in kind to the extrapolations that must be made about the population of galaxies too faint to detect in direct surveys like COLDz.

While we thus believe that the constraints from Figure 7 are those most directly comparable to the plotted previous measurements, there is also value in examining the impact of this L(M) assumption. The space of models at EoR redshifts is not mature enough to justify a full Bayesian treatment of this issue, so we will again attempt to make a rough qualitative estimate based on our literature models.

If we continue to assume a constant, linear conversion between L_CO and M_H2, then we guarantee that $\left\langle {T}_{1-0}\right\rangle \propto {\rho }_{{\rm{H}}2}$. We can therefore convert a Fisher forecast on $\left\langle {T}_{1-0}\right\rangle$ directly to a constraint on ρ_H2. However, $\left\langle {T}_{1-0}\right\rangle$ in the power spectrum measurement is strongly degenerate with the average bias, a degeneracy that COMAP-EoR lacks the sensitivity to break at high redshift. Using our previously assumed flat priors on the power spectrum would overestimate the allowed range of bias values, giving an unreasonably large error bar on ρ_H2.

For the purposes of this estimate, we will assume a Gaussian prior on the average bias with a fractional error of 25%, roughly equivalent to the difference between the bias values for the Li16/Keating20 and Yang22 models. The right panel of Figure 11 shows the molecular gas constraints under this relaxed assumption. We see that the constraints worsen by a factor of a few, though the data remain quite constraining. In a real measurement, we expect to have access to more sophisticated models and more detailed experimental priors, including lower-redshift information from the COMAP pathfinder. We also expect to make use of information beyond the power spectrum such as one-point (Ihle et al. 2019) and cross-correlation (Silva et al. 2021) statistics. It is therefore reasonable to assume that a full measurement would lie somewhere between the two limits shown here.