# Fisher matrix for the angular power spectrum of multi-tracer galaxy surveys

L. Raul Abramo  , João Vitor Dinarte Ferri  , Ian Lucas Tashiro  and Arthur
Loureiro

###### Abstract

Redshift evolution and peculiar velocities break the isotropy of cosmological
surveys with respect to the directions parallel and transverse to the line of
sight, limiting the accuracy of the Fourier representation to small areas and
redshift ranges. In contrast to the Fourier space power spectrum, the full
information about the two-point function of tracers of large-scale structure
is encapsulated in the redshift-dependent angular power spectrum
$C_{\ell}^{ij}(z_{i},z_{j})$ for the tracer species $i$ and $j$ at the
redshift slices $z_{i}$ and $z_{j}$, expressed in harmonic space. In this
paper we derive semi-analytical expressions for the multi-tracer Fisher matrix
of angular power spectra, in real and in redshift space, which are exact in
the linear regime of structure formation. Our expressions can be used to
forecast the constraining power of galaxy surveys with many tracers and a
large number of redshift slices, for which the derivation of the Fisher matrix
from numerically evaluated covariance matrices may not be feasible or
practical.

##  1 Introduction

Small-area, shallow astrophysical surveys are increasingly giving way to
large-area, deep surveys that can map the distribution of galaxies and other
tracers of large-scale structure over large fractions of our past light-cone.
This process, which started with the Sloan Digital Sky Survey [1], is now
being driven by surveys such as DESI [2], PFS [3], Euclid [4], the Vera Rubin
Observatory Legacy Survey of Space and Time [5], J-PAS [6], and SKA [7]. These
surveys will go beyond just measuring the scale of baryon acoustic
oscillations, providing exciting new tests of modified gravity [8, 9],
primordial non-Gaussianities [10], relativistic effects [11, 12, 13, 14], as
well as disentangling the cosmological mechanisms that operate at ultra-large
scales [15, 16, 17].

It is therefore critical that we understand to what extent future surveys will
be able to capture the information contained on these largest observable
scales, both in terms of the vast redshift ranges that are now accessible to
our instruments, but also in terms of large angular separations that will be
covered by those surveys. This means, in particular, that a traditional
approach based on the Fourier-space power spectrum cannot possibly provide an
accurate description: in its simplest form, it assumes that the sky is flat,
which is equivalent to assuming that the survey is inscribed by a Cartesian
box placed at a large distance from the observer. More sophisticated
approximation schemes that try to take into account the curvature of the sky
are able to partially offset this problem [18, 19, 20], but they also start to
fail at intermediate scales. We are then led almost inevitably to a spherical
harmonic description of the matter distribution in redshift space [21, 22, 23,
24, 25], which should retain information from large angular scales and all
redshift ranges, as well as taking into account the many different tracers of
large-scale structure: galaxies [26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36,
37], quasars [38], as well as radio and HI intensity mapping [39, 40]. The
description of weak lensing in terms of a spherical harmonic decomposition and
angular power spectra has also become increasingly popular [41, 42, 43, 44].

One of the observables that can encapsulate all the information of our past
light-cone is the angular power spectrum $C^{ij}_{\ell}(x_{i},x_{j})$, which
is the harmonic component ($\ell$) of the clustering between the tracer $i$ at
the redshift slice $z_{i}$ at radius $x_{i}=\chi(z_{i})$ and the tracer $j$ at
the redshift slice $z_{j}$ and radius $x_{j}=\chi(z_{j})$. However, given the
large number of degrees of freedom in these angular power spectra, it can be
very difficult to estimate the constraining power of large surveys on the
basis of these observables. Despite recent breakthroughs in numerical
implementations of the harmonic approach [45, 46, 47, 48], analytic and semi-
analytic tools are still needed. In particular, spectroscopic and intensity
mapping surveys allow us to split the observations into a large number of thin
redshift slices, which then become increasingly correlated. By using thinner
redshift slices we avoid smearing the maps of the tracers over the radial
direction, but this creates a problem insofar as the covariance matrix for the
angular power spectra grows fast with the number of slices: e.g., even for a
single tracer, splitting a survey into 50 redshift slices leads to a
$(1275)^{2}$ covariance matrix for each multipole $\ell$. Computing such an
object with sufficient accuracy may be fraught with difficulties, especially
given that what we actually need is either the inverse of that covariance
(i.e., the Fisher matrix) or its Cholesky decomposition. From a different
perspective, with the increasing popularity of Bayesian Hierarchical Models
for Field-Level Inference[49, 50, 51, 52, 53, 54, 55, 56], analytic
expressions for the covariance matrix and its inverse have an important
application as the mass matrix or step-size proxies in a Hamiltonian Monte-
Carlo framework [57, 58], allowing for a significant increase in efficiency
for these samplers. As is the case for any Fisher analysis, we do assume
Gaussianity (since that is what enables us to write the covariance of the
power spectra in terms of a 4-point function), but our semi-analytical results
are able to capture the full time evolution along the past light-cone. In this
paper we derive a semi-analytical expression for the multi-tracer Fisher
matrix of the angular power spectra. Our results will pave the way to
forecasting accurately the constraining power of surveys on ultra-large
scales, and to estimate the uncertainties in the angular power spectra. In
particular, since we do not employ the Limber approximation [59], we are able
to obtain an analytical expression for the radial mixing matrix that
determines the coupling between radial modes in the angular power spectra.
Finally, we should stress that in this paper we do not resort (at least not
explicitly) to the spherical Fourier-Bessel transform of the density field
[23, 60, 45, 48].

This paper is organized as follows. In Section 2 we review general expressions
for the covariance and Fisher matrices, including the case when there are
multiple tracers, and including their cross-correlations – always under the
assumption of Gaussianity. In Section 3 we show how to write the data
covariance (i.e., the observed angular power spectra, including shot noise)
and how it is possible to invert that object – which is the key step that
allows us to compute the Fisher matrix. Still in Section 3 we present an exact
solution for a toy model, that also serves as a check for our formulae. In
Section 4 we obtain the data covariance in redshift space, and present a semi-
analytical expression for its inverse – the Fisher matrix. At the end of that
Section we also discuss in which circumstances one would employ our semi-
analytical methods, as opposed to direct numerical computations in terms of
mocks or simulations. We conclude with a discussion in Section 5. For the sake
of clarity and readability, some technical details are left for the three
appendices.

##  2 Fisher matrix: general expressions

The fundamental degrees of freedom in a survey are the positions of the
tracers (galaxies, halos or other point-like objects that follow the
underlying matter distribution). When we measure the number densities
$n^{i}(\vec{x})$ of a tracer species $i$, over some volume around the position
$\vec{x}$, that number reflects the mean density of those tracers,
$\bar{n}^{i}(\vec{x})$, as well as the fluctuations $\delta
n^{i}=n^{i}-\bar{n}^{i}$. From these observables we compute the main object
that carries information about cosmology, the data (or “pixel”) covariance:

$\Gamma^{ij}(\vec{x},\vec{y})=\langle\delta n^{i}(\vec{x})\,\delta
n^{j}(\vec{y})\rangle=\bar{n}^{i}(\vec{x})\,\bar{n}^{j}(\vec{y})\,\xi^{ij}(\vec{x},\vec{y})+\bar{n}^{i}(\vec{x})\delta^{ij}\delta_{D}(\vec{x}-\vec{y})\;,$
(2.1)

where $\xi^{ij}(\vec{x},\vec{y})$ is the 2-point correlation function in
configuration space, and the last term is shot noise, which we assume here to
follow Poisson statistics. The multi-tracer 2-point correlation function is
generally assumed to be related to the matter correlation function,
$\xi^{(m)}(\vec{x},\vec{y})$, through some knowable relations such as tracer
bias, redshift-space distortions [61], etc. In real space (i.e., excluding
redshift-space distortion), the matter two-point correlation function can be
written in terms of the matter power spectrum as:

$\xi^{(m)}(\vec{x},\vec{y})=\xi^{(m)}(|\vec{x}-\vec{y}|)=\int\frac{d^{3}k}{(2\pi)^{3}}\,e^{-i\vec{k}\cdot(\vec{x}-\vec{y})}\,P^{(m)}(k)\;.$
(2.2)

Let’s say we want to estimate the constraints on some set of parameters
$\theta^{\mu}$, which can be derived from the two-point correlation functions.
We start by constructing the Fisher information matrix for that set of
parameters [62]:

$\displaystyle F[\theta^{\mu},\theta^{\nu}]=F^{\mu\nu}$ $\displaystyle=$
$\displaystyle\frac{1}{2}{\rm
Tr}\left[\frac{\partial\Gamma}{\partial\theta^{\mu}}\Gamma^{-1}\frac{\partial\Gamma}{\partial\theta^{\nu}}\Gamma^{-1}\right]$
$\displaystyle=$ $\displaystyle\frac{1}{2}\sum_{iji^{\prime}j^{\prime}}\int
d^{3}x_{i}\int d^{3}x_{j}\int d^{3}x_{i^{\prime}}\int d^{3}x_{j^{\prime}}$
$\displaystyle\times\left\\{\frac{\partial\Gamma^{ij}(\vec{x}_{i},\vec{x}_{j})}{\partial\theta^{\mu}}[\Gamma^{ji^{\prime}}(\vec{x}_{j},\vec{x}_{i^{\prime}})]^{-1}\frac{\partial\Gamma^{i^{\prime}j^{\prime}}(\vec{x}_{i^{\prime}},\vec{x}_{j^{\prime}})}{\partial\theta^{\nu}}[\Gamma^{j^{\prime}i}(\vec{x}_{j^{\prime}},\vec{x}_{i})]^{-1}\right\\}\;,$

where the compact notation for the trace in the definition of the first line
denotes both a sum over all tracers, as well as spatial integrations over all
positions. This expression also gives us the first hints that the tracer
indices appear linked to the spatial positions where we measure those same
tracers: $i\leftrightarrow\vec{x}_{i}$, $j\leftrightarrow\vec{x}_{j}$, and so
on and so forth. This means that we can make use, whenever convenient, of a
compact notation whereby the sum over a tracer index also denotes an integral
over its corresponding positions in space.

Finally, after computing the Fisher matrix we can invert it to find the
covariance matrix for the parameters $\theta^{\mu}$:

${\rm Cov}[\theta^{\mu},\theta^{\nu}]={\rm
Cov}_{\mu\nu}=\left[F_{\mu\nu}\right]^{-1}\;,$ (2.4)

which is an estimate of the marginalized uncertainties and degeneracies
between the parameters.

###  2.1 Fisher matrix in Fourier space

Working directly in Fourier space allows us to introduce some notions and
notation that will become useful later when we compute the multi-tracer Fisher
matrix in harmonic space. For simplicity, in this section we will consider the
density contrast for the tracers,
$\delta^{i}=(n^{i}-\bar{n}^{i})/\bar{n}^{i}$, as our data set. The Fourier
mode $\vec{k}$ of the density contrast can be expressed as:

$d^{i}_{a}(\vec{k})=\\{\tilde{\delta}^{i}(\vec{k})\,,\,\tilde{\delta}^{i*}(\vec{k})\\}\;,$
(2.5)

where $i=1,2,\ldots,N_{t}$ denotes the tracer, and $a=1,2$ for the Fourier
mode and its complex conjugate, respectively. The data covariance is then:

$\langle
d^{i}_{a}(\vec{k})d^{j}_{b}(\vec{k}{}^{\prime})\rangle=D_{ab}\,\Gamma^{ij}(\vec{k},\vec{k}{}^{\prime})=\Gamma_{ab}^{ij}(\vec{k},\vec{k}{}^{\prime})\;,$
(2.6)

where $D_{ab}=1-\delta_{ab}$. The data covariance in Fourier space is simply
the observed power spectrum, including shot noise:

$\Gamma^{ij}(\vec{k},\vec{k}{}^{\prime})=\delta_{\vec{k}\,\vec{k}{}^{\prime}}\left(P^{ij}+\frac{\delta_{ij}}{\bar{n}^{i}}\right)\;,$
(2.7)

where $P^{ij}(\vec{k})$ is the power spectrum for the tracers $i,j$.

The Fisher matrix for the set of parameters $\theta^{\mu}$ is given by the
generalized trace:

$F_{\mu\nu}=\frac{1}{4}\sum_{k}V\tilde{V}_{k}\sum_{iji^{\prime}j^{\prime}}\sum_{aba^{\prime}b^{\prime}}\frac{\partial\,\Gamma_{ab}^{ij}}{\partial\theta^{\mu}}\left[\Gamma_{ba^{\prime}}^{ji^{\prime}}\right]^{-1}\frac{\partial\,\Gamma_{a^{\prime}b^{\prime}}^{i^{\prime}j^{\prime}}}{\partial\theta^{\nu}}\left[\Gamma_{b^{\prime}a}^{j^{\prime}i}\right]^{-1}\;,$
(2.8)

where $V$ is the survey volume, $\tilde{V}_{k}$ is the volume in Fourier space
of the bandpowers (Fourier bins) $k$, and the additional factor of $1/2$ in
Eq. (2.8) is due to the fact that our degrees of freedom in Eq. (2.5) count
the Fourier modes twice.

Now, using the fact that the data covariance is separable,
$\left[\Gamma_{ab}^{ij}\right]^{-1}=\left[D_{ab}\right]^{-1}\left[\Gamma^{ij}\right]^{-1}$,
and using that $\left[D_{ab}\right]^{-1}=D_{ab}$, we obtain that
$\sum_{abcd}D_{ab}D_{bc}D_{cd}D_{da}=\sum_{ac}\delta_{ac}\delta_{ca}=2$,
hence:

$F_{\mu\nu}=\frac{1}{2}\sum_{k}V\tilde{V}_{k}\sum_{iji^{\prime}j^{\prime}}\frac{\partial\,\Gamma^{ij}}{\partial\theta^{\mu}}\left[\Gamma^{ji^{\prime}}\right]^{-1}\frac{\partial\,\Gamma^{i^{\prime}j^{\prime}}}{\partial\theta^{\nu}}\left[\Gamma^{j^{\prime}i}\right]^{-1}\;,$
(2.9)

In Fourier space the inverse of the data covariance has a trivial expression,
namely:

$\left[\Gamma^{ij}\right]^{-1}=\bar{n}^{i}\,\delta_{ij}-\bar{n}^{i}\,\frac{P^{ij}}{1+{\cal{P}}}\,\bar{n}^{j}\;,$

with ${\cal{P}}=\sum_{i}\bar{n}^{i}\,P^{ii}$.

We now set the parameters $\theta^{\mu}$ to be the auto- and cross-spectra of
the tracers evaluated at some bandpower, $P^{ij}(k)$. But before we compute
that Fisher matrix, it is worth noting that these spectra (as well as the
corresponding data covariance) are symmetric, $P^{ij}=P^{ji}$. Therefore, if
we simply identify the parameters of the Fisher matrix as those spectra, we
would be counting the cross-spectra twice. We handle this double-counting by
defining non-degenerate spectra as:

$P^{[ij]}=\left\\{\begin{array}[]{cc}P^{ij}&i\leq j\\\
0&i>j\end{array}\right.\;.$ (2.10)

Conversely, we can rewrite the spectra in terms of the non-degenerate spectra
as:

$P^{ij}=P^{[ij]}+P^{[ji]}-\delta_{ij}P^{ii}\;,$ (2.11)

and a similar expression for the data covariance. We can now evaluate the
partial derivatives assuming that the parameters are the non-denegerate
spectra:

$\frac{\partial\Gamma^{ij}(k)}{\partial
P^{[i^{\prime}j^{\prime}]}(k^{\prime})}=\delta_{k,k^{\prime}}\left[\delta_{ii^{\prime}}\delta_{jj^{\prime}}+\delta_{ij^{\prime}}\delta_{ji^{\prime}}-\delta_{ij}\delta_{ji^{\prime}}\delta_{i^{\prime}j^{\prime}}\delta_{j^{\prime}i}\right]\,\equiv\,\delta_{k,k^{\prime}}\delta_{[ij],[i^{\prime}j^{\prime}]}\;.$
(2.12)

Substituting this identity into Eq. (2.9) results in a Fisher matrix which is
diagonal in the bandpowers, and can be expressed as [63]:

$F[P^{[ij]},P^{[i^{\prime}j^{\prime}]}]=F^{[ij],[i^{\prime}j^{\prime}]}=\frac{V\tilde{V}_{k}}{4}\,\left(2-\delta_{ij}\right)\left(2-\delta_{i^{\prime}j^{\prime}}\right)\,\left([\Gamma^{ii^{\prime}}]^{-1}[\Gamma^{jj^{\prime}}]^{-1}+[\Gamma^{ij^{\prime}}]^{-1}[\Gamma^{ji^{\prime}}]^{-1}\right)\,.$
(2.13)

Eq. (2.13) makes it clear that, in order to compute the Fisher matrix, we need
first to invert the data covariance $\Gamma^{ij}$. The reason we compute the
Fisher matrix for the power spectrum, but not for the correlation function, is
that we can, under some approximations, invert the Fourier space data
covariance [64, 60], whereas it is unfeasible to invert the configuration-
space data covariance. As we will show in the next Sections, it is possible to
invert the harmonic space data covariance as well, in real and in redshift
space. Ultimately, this is what has allowed us to compute the Fisher matrix
for the angular power spectrum.

Using this notation it is straightforward to show that the inverse of the
Fisher matrix of Eq. (2.13) is indeed the familiar expression for the
covariance of the spectra which follows from the 4-point function under the
assumption of Gaussianity:

${\rm Cov}[P^{[ij]},P^{[i^{\prime}j^{\prime}]}]={\rm
Cov}^{[ij],[i^{\prime}j^{\prime}]}=\frac{1}{V\tilde{V}_{k}}\,\left(\Gamma^{ii^{\prime}}\Gamma^{jj^{\prime}}+\Gamma^{ij^{\prime}}\Gamma^{ji^{\prime}}\right)\,.$
(2.14)

In order to show that Eq. (2.14) is in fact the inverse of Eq. (2.13), the
following identity is useful:

$\displaystyle
2\sum_{[mn]}[\Gamma^{im}]^{-1}\Gamma^{mi^{\prime}}[\Gamma^{jn}]^{-1}\Gamma^{nj^{\prime}}$
$\displaystyle=$
$\displaystyle\sum_{mn}[\Gamma^{im}]^{-1}\Gamma^{mi^{\prime}}[\Gamma^{jn}]^{-1}\Gamma^{nj^{\prime}}+\sum_{m}[\Gamma^{im}]^{-1}\Gamma^{mi^{\prime}}[\Gamma^{jm}]^{-1}\Gamma^{mj^{\prime}}$
$\displaystyle=$
$\displaystyle\delta_{ii^{\prime}}\delta_{jj^{\prime}}+\sum_{m}[\Gamma^{im}]^{-1}\Gamma^{mi^{\prime}}[\Gamma^{jm}]^{-1}\Gamma^{mj^{\prime}}\;.$

With the help of this expression it is then straightforward to show that:

$\sum_{[mn]}F^{[ij],[mn]}\,{\rm
Cov}^{[mn],[i^{\prime}j^{\prime}]}=\delta_{[ij],[i^{\prime}j^{\prime}]}\,,$
(2.15)

where we used the same notation as in Eq. (2.12).

###  2.2 Correlations in harmonic space: the angular power spectrum

The problem with the Fourier expressions found in the last Section is that,
first, we do not observe the Universe in snapshots ($t=constant$
hypersurfaces), but over the past light-cone. And second, the surfaces of
constant time in that light-cone are not flat, but 2D spherical shells. This
means that, in order not to miss any physics that depends on evolution,
redshift-space effects, or that is only manifested on the largest observable
scales, we need a description in terms of the distance over the light cone
(or, equivalently, redshift) and spherical harmonics. This representation can
be used both for the density fields (which are scalars under rotations) [65,
23, 45] as well as for higher-spin fields such as cosmic shear or flexions
[66].

We wish to express the number of galaxies of type $i$ that occupy some
spherical shell, with radial distances inside the interval
$x\in[\bar{x}_{i}-\Delta x/2,\bar{x}_{i}+\Delta x/2]$, in terms of an
expansion over spherical harmonic functions. We then define:

$N^{i}_{\ell m}(\bar{x}_{i})=\int_{\bar{x}_{i}}d^{3}x\,Y_{\ell
m}^{*}(\hat{x})\,n^{i}(\vec{x})=\int_{\bar{x}_{i}-\Delta
x/2}^{\bar{x}_{i}+\Delta x/2}dx\,x^{2}\,\int d^{2}\hat{x}\,Y_{\ell
m}^{*}(\hat{x})\,n^{i}(\vec{x})\;,$ (2.16)

where $\hat{x}=\vec{x}/x$. It is useful to define the angular density of
tracers in the radial bin $\bar{x}_{i}$, or counts in shells, as:

$\displaystyle N^{i}_{\bar{x}_{i}}(\hat{x})$ $\displaystyle=$
$\displaystyle\int_{\bar{x}_{i}}dx\,x^{2}\,n^{i}(x\,\hat{x})=\int_{\bar{x}_{i}}dx\,x^{2}\,\sum_{p=1}^{N^{i}_{Tot}}\delta_{D}(\vec{x}{\,}_{i}^{p}-x\,\hat{x})$
$\displaystyle=$
$\displaystyle\int_{\bar{x}_{i}}dx\,x^{2}\,\sum_{p=1}^{N^{i}_{Tot}}\,\frac{\delta_{D}(x^{p}_{i}-x)}{x^{2}}\,\delta_{D}(\hat{x}^{p}_{i}-\hat{x})=\sum_{q=1}^{N^{i}_{\bar{x}_{i},Tot}}\delta_{D}(\hat{x}^{q}_{i}-\hat{x})\;,$

where $N^{i}_{Tot}$ denotes the total number of tracers of type $i$, each one
at a position $\vec{x}{\,}_{i}^{p}$, and $N^{i}_{\bar{x}_{i},Tot}$ denotes the
number of those tracers that fall inside the radial bin $\bar{x}_{i}$. Notice
that we have expressed the Dirac delta function for the angular coordinates as
$\delta_{D}(\hat{x}-\hat{x}{}^{\prime})=\delta_{D}(\cos\theta-\cos\theta^{\prime})\,\delta_{D}(\varphi-\varphi^{\prime})$,
and that the total number of tracers in the radial bin can also be found by
integrating the angular density of tracers over the solid angle,
$N^{i}_{\bar{x}_{i},Tot}=\int d^{2}\hat{x}\,N^{i}_{\bar{x}_{i}}$. With these
definitions we can rewrite Eq. (2.16) as:

$N^{i}_{\ell m}(\bar{x}_{i})=\int d^{2}\hat{x}\,Y_{\ell
m}^{*}(\hat{x})\,N^{i}_{\bar{x}}(\hat{x})=\sum_{q=1}^{N^{i}_{\bar{x}_{i},Tot}}Y_{\ell
m}^{*}(\hat{x}^{q}_{i})\;,$ (2.18)

although we will not make use of the last expression above. Notice that, with
these definitions, and assuming a full sky, the monopole of the number of
galaxies in the bin $\bar{x}$ is
$N^{i}_{00}(\bar{x}_{i})=N^{i}_{\bar{x}_{i},Tot}/\sqrt{4\pi}$.

Let’s now define the harmonic space data covariance as the expectation value
of the differential counts of tracers in the shells of radii $\bar{x}$ and
$\bar{y}$:

$\displaystyle\left\langle\delta N^{i}_{\ell m}(\bar{x}_{i})\,\delta
N^{j*}_{\ell^{\prime}m^{\prime}}(\bar{x}_{j})\right\rangle$ $\displaystyle=$
$\displaystyle\Gamma^{ij}_{\ell}(\bar{x}_{i},\bar{x}_{j})\,\delta_{\ell\,\ell^{\prime}}\,\delta_{m\,m^{\prime}}$
$\displaystyle=$
$\displaystyle\int_{\bar{x}_{i}}dx_{i}\,x_{i}^{2}\int_{\bar{x}_{j}}dx_{j}\,x_{j}^{2}\int
d^{2}\hat{x}_{i}\int d^{2}\hat{x}_{j}\,Y_{\ell
m}^{*}(\hat{x}_{i})Y_{\ell^{\prime}m^{\prime}}(\hat{x}_{j})\,\Gamma^{ij}(\vec{x}_{i},\vec{x}_{j})\;,$

where symmetry under rotations leads to the diagonal structure of the angular
power spectrum111The couplings between different $(\ell,m)$ generated by masks
or non-trivial angular selection functions [23, 28] will be the subject of a
future paper. We also assume in this paper that the radial selection function
is constant inside each radial bin, and that the radial bins do not overlap..
Notice that the harmonic data covariance above reduces to the angular power
spectrum only after subtracting shot noise,
$\Gamma_{\ell}^{ij}(\bar{x}_{i},\bar{x}_{j})\to
C_{\ell}^{ij}(\bar{x}_{i},\bar{x}_{j})$.

The expression above relates the data covariance in configuration space,
$\Gamma^{ij}(\vec{x}_{i},\vec{x}_{j})$, to the harmonic data covariance,
$\Gamma_{\ell}^{ij}(\bar{x}_{i},\bar{x}_{j})$, and it can be inverted in the
sense:

$\displaystyle\Gamma^{ij}(\hat{x}_{i}\,\bar{x}_{i},\hat{x}_{j}\,\bar{x}_{j})$
$\displaystyle=$ $\displaystyle\frac{4\pi}{\Delta
V_{\bar{x}_{i}}}\frac{4\pi}{\Delta
V_{\bar{x}_{j}}}\int_{\bar{x}_{i}}dx_{i}\,x_{i}^{2}\int_{\bar{x}_{j}}dx_{j}\,x_{j}^{2}\,\Gamma^{ij}(\vec{x}_{i},\vec{x}_{j})$
$\displaystyle=$ $\displaystyle\frac{4\pi}{\Delta
V_{\bar{x}_{i}}}\frac{4\pi}{\Delta V_{\bar{x}_{j}}}\sum_{\ell m}Y_{\ell
m}(\hat{x}_{i})\,Y_{\ell
m}^{*}(\hat{x}_{j})\,\Gamma^{ij}_{\ell}(\bar{x}_{i},\bar{x}_{j})$
$\displaystyle=$ $\displaystyle\frac{4\pi}{\Delta
V_{\bar{x}_{i}}}\frac{4\pi}{\Delta
V_{\bar{x}_{j}}}\sum_{\ell}\frac{2\ell+1}{4\pi}{\cal{L}}_{\ell}(\hat{x}_{i}\cdot\hat{x}_{j})\,\Gamma^{ij}_{\ell}(\bar{x}_{i},\bar{x}_{j})\;,$

where the volume of the spherical shell around $\bar{x}$ is given by $\Delta
V_{\bar{x}}=\int_{\bar{x}}d^{3}x=4\pi\,\bar{x}^{2}\,\Delta x$. The last line
serves to make it explicit the fact that the correlation function in
configuration space depends on the angles only through the combination
$\mu_{ij}=\hat{x}_{i}\cdot\hat{x}_{j}$ that appear in the argument of the
Legendre polynomial ${\cal{L}}_{\ell}$. Notice also that, in contrast to the
configuration-space covariance, which was defined in terms of the number
densities $n^{i}$, the harmonic-space covariance of Eq. (2.2) is defined in
terms of the counts in shells. As a result, whereas the configuration space
data covariance, $\Gamma^{ij}(\vec{x}_{i},\vec{x}_{j})$, has dimensions of
1/(volume)2, the harmonic space data covariance in spherical shells,
$\Gamma^{ij}_{\ell}(\bar{x}_{i},\bar{x}_{j})$, is adimensional.

Let’s say that we are able to find an expression for the inverse of the data
covariance in harmonic space, such that:

$\sum_{j}\sum_{\bar{y}}\left[\Gamma^{ij}_{\ell}(\bar{x},\bar{x}_{j})\right]^{-1}\Gamma^{ji^{\prime}}_{\ell}(\bar{x}_{j},\bar{x}_{i^{\prime}})=\delta_{ii^{\prime}}\delta_{\bar{x}_{i}\,\bar{x}_{i^{\prime}}}\;.$
(2.21)

Given such an expression, we could define the inverse of the data covariance
in configuration space in analogy to Eq. (2.2), as:

$\left[\Gamma^{ij}(\hat{x}_{i}\,\bar{x}_{i},\hat{x}_{j}\,\bar{x}_{j})\right]^{-1}=\frac{\Delta
V_{\bar{x}_{i}}}{4\pi}\frac{\Delta V_{\bar{x}_{j}}}{4\pi}\sum_{\ell m}Y_{\ell
m}(\hat{x}_{i})\,Y_{\ell
m}^{*}(\hat{x}_{j})\left[\Gamma^{ij}_{\ell}(\bar{x}_{i},\bar{x}_{j})\right]^{-1}\;.$
(2.22)

One can then take the limit back to the continuum by making the volume
elements infinitesimally small. It is then trivial to show, with the help of
Eqs. (2.2)-(2.21), that the expression above is indeed the inverse of the
covariance:

$\sum_{j}\int
d^{3}x_{j}\left[\Gamma^{ij}(\vec{x}_{i},\vec{x}_{j})\right]^{-1}\Gamma^{ji^{\prime}}(\vec{x}_{j},\vec{x}_{i^{\prime}})=\delta_{ii^{\prime}}\delta_{D}(\vec{x}_{i}-\vec{x}_{i^{\prime}})\;,$
(2.23)

or, equivalently, in finite spherical shells:

$\sum_{j}\sum_{\bar{x}_{j}}\Delta
V_{\bar{x}_{j}}\left[\Gamma^{ij}(\hat{x}_{i}\,\bar{x}_{i},\hat{x}_{j}\,\bar{x}_{j})\right]^{-1}\Gamma^{ji^{\prime}}(\hat{x}_{j}\,\bar{x}_{j},\hat{x}_{i^{\prime}}\,\bar{x}_{i^{\prime}})=\delta_{ii^{\prime}}\frac{\delta_{\bar{x}_{i}\,\bar{x}_{i^{\prime}}}}{\Delta
V_{\bar{x}_{i}}}\,\delta_{D}(\hat{x}_{i}-\hat{x}_{i^{\prime}})\;,$ (2.24)

where we made use of the identity:

$\sum_{\ell m}Y_{\ell m}(\hat{x})Y_{\ell
m}^{*}(\hat{y})=\delta_{D}(\hat{x}-\hat{y})\;.$

We can now go back to the expression for the Fisher matrix given in Eq. (2),
and substitute the covariances in configuration space as well as their
inverses, Eqs. (2.2) and (2.22). After a bit of algebra we obtain:

$\displaystyle F^{\mu\nu}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{\ell}(2\ell+1)\sum_{i,\bar{x}_{i}}\sum_{j,\bar{x}_{j}}\sum_{i^{\prime},\bar{x}_{i^{\prime}}}\sum_{j,\bar{x}_{j^{\prime}}}\frac{\Delta
V_{\bar{x}_{i}}}{4\pi}\frac{\Delta V_{\bar{x}_{j}}}{4\pi}\frac{\Delta
V_{\bar{x}_{i^{\prime}}}}{4\pi}\frac{\Delta V_{\bar{x}_{j^{\prime}}}}{4\pi}$
$\displaystyle\times\frac{\partial\Gamma^{ij}_{\ell}(\bar{x}_{i},\bar{x}_{j})}{\partial\theta^{\mu}}\left[\Gamma^{ji^{\prime}}_{\ell}(\bar{x}_{j},\bar{x}_{i^{\prime}})\right]^{-1}\frac{\partial\Gamma^{i^{\prime}j^{\prime}}_{\ell}(\bar{x}_{i^{\prime}},\bar{x}_{j^{\prime}})}{\partial\theta^{\nu}}\left[\Gamma^{j^{\prime}i}_{\ell}(\bar{x}_{j^{\prime}},\bar{x}_{i})\right]^{-1}\;.$

In deriving the expression above, the angular integrals cancel all the
dependence on the spherical harmonic functions, such that only one sum over
$(\ell,m)$ remains, which leads to the factor $\sum_{m=-\ell}^{\ell}=2\ell+1$
since none of the terms in the final expression depends on $m$.

What the result above shows is that the Fisher matrix can be decomposed in
linearly independent spherical harmonic components:

$F^{\mu\nu}=\sum_{\ell}F^{\mu\nu}_{\ell}\;,$ (2.26)

where we define:

$F_{\ell}^{\mu\nu}=\frac{2\ell+1}{2}\,{\rm
Tr}\left\\{\frac{\partial\Gamma^{ij}_{\ell}}{\partial\theta^{\mu}}\left[\Gamma^{ji^{\prime}}_{\ell}\right]^{-1}\frac{\partial\Gamma^{i^{\prime}j^{\prime}}_{\ell}}{\partial\theta^{\nu}}\left[\Gamma^{j^{\prime}i}_{\ell}\right]^{-1}\right\\}\;,$
(2.27)

with the trace denoting, in short-hand notation, both sums over tracer indices
as well as integrals over the corresponding radii, $\sum_{\bar{x}_{i}}\Delta
V_{\bar{x}_{i}}/4\pi(\cdots)\to\int dx_{i}\,x_{i}^{2}(\cdots)$. Since we will
choose our parameters $\theta^{\mu}\to C^{ij}_{\ell}$, and $\partial
C^{ij}_{\ell}/\partial
C^{i^{\prime}j^{\prime}}_{\ell^{\prime}}\sim\delta_{\ell\ell^{\prime}}$, the
Fisher matrix for the angular power spectrum is given by Eq. (2.27). With
partial sky coverage $f_{sky}=\Delta\Omega/4\pi$, where $\Delta\Omega$ is the
survey angular area, and collecting different $\ell$’s inside a bin
$\bar{\ell}$, the expression becomes [67]:

$F_{\bar{\ell}}^{\mu\nu}=\frac{f_{sky}}{2}\,\sum_{\ell\in\bar{\ell}}\,(2\ell+1)\,{\rm
Tr}\left\\{\frac{\partial\Gamma^{ij}_{\ell}}{\partial\theta^{\mu}}\left[\Gamma^{ji^{\prime}}_{\ell}\right]^{-1}\frac{\partial\Gamma^{i^{\prime}j^{\prime}}_{\ell}}{\partial\theta^{\nu}}\left[\Gamma^{j^{\prime}i}_{\ell}\right]^{-1}\right\\}\;.$
(2.28)

###  2.3 Fisher matrix in harmonic space

Let’s now say that we are able to invert the harmonic data covariance (i.e.,
the observed angular power spectrum), in real or in redshift space. We can
then return to Eq. (2.27) and start writing the Fisher matrix for the
observables of interest, which in our case are the angular power spectra,
$\theta^{\mu}\to C^{ij}_{\ell}(\bar{x},\bar{y})$.

As was the case for the power spectrum in Fourier space, we must be aware of
the potential for double counting the degrees of freedom. For that reason, in
analogy with Eq. (2.10) we define:

$C^{[ij]}_{\ell}(\bar{x}_{i},\bar{x}_{j})=\left\\{\begin{array}[]{cc}C^{ij}_{\ell}(\bar{x}_{i},\bar{x}_{j})&i\leq
j\\\ 0&i>j\end{array}\right.\;,$ (2.29)

and

$C_{\ell}^{ij}(\bar{x}_{i},\bar{x}_{j})=C_{\ell}^{[ij]}(\bar{x}_{i},\bar{x}_{j})+C_{\ell}^{[ji]}(\bar{x}_{j},\bar{x}_{i})-\delta_{ij}C^{ii}(\bar{x}_{i},\bar{x}_{j})\;,$
(2.30)

where the auto-correlation in the last term is symmetric in the spatial
indices by definition,
$C_{\ell}^{ii}(\bar{x}_{i},\bar{x}_{j})=C_{\ell}^{ii}(\bar{x}_{j},\bar{x}_{i})$.

The partial derivatives in the Fisher matrix in the present case are given by
the analog of Eq. (2.12):

$\frac{\partial C_{\ell}^{ij}(\bar{x}_{i},\bar{x}_{j})}{\partial
C_{\ell}^{[i^{\prime}j^{\prime}]}(\bar{x}_{i^{\prime}},\bar{x}_{j^{\prime}})}=\delta_{ii^{\prime}}\delta_{jj^{\prime}}\,\delta_{\bar{x}_{i}\bar{x}_{i^{\prime}}}\delta_{\bar{x}_{j}\bar{x}_{j^{\prime}}}+\delta_{ij^{\prime}}\delta_{ji^{\prime}}\,\delta_{\bar{x}_{i}\bar{x}_{j^{\prime}}}\delta_{\bar{x}_{j}\bar{x}_{i^{\prime}}}-\delta_{ij}\delta_{ji^{\prime}}\delta_{i^{\prime}j^{\prime}}\delta_{j^{\prime}i}\,\delta_{\bar{x}_{i}\bar{x}_{j}}\delta_{\bar{x}_{j}\bar{x}_{i^{\prime}}}\delta_{\bar{x}_{i^{\prime}}\bar{x}_{j^{\prime}}}\delta_{\bar{x}_{j^{\prime}}\bar{x}_{i}}\;.$

At this point it should have become abundantly clear that the tracer index
always appears tied up with the radial bin: $\\{i,\bar{x}_{i}\\}$,
$\\{j,\bar{x}_{j}\\}$, etc. Hence, we will often use this association in order
to shorten our notation, by implicitly assuming that the two are connected. In
particular, we can define
$\bar{\delta}_{ii^{\prime}}=\delta_{ii^{\prime}}\delta_{\bar{x}_{i}\bar{x}_{i^{\prime}}}$,
leading to:

$\frac{\partial C_{\ell}^{ij}}{\partial
C_{\ell}^{[i^{\prime}j^{\prime}]}}=\bar{\delta}_{ii^{\prime}}\bar{\delta}_{jj^{\prime}}+\bar{\delta}_{ij^{\prime}}\bar{\delta}_{ji^{\prime}}-\bar{\delta}_{ij}\bar{\delta}_{ji^{\prime}}\bar{\delta}_{i^{\prime}j^{\prime}}\bar{\delta}_{j^{\prime}i}\;.$
(2.31)

Substituting this expression into Eq. (2.27) leads us to:

$\displaystyle F[C^{[ij]}_{\ell},C^{[i^{\prime}j^{\prime}]}_{\ell}]$
$\displaystyle=$
$\displaystyle\frac{2\ell+1}{4}(2-\bar{\delta}_{ij})(2-\bar{\delta}_{i^{\prime}j^{\prime}})\left\\{[{\Gamma}^{ii^{\prime}}_{\ell}]^{-1}[{\Gamma}^{jj^{\prime}}_{\ell}]^{-1}+[{\Gamma}^{ij^{\prime}}_{\ell}]^{-1}[{\Gamma}^{i^{\prime}j}_{\ell}]^{-1}\right\\}\;,$
(2.32)

in full analogy with the result in Fourier space. And again, in exact
correspondence with Fourier space, the covariance matrix in harmonic space is
given by:

$\displaystyle{\rm Cov}[C^{[ij]}_{\ell},C^{[i^{\prime}j^{\prime}]}_{\ell}]$
$\displaystyle=$
$\displaystyle\frac{1}{2\ell+1}\left[{\Gamma}^{ii^{\prime}}_{\ell}{\Gamma}^{jj^{\prime}}_{\ell}+{\Gamma}^{ij^{\prime}}_{\ell}{\Gamma}^{i^{\prime}j}_{\ell}\right]\;.$
(2.33)

##  3 The harmonic Fisher matrix in real space

In this Section we will derive exact expressions for the Fisher and covariance
matrices of the angular power spectrum for two redshift slices,
$C_{\ell}^{ij}(\bar{x}_{i},\bar{x}_{j})$, under the assumption of linear
biasing for the tracers, but without any redshift-space distortions. In the
Section 4 we will generalize this result to redshift space, but still under
the assumption of linear biasing and using the Kaiser approximation [61] –
i.e., assuming that redshift-space distortions are in the linear regime of
structure formation.

###  3.1 The data covariance in harmonic space

Let’s go back to the covariance in harmonic space, Eq. (2.2), and write it in
the form:

$\displaystyle\Gamma^{ij}_{\ell}$ $\displaystyle=$
$\displaystyle\int_{\bar{x}_{i}}dx_{i}\,x_{i}^{2}\int_{\bar{x}_{j}}dx_{j}\,x_{j}^{2}\int
d^{2}\hat{x}_{i}\int d^{2}\hat{x}_{j}\,Y_{\ell m}(\hat{x}_{i})Y_{\ell
m}^{*}(\hat{x}_{j})\left\langle\delta n^{i}(\vec{x}_{i})\delta
n^{j}(\vec{x}_{j})\right\rangle$ $\displaystyle=$
$\displaystyle\int_{\bar{x}_{i}}dx_{i}\,x_{i}^{2}\int_{\bar{x}_{j}}dx_{j}\,x_{j}^{2}\int
d^{2}\hat{x}_{i}\int d^{2}\hat{x}_{j}\,Y_{\ell m}(\hat{x}_{i})Y_{\ell
m}^{*}(\hat{x}_{j})\left[\bar{n}^{i}(\vec{x}_{i})\,\bar{n}^{j}(\vec{x}_{j})\,\xi^{ij}(\vec{x}_{i},\vec{x}_{j})\right.$
$\displaystyle+\left.\delta_{ij}\,\bar{n}^{i}(\vec{x}_{i})\,\delta_{D}(\vec{x}_{i}-\vec{x}_{j})\right]\;,$

The first term in this expression is the correlation function of the counts in
shells, in harmonic space – i.e., the angular power spectrum. If the number
densities of tracers in the shells vary as a function of the angular position
due to masks or selection functions, then the identity above becomes an
angular convolution between the correlation function and those masks and
selection functions. For simplicity, here we will assume that the mean numbers
of tracers is nearly constant, so that it can be pulled outside the integral
over the radial bins:

$C^{ij}_{\ell}(\bar{x}_{i},\bar{x}_{j})=\bar{n}^{i}(\bar{x}_{i})\bar{n}^{j}(\bar{x}_{j})\int_{\bar{x}_{i}}dx_{i}\,x_{i}^{2}\int_{\bar{x}_{j}}dx_{j}\,x_{j}^{2}\int
d^{2}\hat{x}_{i}\,Y_{\ell m}(\hat{x}_{i})\int d^{2}\hat{x}_{j}\,Y^{*}_{\ell
m}(\hat{x}_{j})\,\xi^{ij}(\vec{x}_{i},\vec{x}_{j})\,.$ (3.2)

Notice also that under the assumption of a trivial angular selection function,
the shot noise term in Eq. (3.1) reduces to
$\delta_{ij}\,\bar{N}^{i}_{\bar{x}_{i}}=\delta_{ij}\,\bar{n}^{i}_{\bar{x}_{i}}\,\Delta
V_{\bar{x}_{i}}$.

Now comes a simple but crucial step that is often overlooked: if we assume
that the correlation function is computed at a $t=\,$constant hypersurface,
and we assume homogeneity, then we can write that, for all
$|\vec{x}_{i}|\in\bar{x}_{i}$ and all $|\vec{x}_{j}|\in\bar{x}_{j}$ we have:

$\displaystyle\xi^{ij}(\vec{x}_{i},\vec{x}_{j})$ $\displaystyle=$
$\displaystyle\int\frac{d^{3}k}{(2\pi)^{3}}\,e^{i\vec{k}(\vec{x}_{i}-\vec{x}_{j})}P^{ij}(\vec{k})\,.$
(3.3)

Here the reader should bear in mind that the power spectrum is given in terms
of the expression
$\langle\tilde{\delta}^{i}(\vec{k}_{i}|z_{i})\tilde{\delta}^{j*}(\vec{k}_{j}|z_{j})\rangle=(2\pi)^{3}\delta_{D}(\vec{k}_{i}-\vec{k}_{j})P^{ij}(\vec{k}_{i})$,
and therefore the power spectrum includes factors of the tracer bias as well
as growth function. In particular, in real space and in the linear regime we
would have $P^{ij}(\vec{k})\to
b^{i}(z_{i})\,b^{j}(z_{j})\,D(z_{i})\,D(z_{j})\,P_{m}(k|z=0)$.

At this point it is important to remember that the wavenumber $\vec{k}$ in the
integral above has support in the entire 3D volume, and not only in the
spherical shells $\bar{x}_{i}$ and $\bar{x}_{j}$. Now, recall that the two
spherical shells $z_{i}$ (of comoving radius $\bar{x}_{i}$) and $z_{j}$ (of
radius $\bar{x}_{j}$) are not in the same $t=\,$constant hypersurfaces – in
fact, the corresponding observables are on the past light-cone of the
observer. Eq. (3.3) is, therefore, the result of bringing both spherical
shells to the same configuration space volume (and the same $t=\,$constant
hypersurface), and then taking the expectation value of the 2-point function
of the Fourier modes
$\langle\tilde{\delta}^{i}(\vec{k}_{i})\tilde{\delta}^{j*}(\vec{k}_{j})\rangle$
in that hypersurface – which is arbitrary, although it is often assumed to be
the one at $z=0$. The procedure is pictured in figure 1. If we assume that the
number densities and biases corresponding to the annuli $\bar{x}_{i}$ and
$\bar{x}_{j}$ are kept fixed, and we restrict our attention to the linear
regime, then it makes no difference what is the common $t=\,$constant
hypersurface that we chose to compute the expectation values. If, on the other
hand, we allow for non-linearities that change as a function of time due to
the growth of structures, then in principle the $t=\,$constant hypersurface
should be defined carefully – or, in a more general approach, one could
consider using unequal-time correlation functions [68].

Figure 1: The fields corresponding to the tracer $i$ on the shell
$\bar{x}_{i}$ and to the tracer $j$ on the shell $\bar{x}_{j}$ lie on the past
light-cone of the observer (shown on the left). In order to compute
correlation functions (in real or in Fourier space) these fields must be first
brought to the same $t=\,$constant hypersurface (shown on the right).
Homogeneity then allows us to extend the modes inside the annuli $\bar{x}_{i}$
and $\bar{x}_{j}$ to the entire hypersurface, and then to compute expectation
values such as the correlation function $\xi^{ij}$ and the power spectrum
$P^{ij}$. Notice that the quantities $\bar{n}^{i}$, $b^{i}$, $\bar{n}^{j}$ and
$b^{j}$, etc., are kept fixed in the common $t=\,$constant hypersurface.

Coming back to Eq. (3.2), substituting Eq. (3.3) and using the Rayleigh
expansion of a plane wave in spherical harmonics:

$e^{i\vec{k}\cdot\vec{x}}=4\pi\sum_{\ell m}i^{\ell}j_{\ell}(kx)Y_{\ell
m}(\hat{k})Y_{\ell m}^{*}(\hat{x})\;,$ (3.4)

after a bit of trivial algebra we obtain:

$\displaystyle C^{ij}_{\ell}$ $\displaystyle=$
$\displaystyle\bar{n}^{i}\,\bar{n}^{j}\int_{\bar{x}_{i}}dx_{i}\,x_{i}^{2}\int_{\bar{x}_{j}}dx_{j}\,x_{j}^{2}\times\frac{2}{\pi}\int_{0}^{\infty}dk\,k^{2}\,j_{\ell}(kx_{i})\,j_{\ell}(kx_{j})\,P^{ij}(k)\,.$
(3.5)

In the linear regime we can simplify this expression further, using
$P^{ij}(k)=b^{i}b^{j}D^{i}D^{j}P^{(m)}(k)$, where $P^{(m)}(k)$ is the matter
power spectrum and $D^{i}=D(z_{i})$ is the growth function. Taking the biases
of the tracers outside the volume integral (again a safe assumption, if the
slices are thin enough) we obtain:

$\displaystyle C^{ij}_{\ell}$ $\displaystyle=$
$\displaystyle\bar{n}^{i}\bar{n}^{j}\,b^{i}\,b^{j}\,D^{i}\,D^{j}\int_{\bar{x}_{i}}dx_{i}\,x_{i}^{2}\int_{\bar{x}_{j}}dx_{j}\,x_{j}^{2}\times\frac{2}{\pi}\int_{0}^{\infty}dk\,k^{2}\,j_{\ell}(kx_{i})\,j_{\ell}(kx_{j})\,P^{(m)}(k)$
(3.6) $\displaystyle=$
$\displaystyle\bar{n}^{i}\bar{n}^{j}\,b^{i}\,b^{j}\,C^{(m)}_{\ell}(\bar{x}_{i},\bar{x}_{j})\,,$

where $C^{(m)}_{\ell}(\bar{x}_{i},\bar{x}_{j})$ is the matter angular power
spectrum for the two slices. It is also possible to extend this expression to
take into account non-linearities – see, e.g., [69]

Another useful approximation arises if the Bessel functions are nearly
constant inside each radial bin:

$\displaystyle C^{ij}_{\ell}(\bar{x}_{i},\bar{x}_{j})$ $\displaystyle\approx$
$\displaystyle\bar{N}^{i}\bar{N}^{j}\times\frac{2}{\pi}\int_{0}^{\infty}dk\,k^{2}\,j_{\ell}(k\bar{x}_{i})\,j_{\ell}(k\bar{x}_{j})\,P^{ij}(k)\equiv\bar{C}{}^{ij}_{\ell}(\bar{x}_{i},\bar{x}_{j})\,,$
(3.7)

where we defined $\bar{N}^{i}=\bar{n}^{i}(\bar{x}_{i})\,\Delta
V_{\bar{x}_{i}}/4\pi$ as the mean number of tracers per unit solid angle, and
we assumed that the spherical shells are sufficiently thin, such that we can
approximate the argument of the spherical Bessel functions to the mean radius,
$x\to\bar{x}$. We will not rely on the approximation of Eq. (3.7) for now –
although it will turn out to be useful later. For more details on this
approximation we refer the reader to Appendix A.

At this point it is relevant to recall that we are working under the
assumption of linear regime. At $z=0$, scales smaller than about $\Delta
r\lesssim 30\,h^{-1}$ Mpc (i.e., $k\gtrsim 0.2\,h$ Mpc-1) are already affected
by non-linear effects, but at higher redshifts these non-linear scales move in
proportion to the matter growth rate. Therefore, one can think of the minimum
radial bin widths where our results would still remain applicable as $\Delta
r\gtrsim 30\,h^{-1}$ Mpc$/(1+z)$, or $\Delta z\gtrsim 0.01/(1+z)$.

###  3.2 Inverting the data covariance matrix

We now return to the task of inverting the data covariance of Eq. (3.1), in
light of the expression for the correlation function found in Eq. (3.5). Using
the closure relation for spherical Bessel functions,

$\int_{0}^{\infty}dt\,t^{2}\,j_{\ell}(at)\,j_{\ell}(bt)=\frac{\pi}{2}\frac{\delta_{D}(a-b)}{a^{2}}\;,$
(3.8)

allows us to rewrite the harmonic data covariance as:

$\displaystyle\Gamma^{ij}_{\ell}$ $\displaystyle=$
$\displaystyle\bar{n}^{i}\bar{n}^{j}\int_{\bar{x}_{i}}dx_{i}\,x_{i}^{2}\int_{\bar{x}_{j}}dx_{j}\,x_{j}^{2}\times\frac{2}{\pi}\int_{0}^{\infty}dk\,k^{2}\,j_{\ell}(kx_{i})\,j_{\ell}(kx_{j})\,M^{ij}(k)\;,$
(3.9)

where we defined the matrix:

$M^{ij}(k)=P^{ij}(k)+\frac{\delta_{ij}}{\sqrt{\bar{n}^{i}\,\bar{n}^{j}}}\;.$
(3.10)

This last expression indicates how we can go about inverting the harmonic
covariance. But before we present that inverse covariance, consider that, in
the continuum limit, we have:

$\sum_{\bar{x}}\int_{\bar{x}}dx\,f(x)\,g(\bar{x})\to\int dx\,f(x)\,g(x)\;.$

With the help of this limit, as well as Eq. (3.9), it is easy to show that the
inverse data covariance is given by:

$\left[\Gamma^{ji^{\prime}}_{\ell}\right]^{-1}=\frac{1}{\bar{n}^{j}\,\bar{n}^{i^{\prime}}}\times\frac{2}{\pi}\int_{0}^{\infty}dk\,k^{2}\,j_{\ell}(k\bar{x}_{j})\,j_{\ell}(k\bar{x}_{i^{\prime}})\,\left[M^{ji^{\prime}}(k)\right]^{-1}\;,$
(3.11)

which then leads to Eq. (2.21), i.e.:

$\sum_{j}\sum_{\bar{x}_{j}}\,\Gamma^{ij}_{\ell}(\bar{x}_{i},\bar{x}_{j})\,\left[\Gamma^{ji^{\prime}}_{\ell}(\bar{x}_{j},\bar{x}_{i^{\prime}})\right]^{-1}=\delta_{ii^{\prime}}\delta_{\bar{x}_{i}\,\bar{x}_{i^{\prime}}}\;.$
(3.12)

It is important to recall that, according to the discussion in the previous
section, the integration variable above, $\bar{x}_{j}$, corresponds to the
radius in a $t=\,$constant hypersurface, and the index $j$ corresponds to all
tracers in that hypersurface, which means that the tracer number densities
$\bar{n}^{j}$ and biases $b^{j}$ are kept constant – i.e., they do not depend
on $\bar{x}_{j}$.

In the definition of the inverse covariance, Eq. (3.11), the inverse of the
matrix $M^{ij}$ appears. We first rewrite that matrix as:

$\displaystyle M^{ij}(k)$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{\bar{n}^{i}}}\left[\delta_{ij}+{\cal{P}}^{ij}(k)\right]\frac{1}{\sqrt{\bar{n}^{j}}}\;,$
(3.13)

where in the last line we defined
${\cal{P}}^{ij}=\sqrt{\bar{n}^{i}\,\bar{n}^{j}}\,P^{ij}$. Notice also that
this adimensional quantity is separable, in the sense that
${\cal{P}}^{ij}{\cal{P}}^{i^{\prime}j^{\prime}}={\cal{P}}^{ii^{\prime}}{\cal{P}}^{jj^{\prime}}$,
and $\sum_{l}{\cal{P}}^{il}{\cal{P}}^{jl}={\cal{P}}^{ij}{\cal{P}}$, where we
defined the trace ${\cal{P}}=\sum_{i}{\cal{P}}^{ii}$. Due to these properties,
the inverse of the matrix $M$ is given simply by:

$\displaystyle[M^{ij}]^{-1}$ $\displaystyle=$
$\displaystyle\sqrt{\bar{n}^{i}}\left[\delta_{ij}-\frac{{\cal{P}}^{ij}}{1+{\cal{P}}}\right]\sqrt{\bar{n}^{j}}\;.$
(3.14)

At this point it is useful to understand the fundamental reason why we were
able to invert the harmonic data covariance in real space, and which will
later allow us to also invert that covariance including redshift-space
distortions. First, we rewrite the covariance as:

$\displaystyle\Gamma^{ij}_{\ell}(\bar{x}_{i},\bar{x}_{j})$ $\displaystyle=$
$\displaystyle\bar{N}^{i}\bar{N}^{j}\;\frac{2}{\pi}\int_{0}^{\infty}dk_{i}\,k_{i}^{2}\;\frac{2}{\pi}\int_{0}^{\infty}dk_{j}\,k_{j}^{2}\;j_{\ell}(k_{i}\bar{x}_{i})\,j_{\ell}(k_{j}\bar{x}_{j})\,\Gamma^{ij}_{\ell}(k_{i},k_{j})\,,$
(3.15)

where we defined the Fourier conjugate of the data covariance:

$\displaystyle\Gamma^{ij}_{\ell}(k_{i},k_{j})$ $\displaystyle=$
$\displaystyle\frac{\pi}{2}\frac{\delta{(k_{i}-k_{j})}}{k_{i}^{2}}\left[\frac{\delta_{ij}}{\bar{n}^{i}}+P^{ij}(k_{i})\right]\;.$
(3.16)

Since this object is perfectly diagonal in the Fourier modes, it is
straightforward to find its inverse – in stark contrast with the angular power
spectrum as a function of the radii of the spherical shells, which cannot be
directly inverted. Therefore, in Fourier space the problem reduces to finding
the inverse of a matrix in the tracer indices, $M^{ij}$, which is completely
trivial. This argument clarifies what we have done above: we first found the
inverse of Eq. (3.16), and then transformed it back to real space, arriving at
Eq. (3.11). We will use a very similar approach later, in Section 4, when we
tackle the angular power spectrum in redshift space.

###  3.3 An exact solution

Consider a single tracer with number density $\bar{n}$, unit bias ($b=1$), and
matter density fluctuations with a “top-hat” power spectrum
$P(k)=P_{0}\left[\theta_{H}(k-k_{1})-\theta_{H}(k-k_{2})\right]$, where
$\theta_{H}$ is the Heaviside (step) function — i.e., the spectrum is only
non-zero for $k_{1}\leq k\leq k_{2}$. To be clear, in this example we neglect
any time evolution effects as well as redshift-space distortions. The Fourier-
space data covariance is then by:

$\Gamma(k)=V_{0}[1+A\,\theta_{H}(k-k_{1})-A\,\theta_{H}(k-k_{2})]\;,$ (3.17)

where $V_{0}=1/\bar{n}$ corresponds to the shot noise term, and $A$ is an
adimensional parameter, so that the amplitude of the spectrum is
$P_{0}=AV_{0}$.

As shown in Appendix B, the harmonic data covariance corresponding to this
power spectrum is given by the analytical expressions:

$\displaystyle\Gamma_{\ell}(\bar{x},\bar{y})$ $\displaystyle=$
$\displaystyle\bar{N}_{\bar{x}}\bar{N}_{\bar{y}}V_{0}\left\\{\delta_{\bar{x}\,\bar{y}}\frac{4\pi}{\Delta
V_{\bar{x}}}+\,\frac{2}{\pi}A\left[\,k_{1}^{3}\,g_{\ell}(k_{1}\bar{x},k_{1}\bar{y})-k_{2}^{3}\,g_{\ell}(k_{2}\bar{x},k_{2}\bar{y})\right]\right\\}\;,$
(3.18) $\displaystyle\Gamma_{\ell}^{-1}(\bar{x},\bar{y})$ $\displaystyle=$
$\displaystyle V_{0}\left\\{\delta_{\bar{x}\,\bar{y}}\frac{4\pi}{\Delta
V_{\bar{x}}}-\,\frac{2}{\pi}\frac{A}{1+A}\,\left[k_{1}^{3}\,g_{\ell}(k_{1}\bar{x},k_{1}\bar{y})-k_{2}^{3}\,g_{\ell}(k_{2}\bar{x},k_{2}\bar{y})\right]\right\\}\;,$
(3.19)

where $\bar{N}_{\bar{x}}=\bar{n}\,\Delta V_{\bar{x}}/4\pi$ is the number of
tracers per unit solid angle in the radial bin ${\bar{x}}$, and:

$\displaystyle g_{\ell}(x,y)$ $\displaystyle\equiv$
$\displaystyle-\frac{x\,j_{\ell-1}(x)j_{\ell}(y)-y\,j_{\ell-1}(y)j_{\ell}(x)}{x^{2}-y^{2}}$
$\displaystyle=$
$\displaystyle\frac{j_{0}(x-y)-(-1)^{\ell}\,j_{0}(x+y)}{2\,x\,y}-\frac{1}{x\,y}\sum_{n=0}^{n_{\ell}}[2(\ell-1-2n)+1]\,j_{\ell-2n}(x)\,j_{\ell-2n}(y)\;.$

The sum in the second line above is over $n=0,1,\ldots,n_{\ell}=(2\ell-1)/4$,
so that $2(\ell-1-2n)+1\geq 0$, so that the sum stops at $n_{\ell}=(\ell-1)/2$
for odd values of $\ell$, and at $n_{\ell}=\ell/2$ for even values of $\ell$.
The expression of the second line shows explicitly that this function is well-
behaved as $x\to y$, and it is symmetric under $x\leftrightarrow y$, since
$j_{0}(x-y)=j_{0}(y-x)$.

Figure 2: Correlation matrix of the angular power spectrum corresponding to a
top-hat Fourier power spectrum, for the multipole $\ell=10$. We split the
radial coordinate in 25 bins of width $\Delta x=10\,h^{-1}\,{\rm Mpc}$ from
$\bar{x}_{1}=5\,h^{-1}\,{\rm Mpc}$ to $\bar{x}_{25}=245\,h^{-1}\,{\rm Mpc}$,
resulting in a $625\times 625$ covariance matrix. The left panel shows the
analytical (“Theory”) correlation matrix computed with the help of Eq. (2.33).
The right panel shows the sample correlation from 1000 Gaussian simulations.
The insets show a zoom in the same parts of the two matrices, for better
comparison between theory and simulations.

We can use the expressions in Eqs. (3.18)-(3.19) to compare our analytical
expression for the covariance with Guassian simulations. We produced a sample
of 1000 simulation boxes using the “top-hat” spectrum, and computed the sample
mean and sample variance of the angular power spectra. We use 25 redshift
slices corresponding to radial bins of width $\Delta x=10\,h^{-1}\,{\rm Mpc}$
from $\bar{x}=5\,h^{-1}\,{\rm Mpc}$ to $\bar{x}=245\,h^{-1}\,{\rm Mpc}$, and
for brevity we focus here on the multipole $\ell=10$. For more details, see
Appendix B.

In figure 2 we compare the analytical covariance of Eq. (2.33) (left panel)
with the sample covariance (right panel). For visualization purposes we show
the entire $[N_{r}\times N_{r}]\times[N_{r}\times N_{r}]$ matrix, where
$N_{r}$ is the number of radial bins from our mock simulations – i.e., we
compute the whole ${\rm
Cov}[C_{\ell}(\bar{x}_{i},\bar{x}_{j}),C_{\ell}(\bar{x}_{i^{\prime}},\bar{x}_{j^{\prime}})]$,
even if that carries redundancies due to double-counting the cross-spectra.
Moreover, in order to further facilitate visualization, we display the
correlation matrix, ${\rm Corr}[X,X^{\prime}]={\rm
Cov}[X,X^{\prime}]/\sqrt{{\rm Cov}[X,X]{\rm Cov}[X^{\prime},X^{\prime}]}$.
Since we have 25 radial bins, these are $(25^{2})\times(25^{2})=625\times 625$
correlation matrices.

It can be seen from figure 2 that all the main structures of the theoretical
correlation matrix (left panel) are reflected in the sample correlation matrix
constructed from our 1000 Gaussian simulations (right panel). Notice, in
particular, the significant levels of correlation (and, therefore, covariance)
between different different pairs of slices. We conclude by pointing out that
the exercise above has served to validate the expressions found in previous
sections for the angular power spectrum and for the covariance matrix of the
power spectra. As a consequence, we also indirectly checked our expression for
the Fisher matrix of the angular power spectra, Eq. (2.32).

##  4 The harmonic Fisher matrix in redshift space

We will now take into account redshift space distortions [61] in the linear
regime of structure formation, and subsequently compute the Fisher matrix
directly in redshift space.

Analogously to the real space case, the two-point correlation function is
defined here as the expectation value of the counts inside the spherical
shells of radii $\bar{x}$ and $\bar{y}$, but now including the effects of the
peculiar velocities – see, e.g., [22, 24, 69]:

$\displaystyle C_{\ell,s}^{ij}(\bar{x}_{i},\bar{x}_{j})$ $\displaystyle=$
$\displaystyle\;\bar{n}^{i}_{\bar{x}_{i}}\,\bar{n}^{j}_{\bar{x}_{j}}\int_{\bar{x}_{i}}dx_{i}\,x_{i}^{2}\int_{\bar{x}_{j}}dx_{j}\,x_{j}^{2}$
$\displaystyle\quad\times\frac{2}{\pi}\int_{0}^{\infty}dk\,k^{2}\left[j_{\ell}(kx_{i})-\beta^{i}j_{\ell}^{\prime\prime}(kx_{i})\right]\left[j_{\ell}(kx_{j})-\beta^{j}j_{\ell}^{\prime\prime}(kx_{j})\right]P^{ij}(k)\;,$

with $\beta^{i}=f(z_{i})/b^{i}$, where $f(z)$ is the matter growth rate. When
$\beta^{i}\to 0$ the effects of the peculiar velocities are erased, and this
expression reduces to the real space correlation function given by Eq. (3.5).
The harmonic data covariance in redshift space is therefore a simple
generalization of Eq. (3.9):

$\displaystyle\Gamma_{\ell,s}^{ij}(\bar{x}_{i},\bar{x}_{j})$ $\displaystyle=$
$\displaystyle\bar{n}^{i}\bar{n}^{j}\int_{\bar{x}_{i}}dx_{i}\,x_{i}^{2}\,\int_{\bar{x}_{j}}dx_{j}\,x_{j}^{2}\,\frac{2}{\pi}\int_{0}^{\infty}dk\,k^{2}\;\Big{\\{}j_{\ell}(kx_{i})j_{\ell}(kx_{j})\frac{\delta_{ij}}{\sqrt{\bar{n}^{i}\,\bar{n}^{j}}}$
$\displaystyle\;+\left[j_{\ell}(kx_{i})-\beta^{i}j_{\ell}^{\prime\prime}(kx_{i})\right]\left[j_{\ell}(kx_{j})-\beta^{j}j_{\ell}^{\prime\prime}(kx_{j})\right]P^{ij}(k)\Big{\\}}\;.$

Just as was the case in real space, the challenge is to invert this data
covariance. We will show in the next section how to invert Eq. (4) in the same
sense of Eq. (2.21), that is:

$\sum_{j}\sum_{\bar{x}_{j}}\left[\Gamma^{ij}_{\ell,s}(\bar{x}_{i},\bar{x}_{j})\right]^{-1}\Gamma^{ji^{\prime}}_{\ell,s}(\bar{x}_{j},\bar{x}_{i^{\prime}})=\delta_{ii^{\prime}}\delta_{\bar{x}_{i}\,\bar{x}_{i^{\prime}}}\;.$
(4.3)

In order to accomplish this inversion, and in analogy to what was done in real
space, we will make use of analytical solutions to the integrals of products
of spherical Bessel functions that appear in Eq. (4). The term without any
derivatives is simply the closure relation, Eq. (3.8), already used in the
real space case. The challenge lies with the terms involving the second
derivatives.

We start by recalling that the spherical Bessel functions satisfy the
differential equation:

$z^{2}j_{\ell}^{\prime\prime}(z)+2zj_{\ell}^{\prime}(z)+\big{[}z^{2}-\ell(\ell+1)\big{]}j_{\ell}(z)=0\;.$
(4.4)

We can combine this with the recurrence relation:

$j_{\ell}^{\prime}(z)=\frac{\ell}{z}j_{\ell}(z)-j_{\ell+1}(z)\;,$ (4.5)

and rewrite the second derivative of the spherical Bessel function as:

$j_{\ell}^{\prime\prime}(z)=\frac{1}{z^{2}}\bigg{\\{}\big{[}\ell^{2}-\ell-z^{2}\big{]}j_{\ell}(z)+2\,z\,j_{\ell+1}(z)\bigg{\\}}\;.$
(4.6)

At this point it is important to draw a distinction between our path and that
chosen since Fisher, Scharf and Lahav [22] first wrote an expression for the
angular power spectrum in redshift space. After that seminal paper, many
subsequent others, e.g. Ref. [29], expressed the second derivative of the
spherical Bessel function in terms of Bessel functions of different orders:

$j_{\ell}^{\prime\prime}(z)=\frac{\ell(\ell-1)}{(2\ell+1)(2\ell-1)}j_{\ell-2}(z)+\frac{(\ell+1)(\ell+2)}{(2\ell+1)(2\ell+3)}j_{\ell+2}(z)-\frac{(2\ell^{2}+2\ell-1)}{(2\ell+3)(2\ell-1)}j_{\ell}(z)\;.$
(4.7)

We, on the other hand, employ the expression given by Eq. (4.6), which, as
shown in Appendix C, leads to analytical integrals.

There are two relevant integrals, the first one involving the combination of
spherical Bessel functions:

$\int
dk\,k^{2}j_{\ell}(k{x})j_{\ell}^{\prime\prime}(k{y})=-\frac{\pi}{2{x}^{2}}\delta_{D}({x}-{y})+H_{\ell}({x},{y})\;,$
(4.8)

where, as shown in Appendix C:

$H_{\ell}({x},{y})\,=\,\frac{\pi}{2y^{2}}\begin{dcases}\frac{1}{x}\left[\frac{2(2\ell+1)+\ell(\ell-1)}{2(2\ell+1)}\right]&\text{if
}{x}={y}\\\\[12.0pt]
\frac{1}{{x}}\left(\frac{{y}}{{x}}\right)^{\ell}\left[\frac{2(2\ell+1)+\ell(\ell-1)}{(2\ell+1)}\right]&\text{if
}{x}>{y}\\\\[12.0pt]
\frac{1}{{y}}\frac{\ell(\ell-1)}{(2\ell+1)}\left(\frac{{x}}{{y}}\right)^{\ell}&\text{if
}{x}<{y}\;.\end{dcases}$ (4.9)

Although this expression shows a discontinuity, that does not turn out to be a
problem since, as we shall see later on, it appears twice, with interchanged
arguments. In particular, the symmetric combination of this expression appears
when we consider auto-correlations of the same tracer, and in that case the
definition above simplifies to a manifestly continuous expression:

$\frac{1}{2}\left[H_{\ell}({x},{y})+H_{\ell}({y},{x})\right]=\frac{\pi}{2\,r_{>}}\left(\frac{r_{<}}{r_{>}}\right)^{\ell}\left[\frac{2(2\ell+1)r_{<}^{2}+\ell(\ell-1)({x}^{2}+{y}^{2})}{2(2\ell+1){x}^{2}{y}^{2}}\right]\;.$
(4.10)

The second relevant integral involves two second derivatives of spherical
Bessel functions of the same order, and as also shown in Appendix C we have:

$\int_{0}^{\infty}dk\,k^{2}\,j^{\prime\prime}_{\ell}(k{x})j^{\prime\prime}_{\ell}(k{y})=\frac{\pi}{2}\frac{\delta_{D}({x}-{y})}{{x}^{2}}+G_{\ell}({x},{y})\;,$
(4.11)

where:

$\displaystyle G_{\ell}({x},{y})$ $\displaystyle\equiv$
$\displaystyle\frac{\pi}{2}\frac{r_{<}^{\ell}}{r_{>}^{\ell+1}}\left\\{\left[\frac{4}{(2\ell+3)}-2\right]\frac{1}{r_{>}^{2}}\right.$
$\displaystyle+\frac{2\,\ell(\ell-1)}{(2\ell+1)(2\ell+3)}\frac{1}{r_{>}^{2}}-\frac{\ell(\ell-1)}{(2\ell+1)}\left(\frac{1}{r_{>}^{2}}+\frac{1}{r_{<}^{2}}\right)$
$\displaystyle+\frac{2\,\ell(\ell-1)}{r^{2}_{<}}\left[\frac{1}{(2\ell+1)}-\frac{1}{(2\ell+3)}\left(\frac{r_{<}}{r_{>}}\right)^{2}\right]$
$\displaystyle\left.+2^{\ell}[\ell(\ell-1)]^{2}\frac{(\ell+1)!(2\ell-3)!!}{(2\ell+3)!}\frac{1}{r_{<}^{2}}\left[2\ell+3-(2\ell-1)\left(\frac{r_{<}}{r_{>}}\right)^{2}\right]\right\\}\;.$

The relevance of these analytic integrals will become clear in the next
Section.

###  4.1 Inverting the data covariance in redshift space

We are now in a position to show how to invert the expression for the harmonic
data covariance in redshift space. The reader who is uninterested in the
details of this calculation can jump to the final results, which are presented
in Section 4.2.

We start this calculation by rewriting the data covariance:

$\displaystyle\Gamma_{\ell,s}^{ij}(\bar{x}_{i},\bar{x}_{j})$ $\displaystyle=$
$\displaystyle\bar{n}^{i}\,\bar{n}^{j}\int_{\bar{x}_{i}}dx_{i}\,x_{i}^{2}\int_{\bar{x}_{j}}dx_{j}\,x_{j}^{2}\frac{2}{\pi}\int_{0}^{\infty}dk\,k^{2}\Bigg{[}j_{\ell}(kx_{i})j_{\ell}(kx_{j})\frac{\delta_{ij}}{\sqrt{\bar{n}^{i}\,\bar{n}^{j}}}$
$\displaystyle\,+\left[j_{\ell}(kx_{i})-\beta^{i}j_{\ell}^{\prime\prime}(kx_{i})\right]\left[j_{\ell}(kx_{j})-\beta^{j}j_{\ell}^{\prime\prime}(kx_{j})\right]P^{ij}(k)\Bigg{]}$
$\displaystyle\approx$
$\displaystyle\bar{N}^{i}\,\delta_{ij}\,\delta_{\bar{x}_{i}\,\bar{x}_{j}}$
$\displaystyle\,+\,\bar{N}^{i}\bar{N}^{j}\frac{2}{\pi}\int_{0}^{\infty}dk\,k^{2}\left[j_{\ell}(\bar{x}_{i})-\beta^{i}j_{\ell}^{\prime\prime}(k\bar{x}_{i})\right]\left[j_{\ell}(k\bar{x}_{j})-\beta^{j}j_{\ell}^{\prime\prime}(k\bar{x}_{j})\right]P^{ij}(k)\;,$

where $\bar{N}^{i}=\bar{n}^{i}(\bar{x}_{i})\Delta V_{\bar{x}_{i}}/4\pi$
denotes the mean angular number density of tracers $i$ in the radial bin
$\bar{x}_{i}$. Also recall that $\delta_{\bar{x}_{i}\,\bar{x}_{j}}=1$ if the
two shells (radial bins) are the same, $\bar{x}_{i}=\bar{x}_{j}$, and 0
otherwise. The approximation taken from the first to the second line is
similar to that which was discussed in Section 2.1 and was shown in figure 4.
Our results do not rely on this approximation, but we use it here just in
order to shorten our expressions.

We will now simplify the nature of our problem by using the condensed
notation:

$\frac{\pi}{2}\frac{\delta_{D}(k-k^{\prime})}{k^{2}}\rightarrow\delta^{kk^{\prime}}_{F}\qquad,\qquad\frac{2}{\pi}\int
dk\,k^{2}\rightarrow\sum_{[k]}\;,$ (4.14)

and define the normalized data covariance:

$\Gamma^{ij}_{\ell,s}\quad\rightarrow\quad\tilde{\Gamma}^{ij}_{\ell,s}\,\equiv\,\frac{1}{\sqrt{\bar{N}^{i}\,\bar{N}^{j}}}\,\Gamma^{ij}_{\ell,s}\;.$
(4.15)

Now let’s define the auxiliary “dual vectors”:

$\displaystyle V^{i}_{\ell}(k,\bar{x}_{i})$ $\displaystyle\equiv$
$\displaystyle\sqrt{\bar{n}^{i}\,\frac{\Delta
V_{\bar{x}_{i}}}{4\pi}\,P^{ii}(k)}\left[j_{\ell}(k\bar{x}_{i})-\beta^{i}j_{\ell}^{\prime\prime}(k\bar{x}_{i})\right]$
$\displaystyle=$ $\displaystyle\sqrt{\frac{\Delta
V_{\bar{x}_{i}}}{4\pi}\,\mathcal{P}^{ii}(k)}\left[j_{\ell}(k\bar{x}_{i})-\beta^{i}j_{\ell}^{\prime\prime}(k\bar{x}_{i})\right]\;,$

where we used the definition of the clustering strength introduced in Eq.
(3.13), i.e., ${\cal{P}}^{ij}(k)=\sqrt{\bar{n}^{i}\bar{n}^{j}}\,P^{ij}(k)$.
Using the associative property of the spectra
$P^{ii}(k)P^{jj}(k^{\prime})=P^{ij}(k)P^{ij}(k^{\prime})$ (which is valid at
least to first approximation, in the linear regime) allows us to write the
data covariance as:

$\displaystyle\tilde{\Gamma}^{ij}_{\ell,s}(\bar{x}_{i},\bar{x}_{j})$
$\displaystyle=$
$\displaystyle\delta_{ij}\delta_{\bar{x}_{i}\,\bar{x}_{j}}+\sum_{[k]}V^{i}_{\ell}(k,\bar{x}_{i})\,V^{j}_{\ell}(k,\bar{x}_{j})$
$\displaystyle=$
$\displaystyle\delta_{ij}\delta_{\bar{x}_{i}\,\bar{x}_{j}}+\sum_{[k_{i}]}\sum_{[k_{j}]}V^{i}_{\ell}(k_{i},\bar{x}_{i})\,\delta^{k_{i}\,k_{j}}_{F}\,V^{j}_{\ell}(k_{j},\bar{x}_{j})$

We also define, in analogy with the expression above, the normalized harmonic
data covariance in the conjugate (Fourier) representation as:

$\displaystyle\tilde{\Gamma}^{ij}_{\ell,s}(k_{i},k_{j})$ $\displaystyle=$
$\displaystyle\delta_{ij}\delta^{k_{i}\,k_{j}}_{F}+\sum_{\bar{x}}V^{i}_{\ell}(k_{i},\bar{x})\,V^{j}_{\ell}(k_{j},\bar{x})$
$\displaystyle=$
$\displaystyle\delta_{ij}\delta^{k_{i}\,k_{j}}_{F}+\sum_{\bar{x}_{i}}\sum_{\bar{x}_{j}}V^{i}_{\ell}(k_{i},\bar{x}_{i})\,\delta_{\bar{x}_{i}\bar{x}_{j}}V^{j}_{\ell}(k_{j},\bar{x}_{j})\;.$

Just as it was in the case of real space, discussed in Section 3.2, these
expressions are instrumental for the inversion of the harmonic data covariance
in redshift space. In real space the Fourier covariance is perfectly diagonal
[see Eq. (3.16)], but redshift-space distortions introduce some off-diagonal
structures. However, since those off-diagonal structures in Fourier space are
relatively small, the Fourier data covariance is dominated by the diagonal
terms, and it becomes feasible to invert that object. All we have to do then
is to transform the Fourier inverse covariance back to configuration space to
obtain the inverse covariance in terms of the spherical shells.

In order to invert the expression in Eq. (4.1) we can employ the Woodbury
identity:

$\begin{split}M&=A+XBY^{tr}\\\ \Rightarrow\quad
M^{-1}&=A^{-1}-A^{-1}X(B^{-1}+Y^{tr}A^{-1}X)^{-1}Y^{tr}A^{-1}\;,\end{split}$
(4.19)

Substituting
$A\rightarrow\delta_{ij}\delta_{\bar{x}_{i}\,\bar{x}_{j}},\,B\rightarrow\delta_{ij}\delta^{k_{i}\,k_{j}}_{F}$
and $X,Y\rightarrow
V^{i}_{\ell}(k_{i},\bar{x}_{i}),V^{j}_{\ell}(k_{j},\bar{x}_{j})$, we obtain:

$\displaystyle[\tilde{\Gamma}^{ij}_{\ell,s}(\bar{x}_{i},\bar{x}_{j})]^{-1}$
$\displaystyle=$ $\displaystyle\delta_{ij}\delta_{\bar{x}_{i}\,\bar{x}_{j}}$
$\displaystyle-\sum_{[k_{i}]}\sum_{[k_{j}]}V^{i}_{\ell}(k_{i},\bar{x}_{i})\left[\delta_{ij}\delta^{k_{i}\,k_{j}}_{F}+\sum_{\bar{x}}V^{i}_{\ell}(k_{i},\bar{x})V^{j}_{\ell}(k_{j},\bar{x})\right]^{-1}V^{j}_{\ell}(k_{j},\bar{x}_{j})\;.$

Notice that the matrix inside the square brackets is just the Fourier space
harmonic covariance, $\tilde{\Gamma}^{ij}_{\ell}(k_{i},k_{j})$, defined in Eq.
(4.1). Hence, if we can invert $\tilde{\Gamma}^{ij}_{\ell}(k_{i},k_{j})$ then
we can also invert $\Gamma^{ij}_{\ell}(\bar{x}_{i},\bar{x}_{j})$:

$[\tilde{\Gamma}^{ij}_{\ell,s}(\bar{x}_{i},\bar{x}_{j})]^{-1}=\delta_{ij}\delta_{\bar{x}_{i}\,\bar{x}_{j}}-\sum_{[k_{i}]}\sum_{[k_{j}]}V^{i}_{\ell}(k_{i},\bar{x}_{i})\left[\tilde{\Gamma}^{ij}_{\ell}(k_{i},k_{j})\right]^{-1}V^{j}_{\ell}(k_{j},\bar{x}_{j})\;.$
(4.21)

Therefore, we have reduced the problem of inverting the multi-tracer data
covariance to the problem of inverting its conjugate
$\tilde{\Gamma}^{ij}_{\ell}(k_{i},k_{j})$. Returning to Eq. (4.1), we now
write it in full:

$\displaystyle\tilde{\Gamma}^{ij}_{\ell}(k_{i},k_{j})$ $\displaystyle=$
$\displaystyle\delta_{ij}\delta^{k_{i}\,k_{j}}_{F}+\sqrt{\mathcal{P}^{ii}(k_{i})\,\mathcal{P}^{jj}(k_{j})}$
$\displaystyle\times\sum_{\bar{x}}\int_{\bar{x}}dx\,x^{2}\left[j_{\ell}(k_{i}x)-\beta^{i}\,j_{\ell}^{\prime\prime}(k_{i}x)\right]\left[j_{\ell}(k_{j}x)-\beta^{j}\,j_{\ell}^{\prime\prime}(k_{j}x)\right]$
$\displaystyle=$
$\displaystyle\delta^{k_{i}\,k_{j}}_{F}\left[\delta_{ij}+(1+\beta^{i})(1+\beta^{j})\,\sqrt{\mathcal{P}^{ii}(k_{i})\,\mathcal{P}^{jj}(k_{j})}\right]$
$\displaystyle\quad+\sqrt{\mathcal{P}^{ii}(k_{i})\,\mathcal{P}^{jj}(k_{j})}\left[\beta^{i}\,\beta^{j}\,G_{\ell}(k_{i},k_{j})-\beta^{i}\,H_{\ell}(k_{i},k_{j})-\beta^{j}\,H_{\ell}(k_{j},k_{i})\right]\;,$

where in the last line of the equations above we used the fact that
$\sum_{\bar{x}}\int_{\bar{x}}dx=\int_{0}^{\infty}dx$, as well as the
definitions of Eqs. (4.11) and (4.8). Notice that here we used again the fact
that correlations (or expectation values) should be evaluated at some common t
= constant hypersurface (see figure 1 and the discussion in Section 2), and
the sum over radial bins in the expression above, combined with the radial
integrals inside the bins, leads to an integration over the entire range of
the radial coordinate ($x$) on that hypersurface, while the quantities
$\bar{n}^{i}$, $b^{i}$, $\beta^{i}=b^{i}/f(z_{i})$, etc., are kept fixed in
that expression. As a result, the exact same expressions $G_{\ell}$ and
$H_{\ell}$ of Eqs. (4.11) and (4.8) that appeared in configuration space show
up again, but now in Fourier space. Notice also that the Fourier harmonic
covariance above is an entirely different object compared with the angular
power spectrum that follows from taking the expectation value of the spherical
modes in a Fourier-Bessel expansion over the past light cone: in our case, the
data covariance has support on the common $t=\,$ constant hypersurface where
the expectation values are evaluated, not on the past light cone.

At this point it is convenient to collect the terms in the last line of Eq.
(4.1) in the Fourier (or, more accurately, radial) mixing matrix:

$T^{ij}_{\ell}(k_{i},k_{j})\equiv\sqrt{{\cal{P}}^{ii}(k_{i}){\cal{P}}^{jj}(k_{j})}\left[\beta^{i}\beta^{j}\,G_{\ell}(k_{i},k_{j})-\beta^{i}H_{\ell}(k_{i},k_{j})-\beta^{j}H_{\ell}(k_{j},k_{i})\right]\;.$
(4.23)

Eq. (4.1) also suggests that we should define the redshift-space clustering
strength:

$\mathcal{P}^{ij}_{(s)}\,\equiv\,(1+\beta^{i})(1+\beta^{j})\mathcal{P}^{ij}\;.$
(4.24)

Recall that, in the linear regime, the clustering strength can be factored,
${\cal{P}}^{ii}(k_{i}){\cal{P}}^{jj}(k_{j})={\cal{P}}^{ij}(k_{i}){\cal{P}}^{ij}(k_{j})$,
and the definition above implies that the same applies for
$\mathcal{P}^{ij}_{(s)}$. With the help of these expressions we can write
finally a simple, analytical expression for the harmonic data covariance in
redshift space:

$\displaystyle\tilde{\Gamma}_{\ell}^{ij}(k_{i},k_{j})$ $\displaystyle=$
$\displaystyle\delta^{k_{i}\,k_{j}}_{F}\left[\delta_{ij}+\mathcal{P}^{ij}_{(s)}\right]+T^{ij}_{\ell}(k_{i},k_{j})$
(4.25) $\displaystyle=$
$\displaystyle\delta^{k_{i}\,k_{j}}_{F}\left[\delta_{ij}+\mathcal{P}^{ij}_{(s)}(k_{i})+\lambda_{\ell}^{ij}(k_{i})\right]+\tilde{T}^{ij}_{\ell}(k_{i},k_{j})$
$\displaystyle=$
$\displaystyle\delta^{k_{i}\,k_{j}}_{F}\mathcal{S}^{ij}_{\ell}(k_{i})+\tilde{T}^{ij}_{\ell}(k_{i},k_{j})\;,$

where $\lambda^{ij}_{\ell}(k)\sim T^{ij}_{\ell}(k_{i},k_{i})$ carries the
diagonal of the radial mixing matrix, and $\tilde{T}^{ij}_{\ell}(k_{i},k_{j})$
has only off-diagonal ($k_{i}\neq k_{j}$) terms, and we defined:

$\mathcal{S}^{ij}_{\ell}(k)=\delta^{ij}+\lambda^{ij}_{\ell}+\mathcal{P}^{ij}_{(s)}\;.$
(4.26)

Notice that both $\mathcal{P}^{ij}$ and $\lambda^{ij}_{\ell}$ are adimensional
– to see that this is the case for the latter, see Eq. (4.30).

Figure 3: Radial (Fourier) mixing matrix $T_{\ell}(k_{i},k_{j})$, normalized
by its diagonal for visualization purposes. Top: mixing matrix for $\ell=10$
(left) and $\ell=80$ (right). Bottom: rows of the normalized mixing matrices
(denoted as the dashed lines of the two top plots), with fixed $k_{j}=0.04$,
0.08, 0.16, 0.32 and 0.64 $h$ Mpc-1 (from left to right, respectively: dark
red, red, orange, light blue and blue). Notice that for $k_{i}=k_{j}$ the
normalized mixing matrix is always equal to 1. The different lines
corresponds, from lightest to heaviest, to the multipoles $\ell=10$, 20, 40
and 80.

At this point it is useful to examine the structure of the radial mixing
matrix $T_{\ell}^{ij}(k_{i},k_{j})$ defined above, in Eq. (4.23). As it turns
out, this matrix is, to an excellent approximation, dominated by the diagonal
terms $k_{i}=k_{j}$. In the top row of figure 3 we plot that matrix
(normalized by the diagonal for visualization purposes) in the case of a
single tracer, for the multipoles $\ell=10$ (left panel) and $\ell=80$
(right), and using $\beta=1$. These two plots show clearly that the matrix
elements decay as we move away from the diagonal. In the lower plot of figure
3 we show the lines (of varying $k_{i}$) of the (normalized) radial mixing
matrices with $\ell=10,20,40$ and $80$, for fixed values $k_{j}=0.04$, $0.08$,
$0.16$, $0.32$ and $0.64\,h$ Mpc-1. The log-scale of this bottom plot makes it
clear that off-diagonal matrix elements decay as a power-law that depends on
$\ell$. It is crucial to note that the mixing matrix normalized by the
diagonal,
$T^{ij}_{\ell}(k_{i},k_{j})/\sqrt{T^{ii}_{\ell}(k_{i},k_{i})\,T^{jj}_{\ell}(k_{j},k_{j})}$,
has no dependence at all on the shape of the power spectrum or on the number
of tracers: it is determined entirely by the functions $H_{\ell}$ and
$G_{\ell}$, as well as the amplitudes, $\beta^{i}=f(\bar{x}_{i})/b^{i}$. In
conclusion, independently of the shape of the spectrum, and for any number of
tracers, for all but the lowest multipoles it should be a very good
approximation to discard the non-diagonal terms of the matrix
$T_{\ell}^{ij}(k_{i},k_{j})$ – i.e., to neglect
$\tilde{T}^{ij}_{\ell}(k_{i},k_{j})$ altogether.

We can now proceed to computing the inverse of the data covariance in harmonic
Fourier space given in Eq. (4.25). First, notice that:

$(1+\mathcal{P}+\lambda)^{-1}=(1+\mathcal{P})^{-1}\left[1+(1+\mathcal{P})^{-1}\lambda\right]^{-1}\;.$

In a similar way to what was done in real space, in Section 3.2, using the
Woodbury identity leads to the inversion:

$\left[\delta^{ij}+\mathcal{P}^{ij}_{(s)}\right]^{-1}=\delta^{ij}-\frac{\mathcal{P}^{ij}_{(s)}}{1+\mathcal{P}_{(s)}}\;,$
(4.27)

where $\mathcal{P}_{(s)}=\sum_{i}\mathcal{P}^{ii}_{(s)}$. The inverse of
$\mathcal{S}^{ij}_{\ell}$ is therefore given by:

$\displaystyle[\mathcal{S}^{ij}_{\ell}]^{-1}$ $\displaystyle=$
$\displaystyle\sum_{l}\left[\delta^{il}-\frac{\mathcal{P}^{il}_{(s)}}{1+\mathcal{P}_{(s)}}\right]\left[\delta^{lj}+\sum_{n}\left(\delta^{ln}-\frac{\mathcal{P}^{ln}_{(s)}}{1+\mathcal{P}_{(s)}}\right)\lambda^{nj}_{\ell}\right]^{-1}\;.$
(4.28)

In practical applications we would use finite Fourier bins, and it is
therefore important to understand the relative importance of the redshift-
space corrections $\lambda^{ij}_{\ell}$. From Eq. (4.25) it is clear that the
diagonal of the data covariance carries a Dirac delta function – recall that
in our condensed notation
$\delta^{k_{i}k_{j}}_{F}=(\pi/2)\,\delta_{D}(k_{i}-k_{j})/k_{i}^{2}$. Using
discrete Fourier bins of width $\Delta k$ we have
$\delta_{D}(k_{i}-k_{j})\to\delta_{k_{i}k_{j}}/\Delta k$, and the covariance
is then given by:

$\tilde{\Gamma}_{\ell}^{ij}(k_{i},k_{j})\to\frac{\pi}{2}\frac{\delta_{k_{i}k_{j}}}{k_{i}^{2}\Delta
k}\mathcal{S}^{ij}_{\ell}(k_{i})+\tilde{T}^{ij}_{\ell}(k_{i},k_{j})\;,$ (4.29)

where now we identify:

$\lambda^{ij}_{\ell}(k_{i})\to\frac{2}{\pi}\,k_{i}^{2}\,\Delta
k\,T^{ij}_{\ell}(k_{i},k_{i})\;.$ (4.30)

Since the mixing matrix has dimensions of volume, the diagonal term
$\lambda^{ij}_{\ell}$ is adimensional. The expression above shows that these
diagonal terms are proportional to the Fourier bin width, while from Eq.
(4.23) and the definitions of $H_{\ell}$ and $G_{\ell}$, Eqs. (4) and (4.9),
it is straightforward to see that, up to factors of the redshift-space
parameter $\beta$,
$T_{\ell}^{ij}(k,k)\sim\alpha_{\ell}\,\mathcal{P}^{ij}(k)/k^{3}$, where
$\alpha_{\ell}$ is a coefficient that could be very large for $\ell\gg 1$.
This means that the term
$\lambda^{ij}_{\ell}\sim\alpha_{\ell}\,\mathcal{P}^{ij}(k)\,\Delta k/k$ is in
general not negligible. If in a particular case $\lambda^{ij}_{\ell}$ happens
to be small compared with $\delta^{ij}$ and $\mathcal{P}_{s}^{ij}$ [see Eq.
(4.26)], then the inverse of the matrix $\mathcal{S}^{ij}_{\ell}$, Eq. (4.28),
could be written in terms of a series expansion:

$\displaystyle[\mathcal{S}^{ij}_{\ell}]^{-1}$ $\displaystyle\approx$
$\displaystyle\delta^{ij}-\frac{\mathcal{P}^{ij}_{(s)}}{1+\mathcal{P}_{(s)}}-\sum_{i^{\prime},j^{\prime}}\left(\delta^{ii^{\prime}}-\frac{\mathcal{P}^{ii^{\prime}}_{(s)}}{1+\mathcal{P}_{(s)}}\right)\lambda^{i^{\prime}j^{\prime}}_{\ell}\left(\delta^{j^{\prime}j}-\frac{\mathcal{P}^{j^{\prime}j}_{(s)}}{1+\mathcal{P}_{(s)}}\right)+{\cal{O}}(\lambda^{2})\;.$
(4.31)

In any case, $\mathcal{S}^{ij}_{\ell}(k)$ is a simple $N_{t}\times N_{t}$
object (where $N_{t}$ is the number of tracers), which has a trivial inverse
for each value of $k$ and for each multipole $\ell$.

Once the inverse $[\mathcal{S}^{ij}_{\ell}(k)]^{-1}$ has been determined, we
can finally obtain the inverse covariance. Since the off-diagonal matrix
$\tilde{T}^{ij}(k_{i},k_{j})$ is small compared with the diagonal part of the
data covariance of Eq. (4.25), we can invert that expression in terms of a
series expansion around $\tilde{T}^{ij}(k_{i},k_{j})$:

$[\tilde{\Gamma}_{\ell}^{ij}(k_{i},k_{j})]^{-1}\approx\delta^{k_{i}k_{j}}_{F}[\mathcal{S}^{ij}_{\ell}(k_{i})]^{-1}-\sum_{i^{\prime},j^{\prime}}[\mathcal{S}^{ii^{\prime}}_{\ell}(k_{i})]^{-1}\,\tilde{T}^{i^{\prime}j^{\prime}}_{\ell}(k_{i},k_{j})\,[\mathcal{S}_{\ell}^{j^{\prime}j}(k_{j})]^{-1}+{\cal{O}}(\tilde{T}^{2})\;.$
(4.32)

As discussed above, it is often sufficient to consider only the first term in
the expression above, neglecting the off-diagonal matrix
$\tilde{T}_{\ell}^{ij}$, which leads to corrections of order
$\sim\,10^{-2}-10^{-3}$, depending on the multipole. This is a manifestation
of the fact that the correlation function is much more diagonal in Fourier
space than it is in real space, and is what ultimately allows us to invert the
harmonic data covariance.

Before we go back to the standard notation for the data covariance, it is
useful to assess the accuracy of our approximations, when we take into account
the off-diagonal matrix $\tilde{T}^{ij}_{\ell}(k_{i},k_{j})$ to first order,
as given in Eq. (4.32). Recall that the inverse must be be such that:

$\displaystyle\sum_{j}\sum_{[k_{j}]}\left[\tilde{\Gamma}_{\ell}^{ij}(k_{i},k_{j})\right]^{-1}\tilde{\Gamma}^{ji^{\prime}}_{\ell}(k_{j},k_{i^{\prime}})$
$\displaystyle=$
$\displaystyle\sum_{j}\frac{2}{\pi}\int_{0}^{\infty}dk_{j}\,k_{j}^{2}\,\left[\tilde{\Gamma}^{ij}_{\ell}(k_{i},k_{j})\right]^{-1}\tilde{\Gamma}^{ji^{\prime}}_{\ell}(k_{j},k_{i^{\prime}})$
(4.33) $\displaystyle=$
$\displaystyle\delta^{k_{i}\,k_{i^{\prime}}}_{F}\delta_{ii^{\prime}}=\frac{\pi}{2}\frac{\delta_{D}(k_{i}-k_{i^{\prime}})}{k_{i}^{2}}\delta_{ii^{\prime}}\;.$

In order to check the accuracy of the approximation in Eq. (4.32) we
considered the case of a single tracer, such that:

$\mathcal{S}_{\ell}=1+\lambda_{\ell}+\mathcal{P}_{(s)}\quad,\quad\mathcal{S}^{-1}_{\ell}=\frac{1}{1+\lambda_{\ell}+\mathcal{P}_{(s)}}\;.$

Substituting Eq. (4.32) into Eq. (4.33) then leads to a residue which is of
second order in $\tilde{T}$. We computed numerically these
$\mathcal{O}(\tilde{T}^{2})$ corrections in discrete Fourier bins assuming
realistic scenarios for the matter power spectrum and for the redshift-space
distortions, and the result is that these extra terms are negligible,
typically of the order $\sim\,10^{-6}$ compared with the lower-order terms in
the expansion. We checked the accuracy of this approximation for $\ell\lesssim
100$ and scales $k\lesssim 1\;h\;{\rm Mpc}^{-1}$.

To summarize the results of this Section, we found that the inverse of the
normalized harmonic data covariance is given by:

$\displaystyle\left[\tilde{\Gamma}^{ij}_{\ell,s}(\bar{x}_{i},\bar{x}_{j})\right]^{-1}$
$\displaystyle\simeq$
$\displaystyle\delta_{ij}\delta_{\bar{x}_{i}\bar{x}_{j}}$ $\displaystyle-$
$\displaystyle\left(\frac{2}{\pi}\right)^{2}\frac{\Delta
V_{\bar{x}_{i}}}{4\pi}\frac{\Delta V_{\bar{x}_{j}}}{4\pi}\int
dk_{i}\,k_{i}^{2}\int
dk_{j}\,k_{j}^{2}\,\sqrt{\mathcal{P}^{ii}(k_{i})\,\mathcal{P}^{jj}(k_{j})}$
$\displaystyle\qquad\times\left[j_{\ell}(k_{i}\bar{x}_{i})-\beta^{i}\,j_{\ell}^{\prime\prime}(k_{i}\bar{x}_{i})\right]\left[j_{\ell}(k_{j}\bar{x}_{j})-\beta^{j}\,j_{\ell}^{\prime\prime}(k_{j}\bar{x}_{j})\right]$
$\displaystyle\qquad\times\left[\frac{\pi}{2}\frac{\delta(k_{i}-k_{j})}{k_{i}^{2}}[\mathcal{S}^{ij}(k_{i})]^{-1}-\sum_{i^{\prime},j^{\prime}}[\mathcal{S}^{ii^{\prime}}(k_{i})]^{-1}\tilde{T}^{i^{\prime}j^{\prime}}_{\ell}(k_{i},k_{j})[\mathcal{S}^{j^{\prime}j}(k_{j})]^{-1}\right]\;.$

All that is left to do is to plug the normalization factor of Eq. (4.15) back
in this expression.

###  4.2 Summary of the main results

In the previous Section we found that the inverse of the configuration-space
harmonic data covariance including redshift-space distortions, Eq. (3.11), is
given by:

$\displaystyle[{\Gamma}^{ij}_{\ell,s}(\bar{x}_{i},\bar{x}_{j})]^{-1}$
$\displaystyle\simeq$
$\displaystyle\frac{1}{\bar{N}^{i}}\,\delta_{ij}\,\delta_{\bar{x}_{i}\bar{x}_{j}}$
$\displaystyle-\frac{2}{\pi}\int
dk\,k^{2}\;\left[j_{\ell}(k\bar{x}_{i})-\beta^{i}j_{\ell}^{\prime\prime}(k\bar{x}_{i})\right]\left[j_{\ell}(k\bar{x}_{j})-\beta^{j}j_{\ell}^{\prime\prime}(k\bar{x}_{j})\right]\;P^{ij}(k)\,[\mathcal{S}^{ij}(k)]^{-1}$
$\displaystyle+\left(\frac{2}{\pi}\right)^{2}\int dk_{i}\,k_{i}^{2}\int
dk_{j}\,k_{j}^{2}\;\left[j_{\ell}(k_{i}\bar{x}_{i})-\beta^{i}j_{\ell}^{\prime\prime}(k_{i}\bar{x}_{i})\right]\left[j_{\ell}(k_{j}\bar{x}_{j})-\beta^{j}j_{\ell}^{\prime\prime}(k_{j}\bar{x}_{j})\right]$
$\displaystyle\quad\times\sqrt{P^{ii}(k_{i})\,P^{jj}(k_{j})}\sum_{i^{\prime},j^{\prime}}[\mathcal{S}^{ii^{\prime}}(k_{i})]^{-1}\tilde{T}^{i^{\prime}j^{\prime}}_{\ell}(k_{i},k_{j})[\mathcal{S}^{j^{\prime}j}(k_{j})]^{-1}\;+\ldots\;,$

where $\mathcal{S}^{ij}(k)$ was defined in Eq. (4.25), and the radial mixing
matrix $T^{ij}(k_{i},k_{j})$ was defined in Eq. (4.23).

Notice that for $\beta^{i}\to 0$ we also get $T_{\ell}^{ij}\to 0$ (so
$\tilde{T}^{ij}_{\ell}\to 0$ and $\lambda^{i}_{\ell}\to 0$), and we recover
the real-space result of Eq. (3.11). Since $T^{ij}_{\ell}(k_{i},k_{j})$ and
$S^{ij}_{\ell}(k)$ are given in terms of analytical expressions — see Eqs.
(4.23) and (4.26) — in a first approximation the inverse covariance above can
be computed in terms of $3\times N_{t}(N_{t}+1)/2$ integrals involving Bessel
functions, which can be handled either by taking the Limber approximation or
numerically, using, e.g., FFTLog [70, 71, 72]. If a more accurate calculation
is needed, then in addition to those single integrals we must also compute
$3\times N_{t}(N_{t}+1)/2$ double integrals. As discussed in the previous
section, it is often a good approximation to neglect the matrix
$\tilde{T}^{ij}_{\ell}$, keeping only the single Fourier integral in Eq.
(4.2). Although this approximation is excellent for high multipoles
($\ell\gtrsim 10$), for very low multipoles one may need to perform the double
integral in Eq. (4.2).

Finally, after finding at last an expression for the inverse data covariance
in configuration space, including redshift-space distortions, we find
ourselves in a position to go back to the Fisher matrix for the multi-tracer
angular power spectra. Since the general expressions derived in Sections 2 and
3 still hold, the Fisher matrix is given by the same expression as Eq. (2.32):

$F[C^{[ij]}_{\ell,s},C^{[i^{\prime}j^{\prime}]}_{\ell,s}]=\frac{2\ell+1}{4}(2-\bar{\delta}_{ij})(2-\bar{\delta}_{i^{\prime}j^{\prime}})\left\\{[{\Gamma}^{ii^{\prime}}_{\ell,s}]^{-1}[{\Gamma}^{jj^{\prime}}_{\ell,s}]^{-1}+[{\Gamma}^{ij^{\prime}}_{\ell,s}]^{-1}[{\Gamma}^{i^{\prime}j}_{\ell,s}]^{-1}\right\\}\;,$
(4.36)

where the inverse data covariance $\left[\Gamma^{ij}_{\ell,s}\right]^{-1}$ is
given in Eq. (4.2). The inverse of the Fisher matrix is formally identical to
the case in real space, and leads to the multi-tracer covariance matrix for
the angular power spectra in redshift space that generalizes Eq. (2.33):

${\rm
Cov}[C^{[ij]}_{\ell,s},C^{[i^{\prime}j^{\prime}]}_{\ell,s}]=\frac{1}{2\ell+1}\left[{\Gamma}^{ii^{\prime}}_{\ell,s}\,{\Gamma}^{jj^{\prime}}_{\ell,s}+{\Gamma}^{ij^{\prime}}_{\ell,s}\,{\Gamma}^{i^{\prime}j}_{\ell,s}\right]\;.$
(4.37)

Of course, this expression can also be derived directly using the Gaussian
approximation for the 4-point function of the angular power spectrum.

As discussed in this Section, the critical step is the computation of the
inverse data covariances $\Gamma^{-1}$ of Eq. (4.36). But the key question is:
why would one resort to these formulas, as opposed to a direct numerical
computation?

Let’s assume that we have 100 redshift bins and three tracers: as an example,
consider a survey that maps red galaxies, blue galaxies and HI in the interval
$0.5\leq z\leq 1.5$ , in bins of width $\Delta z=0.01$, which means comoving
radius bins of widths $\Delta\chi\sim 30/E(z)\,h^{-}1$ Mpc, where
$E(z)=H(z)/H_{0}$. That means $300\times(299)/2=44850$ pairs of tracers and
redshift slices. Naively, 44850 would be the minimum number of mocks (or
simulations) needed to estimate the sample covariance matrix, and depending on
the nature of the tracers, a numerically stable estimation would need
something of the order of $44850^{2}\sim 2\times 10^{9}$ mocks (see also Fig.
2, in Section 3, for an illustration of the type of matrix one would need to
evaluate). However, there may be substantial redundancy in those matrices, so
perhaps a more reasonable assumption would be around $10^{5}$ mocks for a
stable convergence of the sample covariance matrix. The question then is: how
long would it take to produce that number of mocks/simulations (including, of
course, RSDs), using the best available tools today?

The standard code to compute full-sky mocks of the past light-cone is FLASK
[73], or its updated version GLASS 222Tessori, Loureiro & Joachimi, in prep.,
https://glass.readthedocs.io/en/latest/index.html. We have estimated that it
would take approximately $1.6\times 10^{4}$ CPU hours, as well as $\sim 240$
Tb of storage, in order to produce that number of simulations with the given
volume and assuming full sky, with a resolution of $Nside=512$.

On the other hand, assuming we use $\ell_{min}=2$ to $\ell_{max}=500$
(probably an overestimate, since at these radii the scales $k=\ell_{max}/\chi$
are already well into the non-linear regime), we would need to compute
approximately $6\times 5000$ double integrals in Eq. (4.2). With the help of
2DFFTlog [71] it is possible to compute each one of the Bessel integrals of
those double integrals in approximately 2s using a single core in a desktop
computer – notice that this corresponds to all the radial bins at once for a
given multipole $\ell$. We have three of those integrals in Eq. (4.2), which
means a total of approximately 50 hrs on a single core.

##  5 Discussion and conclusions

In this paper we present a semi-analytical expression for the multi-tracer
Fisher matrix of the redshift-space angular power spectra of counts of tracers
in spherical shells. The key step that enabled us to find this expression was
the fact that, according to Eq. (2.13), the Fisher matrix is given by the
square of the inverse data covariance, and, in the linear regime, it is
possible to write an analytical expression for the inverse of the harmonic
data covariance even when we include redshift-space distortions. The main
results are summarized at the end of Section 4, in Eqs. (4.2) and (4.36).

We have not resorted at any time to the Limber approximation: in fact, one of
our main results (a semi-analytic expression for the inverse harmonic data
covariance) follows from the fact that we were able to compute some of the
radial integrals analytically. It was because of those analytic integrals that
we were able to obtain a Fourier-space radial mixing matrix that determines
how the different spherical shells are correlated – see Eq. (4.23).
Nevertheless, our final semi-analytic expression for the inverse harmonic data
covariance in redshift space, Eq. (4.2), would still need to be evaluated
either numerically, or with the help of an approximation scheme such as the
Limber approximation.

It is worth clarifying once again that our formulas are of interest mainly
when a survey has a large number of thin redshift slices, and/or a large
number of tracers. When we have only a few tracers and a small number of
redshift slices, it is often more economical to simply compute the covariance
of the angular power spectra, and then invert that covariance numerically.
However, with a large number of auto- and cross-spectra, the resulting high-
dimensional covariance may be challenging to compute, since numerical
instabilities can make it very difficult to either invert the covariance
matrix, or to even allow for the computation of its Cholesky decomposition.
Therefore, it is in the regime of many tracers and many redshift slices that
our semi-analytic formulas for the Fisher matrix can become useful.

The Fisher matrix for the angular power spectra has a dual use: on one hand,
it is straightforward to incorporate our results into Fisher codes in order to
forecast the constraining power of future surveys, in particular those for
which very large-scale phenomena are critical, such as primordial non-
Gaussianities and relativistic effects. On the other hand, we can employ the
Fisher matrix as an approximation for the uncertainties in the measurements of
the angular power spectra of counts in shells, which can then be used in
Markov-Chain Monte Carlo (MCMC) explorations of the likelihood in real surveys
or in simulations.

A key application of our formulas, which is in fact what will ultimately
justify their use in certain scenarios, is a full comparison between the
constraints derived using angular power spectra, and the constraints derived
using the Fourier power spectra. In a forthcoming paper we will compare the
constraining power of the two approaches, and the biases in the parameters
that either approaches may introduce.

As our final remarks, we point out that our results were derived under some
key assumptions: first, we used isotropic selection functions
$\bar{n}^{i}(\vec{x})\to\bar{n}^{i}(x)$ that cover the full sky, and we
assumed that the radial bins are non-overlapping, with constant tracer
densities inside each radial bin, $\bar{n}^{i}(x)\to\bar{n}^{i}(\bar{x}_{i})$
for all $x\in\bar{x}_{i}$. Relaxing these conditions introduces couplings
between the modes of the angular power spectra, and in full generality those
would need to be treated using numerical methods. However, with regards to the
radial selection function, since our results apply for arbitrarily thin bins,
one can always construct larger bins by summing over the smaller bins with
some weights, and then taking care of the additional covariance through the
appropriate Jacobian. And concerning a mask (or more generically, an angular
selection function), which lead to coupling between the different multipoles,
this can be partially mitigated by constructing bins in $\ell$. Another result
which allows us to sustain the semi-analytical approach is the fact that the
integrals of Bessel functions of different orders, which appear in the context
of mode-coupling, also show a nearly diagonal behavior [74].

And second, we assumed that redshift-space distortions are in the linear
regime – see Eq. (4). But this is no different from what we already face when
using the usual Fourier power spectrum. Therefore, even though these are
strong assumptions, this work constitute a significant first step towards
simplifying the spherical description of galaxy surveys, which should further
facilitate the cosmological use of the angular power spectra as a tool to
study large-scale structure.

## Acknowledgments

We would like to thank David Alonso, Stefano Camera and José Fonseca for
useful insights and comments on an early draft. We also acknowledge the
financial support of FAPESP (R.A.), CNPq (R.A & I.L.T.) and CAPES (J.V.D.F.).
For the purpose of open access, the authors have applied a Creative Commons
Attribution (CC BY) license to any Author Accepted Manuscript version arising
from this submission.

## References

  * [1] D.G. York, J. Adelman, J. Anderson, John E., S.F. Anderson, J. Annis, N.A. Bahcall et al., _The Sloan Digital Sky Survey: Technical Summary_ , _Astron. J._ 120 (2000) 1579 [astro-ph/0006396].
  * [2] DESI Collaboration, A. Aghamousa, J. Aguilar, S. Ahlen, S. Alam, L.E. Allen et al., _The DESI Experiment Part I: Science,Targeting, and Survey Design_ , _arXiv e-prints_ (2016) arXiv:1611.00036 [1611.00036].
  * [3] M. Takada, R.S. Ellis, M. Chiba, J.E. Greene, H. Aihara, N. Arimoto et al., _Extragalactic science, cosmology, and Galactic archaeology with the Subaru Prime Focus Spectrograph_ , _Publ. Astron. Soc. Jap._ 66 (2014) R1 [1206.0737].
  * [4] L. Amendola, S. Appleby, A. Avgoustidis, D. Bacon, T. Baker, M. Baldi et al., _Cosmology and fundamental physics with the Euclid satellite_ , _Living Reviews in Relativity_ 21 (2018) 2 [1606.00180].
  * [5] Ž. Ivezić, S.M. Kahn, J.A. Tyson, B. Abel, E. Acosta, R. Allsman et al., _LSST: From Science Drivers to Reference Design and Anticipated Data Products_ , _Astrophys. J._ 873 (2019) 111 [0805.2366].
  * [6] N. Benitez, R. Dupke, M. Moles, L. Sodre, J. Cenarro, A. Marin-Franch et al., _J-PAS: The Javalambre-Physics of the Accelerated Universe Astrophysical Survey_ , _arXiv e-prints_ (2014) arXiv:1403.5237 [1403.5237].
  * [7] Square Kilometre Array Cosmology Science Working Group, D.J. Bacon, R.A. Battye, P. Bull, S. Camera, P.G. Ferreira et al., _Cosmology with Phase 1 of the Square Kilometre Array Red Book 2018: Technical specifications and performance forecasts_ , _Publ. Astron. Soc. Austr._ 37 (2020) e007 [1811.02743].
  * [8] T. Clifton, P.G. Ferreira, A. Padilla and C. Skordis, _Modified gravity and cosmology_ , _Phys. Rep._ 513 (2012) 1 [1106.2476].
  * [9] D. Munshi, G. Pratten, P. Valageas, P. Coles and P. Brax, _Galaxy clustering in 3D and modified gravity theories_ , _Mon. Not. R. Astron. Soc._ 456 (2016) 1627 [1508.00583].
  * [10] N. Dalal, O. Doré, D. Huterer and A. Shirokov, _Imprints of primordial non-Gaussianities on large-scale structure: Scale-dependent bias and abundance of virialized objects_ , _Phys. Review D_ 77 (2008) 123514 [0710.4560].
  * [11] J. Yoo, _General relativistic description of the observed galaxy power spectrum: Do we understand what we measure?_ , _Phys. Review D_ 82 (2010) 083508 [1009.3021].
  * [12] A. Challinor and A. Lewis, _Linear power spectrum of observed source number counts_ , _Phys. Review D_ 84 (2011) 043516 [1105.5292].
  * [13] C. Bonvin and R. Durrer, _What galaxy surveys really measure_ , _Phys. Review D_ 84 (2011) 063505 [1105.5280].
  * [14] J. Yoo, N. Hamaus, U. Seljak and M. Zaldarriaga, _Going beyond the Kaiser redshift-space distortion formula: A full general relativistic account of the effects and their detectability in galaxy clustering_ , _Phys. Review D_ 86 (2012) 063514 [1109.0998].
  * [15] M. Bruni, R. Crittenden, K. Koyama, R. Maartens, C. Pitrou and D. Wands, _Disentangling non-Gaussianity, bias, and general relativistic effects in the galaxy distribution_ , _Phys. Review D_ 85 (2012) 041301 [1106.3999].
  * [16] L.R. Abramo and D. Bertacca, _Disentangling the effects of Doppler velocity and primordial non-Gaussianity in galaxy power spectra_ , _Phys. Review D_ 96 (2017) 123535 [1706.01834].
  * [17] C. Guandalin, J. Adamek, P. Bull, C. Clarkson, L.R. Abramo and L. Coates, _Observing relativistic features in large-scale structure surveys - I. Multipoles of the power spectrum_ , _Mon. Not. R. Astron. Soc._ 501 (2021) 2547 [2009.02284].
  * [18] K. Yamamoto, M. Nakamichi, A. Kamino, B.A. Bassett and H. Nishioka, _A Measurement of the Quadrupole Power Spectrum in the Clustering of the 2dF QSO Survey_ , _Publ. Astron. Soc. Jap._ 58 (2006) 93 [astro-ph/0505115].
  * [19] D. Bianchi, H. Gil-Marín, R. Ruggeri and W.J. Percival, _Measuring line-of-sight-dependent Fourier-space clustering using FFTs_ , _Mon. Not. R. Astron. Soc._ 453 (2015) L11 [1505.05341].
  * [20] R. Scoccimarro, _Fast estimators for redshift-space clustering_ , _Phys. Review D_ 92 (2015) 083532 [1506.02729].
  * [21] P.J.E. Peebles, _Statistical Analysis of Catalogs of Extragalactic Objects. I. Theory_ , _Astrophys. J._ 185 (1973) 413.
  * [22] K.B. Fisher, C.A. Scharf and O. Lahav, _A spherical harmonic approach to redshift distortion and a measurement of Omega(0) from the 1.2-Jy IRAS Redshift Survey_ , _Mon. Not. R. Astron. Soc._ 266 (1994) 219 [astro-ph/9309027].
  * [23] A.F. Heavens and A.N. Taylor, _A spherical harmonic analysis of redshift space_ , _Mon. Not. R. Astron. Soc._ 275 (1995) 483 [astro-ph/9409027].
  * [24] A.J.S. Hamilton and M. Culhane, _Spherical redshift distortions_ , _Mon. Not. R. Astron. Soc._ 278 (1996) 73 [astro-ph/9507021].
  * [25] F.A. Dahlen and F.J. Simons, _Spectral estimation on a sphere in geophysics and cosmology_ , _Geophysical Journal International_ 174 (2008) 774 [0705.3083].
  * [26] C. Scharf, Y. Hoffman, O. Lahav and D. Lynden-Bell, _Spherical harmonic analysis of IRAS galaxies : implications for the GreatAttractor and cold dark matter._ , _Mon. Not. R. Astron. Soc._ 256 (1992) 229.
  * [27] D. Huterer, L. Knox and R.C. Nichol, _The Angular Power Spectrum of Edinburgh/Durham Southern Galaxy Catalogue Galaxies_ , _Astrophys. J._ 555 (2001) 547 [astro-ph/0011069].
  * [28] C. Blake, P.G. Ferreira and J. Borrill, _The angular power spectrum of NVSS radio galaxies_ , _Mon. Not. R. Astron. Soc._ 351 (2004) 923 [astro-ph/0404085].
  * [29] N. Padmanabhan, D.J. Schlegel, U. Seljak, A. Makarov, N.A. Bahcall, M.R. Blanton et al., _The clustering of luminous red galaxies in the Sloan Digital Sky Survey imaging data_ , _Mon. Not. R. Astron. Soc._ 378 (2007) 852 [astro-ph/0605302].
  * [30] C. Blake, A. Collister, S. Bridle and O. Lahav, _Cosmological baryonic and matter densities from 600000 SDSS luminous red galaxies with photometric redshifts_ , _Mon. Not. R. Astron. Soc._ 374 (2007) 1527 [astro-ph/0605303].
  * [31] S. Ho, A. Cuesta, H.-J. Seo, R. de Putter, A.J. Ross, M. White et al., _Clustering of Sloan Digital Sky Survey III Photometric Luminous Galaxies: The Measurement, Systematics, and Cosmological Implications_ , _Astrophys. J._ 761 (2012) 14 [1201.2137].
  * [32] J. Asorey, M. Crocce, E. Gaztañaga and A. Lewis, _Recovering 3D clustering information with angular correlations_ , _Mon. Not. R. Astron. Soc._ 427 (2012) 1891 [1207.6487].
  * [33] A. Balaguera-Antolínez, M. Bilicki, E. Branchini and A. Postiglione, _Extracting cosmological information from the angular power spectrum of the 2MASS Photometric Redshift catalogue_ , _Mon. Not. R. Astron. Soc._ 476 (2018) 1050 [1711.04583].
  * [34] H.S. Xavier, M.V. Costa-Duarte, A. Balaguera-Antolínez and M. Bilicki, _All-sky angular power spectra from cleaned WISE×SuperCOSMOS galaxy number counts_ , _Journal of Cosmology and Astroparticle Physics_ 2019 (2019) 037 [1812.08182].
  * [35] H. Camacho, N. Kokron, F. Andrade-Oliveira, R. Rosenfeld, M. Lima, F. Lacasa et al., _Dark Energy Survey Year 1 results: measurement of the galaxy angular power spectrum_ , _Mon. Not. R. Astron. Soc._ 487 (2019) 3870 [1807.10163].
  * [36] A. Loureiro, B. Moraes, F.B. Abdalla, A. Cuceu, M. McLeod, L. Whiteway et al., _Cosmological measurements from angular power spectra analysis of BOSS DR12 tomography_ , _Mon. Not. R. Astron. Soc._ 485 (2019) 326 [1809.07204].
  * [37] F. Andrade-Oliveira, H. Camacho, L. Faga, R. Gomes, R. Rosenfeld, A. Troja et al., _Galaxy clustering in harmonic space from the dark energy survey year 1 data: compatibility with real-space results_ , _Mon. Not. R. Astron. Soc._ 505 (2021) 5714 [2103.14190].
  * [38] B. Leistedt, H.V. Peiris, D.J. Mortlock, A. Benoit-Lévy and A. Pontzen, _Estimating the large-scale angular power spectrum in the presence of systematics: a case study of Sloan Digital Sky Survey quasars_ , _Mon. Not. R. Astron. Soc._ 435 (2013) 1857 [1306.0005].
  * [39] L.D. Ferramacho, M.G. Santos, M.J. Jarvis and S. Camera, _Radio galaxy populations and the multitracer technique: pushing the limits on primordial non-Gaussianity_ , _Mon. Not. R. Astron. Soc._ 442 (2014) 2511 [1402.2290].
  * [40] J. Fonseca, S. Camera, M.G. Santos and R. Maartens, _Hunting Down Horizon-scale Effects with Multi-wavelength Surveys_ , _Astrophys. J. Lett._ 812 (2015) L22 [1507.04605].
  * [41] D. Huterer, _Weak lensing and dark energy_ , _Phys. Review D_ 65 (2002) 063001 [astro-ph/0106399].
  * [42] C. Hikage, M. Takada, T. Hamana and D. Spergel, _Shear power spectrum reconstruction using the pseudo-spectrum method_ , _Mon. Not. R. Astron. Soc._ 412 (2011) 65 [1004.3542].
  * [43] A. Loureiro, L. Whittaker, A. Spurio Mancini, B. Joachimi, A. Cuceu, M. Asgari et al., _KiDS & Euclid: Cosmological implications of a pseudo angular power spectrum analysis of KiDS-1000 cosmic shear tomography_, _arXiv e-prints_ (2021) arXiv:2110.06947 [2110.06947].
  * [44] R.E. Upham, M.L. Brown, L. Whittaker, A. Amara, N. Auricchio, D. Bonino et al., _Euclid: Covariance of weak lensing pseudo- $C_{\ell}$ estimates. Calculation, comparison to simulations, and dependence on survey geometry_, _arXiv e-prints_ (2021) arXiv:2112.07341 [2112.07341].
  * [45] B. Leistedt, A. Rassat, A. Réfrégier and J.L. Starck, _3DEX: a code for fast spherical Fourier-Bessel decomposition of 3D surveys_ , _Astron. Astrophys._ 540 (2012) A60 [1111.3591].
  * [46] M. Reinecke and D.S. Seljebotn, _Libsharp - spherical harmonic transforms revisited_ , _Astron. Astrophys._ 554 (2013) A112 [1303.4945].
  * [47] D. Alonso, J. Sanchez, A. Slosar and LSST Dark Energy Science Collaboration, _A unified pseudo-C ℓ framework_, _Mon. Not. R. Astron. Soc._ 484 (2019) 4127 [1809.09603].
  * [48] H.S. Grasshorn Gebhardt and O. Doré, _Fabulous code for spherical Fourier-Bessel decomposition_ , _Phys. Review D_ 104 (2021) 123548 [2102.10079].
  * [49] J. Jasche and B.D. Wandelt, _Bayesian physical reconstruction of initial conditions from large-scale structure surveys_ , _Mon. Not. R. Astron. Soc._ 432 (2013) 894 [1203.3639].
  * [50] J. Alsing, A. Heavens, A.H. Jaffe, A. Kiessling, B. Wandelt and T. Hoffmann, _Hierarchical cosmic shear power spectrum inference_ , _Mon. Not. R. Astron. Soc._ 455 (2016) 4452 [1505.07840].
  * [51] J. Alsing, A. Heavens and A.H. Jaffe, _Cosmological parameters, shear maps and power spectra from CFHTLenS using Bayesian hierarchical inference_ , _Mon. Not. R. Astron. Soc._ 466 (2017) 3272 [1607.00008].
  * [52] D.K. Ramanah, G. Lavaux, J. Jasche and B.D. Wandelt, _Cosmological inference from Bayesian forward modelling of deep galaxy redshift surveys_ , _Astron. Astrophys._ 621 (2019) A69 [1808.07496].
  * [53] G. Lavaux, J. Jasche and F. Leclercq, _Systematic-free inference of the cosmic matter density field from SDSS3-BOSS data_ , _arXiv e-prints_ (2019) arXiv:1909.06396 [1909.06396].
  * [54] N. Porqueres, O. Hahn, J. Jasche and G. Lavaux, _A hierarchical field-level inference approach to reconstruction from sparse Lyman- $\alpha$ forest data_, _Astron. Astrophys._ 642 (2020) A139 [2005.12928].
  * [55] F. Leclercq and A. Heavens, _On the accuracy and precision of correlation functions and field-level inference in cosmology_ , _Mon. Not. R. Astron. Soc._ 506 (2021) L85 [2103.04158].
  * [56] N. Porqueres, A. Heavens, D. Mortlock and G. Lavaux, _Lifting weak lensing degeneracies with a field-based likelihood_ , _Mon. Not. R. Astron. Soc._ 509 (2022) 3194 [2108.04825].
  * [57] J. Jasche and F.S. Kitaura, _Fast Hamiltonian sampling for large-scale structure inference_ , _Mon. Not. R. Astron. Soc._ 407 (2010) 29 [0911.2496].
  * [58] M. Betancourt, _A Conceptual Introduction to Hamiltonian Monte Carlo_ , _arXiv e-prints_ (2017) arXiv:1701.02434 [1701.02434].
  * [59] M. LoVerde and N. Afshordi, _Extended Limber approximation_ , _Phys. Review D_ 78 (2008) 123506 [0809.5112].
  * [60] L.R. Abramo, P.H. Reimberg and H.S. Xavier, _CMB in a box: Causal structure and the Fourier-Bessel expansion_ , _Phys. Review D_ 82 (2010) 043510 [1005.0563].
  * [61] N. Kaiser, _Clustering in real space and in redshift space_ , _Mon. Not. R. Astron. Soc._ 227 (1987) 1.
  * [62] M. Tegmark, A.N. Taylor and A.F. Heavens, _Karhunen-loeve eigenvalue problems in cosmology: How should we tackle large data sets?_ , _The Astrophysical Journal_ 480 (1997) 22.
  * [63] L.R. Abramo, J.V.D. Ferri and I.L. Tashiro, _Fisher matrix for multiple tracers: the information in the cross-spectra_ , _arXiv e-prints_ (2021) arXiv:2112.01812 [2112.01812].
  * [64] A.J.S. Hamilton, _Towards optimal measurement of power spectra - I. Minimum variance pair weighting and the Fisher matrix_ , _Mon. Not. R. Astron. Soc._ 289 (1997) 285 [astro-ph/9701008].
  * [65] J. Binney and T. Quinn, _Gaussian random fields in spherical coordinates_ , _Mon. Not. R. Astron. Soc._ 249 (1991) 678.
  * [66] P.G. Castro, A.F. Heavens and T.D. Kitching, _Weak lensing analysis in three dimensions_ , _Phys. Review D_ 72 (2005) 023516 [astro-ph/0503479].
  * [67] W. Hu and B. Jain, _Joint galaxy-lensing observables and the dark energy_ , _Phys. Review D_ 70 (2004) 043009 [astro-ph/0312395].
  * [68] T.D. Kitching and A.F. Heavens, _Unequal-time correlators for cosmology_ , _Phys. Review D_ 95 (2017) 063522 [1612.00770].
  * [69] H.S. Grasshorn Gebhardt and D. Jeong, _Nonlinear redshift-space distortions in the harmonic-space galaxy power spectrum_ , _Phys. Review D_ 102 (2020) 083521 [2008.08706].
  * [70] A.J.S. Hamilton, _Uncorrelated modes of the non-linear power spectrum_ , _Mon. Not. R. Astron. Soc._ 312 (2000) 257 [astro-ph/9905191].
  * [71] X. Fang, T. Eifler and E. Krause, _2d-fftlog: efficient computation of real-space covariance matrices for galaxy clustering and weak lensing_ , _Monthly Notices of the Royal Astronomical Society_ 497 (2020) 2699.
  * [72] X. Fang, E. Krause, T. Eifler and N. MacCrann, _Beyond Limber: Efficient computation of angular power spectra for galaxy clustering and weak lensing_ , _Journal of Cosmology and Astroparticle Physics_ 2020 (2020) 010.
  * [73] H.S. Xavier, F.B. Abdalla and B. Joachimi, _Improving lognormal models for cosmological fields_ , _Mon. Not. R. Astron. Soc._ 459 (2016) 3693 [1602.08503].
  * [74] L.C. Maximon, _On the evaluation of the integral over the product of two spherical bessel functions_ , _Journal of Mathematical Physics_ 32 (1991) 642 [https://doi.org/10.1063/1.529405].
  * [75] J.K. Bloomfield, S.H. Face and Z. Moss, _Indefinite integrals of spherical bessel functions_ , _arXiv preprint arXiv:1703.06428_ (2017) .
  * [76] G.N. Watson, _A Treatise on the Theory of Bessel Functions_ , Cambridge University Press, Cambridge, England (1944).

##  Appendix A Appendix: Radial averaging of the spherical Bessel functions

Our aim here is to explore the conditions under which the approximation of Eq.
(3.7) can be employed. First of all, notice that the Fourier integration is of
the form:

$\int d(\ln k)j_{\ell}(kx_{i})j_{\ell}(kx_{j})\,k^{3}P(k)\;,$

and at low redshifts the adimensional spectrum $k^{3}P(k)$ differs
significantly from zero for $k\gtrsim 10^{-2}\,[h$ Mpc${}^{-1}]$ . On the
other hand, the comoving radii that are of most relevance for the large-scale
structures are typically $x\gtrsim 10^{3}\;h^{-1}$ Mpc. Therefore, the
arguments of the spherical Bessel functions are most relevant for $kx\gtrsim
10$.

At this point, recall the asymptotic limits of the spherical Bessel functions:

$\displaystyle\lim_{t\ll\ell}\,j_{\ell}(t)$ $\displaystyle=$
$\displaystyle\frac{t^{\ell}}{(2\ell+1)!!}-\frac{t^{2}}{2}\frac{t^{\ell}}{(2\ell+3)!!}+\ldots\;,$
(A.1) $\displaystyle\lim_{t\gg\ell}\,j_{\ell}(t)$ $\displaystyle=$
$\displaystyle\frac{1}{t}\sin\left(t-\frac{\ell\pi}{2}\right)-\frac{\ell(\ell+1)}{2t^{2}}\cos\left(t-\frac{\ell\pi}{2}\right)+\ldots\;.$
(A.2)

This means that, for $\ell\gtrsim 10$ the Bessel function varies slowly at
first, according to Eq. (A.1), but when $kx$ becomes of the same order as
$\ell$ it starts to oscillate, eventually reaching the regime of fast
oscillations given by Eq. (A.2). For $\ell\lesssim 10$, on the other hand, the
Bessel function rapidly reaches the fast oscillating regime of Eq. (A.2),
which further improves the convergence of the spatial integrals in Eq. (3.5).

In figure 4 we show the accuracy of the approximation of Eq. (3.7), in terms
of the residue ${\rm
Res}=C_{\ell}(\bar{x}_{i},\bar{x}_{j})/\bar{C}_{\ell}(\bar{x}_{i},\bar{x}_{j})-1$,
for two values of the radii, $\bar{x}_{i}=\bar{x}_{j}=10^{3}\;h^{-1}$ Mpc
(left panel) and $\bar{x}_{i}=\bar{x}_{j}=5\times 10^{3}\;h^{-1}$ Mpc (right
panel), and for several radial bin widths between $\Delta x=5$ and $\Delta
x=25\;h^{-1}$ Mpc (for $\bar{x}_{i}\neq\bar{x}_{j}$ the amplitude of the
angular power spectrum falls very fast with $|\bar{x}_{i}-\bar{x}_{j}|$, but
the general trends remain the same). The power spectrum we used for this
exercise was computed using standard cosmological parameters, at $z=0$ –
notice that it is only the shape of the spectrum, not its amplitude, that
matters for the accuracy of the approximation.

Figure 4: Residue of the approximation in Eq. (3.7). The left and right panels
show the fractional errors for $\bar{x}_{i}=\bar{x}_{j}=10^{3}$ and
$\bar{x}_{i}=\bar{x}_{j}=5.\times 10^{3}\;h^{-1}$ Mpc, respectively. The
curves, from top to bottom, correspond to radial bin widths of 5, 10, 15, 20
and 25 $h^{-1}$ Mpc. For this exercise we used a $z=0$ matter power spectrum
with standard cosmological parameters.

It is clear from figure 4 that the approximation is more accurate for small
values of $\ell$ as well as higher radii – although even for the smallest bin
width there is a small bias for all values of $\ell$, which comes from
smearing the Fourier modes inside the radial bins. The accuracy of Eq. (3.7)
becomes worse for progressively smaller radii, but we checked that for
$\bar{x}_{i}=\bar{x}_{j}=500\;h^{-1}$ Mpc, with $\Delta x=5\;h^{-1}$ Mpc the
approximation is accurate to within $\lesssim$5% for $\ell<100$. Conversely,
these results imply that for larger radii we can indeed afford to use wider
radial bins.

##  Appendix B Appendix: Toy model for the covariance and Fisher matrices

###  B.1 Analytical solutions for the angular power spectrum

Consider the Fourier-space data covariance given, in condensed notation, by:

$\Gamma(k)=V_{0}[1+\epsilon f(k)]\;,$ (B.1)

where $\epsilon$ is a small parameter, and $f(k)$ is any adimensional function
of $k$ – in other words, shot noise is given by $V_{0}$ and the power spectrum
is $P(k)=\epsilon V_{0}f(k)$. Assuming for simplicity that we have a single
tracer with mean number density $\bar{n}$, and that its bias is $b=1$, the
harmonic space data covariance $\Gamma_{\ell}(\bar{x},\bar{y})$ is given by:

$\displaystyle\Gamma_{\ell}(\bar{x},\bar{y})$ $\displaystyle=$
$\displaystyle\bar{n}^{2}\int_{\bar{x}}dx\,x^{2}\int_{\bar{y}}dy\,y^{2}\frac{2}{\pi}\int_{0}^{\infty}dk\,k^{2}j_{\ell}(kx)j_{\ell}(ky)V_{0}[1+\epsilon
f(k)]$ (B.2) $\displaystyle\simeq$
$\displaystyle\bar{N}_{\bar{x}}\,\bar{N}_{\bar{y}}\,V_{0}\,\left[\delta_{\bar{x}\,\bar{y}}\frac{4\pi}{\Delta
V_{\bar{x}}}+\epsilon\,\frac{2}{\pi}\int_{0}^{\infty}dk\,k^{2}j_{\ell}(k\bar{x})j_{\ell}(k\bar{y})f(k)\right]\;,$

where recall that $\bar{N}_{\bar{x}}=\bar{n}\,\Delta V_{\bar{x}}/4\pi$ is the
number of tracers per unit solid angle in the radial bin ${\bar{x}}$, and in
the passage to the second line we assumed that the widths of the radial bins
are extremely small.

The inverse of the harmonic covariance is now given by:

$\displaystyle\Gamma_{\ell}^{-1}(\bar{x},\bar{y})$ $\displaystyle=$
$\displaystyle\frac{1}{\bar{n}^{2}}\,\frac{2}{\pi}\int_{0}^{\infty}dk\,k^{2}j_{\ell}(k\bar{x})j_{\ell}(k\bar{y})\frac{1}{V_{0}[1+\epsilon
f(k)]}$ $\displaystyle\simeq$
$\displaystyle\frac{1}{\bar{n}^{2}}\,\frac{1}{V_{0}}\left[\delta_{\bar{x}\,\bar{y}}\frac{4\pi}{\Delta
V_{\bar{x}}}-\,\epsilon\,\frac{2}{\pi}\int_{0}^{\infty}dk\,k^{2}j_{\ell}(k\bar{x})j_{\ell}(k\bar{y})f(k)+{\cal{O}}(\epsilon^{2})\right]\;.$

Any function $f(k)$ that has an exact integral in the expression above could
be used here to make an approximate calculation. E.g., a power-law such as
$f(k)=1/k$ has an exact integral.

There is another interesting toy model that we can build using the following
function:

$\displaystyle 1+\epsilon
f(k)=1+\epsilon\frac{1}{1+\left(\frac{k}{k_{0}}\right)^{n}}=\frac{1+\epsilon+\left(\frac{k}{k_{0}}\right)^{n}}{1+\left(\frac{k}{k_{0}}\right)^{n}}\;.$
(B.4)

Notice that the inverse of the spectrum, which appears in the inverse
covariance, is now basically:

$\displaystyle\frac{1}{1+\epsilon
f(k)}=\frac{1+\left(\frac{k}{k_{0}}\right)^{n}}{1+\epsilon+\left(\frac{k}{k_{0}}\right)^{n}}=\frac{\frac{1}{1+\epsilon}+\frac{1}{1+\epsilon}\left(\frac{k}{k_{0}}\right)^{n}}{1+\frac{1}{1+\epsilon}\left(\frac{k}{k_{0}}\right)^{n}}\;.$
(B.5)

Therefore, defining $1+\epsilon^{\prime}\equiv 1/(1+\epsilon)$, and
$k_{0}^{\prime}\equiv k_{0}(1+\epsilon)^{1/n}$, we have that the Fourier-space
correlation function and its inverse are given by basically the same function:

$\displaystyle\Gamma(k)$ $\displaystyle=$ $\displaystyle
V_{0}\left[1+\epsilon\frac{1}{1+\left(\frac{k}{k_{0}}\right)^{n}}\right]$
(B.6) $\displaystyle\Gamma^{-1}(k)$ $\displaystyle=$ $\displaystyle
V_{0}^{-1}\left[1+\epsilon^{\prime}\frac{1}{1+\left(\frac{k}{k_{0}^{\prime}}\right)^{n}}\right]\;.$
(B.7)

Therefore, we find that the covariance and the inverse covariance have terms
with the same functional form:

$\displaystyle\Gamma_{\ell}(\bar{x},\bar{y})$ $\displaystyle=$
$\displaystyle\bar{N}_{\bar{x}}\bar{N}_{\bar{y}}V_{0}\left[\delta_{\bar{x}\,\bar{y}}\frac{4\pi}{\Delta
V_{\bar{x}}}+\epsilon\,\frac{2}{\pi}\int_{0}^{\infty}dk\,k^{2}j_{\ell}(k\bar{x})j_{\ell}(k\bar{y})\frac{1}{1+(k/k_{0}^{\prime})^{n}}\right]$
(B.8) $\displaystyle\Gamma_{\ell}^{-1}(\bar{x},\bar{y})$ $\displaystyle\simeq$
$\displaystyle\frac{1}{\bar{n}^{2}}\,\frac{1}{V_{0}}\left[\delta_{\bar{x}\,\bar{y}}\frac{4\pi}{\Delta
V_{\bar{x}}}-\,\epsilon^{\prime}\,\frac{2}{\pi}\int_{0}^{\infty}dk\,k^{2}j_{\ell}(k\bar{x})j_{\ell}(k\bar{y})\frac{1}{1+(k/k_{0}^{\prime})^{n}}\right]\;.$
(B.9)

In this way, we only need to integrate a single functional form for the two
expressions.

An interesting limit happens when we take $n\to+\infty$: in that case, we
obtain a step function, and the integral becomes limited to $k\leq k_{0}$ for
both the spectrum and its inverse (since $k_{0}^{\prime}\to k_{0}$ in that
limit). We then have:

$\displaystyle\Gamma_{\ell}(\bar{x},\bar{y})$ $\displaystyle=$
$\displaystyle\bar{N}_{\bar{x}}\bar{N}_{\bar{y}}V_{0}\left[\delta_{\bar{x}\,\bar{y}}\frac{4\pi}{\Delta
V_{\bar{x}}}+\epsilon\,\frac{2}{\pi}\int_{0}^{k_{0}}dk\,k^{2}j_{\ell}(k\bar{x})j_{\ell}(k\bar{y})\right]\;,$
(B.10) $\displaystyle\Gamma_{\ell}^{-1}(\bar{x},\bar{y})$ $\displaystyle=$
$\displaystyle\frac{1}{\bar{n}^{2}}\,\frac{1}{V_{0}}\left[\delta_{\bar{x}\,\bar{y}}\frac{4\pi}{\Delta
V_{\bar{x}}}-\,\epsilon^{\prime}\,\frac{2}{\pi}\int_{0}^{k_{0}}dk\,k^{2}j_{\ell}(k\bar{x})j_{\ell}(k\bar{y})\right]\;,$
(B.11)

where remember that $\epsilon^{\prime}\equiv
1/(1+\epsilon)-1=-\epsilon/(1+\epsilon)$. Now, the integral above has an exact
solution [75]:

$\int_{0}^{1}dz\,z^{2}\,j_{\ell}(zx)j_{\ell}(zy)=\frac{y\,j_{\ell-1}(y)\,j_{\ell}(x)-x\,j_{\ell-1}(x)\,j_{\ell}(y)}{x^{2}-y^{2}}\;\equiv\;g_{\ell}(x,y)\;,$
(B.12)

where we defined:

$\displaystyle g_{\ell}(x,y)$ $\displaystyle\equiv$
$\displaystyle-\frac{x\,j_{\ell-1}(x)j_{\ell}(y)-y\,j_{\ell-1}(y)j_{\ell}(x)}{x^{2}-y^{2}}$
$\displaystyle=$
$\displaystyle\frac{j_{0}(x-y)-(-1)^{\ell}\,j_{0}(x+y)}{2\,x\,y}-\frac{1}{x\,y}\sum_{n=0}^{n_{\ell}}[2(\ell-1-2n)+1]\,j_{\ell-2n}(x)\,j_{\ell-2n}(y)\;.$

The sum in the second line above is over $n=0,1,\ldots,n_{\ell}=(2\ell-1)/4$,
so that $2(\ell-1-2n)+1\geq 0$. The sum stops at $n_{\ell}=(\ell-1)/2$ for odd
values of $\ell$, and at $n_{\ell}=\ell/2$ for even values of $\ell$. This
expression in the second line shows explicitly that this term is well-behaved
as $x\to y$. It is also symmetric under $x\leftrightarrow y$, since
$j_{0}(x-y)=j_{0}(y-x)$ is an even function of its argument.

Therefore, we obtain that:

$\int_{0}^{k_{0}}dk\,k^{2}\,j_{\ell}(kx)j_{\ell}(ky)=k_{0}^{3}\,g_{\ell}(k_{0}x,k_{0}y)\;.$
(B.14)

Substituting this back into Eqs. (B.10) and (B.11) we get:

$\displaystyle\Gamma_{\ell}(\bar{x},\bar{y})$ $\displaystyle=$
$\displaystyle\bar{N}_{\bar{x}}\bar{N}_{\bar{y}}V_{0}\left[\delta_{\bar{x}\,\bar{y}}\frac{4\pi}{\Delta
V_{\bar{x}}}+\,\epsilon\,\frac{2}{\pi}\,k_{0}^{3}\,g_{\ell}(k_{0}\bar{x},k_{0}\bar{y})\right]\;,$
(B.15) $\displaystyle\Gamma_{\ell}^{-1}(\bar{x},\bar{y})$ $\displaystyle=$
$\displaystyle\frac{1}{\bar{n}^{2}}\,\frac{1}{V_{0}}\left[\delta_{\bar{x}\,\bar{y}}\frac{4\pi}{\Delta
V_{\bar{x}}}-\,\frac{\epsilon}{1+\epsilon}\,\frac{2}{\pi}\,k_{0}^{3}\,g_{\ell}(k_{0}\bar{x},k_{0}\bar{y})\right]\;.$
(B.16)

It can be shown, after quite a bit of algebra, that these expressions indeed
satisfy the defining condition for the inverse, Eq. (3.12).

We can in fact go one step further, and combine two "steps" like the one we
used above. Let’s define here the step-functions
$\theta_{1}=\theta_{H}(k-k_{1})$ and $\theta_{2}=\theta_{H}(k-k_{2})$. By
ordering the k’s, $k_{2}>k_{1}$, we get that $\theta_{1}^{2}=\theta_{1}$,
$\theta_{2}^{2}=\theta_{2}$, and $\theta_{1}\theta_{2}=\theta_{2}$. We then
define the Fourier-space correlation function as:

$\Gamma(k)=V_{0}\left(1+\epsilon_{1}\theta_{1}+\epsilon_{2}\theta_{2}\right)\;.$
(B.17)

The inverse of that correlation function is now given by:

$\Gamma^{-1}(k)=V_{0}^{-1}\left[1-\frac{\epsilon_{1}}{1+\epsilon_{1}}\theta_{1}-\frac{\epsilon_{2}}{(1+\epsilon_{1})(1+\epsilon_{1}+\epsilon_{2})}\theta_{2}\right]\;.$
(B.18)

It is easy to see how this expression can be used to generalize Eqs.
(B.15-B.16): it becomes a sum over terms like the ones we have for the case of
the single step-function.

We can in fact keep adding steps to the spectrum and its inverse, and build up
any solution that we want. Using the expressions above we get:

$\displaystyle\Gamma^{\ell}(\bar{x},\bar{y})$ $\displaystyle=$
$\displaystyle\bar{N}_{\bar{x}}\bar{N}_{\bar{y}}V_{0}\left[\delta_{\bar{x}\,\bar{y}}\frac{4\pi}{\Delta
V_{\bar{x}}}+\,\frac{2}{\pi}\,\epsilon_{1}\,k_{1}^{3}\,g_{\ell}(k_{1}\bar{x},k_{1}\bar{y})+\,\frac{2}{\pi}\,\epsilon_{2}\,k_{2}^{3}\,g_{\ell}(k_{2}\bar{x},k_{2}\bar{y})+\ldots\right]\;,$
(B.19) $\displaystyle\Gamma_{\ell}^{-1}(\bar{x},\bar{y})$ $\displaystyle=$
$\displaystyle\frac{1}{\bar{n}^{2}}\,\frac{1}{V_{0}}\left[\delta_{\bar{x}\,\bar{y}}\frac{4\pi}{\Delta
V_{\bar{x}}}-\,\frac{2}{\pi}\,\frac{\epsilon_{1}}{1+\epsilon_{1}}\,k_{1}^{3}\,g_{\ell}(k_{1}\bar{x},k_{1}\bar{y})\right.$
$\displaystyle\quad\quad\quad\left.-\,\frac{2}{\pi}\,\frac{\epsilon_{2}}{(1+\epsilon_{1})(1+\epsilon_{1}+\epsilon_{2})}\,k_{2}^{3}\,g_{\ell}(k_{2}\bar{x},k_{2}\bar{y})+\ldots\right]\;.$

We can make a connection of this last equation with Eq. (B.1), where remember
that $f(k)\to\bar{n}P(k)$ in that equation. We can write it as:

$\displaystyle\Gamma_{\ell}^{-1}(\bar{x},\bar{y})$ $\displaystyle=$
$\displaystyle\frac{2}{\pi}\int_{0}^{\infty}dk\,k^{2}j_{\ell}(k\bar{x})j_{\ell}(k\bar{y})\frac{1}{V_{0}}\left[1-\frac{f(k)}{1+f(k)}\right]$
$\displaystyle=$
$\displaystyle\frac{1}{\bar{n}^{2}}V_{0}^{-1}\left[\delta_{\bar{x}\,\bar{y}}\frac{4\pi}{\Delta
V_{\bar{x}}}-\frac{2}{\pi}\int_{0}^{\infty}dk\,k^{2}j_{\ell}(k\bar{x})j_{\ell}(k\bar{y})\frac{f(k)}{1+f(k)}\right]$
$\displaystyle=$
$\displaystyle\frac{1}{\bar{n}^{2}}V_{0}^{-1}\left[\delta_{\bar{x}\,\bar{y}}\frac{4\pi}{\Delta
V_{\bar{x}}}+\frac{2}{\pi}\int_{0}^{\infty}dk\,k^{3}g_{\ell}(k\bar{x},k\bar{y})\frac{d}{dk}\left(\frac{f}{1+f}\right)\right]\;,$

where the last equality follows from the expression in terms of the “steps”.
So, in some sense we did an integration by parts to arrive at this last
equality.

But let’s go back to Eqs. (B.19-B.1), and consider a “top-hat” function,
$P=\epsilon V_{0}$ for $k_{1}\leq k\leq k_{2}$, and $P=0$ outside of this
interval. This spectral shape can be easily implemented by taking
$\epsilon_{1}\to\epsilon$, $\epsilon_{2}\to-\epsilon$, and we obtain:

$\displaystyle\Gamma^{\ell}(\bar{x},\bar{y})$ $\displaystyle=$
$\displaystyle\bar{N}_{\bar{x}}\bar{N}_{\bar{y}}V_{0}\left\\{\delta_{\bar{x}\,\bar{y}}\frac{4\pi}{\Delta
V_{\bar{x}}}+\,\frac{2}{\pi}\epsilon\left[\,k_{1}^{3}\,g_{\ell}(k_{1}\bar{x},k_{1}\bar{y})-k_{2}^{3}\,g_{\ell}(k_{2}\bar{x},k_{2}\bar{y})\right]\right\\}\;,$
(B.22) $\displaystyle\Gamma_{\ell}^{-1}(\bar{x},\bar{y})$ $\displaystyle=$
$\displaystyle\frac{1}{\bar{n}^{2}}V_{0}^{-1}\left\\{\delta_{\bar{x}\,\bar{y}}\frac{4\pi}{\Delta
V_{\bar{x}}}-\,\frac{2}{\pi}\frac{\epsilon}{1+\epsilon}\,\left[k_{1}^{3}\,g_{\ell}(k_{1}\bar{x},k_{1}\bar{y})-k_{2}^{3}\,g_{\ell}(k_{2}\bar{x},k_{2}\bar{y})\right]\right\\}\;.$
(B.23)

This is, therefore, an exact expression that can be used to compare analytic
and numerical results – and notice that $\epsilon$ does not have to be small.

###  B.2 Comparing the analytical and numerical solutions

The results of the previous section allow us to make an interesting check
using the “top-hat” power Spectrum of Eqs. (B.22-B.23). This will allow us to
compute not only the data covariance, but the covariance matrix for the
angular power spectrum as well, and compare the exact result with a
simulation.

In order to check the validity of our formulas we simulate Gaussian random
fields in a box of $L=500$ h-1 Mpc, with cells of $\Delta L=5$ h-1 Mpc on the
side – i.e., we have a box with $100^{3}$ cells. Our top-hat power spectrum is
given by:

$P_{m}(k)=\epsilon\,[\theta(k-k_{1})-\theta(k-k_{2})]\,,$ (B.24)

with $\epsilon=100$ and:

$\begin{gathered}k_{1}=0.2\,h\,{\rm Mpc}^{-1}\\\ k_{2}=0.3\,h\,{\rm
Mpc}^{-1}\;.\end{gathered}$ (B.25)

In our simulations, for each box we draw Gaussian samples of the Fourier modes
$\tilde{\delta}_{m}(\vec{k})$ of the density field from a probability
distribution function given by the power spectrum:

$P_{m}(|\tilde{\delta}|)\sim
e^{-\frac{1}{2}\frac{|\tilde{\delta}_{m}|^{2}}{P(k)}}.$ (B.26)

Once the Fourier modes have been sampled in the corresponding grid, we take
the inverse Fourier transform in order to obtain the matter density contrast
field $\delta(\vec{r})$. We then assume a mean density of objects in our box
in order to populate each cell with a certain number of tracers, given a bias
for that tracer:

$\begin{gathered}N(\vec{r})\to\bar{N}\,[1+b\,\delta_{m}(\vec{r})]\;\longrightarrow\;\delta(\vec{r})=\frac{N(\vec{r})-\bar{N}}{\bar{N}}\;,\end{gathered}$
(B.27)

where for simplicity we took the tracer bias $b=1$, and used $\bar{N}=10^{3}$
objects per cell. For this exercise we also used a low amplitude for the
spectrum, and sigma-clipping for the very rare instances where cells end up
with $1+b\,\delta_{m}(\vec{r})<0$. The actual number of objects in each cell
is finally obtained by sampling the expression above from a Poisson
distribution.

Once the density field of the tracers $\delta(\vec{x})$ has been determined,
we can compute its Fourier transform and take the expectation value, with the
usual result that:

$\langle\tilde{\delta}(\vec{k})\,\tilde{\delta}^{*}(\vec{k}^{\prime})\rangle=L^{3}\,\delta_{\vec{k}\,\vec{k}^{\prime}}\,\bigg{[}P_{m}(k)+\frac{1}{\bar{n}}\bigg{]}\;.$
(B.28)

Figure 5: Harmonic data covariance $\Gamma_{\ell}(\bar{x}_{i},\bar{x}_{j})$
derived from a top-hat Fourier power spectrum (see text), for $\ell=10$. The
upper panels shows the harmonic data covariance
$\Gamma_{\ell}(\bar{x}_{i},\bar{x}_{j})$. In this example we specified the
number of particles such that shot noise is negligible for
$\bar{x}_{i}\,,\bar{x}_{j}\gtrsim 30\;h^{-1}$ Mpc, i.e., on large scales the
data covariance is identical with the angular power spectrum
$\Gamma_{\ell}(\bar{x}_{i},\bar{x}_{j})\to C_{\ell}(\bar{x}_{i},\bar{x}_{j})$.
The plot on the top left shows a high-resolution image of the harmonic
covariance computed with the help of Eq. (B.31), normalized by its diagonal.
The central panel on the top row shows the result of averaging the high-
resolution covariance in radial bins of $10\;h^{-1}$ Mpc. The upper right
panel shows the normalized covariance in radial bins obtained from the
simulations, also normalized by the diagonal. On the bottom panel we plot the
diagonal of the two matrices. The shaded regions indicate the variance around
the mean: in the case of the simulation (sample mean: blue line), the sample
variance was used, and in the case of the analytical result (orange line), the
expression of Eq. (2.33) was used. The steep rise of the curves as $\bar{x}\to
0$ is due to shot noise, which falls as $\sim\bar{x}^{-2}$.

The next step is to compute the spherical harmonic modes $\delta_{\ell m}$ of
the tracer:

$\delta_{\ell
m}(\bar{x})=\int_{\bar{x}}d^{2}\hat{x}\,\delta(\vec{x})Y^{*}_{\ell
m}(\hat{x})$ (B.29)

where the integral over the solid angle becomes a sum over all cells within
some radial width width, $[\bar{x}-\Delta x/2,\bar{x}+\Delta x/2]$. In the
discretization it is important to enforce the normalization by the effective
solid angle subtended by the cells, $\int_{\bar{x}}d^{2}\hat{x}/4\pi$.

Taking $N_{s}$ samples of the simulated density contrast $\delta_{\ell m}$, we
can compute the tracer angular power spectrum:

$\Gamma_{\ell,\,Sim}(\bar{x},\bar{y})=\dfrac{\bar{N}^{2}}{2\ell+1}\sum_{m=-\ell}^{\ell}\langle\delta_{\ell
m}(\bar{x})\,\delta^{*}_{\ell m}(\bar{y})\rangle_{S}\;,$ (B.30)

where $\langle\ldots\rangle_{S}$ is the sample average.

Our goal here is to compare the angular power spectrum of our simulations
$C_{\ell,\,Sim}(\bar{x},\bar{y})$ to the analytical (or theoretical) result of
Eq. (B.22), which we will denote by $C_{\ell,\,Th}$. In order to draw a fair
comparison, the analytical data covariance is averaged in each shell:

$C_{\ell,\,Th}(\bar{x},\bar{y})=\bar{N}^{2}_{\bar{x}}\,\frac{\delta_{\bar{x},\bar{y}}\,}{\bar{x}^{2}\,\Delta
x}+\,\frac{2}{\pi}\,\epsilon\big{\langle}\left[\,k_{1}^{3}\,g_{\ell}(k_{1}x,k_{1}y)-k_{2}^{3}\,g_{\ell}(k_{2}x,k_{2}y)\right]\big{\rangle}_{\bar{x},\bar{y}}$
(B.31)

where the brackets $\langle\,\ldots\rangle_{\bar{x},\bar{y}}$ denote a mean
over the radial bins.

In order to compare our analytical expressions with simulations we take the
multipole $\ell=10$, and use radial bins of width $\Delta x=10\,h^{-1}\,{\rm
Mpc}$ from $\bar{x}=5\,h^{-1}\,{\rm Mpc}$ to $\bar{x}=245\,h^{-1}\,{\rm Mpc}$.
We produced 1000 boxes following the prescriptions described above, and
computed the sample mean and sample variance of the angular power spectra. The
comparison between the simulations and the theoretical expression is shown in
figure 5: the blue and orange solid lines correspond, respectively, to the
sample mean $\langle C_{\ell,Sim}(\bar{x},\bar{x}{\,}^{\prime})\rangle$ and
the theory data covariance from Eq. (B.31) (orange). The small difference
between the simulations and the theoretical values is concentrated at small
radii, where we indeed expect some smearing of the signal due to poor
sampling. For larger radii the accuracy in the shape of the angular power
spectrum improves, apart from a small bias which also diminishes as we take
smaller radial bins. Moreover, the orange shaded region denotes the sample
variance from the 1000 simulations, while the blue shaded region corresponds
to the covariance of the angular power spectrum calculated from Eq. (2.33),
where we used the $C_{\ell}(\bar{x},\bar{x}{\,}^{\prime})$ from Eq. (B.31).

##  Appendix C Appendix: Some integrals of spherical Bessel functions

Let’s first solve the mixed integral, involving one bare and one second
derivative of a spherical Bessel function, both of the same order $\ell$:

$A=\int dk\,k^{2}\,j_{\ell}(kx)j^{\prime\prime}_{\ell}(ky)\;.$ (C.1)

Using Eq. (4.6) we split this into three integrals:

$\displaystyle A$
$\displaystyle=\int\frac{dk\,k^{2}\,j_{\ell}(kx)}{(k\bar{y})^{2}}\bigg{\\{}\big{[}\ell^{2}-\ell-k^{2}y^{2}\big{]}j_{\ell}(ky)+2ky\,j_{\ell+1}(ky)\bigg{\\}}$
(C.2) $\displaystyle=\frac{(\ell^{2}-\ell)}{y^{2}}\underbrace{\int
dk\,j_{\ell}(kx)j_{\ell}(ky)}_{A_{1}}\,-\underbrace{\int
dk\,k^{2}\,j_{\ell}(kx)j_{\ell}(ky)}_{A_{2}}\,+\,\frac{2}{y}\underbrace{\int
dk\,k\,j_{\ell}(kx)j_{\ell+1}(ky)}_{A_{3}}\;.$ (C.3)

For the first term, we use the following well-known result:

$A_{1}=\int
dk\,j_{\ell}(kx)j_{\ell}(ky)=\frac{\pi}{2}\frac{1}{(2\ell+1)r_{>}}\left(\frac{r_{<}}{r_{>}}\right)^{\ell}\;,$
(C.4)

where $r_{>}$ and $r_{<}$ denote, respectively, the larger and smaller value
between $x$ and $y$. The second term is, once again, simply the closure
relation:

$A_{2}=\int_{0}^{\infty}dk\,k^{2}\,j_{\ell}(kx)j_{\ell}(ky)=\frac{\pi}{2}\frac{\delta_{D}(x-y)}{x^{2}}\;.$
(C.5)

And as for the last integral, we first express the spherical Bessel function
in terms of the usual the Bessel function of the first kind:

$j_{\ell}(z)=\sqrt{\frac{\pi}{2z}}J_{\ell+1/2}(z)\;,$ (C.6)

so we can rewrite that term as:

$A_{3}=\int
dk\,k\,j_{\ell}(kx)j_{\ell+1}(ky)=\frac{\pi}{2}\,\frac{1}{\sqrt{xy}}\int
dk\,J_{\ell+1/2}(kx)\,J_{\ell+3/2}(ky)\;.$ (C.7)

This integral is also known in the literature (see, e.g., [76]), and yields:

$A_{3}=\int
dk\,k\,j_{\ell}(kx)j_{\ell+1}(ky)=\frac{\pi}{2}\times\begin{dcases}\frac{1}{2x^{2}}&\text{if
}x=y\\\\[12.0pt] 0&\text{if }x>y\\\\[12.0pt]
\frac{x^{\ell}}{y^{\ell+2}}&\text{if }x<y\;.\end{dcases}$ (C.8)

Plugging these results back into Eq. (C.2), we obtain:

$\int
dk\,k^{2}j_{\ell}(k{x})j_{\ell}^{\prime\prime}(k{y})=-\frac{\pi}{2{x}^{2}}\delta_{D}({x}-{y})+H_{\ell}({x},{y})\;,$
(C.9)

where:

$H_{\ell}({x},{y})\,=\,\frac{\pi}{2y^{2}}\begin{dcases}\frac{1}{x}\left[\frac{2(2\ell+1)+\ell(\ell-1)}{2(2\ell+1)}\right]&\text{if
}{x}={y}\\\\[12.0pt]
\frac{1}{{x}}\left(\frac{{y}}{{x}}\right)^{\ell}\left[\frac{2(2\ell+1)+\ell(\ell-1)}{(2\ell+1)}\right]&\text{if
}{x}>{y}\\\\[12.0pt]
\frac{1}{{y}}\frac{\ell(\ell-1)}{(2\ell+1)}\left(\frac{{x}}{{y}}\right)^{\ell}&\text{if
}{x}<{y}\;.\end{dcases}$ (C.10)

The next step is to find an analytical expression for the integral which
contains two second derivatives of spherical Bessel functions:

$\int dk\,k^{2}j_{\ell}^{\prime\prime}(k{x})j_{\ell}^{\prime\prime}(k{y})\;.$

In similar fashion to what we did for Eq. (C.1), we rewrite the derivatives
and split this integral into nine parts:

$\begin{gathered}\int\frac{dk}{k^{2}{x}^{2}{y}^{2}}\left\\{\left[\ell^{2}-\ell-k^{2}{y}^{2}\right]\,j_{\ell}(k{y})+2k{y}\,j_{\ell+1}(k{y})\right\\}\left\\{\left[\ell^{2}-\ell-k^{2}{x}^{2}\right]\,j_{\ell}(k{x})+2k{x}\,j_{\ell+1}(k{x})\right\\}\\\
=\frac{(\ell^{2}-\ell)^{2}}{({x}{y})^{2}}\underbrace{\int\frac{dk}{k^{2}}\,j_{\ell}(k{x})j_{\ell}(k{y})}_{B_{1}}\,-\,\frac{(\ell^{2}-\ell)}{{y}^{2}}\underbrace{\int
dk\,j_{\ell}(k{x})j_{\ell}(k{y})}_{B_{2}}\,+\,2\frac{(\ell^{2}-\ell)}{{x}{y}^{2}}\underbrace{\int\frac{dk}{k}j_{\ell+1}(k{x})j_{\ell}(k{y})}_{B_{3}}\\\
\,\,-\frac{(\ell^{2}-\ell)}{{x}^{2}}\underbrace{\int
dk\,j_{\ell}(k{x})j_{\ell}(k{y})}_{B_{4}}\quad+\quad\underbrace{\int
dk\,k^{2}\,j_{\ell}(k{x})j_{\ell}(k{y})}_{B_{5}}\quad-\quad\frac{2}{{x}}\underbrace{\int
dk\,k\,j_{\ell+1}(k{x})j_{\ell}(k{y})}_{B_{6}}\\\
\quad+2\frac{(\ell^{2}-\ell)}{{y}{x}^{2}}\underbrace{\int\frac{dk}{k}j_{\ell}(k{x})j_{\ell+1}(k{y})}_{B_{7}}\,-\,\frac{2}{{y}}\underbrace{\int
dk\,k\,j_{\ell}(k{x})j_{\ell+1}(k{y})}_{B_{8}}\,+\,\frac{4}{{x}{y}}\underbrace{\int
dk\,j_{\ell+1}(k{x})j_{\ell+1}(k{y})}_{B_{9}}\end{gathered}$

Most of the integrals above are known, except for $B_{1}$ and $B_{3}$ (which
is identical to $B_{7}$). Although those two integrals have analytical
formulas in terms of hypergeometric functions [76], we have found more
suitable expressions for these integrals when $\ell$ assumes integer values.
We start by rewriting the integral in $B_{1}$ as:

$B_{1}=r_{>}\int_{0}^{\infty}\frac{dq}{q^{2}}\,j_{\ell}(q)j_{\ell}\left(q\frac{r_{<}}{r_{>}}\right)$
(C.11)

with $q=kr_{>}$, and after many manipulations we find the expression:

$B_{1}=\frac{\pi}{2}\frac{2^{\ell}(\ell+1)!(2\ell-3)!!}{(2\ell+3)!}r_{>}\left(\frac{r_{<}}{r_{>}}\right)^{\ell}\left[2\ell+3-(2\ell-1)\left(\frac{r_{<}}{r_{>}}\right)^{2}\right]\;.$
(C.12)

Similarly, for the integral $B_{3}$ after some algebra we obtain the
analytical expression:

$B_{3}=\frac{\pi}{2}\;\begin{dcases}\frac{1}{(2\ell+1)(2\ell+3)}&\text{if
}{x}={y}\\\\[12.0pt]
\frac{1}{(2\ell+1)(2\ell+3)}\left(\frac{{x}}{{y}}\right)^{\ell+1}&\text{if
}{x}<{y}\;,\\\\[12.0pt]
\frac{1}{2}\left(\frac{{y}}{{x}}\right)^{\ell}\left[\frac{1}{2\ell+1}-\frac{1}{2\ell+3}\left(\frac{{y}}{{x}}\right)^{2}\right]&\text{if
}{x}>{y}\;.\end{dcases}$ (C.13)

We checked these expressions against numerical integrations with the help of
FFTLog [70, 71, 72], from $\ell=2$ up to $\ell=50$, and verified that they are
correct. We were unable to ascertain whether these closed forms of the
integrals $B_{1}$ and $B_{3}$ were previously known.

Collecting these results, together with Eqs. (C.4), (C.5) and (C.8), we can
finally write an analytical expression for the integral: