A User's Guide to Extracting Cosmological Information from Line-Intensity Maps

TL;DR

This paper develops a comprehensive framework for analyzing line-intensity mapping data to extract cosmological information while accounting for astrophysical degeneracies, improving measurement precision of key parameters.

Contribution

It introduces a detailed model of the LIM power spectrum including anisotropies and degeneracy mitigation, enhancing the robustness of cosmological inference from LIM observations.

Findings

Including the hexadecapole improves BAO scale precision by 10-60%.

Fewer, wider bins are optimal for certain cosmological parameters.

The formalism enables robust cosmological constraints despite astrophysical uncertainties.

Abstract

Line-intensity mapping (LIM) provides a promising way to probe cosmology, reionization and galaxy evolution. However, its sensitivity to cosmology and astrophysics at the same time is also a nuisance. Here we develop a comprehensive framework for modelling the LIM power spectrum, which includes redshift space distortions and the Alcock-Paczynski effect. We then identify and isolate degeneracies with astrophysics so that they can be marginalized over. We study the gains of using the multipole expansion of the anisotropic power spectrum, providing an accurate analytic expression for their covariance, and find a 10%-60% increase in the precision of the baryon acoustic oscillation scale measurements when including the hexadecapole in the analysis. We discuss different observational strategies when targeting other cosmological parameters, such as the sum of neutrino masses or primordial non-Gaussianity, finding that fewer and wider bins are typically more optimal. Overall, our formalism facilitates an optimal extraction of cosmological constraints robust to astrophysics.

Tables

Table 1. Table 1: Forecasted 68% confidence-level marginalized relative constraints using our generalized experiment on the BAO parameters and ⟨ T ⟩ ​ f ​ σ 8 delimited-⟨⟩ 𝑇 𝑓 subscript 𝜎 8 \langle T\rangle f\sigma_{8} , expressed as percentages. We compare results for only the monopole and quadrupole to the case of adding the hexadecapole, leading to marked improvement.

$z$	${\sigma }_{\mathrm{rel}}\left({\alpha }_{⟂}\right)(\%)$		${\sigma }_{\mathrm{rel}}\left({\alpha }_{\parallel }\right)(\%)$		${\sigma }_{\mathrm{rel}}\left(⟨T⟩f{\sigma }_8\right)(\%)$
	$\ell \le 2$	$\ell \le 4$	$\ell \le 2$	$\ell \le 4$	$\ell \le 2$	$\ell \le 4$
2.73	1.5	1.1	2.0	1.3	3.4	2.1
4.01	1.4	1.1	2.0	1.4	3.5	2.5
5.30	1.7	1.5	2.4	1.9	3.8	3.0
6.58	3.8	3.4	4.8	4.1	5.9	5.1
7.87	10.6	9.5	12.2	11.0	9.5	8.6

Equations

$$\begin{array}{rl}& T(z)=\frac{c^3(1+z{)}^2}{8\pi k_B{\nu }^{3}H(z)}{\rho }_L(z),\\ & I(z)=\frac{c}{4\pi \nu H(z)}{\rho }_L(z),\end{array}$$

$$⟨{\rho }_L⟩(z)=\int \mathrm{d}ML(M,z)\frac{\mathrm{d}n}{\mathrm{d}M}(M,z).$$

$$b(z)=\frac{\int \mathrm{d}ML(M,z)b_{\mathrm{h}}(M,z)\frac{\mathrm{d}n}{\mathrm{d}M}(M,z)}{\int \mathrm{d}ML(M,z)\frac{\mathrm{d}n}{\mathrm{d}M}(M,z)}.$$

$$F_{\mathrm{RSD}}(k,\mu ,z)=\left(1+\frac{f(z)}{b(z)}{\mu }^2\right)\frac{1}{1+0.5{\left(k\mu {\sigma }_{\mathrm{FoG}}\right)}^2},$$

$$\begin{array}{rl}& P(k,\mu ,z)=P_{\mathrm{clust}}(k,\mu ,z)+P_{\mathrm{shot}}(z);\\ & P_{\mathrm{clust}}(k,\mu ,z)=⟨T{⟩}^2(z)b^2(z)F_{\mathrm{RSD}}^2(k,\mu ,z)P_{\mathrm{m}}(k,z);\\ & P_{\mathrm{shot}}(z)={\left(\frac{c^3(1+z{)}^2}{8\pi k_B{\nu }^{3}H(z)}\right)}^2\int dM{L}^2(M,z)\frac{\mathrm{d}n}{\mathrm{d}M},\end{array}$$

$${\alpha }_{\perp }=\frac{D_A(z)/r_{\mathrm{s}}}{{\left(D_A(z)/r_{\mathrm{s}}\right)}^{\mathrm{fid}}},{\alpha }_{\parallel }=\frac{{\left(H(z)r_{\mathrm{s}}\right)}^{\mathrm{fid}}}{H(z)r_{\mathrm{s}}},$$

$$\begin{array}{rl}& k^{\mathrm{true}}=\frac{k^{\mathrm{meas}}}{{\alpha }_{\perp }}{\left[1+{\left({\mu }^{\mathrm{meas}}\right)}^2\left(F_{\mathrm{AP}}^{-2}-1\right)\right]}^{1/2},\\ & {\mu }^{\mathrm{true}}=\frac{{\mu }^{\mathrm{meas}}}{F_{\mathrm{AP}}}{\left[1+{\left({\mu }^{\mathrm{meas}}\right)}^2\left(F_{\mathrm{AP}}^{-2}-1\right)\right]}^{-1/2},\end{array}$$

$$P(k,\mu )={\left(\frac{⟨T⟩b{\sigma }_8+⟨T⟩f{\sigma }_8{\mu }^2}{1+0.5{\left(k\mu {\sigma }_{\mathrm{FoG}}\right)}^2}\right)}^2\frac{P_m(k,\vec{\varsigma })}{{\sigma }_8^2}+P_{\mathrm{shot}},$$

$$\vec{\theta }=\{{\alpha }_{\perp },{\alpha }_{\parallel },⟨T⟩f{\sigma }_8,⟨T⟩b{\sigma }_8,{\sigma }_{\mathrm{FoG}},P_{\mathrm{shot}},\vec{\varsigma }\}.$$

$$P(k,\mu )={\left(\frac{⟨T⟩b{\sigma }_8\left(1+\beta {\mu }^2\right)}{1+0.5{\left(k\mu {\sigma }_{\mathrm{FoG}}\right)}^2}\right)}^2\frac{P_{\mathrm{m}}(k)}{{\sigma }_8^2}+P_{\mathrm{shot}}.$$

$$\begin{array}{rl}\tilde{P}(k,\mu ,z)& =W(k,\mu ,z)P(k,\mu ,z)=\\ & =W_{\mathrm{vol}}(k,\mu ,z)W_{\mathrm{res}}(k,\mu ,z)P(k,\mu ,z),\end{array}$$

$$\begin{array}{rl}{\sigma }_{\parallel }=& \frac{c(1+z)\delta \nu }{H(z){\nu }_{\mathrm{obs}}};\\ {\sigma }_{\perp }^{\mathrm{dish}}=\chi (z){\sigma }_{\mathrm{beam}},& {\sigma }_{\perp }^{\mathrm{interf}}=\frac{c\chi (z)}{{\nu }_{\mathrm{obs}}{D}_{\max }},\end{array}$$

$$W_{\mathrm{res}}(k,\mu )=\exp \left\{-k^2\left[{\sigma }_{\perp }^2(1-{\mu }^2)+{\sigma }_{\parallel }^2{\mu }^2\right]\right\}.$$

$$V_{\mathrm{field}}=\left[{\chi }^2(z){\Omega }_{\mathrm{field}}\right]\left[\frac{c(1+z{)}^2\Delta \nu }{H(z)\nu }\right].$$

$$\begin{array}{rl}{W}_{\mathrm{vol}}(k,\mu )=& \left(1-\exp \left\{-{\left(\frac{k}{k_{\perp }^{\min }}\right)}^2\left(1-{\mu }^2\right)\right\}\right)\times \\ & \times \left(1-\exp \left\{-{\left(\frac{k}{k_{\parallel }^{\min }}\right)}^2{\mu }^2\right\}\right).\end{array}$$

$$\begin{array}{rl}{\tilde{P}}_{\ell }(k^{\mathrm{meas}})=& \frac{H(z)}{H^{\mathrm{fid}}(z)}{\left(\frac{D_A^{\mathrm{fid}}(z)}{D_A(z)}\right)}^2\frac{2\ell +1}{2}\times \\ & \times {\int }_{-1}^1\mathrm{d}{\mu }^{\mathrm{meas}}\tilde{P}(k^{\mathrm{true}},{\mu }^{\mathrm{true}}){\mathcal{L}}_{\ell }({\mu }^{\mathrm{meas}}),\end{array}$$

$$\begin{array}{rl}{P}_n^{\mathrm{dish}}=& \frac{T_{\mathrm{sys}}^2{V}_{\mathrm{field}}}{\Delta \nu t_{\mathrm{obs}}{N}_{\mathrm{pol}}{N}_{\mathrm{feeds}}{N}_{\mathrm{ant}}};\\ P_n^{\mathrm{interf}}=& \frac{T_{\mathrm{sys}}^2{V}_{\mathrm{field}}{\Omega }_{\mathrm{FOV}}}{\Delta \nu t_{\mathrm{obs}}{N}_{\mathrm{pol}}{N}_{\mathrm{feeds}}{n}_s},\end{array}$$

$$n_s=\frac{2\pi {\nu }_{\mathrm{obs}}^2(D_{\max }^2-D_{\min }^2)}{c^2{N}_{\mathrm{ant}}(N_{\mathrm{ant}}-1)}.$$

$$\begin{array}{rl}& {\tilde{\sigma }}^2(k_i,k_j,\mu )=\frac{1}{\sqrt{N_{\mathrm{modes}}(k_i,\mu )}}\left(\tilde{P}(k_i,\mu )+P_n\right)\times \\ & \times \frac{1}{\sqrt{N_{\mathrm{modes}}(k_j,\mu )}}\left(\tilde{P}(k_j,\mu )+P_n\right){\delta }_{ij}^K,\end{array}$$

$$N_{\mathrm{modes}}(k,\mu )=\frac{k^2\Delta k\Delta \mu }{8{\pi }^2}{V}_{\mathrm{field}},$$

$$\begin{array}{rl}{\tilde{\mathcal{C}}}_{\ell {\ell }^{\prime }}(k_i,k_j)=& \frac{(2\ell +1)(2{\ell }^{\prime }+1)}{2}\times \\ & \times {\int }_{-1}^1\mathrm{d}\mu {\tilde{\sigma }}^2(k_i,k_j,\mu ){\mathcal{L}}_{\ell }(\mu ){\mathcal{L}}_{{\ell }^{\prime }}(\mu ),\end{array}$$

$$\begin{array}{rl}\left[S/N\right(k){]}^2& ={\vec{\Theta }}^T(k){\tilde{\mathcal{C}}}^{-1}(k)\vec{\Theta }(k),\\ {\left[S/N\right]}^2& ={\vec{\Theta }}^T{\tilde{\mathcal{C}}}^{-1}\vec{\Theta },\end{array}$$

$${\chi }^2=\Delta {\vec{\Theta }}^T{\tilde{\mathcal{C}}}^{-1}\Delta \vec{\Theta },$$

$$F_{{\vartheta }_a{\vartheta }_b}=⟨\frac{{\partial }^2\log \mathcal{L}}{\partial {\vartheta }_a\partial {\vartheta }_b}⟩=\left(\frac{\partial {\vec{\Theta }}^T}{\partial {\vartheta }_a}{\tilde{\mathcal{C}}}^{-1}\frac{\partial \vec{\Theta }}{\partial {\vartheta }_b}\right).$$

$$\Phi (x)=\varphi (x)+f_{\mathrm{NL}}\left({\varphi }^2(x)-⟨{\varphi }^2⟩\right).$$

$$\Delta b_h(k,z)=\left[b_h^{\mathrm{G}}(z)-1\right]f_{\mathrm{NL}}{\delta }_{\mathrm{ec}}\frac{3{\Omega }_m{H}_0^2}{c^2{k}^2{T}_m(k)D(z)}.$$

$${\sigma }_{{\vartheta }_{\varsigma }}={\left\{\sum _z{\left[{\left(F_z^{-1}\right)}_{{\vartheta }_{\varsigma }{\vartheta }_{\varsigma }}\right]}^{-1}\right\}}^{-1/2},$$

$$\begin{array}{rl}& P^{XY}=P_{\mathrm{clust}}^{XY}+P_{\mathrm{shot}}^{XY};\\ & P_{\mathrm{clust}}^{XY}=⟨T_X⟩⟨T_Y⟩b_X{b}_Y{F}_{\mathrm{RSD}}^X{F}_{\mathrm{RSD}}^Y{P}_{\mathrm{m}};\\ & P_{\mathrm{shot}}^{XY}={\left(\frac{c^3(1+z{)}^2}{8\pi k_B{\nu }^{3}H(z)}\right)}^2\int dM{L}_X(M)L_Y(M)\frac{\mathrm{d}n}{\mathrm{d}M},\end{array}$$

$$\begin{array}{rl}{P}^{XY}& =\frac{P_{\mathrm{m}}/{\sigma }_8^2}{{\left(1+0.5{\left[k\mu {\sigma }_{\mathrm{FoG}}\right]}^2\right)}^2}\times \\ & \times \left({\tilde{T}}_{XY}^{1/2}{b}_X{\sigma }_8+{\tilde{T}}_{XY}^{1/2}f{\sigma }_8{\mu }^2\right)\times \\ & \times \left({\tilde{T}}_{XY}^{1/2}{b}_Y{\sigma }_8+{\tilde{T}}_{XY}^{1/2}f{\sigma }_8{\mu }^2\right)+P_{\mathrm{shot}}^{XY},\end{array}$$

$$\begin{array}{rl}{\vec{\theta }}_{XY}=\{& {\alpha }_{\perp },{\alpha }_{\parallel },{\tilde{T}}_{XY}^{1/2}f{\sigma }_8,{\tilde{T}}_{XY}^{1/2}{b}_X{\sigma }_8,{\tilde{T}}_{XY}^{1/2}{b}_Y{\sigma }_8,\\ & {\sigma }_{\mathrm{FoG}},P_{\mathrm{shot}}^{XY},\vec{\varsigma }\}.\end{array}$$

Full text

I Introduction

Line-intensity mapping (LIM) Kovetz_IMstatus has recently arisen as a key technique to surpass and expand the remarkable achievements of precision observational cosmology over the past decades. Impressive experimental efforts concentrated mostly on two observables, namely the cosmic microwave background (CMB) and galaxy number counts, have attained percent-level measurements of the standard cosmological model, $\Lambda$-Cold Dark Matter ($\Lambda$CDM), and provided stringent constraints on possible deviations from it (Planck18_pars, ; alam_bossdr12, ). However, the former only allows limited (due to damping on small scales) access to a brief moment in the early history of the Universe, while the latter loses effectiveness as it probes deeper into the Universe (as the discrete sources become too faint to detect).

LIM provides a way to bridge the gap between the CMB and galaxy surveys, and probe the huge swaths of the observable Universe that remain uncharted Karkare:2018sar ; Kovetz_IM2020 . It measures the integrated emission from either atomic or molecular spectral lines originating from all the galaxies—individually detectable or not—as well as from the diffuse intergalactic medium along the line of sight, and generates three-dimensional maps of the targeted volume that in principle contain all the information that can be harvested from the incoming photons.

LIM’s experimental landscape is rapidly progressing. Probably the most studied line to date is the 21-cm spin-flip transition in neutral hydrogen, first detected by Ref. Chang_21cmdetection . A comprehensive suite of experiments is on its way to map this line across the history of the observable Universe, with several experiments targeting the epochs of cosmic dawn and reionization (Haarlem_lofar, ; Tingay_MWA, ; hera, ; Koopmans_SKAlow, ) and others focusing on lower redshifts (Masui_GBT, ; CHIME_pathfinder, ; HIRAX, ; Santos_Meerkat, ; redbook, ; Battye_bingo, ). Meanwhile, growing attention has been given to other lines. Some examples include: carbon monoxide (CO) rotational lines (Lidz2011, ; Pullen2013, ; Breysse2014, ; Li_CO_16, ; Padmanhaban_CO, ), [CII] (Silva_CII, ; Pullen_CII, ; Padmanabhan_CII, ), H$\alpha$ and H$\beta$ (Gong_lines, ; Silva_Halpha, ), oxygen lines (Gong_lines, ) and Lyman-$\alpha$ (Pullen2014, ; Silva_Lyalpha, ). A few of these have already been (at least tentatively) detected at intermediate redshifts (Keating_COPSS, ; Keating_COPSSII, ; Pullen_CII, ; Yang2019, ; Croft_Lyalphadetection, ). A significant effort is now being invested in LIM experiments targeting these lines, with some instruments already observing and others to come online soon Cleary_COMAP ; TIME ; Hetdex ; spherex ; Concerto ; CCAT-prime ; Starfire ; aim-co ; Cooray_CDIM ; OST_origins .

LIM holds unique promise for the study of both cosmology and astrophysics (Kovetz_IMstatus, ; Kovetz_IM2020, ). Since the line emission originates from halos, the intensity of the spectral lines acts as a biased tracer of the underlying density distribution (sourced by the cosmological primordial fluctuations), with the bias depending on the specific line. Meanwhile, the signal is intimately related to various astrophysical processes that take place during reionization and galaxy evolution (Breysse2016, ; Serra2016, ; Breysse2017, ; Wolz2017, ; Lidz2018, ; Padmanabhan2018, ; Breysse2019, ). A primary challenge is therefore to disentangle between the astrophysical and cosmological information contained in the maps.

Previous works have investigated the potential of using LIM observations to study cosmology within the confines of $\Lambda$CDM and beyond it (see e.g., Refs. Bull_21cm ; Raccanelli_ska ; Carucci_wdm21cm ; Pourtsidou_21cmcosmo ; Karkare_IMBAO ; Fonseca_multitracerlines ; Obuljen_21cmcrossothers ; Munoz_chargedDM21cm ; Chen_hiddenvalley ; Moradinezhad_Dizga_COfNL1 ; Moradinezhad_Dizga_COfNL2 ; redbook ; SKA_fundamental ; Castorina_21cmgrowth ; Foreman_lensing21cm ; Schaan_IMlensing ; Jalilvand_21cmgallensing ). However, existing studies are limited in the sense that they do not fully account for or properly model one or more of the following: the observable signal, observing limitations, intrument response, degeneracies between astrophysics and cosmology, or the noise covariance. Ignoring these effects may result in a significant underestimation of the errors and bias of the best-fit parameters even when dealing merely with forecasts (see e.g., Refs. Bellomo_Fisher ; Bernal_Fisher for some examples related to galaxy surveys).

In this paper we aim to address this deficit by providing a general and comprehensive formalism to optimally extract robust cosmological information from the power spectrum of fluctuations in observed line-intensity maps. We focus especially on maximizing the precision of baryon acoustic oscillation (BAO) measurements, but also consider the potential of LIM to constrain extensions to $\Lambda$CDM, such as primordial non-Gaussianity or neutrinos with a sum of masses higher than $0.06\,{\rm eV}$. In particular, we consistently account for the following:

• •

Degeneracies between cosmology and astrophysics. We identify and isolate the cosmological information encoded in the power spectrum that can be extracted given known astrophysical dependencies.

• •

LIM power spectrum anisotropies. We present a model for the power spectrum which includes the anisotropic information from redshift-space distortions and the Alcock-Paczynski effect (see below).

• •

Multipole decomposition of the power spectrum. We show the optimal way to account for these anisotropies, including the full covariance between different multipoles. We demonstrate generically the importance of going to higher multipoles.

• •

Effects of survey volume and instrument response. We show how to accurately model the suppression of the power spectrum on large and small scales due to survey and beam shapes. We discuss the effects of redshift binning and experimental optimization.

In order to remain as generic as possible in our quantitative estimates, we refrain from forecasting for specific experiments and consider a generic LIM instrument which can easily be replaced with any concrete experimental setup. This manuscript (hopefully) results in a useful manual for LIM enthusiasts who wish to precisely quantify the cosmological information accessible in their observed maps.

To perform the analyses this work is based on, we modify the public code LIM111https://github.com/pcbreysse/lim. We analytically compute the anisotropic LIM power spectrum and the corresponding covariance, using outputs from CAMB222https://camb.info/ (CAMB, ) to extract the matter power spectrum and other cosmological quantities, and Pylians333https://github.com/franciscovillaescusa/Pylians (VillaescusaNavarro_pylians, ) to obtain the halo mass function and halo bias. We plan to release an updated version of the LIM package upon publication.

Our calculations are based on several standard assumptions and parameter choices. Throughout this work we adopt the halo mass function and halo bias presented in Ref. Tinker_hmf2010 . We assume, unless otherwise stated, the Planck baseline $\Lambda$CDM model and take the best-fit parameter values from the combined analysis of the full Planck 2018 and galaxy BAO measurements (Planck18_pars, ; alam_bossdr12, ): baryon and dark matter physical densities at $z=0$, $\Omega_{\rm b}h^{2}=0.0224$ and $\Omega_{\rm cdm}h^{2}=0.1193$, respectively; spectral index $n_{s}=0.967$; amplitude $\log\left(A_{s}10^{10}\right)=3.047$ of the primordial power spectrum of scalar perturbations at the pivot scale; Hubble constant $H_{0}=67.67$ km/s/Mpc; a sum of neutrino masses $\sum m_{\nu}=0.06$ eV; and Gaussian primordial perturbations, i.e., $f_{\rm NL}=0$.

This paper is structured as follows. In Section II we review the basics of LIM power spectrum modelling. We highlight the degeneracies between astrophysics and cosmology and provide a formulation which clearly delineates them. Properties of the measured power spectrum are derived in Section III, including the Legendre multipole expansion, and the effects of the window function and the covariance matrix. A generalized LIM experiment is introduced in Section IV, along with the prescription for our Fisher-based forecasts. In Section V we investigate the extraction of cosmological information from the LIM power spectrum, demonstrating the importance of higher multipoles when measuring BAO and RSD, and describing strategies to improve cosmological constraints, including redshift binning and experimental optimization. We defer to an Appendix the adaptation of our methodology to cross-correlations between different lines or with galaxy surveys (Appendix A), as well as the discussion of how to fold in the instrument response and the signal suppression due to the finite volume surveyed, and the effects of foregrounds and line interlopers (Appendix B). We present our conclusions in Section VI.

II Modelling the IM Power Spectrum

II.1 From line luminosity to the power spectrum

The brightness temperature $T$ and specific intensity $I$ of a given radiation source are equivalent quantities that depend on the expected luminosity density $\rho_{L}$ of a spectral line with rest frame frequency $\nu$ at redshift $z$:

[TABLE]

where $c$ is the speed of light, $k_{B}$ is the Boltzmann constant and $H$ is the Hubble parameter Lidz2011 . Depending on the frequency band of the experiment, either $T$ or $I$ are conventionally used. Hereinafter, we will use $T$ in order to homogenize the nomenclature, but all the expressions are equally valid if intensity is used instead. Assuming a known relation, $L(M,z)$, between the luminosity of the spectral line and the mass $M$ of the host halo, the expected mean luminosity density can be computed using the halo mass function ${\rm d}n/{\rm d}M(z)$:

[TABLE]

The main statistic we will focus on in this work is the LIM power spectrum, which is given by the Fourier transform of the two-point correlation function of the perturbations of the brightness temperature, denoted by $\delta T\equiv T-\langle T\rangle$. Since spectral lines are sourced in halos, $\delta T$ can be used to trace the halo distribution and, equivalently, the underlying matter density perturbations. At linear order, matter density and brightness temperature perturbations are related by an effective linear bias, which is given in terms of the halo bias $b_{\rm h}$ by

[TABLE]

In practice, $T$ maps are obtained in redshift space, where the observed position along the line of sight is disturbed with respect to real space due to peculiar velocities, producing the so-called redshift-space distortions (RSD). RSD introduce anisotropies in the a-priori isotropic real-space power spectrum. The relation between real and redshift space perturbations ($\delta^{r}$ and $\delta^{s}$, respectively) can be approximated by a linear factor, such that $\delta^{s}=F_{\rm RSD}\delta^{r}$, with $F_{\rm RSD}$ in Fourier space given by:

[TABLE]

where $f(z)$ is the logarithmic derivative of the growth factor (also known as the growth rate), $k$ is the module of the Fourier mode, and $\mu=\hat{k}\cdot\hat{k}_{\parallel}$ is the cosine of the angle between the mode vector $\vec{k}$ and the line-of-sight component, $k_{\parallel}$ (i.e. $\mu\in\left[-1,1\right]$). At linear scales, coherent peculiar velocities boost the power spectrum in redshift space through the Kaiser effect (first term in Eq. (4)) Kaiser1987 . In turn, small-scale velocities suppress the clustering on small scales, an effect known as the fingers of God (second term in Eq. (4)). We use a Lorentzian damping factor whose scale-dependence is driven by the parameter $\sigma_{\rm FoG}$, the value of which is related to the halo velocity dispersion444The fingers-of-God damping can be also modelled with a Gaussian function, yielding similar results. We have checked that the conclusions of this work do not depend on this choice.. We assume a fiducial value of $7\,{\rm Mpc}$. We refer the interested reader to e.g., Ref. Hamilton_RSDreview , for a review on RSD discussing these contributions.

On top of all these effects, there is a scale-independent shot noise contribution to the LIM power spectrum stemming from the discreteness of the source population, which would be present even in the absence of clustering. Thus we can express the anisotropic LIM power spectrum as the sum of clustering and shot noise contributions:

[TABLE]

where $P_{\rm m}$ is the matter power spectrum.

We note that the above expressions assume that each halo contains only a single point-source line emitter at its center. In many cases, there will be wider-scale emission from diffuse gas within a halo, or there will be satellite galaxies which contribute additional intensity to that of the central emitter. Either of these will lead to an additional “one-halo” contribution to the clustering term which traces the halo density profile. We direct the interested reader to Ref. Wolz2019 for a discussion of one-halo effects in intensity mapping surveys. We also leave for future work the effect of anisotropic halo assembly bias Obuljen:2019ukz .

II.2 Dependence on astrophysics and cosmology

The emission of any spectral line is intimately related to astrophysical processes. Therefore, the observed density field is strongly dependent on astrophysics through the line luminosity function. Fortunately, according to Eq. (5), the relation between astrophysics and the LIM power spectrum at linear order is limited to the amplitude of the clustering and shot noise contributions555A detailed study of the nonlinearity and scale dependence of the astrophysical terms in the LIM power spectrum lies beyond the scope of this work.. This means that the only accessible astrophysical information present in the linear LIM power spectrum is contained in the first two moments of $L(M)$666Alternative probes such as the voxel intensity distribution Breysse_VID are more suitable to infer $L(M)$, see Section VI..

Since $\delta T$ traces matter overdensities, LIM observations also carry cosmological information. This information is mostly encoded in $P_{\rm m}$, $F_{\rm RSD}$ and the BAO. For example, as seen in Eq. (4), $F_{\rm RSD}$ depends on the growth factor $f$, and so can be efficiently used to constrain deviations from general relativity.

In order to measure the power spectrum, redshifts need to be transformed into distances. This procedure introduces further anisotropies in the power spectrum if the assumed cosmology does not match the actual one. This is known as the Alcock-Paczynski effect Alcock_paczynski . Radial and transverse distances are distorted in different ways, which can be modelled by introducing rescaling parameters to redefine distances:

[TABLE]

where $D_{A}$ is the angular diameter distance, $r_{\rm s}$ is the sound horizon at radiation drag, and the superscript ‘${\rm fid}$’ denotes the corresponding values in the assumed (fiducial) cosmology. Due to this distortion, the true wave numbers are related to the measured ones in the transverse and line-of-sight directions by $k^{\rm true}_{\perp}=k^{\rm meas}_{\perp}/\alpha_{\perp}$ and $k^{\rm true}_{\parallel}=k^{\rm meas}_{\parallel}/\alpha_{\parallel}$, respectively. These relations can then be expressed in terms of $k$ and $\mu$ (Ballinger_scaling96, ):

[TABLE]

where $F_{\rm AP}\!\equiv\!\alpha_{\parallel}/\alpha_{\perp}$. To correct for the modification of the volumes, the power spectrum is then multiplied by the factor: $H(z)/H^{\rm fid}(z)\times\left(D_{A}^{\rm fid}(z)/D_{A}(z)\right)^{2}$. BAO provide useful cosmic rulers, as they allow high precision measurements of $\alpha_{\perp}$ and $\alpha_{\parallel}$, which can effectively constrain the expansion history of the Universe, see Ref. Bernal_IM_letter . Moreover, the Alcock-Paczynski effect has been proposed as a method to identify and remove line interlopers Lidz_APforeground .

Other cosmological dependences in Eq. (5) arise through the halo bias and the halo mass function, and also due to the halo velocity dispersion (via $\sigma_{\rm FoG}$). However, the cosmological information encoded in these quantities will be degenerate with the highly-uncertain $L(M)$ modeling. Therefore, in our analysis below we will mostly focus on BAO and RSD.

II.3 Explicitly accounting for the degeneracies

As clearly evident, the quantities appearing in Eq. (5) present strong degeneracies. Therefore, after identifying where the cosmological information is encoded, it is useful to reparametrize Eq. (5) to minimize the degeneracies.

The overall shape of $P_{\rm m}$ is not very sensitive to small variations of the cosmological parameters of $\Lambda$CDM. Thus, it is preferable to employ a template for $P_{\rm m}(k,z)$, computed based on the fiducial cosmology (parametrizing its amplitude with $\sigma_{8}$, the root mean square of the density fluctuations within $8\,h{\rm Mpc}^{-1}$), and then measure the terms that relate it to the LIM power spectrum.

Examining Eq. (5), we find that $\sigma_{8}$, $b$ and $\langle T\rangle$ are degenerate within a given single redshift, as well as $\sigma_{8}$, $f\mu^{2}$, and $\langle T\rangle$. These degeneracies can be broken with an external prior on $\langle T\rangle$, or using tomography and interpreting the measurements under a given cosmological model, taking advantage of their different redshift evolutions.

Other ways to break these degeneracies, at least partially, include the combination of the LIM power spectrum and higher-order statistics (see e.g., Ref. GilMarin_bispectrunNbody in the case of halo number counts, and Ref. Sarkar_21cmBispec for the 21cm LIM case); multi-tracer techniques with several spectral lines and galaxy surveys Obuljen_21cmcrossothers ; or the exploitation of the mildly non-linear regime of the LIM power spectrum by modelling the halo clustering with Lagrangian perturbation theory (Castorina_21cmgrowth, ).

Taking all this into account, we group all degenerate parameters and reparameterize the LIM power spectrum from Eq. (5) as:

[TABLE]

where all quantities depend on $z$. From Eq. (8), the set of parameter combinations that can be directly measured from the LIM power spectrum at each independent redshift bin and observed patch of sky is:

[TABLE]

In $\vec{\varsigma}$ we group all parameters that modify the template of the power spectrum used in the analysis, such as any of the standard $\Lambda$CDM parameters, or extensions like primordial non-Gaussianity or neutrino masses.

As explained above, out of all the parameter combinations included in Eq. (9), only $\alpha_{\perp}$, $\alpha_{\parallel}$ and $\langle T\rangle f\sigma_{8}$ contain useful cosmological information (modulo $\langle T\rangle$, in the case of $\langle T\rangle f\sigma_{8}$). Given the difficulty in tracing cosmological information hidden in the other parameter combinations, we prefer to be conservative and consider them simply as nuisance parameters in our analysis.

Lastly, note that we choose this parameterization in order to separate $b$ and $f$. However, the LIM power spectrum can be cast in different ways in order to target other parameter combinations. For example, if a measurement of $\beta=f/b$ is wanted, one could use:

[TABLE]

In this case, the measured quantities would be $\beta$ and $\langle T\rangle b\sigma_{8}$, instead of $\langle T\rangle f\sigma_{8}$ and $\langle T\rangle b\sigma_{8}$.

III Measuring the IM Power Spectrum

III.1 The window function and multipole expansion

LIM experiments have a limited resolution and probe a finite volume, which limits the minimum and maximum accessible scales, respectively. This effectively renders the observed brightness temperature fluctuations smoothed with respect to the true ones. This smoothing can be modelled by convolving the true $T$ map with window functions describing the instrument response and the effects due to surveying a finite volume. This would yield an observer-space power spectrum, given by:

[TABLE]

where $W_{\rm vol}$ and $W_{\rm res}$ respectively are the survey-area and instrument-response window functions in Fourier space.

The line-of-sight resolution of a LIM instrument is given by the width of the frequency channels. The resolution in the plane of the sky depends on whether the experiment uses only the auto-correlation of each of its antennas (i.e., a single-dish approach) or employs interferometric techniques. In the former case, the resolution is determined by the antenna’s beam profile, while in the latter, by the largest baseline of the interferometer, $D_{\rm max}$. The characteristic resolution limits in the radial and transverse directions are given by:

[TABLE]

where $\delta\nu$ is the frequency-channel width, $\nu_{\rm obs}$ is the observed frequency, $\chi(z)$ is the radial comoving distance, and the width of the beam profile is given by $\sigma_{\rm beam}=\theta_{\rm FWHM}/\sqrt{8\log 2}$ (where $\theta_{\rm FWHM}$ is its full width at half maximum). Then, the resolution window function $W_{\rm res}$ that models the instrument response in Fourier space can be computed as (Li_CO_16, ):

[TABLE]

The beam profile, field of view, and frequency band window might not be Gaussian; for instance, the frequency channels are likely to be discrete bins (which would correspond to a top-hat window in real space). In practice, the exact instrument response will be accurately characterized by each experiment. Investigation of more realistic window functions is beyond the scope of this work.

On the other hand, at a given redshift, observing a patch of the sky with solid angle $\Omega_{\rm field}$ over a frequency band $\Delta\nu$ corresponds to a surveyed volume $V_{\rm field}$, which in the absence of complex observation masks is given by:

[TABLE]

The two factors in Eq. (14) are the transverse area and the length of the radial side of the volume observed (i.e., $V_{\rm field}\sim L_{\perp}^{2}L_{\parallel}$). These are the largest scales that can be probed by a single-dish like experiment, hence we define $k_{\parallel}^{\rm min,dish}\equiv 2\pi/L_{\parallel}$ and $k_{\perp}^{\rm min,dish}\equiv 2\pi/L_{\perp}$. While an interferometer shares this limitation in the radial direction, the largest scales that an interferometer can measure are determined by the shortest baseline $D_{\rm min}$, so that $k_{\perp}^{\rm min,interf}\equiv 2\pi\nu_{\rm obs}D_{\rm min}/(c\chi(z))$. Note that $k_{\perp}^{\rm min,interf}$ will always be ultimately limited by $L_{\perp}$. LIM power spectrum measurements beyond these scales will be suppressed by the lack of observable modes.

To account for the loss of modes, we define a volume window function $W_{\rm vol}$, which in Fourier space is given by:

[TABLE]

We use a smoothed window in order to extend the loss of modes to smaller scales than those corresponding to $k_{\perp,\parallel}^{\rm min}$. This mimics the effect of having residuals present in the map after the removal of foregrounds (and line-interlopers), polluting the signal, since large scales are the ones more likely to be affected777One could assume a top-hat window in real space to perform a more aggressive analysis. Complex observational masks may entail more complicated $W_{\rm vol}$, but these effects are very dependent on the particular circumstances of a given experimental setup..

Though we have worked up to this point in $(k,\mu)$ coordinates, we must point out that while the Fourier transform of the observed map is performed using Cartesian coordinates, the line of sight changes with different pointings on the sky. It is therefore not parallel to any Cartesian axis. Therefore, it is not possible to obtain a well defined $\mu$ from the observations, meaning that we cannot directly measure $\tilde{P}(k,\mu)$. Nevertheless, it is possible to directly measure the multipoles of the anisotropic power spectrum using, e.g., the Yamamoto estimator YamamotoEstimator .

Accounting for all the effects described in this section, the multipole expansion of the observed LIM power spectrum is given by:

[TABLE]

where $\mathcal{L}_{\ell}$ is the Legendre polynomial of degree $\ell$. In our analysis below, we will calculate the integral in Eq. (16) numerically. Analytic results (neglecting the fingers-of-God contribution to the RSD and not considering $W_{\rm vol}$) can be found in Ref. Chung_IMPkvoxel . We note that there is a non-vanishing contribution from the shot noise of the LIM power spectrum to multipoles higher than the monopole, due to the effect of the window function.

Most previous works do account for the smearing of the brightness temperature maps, but choose to include this effect in the covariance of the LIM power spectrum rather than in the observed signal. This is equivalent to deconvolving the observed LIM power spectrum by $W(k,\mu)$. As we show in Appendix B, this approach may result in a reduction of the significance of the measured LIM power spectrum multipoles when their covariance is properly modelled888The contribution to the covariance of the monopole is often modelled incorrectly, see Appendix B., unless special treatment is applied to both signal and errors. We explore this avenue and derive a nearly optimal estimator of the multipoles of the true unsmoothed LIM power spectrum. However, given the complexity of applying this approach to observations, we advocate for modelling the smoothed LIM power spectrum rather than deconvolving the observed one.

III.2 The power spectrum covariance matrix

There are three contributions to the LIM power spectrum covariance: (i) sample variance; (ii) instrumental noise; and (iii) residual contamination after removal of foregrounds and line interlopers. Since the magnitude of the effect of foregrounds and line interlopers greatly depends on the spectral line and redshift of interest, we consider only the first two contributions to the covariance and assume that the signal has been cleaned such that the effects included in $W_{\rm vol}$ are sufficient. We refer the reader interested in foreground and line-interloper removal strategies to, e.g., Refs. Oh2003 ; Wang2006 ; Liu2011 ; Lidz_APforeground ; Cheng_foregroundsAP ; Sun_foregrounds ; Breysse_foregrounds . We also assume a Gaussian covariance without mode coupling.

The instrumental noise, set by the practical limitations of the experiment, introduces an artificial floor brightness temperature in the observed map. The noise power spectrum for a given single-dish like experiment or an interferometer is given by (see e.g., Ref. Bull_21cm ):

[TABLE]

where $t_{\rm obs}$ is the total observing time, $N_{\rm feeds}$ is the number of detectors in each antenna, each of them able to measure $N_{\rm pol}=1,2$ polarizations, $T_{\rm sys}$ is the system temperature, $\Omega_{\rm FOV}=c^{2}/(\nu_{\rm obs}D_{\rm dish})^{2}$ is the field of view of an antenna, and $n_{s}$ is the average number density of baselines in the visibility space. Assuming a constant number density of baselines, the latter is given by Bull_21cm :

[TABLE]

Although the anisotropic power spectrum may not be measurable, it is useful to define the covariance per $k$ and $\mu$ bin:

[TABLE]

where $\delta^{K}$ is the Kronecker delta (introduced because we neglect mode coupling). In Eq. (19), the first term in the parentheses corresponds to the sample variance, and $N_{\rm modes}$ denotes the number of modes per bin in $k$ and $\mu$ in the observed field:

[TABLE]

with $\Delta k$ and $\Delta\mu$ referring to the width of the $k$ and $\mu$ bins, respectively.

The covariance matrix of the LIM power spectrum multipoles is comprised of the sub-covariance matrices of each multipole, and those between different multipoles. The sub-covariance matrix for multipoles $\ell$ and $\ell^{\prime}$ is given by:

[TABLE]

where the Delta function we assumed in Eq. (19) makes each sub-covariance matrix diagonal. However, modeling the non-zero covariance between multipoles is essentiall for an unbiased analysis. We refer the interested reader to Ref. Grieb_covmat for a thorough analytic derivation of the covariance of the multipoles of the galaxy power spectrum under the Gaussian assumption.

IV Forecasting for IM Experiments

IV.1 Parameterizing a generic IM experiment

In order to illustrate the potential of our methodology, it is useful to consider a concrete example for a LIM experiment — without loss of generality — so that we can compute the window functions and the covariance. Without specifying the targeted emission line, we choose as our strawperson LIM experiment an ambitious single dish-like instrument with a total frequency band $\Delta\nu/\nu~{}\approx~{}0.2$. We split the corresponding volume into several redshift bins to study how the signal — and the extraction of cosmological information — changes with redshift. In our analysis below we use five non-overlapping, independent redshift bins centered at $z=\{2.7,4.0,5.3,6.6,7.9\}$, such that $\log_{10}\left[\Delta(1+z)\right]=\log_{10}\left[\Delta(\nu/\nu_{\rm obs})\right]=0.1$.

In order to compute the signal, we need to assume some numbers for the astrophysical quantities that are present in Eq. (5). Strictly as an example, we will assume a CO intensity mapping signal described by the model from Ref. Li_CO_16 . In the redshift bins defined above, this model gives mean brightness temperature $\langle T\rangle=\{2.9,3.3,3.5,3.1,2.4\}$, luminosity-averaged bias $b=\{1.5,2.1,2.6,3.1,4.0\}$, and shot noise contribution to the LIM power spectrum $P_{\rm shot}=\{354,331,214,114,59\}$ (Mpc/$h$)${}^{3}\mu$K2. Note that we are neglecting the evolution of these quantities within each bin, simply using the value at the central redshift each time.

As an analogy with the futuristic experiment envisioned in Ref. Bernal_IM_letter , we consider an array of single dish antennas with total $N_{\rm pol}N_{\rm feeds}N_{\rm ant}t_{\rm obs}/T_{\rm sys}^{2}=10500$ h/K2 and $18900$ h/K2 for the first two redshift bins, respectively, and $25\times 10^{3}$ h/K2 for the last three redshift bins. The value of $N_{\rm pol}N_{\rm feeds}N_{\rm ant}t_{\rm obs}/T_{\rm sys}^{2}$ changes with redshift because we assume that $T_{\rm sys}$ depends on $\nu_{\rm obs}$ until it saturates: $T_{\rm sys}={\rm max}\left[20,\,\nu_{\rm obs}\left({\rm K}/{\rm GHz}\right)\right]$ (see e.g., Prestage_GBT ; Murphy_nextVLA ). Furthermore, we assume a spectral resolution $\nu_{\rm obs}/\delta\nu=\{15450,11500,9150,7600,6500\}$, and $\theta_{\rm FWHM}=4$ arcmin. Taking into account these experiment specifications, we choose $\Omega_{\rm field}=1000$ deg2 to maximize the significance of the measurement around the scales of the BAO. For this experiment, we would have $\sigma_{\perp}<\sigma_{\parallel}$ and $L_{\perp}>L_{\parallel}$.

We emphasize that our choice of this particular CO model is only an example. Our formalism is general and applicable to any line and experiment, and the quantities we compute here can easily be calculated for any other model the reader might have in mind.

In Fig. 1 we demonstrate the effect of the LIM power spectrum smoothing, discussed in Section III, for the specific case of the first redshift bin of our strawperson experiment, comparing $P_{\ell}$ and $\tilde{P}_{\ell}$. While for the monopole, $W$ just suppresses the power spectrum at high and low $k$, the effect on the quadrupole, and especially on the hexadecapole, is much more significant. The effect largely depends on the ratio between $\sigma_{\perp}$ and $\sigma_{\parallel}$, and between $k_{\perp}^{\rm min}$ and $k_{\parallel}^{\rm min}$. Therefore, the effect of $W$ on the power spectrum is of course strongly dependent on the experiment. This was similarly demonstrated in Ref. Chung_IMPkvoxel , though the survey-volume window effect was not included there.

We also show how the LIM power spectrum multipoles change when each of the parameters of Eq. (9) is $5\%$ larger than the fiducial values. We can see that the derivatives of the quadrupole and hexadecapole with respect to the parameters considerably change whether we consider $P$ or $\tilde{P}$. This figure provides an intuition of the dependence of the power spectrum on each of the considered free parameters, both for the true and the measured LIM power spectra, and may guide the design of future experiments depending on the parameter targeted.

IV.2 Likelihood and Fisher Matrix

If we construct the vectors $\vec{\Theta}(k)=\left[\tilde{P}_{0}(k),\tilde{P}_{2}(k),...\right]$, and $\vec{\Theta}=\left[\tilde{P}_{0}(k_{0}),\tilde{P}_{0}(k_{1}),...,\tilde{P}_{2}(k_{0}),\tilde{P}_{2}(k_{1}),...\right]$, we can compute S/N$(k)$, the total signal-to-noise ratio per $k$ bin, as well as the total S/N summed over all bins as:

[TABLE]

where the superscript T denotes the transpose operator and $\tilde{\mathcal{C}}(k_{i})$ is a $N_{\ell}\times N_{\ell}$ matrix (with $N_{\ell}$ being the number of multipoles included in the analysis) made up of the corresponding elements of $\mathcal{C}_{\ell\ell^{\prime}}(k_{i},k_{i})$. Similarly, we can compute the $\chi^{2}$ of the LIM power spectrum multipoles for a given experiment, using:

[TABLE]

where $\Delta\vec{\Theta}$ is the difference between the model prediction and the actual measurements.

We will use the Fisher matrix formalism to forecast constraints from our generalized LIM experiment (Fisher:1935, ; Tegmark_fisher97, ). The Fisher matrix is the average of the second partial derivatives of the logarithm of the likelihood, $\log\mathcal{L}$, around the best fit (or assumed fiducial model). Eq. (23) can be adapted to form the Fisher matrix by replacing $\Delta\vec{\Theta}$ by the corresponding derivatives. The Fisher matrix element corresponding to the parameters $\vartheta_{a}$ and $\vartheta_{b}$ is:

[TABLE]

V Extracting Cosmology from LIM

Based on the previous sections, we are now equipped to calculate quantitative estimates of the potential of LIM experiments to set robust cosmological constraints. This is a good way to evaluate the methodology presented above. We will use our generic LIM experiment and forecast constraints on $\alpha_{\perp}$, $\alpha_{\parallel}$ and $\langle T\rangle f\sigma_{8}$, as well as on the sum of neutrino masses and primordial non-Gaussianity.

V.1 The importance of going to higher multipoles

In order to exploit the anisotropic BAO, it is necessary to measure at least the monopole and quadrupole of the LIM power spectrum, Eq. (16). This is equivalent to constraining $D_{A}/r_{\rm s}$ and $Hr_{\rm s}$, rather than a combination of them if only the isotropic signal is measured, with the obvious benefits that it entails Padmanabhan_anisoBAO . However, while measuring the monopole and quadrupole is relatively easy, given their reasonably high S/N, detecting the hexadecapole is challenging, since the signal is typically well below the noise. The situation is even more pessimistic for $\ell>4$.

Neglecting the fingers-of-God effect and $W_{\rm vol}$, Ref. Chung_IMPkvoxel finds that including the hexadecapole ($\ell=4$) does not significantly improve the estimated constraints on the astrophysical parameters. Nevertheless, focusing on cosmology and taking into account that the Alcock-Paczynski effect is anisotropic, we find that including the hexadecapole might be useful. Even if the hexadecapole is not detectable, it might help to break degeneracies between the BAO and RSD parameters, since its amplitude is mostly independent of the luminosity-averaged bias, but it does depend greatly on the $\alpha_{\perp}$ and $\alpha_{\parallel}$ configuration. For instance, including the hexadecapole in the analysis of the quasar power spectrum in eBOSS resulted in smaller degeneracies between $\alpha_{\parallel}$ and $\alpha_{\perp}$, and between $\alpha_{\parallel}$ and $f\sigma_{8}$, compared with the case without including the hexadecapole (GilMarin_qsoeBOSS, ).

Moreover, the hexadecapole can be measured and estimated at the same time as the lower order multipoles, without requiring additional computing time. However, robust measurements of the hexadecapole impose stronger requirements on experiments, given its anisotropy. In order to assess the benefits of adding the hexadecapole to the LIM power spectrum analysis, we compare the projected results with and without its inclusion.

We first evaluate the straightforward gain of including more multipoles than just the monopole: we compare the S/N$(k)$ of the LIM power spectrum using only the monopole and then consecutively adding the quadrupole and hexadecapole. The results given our strawperson experiment in each of the five redshift bins are shown in Fig. 2. We can see that adding the quadrupole or the hexadecapole does not significantly increase the total S/N per $k$ bin. Note that this does not mean that the S/N of $\tilde{P}_{2}$ or $\tilde{P}_{4}$ is low, since the total S/N is not the sum of the individual significances of each multipole; the covariance between multipoles must be included (see Eq. (22)).

It is also evident that the signal-to-noise decreases with redshift. This is because the amplitude of the LIM power spectrum decreases with redshift faster than the noise level over the range considered here. In addition, the window function suppresses the LIM power spectrum at wider $k$ ranges. While the behavior in Fig. 2 depends on the specific experiment due to the covariance and the anisotropy of $W$, the qualitative result applies generally.

Nonetheless, adding multipoles beyond the monopole helps to break degeneracies between parameters. The correlation matrices for each redshift bin after marginalizing over nuisance parameters are shown in Fig. 3. We compare the results with and without the hexadecapole. This figure shows a reduction of the correlation between the parameters when the hexadecapole is included with respect to the case in which it is not. Note that this reduction is smaller for larger redshifts (where the global S/N, and especially that of the hexadecapole, is smaller).

Lower absolute correlations have a positive impact on the final marginalized constraints. Table 1 reports the forecasted 68% confidence-level marginalized precision of the measurements of the BAO rescaling parameters and the parameter combination $\langle T\rangle f\sigma_{8}$, for each of the redshift bins of our general experiment. The fiducial values of $\alpha_{\perp}$ and $\alpha_{\parallel}$ are 1 for all redshifts, the corresponding fiducial values of $\langle T\rangle f\sigma_{8}$ are $\left\{0.78,\,0.67,\,0.56,\,0.43,\,0.29\right\}\,\mu$K for increasing $z$ (assuming the $\langle T\rangle$ values reported in Section IV.1).

The worsening of the forecasted constraints with redshift is expected, given the decreasing S/N (as shown in Fig. 2). Interestingly, the marginalized constraints on the BAO rescaling parameters are between 11% and 62% stronger when including the hexadecapole. The improvement obtained by including the hexadecapole decreases with redshift, as expected from the reduction of the difference between correlations with and without the hexadexapole shown in Fig. 3. Although measuring even higher multipoles might yield further gain, the improvement would be marginal, and the level of observational systematics would be too high to consider these measurements reliable. This is why we limit our study to $\ell\leq 4$.

V.2 Redshift binning and experimental optimization

When aiming to measure BAO and RSD, it is preferable to bin in redshift as much as possible, so that the evolution of the expansion of the Universe and growth of structure is better constrained. However, if the target is a parameter that does not change with redshift, wider redshift bins are more optimal, since that reduces the covariance of the LIM power spectrum (see Section III.2). Therefore, the analysis should always be adapted to the target parameter. We illustrate this fact by forecasting constraints on the sum $\sum m_{\nu}$ of the neutrino masses (when allowed to be larger than $0.06\,{\rm eV}$), and on deviations from primordial Gaussian perturbations.

Any imprint caused by primordial deviations from Gaussian initial conditions would have been preserved in the ultra-large scales of the matter power spectrum, which have remained outside the horizon since inflation. In the local limit, primordial non-Gaussianity can be parametrized with $f_{\rm NL}$ as the amplitude of the local quadratic contribution of a single Gaussian random field $\phi$ to the Bardeen potential $\Phi$ 999Here we assume the convention of large scale structure rather than the one for CMB ($f_{\rm NL}^{\rm LSS}\approx 1.3f_{\rm NL}^{\rm CMB}$ (Xia_PNG, )).:

[TABLE]

The skewness introduced in the density probability distribution by local primordial non-Gaussianity increases the number of massive objects, hence introducing a scale dependence on the halo bias (Matarrese_png00, ; Dalal_png07, ; Matarrese_png08, ; Desjacques_png, ). Denoting the Gaussian halo bias with $b^{\rm G}_{h}$, the total halo bias appearing in Eq. (3) is given by $b_{h}(k,z)=b^{\rm G}_{h}(z)+\Delta b_{h}(k,z)$, where:

[TABLE]

Here $\Omega_{m}$ is the matter density parameter at $z=0$, $\delta_{\rm ec}=1.68$ is the critical value of the matter overdensity for ellipsoidal collapse, $D(z)$ is the linear growth factor (normalized to 1 at $z=0$) and $T_{m}(k)$ is the matter transfer function (which is approximately $\sim 1$ at large scales).

The extensions we consider to $\Lambda$CDM therefore have $\sum m_{\nu}$ and $f_{\rm NL}$ as extra parameters, respectively, which we include in $\vec{\varsigma}$ in our Eq. (9), since they modify the power spectrum template. Naturally, varying either $\sum m_{\nu}$ or $f_{\rm NL}$ would also slightly modify the halo mass function and therefore alter the LIM power spectrum via modifications of $\langle T\rangle$, $P_{\rm shot}$ and $b$. We do not model this dependence, but our findings should not affected by this: the halo mass function (minimally) affects the amplitudes of the $P_{\rm clust}$ and $P_{\rm shot}$ terms, which we marginalize over.

Moreover, we should note that more massive neutrinos induce a slight scale dependence on the halo bias that can be efficiently removed if the power spectrum of baryons and cold dark matter (instead of the power spectrum of all matter, including neutrinos) is used (see e.g., Raccanelli_neutrino_worry ; Vagnozzi_neutrino_uncorrected ; Valcin_beHappy ). Since we do not aim here for detailed constraints, but for an illustration of the survey optimization, we ignore these effects and leave their study to future work.

We compare the performance using the redshift binning proposed in Section IV.1 with those of a single redshift bin centered at $z=3.80$ and covering the whole frequency band of the experiment. In order to maximize the coverage at large scales, we use logarithmic binning in $k$ in both cases. All measured quantities present in Eq. (9), except for $\vec{\varsigma}$ in some cases (e.g., $\sum m_{\nu}$ or $f_{\rm NL}$), are different in each redshift bin. Therefore, the final marginalized constraints on $\vec{\varsigma}$ obtained using several redshift bins is given by the result of the combination of each individual marginalized constraint obtained from each redshift bin:

[TABLE]

where $F_{z}$ is the Fisher matrix of a given redshift bin.

Fig. 4 shows a comparison of the forecasted marginalized 68% confidence-level constraints on $f_{\rm NL}$ and $\sum m_{\nu}$ from our generalized experiment using one and five redshift bins. It also shows the dependence of the constraints on $\Omega_{\rm field}$ and $t_{\rm obs}N_{\rm feeds}/T_{\rm sys}^{2}$. While the former affects the volume probed (which determines $N_{\rm modes}$ and $W_{\rm vol}$), the latter only affects the amplitude of $P_{n}$. We restrict our investigation of the optimization of the survey to the variation of these parameters (for a comprehensive study of survey optimization to constrain $f_{\rm NL}$, see Ref. Moradinezhad_Dizga_COfNL2 ).

Having five redshift bins reduces $L_{\parallel}$ in each of them, which suppresses the power spectrum up to higher $k_{\parallel}^{\rm min}$ values, through $W_{\rm vol}$. This is the main reason that our forecasted errors on $f_{\rm NL}$ and $\sum m_{\nu}$ using a single redshift bin are approximately a factor of two better. For high values of $\Omega_{\rm field}$, sample variance is smaller than the instrumental noise, which is why the constraints improve so much when $t_{\rm obs}N_{\rm feeds}/T_{\rm sys}^{2}$ increases. The improvement is smaller for low $\Omega_{\rm field}$ since the sample variance contribution to the error dominates and reducing the instrumental noise does not significantly reduce the total covariance. Since the imprints of $f_{\rm NL}$ lie on large scales, this effect is more critical than for $\sum m_{\nu}$, which affects a wider range of scales. Also, this effect is not very strong with five redshift bins, since the $\tilde{P}$ suppression limits the exploitation the very large line-of-sight scales.

V.3 Further improvements

Besides the strategies explored above, the extraction of cosmological information from the LIM power spectrum can be further improved. One possibility would be to assign an arbitrary normalized weight, $w$, to each of the observed voxels in configuration space. Then, $w(\vec{r})$ can be chosen to maximize the S/N of the measured IM power spectrum. With the introduction of these weights, a minimum-variance estimator can be derived, as proposed in Ref. (FKPweights, ) for galaxy surveys. An adaptation to LIM can be found in Ref. Blake_Pkmod .

Different redshift binning strategies or overlapping bins may be more optimal than the prescription used above. This will depend on the main objectives of the survey, as argued in Section V.2. However, if overlapping bins are to be used, they would not be independent anymore, and the corresponding covariance will have to be taken into account. Moreover, the actual optimal frequency band of the experiment depends on the targeted parameters.

Finally, in this work we have enforced the redshift bins to be narrow enough such that the redshift evolution of the large-scale brightness temperature fluctuations does not change significantly (except in Section V.2). However, this can be incorporated into the analysis using redshift-weighting techniques (Zhu_redshiftweight, ): redshift-dependent weights are chosen in order to minimize the projected error on the target cosmological or astrophysical parameter using a Fisher matrix forecast, so that the constraining power is maximized. These techniques can become especially useful when wide redshift bins are needed, e.g. for primordial non-Gaussianity, see Ref. Mueller_fNL .

VI Conclusions

LIM techniques have attracted substantial attention due to their untapped potential to constrain astrophysics and cosmology by probing huge volumes of the observable Universe which are beyond the reach of other methods. Nonetheless, the versatility of LIM can also be considered a nuisance, since cosmological information is always intrinsically degenerate with astrophysics. These degeneracies, as well as the effect of subtle contributions to the observed LIM power spectrum and its covariance, as well as various assumptions entailed in the analyses, are too commonly ignored or only partially treated.

The aim of this work has been to provide a comprehensive and general framework to optimally exploit the cosmological information encoded in the LIM power spectrum. We have presented a reparameterization of it to identify and isolate the degeneracies between cosmology and astrophysics, including redshift-space distortions and the Alcock-Paczynski effect. Furthermore, we have introduced and advocated for the use of the multipole expansion of the LIM power spectrum. We also derived an accurate analytic covariance for the multipoles, and discussed common errors in previous analyses.

Using a generalized experiment, and focusing on baryon acoustic oscillations and redshift-space distortions measurements, we showed that adding the hexadecapole to the monopole and quadrupole of the LIM power spectrum in the analysis returns a $10\%-60\%$ improvement in the forecasted constraints. This is mainly due to the reduction of the correlation between the parameters, as the S/N of hexadecapole measurements is usually low. However, redshift-space distortions as inferred from the LIM power spectrum are completely degenerate with $\langle T\rangle$. This degeneracy may be partially broken using mildly non-linear scales and perturbation theory, as proposed by Ref. Castorina_21cmgrowth . Other possible strategies include the combination with other cosmological probes Obuljen_21cmcrossothers or joint analyses of two-point and higher order statistics measurements GilMarin_bispectrunNbody ; Sarkar_21cmBispec .

We also investigated different survey strategies depending on the target of the experiment. As an illustration, we compared forecasts for constraints on $\sum m_{\nu}$ and local primordial non-Gaussianity using one or multiple redshift bins over the same total volume, finding that in both cases using one redshift bin returns stronger constraints. We also explored the improvement of the constraints by decreasing the instrumental noise and increasing the size of the patch of the sky observed. While the constraints on $\sum m_{\nu}$ do not significantly improve when varying $\Omega_{\rm field}$, the size of the patch of the sky observed is critical to constrain primordial non-Gaussianity, since its signatures are dominant at large scales (however, increasing $\Omega_{\rm field}$ beyond $\sim 4000-6000$ deg2 does not improve significantly the constraints, since in this case the instrumental error dominates the error budget at large scales).

Finally, we briefly discussed further improvements to the LIM power spectrum analyses to be studied in future works, especially to be applied to simulated and real observations, rather than theoretical work. In the Appendices, we provide an adaptation of our proposed framework to cross-power spectra between spectral lines or between LIM and galaxy surveys, and derive a close-to-optimal estimator for the true LIM power spectrum (i.e., without smoothing), and compare different approaches.

Though we have focused here on power spectrum analyses, line-intensity maps contain a significant amount of information beyond their power spectra. Line-intensity fluctuations are generated by highly complex gas dynamics on sub-galactic scales, which gives rise to a significantly non-Gaussian intensity field. Alternative summary statistics, such as higher order correlations or the voxel intensity distribution Breysse_VID , are then needed in order to exhaust the information encoded in LIM observations. Ideally, these summary statistics will be measured and analyzed with the power spectrum in tandem, accounting for their correlation, as first explored in Ref. (Ihle_VID-PS, ).

The methodology developed in this work should find ample opportunities for implementation. As a compelling example, we demonstrate in a companion paper, Ref. (Bernal_IM_letter, ), that LIM can be used to efficiently probe the expansion history of the universe up to extremely high redshifts ($z\lesssim 9$), possibly weighing in on the growing Hubble tension Freedman:2017yms ; RiessH0_19 ; Wong:2019kwg and suggested models to explain it (including models of evolving dark energy, exotic models of dark-matter, dark-matter–dark-energy interaction, modified gravity, etc., see e.g. Refs. Karwal:2016vyq ; Raveri_DEMG ; Kyriakos_dcdm ; Poulin_EDE ; Lin:2019qug ; Kreisch_selfnu ; Yang:2019nhz ; Archidiacono:2019wdp ; Hooper:2019gtx ; Vattis:2019efj ; Pandey:2019plg ; Pan:2019gop ; Desmond:2019ygn ; Agrawal:2019lmo ; Vagnozzi:2019ezj ). It can also be used to constrain model-independent expansion histories of the Universe BernalH0 ; Poulin_H0 ; Joudaki:2017zhq at the few-percent level. Our formalism could also be adapted to measure the velocity-induced acoustic oscillations (Munoz_vao_meth, ), recently proposed in Ref. Munoz_vao as a standard ruler at cosmic dawn.

We are eager for the next generation of LIM observations to be available for cosmological analyses, and hope that the general framework reported in this manuscript can guide its precise and robust exploitation.

Acknowledgements.

Funding for this work was partially provided by the Spanish MINECO under projects AYA2014-58747-P AEI/FEDER UE and MDM-2014-0369 of ICCUB (Unidad de Excelencia Maria de Maeztu). JLB is supported by the Spanish MINECO under grant BES-2015-071307, co-funded by the ESF, and thanks Johns Hopkins University for hospitality during early stages of this work. HGM acknowledges the support from ”la Caixa” Foundation (ID 100010434) which code LCF/BQ/PI18/11630024”.

Appendix A Formalism for cross-correlations

A.1 The IM cross-power spectrum

Let us consider two generic brightness-temperature maps of two different spectral lines, denoted by $X$, $Y$, respectively. In this case, the Kaiser effect present in the RSD factor is different for each tracer. For instance, in $F^{X}_{\rm RSD}\propto\left(1+f\mu^{2}/b_{X}\right)$, with $b_{X}$ being the luminosity-averaged bias for the line X. Then, the LIM cross-power spectra of the $X$ and $Y$ lines is given by:

[TABLE]

where we have dropped the explicit notation regarding dependence on $k$, $\mu$ and $z$ (as will be done hereafter) for the sake of readability (See Refs. Wolz2017 ; Breysse2019 for a derivation). Following the same arguments as in Section II, we express the cross-power spectrum between two different lines as:

[TABLE]

where $\tilde{T}_{XY}=\langle T_{X}\rangle\langle T_{Y}\rangle$. Since we have two different biases and two different brightness temperatures (one per line), the measurable combinations of parameters are:

[TABLE]

The covariance of the LIM cross-power spectrum is also slightly different. Considering, without loss of generality, the case in which each line is observed by a different experiment, the covariance per $\mu$ and $k$ bin of the cross-power spectra of two different lines is (neglecting mode coupling) given by:

[TABLE]

where $\tilde{\sigma}_{X}$ and $\tilde{\sigma}_{Y}$ are the square root of the corresponding covariances in Eq. (19).

Often, if the cross-power spectrum of two tracers can be measured, the corresponding auto-power spectra can be as well. In this case, when comparing the model to observations, all this information needs to be taken into account. In this case, Eqs. (22), (23), and (24) are still correct, but both the data vector and the covariance matrix need to be changed. The former will include both auto- and cross-power spectra, hence $\vec{\Theta}$ would be the concatenation of $\vec{\Theta}_{XX}$, $\vec{\Theta}_{XY}$ and $\vec{\Theta}_{YY}$. In turn, the covariance matrix will be formed by four square blocks of the same size, where the diagonal blocks would be $\tilde{\mathcal{C}}_{XX}$ and $\tilde{\mathcal{C}}_{YY}$, and the off-diagonal blocks, $\tilde{\mathcal{C}}_{XY}$.

A.2 Cross-correlation with galaxy surveys

LIM observations can also be cross-correlated with galaxy number counts. Denoting the galaxy catalog and related quantities with subscript/superscript $g$, we have an RSD factor $F^{g}_{\rm RSD}\propto\left(1+f\mu^{2}/b_{g}\right)$. In this case, the cross-power spectrum of galaxy number counts and LIM is:

[TABLE]

where $\langle\rho_{L}^{X}\rangle_{g}$ is the expected luminosity density sourced only from the galaxies belonging to the galaxy catalog used, and $n_{g}$ is the number density of such galaxies Wolz2017 ; Breysse2019 . Note the difference in the shot noise terms in Eq. (32) with respect to Eqs. (5) and (28). This is because contributions to $P_{\rm shot}^{Xg}$ come only from the locations occupied by the galaxies targeted by the galaxy survey, and the shot noise between the galaxy distribution and the luminosity sourced elsewhere vanishes. We also assume that the shot noise of the galaxy power spectrum is Poissonian (i.e., $P_{\rm shot}^{gg}=1/n_{g}$). Nonetheless, clustering and halo exclusion may introduce deviations from a Poissonian shot noise. This non-Poissonian contribution can change the amplitude of the shot noise and even induces a small scale dependence (Ginzburg_halomodel, ; Schmittfull_tracers, ).

Similarly to Eq. (29), we prefer to express the cross-power spectrum of LIM and the galaxy number counts as:

[TABLE]

In this case, the parameter combinations measured would be:

[TABLE]

This parameterization already accounts for the possible variation of the amplitude of $P_{\rm shot}^{Xg}$ due to the non-Poissonian contributions mentioned above, since we have marginalized over $P^{Xg}_{\rm shot}$. Note that if the goal is to measure $\langle\rho_{L}^{X}\rangle_{g}$, the non-Poissonian contribution to the shot noise needs to be explicitly modelled. We neglect the potential scale dependence introduced, which is a good approximation at this stage. As in the auto-spectrum case, one may be able to access some of the degenerate astrophysical information in this cross-spectrum through a map’s one-point statistics, in this case by conditioning the LIM statistics on those of the galaxy catalog Breysse2019b .

The covariance per $\mu$ and $k$ bin of the cross-power spectra of one line and galaxy number counts can be computed (neglecting mode coupling) as:

[TABLE]

where $\left(P_{gg}+1/n_{g}\right)/\sqrt{N_{\rm modes}}$ is the square root of the covariance of the galaxy power spectrum (when Poissonian shot noise is assumed).

Appendix B True vs. observed power spectra

In most of the literature, the effect of the instrument response, finite-volume surveyed and remaining residuals after foreground and interlopers removal is included in the error budget of the LIM power spectrum measurements, rather than incorporate it in the signal, as we do here. This implies that the smearing of $\delta T$ must be removed from the data (i.e., a deconvolution of the window function is required), and therefore the true LIM power spectrum would be estimated in the analysis.

The derivation of the observables, degeneracies and covariance is equivalent to the one presented in Sections II and III.2, with two exceptions. On one hand, the summary statistic would be $P$ instead of $\tilde{P}$, which affects the power spectrum multipoles (Eq. (16)), as shown in Fig. 1. On the other hand, since the noise power spectrum $P_{n}$, depicted in Eq. (17), corresponds to the noise of the observed temperature fluctuations, the inverse of the window needs to be applied to $P_{n}$ (i.e., the same operation to obtain the true power spectrum from the observations). In this case, the variance of $P$ per $k$ and $\mu$ bin (neglecting mode coupling) is:

[TABLE]

Eq. (36) has an immediate consequence on the S/N of the power spectrum multipoles. When computing $C_{\ell\ell^{\prime}}$ (following Eq. (21), but using $\sigma^{2}$ instead of $\tilde{\sigma}^{2}$), the integral over $\mu$ tends to asymptote to infinity beyond the first resolution limit, i.e., $k\gtrsim\min\left(k^{\rm max}_{\perp},k^{\rm max}_{\parallel}\right)$, due to the presence of $W$ in the denominator. The effect on $\mathcal{C}_{00}$ is shown in Fig. 5, considering only $W_{\rm res}$. This behavior was unnoticed in previous work because either the anistropic power spectrum was directly used (rather than the multipoles), or due to an incorrect computation of the covariance. In many previous works, the thermal-noise contribution to the monopole error is given as proportional to $1/\int d\mu W(k,\mu)$, when in fact it should be proportional to $\int d\mu 1/W(k,\mu)$. Given the exponential behavior of $W$, $W_{0}^{-1}$ differs substantially from $(W^{-1})_{0}$, where the subindex 0 denotes the monopole.

However, an optimal estimator of the multipoles of $P$ should contain the same information that is contained in $\tilde{P}$, if $W$ is accurately modelled. In the following we discuss a next-to-optimal estimator of the multipoles of the true power spectrum, $\hat{P}$. Let us consider that we could assign a weight, $w(\vec{k})$, to each mode of the brightness temperature perturbations in Fourier space, $\delta_{k}T$, such that the temperature fluctuations, $\mathcal{F}$, where $P=\langle\lvert\mathcal{F}\lvert^{2}\rangle$, becomes:

[TABLE]

where $A=V_{k}^{-1}\int{\rm d}^{3}\vec{k}w^{2}$, $V_{k}$ is the volume in Fourier space, and all quantities depend on $k$ and $\mu$. Hereinafter we will not show the explicit notation for the dependence for the sake of simplicity and readability.

Applying the weights also to the noise power spectrum, we obtain $\hat{P}=w^{2}P/A$ and $\hat{P}_{n}=w^{2}P_{n}/A$. If we use these estimators for Eq. (36), the covariance of the multipoles of $\hat{P}$ is given by:

[TABLE]

where $\mathcal{P}=(P+P_{n}/W)^{2}$. In order to obtain the minimum variance using this framework, we impose the stability of the covariance under small changes in $w$. Since most of the S/N comes from the monopole (see Figure 2), let us consider the minimization of the average of the covariance of the monopole in a shell in Fourier space. Imposing $w=w_{0}+\delta w$:

[TABLE]

where we have expanded over $\delta w/w_{0}$ up to second order. The stability condition of $\hat{\mathcal{C}}_{00}$ is fulfilled if:

[TABLE]

We compare the three modellings of the LIM power spectrum and the covariances discussed in this appendix in Fig. 5, only accounting for $W_{\rm res}$ and neglecting the effects of $W_{\rm vol}$ for the sake of clarity in this illustration. We show the monopole of the measured LIM power spectrum, $\tilde{P}$, of the true LIM power spectrum, $P$, and applying the weights of Eq. (40) to the true LIM power spectrum, $\hat{P}$, as well as the square root of the corresponding covariances. As can be seen, using $P$ is very limited, since the S/N tends to zero at the first resolution limit. On the contrary, $\hat{P}$ shows a slightly more optimal behaviour than $\tilde{P}$. However, as discussed in Section III, measuring $\delta_{k}T$ will be extremely difficult when the plane-parallel approximation breaks down. Therefore, we advocate for the use of $\tilde{P}$ and $\tilde{\mathcal{C}}_{\rm\ell\ell^{\prime}}$, since their measurement is more straightforward.