1. Theory

Authors:

Nicolas Tessore, Lucia F. de la Bella

1.1. Projections in cosmic lensing

Our aim is to compute the angular correlation function \(w(\theta)\) of the cosmic convergence field,

(1.1)\[ w(\theta) = \frac{9H_0^4\Omega_m^2}{4c^4} \iint_{0}^{\infty} \! \frac{q_1(x_1) \, x_1}{a(x_1)} \, \frac{q_2(x_2) \, x_2}{a(x_2)} \, \xi(r_{12}; t_1, t_2) \, dx_1 \, dx_2 \;,\]

where \(H_0\) and \(\Omega_m\) are the cosmological parameters, \(x\) is comoving distance, \(a(x)\) is the scale factor corresponding to \(x\), \(q_1(x)\) and \(q_2(x)\) are the lensing efficiencies given by

(1.2)\[ q_i(x) = \int_{x}^{\infty} \! \frac{x' - x}{x'} \, n_i(x') \, dx'\]

for the distributions \(n_1(x)\) and \(n_2(x)\) of observed sources, and \(\xi(r_{12}; t_1, t_2)\) is the unequal-time matter correlation function for the separation

(1.3)\[r_{12} = \sqrt{x_1^2 + x_2^2 - 2x_1x_2 \cos\theta}\]

at cosmic times \(t_1\) and \(t_2\) corresponding to \(x_1\) and \(x_2\). For simplicity, we have assumed a flat universe in writing (1.1) and (1.2); however, everything that follows applies identically to nonflat geometries.

Other important two-point functions in weak lensing are galaxy clustering and galaxy-galaxy lensing. Galaxy clustering quantifies correlations between galaxy number density fields

(1.4)\[ w_{gg}(\theta) = \frac{1}{c^2} \iint_{0}^{\infty} \! b(x_1) n_1(x_1)H(x_1) \, b(x_2) n_2(x_2)H(x_2) \, \xi(r_{12}; t_1, t_2) \, dx_1 \, dx_2 \;,\]

where \(b(x)\) is the bias parameter and \(H(x)\) the Hubble parameter at a redshift corresponding to a comoving distance \(x\).

Galaxy-galaxy lensing quantifies the correlation between the shape of background (or source) and foreground (lens) galaxy number density. In the weak lensing regime, the observed galaxy shape is the sum of an intrinsic (unlensed component) and a shear due to gravitational lensing. For simplicity, we only consider the shear component:

(1.5)\[ w_{\kappa g}(\theta) = \frac{3H_0^2\Omega_m}{2c^3} \iint_{0}^{\infty} \! \frac{q_1(x_1) \, x_1}{a(x_1)} \, b(x_2) n_2(x_2)H(x_2) \, \xi(r_{12}; t_1, t_2) dx_1 \, dx_2 \;.\]

Generic form:

(1.6)\[ w(\theta) = \iint_{0}^{\infty} \! f_1(x_1) \, f_2(x_2) \, \xi(r_{12}; t_1, t_2) \, dx_1 \, dx_2 \;.\]

Transformation:

(1.7)\[ \xi(r; t_1, t_2) = \frac{1}{2\pi^2} \int_{0}^{\infty} \! P(k; t_1, t_2) \, \frac{\sin kr}{kr} \, k^2 \, dk\]

Addition theorem for spherical Bessel functions ([1], 10.1.45):

(1.8)\[ \frac{\sin kr_{12}}{kr_{12}} = \sum_{l} (2l + 1) \, j_l(kx_1) \, j_l(kx_2) \, P_l(\cos\theta) \;.\]

Inserting (1.7) and (1.8) into (1.6), and exchanging the order of summation and integration, yields the relation

(1.9)\[ w(\theta) = \sum_{l} \frac{2l + 1}{4\pi} \, C_l \, P_l(\cos\theta)\]

between angular correlation function \(w(\theta)\) and the angular power spectrum

(1.10)\[ C_l = \frac{2}{\pi} \iiint_{0}^{\infty} \! f_1(x_1) \, f_2(x_2) \, P(k; t_1, t_2) \, j_l(kx_1) \, j_l(kx_2) \, k^2 \, dk \, dx_1 \, dx_2 \;.\]

In practice, it is not cosmic convergence but cosmic shear that is observable. The two-point statistics are related through their respective angular power spectra \(C_l^{\kappa\kappa}\) and \(C_l^{\gamma\gamma}\), with

(1.11)\[C_l^{\gamma\gamma} = \frac{(l-1) \, (l+2)}{l \, (l+1)} \, C_l^{\kappa\kappa} \;.\]

The results we obtain for cosmic convergence are therefore readily applied to cosmic shear.

1.2. How to compute the projection

Evaluating the angular power spectrum is far from straightforward. The reason is that the integral (1.10) contains the product of two highly oscillatory spherical Bessel functions \(j_l(kx_1) \, j_l(kx_2)\). However, this is not dissimilar to the well-known problem of integrating a spherical function against a pair of highly oscillatory spherical harmonics \(Y_{lm} \, Y_{l'm'}\). Encouragingly, the integral over the sphere is routinely evaluated, by first expanding the function to be integrated into spherical harmonics, and subsequently using Gaunt’s triple integral,

(1.12)\[\int_{S^2} \! Y_{l_1m_1}(\hat{n}) \, Y_{l_2m_2}(\hat{n}) \, Y_{l_3m_3}(\hat{n}) \, d\hat{n} = Y_{l_1l_2l_3m_1m_2m_3} \;,\]

where the Gaunt coefficient \(Y_{l_1l_2l_3m_1m_2m_3}\) is readily computed. We posit that a similar result would represent the most useful analytical solution for the angular power spectrum integral (1.10). It turns out that we can derive the general form of such a solution, if it exists, without performing any actual calculations.

To this end, let us assume for a moment that there exists a set of basis functions \(\tilde{\jmath}_{l'}(k)\) such that i) the unequal-time power spectrum can be expanded in this basis,

(1.13)\[ P(k; t_1, t_2) = \sum_{l'} a_{l'}(t_1, t_2) \, \tilde{\jmath}_{l'}(k) \;,\]

where the modes \(a_{l'}(t_1, t_2)\) of the expansion are necessarily time-dependent, and ii) that the overlap integral between two spherical Bessel functions and the basis functions reduces to coefficients \(\mathfrak{J}_{ll'}\) which can be evaluated,

(1.14)\[ \int_{0}^{\infty} \! j_l(kx_1) \, j_l(kx_2) \, \tilde{\jmath}_{l'}(k) \, k^2 \, dk = \frac{\pi}{2} \, \mathfrak{J}_{ll'}(x_1, x_2) \;.\]

Since \(x_1\) and \(x_2\) appear on the left-hand side of (1.14) as independent variables, the coefficients on the right-hand side must necessarily be functions \(\mathfrak{J}_{ll'}(x_1, x_2)\). Let us finally assume the best-case scenario in which the matrix of functions \(\mathfrak{J}_{ll'}\) is diagonal. Putting together our hypothetical analytical solution (1.13) and (1.14), we find that the angular power spectrum reduces at most to an integral

(1.15)\[ C_l = \iint_{0}^{\infty} \! f_1(x_1) \, f_2(x_2) \, a_l(t_1, t_2) \, \mathfrak{J}_{ll}(x_1, x_2) \, dx_1 \, dx_2 \;,\]

which is precisely the same form as (1.6), the unequal-time angular correlation function!

This short exercise shows that even if there was a convenient analytical solution, similar to the case of the spherical harmonics, the remaining integral (1.15) would still be at least as difficult to compute as the unequal-time angular correlation function.

1.3. Exact projection in real space

The preceding section has shown that the most convenient approach to the exact unequal-time projection of angular correlations is via the real-space integral,

(1.16)\[ w(\theta) = \iint_{0}^{\infty} \! f_1(x_1) \, f_2(x_2) \, \xi(r_{12}; t_1, t_2) \, dx_1 \, dx_2 \;.\]

Evaluating this double integral should, at least in principle, be straightforward. The filter functions \(f_1\) and \(f_2\) are determined by observations, and thus known over a fixed grid of points \(x_1\) and \(x_2\). The filter grid hence provides a natural resolution for numerical integration and, if the remaining factor \(\xi(r_{12}; t_1, t_2)\) can be computed, the value of (1.16) is found by any suitable quadrature scheme, e.g. the trapezoidal rule if the filters are finely enough resolved. In practice, this is exactly how we perform the integration, with one additional complication that we describe shortly.

But first, we compute the matter correlation function \(\xi(r)\) from the matter power spectrum \(P(k)\), which usually is the fundamental input. This is done by expressing the integral relationship (1.7) between the two functions in terms of the Bessel function \(J_{1/2}\),

(1.17)\[ \xi(r; t_1, t_2) = \frac{1}{(2\pi)^{3/2}} \int_{0}^{\infty} \! P(k; t_1, t_2) \, \frac{J_{1/2}(kr)}{\sqrt{kr}} \, k^2 \, dk \;.\]

Integrals of this form can be efficiently evaluated over a logarithmic range of \(r\) values with the FFTLog algorithm [2].

Having obtained the function \(\xi(r; t_1, t_2)\), we need to integrate it against the filter functions \(f_1(x_1)\) and \(f_2(x_2)\). These are usually obtained, either directly or indirectly, from observations, and thus given on a fixed grid of \(x_1\) and \(x_2\) values. Furthermore, the correlations \(\xi(r; t_1, t_2)\) change smoothly with \(r\), \(t_1\), and \(t_2\). Therefore, the numerical integration appears indeed straightforward: Use a fixed-grid quadrature scheme over the \((x_1, x_2)\) grid and interpolate \(\xi(r; t_1, t_2)\) to the points \(r_{12}\), \(t(x_1)\) and \(t(x_2)\) of the filter grid.

In practice, this naive approach sometimes suffers from a catastrophic loss of precision at small scales. The reason is that, although the function \(\xi(r; t_1, t_2)\) does change smoothly with all of \(r\), \(t_1\), and \(t_2\), the resolution of the filter grid \(x_1, x_2\) can nevertheless be insufficient to properly resolve \(\xi(r; t_1, t_2)\) over the entire function domain. To see this, we construct a change of variables from \(x_1, x_2\) to a polar coordinate system where \(r = \sqrt{x_1^2 + x_2^2 - 2x_1x_2 \cos\theta}\) is the radial coordinate. To find the angular coordinate, we only have to write \(r^2\) as a sum of squares; a symmetric choice is

(1.18)\[r^2 = \Bigl[(x_1 - x_2) \, \sqrt{(1 + \cos\theta)/2}\Bigr]^2 + \Bigl[(x_1 + x_2) \, \sqrt{(1 - \cos\theta)/2}\Bigr]^2 \;.\]

We therefore introduce the angle \(\alpha\) as \(\tan(\alpha) = \sqrt{\frac{1 - \cos\theta}{1 + \cos\theta}} \frac{x_1 + x_2}{x_1 - x_2}\). Conversely, the distances \(x_1, x_2\) for given \(r, \alpha\) are

(1.19)\[ x_{1,2} = \sqrt{\frac{1 + \cos(\theta)}{2}} \, \frac{r \sin\alpha}{\sin\theta} \pm \sqrt{\frac{1 - \cos\theta}{2}} \, \frac{r \cos\alpha}{\sin\theta} \;.\]

We can now overlay a regular \((x_1, x_2)\) grid with \((x_1, x_2)\) ellipses of constant \(r\), given by (1.19). The result is shown in Fig. 1.1. We see that the ellipses below a certain scale \(r\) fall entirely between the first grid square, and are hence invisible to the fixed-grid quadrature. Since the correlation function generally shows power law behaviour, these lost contributions to the integral are significant, since they come from where \(\xi(r)\) is largest.

Fig. 1.1 The different grids for the numerical quadrature.

Fortunately, Fig. 1.1 suggests a solution to the problem. To account for every ellipse, i.e. every radius \(r\) at which \(xi(r, t_1, t_2)\) is known, it suffices to go through the \(x_1\) and compute extra points \(x_2\) as the solution for given \(x_1\) and \(r\) (i.e. the intersections of ellipses and vertical lines in Fig. 1.1).

1.4. Limber’s approximation

Write the angular correlation function (1.6) as the integral over the mean radial distance \(x = (x_1 + x_2)/2\) and the radial separation \(R = x_1 - x_2\):

(1.20)\[ w(\theta) = \int_{0}^{\infty} \! \int_{-2x}^{2x} \! f_1(x+R/2) \, f_2(x-R/2) \, \xi(r_{12}; t_1, t_2) \, dR \, dx \;,\]

where the distance between the points in terms of \(x\) and \(r\) is now

(1.21)\[r_{12} = \sqrt{2 x^2 \, (1-\cos\theta) + R^2 \, (1 + \cos\theta)/2} \;.\]

Limber [3][4] introduced an approximation for the integral (1.20) using the assumption i) that the filters and correlation function change slowly and can be approximated by their midpoint values,

(1.22)\[f_1(x+R/2) \, f_2(x-R/2) \, \xi(r_{12}; t_1, t_2) \approx f_1(x) \, f_2(x) \, \xi(r_{12}; t) \;,\]

where \(t = t(x)\); ii) that the angle \(\theta\) between the points is small, \(\theta \ll 1\), so that the distance \(r_{12}\) can be approximated as

(1.23)\[r_{12} \approx \sqrt{x^2\theta^2 + R^2} \;;\]

and iii) that the integral over \(R\) can be extended over the entire real line. Limber’s approximation for the correlation function (1.20) is thus

(1.24)\[ w_{\rm L}(\theta) = \int_{0}^{\infty} \! f_1(x) \, f_2(x) \, \xi_{\rm L}(x\theta; t) \, x\theta \, dx \;,\]

where \(\xi_{\rm L}\) is Limber’s matter correlation function, defined as

(1.25)\[ \xi_{\rm L}(r; t) = \frac{1}{r} \, \int_{-\infty}^{\infty} \! \xi(\sqrt{r^2 + R^2}; t) \, dR \;.\]

If \(\xi_{\rm L}\) is known, the angular correlation function (1.24) is a single integral over a slowly changing combined filter function \(f_{12}(x) = f_1(x) \, f_2(x)\). If the filter is determined by observations, the resolution of the integral fixed, and it can be evaluated by any suitable method.

It remains to compute \(\xi_{\rm L}\) from the integral (1.25). Using the relation (1.7) to express the matter correlation function \(\xi\) as an integral over the matter power spectrum and exchanging the order of the integrals yields a representation of the Bessel function \(J_0\),

(1.26)\[\frac{1}{\pi} \int_{-\infty}^{\infty} \! \frac{\sin(\sqrt{x^2 + y^2})}{\sqrt{x^2 + y^2}} \, dy = J_0(x) \;.\]

We thus obtain an expression for Limber’s matter correlation function in terms of the matter power spectrum,

(1.27)\[ \xi_{\rm L}(r; t) = \frac{1}{2\pi} \int_{0}^{\infty} \! P(k; t) \, \frac{J_0(kr)}{kr} \, k^2 \, dk \;.\]

This is of similar form as the unequal-time matter correlation function (1.17), and we can express both in the generic form

(1.28)\[ \xi_{\mu}(r; \ldots) = \frac{1}{(2\pi)^{1+\mu}} \int_{0}^{\infty} \! P(k; \ldots) \, \frac{J_\mu(kr)}{(kr)^{1-\mu}} \, k^2 \, dk \;,\]

setting \(\mu = 0\) in the Limber case, and \(\mu = 1/2\) in the exact case. In practice, this allows us to use a single generic implementation of the FFTLog algorithm to compute either the unequal-time matter correlation function (1.17) or Limber’s matter correlation function (1.27).

1.5. Angular power spectrum

The angular correlation function \(w(\theta)\) of a scalar field is related to its angular power spectrum \(C_l\) by the sum (1.9). The inverse relation is the integral [5]

(1.29)\[C_l = 2\pi \int_{0}^{\pi} \! w(\theta) \, P_l(\cos\theta) \sin(\theta) \, d\theta \;.\]

Instead of a numerical integration of the angular power spectrum, we perform the following procedure.

Given a set of angles \(\theta_1, \theta_2, \ldots\), the computed angular correlation function forms the vector \(w = (w_k)\) with components \(w_k = w(\theta_k)\). Let \(A = (a_{kl})\) be the matrix with entries \(a_{kl} = (2l + 1)/(4\pi) \, P_l(\cos\theta_k)\) up to some maximum number \(l_{\max}\). The truncated sum (1.9) can hence be written

(1.30)\[ w_k = \sum_{l=0}^{l_{\max}} a_{kl} \, C_l\]

or, in matrix form, \(w = Ac\), if \(c = (C_l)\) is the vector of angular power spectrum entries.

Hence, to obtain \(C_l\) for \(l \le l_{\max}\) from \(w(\theta)\), we compute sufficiently many values \(w_k\), and use a least squares solution of the matrix equation.

1.6. References

1

M. Abramowitz and I. A. Stegun. Handbook of Mathematical Functions. National Bureau of Standards, 1972.

2

A. J. S. Hamilton. Uncorrelated modes of the non-linear power spectrum. MNRAS, 312(2):257–284, February 2000. arXiv:astro-ph/9905191, doi:10.1046/j.1365-8711.2000.03071.x.

3

D. Nelson Limber. The Analysis of Counts of the Extragalactic Nebulae in Terms of a Fluctuating Density Field. ApJ, 117:134, January 1953. doi:10.1086/145672.

4

D. Nelson Limber. The Analysis of Counts of the Extragalactic Nebulae in Terms of a Fluctuating Density Field. II. ApJ, 119:655, May 1954. doi:10.1086/145870.

5

Nicolas Tessore. The Spectral Representation of Homogeneous Spin-Weighted Random Fields on the Sphere. arXiv e-prints, pages arXiv:1904.09973, April 2019. arXiv:1904.09973.