ULTRA-LARGE-SCALE COSMOLOGY IN NEXT-GENERATION EXPERIMENTS WITH SINGLE TRACERS

Consider a set of sources with comoving number density ${n}_{s}(\eta ,{\boldsymbol{x}})$ (as measured in their own rest frame) and 4-velocity ${u}_{s}^{\mu }.$ These sources emit photons with a wave vector ${k}^{\mu }\equiv {{dx}}^{\mu }/d\lambda$ (λ is an affine parameter of the photon geodesic) and rest-frame energy ${k}_{\mu }{u}_{s}^{\mu }.$ During an interval $d\lambda$ of the affine parameter, the photons cover a volume ${{dA}}_{e}({k}_{\mu }\;{u}_{s}^{\mu })d\lambda ,$ where dA_e is the invariant area of the wavefront corresponding to the observed solid angle $d{{\rm{\Omega }}}_{o}.$ Throughout, we have labelled quantities measured in the emitter's and observer's frames with subscripts e and o respectively.

The total number count in a redshift interval dz corresponding to $d\lambda$ is therefore

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{{dz}\;d{{\rm{\Omega }}}_{o}}={n}_{s}\;\displaystyle \frac{{{dA}}_{e}}{d{{\rm{\Omega }}}_{o}}({k}_{\mu }\;{u}_{s}^{\mu })\displaystyle \frac{d\lambda }{{dz}}.\end{eqnarray} \tag{ 1 }$$

Each of these terms is straightforward to compute in the absence of perturbations:

$$\begin{eqnarray*}{n}_{s}(\eta ,{\boldsymbol{x}}) & = & {\bar{n}}_{s}(\eta (z)),\\ \displaystyle \frac{{{dA}}_{e}}{d{{\rm{\Omega }}}_{o}} & = & {a}^{2}(\eta (z)){r}^{2}(z),\\ ({k}_{\mu }\;{u}_{s}^{\mu })\displaystyle \frac{d\lambda }{{dz}} & = & \displaystyle \frac{a(\eta (z))}{H(\eta (z))}.\end{eqnarray*}$$

Here $\eta (z)$ is the background conformal time at redshift z, r(z) is the background comoving angular diameter distance, $a(\eta (z))=1/(1+z)$ is the scale factor, and $H\equiv \dot{a}/a$ is the expansion rate. For the rest of this work we will assume a flat background cosmology, so that radial (χ) and angular distances are the same.

In the presence of inhomogeneities, all of these quantities are perturbed with respect to their background values at redshift z, and in general we can write

$$\begin{eqnarray}&&\eta (z,\hat{{\boldsymbol{n}}})\equiv \eta (z)+\delta \eta ,\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{n}_{s}(z,\hat{{\boldsymbol{n}}})={\bar{n}}_{s}(\eta (z))\left[1+{\delta }_{n}+\displaystyle \frac{\partial \mathrm{ln}{\bar{n}}_{s}}{\partial \eta }\delta \eta \right],\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{{dA}}_{e}}{d{{\rm{\Omega }}}_{o}}\equiv {a}^{2}(\eta (z)){\chi }^{2}(z)[1+2{\delta }_{\perp }],\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&({k}_{\mu }\;{u}_{s}^{\mu })\displaystyle \frac{d\lambda }{{dz}}\equiv \displaystyle \frac{a(\eta (z))}{H(\eta (z))}[1+{\delta }_{\parallel }],\end{eqnarray} \tag{ 5 }$$

where ${\delta }_{n}$ is the perturbation to n_s, and we have defined the perturbations to the conformal time, $\delta \eta ,$ transverse distance, ${\delta }_{\perp },$ and radial distance, ${\delta }_{\parallel }.$

One extra detail must be taken into account: not all sources are equally bright, and sources will generally be distributed according to a particular luminosity function, ${n}_{s}(\eta ,\mathrm{ln}L,{\boldsymbol{x}}),$ which we define as the density of sources in a logarithmic interval of luminosity:

$$\begin{eqnarray}&&{n}_{s}\equiv \displaystyle \frac{d(\#\;\mathrm{sources})}{{dVd}\mathrm{ln}L}.\end{eqnarray} \tag{ 6 }$$

Only sources with a flux (observed power per unit detector area) above a given detection cut, ${F}_{{\rm{cut}}},$ will be detected. Flux and luminosity are related by an inverse-square law in angular distance, so perturbations to the angular diameter distance will affect the measured flux. Linearizing with respect to these perturbations gives

$$\begin{eqnarray}&&F(z,\hat{{\boldsymbol{n}}})=\displaystyle \frac{L}{4\pi {(1+z)}^{4}\;{a}^{2}(\eta (z)){\chi }^{2}(z)}[1-2{\delta }_{\perp }].\end{eqnarray} \tag{ 7 }$$

At a given redshift and flux cut, we will only observe sources with luminosities above a threshold ${L}_{{\rm{cut}}},$ related to ${F}_{{\rm{cut}}}$ by the previous equation. In order to take this into account, we must replace n_s by the cumulative luminosity function,

$$\begin{eqnarray}&&{\mathcal{N}}(\eta ,{\boldsymbol{x}},\gt \mathrm{ln}L)\equiv {\displaystyle \int }_{\mathrm{ln}L}^{\infty }d\mathrm{ln}L^{\prime} \;{n}_{s}(\eta ,{\boldsymbol{x}},\mathrm{ln}L^{\prime} ),\end{eqnarray} \tag{ 8 }$$

so that Equation (3) becomes

$$\begin{eqnarray*}&&{\mathcal{N}}(z,{\boldsymbol{n}},{F}_{{\rm{cut}}})=\bar{{\mathcal{N}}}\left[1+{\delta }_{{\mathcal{N}}}+\displaystyle \frac{\partial \mathrm{ln}\bar{{\mathcal{N}}}}{\partial \eta }\;\delta \eta -2\displaystyle \frac{{\bar{n}}_{s}}{\bar{{\mathcal{N}}}}\;{\delta }_{\perp }\right].\end{eqnarray*}$$

We have shortened our notation such that

$$\begin{eqnarray}&&\bar{{\mathcal{N}}}\equiv \bar{{\mathcal{N}}}(\eta (z),\gt \mathrm{ln}{\bar{L}}_{{\rm{cut}}}(z,{F}_{{\rm{cut}}}))\end{eqnarray} \tag{ 9 }$$

(and likewise for ${\bar{n}}_{s}$), and have overlined all quantities evaluated in the background.

The full linear expression for the source number counts can finally be written as

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{{dzd}{{\rm{\Omega }}}_{o}}=\bar{{\mathcal{N}}}\;\displaystyle \frac{{a}^{3}(z)}{H(z)}\;{\chi }^{2}(z)[1+{{\rm{\Delta }}}_{N}(z,\hat{{\boldsymbol{n}}})],\end{eqnarray} \tag{ 10 }$$

where the perturbation is given by

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{N}(z,\hat{{\boldsymbol{n}}})={\delta }_{{\mathcal{N}}}+\displaystyle \frac{\partial \mathrm{ln}\bar{{\mathcal{N}}}}{\partial \eta }\;\delta \eta +{\delta }_{\parallel }+2{\delta }_{\perp }\left[1-\displaystyle \frac{{\bar{n}}_{s}}{\bar{{\mathcal{N}}}}\right].\end{eqnarray} \tag{ 11 }$$

In order to simplify the notation, from now on we will refer to the observed background number of sources found per unit redshift and solid angle simply as $\bar{N}(z),$ i.e.,

$$\begin{eqnarray}&&\bar{N}(z)\equiv \bar{{\mathcal{N}}}\;\displaystyle \frac{{a}^{3}(z)}{H(z)}\;{\chi }^{2}(z).\end{eqnarray} \tag{ 12 }$$

The terms $\delta \eta ,\;$ ${\delta }_{\parallel },$ and ${\delta }_{\perp }$ can be related to the metric, density, and velocity perturbations by solving the geodesic equation for photons in any gauge. In the conformal Newtonian gauge, defined by the line element

$$\begin{eqnarray}&&{{ds}}^{2}={-a}^{2}(\eta )\left[(1+2\psi )d{\eta }^{2}-(1-2\phi ){\delta }_{{ij}}{{dx}}^{i}{{dx}}^{j}\right],\end{eqnarray} \tag{ 13 }$$

these perturbations read

$$\begin{eqnarray}&&{\mathcal{H}}\delta \eta =-\psi +\int (\phi ^{\prime} +\psi ^{\prime} )d\eta +{v}_{r}\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{\delta }_{\parallel }=\left[1-\displaystyle \frac{{\mathcal{H}}^{\prime} }{{{\mathcal{H}}}^{2}}\right]{\mathcal{H}}\delta \eta +\psi +{v}_{r}\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&\quad +\displaystyle \frac{1}{{\mathcal{H}}}\left[-\displaystyle \frac{d\psi }{d\eta }+\phi ^{\prime} +\psi ^{\prime} +\displaystyle \frac{{{dv}}_{r}}{d\eta }\right]\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{\delta }_{\perp }={\mathcal{H}}\delta \eta -\displaystyle \frac{1}{\chi }\left[\delta \eta +\int (\phi +\psi )d\eta \right]-\phi -\kappa \end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&\kappa \equiv \displaystyle \frac{1}{2}{\displaystyle \int }_{0}^{\chi }\displaystyle \frac{\chi -\chi ^{\prime} }{\chi \chi ^{\prime} }{{\rm{\nabla }}}_{{\rm{\Omega }}}^{2}(\phi +\psi )d\chi ^{\prime} ,\end{eqnarray} \tag{ 18 }$$

where ${v}_{r}\equiv \hat{{\boldsymbol{n}}}\cdot {{\boldsymbol{v}}}_{s}$ is the radial peculiar velocity of the sources, ${{\rm{\nabla }}}_{{\rm{\Omega }}}^{2}$ is the Laplacian on the unit sphere, and κ is the lensing convergence. Note that we have denoted all partial derivatives with respect to conformal time as ${\partial }_{\eta }b\equiv b^{\prime}$ (and ${\mathcal{H}}\equiv a^{\prime} /a$), and that the operator $d/d\eta$ denotes a total lightcone derivative along the unperturbed photon trajectory,

$$\begin{eqnarray}&&\displaystyle \frac{{db}}{d\eta }\equiv \displaystyle \frac{d}{d\eta }[b(\eta ,{\boldsymbol{x}}=({\eta }_{0}-\eta )\hat{{\boldsymbol{n}}})],\end{eqnarray} \tag{ 18 }$$

where ${\eta }_{0}$ is the age of the universe. Likewise, all integrals shown in the equations above must be understood as lightcone integrals along the same trajectory.

Besides source number counts, another promising observational tool for studying large-scale structure is a technique known as intensity mapping. The technical details of this method are discussed in Section 4.1, but we will describe the relevant relativistic effects here (see also Hall et al. 2013).

In intensity mapping, the observable used to trace the matter density is the intensity received from a line-emitting medium integrated over a patch of the sky (i.e., the total power measured in a frequency interval per unit detector area and observed solid angle). We assume that this line emission is caused by some well-defined transition line, and can therefore be used to recover the redshift of the source by comparing the observed frequency with the known rest-frame one. In the rest frame of a set of line-emitting sources, the emissivity is

$$\begin{eqnarray}&&\displaystyle \frac{{{dE}}_{e}}{{{dt}}_{e}d{\nu }_{e}d{{\rm{\Omega }}}_{e}{{dV}}_{e}}=\displaystyle \frac{{\hslash }}{2}{A}_{21}{\nu }_{e}\varphi ({\nu }_{e})\displaystyle \frac{{x}_{2}\;{\rho }_{a}}{{m}_{a}},\end{eqnarray} \tag{ 20 }$$

where A₂₁ is the Einstein coefficient for the transition, ${\rho }_{a}$ is the comoving density of the emitting gas, m_a is its atomic mass, x₂ is the (number) fraction of the gas in the excited state, and $\varphi (\nu )$ is the line profile (normalized to unity when integrated over all frequencies).

As shown in Section 2.1.1, the volume covered by the emitted photons in an affine parameter interval $d\lambda$ is

$$\begin{eqnarray}&&{{dV}}_{e}={{dA}}_{e}({k}_{\mu }\;{u}_{s}^{\mu })d\lambda .\end{eqnarray} \tag{ 21 }$$

Assuming that no absorption or scattering of the emitted photons occurs, the emissivity can be related to the measured intensity by

$$\begin{eqnarray*}I({\nu }_{o},\hat{{\boldsymbol{n}}}) & \equiv & \displaystyle \frac{{{dE}}_{o}}{{{dt}}_{o}{{dA}}_{o}d{{\rm{\Omega }}}_{o}d{\nu }_{o}}\\ & = & \displaystyle \frac{{\hslash }{\nu }_{o}{A}_{21}{x}_{2}{\rho }_{a}}{2{m}_{a}}\varphi ({\nu }_{e})({k}_{\mu }{u}_{s}^{\mu })d\lambda \displaystyle \frac{{{dA}}_{e}d{{\rm{\Omega }}}_{e}}{{{dA}}_{o}{{\rm{\Omega }}}_{o}}\displaystyle \frac{d{\nu }_{e}{{dt}}_{e}}{d{\nu }_{o}{{dt}}_{o}}.\end{eqnarray*}$$

The frequencies and time intervals in both frames are directly related through the redshift z, as are angles and invariant areas (from Etherington's reciprocity relation),

$$\begin{eqnarray}&&\displaystyle \frac{{{dt}}_{e}d{\nu }_{e}}{{{dt}}_{o}d{\nu }_{o}}=1,\quad \displaystyle \frac{{{dA}}_{e}d{{\rm{\Omega }}}_{e}}{{{dA}}_{o}d{{\rm{\Omega }}}_{o}}=\displaystyle \frac{1}{{(1+z)}^{2}}.\end{eqnarray} \tag{ 22 }$$

Using these relations, and assuming that observations will take place on frequency intervals ${\rm{\Delta }}{\nu }_{o}$ much larger than the line width, we finally obtain the relation

$$\begin{eqnarray}&&I({\nu }_{o},\hat{{\boldsymbol{n}}})=\displaystyle \frac{{\hslash }{A}_{21}{\nu }_{21}{x}_{2}}{2{m}_{a}{(1+z)}^{2}}\;{\rho }_{a}({k}_{\mu }\;{u}_{s}^{\mu })\displaystyle \frac{d\lambda }{{dz}},\end{eqnarray} \tag{ 23 }$$

where ${\nu }_{21}$ is the rest-frame line frequency. We can see that this is equivalent to Equation (1) for number counts, except for the factor of the angular diameter distance, ${{dA}}_{e}/d{{\rm{\Omega }}}_{o}.$ This is because the observable in intensity mapping is not the total number of objects in a given patch of sky, but the combined emitted light from the same patch. Since luminosities and angular distances are affected in the same way by lightcone effects, these cancel exactly for intensity mapping.

Expanding both ${\rho }_{a}$ and $({k}_{\mu }\;{u}_{s}^{\mu })\frac{d\lambda }{{dz}}$ to linear order, we can therefore compute all of the linear perturbations to the observed intensity mapping signal,

$$\begin{eqnarray}I(\nu ,\hat{{\boldsymbol{n}}}) & = & \bar{I}(\nu )[1+{{\rm{\Delta }}}_{I}]\\ & & \equiv \bar{I}(\nu )\left[1+{\delta }_{a}+\displaystyle \frac{\partial \mathrm{ln}{\bar{\rho }}_{a}}{d\eta }\delta \eta +{\delta }_{\parallel }\right],\end{eqnarray} \tag{ 24 }$$

where ${\delta }_{a}$ is the intrinsic perturbation of the emitting gas density. By comparing this with Equation (11), we can see that the linear perturbation for intensity mapping is equivalent to the perturbation of the number counts for a population of sources with a particular form of the luminosity function, such that the number of sources observed above a given flux is proportional to the luminosity associated with that flux,

$$\begin{eqnarray}&&{\bar{{\mathcal{N}}}}_{{\rm{IM}}}(\gt L)\propto L.\end{eqnarray} \tag{ 25 }$$

The somewhat unfortunate consequence of this result is that there are no linear perturbations to angular distances for intensity mapping,⁷ which could potentially reduce the amount of cosmological information that can be extracted from this probe.

The background term in Equation (24) is commonly expressed in terms of antenna temperature, defined through the Rayleigh–Jeans relation for a blackbody emitter $T(\nu )=I(\nu ){c}^{2}/(2{k}_{{\rm{B}}}{\nu }^{2}),$ where k_B is the Boltzmann constant. In terms of background quantities, the homogeneous intensity mapping signal is

$$\begin{eqnarray}&&\bar{T}(z)=\displaystyle \frac{3{\hslash }{A}_{21}{x}_{2}{c}^{2}}{32\pi {{Gk}}_{{\rm{B}}}{m}_{a}{\nu }_{21}^{2}}\displaystyle \frac{{H}_{0}^{2}\;{{\rm{\Omega }}}_{b,0}\;{x}_{a}(z){(1+z)}^{2}}{H(z)},\end{eqnarray} \tag{ 26 }$$

where x_a(z) is the fraction of baryons made up by the line-emitting species under study.

Until now, we have not related the intrinsic perturbation in the number density of sources to the perturbations of the energy–momentum tensor. Assuming that galaxies form in dark matter haloes, which themselves form preferentially in high-density regions, one would expect the halo (or galaxy) number overdensity to trace the fluctuations in the overall matter density on large (linear) scales with a simple linear bias factor, ${\delta }_{{\rm{halo}}}\simeq {b}_{{\rm{halo}}}\;{\delta }_{M}.$ This bias is a central piece of the halo model of structure formation (Mo & White 1996; Peacock & Smith 2000), the validity of which has been extensively tested against numerical simulations (Cole et al. 2008). Although the linear bias is expected to be scale-dependent on nonlinear and mildly nonlinear scales, where nonlinear and stochastic bias terms could also be important, we are mainly interested in large-scale observables in this paper, where the approximation of a linear, scale-independent bias should be valid (although see Section 2.2). This bias will depend on redshift and luminosity, however (e.g., more luminous, and therefore rarer, objects are expected to be more highly biased).

Since it is not possible to unambiguously define the matter overdensity ${\delta }_{M}$ in a gauge-invariant way in a general-relativistic context, a subtle point is the choice of overdensity field on which the bias relation is applied. In this work we take the point of view of Challinor & Lewis (2011), Baldauf et al. (2011), Jeong et al. (2012), and Bruni et al. (2012), and argue that, since the process of galaxy formation is due to local physics, and since we expect our sources to follow the same velocity field as the dark matter, the bias relation should be applied in the synchronous comoving gauge. Note that it is also the comoving gauge perturbation that appears in the Poisson equation. A more complete discussion of this argument, which can also be extended to the case of primordial non-Gaussianity, can be found in Baldauf et al. (2011).

It follows that the intrinsic perturbation to the number density of sources, ${\delta }_{N},$ in the Newtonian gauge—our choice for this work—is related to the matter overdensity in the synchronous comoving gauge, ${\delta }_{M,{\rm{syn}}},$ through

$$\begin{eqnarray}&&{\delta }_{{\rm{N}}}=b(L,z,k){\delta }_{M,{\rm{syn}}}+\displaystyle \frac{\partial \mathrm{ln}\bar{{\mathcal{N}}}}{d\eta }\displaystyle \frac{v}{k},\end{eqnarray} \tag{ 27 }$$

where v is the peculiar velocity in the Newtonian gauge, and we have allowed the bias, b, to be scale-dependent, in anticipation of the discussion in Section 2.2.

The amplitudes of the perturbations to the conformal time and transverse distances depend explicitly on the derivatives of the luminosity function of the source population with respect to luminosity and time (see Equation (11)). It has become common to express these derivatives in terms of the so-called magnification bias, $s(\eta ),$ and evolution bias, ${f}_{{\rm{evo}}}(\eta ),$ defined as

$$\begin{eqnarray}&&s(\eta )\equiv \displaystyle \frac{5}{2}\displaystyle \frac{{\bar{n}}_{s}(\eta ,\mathrm{ln}{\bar{L}}_{{\rm{cut}}})}{\bar{{\mathcal{N}}}(\eta ,\gt \mathrm{ln}{\bar{L}}_{{\rm{cut}}})},\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&{f}_{{\rm{evo}}}(\eta )\equiv \displaystyle \frac{\partial \mathrm{ln}[{a}^{3}\bar{{\mathcal{N}}}(\eta ,\gt \mathrm{ln}{\bar{L}}_{{\rm{cut}}})]}{\partial \mathrm{ln}a}.\end{eqnarray} \tag{ 29 }$$

As with the clustering bias, the values of s and ${f}_{{\rm{evo}}}$ depend on the source population under study, so must be modeled correctly in order to maximize the information that can be extracted from any clustering analysis. While b(z) must be determined directly from clustering statistics, it is possible to estimate s and ${f}_{{\rm{evo}}}$ directly from the overall number counts of sources as a function of redshift and magnitude. Let $\bar{N}(z,\lt {m}_{*})$ be the cumulative number of sources with magnitude m brighter than ${m}_{*},$ per unit solid angle and redshift interval. $\bar{N}$ is related to the luminosity function, ${\bar{n}}_{s},$ through

$$\begin{eqnarray}&&\bar{N}(z,\lt {m}_{*})=\displaystyle \frac{c\;{\chi }^{2}(z)}{{(1+z)}^{3}\;H(z)}{\displaystyle \int }_{\mathrm{ln}{L}_{*}}^{\infty }{\bar{n}}_{s}(\eta ,\mathrm{ln}L)d\mathrm{ln}L,\end{eqnarray} \tag{ 30 }$$

where the threshold luminosity L_* is ${L}_{*}=4\pi {(1+z)}^{2}{\chi }^{2}(z){F}_{*},$ and fluxes and magnitudes are related through

$$\begin{eqnarray}&&m=-\displaystyle \frac{5}{2}{\mathrm{log}}_{10}\left[\displaystyle \frac{F}{{F}_{0}}\right].\end{eqnarray} \tag{ 31 }$$

Note that we have neglected evolution and k-corrections.

Using the definitions of ${f}_{{\rm{evo}}}$ and s in Equations (28) and (29), these quantities can be related to the derivatives of $\bar{N}$ with respect to z and m_* by

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\mathrm{log}}_{10}\bar{N}}{\partial {m}_{*}}=s,\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\mathrm{log}}_{10}\bar{N}}{\partial {\mathrm{log}}_{10}(1+z)}=\displaystyle \frac{(2-5\;s)}{\chi \;{aH}}-5\;s+\displaystyle \frac{H^{\prime} }{{{aH}}^{2}}-{f}_{{\rm{evo}}}.\end{eqnarray} \tag{ 33 }$$

Note that, in order to use these relations to estimate s and ${f}_{{\rm{evo}}}$, it is necessary to have full redshift information about the source distribution. While this is available by default for spectroscopic surveys, determining the redshift distribution becomes more involved for photometric and radio continuum surveys. This is nevertheless a necessary task if these probes are to be usable for cosmological studies, where, e.g., the redshift and photometric redshift distributions must be correctly modeled. In any case, the uncertainties on s and ${f}_{{\rm{evo}}}$ will tend to grow toward large z, and must therefore be taken into account in any cosmological analysis.

As we described in Section 2.1.2, the case of intensity mapping is slightly different. In this case, perturbations to the angular distance cancel, which is equivalent to setting the magnification bias to the critical value ${s}_{{\rm{IM}}}=2/5$. ${f}_{{\rm{evo}}}$ can be determined directly from the redshift dependence of the background brightness temperature, $\bar{T}(z)$ (Equation (26)).

The most informative observable regarding the clustering of astrophysical sources is their two-point correlation, $\langle {{\rm{\Delta }}}_{{\mathcal{N}}}({z}_{1},{\hat{{\boldsymbol{n}}}}_{1}){{\rm{\Delta }}}_{{\mathcal{N}}}({z}_{2},{\hat{{\boldsymbol{n}}}}_{2})\rangle .$ The perturbation ${{\rm{\Delta }}}_{{\mathcal{N}}}$ can be expressed in terms of spherical harmonic coefficients,

$$\begin{eqnarray}&&{a}_{{\ell }m}(z)\equiv \int d\hat{{\boldsymbol{n}}}\;{{\rm{\Delta }}}_{{\mathcal{N}}}(z,\hat{{\boldsymbol{n}}}){Y}_{{\ell }m}(\hat{{\boldsymbol{n}}}),\end{eqnarray} \tag{ 34 }$$

where ${Y}_{{\ell }m}(\hat{{\boldsymbol{n}}})$ are the spherical harmonics. The clustering of number counts can then be studied through the angular power spectrum, defined by the correlation

$$\begin{eqnarray}&&\langle {a}_{{\ell }m}({z}_{1}){a}_{{\ell }^{\prime} m^{\prime} }^{*}({z}_{2})\rangle \equiv {\delta }_{{\ell }{\ell }^{\prime} }{\delta }_{{mm}^{\prime} }\;{C}_{{\ell }}({z}_{1},{z}_{2}),\end{eqnarray} \tag{ 35 }$$

where angle brackets denote an ensemble average.

In practice, ${\rm{\Delta }}(z,\hat{{\boldsymbol{n}}})$ is not measured in infinitesimal intervals of z, but by averaging over a set of finite radial bins, which we will label here by a Latin index, i. The observed anisotropy in bin i is

$$\begin{eqnarray}&&{a}_{{\ell }m}^{i}\equiv \int {dz}\;{W}_{i}(z){\rm{\Delta }}(z,\hat{{\boldsymbol{n}}}),\end{eqnarray} \tag{ 36 }$$

where the window function W_i is normalized to 1 when integrated over redshift. The shape of W_i is determined by both the background redshift distribution of observed sources, $\bar{N}(z),$ and the probability that a source at redshift z will be included in the ith bin ${p}_{i}(z),$ so that

$$\begin{eqnarray}&&{W}_{i}(z)\propto \bar{N}(z){p}_{i}(z).\end{eqnarray} \tag{ 37 }$$

Using this, one can show that the cross-spectrum between two bins can be written as (Di Dio et al. 2013)

$$\begin{eqnarray}&&{C}_{{\ell }}^{{ij}}=4\pi {\displaystyle \int }_{0}^{\infty }\displaystyle \frac{{dk}}{k}{\mathcal{P}}(k){{\rm{\Delta }}}_{{\ell }}^{i}(k){{\rm{\Delta }}}_{{\ell }}^{j}(k),\end{eqnarray} \tag{ 38 }$$

where ${\mathcal{P}}(k)$ is the dimensionless primordial power spectrum, which is assumed to take the form ${\mathcal{P}}(k)={A}_{s}{(k/{k}_{0})}^{{n}_{s}-1},$ and ${{\rm{\Delta }}}_{{\ell }}^{i}(k)$ contains the transfer functions of the terms that contribute to the anisotropy in bin i, in Fourier space and projected on the sky. Expanding the various contributions to Equation (11), the functions ${{\rm{\Delta }}}_{{\ell }}^{i}(k)$ can be written as a sum of 10 terms corresponding to different physical effects (Di Dio et al. 2013):

$$\begin{eqnarray}{{\rm{\Delta }}}_{{\ell }}^{{\rm{D}},i}(k) & \equiv & \displaystyle \int d\eta \;b\;{\tilde{W}}_{i}\;{\delta }_{M,{\rm{syn}}}(k,\eta ){j}_{{\ell }}(k\chi (\eta )),\\ {{\rm{\Delta }}}_{{\ell }}^{{\rm{RSD}},i}(k) & \equiv & \displaystyle \int d\eta {({aH})}^{-1}{\tilde{W}}_{i}(\eta )\theta (k,\eta ){j}_{{\ell }}^{\prime\prime }(k\chi (\eta )),\end{eqnarray} \tag{ 39 }$$

$$\begin{eqnarray}{{\rm{\Delta }}}_{{\ell }}^{{\rm{L}},i}(k) & \equiv & {\ell }({\ell }+1)\displaystyle \int d\eta \;{\tilde{W}}_{i}^{{\rm{L}}}(\eta )(\phi +\psi )(k,\eta ){j}_{{\ell }}(k\chi (\eta )),\\ {{\rm{\Delta }}}_{{\ell }}^{{\rm{V1}},i}(k) & \equiv & \displaystyle \int d\eta ({f}_{{\rm{evo}}}-3){aH}\;{\tilde{W}}_{i}(\eta )\displaystyle \frac{\theta (k,\eta )}{{k}^{2}}\;{j}_{{\ell }}(k\chi (\eta )),\end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray}{{\rm{\Delta }}}_{{\ell }}^{{\rm{V2}},i}(k) & \equiv & \displaystyle \int d\eta \left(1+\displaystyle \frac{H^{\prime} }{{{aH}}^{2}}+\displaystyle \frac{2-5\;s}{\chi \;{aH}}+5\;s-{f}_{{\rm{evo}}}\right)\\ & & \times {\tilde{W}}_{i}(\eta )\displaystyle \frac{\theta (k,\eta )}{k}\;{j}_{{\ell }}^{\prime }(k\chi (\eta )),\end{eqnarray} \tag{ 41 }$$

$$\begin{eqnarray}{{\rm{\Delta }}}_{{\ell }}^{{\rm{P1}},i}(k) & \equiv & \displaystyle \int d\eta \left(2+\displaystyle \frac{H^{\prime} }{{{aH}}^{2}}+\displaystyle \frac{2-5\;s}{\chi \;{aH}}+5\;s-{f}_{{\rm{evo}}}\right)\\ & & \times {\tilde{W}}_{i}(\eta )\psi (k,\eta ){j}_{{\ell }}(k\chi (\eta )),\end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray}{{\rm{\Delta }}}_{{\ell }}^{{\rm{P2}},i}(k) & \equiv & \displaystyle \int d\eta \;(5\;s-2){\tilde{W}}_{i}(\eta )\phi (k,\eta ){j}_{{\ell }}(k\chi (\eta )),\\ {{\rm{\Delta }}}_{{\ell }}^{{\rm{P3}},i}(k) & \equiv & \displaystyle \int d\eta \;{({aH})}^{-1}{\tilde{W}}_{i}(\eta )\phi ^{\prime} (k,\eta ){j}_{{\ell }}(k\chi (\eta )),\end{eqnarray} \tag{ 43 }$$

$$\begin{eqnarray}{{\rm{\Delta }}}_{{\ell }}^{{\rm{P4}},i}(k) & \equiv & \displaystyle \int d\eta \;{\tilde{W}}_{i}^{{\rm{P4}}}(\eta )(\phi +\psi )(k,\eta ){j}_{{\ell }}(k\chi (\eta )),\\ {{\rm{\Delta }}}_{{\ell }}^{{\rm{ISW}},i}(k) & \equiv & \displaystyle \int d\eta \;{\tilde{W}}_{i}^{{\rm{ISW}}}(\eta )(\phi +\psi )^{\prime} (k,\eta ){j}_{{\ell }}(k\chi (\eta )),\end{eqnarray} \tag{ 44 }$$

where we have defined the window functions

$$\begin{eqnarray}{\tilde{W}}_{i}(\eta (z)) & \equiv & {W}_{i}(z){\left(\displaystyle \frac{d\eta }{{dz}}\right)}^{-1},\\ {\tilde{W}}_{i}^{{\rm{L}}}(\eta ) & \equiv & {\displaystyle \int }_{0}^{\eta }d\eta ^{\prime} {\tilde{W}}_{i}(\eta ^{\prime} )\displaystyle \frac{2-5\;s(\eta ^{\prime} )}{2}\displaystyle \frac{\chi (\eta )-\chi (\eta ^{\prime} )}{\chi (\eta )\chi (\eta ^{\prime} )},\end{eqnarray} \tag{ 45 }$$

$$\begin{eqnarray}{\tilde{W}}_{i}^{{\rm{P4}}}(\eta ) & \equiv & {\displaystyle \int }_{0}^{\eta }d\eta ^{\prime} {\tilde{W}}_{i}(\eta ^{\prime} )\displaystyle \frac{2-5\;s}{\chi },\\ {\tilde{W}}_{i}^{{\rm{ISW}}}(\eta ) & \equiv & {\displaystyle \int }_{0}^{\eta }d\eta ^{\prime} {\tilde{W}}_{i}(\eta ^{\prime} )\\ & & \times {\left(1+\displaystyle \frac{H^{\prime} }{{{aH}}^{2}}+\displaystyle \frac{2-5\;s}{\chi \;{aH}}+5\;s-{f}_{{\rm{evo}}}\right)}_{\eta ^{\prime} }.\end{eqnarray} \tag{ 46 }$$

In these equations, the quantities ${\delta }_{M,{\rm{syn}}}(k,\eta ),\;$ $\theta (k,\eta ),\;$ $\psi (k,\eta ),$ and $\phi (k,\eta )$ are the transfer functions for the perturbation in the synchronous comoving gauge matter density, the divergence of the peculiar velocity, and the two metric potentials respectively.

Each term is sourced by a different physical effect. ${{\rm{\Delta }}}_{{\ell }}^{{\rm{D}}}$ corresponds to the intrinsic perturbation in the comoving number density of sources, which is the dominant contribution in most cases, and is the only term that has traditionally been taken into account when RSDs and lensing can be neglected. ${{\rm{\Delta }}}_{{\ell }}^{{\rm{RSD}}}$ is the usual RSD term corresponding to the Kaiser effect, due to the deformation of the Lagrangian volume in redshift space. ${{\rm{\Delta }}}_{{\ell }}^{{\rm{L}}}$ is the lensing convergence term, caused by the deformation of the Lagrangian volume in the transverse directions due to weak lensing. The terms ${{\rm{\Delta }}}_{{\ell }}^{{\rm{V1}}}$ and ${{\rm{\Delta }}}_{{\ell }}^{{\rm{V2}}}$ are extra RSD contributions that come from evaluating the background terms at a redshift perturbed by the Doppler effect. The remaining terms correspond to the same effect, but for redshift perturbations caused by gravitational redshifting instead of peculiar velocities. In particular, ${{\rm{\Delta }}}_{{\ell }}^{{\rm{ISW}}}$ is the analogue of the integrated Sachs–Wolfe (ISW) effect (Sachs & Wolfe 1967) for number counts.

Of these terms, the first three give the largest contribution to the total clustering anisotropy, so only these have traditionally been included in clustering analyses. The remaining terms are mainly relevant on super-horizon scales at the position of the sources, and even on those scales their amplitude is significantly smaller than the first three (see Figure 1). Nevertheless, these terms contain useful information that could potentially be used, for example, to constrain different theories of gravity (Lombriser et al. 2013; T. Baker & P. Bull 2015, in preparation). One of the aims of this paper is to forecast the detectability of these terms by future experiments. In order to do so, we have defined an effective parameter, ${\epsilon }_{{\rm{GR}}},$ which multiplies the terms ${{\rm{\Delta }}}_{{\ell }}^{{\rm{V1,2}}},{{\rm{\Delta }}}_{{\ell }}^{{\rm{P}}1-4},\;{\rm{and}}\;{{\rm{\Delta }}}_{{\ell }}^{{\rm{ISW}}}$ and has a fiducial value of 1. ${\epsilon }_{{\rm{GR}}}$ therefore parameterizes the amplitude of the relativistic corrections to the clustering of sources.

Zoom In Zoom Out Reset image size

Even though the origin of the lensing convergence term (${{\rm{\Delta }}}_{{\ell }}^{{\rm{L}}}$) is clearly general-relativistic, we have not included it under the umbrella of ${\epsilon }_{{\rm{GR}}}$ for two main reasons: first, we would like to focus on yet-undetected effects, and lensing magnification has already been detected by cross-correlating pairs of distant tracers (Scranton et al. 2005; Hildebrandt et al. 2009). Second, in this work we aim to identify possibly detectable observables on ultra-large scales, but the lensing term has a non-negligible effect on small angular scales. Nevertheless, for completeness we have also forecasted for the detectability of lensing magnification by defining an effective amplitude, ${\epsilon }_{{\rm{WL}}},$ multiplying ${{\rm{\Delta }}}_{{\ell }}^{{\rm{L}}}.$ In keeping with the main aim of this paper, however, we will only produce forecasts for this parameter based on its effects on the largest angular scales (lowest multipoles).

Thus, to clarify the terminology used here, we will refer to the terms ${{\rm{\Delta }}}^{{\rm{V1,2}}},$ ${{\rm{\Delta }}}^{{\rm{P}}1-4}$, and ${{\rm{\Delta }}}^{{\rm{ISW}}}$ as "GR effects" or "GR terms" and to ${{\rm{\Delta }}}^{{\rm{L}}}$ as the "lensing term," even if the nature of the latter is clearly relativistic.

As we have argued above, the two main sources of extra power on ultra-large scales are primordial non-Gaussianity and relativistic corrections. In order to study the detectability of the latter, we can therefore treat only ${f}_{{\rm{NL}}}$ and ${\epsilon }_{{\rm{GR}}}$ as free parameters, and fix the rest to their fiducial values. While this procedure would clearly yield an optimistic prediction of the actual constraint on ${\epsilon }_{{\rm{GR}}},$ it mimics what a survey attempting a first detection of any new effect would do: fix all non-degenerate parameters to their best-fit values, and fit for the amplitude of the terms related to the new effect. If, in doing this, the S/N on the amplitude of the effect is smaller than unity, then there is no point in even considering the covariance with other parameters.

We applied this procedure for both SKA1-MID and a cosmic variance-limited survey (${f}_{{\rm{sky}}}=1,{N}_{{\ell }}^{{ij}}=0$), obtaining the following result:

$$\begin{eqnarray}&&\mathrm{SKA}1-\mathrm{MID}\;\longrightarrow \;\sigma ({\epsilon }_{{\rm{GR}}})=2.75,\end{eqnarray} \tag{ 53 }$$

$$\begin{eqnarray}&&\mathrm{CV}-\mathrm{limited}\;\longrightarrow \;\sigma ({\epsilon }_{{\rm{GR}}})=1.97.\end{eqnarray} \tag{ 54 }$$

Thus, even in the best-case scenario, it is not possible to detect the effect of the relativistic corrections on the brightness temperature fluctuations.

This may seem like a striking result, as one of the relativistic effects is the equivalent of the CMB ISW effect for large-scale structure (${{\rm{\Delta }}}^{{\rm{ISW}}}$ in Equation (44)), and the CMB ISW has already been detected above $4\sigma$ by cross-correlating CMB maps with different LSS data sets (Giannantonio et al. 2008). In simplistic terms, intensity mapping surveys can be thought of as a set of uncorrelated "CMB" maps at different redshifts, so it is legitimate to ask why a similarly significant measurement is not possible in this case.

This can be explained in terms of the clustering variance of both data sets. Consider an attempt to measure the ISW effect by cross-correlating two data sets, one at high redshift (e.g., the CMB or a high-z H i intensity map), which we label here by a superscript h, and another at low redshift (e.g., a galaxy survey or low-z intensity map), which we label by g. Assuming that the ISW is the only term that could give rise to a significant cross-correlation between the two data sets, the signal would be given by the cross-power spectrum,

$$\begin{eqnarray}&&{S}_{{\rm{ISW}}}={C}_{{\ell }}^{{hg}}.\end{eqnarray} \tag{ 55 }$$

Neglecting any instrumental or shot noise, and assuming full-sky coverage and Gaussian statistics, the noise is purely given by the sample variance,

$$\begin{eqnarray}&&{N}_{{\rm{ISW}}}\simeq \sqrt{\displaystyle \frac{2}{2{\ell }+1}\;{C}_{{\ell }}^{{gg}}\;{C}_{{\ell }}^{{hh}}},\end{eqnarray} \tag{ 56 }$$

where we have assumed that ${C}_{{\ell }}^{{hg}}\ll \sqrt{{C}_{{\ell }}^{{gg}}\;{C}_{{\ell }}^{{hh}}}.$ Except for factors of order unity, the amplitude of the signal depends only on the value of $\phi ^{\prime} +\psi ^{\prime}$ at the redshift of g, and not on the nature of the high-redshift sample, so it will be roughly the same for both a high-z intensity mapping bin and the CMB. The difference in S/N between the two cases must therefore depend primarily on the amplitude of the noise, which differs by the ratio of ${C}_{{\ell }}^{{hh}}$ for the two cases. Since perturbations have grown significantly since ${z}_{{\rm{CMB}}}\sim 1100,$ we can expect the power spectrum of the intensity mapping to have a much larger amplitude, ${C}_{{\ell }}^{{hh},{\rm{IM}}}\gg {C}_{{\ell }}^{{hh},{\rm{CMB}}},$ which would explain the difficulty of achieving a good S/N in this case. This is shown explicitly in Figure 4: even at the highest redshift we considered, the power spectrum of the intensity mapping is four orders of magnitude larger than that of the CMB.

Zoom In Zoom Out Reset image size

Another effect conspiring against the detection of the relativistic terms in an intensity mapping survey is the fact that the perturbations on transverse scales cancel exactly (i.e., $s(z)\equiv 2/5$), as previously discussed. This further suppresses the overall amplitude of the relativistic effects, and is the reason why we do not present forecasts for ${\epsilon }_{{\rm{WL}}}$ in this case.

The uncertainty on the level of primordial non-Gaussianity measured by an intensity mapping experiment will depend on our prior knowledge of the free parameters of the model. Here we have imposed a (non-diagonal) Gaussian prior on the cosmological parameters $\{{{\rm{\Omega }}}_{M},{f}_{b},h,{w}_{0},{n}_{s},{A}_{s}\}$ using a prior covariance matrix estimated from the appropriate Planck 2015 MCMC chains for our set of parameters (Planck Collaboration 2015a). As we have argued, it is reasonable to assume that, by the time the SKA attempts to measure ${f}_{{\rm{NL}}}$, prior information will be available regarding the H i clustering bias (e.g., from experiments on smaller scales) and evolution bias (e.g., from external measurements of ${T}_{{\rm{H}}\;{\rm{I}}}(z)$).

Before we assume any specific priors for these parameters, it is worth studying their effect on $\sigma ({f}_{{\rm{NL}}})$ in order to quantify how good prior measurements will need to be in order to optimize the constraints on ${f}_{{\rm{NL}}}$. We first studied the level of degeneracy between ${f}_{{\rm{NL}}}$ and b(z) by estimating the value of $\sigma ({f}_{{\rm{NL}}})$ assuming a relative Gaussian prior on $b(z),$ constant across the whole redshift range. The result is shown in the left panel of Figure 5. We observe that a mild improvement on $\sigma ({f}_{{\rm{NL}}})$ ($\sim 10\%$ for SKA1-MID and $\sim 20\%$ for a CV-limited survey) can be achieved only for extremely accurate prior measurements of the clustering bias (${\rm{\Delta }}b/b\lesssim {10}^{-3}$). Since it would not be realistic to expect such a tremendous accuracy, we adopted a fiducial prior on b of $10\%,$ more compatible with current measurements of the bias of neutral hydrogen (Masui et al. 2013).

Zoom In Zoom Out Reset image size

Using this fiducial prior, we then explored the degeneracy between ${f}_{{\rm{NL}}}$ and ${f}_{{\rm{evo}}}$ by studying the dependence of $\sigma ({f}_{{\rm{NL}}})$ on the prior uncertainty, ${\rm{\Delta }}{f}_{{\rm{evo}}}$. The result is shown in the right panel of Figure 5. While $\sigma ({f}_{{\rm{NL}}})$ increases only slightly when factoring in the uncertainties on the clustering bias, we observe a much larger increase (e.g., by a factor 2.4 for SKA1-MID) when we assume no knowledge about the evolution bias of the sample at all. This suggests a much stronger degeneracy between ${f}_{{\rm{NL}}}$ and ${f}_{{\rm{evo}}}$, which could affect any attempt at measuring ${f}_{{\rm{NL}}}$ with LSS probes (note that a similar issue was reported by Camera et al. 2015b).

The source of the degeneracy can be understood by comparing the scale dependence of the terms in Equations (39)–(44) that are proportional to ${f}_{{\rm{evo}}}$ and ${f}_{{\rm{NL}}}$ respectively. Primordial non-Gaussianity introduces a term, included in ${{\rm{\Delta }}}_{{\ell }}^{{\rm{D}}},$ with a k-dependence of the form $\propto \delta ({\boldsymbol{k}}){j}_{{\ell }}(k\chi )/({k}^{2}\;T(k))\simeq \delta ({\boldsymbol{k}}){j}_{{\ell }}(k\chi )/{k}^{2},$ where the second equality holds on ultra-large scales. The evolution bias, on the other hand, multiplies four different terms:

• 1.  
  ${{\rm{\Delta }}}_{{\ell }}^{{\rm{P1}}},$ proportional to $\psi ({\boldsymbol{k}}){j}_{{\ell }}(k\chi )\propto \delta ({\boldsymbol{k}}){j}_{{\ell }}(k\chi )/{k}^{2}.$
• 2.  
  ${{\rm{\Delta }}}_{{\ell }}^{{\rm{V1}}},$ which is proportional to $\theta ({\boldsymbol{k}}){j}_{{\ell }}(k\chi )/{k}^{2}\propto \delta ({\boldsymbol{k}}){j}_{{\ell }}(k\chi )/{k}^{2}.$
• 3.  
  ${{\rm{\Delta }}}_{{\ell }}^{{\rm{V2}}},$ proportional to $\theta ({\boldsymbol{k}}){j}_{{\ell }}^{\prime }(k\chi )/k.$ For sufficiently large ℓ ($\gtrsim 5$), this term is also proportional to $\delta ({\boldsymbol{k}}){j}_{{\ell }}(k\chi )/{k}^{2}.$
• 4.  
  ${{\rm{\Delta }}}_{{\ell }}^{{\rm{ISW}}},$ the ISW term, which involves a much wider window function covering the full photon path from the source.

Thus, on large scales, three out of the four terms involving ${f}_{{\rm{evo}}}$ have the same scale dependence as the ${f}_{{\rm{NL}}}$ term, which explains the result found above. Fortunately, as can be seen in Figure 5, the effect of this degeneracy disappears if we can assume a relatively loose prior on ${f}_{{\rm{evo}}}$ of ${\rm{\Delta }}{f}_{{\rm{evo}}}\lesssim 1.$ That is, if we parameterize the evolution of background density of H i in the universe as ${\rho }_{{\rm{H}}\;{\rm{I}}}(a)\propto {a}^{\alpha },$ the slope α should then be measured with an error ${\rm{\Delta }}\alpha \lesssim 1.$ In what follows we assume that such an accuracy will be available from external measurements of ${T}_{{\rm{H}}\;{\rm{I}}}(z),$ and impose a Gaussian prior of ${\rm{\Delta }}{f}_{{\rm{evo}}}=1.$

For our fiducial set of priors (Planck CMB priors for the cosmological parameters, ${\rm{\Delta }}b/b=0.1,$ and ${\rm{\Delta }}{f}_{{\rm{evo}}}=1$), the final constraints on ${f}_{{\rm{NL}}}$ for SKA1-MID and for a cosmic variance-limited survey are given in Table 3. We also include results for an IM survey with a higher fiducial bias, which we discuss in Section 5.

Table 3.  Forecasted Constraints on ${f}_{{\rm{NL}}}$, ${\epsilon }_{{\rm{GR}}},$ and ${\epsilon }_{{\rm{WL}}}$.

Experiment Type	Experiment	$\sigma ({f}_{{\rm{NL}}})$	$\sigma ({\epsilon }_{{\rm{GR}}})$	$\sigma ({\epsilon }_{{\rm{WL}}})$
Intensity mapping	SKA1-MID	3.01	2.75	⋯
	(w. $1.5\times$ bias)	0.90	1.90	⋯
	CV-limited	1.68	1.97	⋯
Continuum survey	${S}_{{\rm{cut}}}=10\;\mu {\rm{Jy}}$	18.5	26.7	1.90
	${S}_{{\rm{cut}}}=5\;\mu {\rm{Jy}}$	16.0	24.6	0.57
	${S}_{{\rm{cut}}}=1\;\mu {\rm{Jy}}$	11.8	17.1	0.25
Spectroscopic survey	Hα survey	6.64	2.57	0.19
	CV-limited	3.02	1.35	0.10
Photometric survey	LSST-red	4.32	1.41	0.14
	LSST-full	1.71	2.33	0.04

Download table as:  ASCIITypeset image

One of the most important observational challenges for intensity mapping is the presence of galactic and extragalactic radio foregrounds (e.g., galactic synchrotron emission and extragalactic continuum radio sources) with amplitudes several orders of magnitude larger than the cosmological H i signal. The potential bias and extra variance induced in the measured signal by the process of foreground removal must be correctly taken into account in any analysis.

The strategy underlying most foreground cleaning methods is to exploit the very different frequency structures of signal and foregrounds. Most foreground signals have a very smooth frequency dependence, while the cosmological signal traces the stochastic fluctuations in the matter density, and is therefore much "noisier" in the radial (frequency) direction. Broadly speaking, most cleaning methods try to remove the foregrounds by fitting and subtracting a set of smooth functions of the frequency from the combined foreground + cosmological signal. See Alonso et al. (2015) for a description and comparison of different methods.

Since foregrounds are smooth in frequency, and frequency is a proxy for radial distances for an IM experiment, we can expect the foreground-cleaned maps to be dominated by systematics on large radial scales. These scales must then be omitted from the analysis, which reduces the sensitivity of an experiment to ${f}_{{\rm{NL}}}$. In order to understand the importance of this effect, we have studied the dependence of $\sigma ({f}_{{\rm{NL}}})$ on the maximum radial separation between redshift bins included in the computation of the Fisher matrix, i.e., we set to zero the off-diagonal elements of ${{\mathsf{C}}}_{{\ell }}$ corresponding to pairs of bins separated by more than some radial separation ${\rm{\Delta }}\chi .$

The results are shown in Figure 6 for SKA1-MID and for a cosmic variance-limited survey. Reducing the range of the cross-correlations included in the analysis can degrade the sensitivity to ${f}_{{\rm{NL}}}$ significantly, enlarging the errors by up to a factor ∼3.3 in the case of SKA1-MID.

Zoom In Zoom Out Reset image size

Given that measurements of ultra-large-scale modes will typically need to be done in single-dish (autocorrelation) mode, one also needs to worry about the stability of the instrument and the observational strategy. Much as in CMB experiments, long-term noise drifts (the "$1/f$" noise) will lead to striping in the maps, i.e., a coherent set of large-angle features that have been artificially projected on the sky. A sensible choice of scan strategy that leads to appreciable cross-linking between the scans can mitigate the effect, but there will always be a residual large-angle contaminant. Again, we can model this effect by not including the very large angular modes in the analysis (tantamount to assuming that they are filtered out by the destriping process). Figure 7 gives some idea of the impact of this effect on the constraints: a severe cut in the large-angle data significantly degrades any attempt to detect large-scale features.

Zoom In Zoom Out Reset image size

Relativistic effects in the continuum angular power spectrum were considered in Maartens et al. (2013) and Chen & Schwarz (2015). As in the previous section, we first present the most optimistic forecasts for ${\epsilon }_{{\rm{GR}}}$ by marginalizing only over ${f}_{{\rm{NL}}}$. If the uncertainty on ${\epsilon }_{{\rm{GR}}}$ is larger than unity (i.e., no detection of GR effects) in this optimistic case, there is no point in exploring more realistic scenarios.

There are two main differences between intensity mapping and continuum surveys in terms of the quantities that can affect these forecasts. First of all, we can expect the lack of radial information in continuum surveys to considerably degrade the constraints on most parameters. Perturbations to transverse scales will affect the observed clustering of radio sources in this case, however (i.e., $s(z)\ne 2/5,$ unlike for intensity mapping), which will enhance the amplitude of the relativistic terms.

The constraints on ${\epsilon }_{{\rm{GR}}}$ found for the three flux limits that we considered are shown in Table 3. Due to the lack of radial information, there is no hope of detecting GR effects using only a single-tracer analysis, in spite of the enormous volume probed; $\sigma ({\epsilon }_{{\rm{GR}}})=17.1$ even for the deepest ($1\;\mu {\rm{Jy}}$) survey.

We have also produced forecasts for the detectability of the weak lensing term on large scales for continuum surveys, following the same logic used in the case of ${\epsilon }_{{\rm{GR}}}$ (i.e., we keep all other parameters fixed, except for ${f}_{{\rm{NL}}}$). In order to pick up only the large-scale lensing contribution we also used a more stringent value of ${{\ell }}_{{\rm{max}}}=100.$ The results are summarized in Table 3: a continuum survey with a flux limit of $1\;\mu {\rm{Jy}}$ would be able to clearly detect the large-scale lensing effect above $\sim 4\sigma ,$ although the level of this detection would be below $2\sigma$ for ${S}_{{\rm{cut}}}=5\;\mu {\rm{Jy}}.$ No detection would be possible for a flux cut of $10\;\mu {\rm{Jy}}.$

In a continuum survey, the available information is compressed into only a small amount of data—the angular clustering statistics of radio sources—due to the lack of any sensitivity to radial modes. We can therefore expect an even larger degree of degeneracy between the various cosmological and nuisance parameters than for intensity mapping. As before, we assume Planck CMB priors on all cosmological parameters (except ${f}_{{\rm{NL}}}$), and start our discussion of the primordial non-Gaussianity forecasts by exploring the effect of prior information about the bias functions on the ${f}_{{\rm{NL}}}$ constraints. The results of this analysis are shown in Figure 8.

Zoom In Zoom Out Reset image size

The top panel of Figure 8 shows the forecast uncertainty on ${f}_{{\rm{NL}}}$ as a function of a constant relative Gaussian prior for b. In contrast with the situation for intensity mapping, the constraints are much more sensitive to the prior on the clustering bias. The main reason for this is that, for IM experiments, the availability of redshift information helps to break the degeneracy between the parameters through the scale dependence of the ${f}_{{\rm{NL}}}$ term along the line of sight. We find that a $\sim 10\%$ error on b(z) would be sufficient to minimize the uncertainty on ${f}_{{\rm{NL}}}$. Since this is compatible with previous smaller-scale observations (see Lindsay et al. 2014), we chose this value as our fiducial prior on the clustering bias.

The effect of a prior on the magnification bias is shown in the bottom left panel of Figure 8. We observe a similar degeneracy with ${f}_{{\rm{NL}}}$, again mainly due to the lack of redshift information, which can only be mitigated by prior measurements of s(z) with an error better than ${\rm{\Delta }}s\approx \pm 0.1.$ We have assumed that such an accuracy would be achievable using the magnitude–redshift distribution of sources, although this could be an optimistic assumption at the highest redshifts. Finally, as in the case of intensity mapping, we observe a significant degradation in the uncertainty on ${f}_{{\rm{NL}}}$ if we assume no knowledge about the evolution bias of the sample (bottom right panel in Figure 8). This is again due to the degeneracy in the scale dependence of the terms corresponding to both quantities (see Section 4.1.2), and can be mitigated by measuring ${f}_{{\rm{evo}}}$ to an accuracy of better than ${\rm{\Delta }}{f}_{{\rm{evo}}}\lesssim 10.$

In view of this analysis, our final fiducial set of priors on the bias functions is ${\rm{\Delta }}b/b=0.1,\;$ ${\rm{\Delta }}s=0.1,$ and ${\rm{\Delta }}{f}_{{\rm{evo}}}=1.$ Forecasts for ${f}_{{\rm{NL}}}$ with these priors are listed in Table 3 and, with $\sigma ({f}_{{\rm{NL}}})=16$ for the 5 μJy sample, are compatible with the results of Ferramacho et al. (2014) for their combined sample. This is far worse than the constraints possible with multiple tracers, however.

The different frequency range and observational techniques involved in radio astronomy give rise to potentially very different sources of systematics for continuum surveys compared with optical and near-infrared surveys. To begin with, the diffuse nature of galactic synchrotron emission (the largest galactic foreground at radio frequencies) makes it virtually transparent to the long interferometer baselines needed for a survey aiming to resolve individual sources, and hence the problem of galactic foregrounds is greatly ameliorated. On the other hand, in order to produce a full-sky catalog, mosaicing of the individual pointings must be implemented. If the mosaicing pattern and correlations in the noise properties between pointings are not fully understood, they could introduce systematic deviations on large angular scales. Ionospheric effects will also be relevant at low frequencies, although this should not be a problem in the SKA1-MID frequency range. Bright point sources would also need to be masked in a non-trivial way, due to dynamical range issues causing increased noise in the far beam sidelobes. The extent to which this would affect ultra-large scales is instrument-dependent, however. Avoiding these systematics might again entail removing the smallest multipoles of the power spectrum from the analysis; we show how the constraints on ${f}_{{\rm{NL}}}$ depend on the minimum multipole ${{\ell }}_{{\rm{min}}}$ in Figure 9.

Zoom In Zoom Out Reset image size

The most optimistic forecast for ${\epsilon }_{{\rm{GR}}},$ marginalizing over ${f}_{{\rm{NL}}}$ only, yields $\sigma ({\epsilon }_{{\rm{GR}}})=2.6$ for the Hα survey's narrowest redshift binning (${\rm{\Delta }}z=0.025$). This result is insensitive to bin width, increasing only slightly to 2.7 for the widest binning (${\rm{\Delta }}z=0.1$); the additional information gained by decreasing the bin width is mostly confined to small scales, where the relativistic effects are essentially negligible. The correlation between ${\epsilon }_{{\rm{GR}}}$ and ${f}_{{\rm{NL}}}$ is very weak, and there is no change in the constraint whether ${f}_{{\rm{NL}}}$ is marginalized or fixed. The survey is also quite close to its ideal (sample variance-limited) performance, with $\sigma ({\epsilon }_{{\rm{GR}}})$ improving only slightly to 2.2 in the limit ${N}_{{\ell }}\to 0$ for ${\rm{\Delta }}z=0.025.$ These results, together with the constraints corresponding to a cosmic variance-limited results (assuming ${N}_{{\ell }}=0$ and ${f}_{{\rm{sky}}}=1$) are summarized in Table 3, and have a qualitatively similar behavior.

As with the two previous surveys, then, the relativistic effects are undetectable. This is despite the relatively high magnification bias of the Hα galaxies, which boosts the size of some of the relativistic correction terms. The Hα survey's sky coverage and maximum redshift are smaller than for the other surveys though, which weakens its constraining power. An experiment with the same specifications as the Hα survey but covering twice the area (30,000 deg²) would give $\sigma ({\epsilon }_{{\rm{GR}}})=1.8$ (compared with 1.4 in the CV-limited case), which is still not enough to gain a detection.

Nevertheless, the Hα surveys's forecast constraint of $\sigma ({\epsilon }_{{\rm{GR}}})=2.6$ is the best so far, and the CV-limited figure of 1.36 is markedly better than SKA1-MID's value of 1.97, despite the IM survey having a significantly wider redshift range. The enhanced performance of the spectroscopic survey over intensity mapping is primarily due to the different behavior of the bias functions, particularly $s(z),$ which caused many of the relativistic effects to cancel for the IM survey.

As before, we also forecast for the detectability of the large-scale lensing effect, parameterized by ${\epsilon }_{{\rm{WL}}}.$ After restricting to modes ${\ell }\leqslant 100$ and marginalizing only over ${f}_{{\rm{NL}}}$, we find $\sigma ({\epsilon }_{{\rm{WL}}})=0.19$ for all three choices of redshift bin width—a strong detection (see Table 3).

As with the previous two probes, our ${f}_{{\rm{NL}}}$ forecasts include a Planck CMB prior and priors on the bias functions (${\rm{\Delta }}b/b=0.1,\;$ ${\rm{\Delta }}s=1,\;$ ${\rm{\Delta }}{f}_{{\rm{evo}}}=1$). The results are shown in Table 3 for the ${\rm{\Delta }}z=0.025$ redshift binning.

The forecast constraint from the Hα survey is $\sigma ({f}_{{\rm{NL}}})=6.8,$ which is worse than the intensity mapping survey by a factor of ∼2. While the Hα survey has a consistently higher bias (which enhances the non-Gaussian bias signal, $\propto b-1$), it covers a narrower redshift range and smaller area than the IM survey, so ultimately loses out when the higher-redshift bins of the IM survey are taken into account (see Figure 16, below). The difference in performance remains in the CV-limited case, again mostly due to the wider redshift range of the IM survey.

The Hα constraint degrades only slightly to $\sigma ({f}_{{\rm{NL}}})=7.2$ for the widest redshift binning, ${\rm{\Delta }}z=0.1.$ Similarly, it is only weakly sensitive to ${{\ell }}_{{\rm{max}}},$ improving from $\sigma ({f}_{{\rm{NL}}})=7.2$ for ${{\ell }}_{{\rm{max}}}=200$ to 6.6 for ${{\ell }}_{{\rm{max}}}=1000$ (both for ${\rm{\Delta }}z=0.025$). The addition of significantly more small-scale information in both the radial and transverse directions is therefore only mildly beneficial.

As with the IM survey, there is a reasonably strong degeneracy between ${f}_{{\rm{NL}}}$ and ${f}_{{\rm{evo}}}$, predominantly for the highest-redshift nuisance parameter bin. Figure 10 shows the effect of changing the prior on ${f}_{{\rm{evo}}}$—an ${\mathcal{O}}(1)$ prior is sufficient to completely break the degeneracy. The results are insensitive to the prior on the magnification bias, and there is no gain to be had from tightening the bias prior until a very low level of ${\rm{\Delta }}b/b\lesssim 1\%$ is reached.

Zoom In Zoom Out Reset image size

Spectroscopic surveys are prone to systematic effects on large angular scales. Redshift surveys commonly consist of samples from several non-contiguous fields, surveyed during different observing seasons and possibly even with different instruments, which can make it tricky to patch them together into a single coherent survey volume. They may also suffer from the problem of not having a homogeneous magnitude limit, i.e., the magnitude cuts vary, and cannot easily be mapped onto a substantial and complete 3D volume of the sky.

Galactic extinction is a dominant source of systematic error on large angular scales. Dust in our galaxy changes the overall true flux cut of the survey, and introduces number density fluctuations that vary with the shape of the Galaxy. This must be corrected for (e.g., by fitting an extinction template) to avoid biasing the inferred large-scale power. Variations in airmass and seeing also affect the number of photons reaching the detector in a way that is correlated with the elevation of the telescope. This induces additional dispersion in the magnitude of the measured galaxies, and the ability to distinguish them from stars.

Stars themselves are problematic. For example, a small fraction of the observed sources may in fact be misidentified stars that contaminate the galaxy sample. Stars also obscure regions of the sky of the order of the size of the point-spread function, which reduces the observed density of galaxies nearby. While this should be a small effect for an individual star, which will mask an area of $\sim {10}^{-6}$ of a degree, the total obscured area can be significant given the large stellar density (which grows toward the galactic plane). In fact, stellar contamination was found to be a dominant source of systematic error in the recent analysis of the BOSS data (Ross et al. 2011), where it introduced a significant bias in the measurement of the correlation function on large angular scales if left uncorrected. This bias was well above the statistical uncertainty, to the extent that the correlation functions measured in the northern and southern Galactic hemispheres of the survey were inconsistent with one other.

While many of these effects can lead to fluctuations in the number density of galaxies as a function of redshift, inducing systematics along the radial direction, the dominant effect is on large angular scales, reducing the effective area of the survey and hampering accurate recovery of the lowest ℓ modes. Figure 11 shows how the forecast constraints on ${\epsilon }_{{\rm{GR}}}$ and ${f}_{{\rm{NL}}}$ depend on the minimum recoverable ℓ mode of the survey. There is a rapid loss of information on both parameters as ${{\ell }}_{{\rm{min}}}$ increases, with $\sigma ({f}_{{\rm{NL}}})$ doubling from 6.8 for ${{\ell }}_{{\rm{min}}}=2$ to around 13 at ${{\ell }}_{{\rm{min}}}=10.$ The degradation is similar for $\sigma ({\epsilon }_{{\rm{GR}}}),$ which also doubles in the same range. As such, future spectroscopic surveys will likely need excellent control over large-scale systematics if they are to be used to constrain ${f}_{{\rm{NL}}}$.

Zoom In Zoom Out Reset image size

As with the previous experiments, we start by exploring the possibility of detecting the contribution of the relativistic terms to the clustering of LSST galaxies in the best-case scenario, by marginalizing over only ${f}_{{\rm{NL}}}$ while keeping all other cosmological parameters fixed to their fiducial values. Two main differences with respect to the previous tracers give some hope for detecting ${\epsilon }_{{\rm{GR}}}$ with LSST. First of all, the sharp decay in the number density of red galaxies can enhance the amplitude of the relativistic terms thanks to the large value of ${f}_{{\rm{evo}}}$. Also, LSST covers a wider survey area and redshift range than spectroscopic surveys, so has access to larger scales.

The results are summarized in Table 3. Even though the higher value of ${f}_{{\rm{evo}}}$ for red galaxies helps to decrease the forecast uncertainty on ${\epsilon }_{{\rm{GR}}},$ it is still impossible to detect relativistic effects using either sample; $\sigma ({\epsilon }_{{\rm{GR}}})=1.4$ and 2.3 for the red and full samples respectively. The red sample nevertheless produces the best constraint on ${\epsilon }_{{\rm{GR}}}$ of any of the surveys considered above.

As before, we predict the detectability of the large-scale magnification lensing by marginalizing only over ${f}_{{\rm{NL}}}$ and using a small-scale cutoff ${{\ell }}_{{\rm{max}}}=100.$ Our results (see Table 3) show that this effect should be clearly detectable (well above $5\sigma$) for both the "red" and "full" samples.

Even though photometric redshifts erase most of the clustering signal on all but the largest radial scales, they are still sufficient to enable a tomographic analysis of galaxy clustering to be performed. This helps immensely in breaking many of the degeneracies reported for continuum surveys (which retain essentially no radial information).

We studied the importance of breaking these degeneracies by again calculating $\sigma ({f}_{{\rm{NL}}})$ as a function of the priors on the bias parameters, finding that the forecasted uncertainty is almost completely insensitive to any priors on the clustering and magnification biases, b(z) and $s(z).$ For the reasons outlined in Section 4.1.2, this behavior does not follow for the evolution bias, so it is useful to explore the prior constraints on ${f}_{{\rm{evo}}}$ that are required in order to optimize the measurement of ${f}_{{\rm{NL}}}$. Figure 12 shows the dependence of $\sigma ({f}_{{\rm{NL}}})$ on a constant Gaussian prior imposed on ${f}_{{\rm{evo}}}$. The degeneracy between the parameters can be largely mitigated by measuring the evolution bias with an accuracy of ${\rm{\Delta }}{f}_{{\rm{evo}}}\lesssim 1.$

Zoom In Zoom Out Reset image size

As in the previous cases, we produced our final forecasts for ${f}_{{\rm{NL}}}$ by assuming Planck CMB priors for the cosmological parameters, a 10% uncertainty on the clustering bias, and priors of ${\rm{\Delta }}s=1$ and ${\rm{\Delta }}{f}_{{\rm{evo}}}=1.$ The final results are summarized in Table 3. LSST should be able to impose very tight constraints of $\sigma ({f}_{{\rm{NL}}})\simeq 1.7$ using galaxy clustering autocorrelations (single-tracer) alone.

Most of the sources of systematics that affect photometric redshift surveys are exactly the same as for their spectroscopic counterparts: galactic extinction, variations in sky brightness, seeing, and stellar contamination (due to both stars affecting the local observed number density of galaxies and stars erroneously being included in the galaxy sample). All of these effects can potentially contaminate the signal measured on large angular scales. Figure 13 shows the degradation in the constraints on ${f}_{{\rm{NL}}}$ when the largest scales are omitted in the analysis. Limiting ourselves to scales ${\ell }\geqslant 10$ would increase our best-case error bars by over 50%, to $\sigma ({f}_{{\rm{NL}}})\sim 2.8.$

Zoom In Zoom Out Reset image size

The use of photometric redshifts also carries its own systematic effects. In order to obtain a reliable estimate of the power spectrum of the galaxy density field that we can use to constrain large scales, it is necessary to have a sufficiently accurate model of the window function for every redshift bin. Unless a spectroscopic subsample drawn from the same distribution as the photometric one is available, which is rarely the case, deriving a correct model for the true redshift distribution $\bar{N}(z)$ is a challenging task, although it has been noted that this issue could be ameliorated by cross-correlating the photometric sample with any spectroscopic survey (Newman 2008). The presence of photometric redshift outliers can also modify the tails of the photo-z distribution, which affects the shape of the redshift window functions. The level to which this effect is problematic will depend on how accurately the photo-z probability density function (pdf) can be characterized.

For an IM autocorrelation experiment, the simplest case is to assume that the noise is uncorrelated between different frequency channels, and has a white noise power spectrum:

$$\begin{eqnarray}&&{N}_{{\ell }}^{{ij}}={\delta }_{{ij}}{\sigma }_{{\rm{sr}}}^{2}.\end{eqnarray} \tag{ 57 }$$

Here ${\sigma }_{{\rm{sr}}}^{2}$ is the noise variance per steradian, and can be calculated as follows: the noise per pointing can be estimated as the rms temperature fluctuation of the system, ${T}_{{\rm{sys}}},$ scaled by the number of independent samples measured (given by $\delta \nu \;{t}_{{\rm{p}}},$ where $\delta \nu$ is the frequency channel width and ${t}_{{\rm{p}}}$ is the integration time per pointing). ${t}_{{\rm{p}}}$ can be approximated by ${t}_{{\rm{tot}}}\;{\rm{\Delta }}{\rm{\Omega }}/(4\pi {f}_{{\rm{sky}}}),$ where ${f}_{{\rm{sky}}}$ is the surveyed fraction of the sky, ${\rm{\Delta }}{\rm{\Omega }}$ is the solid angle covered in each pointing, and ${t}_{{\rm{tot}}}$ is the total survey time. Finally, scaling this by the total number of dishes in the experiment, we obtain the power spectrum

$$\begin{eqnarray}&&{N}_{{\ell }}^{{ij}}={\delta }_{{ij}}\displaystyle \frac{{T}_{{\rm{sys}}}^{2}({\nu }_{i})4\pi \;{f}_{{\rm{sky}}}}{\delta \nu \;{t}_{{\rm{tot}}}\;{N}_{{\rm{dish}}}}.\end{eqnarray} \tag{ 58 }$$

Note that the variance per pointing has been multiplied by ${\rm{\Delta }}{\rm{\Omega }}$ to obtain the variance per steradian, which cancels the dependence on ${\rm{\Delta }}{\rm{\Omega }}.$ The system temperature receives two contributions, ${T}_{{\rm{sys}}}={T}_{{\rm{sky}}}+{T}_{{\rm{inst}}},$ due to atmospheric and background radio emission (${T}_{{\rm{sky}}}\simeq 60\;{\rm{K}}\times {(\nu /300\;{\rm{MHz}})}^{-2.5}$) and instrumental noise (${T}_{{\rm{inst}}}$).

Finally, it is worth noting that two different conventions have been adopted in the literature regarding the effect of the beam that defines the angular resolution of the experiment. The difference is in interpreting the beam as smoothing the signal on scales beyond the resolution, or as enhancing the noise at those same scales. We use the former, so that the model for the total observed power spectrum is

$$\begin{eqnarray}&&{C}_{{\ell }}^{{ij}}={C}_{{\ell }}^{S,{ij}}\;{B}_{{\ell }}^{i}\;{B}_{{\ell }}^{j}+{N}_{{\ell }}^{{ij}},\end{eqnarray} \tag{ 59 }$$

where ${B}_{{\ell }}^{i}$ is the harmonic transform of the instrumental beam in the ith frequency bin. We have assumed that the beams are Gaussian, ${B}_{{\ell }}^{i}=\mathrm{exp}(-{\ell }({\ell }+1){\theta }_{{\rm{B}}}^{2}/2),$ where ${\theta }_{{\rm{B}}}$ is related to the beam FWHM through ${\theta }_{{\rm{FWHM}}}=2\sqrt{2\mathrm{ln}2}\;{\theta }_{{\rm{B}}}.$ The beam width can be related to the dish diameter approximately as ${\theta }_{{\rm{FWHM}}}\simeq c/(\nu \;{D}_{{\rm{dish}}}).$

In our forecasts for SKA1-MID, we used the instrumental parameters ${T}_{{\rm{inst}}}=25\;K$, ${f}_{{\rm{sky}}}=0.75$, ${t}_{{\rm{total}}}={10}^{4}\;{\rm{h}}$, ${D}_{{\rm{dish}}}=15\;{\rm{m}},$ and ${N}_{{\rm{dish}}}=254,$ and assumed a minimum frequency of $350\;{\rm{MHz}},$ corresponding to ${z}_{{\rm{max}}}\simeq 3.$

Modern radio receivers have very high frequency resolution (e.g., $\delta \nu \sim 0.1\;{\rm{MHz}}$), so H i intensity mapping experiments should be able to resolve radial structures on scales much smaller than those relevant for cosmology. We are therefore free to choose the width and shape of the redshift bins used for the cosmological analysis. In order to avoid inhomogeneous coverage of radial scales, we divide the total frequency band into frequency bins of varying width ${\rm{\Delta }}\nu (\nu )$ such that the corresponding comoving size ${\rm{\Delta }}\chi$ is held constant. We estimated the minimum number of frequency bins needed for the constraints on ${f}_{{\rm{NL}}}$ and ${\epsilon }_{{\rm{GR}}}$ to converge, finding that at least 100 bins were necessary. This corresponds to a radial width of ${\rm{\Delta }}\chi \simeq 44\;{\rm{Mpc}}/h.$

We model the clustering, magnification, and evolution biases for H i using an approach based on the halo model. We first assume that a one-to-one relationship exists between halo mass and H i mass, ${M}_{{\rm{H}}\;{\rm{I}}}={M}_{{\rm{H}}\;{\rm{I}}}(M,z).$ The density and clustering bias can then be computed as

$$\begin{eqnarray}&&{\rho }_{{\rm{H}}\;{\rm{I}}}(z)={\displaystyle \int }_{{M}_{{\rm{min}}}}^{{M}_{{\rm{max}}}}{dM}\;n(M,z){M}_{{\rm{H}}\;{\rm{I}}}(M,z),\end{eqnarray} \tag{ 60 }$$

$$\begin{eqnarray}&&{b}_{{\rm{H}}\;{\rm{I}}}(z)={\displaystyle \int }_{{M}_{{\rm{min}}}}^{{M}_{{\rm{max}}}}{dM}\;n(M,z)b(M,z)\displaystyle \frac{{M}_{{\rm{H}}\;{\rm{I}}}(M,z)}{{\rho }_{{\rm{H}}\;{\rm{I}}}(z)},\end{eqnarray} \tag{ 61 }$$

where $n(M,z)$ is the halo mass function (comoving number density per unit mass) and $b(M,z)$ is the halo bias. The background brightness temperature can then be computed in terms of ${\rho }_{{\rm{H}}\;{\rm{I}}}$ using Equation (26).

As described in Section 2.1.2, the transverse distance perturbations cancel out for intensity mapping, so that ${s}_{{\rm{H}}\;{\rm{I}}}(z)=2/5$ exactly. Furthermore, since we observe the emission from all of the H i in each patch of the sky, the evolution bias can be computed directly by differentiating ${\rho }_{{\rm{H}}\;{\rm{I}}}$ with respect to z,

$$\begin{eqnarray}&&{f}_{{\rm{evo}}}=-\displaystyle \frac{d\;\mathrm{log}[{\rho }_{{\rm{H}}\;{\rm{I}}}(z)]}{d\;\mathrm{log}(1+z)}.\end{eqnarray} \tag{ 62 }$$

All that remains is to specify the function ${M}_{{\rm{H}}\;{\rm{I}}}(M,z).$ As in Bull et al. (2015), we assume a power-law relation ${M}_{{\rm{H}}\;{\rm{I}}}(M,z)\propto {M}^{\alpha }$ with an exponent $\alpha \simeq 0.6,$ and with the normalization set by constraints on ${{\rm{\Omega }}}_{{\rm{H}}\;{\rm{I}}}$ at z = 0.8 from Switzer et al. (2013).

The models for the signal and noise power spectra are very simple compared with the other probes, as the galaxy sample is distributed in a single redshift bin, with a window function given by the redshift distribution,

$$\begin{eqnarray}&&W(z)\propto \bar{N}(z).\end{eqnarray} \tag{ 63 }$$

We estimate $\bar{N}(z)$ for the radio sources from empirical estimates of the luminosity functions of the main radio populations, since these also contain the necessary information to estimate the magnification and evolution biases. We consider four main radio galaxy types: star-forming galaxies (SF), starbursts (SB), radio-quiet quasars (RQQ), and Faranoff–Riley type I active galactic nuclei (AGNs) (FRI).¹² The luminosity functions for each population were computed following the prescriptions of Wilman et al. (2008).

We will now outline the procedure used to calculate the redshift distribution, s, and ${f}_{{\rm{evo}}}$ in all cases, and refer the reader to Wilman et al. (2008) and the references Yun et al. (2001) (SF and SB), Ueda et al. (2003) (RQQ), and Willott et al. (2001) (FRI) for details on the observations that the luminosity functions are based on. These details are also summarized in Appendix C. As with optical and IR surveys (see Sections 4.3 and 4.4), the k-correction to the flux measured in a given band is also needed in order to accurately estimate the observed number counts. This can be done for radio sources by assuming a particular SED for each population. For this, we again used the models from Wilman et al. (2008).

In its rest frame, a radio source has a luminosity per unit frequency given by

$$\begin{eqnarray}&&{L}_{\nu }\equiv \displaystyle \frac{{{dE}}_{e}}{{{dt}}_{e}d{\nu }_{e}}={L}_{{\nu }_{*}}\displaystyle \frac{\varphi (\nu )}{\varphi ({\nu }_{*})},\end{eqnarray} \tag{ 64 }$$

where ${\nu }_{*}$ is a pivot frequency and $\varphi (\nu )$ is the source SED. This is related to the flux per unit frequency measured by the observer by

$$\begin{eqnarray}&&{S}_{\nu }\equiv \displaystyle \frac{{{dE}}_{o}}{{{dt}}_{o}d{\nu }_{o}{{dA}}_{o}}=\displaystyle \frac{{L}_{\nu (1+z)}}{4\pi {\chi }^{2}(z)(1+z)}(1-2{\delta }_{\perp }).\end{eqnarray} \tag{ 65 }$$

Consider a radio survey in a given frequency band ($\nu \in [{\nu }_{0},{\nu }_{f}]$). The average flux density measured for one source is defined as

$$\begin{eqnarray}&&\bar{S}({\nu }_{0},{\nu }_{f})={\displaystyle \int }_{{\nu }_{0}}^{{\nu }_{f}}{S}_{\nu }\displaystyle \frac{d\nu }{{\nu }_{f}-{\nu }_{0}}.\end{eqnarray} \tag{ 66 }$$

A given source will be detected if its average flux is above the detection limit ${S}_{{\rm{cut}}},$ and therefore all sources with a pivot luminosity above a minimum value ${L}_{{\nu }_{*}}^{{\rm{cut}}}\equiv {\bar{L}}_{{\nu }_{*}}^{{\rm{cut}}}(1+2{\delta }_{\perp })$ will be included in the sample, where the average threshold luminosity is

$$\begin{eqnarray}&&{\bar{L}}_{{\nu }_{*}}^{{\rm{cut}}}(z,{S}_{{\rm{cut}}})=\displaystyle \frac{4\pi {S}_{{\rm{cut}}}{\chi }^{2}(z)(1+z)\varphi ({\nu }_{*})}{\bar{\varphi }({\nu }_{0}(1+z),{\nu }_{f}(1+z))},\end{eqnarray} \tag{ 67 }$$

and the average SED in the observed band is

$$\begin{eqnarray}&&\bar{\varphi }({\nu }_{1},{\nu }_{2})\equiv {\displaystyle \int }_{{\nu }_{1}}^{{\nu }_{2}}\varphi (\nu )\displaystyle \frac{d\nu }{{\nu }_{2}-{\nu }_{1}}.\end{eqnarray} \tag{ 68 }$$

Given a pivot frequency ${\nu }_{*},$ the luminosity function ${\bar{n}}_{s}(z,\mathrm{ln}{L}_{{\nu }_{*}})$ at that frequency, and a characteristic SED $\varphi (\nu ),$ the redshift distribution of sources can be computed as

$$\begin{eqnarray}&&\bar{N}(z)=\displaystyle \frac{c\;{\chi }^{2}(z)}{{(1+z)}^{3}\;H(z)}\bar{{\mathcal{N}}}(z,{\gt \bar{L}}_{{\nu }_{*}}^{{\rm{cut}}}),\end{eqnarray} \tag{ 69 }$$

where $\bar{{\mathcal{N}}}$ is defined as in Equation (8). s(z) and ${f}_{{\rm{evo}}}(z)$ are then calculated from $\bar{{\mathcal{N}}}$ using Equations (28) and (29).

For the clustering bias, we follow the same approach used in Wilman et al. (2008) and assign a fixed halo mass to each population. The corresponding bias is then found as the halo-model bias for that mass as a function of redshift. For this we parameterize the halo bias as in Sheth & Tormen (1999). We have also explored the approach followed in Ferramacho et al. (2014), where each population is given a distribution of halo masses rather than a fixed one, and the bias is found by averaging over that distribution. No significant differences were found between the approaches, so we use the first, simpler one.

Once the redshift distribution and bias functions have been calculated for each population, we compute them for the combined sample as a weighed average of the individual ones,

$$\begin{eqnarray}&&{\bar{N}}_{{\rm{tot}}}(z)=\displaystyle \sum _{a}{\bar{N}}_{a}(z),\end{eqnarray} \tag{ 70 }$$

$$\begin{eqnarray}&&{(b,s,{f}_{{\rm{evo}}})}_{{\rm{tot}}}(z)=\displaystyle \sum _{a}{(b,s,{f}_{{\rm{evo}}})}_{a}(z)\displaystyle \frac{{\bar{N}}_{a}(z)}{{\bar{N}}_{{\rm{tot}}}(z)},\end{eqnarray} \tag{ 71 }$$

where a labels the population. Separating the different populations is a very costly observational task. For our purposes, the main benefit of doing this is to allow the use of the multi-tracer technique to circumvent cosmic variance. Since we have postponed the multi-tracer analysis for future work, we will only report our forecasts here for the combined sample of radio sources. Figure 18 shows the redshift distribution for the combined sample for the three different detection limits considered here. The luminosity functions, redshift distributions, and bias can be obtained using a code that we have made publicly available.¹³ This provides an easy way to obtain number counts and power spectra without needing to query the full simulation of Wilman et al. (2008). Our results are consistent with the simulation except for the total number counts of star-forming galaxies, where our numbers are a factor of $2.5\times$ higher. This is consistent with what is described in Jarvis et al. (2015) and references therein, however.

Zoom In Zoom Out Reset image size

The most relevant source of noise in a clustering analysis of discrete sources is shot noise, determined by the number density of sources. The angular number density of radio sources in the sample, ${\bar{N}}_{{\rm{\Omega }}},$ is determined by integrating the redshift distribution

$$\begin{eqnarray}&&{\bar{N}}_{{\rm{\Omega }}}\equiv {\displaystyle \int }_{0}^{\infty }{dz}\;\bar{N}(z),\end{eqnarray} \tag{ 72 }$$

and the noise power spectrum is given by

$$\begin{eqnarray}&&{N}_{{\ell }}={\bar{N}}_{{\rm{\Omega }}}^{-1}.\end{eqnarray} \tag{ 73 }$$

Spectroscopic galaxy surveys measure Equation (11), the perturbation to flux-limited number counts ${{\rm{\Delta }}}_{{\rm{N}}}(z,\hat{{\boldsymbol{n}}},\gt \mathrm{log}{F}_{{\rm{cut}}}).$ This depends on the clustering, magnification, and evolution bias functions for the source population. The latter two can be derived from the background luminosity function of the sources, ${n}_{s}(z,\mathrm{log}L),$ given the flux limit and efficiency of the survey. Number counts are subject to a Poisson shot noise term that depends on the number density of sources, which can also be obtained from the luminosity function.

We follow the current set of public specifications for a realistic Hα spectroscopic survey, described in Amendola et al. (2013). The number counts presented in Table 1.3 of that paper can be approximately reproduced by using the best-fit Hα Schechter luminosity function found by Geach et al. (2010),

$$\begin{eqnarray}&&\begin{array}{l}n(z,L\gt {L}_{{\rm{cut}}})={\displaystyle \int }_{{x}_{{\rm{cut}}}}^{\infty }\epsilon \;{\varphi }^{*}\;{x}^{\alpha }{e}^{-x}\;{dx};\\ \quad x\equiv L{/L}^{*}(z),\\ {L}^{*}(z)=\;\left\{\begin{array}{l}5.1\times {10}^{41}{(1+z)}^{3.1}\;\;\;(z\lt 1.3)\\ 6.8\times {10}^{42}\;\;\;(1.3\lt z\lt 2.2)\end{array}\right.\mathrm{erg}\;{{\rm{s}}}^{-1},\end{array}\end{eqnarray}$$

where L is the Hα line luminosity and we have assumed a flux limit of ${F}_{{\rm{cut}}}=3\times {10}^{-16}$ erg s⁻¹ cm⁻², an efficiency of $\epsilon =0.45,$ a faint-end slope of $\alpha =-1.35,$ and comoving number density normalization ${\phi }^{*}=1.37\times {10}^{-3}$ Mpc⁻³. The predicted galaxy number density as a function of redshift is shown in Figure 19. There are significant uncertainties in this model, which we account for by marginalizing over a set of bias function nuisance parameters, as explained in Section 3. Finally, for the clustering bias we use the simplified prescription from Amendola et al. (2013), $b(z)=\sqrt{1+z},$ which we also subject to the nuisance parameterization.

Zoom In Zoom Out Reset image size

The assumed bias functions for the Hα survey were shown in Figure 2. The evolution bias is large and negative, and grows more negative with redshift until $z=1.3,$ where there is a discontinuity in the ${L}^{*}(z)$ model, beyond which ${L}^{*}(z)={\rm{const.}}$ and so ${f}_{{\rm{evo}}}=0.$ This is not particularly realistic, but is the best that can be done until updated constraints on the high-redshift Hα luminosity function become available. The magnification bias grows rapidly with redshift, and deviates significantly from s = 0.4 over most of the range, meaning that there is little or no cancellation of the transverse scale perturbation as there was in the intensity mapping survey (although this is contingent on the uncertain behavior of the luminosity function model for $z\gt 1.3$). The clustering bias is relatively close to unity across the entire redshift range, slowly evolving from a minimum of $b\approx 1.3$ at z = 0.65 to a maximum of $b\approx 1.8$ at $z\approx 2.$

We consider three constant-width redshift binnings for the selection function over the interval $0.65\leqslant z\leqslant 2.05:$ ${\rm{\Delta }}z=(0.025,0.05,0.10),$ yielding (56, 28, 14) bins respectively. The target redshift uncertainty for, e.g., Euclid is ${\sigma }_{z}\leqslant 0.001\;(1+z)$ (Laureijs et al. 2011), which is always significantly smaller than the narrowest bin width. As such, we assume a uniform (top-hat) selection function, weighted by the source redshift distribution, $\bar{N}(z).$

Our noise model assumes that only shot noise is relevant, and that effects such as spectroscopy failures and point source masking have been taken into account in the survey efficiency, . The shot noise angular power spectrum is

$$\begin{eqnarray}&&{N}_{{\ell }}^{{ij}}={\delta }_{{ij}}{/n}_{i};\;\;\;{n}_{i}\equiv {\displaystyle \int }_{{z}_{i}}\bar{N}(z){dz}.\end{eqnarray} \tag{ 74 }$$

with n_i measured in units of sr⁻¹. In our forecasts, we assumed that a wide range of multipoles can be recovered, $2\leqslant {\ell }\leqslant 1000,$ with no cuts at high or low ℓ due to systematics, nonlinear effects, and so on (this assumption was relaxed in Section 4.3.3).

In order to compute $\bar{N}(z),\;$ $s(z),$ and ${f}_{{\rm{evo}}}(z)$ for our two samples, we need an estimate of the luminosity function for both red and blue galaxies, preferably in the r-band, for which the LSST specifications are provided. We describe the method used for this task here.

For the red sample we follow a method similar to that used by Joachimi et al. (2011). First, an estimate of the B-band luminosity function for red galaxies is obtained from Faber et al. (2007) as a Schechter function with constant slope $\alpha =-0.5,$ and z-dependent ${\phi }_{*}$ and ${M}_{*},$ measured in a number of redshift bins in the interval $z\in (0.2,1.2).$ We extrapolate the luminosity function to higher/lower redshifts by fitting the values of these parameters, measured by Faber et al. (2007), to the models

$$\begin{eqnarray}&&{M}_{*}(z)={M}_{0}+{M}_{1}\;z,\end{eqnarray} \tag{ 75 }$$

$$\begin{eqnarray}&&{\phi }_{*}(z)=\displaystyle \frac{{\phi }_{0}}{1+{(z/{z}_{0})}^{a}}\;[{10}^{-3}\;{{\rm{Mpc}}}^{-3}]\end{eqnarray} \tag{ 76 }$$

with ${M}_{0}=-20.6,\;$ ${M}_{1}=-0.49,\;$ ${\phi }_{0}=1.82,\;$z₀ = 1.04, and a = 7.17. In order to translate this into an r-band luminosity function, we use $B-r=1.32$ (Fukugita et al. 1995) and assume that B − r does not evolve significantly for the red sample in the redshift range under study.

For the full sample, we use the $r^{\prime}$-band luminosity function found by Gabasch et al. (2006) and approximate $r^{\prime} \simeq r.$ This is again given as a Schechter function with constant slope $\alpha =-1.33$ and z-dependent M_* and ${\phi }_{*},$ for which we have used the following parameterizations:

$$\begin{eqnarray}&&{M}_{*}(z)={M}_{0}+a\;\mathrm{ln}(1+z),\end{eqnarray} \tag{ 77 }$$

$$\begin{eqnarray}&&{\phi }_{*}(z)=({\phi }_{0}+{\phi }_{1}\;z+{\phi }_{2}\;{z}^{2})[{10}^{-3}\;{{\rm{Mpc}}}^{-3}]\end{eqnarray} \tag{ 78 }$$

with ${M}_{0}=-21.49,\;$ $a=-1.25,\;$ ${\phi }_{0}=2.59,\;$ ${\phi }_{1}=-0.136,$ and ${\phi }_{2}=-0.081.$ The luminosity function for blue galaxies is then estimated as the difference between those of the full and red samples.

An absolute magnitude, M, measured in a given rest-frame band for a galaxy at redshift z is related to the apparent magnitude in the observer-frame band, m, by

$$\begin{eqnarray}M & = & m-25-5\;{\mathrm{log}}_{10}\left[\displaystyle \frac{{d}_{L}(z)}{1\;{\rm{Mpc}}\;{h}^{-1}}\right]\\ & & +{\mathrm{log}}_{10}h-k(z),\end{eqnarray} \tag{ 79 }$$

where d_L(z) is the luminosity distance and k(z) is the k-correction corresponding to that galaxy's SED redshifted to z. We estimated k(z) for red and blue galaxies by running the code kcorrect (Blanton & Roweis 2007) on the spectra of an elliptical galaxy and a barred spiral galaxy (Sbc) respectively, as measured by Coleman et al. (1980). We approximate and extrapolate these k-corrections as ${k}_{{\rm{red}}}(z)\sim 2.5\;z$ and ${k}_{{\rm{blue}}}\sim 1.5\;z.$ We verified that these parameterizations are compatible with the k-corrections shown in Fukugita et al. (1995) for both types, and also that the final redshift distributions did not vary significantly when other functional forms are used. For a given magnitude limit ${m}_{{\rm{cut}}},$ we use Equation (79) to compute the corresponding luminosity cut ${\bar{L}}_{{\rm{cut}}}$ as a function of redshift. The redshift distribution, magnification bias, and evolution bias for each population are then estimated using Equations (28)–(30), respectively, and with the luminosity functions described above.

Regarding the clustering bias, for the full sample we use the parameterization ${b}_{{\rm{full}}}(z)=1+0.84\;z,$ based on the simulations of Weinberg et al. (2004) and quoted in the LSST science book (LSST Collaboration et al. 2009). Red galaxies should have a larger bias, which we parameterize as ${b}_{{\rm{red}}}(z)\simeq 1+z.$ This parameterization is compatible with bias measurements at redshifts $z\lt 1$ (e.g., Coil et al. 2008). Since the red population dies off at $z\sim 1.4,$ this extrapolation should not significantly influence our final result, especially as we ultimately marginalize over $b(z).$

The redshift distributions for the two samples considered here are shown in Figure 20. For this figure, as well as in our forecasts, we assume a magnitude limit of $i=25.3,$ corresponding to $r\sim 26$ for typical galaxy colors, as quoted in LSST Collaboration et al. (2009). According to the models used here, the total number of galaxies observed by LSST should be ∼40 per arcmin², in qualitative agreement with previous results (Ilbert et al. 2006; LSST Collaboration et al. 2009).

Zoom In Zoom Out Reset image size

The quality of a given photo-z algorithm is normally quoted in terms of its rms error, ${\sigma }_{z}^{2}\equiv \langle {({z}_{{\rm{photo}}}-{z}_{{\rm{true}}})}^{2}\rangle ,$ which is typically parameterized as

$$\begin{eqnarray}&&{\sigma }_{z}(z)={\sigma }_{0}(1+z).\end{eqnarray} \tag{ 80 }$$

The photometric redshift requirement for the LSST gold sample, as quoted in LSST Collaboration et al. (2009), is ${\sigma }_{0}\lt 0.05,$ with a goal of 0.02. As described above, since the spectral properties of red galaxies make their photometric redshifts more accurate than those of the blue population, we have assumed the following photometric redshift uncertainties for the two samples:

$$\begin{eqnarray}&&{\sigma }_{0}^{{\rm{red}}}=0.02,{\sigma }_{0}^{{\rm{full}}}=0.05.\end{eqnarray} \tag{ 81 }$$

We have also assumed that the data will be analyzed by dividing the sample into a number of photo-z bins. Let z₀ⁱ and z_fⁱ be the limits of the ith bin. The window function in this bin must trace the true z distribution of galaxies within it, and is therefore given by the product of the overall redshift distribution and the photo-z probability distribution, integrated over the bin:

$$\begin{eqnarray}&&{W}^{i}(z)\propto \bar{N}(z){w}^{i}(z),\end{eqnarray} \tag{ 82 }$$

$$\begin{eqnarray}&&{w}^{i}(z)={\displaystyle \int }_{{z}_{0}^{i}}^{{z}_{f}^{i}}{{dz}}_{p}\;p({z}_{p}| z),\end{eqnarray} \tag{ 83 }$$

where $p({z}_{p}| z)$ is the photo-z pdf. Assuming a Gaussian photo-z distribution, we can write the window functions ${w}^{i}(z)$ analytically as

$$\begin{eqnarray}&&{w}^{i}(z)=\displaystyle \frac{1}{2}\left[{\rm{erf}}\left(\displaystyle \frac{z-{z}_{0}^{i}}{\sqrt{2}{\sigma }_{z}}\right)-{\rm{erf}}\left(\displaystyle \frac{z-{z}_{f}^{i}}{\sqrt{2}{\sigma }_{z}}\right)\right].\end{eqnarray} \tag{ 84 }$$

The tails of ${w}^{i}(z)$ correlate different redshift bins to a much larger degree than the intrinsic correlations due to gravitational clustering do, which is an expression of the loss of information on radial scales. In order to reduce this correlation and avoid redundant calculations, the widths of the redshift bins are usually defined to be of the order of ${\sigma }_{z}.$ In this work we have chosen to define the width of our bins to be three times the photo-z dispersion at the bin center. For the values of ${\sigma }_{z}$ assumed here, this results in 15 bins for the red sample and 9 bins for the full sample. The window functions for the bins are shown in Figure 20.

Finally, as in the case of spectroscopic and continuum surveys, the main source of statistical noise in the measurement of clustering anisotropies is shot noise. The noise bias term for photometric surveys is thus given by

$$\begin{eqnarray}&&{N}_{{\ell }}^{{ij}}=\displaystyle \frac{{\delta }_{{ij}}}{{n}^{i}\;{n}^{j}},\end{eqnarray} \tag{ 85 }$$

where

$$\begin{eqnarray}&&{n}^{i}\equiv {\displaystyle \int }_{0}^{\infty }{dz}\;\bar{N}(z){w}^{i}(z),\end{eqnarray} \tag{ 86 }$$

Note that this reduces to Equation (74) in the case of top-hat windows (i.e., ${\sigma }_{z}\to 0$).