Joint Power Spectrum and Voxel Intensity Distribution Forecast on the CO Luminosity Function with COMAP

The simulations used in this work consist of two components only, namely the target CO signal and random white noise with properties corresponding to the parameters described above. Explicitly, the noise per voxel is given by

$$\begin{eqnarray}&&{\sigma }_{T}=\displaystyle \frac{{T}_{\mathrm{sys}}}{\sqrt{\tau \,\delta \nu }}=\displaystyle \frac{{T}_{\mathrm{sys}}\sqrt{{N}_{\mathrm{pixels}}}}{\sqrt{{\tau }_{\mathrm{tot}}{e}_{\mathrm{obs}}{N}_{\mathrm{feeds}}\delta \nu }},\end{eqnarray} \tag{ 1 }$$

where T_sys is the system temperature, τ is the observation time per pixel, τ_tot is the total observation time, e_obs is the observation efficiency, N_feeds is the number of feeds, N_pixels is the number of pixels, and δν is the frequency resolution. This gives us σ_T ≈ 11 μK and 8 μK for the COMAP1 and COMAP2 phases, respectively. For simplicity we assume that the noise is evenly distributed over all voxels.

A voxel is the 3D equivalent of a pixel. Two of the dimensions correspond to a regular pixel on the sky, while the third dimension corresponds to a small range of redshifts from where line emission would redshift into a given frequency bin of our instrument.

Both instrumental systematics and astrophysical foreground contamination are neglected in the following. However, since our estimator is inherently simulation based, these effects can be added at a later stage when a sufficiently realistic instrument model is available. For discussion of foreground contamination in similar line intensity surveys see, e.g., da Cunha et al. (2013), Breysse et al. (2015, 2017), and Chung et al. (2017).

The signal component is based on the peak patch DM halo approach described by Bond & Myers (1996) and Stein et al. (2019), coupled to the L_CO(M_halo) model presented by Li et al. (2016). Additionally, we adopt the same cosmological parameters as the Li et al. (2016) analysis for the DM simulations: Ω_m = 0.286, Ω_Λ = 0.714, Ω_b = 0.047, h = 0.7, σ₈ = 0.82, and n_s = 0.96.

The DM simulations in this paper were created using the peak patch method of Bond & Myers (1996) and Stein et al. (2019). To cover the full redshift range of the COMAP experiment we simulated a volume of (1140 Mpc)³ using a particle-mesh resolution of N_cells = 4096³. Projecting this onto the sky results in a 96 × 96 field of view covering the redshift range 2.4 < z < 3.4, with a minimum DM halo mass of 2.5 × 10¹⁰ M_⊙.

The resulting halo catalog contains roughly 54 million halos, each with a position, velocity, and mass. The peak patch method can simulate continuous light cones on-the-fly, so stitching snapshots together was not required to create the light cone. Although peak patch simulations result in quite accurate halo masses, the DM halo catalogs were additionally mass corrected by abundance matching along the light cone to Tinker et al. (2008) to ensure statistically the same mass function as the simulations used in the Li et al. (2016) analysis. For a detailed study of the clustering properties of peak patch simulations and other approximate methods, see Lippich et al. (2019), Blot et al. (2018), and Colavincenzo et al. (2019).

A single run required 900 s of computation time on 2048 Intel Xeon EE540 2.53 GHz CPU cores of the Scinet-GPC cluster, with a memory footprint of ≃2.4 TB. This efficiency of the peak patch method allowed for 161 independent realizations of the full 1140 Mpc, ${N}_{\mathrm{cells}}={4096}^{3}$ volume, taking a total computation time of only ∼82,000 CPU hours, over three orders of magnitude faster when compared to an N-body method of equivalent size.

There are many approaches in the literature for estimating the expected CO signal based on DM halos (e.g., Righi et al. 2008; Obreschkow et al. 2009; Visbal & Loeb 2010; Carilli 2011; Gong et al. 2011; Lidz et al. 2011; Fu et al. 2012; Carilli & Walter 2013; Pullen et al. 2013; Breysse et al. 2014; Greve et al. 2014; Mashian et al. 2015; Li et al. 2016; Padmanabhan 2018), with resulting estimates of the CO luminosities spanning roughly an order of magnitude.

Here we adopt the model described by Li et al. (2016) to convert from simulated light cones populated with DM halos to observed CO brightness temperature. This model is defined by a set of parametric relations between DM halo masses, star formation rates (SFR), infrared (IR) luminosities, L_IR, and CO luminosities, L_CO.

The model uses the results from Behroozi et al. (2013a, 2013b) to obtain average SFR from DM halo masses and adds an additional log-normal scatter on top of the average, determined by σ_SFR. IR luminosities are then obtained through the relation

$$\begin{eqnarray}&&\mathrm{SFR}={\delta }_{\mathrm{MF}}\times {10}^{-10}\,{L}_{\mathrm{IR}}.\end{eqnarray} \tag{ 2 }$$

Further, to obtain CO luminosities, the relation

$$\begin{eqnarray}&&\mathrm{log}{L}_{\mathrm{IR}}=\alpha \mathrm{log}{L}_{\mathrm{CO}}^{{\prime} }+\beta \end{eqnarray} \tag{ 3 }$$

is used before a second round of log-normal scatter is added, determined by the parameter ${\sigma }_{{L}_{\mathrm{CO}}}$.

This gives us a L_CO(M_halo) model with five free parameters, $\theta =\{{\sigma }_{\mathrm{SFR}},\mathrm{log}{\delta }_{\mathrm{MF}},\alpha ,\beta ,{\sigma }_{{L}_{\mathrm{CO}}}\}$. The relation between L_CO and M_halo for our fiducial model parameters is shown in Figure 1. For more discussion of the physical and observational motivation for this model, see the original paper, Li et al. (2016).

Zoom In Zoom Out Reset image size

This model is applied to each DM halo separately and we create high-resolution maps from the resulting CO luminosities by converting the total luminosity in a given voxel into brightness temperature. These maps were created using the publicly available limlam_mocker code.¹¹ Next, we convolve these maps with the (Gaussian) instrumental beam profile, degrade to the low-resolution voxel size used in the analysis, and, finally, we add Gaussian uncorrelated noise with standard deviation σ_T, as specified above.

The estimated power spectrum, P(k_i), is calculated simply by taking the 3D Fourier transform of the temperature cube, binning the absolute squared values of the Fourier coefficients according to the magnitude of corresponding wave number k, and averaging over all the contributions within each bin. For a Gaussian map, the Fourier components within each bin follow a perfect normal distribution with mean zero and variance given by the value of the power spectrum. For a non-Gaussian field, the distribution of the Fourier components is more complicated, and thus the power spectrum does not contain all the statistical information in the map. We expect the CO signal to form a highly non-Gaussian map, therefore, in this paper we simply consider the power spectrum as a useful observable that carries some, but far from all, of the statistical information in the map.

As an observable, the power spectrum needs to be accompanied by a covariance matrix ${\xi }_{{ij}}^{P}\equiv \mathrm{Cov}(P({k}_{i}),P({k}_{j}))$ in the analysis, since there are correlations between the power spectrum at different k values.

We consider the VID as another observable, complementary to the PS and more closely related to the luminosity function.

Unlike in many other works on ${ \mathcal P }(D)$ analysis (e.g., Lee et al. 2009; Glenn et al. 2010; Vernstrom et al. 2014; Breysse et al. 2017; Leicht et al. 2019), we do not try to estimate the VID analytically, rather we estimate it based on simulations. This allows us to fully take into account the effects of the beam, clustering, and covariance between temperature bins in a very straightforward manner.

We consider two contributions to the VID, namely the CO signal itself and the instrumental noise. Together they result in the the full probability distribution of voxel temperatures, ${ \mathcal P }(T)$, where T is the observed brightness temperature from a voxel. Since we assume the noise to be uniformly distributed over all voxels in the observed field and assume that the CO signal itself is statistically homogeneous and isotropic, the total probability distribution, ${ \mathcal P }(T)$, is the same across all voxels.

The basic observable related to the VID are the temperature bin counts (i.e., the histogram of voxel temperatures), B_i. The expectation value of these are given by the VID itself,

$$\begin{eqnarray}&&\langle {B}_{i}\rangle ={N}_{\mathrm{vox}}{\int }_{{T}_{i}}^{{T}_{i+1}}{ \mathcal P }(T){dT},\end{eqnarray} \tag{ 4 }$$

where N_vox is the number of voxels observed and B_i is the number of voxels with a temperature between T_i and T_i+1.

If the temperatures of all the voxels that we sample were completely independent, then each of the voxel bins would be approximately independent and would follow a binomial distribution with variance ${\mathrm{Var}}_{\mathrm{ind}}({B}_{i})=\langle {B}_{i}\rangle (1-\langle {B}_{i}\rangle /{N}_{\mathrm{vox}})$. However, even in this ideal case they would not be perfectly independent. We only have a finite number of voxels, and, therefore, if one bin contains a number of voxels above average, then the other bins must have a number lower than average.

In general, the samples will not be independent for many reasons, including correlated sky or noise structures and processing effects, and we therefore need the full covariance matrix between bins, ${\xi }_{{ij}}^{B}\equiv \mathrm{Cov}({B}_{i},{B}_{j})$. This covariance matrix will depend on the DM density field, the CO bias, and the luminosity function, and we will estimate it using simulations.

The main missing component in the above method is definition of a joint power spectrum and VID covariance matrix. By having access to the computationally cheap yet realistic Monte Carlo simulations described above, we can approximate this matrix with simulations. In addition to giving us covariance matrices to do our analysis, this also allows us to check under what conditions the full covariance matrix is necessary and when we can get away with assuming that individual samples are independent.

In this paper, we start with 161 independent simulated light-cone cubes of DM halos, each covering about 96 × 96 on the sky and a frequency range between 26 and 34 GHz, corresponding to redshifts between 2.4 and 3.4. The frequency dimension is divided equally into 512 frequency bins, each spanning δν ≈ 15.6 MHz, corresponding to a redshift resolution of δz ≈ 0.002. Since the COMAP field only spans 15 × 15 on the sky, we sub-divide each of the 96 × 96 light-cone cubes, after beam convolution, into 36 square fields, each covering 15 × 15 , resulting in a total of 5796 semi-independent sky realizations. The final pixelization of these maps is a 22 × 22 grid of square pixels, resulting in a pixel size of $\delta \theta \approx 4\buildrel{\,\prime}\over{.} 1$. To these maps, we add uniformly distributed white noise at the appropriate levels for the COMAP1 and COMAP2 experiment setups described above.

When choosing the pixel size to use for the analysis, we follow Vernstrom et al. (2014). They show that, for ${ \mathcal P }(D)$ analysis, choosing a pixel size to be equal to the FWHM of the beam is a good tradeoff between picking a small pixel size to include the maximal information, and choosing a larger pixel size to reduce the pixel to pixel correlations induced by the beam.

We combine our two observables into a joint one-dimensional vector of the form

$$\begin{eqnarray}&&{d}_{i}=({P}_{{k}_{i}},{B}_{i}),\end{eqnarray} \tag{ 5 }$$

where ${P}_{{k}_{i}}$ is the binned power spectrum and B_i are the temperature bin counts. Let us first consider the ideal case in which all elements in this vector are independent and the Fourier components are approximately Gaussian. In that case we can compute the expected variance, which we will simply call the independent variance, analytically,

$$\begin{eqnarray}&&{\mathrm{Var}}_{\mathrm{ind}}({P}_{{k}_{i}})=\langle {P}_{{k}_{i}}{\rangle }^{2}/{N}_{\mathrm{modes}}({k}_{i}),\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{\mathrm{Var}}_{\mathrm{ind}}({B}_{i})=\langle {B}_{i}\rangle (1-\langle {B}_{i}\rangle /{N}_{\mathrm{vox}})\approx \langle {B}_{i}\rangle ,\end{eqnarray} \tag{ 7 }$$

where N_modes(k_i) is the number of modes in the ith k bin and where we have introduced the notation Var_ind(d_i) for this conditionally independent variance.

With this notation in hand, we define a "pseudo-correlation matrix" as

$$\begin{eqnarray}&&{c}_{{ij}}\equiv \displaystyle \frac{{\xi }_{{ij}}}{\sqrt{{\mathrm{Var}}_{\mathrm{ind}}({d}_{i}){\mathrm{Var}}_{\mathrm{ind}}({d}_{j})}},\end{eqnarray} \tag{ 8 }$$

where, as in Section 3.4, ξ_ij is the full covariance matrix. Note that c_ij is the exact correlation matrix in the limit that Var_ind(d_i) is the true full variance. An important advantage of the pseudo-correlation matrix, however, is the fact that Var_ind(d_i) may be estimated directly from the average data itself, and this is required for our MCMC procedure to be sufficiently fast.

The full covariance matrix ξ is estimated for the model described by Li et al. (2016), adopting the fiducial parameters θ₀, using the set of 5796 simulations described above. However, for the MCMC sampler described in Section 3.4, we actually need the full covariance matrix, corresponding to different model parameters θ, at each step in the Markov chain. Generating the full covariance matrix with the above procedure at each MC step is clearly not computationally feasible and we therefore need to approximate this somehow.

With regard to this last point, we introduce the following proposal: we assume that the full covariance matrix scales, under a change of model parameters from θ₀ to θ, the same way as the independent variance, Var_ind(d_i),

$$\begin{eqnarray}&&{\hat{\xi }}_{{ij}}(\theta )\approx {\xi }_{{ij}}({\theta }_{0})\displaystyle \frac{\sqrt{{\mathrm{Var}}_{\mathrm{ind}}^{\theta }({d}_{i}){\mathrm{Var}}_{\mathrm{ind}}^{\theta }({d}_{j})}}{\sqrt{{\mathrm{Var}}_{\mathrm{ind}}^{{\theta }_{0}}({d}_{i}){\mathrm{Var}}_{\mathrm{ind}}^{{\theta }_{0}}({d}_{j})}},\end{eqnarray} \tag{ 9 }$$

where ${\mathrm{Var}}_{\mathrm{ind}}^{{\theta }_{0}}({d}_{i})$ is the independent variance for the fiducial model and ${\mathrm{Var}}_{\mathrm{ind}}^{\theta }({d}_{i})$ is the independent variance for arbitrary parameters θ. Since this latter function only depends on the average quantities $\langle {d}_{i}\rangle$, it is computationally straightforward to compute ${\hat{\xi }}_{{ij}}(\theta )$ at any position in an MCMC sampler. Note also that ${\hat{\xi }}_{{ij}}(\theta )$ is, by construction, positive definite, as required for a proper covariance matrix.

For a noise-dominated experiment, where all samples are approximately independent, the independent variance, Var_ind(d_i), is the correct variance and Equation (9) is the correct scaling of the covariance matrix. However, we use this scaling as a first approximation even in cases where there is some covariance in the data.

Intuitively, Equation (9) is equivalent to postulating that the pseudo-correlation matrix, c_ij, is approximately constant (i.e., independent of the specific parameters in question). For real-world applications, we recommend testing this assumption explicitly by computing the covariance matrix by brute force simulation for a few extreme parameter combinations drawn from the Markov chains produced during the analysis.

The above prescription applies straightforwardly to single-field observations as, for instance, planned for COMAP1. In contrast, COMAP2 will, under our assumptions, span N = 4 independent but statistically identical fields. Since the mean vector of observables evaluated across those four fields equals the average of the four corresponding independent observable vectors, the full covariance matrix is simply given by the single-field covariance matrix divided by the number of fields:

$$\begin{eqnarray}&&{\xi }_{{ij}}^{N\,\mathrm{field}}=\tfrac{{\xi }_{{ij}}^{1\ \mathrm{field}}}{N}.\end{eqnarray} \tag{ 10 }$$

Note that c_ij then, assuming the fields are of the same size, only depends on the noise level per field, so for a given noise level per field, c_ij is independent of the number of fields.

Finally, we note that the total number of degrees of freedom in our joint PS+VID statistic is in this paper equal to 45, corresponding to 20 power spectrum bins and 25 VID bins. For this number of degrees of freedom, a set of 5796 (semi-independent) simulations provides a very good estimate of all numerically dominant components of the covariance matrix, including both the diagonal and the leading off-diagonal modes, and ξ_ij is well conditioned.

As previously mentioned, we use an MCMC algorithm to sample from the posterior distribution of the L_CO(M_halo) model parameters, $\theta =\{{\sigma }_{\mathrm{SFR}},\mathrm{log}{\delta }_{\mathrm{MF}},\alpha ,\beta ,{\sigma }_{{L}_{\mathrm{CO}}}\}$. This posterior distribution is, as usual, given by Bayes' theorem,

$$\begin{eqnarray}&&P(\theta | d)\propto P({d}_{i}| \theta )P(\theta ),\end{eqnarray} \tag{ 11 }$$

where d represents our compressed data set, $P(d| \theta )$ is the likelihood defined below, and P(θ) is some set of priors. We use the emcee package (Foreman-Mackey et al. 2013) and its implementation of an affine-invariant ensemble MCMC algorithm, with 142 walkers.

We use a burn-in period of 1000 steps, and use the next 1000 steps for the posterior estimation.

We assume a Gaussian likelihood for our observables d_i of the form (up to an additive constant)

$$\begin{eqnarray}&&-2\mathrm{ln}P(d| \theta )=\sum _{{ij}}[{d}_{i}-\langle {d}_{i}\rangle ]{\left({\xi }^{-1}\right)}_{{ij}}[{d}_{j}-\langle {d}_{j}\rangle ]+\mathrm{ln}| \xi | ,\end{eqnarray} \tag{ 12 }$$

where the means $\langle {d}_{i}\rangle$ depend on the model parameters θ, and the covariance matrix ξ_ij is approximated by the expression given in Equation (9). (Note that we do not need to assume that the low-level data are Gaussian, but only that the compressed observables may be well approximated by a multivariate Gaussian distribution. Due to the central limit theorem, this is in practice very often an excellent approximation.)

For both the power spectrum and the low and intermediate temperature VID bins, for which there is a large number of voxel counts per bin, this Gaussian approximation holds to a high degree. However, for the highest VID temperature bins, where there are only a few voxel counts per bin, the discrete nature of the bin count may become relevant and the full binomial distribution should, in principle, be taken into account. However, this effect can also be easily remedied by increasing the bin width, albeit at the cost of a slight loss of information, as is suggested in Vernstrom et al. (2014), and we therefore neglect it in the following, since our primary focus is the dominant Gaussian component of the likelihood. A more thorough analysis may take this issue into account either analytically or by simulations.

An advantage of using a Gaussian likelihood for the VID is that it gives us a straightforward way to take into account the correlations between temperature bins apparent in the covariance matrix, ξ_ij (e.g., in Figure 2). For another approach to building up a ${ \mathcal P }(D)$ likelihood, see Glenn et al. (2010).

Zoom In Zoom Out Reset image size

To estimate $\langle {d}_{i}\rangle$, we compute 10 maps of the survey volume at each step in the MCMC chain using the current model parameters θ with different DM halo realizations (randomly drawn from 252 independent catalogs corresponding to the survey volume). The specific number of realizations, 10 in our case, represents a compromise between minimizing the sample variance in the estimate of $\langle {d}_{i}\rangle$ and maintaining a reasonable computational cost per MC step. Finally, we bin all of the halos in the 10 realizations according to their luminosity and use this histogram to estimate the luminosity function at the current values of θ. This way the MCMC procedure gives us the luminosity function at different points in parameter space, sampled according to the posterior of the model parameters, which we can use to derive constraints on the luminosity function itself.

We adopt the same physically motivated priors as discussed by Li et al. (2016). Specifically, these read

$$\begin{eqnarray}&&P({\sigma }_{\mathrm{SFR}})={ \mathcal N }(0.3,0.1)\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&P(\mathrm{log}{\delta }_{\mathrm{MF}})={ \mathcal N }(0.0,0.3)\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&P(\alpha )={ \mathcal N }(1.17,0.37)\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&P(\beta )={ \mathcal N }(0.21,3.74)\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&P({\sigma }_{{L}_{\mathrm{CO}}})={ \mathcal N }(0.3,0.1),\end{eqnarray} \tag{ 17 }$$

where ${ \mathcal N }(\mu ,\sigma )$ corresponds to a Gaussian distribution with mean μ and standard deviation σ. Additionally, we require the two logarithmic scatter parameters, σ_SFR and ${\sigma }_{{L}_{\mathrm{CO}}}$, to be positive. We choose the mean of all these distributions as the fiducial model, θ₀.

To quantify the importance of joint PS+VID analysis, we perform the above analysis both with each observable separately and with the joint analysis. The main result in this paper may then be formulated in terms of the relative improvement on the CO luminosity function uncertainty derived from the joint analysis to those found in the independent analyses.

When calculating our observables (PS and VID), we assume that our survey volume can be treated as a rectangular grid of voxels with constant co-moving volume. We also neglect the evolution of our observables over redshifts between z = 2.4 and 3.4. That is, we assume that samples from different redshifts are drawn from the same distribution, whether they are power spectrum modes or voxel temperatures. We also assume that the instrument beam is achromatic and is equal to the value at the central frequency. This is of course just an approximation that we make in order for the analysis to be simple. If we were doing experiments with higher signal to noise, we might divide our data into two different redshift regions and do an independent analysis of each region. This could allow us to study the redshift evolution of the observables. For COMAP (1 and 2), however, we are probably best off combining all the data, like we do here, in order to increase the overall signal to noise.

Finally, since COMAP will not measure absolute zero levels, we subtract the mean from all maps. For the power spectrum, this has a negligible impact, as it simply removes one out of N_vox modes. However, it has a significantly higher impact for the VID. Specifically, it makes it much harder to distinguish a potential background of weak sources from noise. Indeed, as shown by Breysse et al. (2017), removing the monopole makes it much harder to detect a possible low luminosity cutoff in the CO luminosity function using the VID.

Figure 2 shows the pseudo-correlation matrices, c_ij, for our two experimental setups, as well as for pure signal alone, for reference. In order to illustrate the effect of the beam, we show covariance matrices from maps both without and with beam smoothing in the left and right columns, respectively.

The first thing to notice is that instrumental noise significantly reduces the numerical values of the normalized covariance matrices, bringing it closer to the independent white noise case for which c_ij = δ_ij. This agrees with intuition, since the noise itself is white and uncorrelated.

Beam smoothing also leads to weaker correlations. This is mainly due to the beam diluting the signal at small scales, where the correlation is otherwise strongest.

Next, we notice that the cross-correlations between the power spectrum and VID are of the same order of magnitude as the correlations internal to each observable itself. Thus, it is essential to account for all these correlations in any joint PS and VID analysis, as is done in the present paper.

Finally, we note that when designing an experiment like COMAP, one of the important trade-offs involves observation time per field. To obtain a fast signal detection it is in general advantageous to observe deep on the smallest possible field. However, this only holds true while the signal-to-noise per voxel is significantly less than unity. When the noise starts to become comparable to the signal, the signal-induced voxel–voxel correlations starts to become important, and the effective uncertainties no longer scale as ${ \mathcal O }(1/\sqrt{\tau })$, where τ is the observation time per pixel. Generally, in such a tradeoff, any significant correlations between different power spectrum modes or voxel temperatures will tend to favor larger survey area or multiple fields, both effectively leading to more independent samples, and thereby higher overall integration efficiency.

We are now ready to present both individual and joint PS+VID constraints on the CO luminosity function, which are summarized in Figure 3 for COMAP1 (left column) and COMAP2 (right column). The top row shows the constraints obtained from the power spectrum alone; the middle row shows the constrains obtained from the VID alone; and the third row shows the constraints from the joint analysis. In each panel, the shaded colored region shows the 95% credibility region from the MCMC samples and the solid line with the same color shows the posterior median. The purple solid line shows the average luminosity function obtained from the mean of all available halo catalogs, and thus represents the ensemble average of our input model. Note that the colored regions correspond to one single realization and the uncertainties therefore contain contributions from instrumental noise, cosmic variance, and sample variance. The agreement between the estimated confidence regions and the ensemble mean is quite satisfactory in all cases, with uncertainties that appear neither too large nor too small.

Zoom In Zoom Out Reset image size

Considering first the individual PS and VID estimates, shown in the top two rows, we see that the two observables are indeed complementary. In particular, the VID primarily constrains the high luminosity end of the luminosity function, while the power spectrum imposes relatively stronger constraints on the low luminosity end. This makes sense intuitively, since the VID is essentially optimized to look for strong outliers above the noise, whereas the power spectrum represents a weighted mean across the full field for each physical scale. It is interesting to note, however, that the VID provides, on average, stronger constraints on the luminosity function than the power spectrum does.

Due to this complementarity, the joint estimator provides the strongest constraints of all. To make this point more explicit, the fourth row compares the uncertainties of the independent power spectrum and VID analyses to the joint constraints. Of course, there is a significant amount of cosmic variance in each of these functions and the precise numerical value of the uncertainty ratio therefore varies significantly with luminosity; but the mean trend is clear: The individual analyses typically result in 20%–70% larger uncertainties than the joint analysis when averaged over luminosities between L_CO = 10⁴–10⁷ L_⊙. Over 10 cosmological realizations, the PS and VID resulted in, on average, 58% and 30% larger uncertainties (in dex) individually, than the joint analysis. This is the main novel result presented in this paper.

Lastly we present the constraints of the model parameters themselves. When doing the MCMC posterior mapping we explore the parameter space of the Li et al. (2016) L_CO(M_Halo) model. Figure 4 shows the posterior distribution for these parameters derived from one realization of the COMAP2 experiment (the same realization as the COMAP2 results in Figure 3).

Zoom In Zoom Out Reset image size

Results for PS, VID, and joint PS+VID analysis are shown in blue, red, and black, respectively. Prior distributions are shown in green. The two curves of each color correspond 68% and 95% credibility regions.

We see that the two parameters that are mainly constrained are α and β, the two parameters from the average L_CO–L_IR relation. These two parameters are fairly degenerate, and the direction in which they are degenerate is given roughly by the line α = −0.1β + 1.19 (Li et al. 2016). In Figure 5, we show the luminosity function for different points on this line. For the figure, the values of ${\sigma }_{\mathrm{SFR}},{\sigma }_{{L}_{\mathrm{CO}}}$, and $\mathrm{log}{\delta }_{\mathrm{MF}}$ are fixed at 0.3, 0.3, and 0.0, respectively. Although the overall signal strength, at least in terms of detectability, is fairly constant along this line, the shape of the luminosity function changes significantly. Lower values of alpha imply a more steep power-law relation between L_CO and L_IR leading to more sources with very high or very low luminosities. We see this as a flattening of the luminosity function. In such a case, a larger fraction of the overall signal will be given by high-luminosity sources.

Zoom In Zoom Out Reset image size

The other parameter that is also slightly constrained is the log-normal scatter parameter from the L_CO–L_IR relation, ${\sigma }_{{L}_{\mathrm{CO}}}$. This parameter is only slightly more constrained compared to the prior, with the highest values of ${\sigma }_{{L}_{\mathrm{CO}}}$ being disfavored. The posterior of the other scatter parameter, σ_SFR, is basically given by the corresponding prior (i.e., this parameter is not very well constrained by the experiment), although, as expected from the fact that the scatter parameters have basically the same effect, we see signs of the degeneracy between them in the posterior.

Interestingly, the normalization parameter in the SFR–L_IR relation, $\mathrm{log}{\delta }_{\mathrm{MF}}$, actually has a posterior that is wider than the prior. This may be because the best fit of this parameter from each of the different patches have an intrinsic scatter larger than the scatter in the prior. We note that we see the same effect in Li et al. (2016; their Figure 7).

From the mean relations in the Li et al. (2016) model, we have log L_CO ∼ −β − log δ_MF. Intuitively, we would then expect $\mathrm{log}{\delta }_{\mathrm{MF}}$ to be completely degenerate with β. However, since the SFR–L_IR is much better constrained by observations than the L_CO–L_IR relation, the prior on $\mathrm{log}{\delta }_{\mathrm{MF}}$ is much tighter than the one on β. The degeneracy thus prevents us from constraining $\mathrm{log}{\delta }_{\mathrm{MF}}\,{\text{}}{until}\,\beta$ is constrained to a comparable level.