3D cosmic shear: Numerical challenges, 3D lensing random fields generation, and Minkowski functionals for cosmological inference

A. Spurio Mancini*

Institut für Theoretische Physik, Universität Heidelberg, Philosophenweg 12, 69120 Heidelberg, Germany

P. L. Taylor

Mullard Space Science Laboratory, University College London,
Holmbury St. Mary, Dorking, Surrey RH5 6NT, United Kingdom

R. Reischke

Astronomisches Rechen-Institut, Zentrum für Astronomie der Universität Heidelberg,
Philosophenweg 12, 69120 Heidelberg, Germany,
Institut für Kernphysik, Karlsruher Institut für Technologie, 76344 Eggenstein-Leopoldshafen, Germany

T. Kitching

Mullard Space Science Laboratory, University College London,
Holmbury St. Mary, Dorking, Surrey RH5 6NT, United Kingdom

V. Pettorino

AIM, CEA, CNRS, Université Paris-Saclay, Université Paris Diderot,
Sorbonne Paris Cité, F-91191 Gif-sur-Yvette, France

B. M. Schäfer

Astronomisches Rechen-Institut, Zentrum für Astronomie der Universität Heidelberg,
Philosophenweg 12, 69120 Heidelberg, Germany

B. Zieser and Ph. M. Merkel

Institut für Theoretische Astrophysik, Zentrum für Astronomie der Universität Heidelberg,
Philosophenweg 12, 69120 Heidelberg, Germany

<figure><img src="image_1.png" /></figure>

(Received 21 August 2018; published 8 November 2018)

Cosmic shear—the weak gravitational lensing effect generated by fluctuations of the gravitational tidal fields of the large-scale structure—is one of the most promising tools for current and future cosmological analyses. The spherical-Bessel decomposition of the cosmic shear field (3D cosmic shear) is one way to maximize the amount of redshift information in a lensing analysis and therefore provides a powerful tool to investigate in particular the growth of cosmic structure that is crucial for dark energy studies. However, the computation of simulated 3D cosmic shear covariance matrices presents numerical difficulties, due to the required integrations over highly oscillatory functions. We present and compare two numerical methods and relative implementations to perform these integrations. We then show how to generate 3D Gaussian random fields on the sky in spherical coordinates, starting from the 3D cosmic shear covariances. To validate our field-generation procedure, we calculate the Minkowski functionals associated with our random fields, compare them with the known expectation values for the Gaussian case and demonstrate parameter inference from Minkowski functionals from a cosmic shear survey. This is a first step towards producing fully 3D Minkowski functionals for a lognormal field in 3D to extract Gaussian and non-Gaussian information from the cosmic shear field, as well as towards the use of Minkowski functionals as a probe of cosmology beyond the commonly used two-point statistics.

DOI: 10.1103/PhysRevD.98.103507

I. INTRODUCTION

The weak gravitational lensing effect is the distortion of
light bundles due to the differential deflection of light rays

<footnote>*Corresponding author.
spuriomancini@thphys.uni-heidelberg.de</footnote>
coming from distant sources [see e.g., [1–3], for reviews on the topic] caused by the gravitational tidal fields generated by the mass distribution between the source and the observer. This effect induces a small change in the ellipticity of a background galaxy’s image, which is referred to as ‘cosmic shear’ when it is caused by the gravitational fields of the large-scale structure of the Universe. Cosmic shear measurements are of a statistical nature, as the lensing effect is not associated with a particular intervening lens, but rather corresponds to small distortions (of the order of 1%) by all potential fluctuations along the line of sight; detecting the extremely faint cosmic shear signal requires averaging over many background galaxies. The statistical properties of the shear field reflect those of the underlying density field by virtue of the gravitational field equations. The cosmic shear field has zero mean; at the level of one-point statistics, cosmological information can be extracted from e.g., peak counts [4–7], while for two-point statistics one looks in configuration space at the angular correlation function of the shear field, or its equivalent in Fourier space, the cosmic shear angular power spectrum. At higher order, cosmic shear can also break degeneracies between the dark sector and neutrinos [8].

Since the first detections in early 2000s [e.g., [9–11]], cosmic shear has flourished into a well-established theoretical framework and important cosmological constraints have already been derived from cosmic shear analyses of different surveys over the last decade, such as the Canada-France-Hawaii Telescope Lensing Survey [CFHTLenS; [12–14]] and the Kilo Degree Survey [KiDS; [15–17]]. Even stronger constraints are expected in the next decade from planned Stage IV surveys such as Euclid<fnref n="1" />[18] and the Large Synoptic Survey Telescope<fnref n="2" />[LSST; [19]], which will use the weak gravitational lensing effect as one of the main cosmological probes to achieve their science goals, most importantly the investigation of the accelerated expansion of the Universe through its influence on the background expansion and cosmic structure growth. Studies to ascertain the nature of the dark energy component of the Universe, to which the acceleration can be ascribed, crucially depend on the sensitivity of the cosmological probes to the growth of cosmic structures. Cosmic shear appears therefore very promising for this purpose, because it is particularly sensitive to the growth.

Extracting information along the radial direction becomes then a necessary requirement for any cosmic shear analysis that aims at increasing the sensitivity to the growth of cosmic structure. For this reason, a standard approach [“tomography,” first introduced in [20]] for analyzing a cosmic shear survey consists in calculating correlations of the lensing signal between different redshift bins, to which the observed galaxies are assigned according

to their estimated (photometric) redshifts, in an attempt to reduce loss of information due to radial averaging.

Alternatively to this tomographic approach, a spherical-Bessel decomposition of the cosmic shear field [“3D cosmic shear”, first introduced in [21]] allows for the inclusion in the analysis of the redshift information on each galaxy within the survey. This fully 3D approach has been recently studied in [22,23] in the context of modified gravity theories. In [22,24], the cosmological sensitivity of the 3D and tomographic approach have been compared. The code used to produce the results in [24] has been released in [25],<fnref n="3" /> where a study of the flat universe approximation has been performed in view of future Stage IV surveys. [24] also discussed alternative weights to the spherical-Bessel ones, which may be easier to compute while still preserving most of the radial information.

The 3D spherical-Bessel formalism presents challenging integrals to evaluate numerically: we highlight them while reviewing the 3D cosmic shear formalism in Sec. II. Subsequently, we describe two different numerical techniques used to evaluate those integrals, namely the numerical recipes underlying the results presented in [22] and [24,25]: we present them in Secs. III A and III B, respectively. In Sec. IV, we present a comparison of 3D cosmic shear covariance matrices obtained with the two numerical methods. Having two completely independent numerical techniques to tackle the 3D cosmic shear integrations, producing results in excellent agreement between them, is useful for a number of future applications that employ the simulated 3D covariance matrices. These include e.g., a cross-correlation analysis of 3D cosmic shear and galaxy clustering [see [26], for a spherical-Bessel analysis of a spectroscopic galaxy clustering survey], or the development of a Bayesian Hierarchical Model for 3D cosmic shear power spectra estimation [see [27,28], for a Bayesian Hierarchical Model for tomography]. We show in Sec. V how to make use of the 3D covariance matrices to generate Gaussian random fields on the sky, showing in particular how to overcome the difficulties arising from the non-diagonality of the covariance matrices in the radial coordinate, which originates from the inhomogeneity of the lensing field along the line of sight. We test the validity of our field-generation procedure in Sec. VI, where we briefly describe and then calculate the Minkowski functionals of our generated Gaussian random fields, and compare them with their expectation values, known analytically in the Gaussian case. We also demonstrate in a first simplified case how a likelihood analysis for cosmological inference can be carried out using the estimated Minkowski functionals. Finally, we draw our conclusions in Sec. VII.

<footnote n="1">¹https://www.euclid-ec.org/.</footnote>

<footnote n="2">²https://www.lsst.org/.</footnote>

<footnote n="3">³The code is available at https://github.com/astro-informatics/GLaSS. For questions on the code, please contact P. Taylor at peterllewelyntaylor@gmail.com.</footnote>
## II. 3D COSMIC SHEAR

We begin by reviewing the formalism for a fully 3D expansion of the shear field based on its spherical-Bessel decomposition, as first introduced in lensing studies by Heavens [21]. Here we follow the notation and conventions of Spurio Mancini et al. [22] and Zieser and Merkel [29] [see also [24,25,30–34]].

Information on the gravitational lensing effect is encoded in the lensing potential, a weighted projection of the gravitational potential along the line of sight. In a standard general relativity context [see [22,23], for an extended formalism valid also for modified gravity theories], the lensing potential $\phi$ is related to the gravitational potential $\Phi$ by

$$ \phi(\chi, \hat{\mathbf{n}}) = \frac{2}{c^2} \int_0^\chi d\chi' \frac{\chi - \chi'}{\chi \chi'} \Phi(\chi, \hat{\mathbf{n}}), \quad (1) $$

where $\chi$ is a comoving distance, and the normalized vector $\hat{\mathbf{n}}$ selects a direction on the sky. Here and throughout the paper spatial flatness will be assumed [for expressions for a nonflat Universe, see [25]], and the integration in Eq. (1) is carried out in Born approximation, i.e., along the unperturbed light path. The shear tensor $\gamma(\chi, \hat{\mathbf{n}})$ is defined as the second $\partial\partial'$-derivative [35,36] of the lensing potential [21,37]

$$ \gamma(\chi, \hat{\mathbf{n}}) = \frac{1}{2} \partial\partial\phi(\chi, \hat{\mathbf{n}}). \quad (2) $$

The $\partial$-derivative acts as a covariant differentiation operator on the celestial sphere and relates quantities of different spin, raising the spin $s$ of a field, a number which characterizes its transformation properties under rotations. Acting twice on $\phi$, the $\partial$ operator relates the scalar (spin-0) lensing potential $\phi$ to the spin-2 shear field $\gamma$. The shear $\gamma$ can be expanded with a choice of basis functions given by a combination of spherical Bessel functions $j_\ell(z)$ [38] and spin 2-weighted spherical harmonics $_2Y_{\ell m}(\hat{\mathbf{n}})$

$$ \gamma(\chi, \hat{\mathbf{n}}) = \sqrt{\frac{2}{\pi}} \sum_{\ell m} \int k^2 dk \gamma_{\ell m}(k) {}_2Y_{\ell m}(\hat{\mathbf{n}}) j_{\ell}(k\chi), \quad (3) $$

where the coefficients $\gamma_{\ell m}(k)$ are given by

$$ \gamma_{\ell m}(k) = \sqrt{\frac{2}{\pi}} \int \chi^2 d\chi \int d\Omega \gamma(\chi, \hat{\mathbf{n}}) j_{\ell}(k\chi) {}_2Y_{\ell m}^*(\hat{\mathbf{n}}). \quad (4) $$

Inserting Eqs. (1) and (2) in (4), and applying a spherical-Bessel expansion to the gravitational potential $\Phi$, we can rewrite $\gamma$ as

$$ \gamma(\chi, \hat{\mathbf{n}}) = \sqrt{\frac{2}{\pi}} \frac{1}{c^2} \int_0^\chi d\chi' \frac{\chi - \chi'}{\chi \chi'} \times \int k^2 dk \sum_{\ell m} \sqrt{\frac{(\ell + 2)!}{(\ell - 2)!}} \times \Phi_{\ell m}(k, \chi') j_\ell(k\chi') {}_2Y_{\ell m}(\hat{\mathbf{n}}). \quad (5) $$

Poisson's equation can be used to link the coefficients in the spherical-Bessel decomposition of the lensing potential to those of the density contrast field $\delta_{\ell m}(k, \chi)$,

$$ \frac{\Phi_{\ell m}(k, \chi)}{c^2} = -\frac{3}{2} \frac{\Omega_m}{(k\chi_H)^2} \frac{\delta_{\ell m}(k, \chi)}{a(\chi)}, \quad (6) $$

with the Hubble radius $\chi_H \equiv c/H_0$. If fluctuations in $\Phi$ are weak, the density field is statistically homogeneous and isotropic, characterized by a power spectrum $P_\delta(k, z, z')$ which is diagonal in harmonic space

$$ \langle \delta_{lm}(k, z) \delta_{\ell' m'}^*(k', z') \rangle = \frac{P_\delta(k, z, z')}{k^2} \delta^D(k-k') \delta_{\ell\ell'}^K \delta_{mm'}^K. \quad (7) $$

Using this, we can relate the covariance of shear modes to the matter power spectrum by

$$ \langle \bar{\gamma}_{lm}(k) \bar{\gamma}_{\ell'm'}^*(k') \rangle = \frac{9 \Omega_m^2 (\ell + 2)!}{16 \pi^4 \chi_H^4 (\ell - 2)!} \times \int \frac{d\tilde{k}}{\tilde{k}^2} G_\ell(k, \tilde{k}) G_{\ell'}(k', \tilde{k}) \delta_{\ell\ell'}^K \delta_{mm'}^K, \quad (8) $$

where

$$ G_{\ell}(k, k') = \int dz n_z(z) F_{\ell}(z, k) U_{\ell}(z, k'), \quad (9) $$

$$ F_{\ell}(z, k) = \int dz_p p(z_p|z) j_{\ell}[k\chi^0(z_p)], \quad (10) $$

$$ U_{\ell}(z, k) = \int_{0}^{\chi(z)} \frac{d\chi'}{a(\chi')} \frac{\chi - \chi'}{\chi \chi'} j_{\ell}(k\chi') \sqrt{P_{\delta}(k, z(\chi))}. \quad (11) $$

$\bar{\gamma}$ are estimates of the shear modes that, in addition to the pure lensing effect, keep into account the redshift distribution of the lensed galaxies $n_z(z)$ and the conditional probability $p(z_p|z)$ of estimating the redshift $z_p$ given the true redshift $z$. The contribution to the total signal coming from sources situated at different distances is governed by the source density $n_z(z)$; through this term, the survey depth affects the strength of the overall signal. Angular variations are assumed to be negligible by considering a uniform source density: the number of sources per steradian and redshift interval is approximated by the mean $n_z(z)/(4\pi)$ across the sky. The influence of incomplete sky
coverage is ignored in this formalism: for applications to a Fisher matrix analysis, for example, the effect of partial sky coverage can be well compensated by a multiplying factor $f_{\text{sky}}$ denoting the fraction of sky spanned by the survey, prepended to the expression for the Fisher matrix [22,32]. A finite field of view can be incorporated in the analysis by considering a suitable window function $W(\hat{n})$ that represents the angular distribution of the sources (e.g., a top hat filter corresponding to a rectangular field of view). For details on the extended formalism to include in the analysis of this inhomogeneous sampling, we refer the reader to e.g., Heavens [21], while e.g., Leistedt et al. [39] consider alternative methods, such as wavelets, to deal with survey geometry in 3D cosmic shear.

Statistical isotropy guarantees that the covariance in Eq. (8) does not depend on the multipole order $m$, while the assumed full sky coverage also prevents mixing of different $\ell$ modes. If the finite field of view is taken into account, statistical isotropy is broken (e.g., by the absence of data across parts of the sky) leading to a coupling of different $\ell$-modes; furthermore, if the field of view is not square, even for a fixed $\ell$ there will be different results for different $m$-modes. The lensing weight function, the redshift errors and the redshift-dependence of the source distribution, instead, always introduce correlations between the amplitudes of the signal on different scales; the covariance matrix then acquires off-diagonal terms, the calculation of which is numerically involved [see [14], $C_{\ell}$ approach in 3D].

The basis of spherical Bessel functions leads to integrals with rapidly oscillatory kernels, which in the inference process have to be solved for a large number of parameter combinations. The $\sqrt{P_{\delta}(k, z(\chi))}$ term comes from an approximation, introduced and qualitatively justified in Castro et al. [37], to calculate unequal-time correlators appearing in the comoving distance integrations by means of a geometric mean $P(k, z, z') \simeq \sqrt{P(k, z)P(k, z')}$ [see also [40]]. This expression simplifies considerably in the linear regime of structure formation, retrieving the one presented in the seminal paper of Heavens [21] where a product of the linear growth factors at different redshifts is present, acting on the matter power spectrum evaluated at the present time.

The noise term for the covariance matrix of the shear modes is given by the intrinsic ellipticity dispersion of source galaxies, as a result of the fact that the observed ellipticity $\epsilon$ is assumed to be the sum of the shear $\gamma$ and the intrinsic ellipticity $\epsilon_S$. The intrinsic ellipticity dispersion is given by $\langle \epsilon_S^2 \rangle = \sigma_{\epsilon}^2$. As given in Heavens et al. [32] and Kitching et al. [33], in the spherical-Bessel formalism (see the Appendix for a complete derivation), this leads to

$$ \langle \gamma_{\ell m}(k) \gamma_{\ell' m'}(k') \rangle_{\text{SN}} = \frac{\sigma_{\epsilon}^2}{2\pi^2} \int dz n_z(z) j_{\ell}(k\chi_0(z)) j_{\ell'}[k'\chi_0(z)] \delta_{\ell\ell'}^K \delta_{mm'}^K, \quad (12) $$

where the subscript SN stands for “shot noise” and we set $\sigma_{\epsilon} = 0.3$. This expression for the noise holds only in absence of intrinsic alignments, i.e., assuming that the intrinsic ellipticities of galaxies are uncorrelated [see [34,41], for a study of intrinsic alignments in 3D cosmic shear].

### III. NUMERICAL IMPLEMENTATION

In this section, we will briefly describe the two methods used to calculate the correlations from Eqs. (8) and (12). While one code implements in C `++` the Levin collocation method [42,43] that makes use of the periodic oscillations of the Bessel functions and has been used to produce the results of Zieser and Merkel [29] and Spurio Mancini et al. [22], the other implements the integrations by matrix multiplications and appropriate use of the Limber approximation [44–46] at high $\ell$. The second code is a Python module, implemented in the code GLASS and used in Taylor et al. [24,25].

#### A. Levin integration

The method presented in [42,43] can be used for efficient evaluation of rapidly oscillatory integrals, once certain conditions are satisfied. The main idea behind the method is to transform the quadrature problem into the solution of a system of linear ordinary differential equations. These are then tackled by collocation, i.e., choosing candidate solutions (polynomials) and a number of points in the domain (called collocation points), and selecting that solution which satisfies the given equations at the collocation points.

As seen in Sec. II the 3D cosmic shear signal and noise [cf. Eqs. (8) and (12)] present several integrals of the form

$$ I_1[h] = \int_{z_1}^{z_2} dz h(z) j_{\ell}(k\chi(z)) \quad (13) $$

or

$$ I_2[h] = \int_{z_1}^{z_2} dz h(z) j_{\ell}(k_1\chi(z)) j_{\ell}(k_2\chi(z)). \quad (14) $$

The comoving distance between two events at redshift $z_1$ and $z_2$ is given by

$$ \chi(z_1, z_2) = \int_{a(z_2)}^{a(z_1)} \frac{cda}{a^2 H(a)} = \chi_H \int_{z_1}^{z_2} \frac{dz}{E(z)}, \quad (15) $$

where $a$ is the scale factor and $\dot{a}/a = H(a) = H_0 E(a)$ is the Hubble function. Rather than redshift integrals, Eqs. (13) and (14) can be rewritten, using $dz = d\chi E[z(\chi)]$, with the comoving distance as the integration variable:

$$ I_1[h] = \frac{1}{\chi_H} \int_{\chi(z_1)}^{\chi(z_2)} d\chi E[z(\chi)] h[z(\chi)] j_{\ell}(k\chi), \quad (16) $$
$$I_2[h] = \frac{1}{\chi_H} \int_{\chi(z_1)}^{\chi(z_2)} d\chi E[z(\chi)] h[z(\chi)] j_{\ell}(k_1\chi) j_{\ell}(k_2\chi). \quad (17)$$

Due to the highly oscillatory nature of the spherical Bessel functions, especially at high $\ell$ or $k$, the numerical solution of these integrals by standard quadrature routines is extremely inaccurate when a large number of zero-crossings occurs in the interval $[\chi(z_1), \chi(z_2)]$, unless an enormous number of points is used to sample the integrand; however, the procedure becomes then exceedingly time consuming, especially if many combinations of $\ell$ and $k$ need to be considered.

Here we describe an alternative method, presented by [42,43], which we use to evaluate our integrals. It is applicable to integrals of the form

$$I[F] = \int_a^b dx F^T(x) w(x) = \int_a^b dx \langle F, w \rangle(x), \quad (18)$$

where $F(x) = [F_1(x), \dots, F_d(x)]^T$ and $w(x) = [w_1(x), \dots, w_d(x)]^T$ are vectors of functions, for which the second equality of Eq. (18) defines a scalar product $\langle \cdot, \cdot \rangle$ and the functions $w_i(x), i = 1, 2, \dots, d$, but not $F_i(x)$, are rapidly oscillatory across the integration domain. A matrix of functions $A(x)$ is defined, such that the derivatives of $w(x)$, denoted by $w'(x)$, fulfil

$$w'(x) = A(x)w(x). \quad (19)$$

The components $A_{iq}(x)$ should not be highly oscillatory. We show below an example of such a matrix for the particular cases given in Eqs. (16) and (17). In the Levin formalism, a vector $p(x)$ is constructed to approximate the integrand in Eq. (18) by

$$\langle p, w \rangle' = \langle p' + A^T p, w \rangle \approx \langle F, w \rangle. \quad (20)$$

The first equality follows from applying the Leibniz rule for derivatives and Eq. (19), with $\langle p, Aw \rangle = \langle A^T p, w \rangle$. If such a vector is found, then the integral in Eq. (18) can be approximated by

$$I[F] \approx \int_a^b dx \langle p, w \rangle'(x) = \langle p, w \rangle(b) - \langle p, w \rangle(a). \quad (21)$$

This can be achieved by demanding that both terms should be equal, $\langle p, w \rangle' = \langle F, w \rangle$, at n collocation points $x_j, j = 1, 2, \dots, n$. The requirement

$$\langle p' + A^T p - F, w \rangle(x_j) = 0, \quad j = 1, \dots, n \quad (22)$$

generally means that the vector $\langle p' + A^T p - F \rangle$ must be orthogonal to $w$ at the points $x_j$, for example by demanding that it should be the null vector:

$$p'(x_j) + A^T(x_j) p(x_j) = F(x_j). \quad (23)$$

Finding a vector $p$ which has this property can be achieved by choosing a set of $n$ linearly independent and differentiable basis functions $u_m(x)$ and writing each component $p_i(x)$ as a linear combination:

$$p_i(x) = c_i^{(m)} u_m(x), \quad i = 1, \dots, d; \quad m = 1, \dots, n. \quad (24)$$

Equation (23) then leads to the following linear system of equations for the $d \times n$ coefficients $c_i^{(m)}$:

$$c_i^{(m)} u'_m(x_j) + A_{qi} c_q^{(m)} u_m(x_j) = F_i(x_j), \quad i, q = 1, \dots, d; \quad j, m = 1, \dots, n. \quad (25)$$

Levin [42] showed how to concretely apply this algorithm to several cases of integrals with highly oscillatory kernels. The performance varies depending on the integrand, but accuracies below $10^{-6}$ can often be achieved with less than 10 collocation points. As suggested by Levin [42], in our implementation we use equidistant collocation points

$$x_j = a + (j - 1) \frac{b - a}{n - 1}, \quad j = 1, \dots, n \quad (26)$$

and choose the $n$ lowest-order polynomials as basis functions:

$$u_m(x) = \left( \frac{x - \frac{a+b}{2}}{b-a} \right)^{m-1}, \quad m = 1, \dots, n. \quad (27)$$

We note that the polynomials with $m > 1$ and the derivatives with $m > 2$ share the root $x = (a + b)/2$: to prevent the linear system of equations from becoming singular, that root should not be used as a collocation point. The factor $1/(b-a)$ is included for numerical reasons: if $b \gg 1$ or $b \ll 1$, the values of polynomials of different order may differ by several orders of magnitude; the normalizing factor guarantees that $|u_m(x)| \le 1$ across the integration domain, in order to regulate the range of the coefficients of the linear system of equations in Eq. (25) and thus the condition of the corresponding matrix. Suitable vectors $w$ for the integrals in Eqs. (16) and (17) can be identified by considering the following recurrence relations for the spherical Bessel functions [38]:

$$\frac{d}{dx} j_{\ell}(x) = j_{\ell-1}(x) - \frac{\ell+1}{x} j_{\ell}(x), \quad (28)$$

$$\frac{d}{dx} j_{\ell-1}(x) = -j_{\ell}(x) + \frac{\ell-1}{x} j_{\ell-1}(x). \quad (29)$$

Rewriting these relations in the form $w' = Aw$, one finds that
$$w(\chi) = \begin{pmatrix} j_{\ell}(k\chi) \\ j_{\ell-1}(k\chi) \end{pmatrix}, \quad A(\chi) = \begin{pmatrix} -\frac{\ell+1}{\chi} & k \\ -k & \frac{\ell-1}{\chi} \end{pmatrix} \quad (30)$$

is a suitable choice for the integral in Eq. (16), with $F(\chi) = \{E[z(\chi)]h[z(\chi)], 0\}^T$. It is easy to verify that neither the entries of the matrix A nor the integral kernels F are rapidly oscillatory. For integrals of the type in Eq. (17), four-dimensional vectors are needed:

$$w(\chi) = \begin{pmatrix} j_{\ell}(k_1\chi)j_{\ell}(k_2\chi) \\ j_{\ell-1}(k_1\chi)j_{\ell}(k_2\chi) \\ j_{\ell}(k_1\chi)j_{\ell-1}(k_2\chi) \\ j_{\ell-1}(k_1\chi)j_{\ell-1}(k_2\chi) \end{pmatrix}, \\ A(\chi) = \begin{pmatrix} -\frac{2(l+1)}{\chi} & k_1 & k_2 & 0 \\ -k_1 & -\frac{2}{\chi} & 0 & k_2 \\ -k_2 & 0 & -\frac{3}{\chi} & k_1 \\ 0 & -k_2 & -k_1 & \frac{2(\ell-1)}{\chi} \end{pmatrix}. \quad (31)$$

Similarly, $F(\chi) = \{E[z(\chi)]h[z(\chi)], 0, 0, 0\}^T$.

## B. GLaSS

The Generalised Lensing and Shear Spectra (GLaSS) code is written in Python and integrated into the modular cosmological package Cosmosis [47]. Cosmological information can be read from an external source as in this work, or directly from the Cosmosis pipeline. More information can be found in [25].

GLaSS is written to compute the lensing spectra for an arbitrary weight function $W_{\ell}[k\chi^0(z_p)]$ which takes the place of the Bessel functions in Eq. (10); see [24] for more details about this generalized spherical-transform. Nevertheless, 3D cosmic shear comes as an in-built run-mode option.

All nested integrals in Eqs. (8)–(11) are computed as matrix multiplications because this is one of the few operations that releases the global interpreter lock in Python allowing parallelization. For example,

$$U_{\ell}(z, k) \approx \sum_{\chi'} A(\chi(z), \chi') B(\chi', k), \quad (32)$$

$A(\chi, \chi') \equiv \Delta\chi' \frac{F_K(\chi, \chi')}{a(\chi')}$, where $\Delta\chi'$ is the spacing between the sampled points in $\chi'$ and $B(\chi, \chi') \equiv j_{\ell}(k\chi')\sqrt{P(k; \chi')}.$

To further speed up computations all Bessel functions data is pre-computed in GLaSS. To save memory, values of the Bessel functions $j_{\ell}(x)$ are stored in a 2D look up table in $\ell$ and $x$ and the $j_{\ell}(k\chi)$ are found as needed. This procedure was first described in [48,49].

While computing the lensing spectra in terms of nested matrix multiplications allows for easy parallelization, this

procedure does not efficiently sample the $z-k$ space as efficiently as the Levin integration. At high-$\ell$ where the Bessel functions oscillate quickly this means the lensing spectra must be evaluated at very high resolutions.

To reduce the resolution at which the lensing spectra must be evaluated, GLaSS takes the extended Limber approximation [50] above $\ell > 100$. This was shown to have negligible impact for stage IV surveys [51]. Taking the Limber approximation, equation (11) can be rewritten as:

$$U_{\ell}(\chi, k) = \frac{F_k(\chi, \nu(k))}{ka(\nu(k))} \sqrt{\frac{\pi}{2(\ell + 1/2)}} \sqrt{P(k, \nu(k))}, \quad (33)$$

where $\nu(k) \equiv \frac{\ell+1/2}{k}$. Meanwhile at low-$\ell$ the Bessel functions oscillate slowly and the nested integrals can be evaluated at lower resolution.

## IV. CODE COMPARISON

In the following, we compare the predictions for the 3D cosmic shear covariance matrices produced with the Levin method and with the algorithm implemented in the GLaSS code. For the code comparison we fix the fiducial cosmological model to a flat cosmology with parameters given in Table I. The source distribution and the redshift error probability need to be the same for the two codes. For the source distribution, we follow Amendola et al. [52] and choose

$$n_z(z) \propto n_0 \left(\frac{\sqrt{2}}{z_m}\right)^3 z^2 \exp\left[-\left(\frac{\sqrt{2}z}{z_m}\right)^{3/2}\right], \quad (34)$$

where $z_m$ is the median redshift of the survey and $n_0$ is the observed redshift-integrated source density. We set $z_m = 0.9, n_0 = 30 \text{ arcmin}^{-2}$. We take the redshift error distribution to be a Gaussian

$$p(z_p|z) = \frac{1}{\sqrt{2\pi}\sigma(z)} \exp\left[-\frac{(z_p - z)^2}{2\sigma^2(z)}\right], \quad (35)$$

with a redshift-dependent dispersion

$$\sigma(z) = \sigma_z(1+z). \quad (36)$$

We first compare the signal-to-noise curve. The cumulative signal-to-noise ratio, summed over the contributions

<table><caption>TABLE I. Values of the cosmological parameters in the fiducial model assumed for the code comparison.</caption><thead><tr><th>Ω<sub>m</sub></th><th>Ω<sub>b</sub></th><th>Ω<sub>r</sub></th><th>Ω<sub>k</sub></th><th>w<sub>0</sub></th><th>w<sub>a</sub></th><th>σ<sub>8</sub></th><th>n<sub>s</sub></th><th>h</th></tr></thead><tbody><tr><td>0.315</td><td>0.0486</td><td>9.187 × 10<sup>-5</sup></td><td>0.0</td><td>-1.0</td><td>0.0</td><td>0.834</td><td>0.962</td><td>0.674</td></tr></tbody></table>
<figure><img src="image_2.png" /><figcaption>FIG. 1. Comparison of the differential signal-to-noise curve [Eq. (37)] as a function of the angular multipole. The two curves have been obtained from the signal and noise parts of the covariance matrices produced with GLaSS (red) and the Levin method (blue).</figcaption></figure>

at different multipoles up to a maximum multipole $\ell$, is defined as

$$ \begin{aligned} \Sigma^2(\le \ell) &= f_{\text{sky}} \sum_{\ell'=\ell_{\min}}^{\ell} \frac{2\ell'+1}{2} \text{Tr}[C_{\ell'}^{-1} S_{\ell'} C_{\ell'}^{-1} S_{\ell'}] \\ &\equiv \sum_{\ell'=\ell_{\min}}^{\ell} \Delta\Sigma^2(\ell'), \end{aligned} \quad (37) $$

where $S$ is the signal covariance [Eq. (8)] only, while $C$ refers to the sum of signal and shot noise, i.e., Eqs. (8)+(12). The signal-to-noise curves produced by both codes are shown in Fig. 1, depicting the differential contributions to the signal-to-noise coming from the different multipoles. The Levin and GLaSS predictions for the signal-to-noise curves show agreement with each other for the multipoles considered, reaching differences below 4%, as evidenced by Fig. 2 where the relative difference between the predictions of the two codes are shown. GLaSS has slightly lower signal-to-noise at high $\ell$. This is because GLaSS is not designed specifically for 3D cosmic shear and the signal-to-noise converges as the resolution of the computation grid is increased. This problem will be exacerbated as one includes higher and higher $\ell$ (e.g., $\ell = 3000$ where there is still useful signal to add) because the Bessel functions oscillate more quickly. For our comparison we used 2000 $k$-modes linearly spaced between $k = 0.005 \text{ h/Mpc}$ and $k = 2.0 \text{ h/Mpc}$, and restricted the comparison to multipoles $\ell \le 1000$; see Taylor et al. [25] for details on how much information is captured by GLaSS at different $k$-resolutions and Taylor et al. [24] for a discussion of the run time at different resolutions.

As a second diagnostic for our comparison, we consider individually the signal and noise contributions to the

<figure><img src="image_3.png" /><figcaption>FIG. 2. Relative difference of the signal-to-noise curve calculated with the GLaSS and Levin method, as a function of the multipole $\ell$.</figcaption></figure>

covariance matrices (Eqs. (8) and (12), respectively) for two different multipoles, $\ell = 100$ and $\ell = 500$. For both signal and noise we compare the elements on the diagonal $C_{\ell}(k, k)$, and plot them respectively in Figs. 3 and 4. In the noise case, we also multiply the curves by $k^2$, to check that they effectively become flat as expected. The predictions show good agreement, with differences of at most a few percent (in the lower $k$ range for the signal, and over the entire $k$ range for the noise), as visible also from Fig. 5, where we show the differences between the codes, normalized to the sum of their predictions. The disagreement in the signal plot towards the higher end of the $k$ range is due to the numerical noise present in the GLaSS computations; however, this discrepancy can be disregarded because the contributions from those $k$-regimes ($k \gtrsim 0.2 \text{ h/Mpc}$ for $\ell = 100$, $k \gtrsim 0.4 \text{ h/Mpc}$ for $\ell = 500$) are many orders of magnitude smaller than the main contributions around the peak of the curves, and also much smaller than contributions from the noise (cf. Fig 4). For $\ell = 100$, the Levin and GLaSS predictions coincide until approximately $k \approx 0.2 \text{ h/Mpc}$: at this point the behavior of the curve for GLaSS starts being dominated by numerical noise, while the Levin signal decreases in a smoother way. The same happens for $\ell = 500$, but the disagreement starts at approximately $k \approx 0.4 \text{ h/Mpc}$. In both cases, however, the signal predictions in those $k$-regimes are at least 3–4 orders of magnitude smaller than the contributions around the peak of the curves, just before and after approximately $k \approx 0.1 \text{ h/Mpc}$, respectively. Importantly, the values of the signal curves for those $k$-regimes are even smaller than the contributions from the noise, which dominates in that regime by many orders of magnitude. This means that for practical purposes we can safely ignore the contributions from those $k$-regimes where the codes are apparently in disagreement in their signal predictions. In Figs. 3 and 4 and the left panels of Fig. 5, we demonstrate this point by shading the regions where the
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

FIG. 3. Comparison of the diagonal elements of the signal part of the covariance matrices [Eq. (8)] for two multipoles $ℓ = 100$ and $ℓ = 500$, produced with GLaSS (solid lines, cyan and red for $ℓ = 100$ and $ℓ = 500$, respectively) and the Levin method (dashed lines, blue and black for $ℓ = 100$ and $ℓ = 500$, respectively). All curves have been plotted without performing any interpolation. We show the same curves using a linear (upper panel) and a logarithmic (bottom panel) scale on the y axis. The differences at higher $k$ ($k \gtrsim 0.2$ h/Mpc for $ℓ = 100$, $k \gtrsim 0.4$ h/Mpc for $ℓ = 500$) arise from the higher numerical noise present in the GLaSS computations in that $k$ regime. However, these contributions are many orders of magnitude smaller than the main contributions around the peaks of the curves, and much smaller than the contributions from the noise (cf. Fig. 4), therefore can be safely neglected. We demonstrate this point in the upper panel by indicating the shaded region for each multipole $ℓ$ where the signal represents a fraction $\le 1/1000$ of the noise: these regions correspond to the $k$-ranges where the GLaSS and Levin predictions for the signal are in apparent disagreement (cf. also Fig. 5). In both panels, the curves for $ℓ = 500$ have been multiplied by a factor 1000 for easier visualization.

signal contribution represents a fraction $\le 1/1000$ of the noise contribution at the same $k$. These regions turn out to be the same where the signal predictions of the two codes disagree, thus demonstrating that this discrepancy can be

<figure><img src="image_6.png" /></figure>

FIG. 4. Comparison of the diagonal elements of the noise part of the covariance matrices [Eq. (12)] for the same two multipoles $ℓ = 100$ and $ℓ = 500$ of Fig. 3, produced with GLaSS (solid, cyan for $ℓ = 100$ and red for $ℓ = 500$) and the Levin method (dashed, blue for $ℓ = 100$ and black for $ℓ = 500$). All curves have been plotted without performing any interpolation and multiplied by a factor $k^2$. In the case $ℓ = 500$, the curves produced by both methods have also been multiplied by a factor 20 for easier visualization.

safely disregarded. In the bottom panel of Fig. 3, we plot the same comparison between the signal predictions produced by both two codes, with a linear scale on the y axis instead of the logarithmic one used in the top panel; this is another way to appreciate how subdominant the contributions coming from the higher end of the $k$-range are with respect to the signal coming from the lower $k$-range.

It is interesting to note that the disagreement is practically only evident in the signal predictions, while the noise part is much less affected. This may be due to the increased number of matrix multiplications that need to be performed in the calculation of the signal with respect to the noise [cf. Eqs. (8) and (12)]. The fact that the number of integrations to carry out for the noise is higher means, in the GLaSS implementation, that more matrix multiplications are required and these are sensitive to the resolution in $k$. Additionally, since in the noise part of the covariance matrix there are no multiplications by Bessel functions, this may suggest that the spikes at high-$k$ in the signal may also be due to the Bessel function resolution breaking down (as explained in Sec. III B, in GLaSS the Bessel functions $j_ℓ(x)$ are precomputed in a look up table in $ℓ$ and $x$). The code implementing the Levin method sources the matter power spectrum from the Einstein-Boltzmann solver Cosmic Linear Anisotropy Solving System [CLASS, 53], while for this code comparison the matter power spectrum used by GLaSS has been sourced from the Code for Anisotropies in the Microwave Background [CAMB, 54]. The CLASS and CAMB codes have been compared in their predictions
<figure><img src="image_7.png" /></figure>

FIG. 5. Differences between the predictions for the signal (left panels) and noise (right panels) contributions to the covariance matrices for multipoles $\ell = 100$ (top panels) and $\ell = 500$ (bottom panels), normalized to their sum. We stress here again that the discrepancies at high $k$ should not be a concern because the $k$-regimes where they originate produce contributions very much subdominant with respect to the peaks of the signal curves, and also with respect to the relevant contributions from the noise (cf. Figs. 3 and 4). In the signal plots, we shade the regions where the signal is a fraction $\le 1/1000$ of the noise (cf. Fig. 3): these regions correspond to the $k$ values where the differences between the two codes are bigger; however, since the signal contributions from these regions are negligible, this discrepancy can be safely ignored.

[55]. Therefore, in comparing the Levin and GLaSS methods, the matter power spectrum has been ruled out as a possible source of discrepancy.

We conclude this section with a note on the performance of the two codes. The code implementing the Levin integration has been developed explicitly for the production of precise 3D cosmic shear cosmological forecasts and has been recently used to this purpose in Spurio Mancini et al. [22]. As one can see from Figs. 3 and 4, the curves produced with the Levin method are very smooth, showing the high precision achieved by the method. This compensates for the relatively low speed of the code, necessary to achieve that precision. GLaSS on the other hand, has not been developed for 3D cosmic shear only; in Taylor et al. [24,25] it is introduced as a means to compute lensing spectra for arbitrary weighting functions and, importantly, for integration within the cosmological module

Cosmosis. This means that speed has been a crucial goal in developing the code and the method used for the matrix multiplications indeed allows for greater speed than the one achieved with the Levin method. However, numerical noise remains higher: to overcome this issue, one would need to increase the resolution at which the matrix multiplications are performed, but this would inevitably imply a slower performance of the code. We conclude that the use of the Levin or the GLaSS method depends on the task to perform: if a high level of precision is required, the Levin method should be preferred, while if speed is a crucial requirement, GLaSS can be a better option. For our purposes in this paper, i.e., the demonstration of a method for generating 3D lensing random fields on the sky and the calculation of Minkowski functionals associated with these fields, both methods are equally valid for the computation of the 3D cosmic shear
covariance matrices, which represent the starting point of
the algorithms described in the following sections.

V. GENERATION OF SPIN-2 RANDOM
FIELDS ON THE SKY

In this section, we show how to generate random fields
on the sky starting from the full 3D cosmic shear covari-
ance matrix. As shown in Sec. II, the full covariance matrix
can be decomposed in C<sub>ℓ</sub>(k, k') for each multipole ℓ, given
that the assumed isotropy of the shear field implies multi-
pole independence ⟨γ<sub>ℓm</sub>(k)γ<sub>ℓ'm'</sub>(k')⟩ = C<sub>ℓ</sub>(k, k')δ<sub>ℓℓ'</sub>δ<sub>mm'</sub>.
We detail our procedure considering for simplicity the
convergence κ, as it is a scalar field and therefore easier to
analyze. The convergence shares essentially the same
covariance matrix with the shear field, each ℓ-block only
being rescaled by a prefactor ℓ(ℓ+1)/(ℓ+2)(ℓ−1) [37] that plays a role
only for the very largest angular scales. The generalization
to the spin-2 case for the shear field simply requires starting
from the original shear covariance matrix and replacing
the transforms from Fourier coefficients to configuration
space with their spin-2 extensions. To demonstrate our field
generation procedure we use the covariance matrices
produced with the Levin integration.

Our aim is to generate modes of the convergence field
in Fourier space κℓm and to transform them back into
configuration space using the HEALPix [56] in-built
function alm2map, in its scalar version for the conver-
gence case (for the shear, one simply needs to activate the
option pol=True that allows the user to deal with spin-2
fields). This way we can obtain samples of the convergence
field in configuration space κ(r, θ, φ), on spherical shells
corresponding to different values of the radius; on each
shell, the field can be discretized on a HEALPix map (an
example of the final result is given in Fig. 6). The procedure
described in the following is similar to the one used in the
code FLASK [57] to generate samples of the density,
convergence and shear fields on redshift slices, starting
from tomographic weak lensing covariance matrices
C<sub>ij</sub>(ℓ), where the indices i and j run over the redshift
slices and the type of field (density, convergence or shear).
In FLASK, the problem of generating correlated random
fields across different redshift slices is dealt with by means
of a Cholesky decomposition of the correlated covariance
matrices C<sub>ij</sub>(ℓ). The Cholesky decomposition rewrites the
covariance as the product of an upper and lower triangular
matrix. Here, the situation is similar in that we also have
correlated multipoles belonging to the different radial
slices; however, the correlation is in terms of the wave
vector k rather than the tomographic/field index i. This
difference originates from the fact that we start from the 3D
cosmic shear covariance matrices C<sub>ℓ</sub>(k, k'), as opposed to
the tomographic C<sub>ij</sub>(ℓ) matrices in FLASK. Additionally,
in FLASK correlations between density, convergence and
shear fields can be considered if the user desires, while here

<figure><img src="image_8.png" /><figcaption>FIG. 6. Convergence field sampled at three different values of the radius $\chi$, with the observer situated in the center. A section of the outer and middle sphere has been removed to facilitate visualization. The lensing covariance matrix which we used for sampling the random field is given by Eq. (8). We consider only contributions from the signal part of the covariance matrix, and use 30$\ell$ modes ranging between 10 and 1000. We use a linear matter power spectrum for the calculation of the covariance.</figcaption></figure>

we concentrate singularly on the generation of convergence
or shear fields and do not consider their cross-correlations.
The fact that the random fields at different wave vectors are
correlated is ultimately due to the fact that the lensing field
is not homogenous along the line of sight, due to the
mode-coupling effect of the lensing kernel, the source
redshift distributions and the redshift error probability
[cf. Eqs. (8)–(11)].

The assumption of statistical isotropy implies that modes
κℓm of the convergence field at different multipoles ℓ and m
can be generated independently. The number of ℓ multi-
poles is in principle infinite; however, practically, there will
be a ℓmax which sets the maximum resolution. We use
ℓmax = 3Nside, where Nside is a HEALPix parameter
describing the resolution of the HEALPix grid [56]. The
choice ℓmax = 3Nside is the same made by Lim and Simon
[58] in their CMB analysis and guarantees that the grid
size is comparable to the smallest angular scale considered,
corresponding to ℓmax. For each ℓ value, m ranges from
−ℓ to +ℓ, so that there are 2ℓ + 1 m values for each
multipole ℓ. However, due to the hermiticity of the con-
vergence field, we actually consider only ℓ + 1 modes
from 0 to ℓ. We employ a Cholesky decomposition of the
covariance matrices to deal with the fact that modes
corresponding to different k values are correlated:

$$
C_{\ell}(k, k') = \sum_{p} T_{\ell}(k, p) T_{\ell}(p, k'), \quad (38)
$$
Algorithm 1. Algorithm for generation of lensing Gaussian random fields on spherical shells.

**input:** Covariance matrix (e.g., for the convergence)

$$ \langle \kappa_{\ell m} \kappa_{\ell' m'} \rangle = C_{\ell}(k, k') \delta_{\ell \ell'} \delta_{mm'} $$

**output:** $\kappa(r, \theta, \phi)$. For each fixed radius $r$ in $r_1, ..., r_{N_\ell}$, create a HEALPix map on discretized $\theta$ and $\phi$

**method:** $\forall r \in [r_1, ..., r_{N_\ell}]$:

$$ \forall \ell \in [0, \ell_{\max}] $$

$$ \text{Cholesky decompose } C_{\ell} = T_{\ell} T_{\ell}^{T}; $$

$$ \forall m \in [0, \ell]: $$

$$ \text{sample } z \sim N(0, I); $$

$$ \kappa_{\ell m}(k) = T_{\ell} z; $$

$$ \kappa_{\ell m}(k) \rightarrow \int dk k^2 j_{\ell}(kr) \rightarrow \kappa_{\ell m}(r); $$

$$ \kappa_{\ell m}(r) \rightarrow \text{HEALPix alm2map} \rightarrow \kappa(r, \theta, \phi) $$

where $T(\ell)$ are (lower) triangular matrices, which we can later use to generate correlated random variables $\kappa_{\ell m}(k)$, e.g., Gaussian distributed,

$$ \kappa_{\ell m}(k) = \sum_{p} T_{\ell}(k, p) n_{\ell m}(k), \qquad (39) $$

where $n_{\ell m}(k)$ are independent, Gaussian distributed random variables with zero mean and unit variance. To obtain our convergence field samples in configuration space, we transform back from Fourier space, first by multiplying by a spherical Bessel function and $k^2$ and integrating over $k$, as indicated by Eq. (3), and then acting with the HEALPix routine `alm2map` to obtain the field samples on a discretized grid in the angular coordinates. This way we are essentially performing in two steps an inverse spherical-Bessel transform.

We summarize schematically our procedure in Algorithm 1. We implemented it in a Python routine, leveraging parallelization on multiple cores with `joblib`. The problem is embarrassingly parallel, since the correlation of the fields on different radii is preserved by the starting cosmic shear covariance matrix, while different realizations of the random fields are completely independent. The fact that the covariance does not depend on $m$, but only on the multipole $\ell$, can be used to speed up calculations, as one needs to perform the Cholesky decomposition only once per each multipole $\ell$, and can then use the decomposition for all $m$'s pertaining to that $\ell$ mode.

We show an example of 3D reconstruction on three slices in redshift (or, equivalently, comoving distance) in Fig. 6.

## VI. MINKOWSKI FUNCTIONALS OF SCALAR FIELDS ON THE SPHERE

In this section, we briefly introduce Minkowski functionals (MFs) and apply them to the generated random fields introduced in Sec. V. For Gaussian random fields MFs can be calculated analytically. We will compare this theoretical prediction with the MFs calculated directly from

the HEALPix maps as a proof of concept. In particular, we will calculate the MFs on spheres of different radii (compare Fig. 6) and estimate the covariance between the different MFs at those radii. Repeating this whole procedure for different starting lensing covariance matrices (e.g., varying each time one parameter), we can then produce a likelihood function dependent on the underlying cosmology.

### A. General definition

Here we define the MFs, concentrating on the aspects that are more interesting for cosmological applications and referring the reader to e.g., Mecke et al. [59] for further mathematical details. In our definitions, we follow the notation of Schmalzing and Gorski [60] and Lim and Simon [58].

MFs are quantities that characterize the morphology of sets in an $n$ dimensional space $\mathcal{M}$. To be considered morphological, a quantity needs to be invariant under translation and rotations; while one could think that there are many such quantities, Hadwiger's theorem [61] states that in an $n$ dimensional space there exist only $n+1$ linearly independent morphological functionals, from which all the others can be derived. These are the so-called Minkowski functionals. On the two-dimensional sphere, $S_2$, there are $2+1=3$ MFs, carrying clear geometrical interpretations and in particular representing the area, circumference and integrated geodesic curvature of an excursion set, i.e., a region where the field exceeds some threshold level. Given a threshold $\nu$ and a smooth scalar field $u$, the excursion set $Q_\nu$ is mathematically defined as

$$ Q_{\nu} = \{x \in \mathcal{M} | u(x) > \nu\}, \qquad (40) $$

while its boundary $\partial Q_\nu$ is given by

$$ \partial Q_{\nu} = \{x \in \mathcal{M} | u(x) = \nu\}. \qquad (41) $$

When considering the 2-sphere $S_2$, the first MF $V_0(\nu)$ can be interpreted as the area of $Q_\nu$,

$$ V_0(\nu) := \frac{1}{4\pi} \int_{S_2} d\Omega \Theta(u - \nu), \qquad (42) $$

where $\Theta$ is the Heaviside function. The total length of the boundary of $Q_\nu$ gives the second MF

$$ V_1(\nu) := \frac{1}{16\pi} \int_{\partial S_2} dl = \frac{1}{16\pi} \int_{S_2} d\Omega \delta(u - \nu) |\nabla u|. \qquad (43) $$

Here $\delta$ is the delta distribution and $|\nabla u|$ is the norm of the gradient of $u$. Finally, the third MF is the integral of the quantity $\kappa$ along the boundary
$$V_2(\nu) := \frac{1}{8\pi^2} \int_{\partial\mathcal{S}_2} dl\kappa = \frac{1}{8\pi^2} \int_{\mathcal{S}_2} d\Omega \delta(u - \nu) |\nabla u|\kappa, \quad (44)$$

where $\kappa$ is the geodesic curvature: this describes how much a curve $\gamma$ is different from a straight line, i.e., from a geodetic. For a normalized tangent, i.e., $|\dot{\gamma}| = 1$, the curvature is defined through

$$\kappa := |\nabla_{\dot{\gamma}} \dot{\gamma}|, \quad (45)$$

where $\nabla_{\dot{\gamma}}$ represents the covariant derivative along the tangent vector $\dot{\gamma}$ of the curve. Schmalzing and Gorski [60] show how to calculate $\kappa$ on a generic manifold, which in the case of $\mathcal{S}_2$ reads

$$\kappa = \frac{2u_{;\theta}u_{;\phi}u_{;\theta\phi} - u_{;\theta}^2u_{;\phi\phi} - u_{;\phi}^2u_{;\theta\theta}}{u_{;\theta}^2 + u_{;\phi}^2}, \quad (46)$$

with the semicolon indicating a covariant derivative.

## B. Numerical calculation of MFs

Numerical estimates of the MFs calculated on our realizations of the lensing fields can be obtained using the software HEALPix [56], as we explain in the following. We first generate full sky maps of e.g., the convergence field, on concentric spherical shells at different radii, starting from the 3D covariance matrices; to this purpose we follow the procedure described in Sec. V. We then calculate numerically the MFs by directly implementing the integrals in Eqs. (42)–(44); our algorithm closely follows the one used in Schmalzing and Gorski [60] and Lim and Simon [58] and is reported in the following.

Given the values $u(x_i)$ of a field on a pixelated map, HEALPix provides useful routines that allow for the calculation of first and second partial derivatives at each pixel in $(\ell, m)$ spherical harmonic space. We use this to obtain the three numerical MFs for $\mathcal{S}_2$, which we label $V_i(i = 0, 1, 2)$, via a sum over all pixels

$$V_i(\nu) := \frac{1}{N_{\text{pix}}} \sum_{j=1}^{N_{\text{pix}}} \mathcal{I}_i(\nu, x_j) \quad (47)$$

of the respective integrands

$$\mathcal{I}_0(\nu, x_j) := \Theta(u - \nu), \quad (48)$$

$$\mathcal{I}_1(\nu, x_j) := \frac{1}{4} \delta(u - \nu) \sqrt{u_{;\theta}^2 + u_{;\phi}^2}, \quad (49)$$

$$\mathcal{I}_2(\nu, x_j) := \frac{1}{2\pi} \delta(u - \nu) \frac{2u_{;\theta}u_{;\phi}u_{;\theta\phi} - u_{;\theta}^2u_{;\phi\phi} - u_{;\phi}^2u_{;\theta\theta}}{u_{;\theta}^2 + u_{;\phi}^2}, \quad (50)$$

where the semicolon indicates again a covariant derivative,

$$\begin{aligned} u_{;\theta} &:= \partial_{\theta}u, & u_{;\phi} &:= \frac{1}{\sin\theta}\partial_{\phi}u, \\ u_{;\phi\phi} &:= \frac{1}{\sin^2\theta}\partial_{\phi}^2u + \frac{\cos\theta}{\sin\theta}\partial_{\theta}u, & & \end{aligned} \quad (51)$$

$$\begin{aligned} u_{;\theta\theta} &:= \partial_{\theta}^2 u, & u_{;\theta\phi} &:= \frac{1}{\sin\theta}\partial_{\theta}\partial_{\phi}u - \frac{\cos\theta}{\sin^2\theta}\partial_{\phi}u, \\ u_{;\theta} &:= \frac{\partial}{\partial\theta}. & & \end{aligned} \quad (52)$$

The integrands $I_1$ and $I_2$ involve the delta function: to approximate this numerically, Schmalzing and Gorski [60] and Lim and Simon [58] use the Heaviside function

$$\delta_N(x) := (\Delta\nu)^{-1}[\Theta(x + \Delta\nu/2) - \Theta(x - \Delta\nu/2)]. \quad (53)$$

This approximation of the delta function produces some numerical noise, which Lim and Simon [58] demonstrate to be due to the delta function discretization rather than some random noise which should disappear averaging over $n_R$ realizations. For our purposes, we do not consider the corrections proposed by Lim and Simon [58] to remove this discretization effect and simply average over many realizations of the field. This is enough for our purposes, as our main goal is to test the field generation procedure rather than using the MFs to study e.g., non-Gaussianity as in Lim and Simon [58] (in which case these corrections should be taken into account).

For Gaussian fields, as the ones we are considering here, the expectation values for the MFs are known analytically and equal to

$$\bar{V}_0^G(\nu) := \langle V_0^G(\nu) \rangle = \frac{1}{2} \left( 1 - \text{erf} \left( \frac{\nu - \mu}{\sqrt{2\sigma}} \right) \right), \quad (54)$$

$$\bar{V}_1^G(\nu) := \langle V_1^G(\nu) \rangle = \frac{1}{8} \sqrt{\frac{\tau}{\sigma}} \exp \left( - \frac{(\nu - \mu)^2}{2\sigma} \right), \quad (55)$$

$$\bar{V}_2^G(\nu) := \langle V_2^G(\nu) \rangle = \frac{1}{(2\pi)^{3/2}} \frac{\tau \nu - \mu}{\sigma \sqrt{\sigma}} \exp \left( - \frac{(\nu - \mu)^2}{2\sigma} \right). \quad (56)$$

Therefore we can compare our numerical estimates with the theoretical expectation values as a check for the validity of our field generation procedure. We perform this comparison in Figs. 7–9, where we overplot our numerical estimates and their expectation values. We consider all three MFs and show the comparison for five values of the radii, corresponding to five concentric shells over which we generate our lensing field. We calculate our MFs over a set of thresholds that always ranges between $-4\sqrt{\sigma}$ and
<figure><img src="image_9.png" /><figcaption>FIG. 7. Numerical estimations of the first MF $V_0^G$ (dots), calculated on our generated Gaussian fields at different values of the radius (represented by different colors), compared with the theoretical predictions given by Eq. (54) (joined by lines), as a function of the threshold $\nu$. The range of the thresholds always varies between $-4\sqrt{\sigma}$ and $+4\sqrt{\sigma}$, where $\sigma$ is the (average) variance of the lensing field at a fixed radius.</figcaption></figure>

$+4\sqrt{\sigma}$, where $\sigma$ is the variance of the lensing field at a certain radius.

The error bars associated with our numerical estimates of the MFs are taken as the square root of the diagonal elements of the covariance matrix of the MFs, computed as

$$ \mathrm{Cov}_{ij} = \frac{1}{n_R - 1} \sum_{m=1}^{n_R} (V_i^m - \langle V_i \rangle)(V_j^m - \langle V_j \rangle), \quad i, j = 0, \dots, 3 \cdot n_\chi \cdot n_\nu \quad (57) $$

<figure><img src="image_10.png" /><figcaption>FIG. 8. Numerical estimations of the second MF $V_1^G$ (dots), calculated on our generated Gaussian fields at different values of the radius (represented by different colors), compared with the theoretical predictions. The color scheme is the same as in Fig. 7.</figcaption></figure>

where the indices $i, j$ run over the type of Minkowski functional (the three MFs $V_0, V_1, V_2$), the number of radii $n_\chi$ and the number of thresholds $n_\nu$. $\langle V_i \rangle$ denotes the mean of the MFs over all realizations $n_R$, $\langle V_i \rangle = \frac{1}{n_R} \sum_{m=1}^{n_R} V_i^m$. An example of this matrix is presented in Fig. 10. We consider the covariance between all three MFs ($V_0, V_1, V_2$, each of them as a function of the threshold $\nu$, ranging from $\nu_1$ to $\nu_{max}$), and we include the correlations between MFs belonging to each different radius (labeled by different $\chi$ value, from $\chi_1$ to $\chi_{max}$). The values of the radii are the same used for Figs. 7–9. We stress here that the error bars depicted in Figs. 7–9, associated with the MFs calculated for each value of the threshold are not independent. Also, since the realizations of the random fields on different radial shells are not statistically independent, the MFs on different radii are not independent either. To give a flavor of the correlations between the MFs, in Fig. 11 we show the elements of the correlation matrix, i.e., the Pearson correlation coefficients $r_{ij}$ calculated as

$$ r_{ij} = \frac{\mathrm{Cov}_{ij}}{\sqrt{\mathrm{Cov}_{ii}}\sqrt{\mathrm{Cov}_{jj}}}, \quad i, j = 0, \dots, 3 \cdot n_\chi \cdot n_\nu \quad (58) $$

from the covariance matrix $\mathrm{Cov}_{ij}$ of our numerical estimates of the MFs, Eq. (57). It follows from the Cauchy-Schwarz inequality that $-1 \le r_{ij} \le 1$. A value of $r_{ij} = 1$ indicates a perfect linear correlation between the two variables $i$ and $j$; a common interpretation is that in this case all data points in a sample lie on a straight line. This is also true if $r_{ij} = -1$, but the slope of the line is negative. A vanishing correlation coefficient implies that there is no linear correlation. If the correlation coefficient is positive, deviations of both variables from the mean tend to have

<figure><img src="image_11.png" /><figcaption>FIG. 9. Numerical estimations of the third MF $V_2^G$ (dots), calculated on our generated Gaussian fields at different values of the radius (represented by different colors), compared with the theoretical predictions. The color scheme is the same as in Fig. 7.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>FIG. 10. Covariance matrix between different MFs at different radii. We consider correlations between all three MFs $V_0, V_1, V_2$, all functions of the threshold $\nu$ (ranging from $\nu_1$ to $\nu_{max}$), as calculated at different radii (labeled by different $\chi$ values, ranging from $\chi_1$ to $\chi_{max}$ and specifically equal to 1000, 2000, 3000, 4000 and 5000 Mpc/h, as in Figs. 7–9). In the matrix, we indicate the block submatrices that represent the covariance between the three MFs. We used a logarithmic scale for both positive and negative values to highlight the many orders of magnitude spanned by the entries of the matrix and the different contributions given by the three MFs.</figcaption></figure>

the same sign, whereas opposite signs lead to a negative correlation coefficient.

As expected, we notice in particular a strong anti-correlation for $V_0$ centered around $\nu = 0$, as was expected by looking at Fig. 7. The same is true for $V_2$ (cf. Fig. 9), while $V_1$ is strongly positively correlated (cf. Fig. 8). This high amount of (anti)correlation suggests that in the Gaussian case analyzed here it is not necessary to consider a very high number of threshold values; however, this may not be true in the non-Gaussian case, where a higher resolution in the threshold values may be important to identify non-Gaussian features. Crucially important is, in all cases, a sufficient resolution in the HEALPix maps used at the beginning for the generation of the random fields, and later for the calculation of the MFs (in our estimates, we used the HEALPix parameters $N_{side} = 256$ and $\ell_{max} = 3N_{side} = 768$). This affects considerably the speed of the numerical implementation of these computations; however as mentioned earlier in Sec. V, the generation of random fields and, separately, the calculation of the MFs (both happening at each realization and at each radius) are embarrassingly parallel problems; this can be leveraged in practical implementations by employing

parallelization across multiple cores and nodes, without the need to worry about interprocess communication.

### C. Inference from Minkowski functionals of Gaussian fields

Introduced in cosmology by Mecke *et al.* [59], the main applications of MFs so far have been as probes of primordial non-Gaussianities [60,62,63], widely used in two and three dimensions, for instance on WMAP CMB data [64], Planck CMB data [65–69] and on the SDSS galaxy catalogue [70,71]. In the CMB case, MFs are suboptimal estimators of primordial non-Gaussianity parameters, while it has been shown that polyspectra provide minimum error bars for weak levels of non-Gaussianity [72]. Nevertheless, MFs constitute an attractive alternative to an analysis with polyspectra for a number of reasons: firstly, contrary to the bispectrum, MFs are defined in configuration rather than in Fourier space, so that a robust implementation for MFs becomes in practice easier to achieve; secondly, MFs are sensitive to the full hierarchy of higher order correlations, instead of third order only, and can provide additional information on all the nonlinear
<figure><img src="image_13.png" /><figcaption>FIG. 11. Correlation matrix between different MFs at different radii; the matrix entries represent the Pearson correlation coefficient, obtained from the covariance matrix entries (the same plotted in Fig. 10) following Eq. (58). We consider correlations between all three MFs $V_0, V_1, V_2$, all functions of the threshold $\nu$ (ranging from $\nu_1$ to $\nu_{max}$), as calculated at different radii (labeled by different $\chi$ values, ranging from $\chi_1$ to $\chi_{max}$ and specifically equal to 1000, 2000, 3000, 4000 and 5000 Mpc/h, as in Figs. 7–9). In the matrix, we indicate the block submatrices that represent the correlation between the three MFs.</figcaption></figure>

coupling parameters $f_{NL}, g_{NL}, \dots$ which appear in the perturbative development of the primordial curvature perturbation [73,74]; additionally, MFs can be analytically determined for Gaussian random fields; lastly, they are additive, which makes accounting for complicated survey geometries much easier compared to estimators of polyspectra.

In this work, we propose (for the first time, as to our knowledge) MFs as an alternative probe of Gaussianity, in addition to non-Gaussianity, in the sense specified in the following. We show how, assuming our MFs to be Gaussian distributed, we can use the MFs to probe the cosmology dependence of the fields realizations. This can be leveraged in future work to develop a full cosmological inference process based on the MFs calculated on lensing fields, of which we provide a first example here.

From a Bayesian perspective, assuming that our likelihood $L(V_i|\Omega)$ (the probability of having MFs $V_i$ given the cosmological parameters $\Omega$) is Gaussian is equivalent, considering a flat prior $p(\Omega)$ on the cosmological parameters $\Omega$, to having a Gaussian posterior $p(\Omega|V_i)$, since by virtue of Bayes theorem

$$p(V_i|\Omega) \propto L(\Omega|V_i)p(\Omega). \quad (59)$$

It follows that we are allowed to consider the likelihood and the posterior equivalently. In the Gaussian case, defining $\mathcal{L} = -\ln L$ and ignoring additive constants, we have that $-2\mathcal{L} = \chi^2$, where the $\chi^2$ can be evaluated as

$$\chi^2(\theta) = \sum_{i,j=1}^{n_v \cdot n_\chi} \langle (V_i(\theta) - \langle V_i(\theta_0) \rangle) \rangle \mathrm{cov}_{ij}^{-1}(\theta_0) \times \langle (V_j(\theta) - \langle V_j(\theta_0) \rangle) \rangle, \quad (60)$$

where the averages are performed over the number of realizations $n_R$, while the indices $i, j$ run over the length of our data vector, i.e., we consider the MFs evaluated at all the $n_\nu$ thresholds and all the $n_\chi$ radii. The MFs depend on the cosmological parameters and so does the covariance matrix; for the calculation of the chi-square, we use the inverse evaluated at the fiducial model $\theta_0$.

We calculate the $\chi^2$ statistics with MFs obtained from the realizations of the lensing random fields (we consider the convergence in this example) at different values of one cosmological parameter, for simplicity. We consider 11 values of $\Omega_m$, ranging from 0.25 to 0.35 in equidistant intervals of 0.01 centered on the fiducial value of 0.3. For
each of the $\Omega_m$ values we produce our 3D cosmic shear covariance matrix following the equations in Sec. II, with either the Levin or the GLASS method. Once the full lensing covariance matrix is available, we use it to generate, according to the procedure described in Sec. V, $n_R$ realizations of the convergence field at $r_{N_\chi}$ values of the radius in configuration space. On each shell and for each realization we also calculate the associated MFs, and store them in memory. Subsequently, we use them to build the full covariance matrix, exactly as the one shown in the previous subsection; however, this time we will have one covariance matrix of the MFs for each starting value of $\Omega_m$. Inverting the covariance corresponding to our fiducial value $\Omega_m = 0.3$, we can then use it to calculate the $\chi^2$ following Eq. (60).

The calculation of this inverse covariance matrix poses a numerical problem, in that its entries are very small and standard methods such as Gaussian elimination fail in producing a sensible inverse. We use therefore a Moore-Penrose pseudo-inverse matrix [75,76], after checking that it effectively produces an inverse covariance matrix that, multiplied by the covariance, gives back the identity matrix to within numerical precision.

We calculate the $\chi^2$ isolating the different MFs in our data vector. This implies isolating from the full covariance matrix the relevant sub-blocks for the autocorrelation of $V_0, V_1$ and $V_2$ (which we will in the following schematically indicate with $\langle V_0, V_0 \rangle, \langle V_1, V_1 \rangle, \langle V_2, V_2 \rangle$, or $\text{Cov}(V_0, V_0), \text{Cov}(V_1, V_1), \text{Cov}(V_2, V_2)$). These sub-blocks can be visualized by looking at the corresponding sub-blocks in the covariance matrix plotted in Fig. 10 (e.g., the correlation $\langle V_0, V_0 \rangle$ isolates the top left corner block); in Fig. 12 we plot the $\chi^2$ curves obtained with the three MFs. We notice how the $\chi^2$ increases going from $V_0$ to $V_2$.

<figure><img src="image_14.png" /><figcaption>FIG. 12. $\chi^2$ obtained considering the covariance of different combinations of MFs, i.e., considering the three MFs singularly $(\langle V_0, V_0 \rangle \text{ (blue)}, \langle V_1, V_1 \rangle \text{ (green)} \text{ and } \langle V_2, V_2 \rangle \text{ (red)})$. Our fiducial model is represented by the choice $\Omega_m = 0.3$.</figcaption></figure>

VII. CONCLUSIONS

3D cosmic shear constitutes an alternative to a traditional tomographic analysis of a cosmic shear survey. The spherical-Bessel expansion of the shear field at the core of its formalism maximizes the amount of redshift information; however, the calculation of the covariance matrices presents numerical difficulties due to the numerous integrations over highly oscillatory functions.

In this paper, we described and compared two methods for the calculation of simulated 3D cosmic shear covariance matrices. While the first method implements the Levin technique for integration of the periodic oscillations of the Bessel functions, the second method, implemented in the code GLASS, tackles the integrations by matrix multiplications and appropriate use of the Limber approximation.

We first compared the predictions of the two codes in terms of covariance matrices and found excellent agreement. We compared the output of the codes both in terms of the total signal-to-noise ratio and the single contributions to the covariance matrices $C_\ell$, for two different values of the multipole $\ell$, for both the signal and noise parts.

Once tested the accuracy of the predictions for the covariance matrices, we used the simulated matrices to generate Gaussian lensing fields on the sky. The procedure we described, based on a Cholesky decomposition of the $C_\ell$ matrices, allowed us to generate correlated Gaussian fields at different slices in comoving distance. We remark here that in our formalism for 3D cosmic shear we ignored complications arising from survey masks, assuming full sky coverage. Masked data can be readily accounted for by applying a mixing matrix/pseudo-$C_\ell$ like formalism. For 3D cosmic shear this is described in Kitching et al. [14], where it was applied to data. Both methods for the computation of 3D cosmic shear covariance matrices studied in this paper can be trivially modified to include such an effect. In our work, we did not investigate masked data and the ability of a pseudo-$C_\ell$ method to produce power spectra that account for a mask; however, thanks to our method to generate 3D cosmic shear data, we could now test this pseudo-$C_\ell$ approach by masking the simulated data and we will investigate this in a future paper.

The generation of normal and lognormal fields (the latter being easily obtainable from the former, by exponentiation of the Gaussian maps) can be used in future work to compute a realistic covariance matrix for a full 3D cosmic shear likelihood analysis. This should improve upon e.g., the CFHTLenS analysis for 3D cosmic shear [14], where a covariance implementation similar to GLASS was used. Kitching et al. [14] constructed a likelihood, in which the parameter dependency was in the covariance rather than the mean shear transform coefficients. This could be improved by having a likelihood in which the covariance is used as the mean and the 4-point covariance of the covariance used.
The generation of normal and lognormal random fields, starting from the 3D cosmic shear covariance matrices, also constitutes the first step for the development of a Bayesian Hierarchical Model for 3D cosmic shear power spectra estimation [following e.g., the work of [27,28], and extending it to a spherical-Bessel formalism], which can be investigated in future work.

In this work, we tested our random field generation procedure by calculating Minkowski functionals associated with our Gaussian random fields and comparing them with their known expectation values. We found good agreement between our numerical estimates and their theoretical expectation values. We calculated our Minkowski functionals separately on each spherical shell; however, we stress here that the realizations of the random fields on different radial shells are not statistically independent, as one can appreciate from the correlation matrix presented in Fig. 11. Future work should concentrate on estimating the full correlation between the Minkowski functionals at different values of the radii, implementing a fully three-dimensional approach for their calculation [see e.g., [77–80], for examples of Minkowski functionals in 3D]. Producing fully 3D Minkowski functionals for a lognormal field in 3D can be used in particular to extract non-Gaussian information from the shear field.

Finally, we showed how Minkowski functionals can also be used to extract Gaussian information by means of a likelihood analysis. We show an example of this in Fig. 12, where we plot the $\chi^2$ obtained from the covariance of the different Minkowski functionals as a function of the varying cosmological parameter $\Omega_m$. This is a first example of a full cosmological inference process, making use of the Minkowski functionals, that we plan to develop in future work.

## ACKNOWLEDGMENTS

A. S. M. and R. R. acknowledge financial support from the graduate college “Astrophysics of cosmological probes of gravity” by Landesgraduiertenakademie Baden-Württemberg. P. L. T. is supported by the UK Science and Technology Facilities Council. T. K. is supported by a Royal Society University Research Fellowship. The authors thank the German DFG Excellence Initiative for funding a mobility program through the Grant “Gravity on the largest scales and the cosmic large-scale structure” within which this work originated.

## APPENDIX: SHOT NOISE IN 3D COSMIC SHEAR

Here we derive explicitly the expression for the shot noise contribution to the 3D cosmic shear covariance matrix, as given by Eq. (12). The noise contribution is

present in all the literature on 3D cosmic shear [see e.g., the seminal papers [21,32]]; however, in the context of our code comparison for 3D cosmic shear, it may arise difficulties due to the different conventions used for the spherical-Bessel formalism. We refer the reader also to Appendix C in Lanusse et al. [26], where a derivation of the shot noise term for 3D galaxy clustering is presented.

Shot noise arises by discretizing the survey in cells that either contain one or zero galaxies [81]. We will keep the discussion more general here, for a random field $f(\vec{x})$ that is discretized on our series of cells labeled by index $i$. We will later specialize to our intrinsic ellipticity field. $n_i$ represents the occupation number of the cell and $f_i$ the value of the field in cell $i$:

$$f(\vec{x}) = \sum_{i} \delta(\vec{x} - \vec{x}_{i}) n_{i} f_{i}. \qquad (A1)$$

We calculate the correlation (where $V$ is the “volume” factor for our field):

$$\langle f(\vec{x})f(\vec{x}) \rangle = \sum_{i,j} \langle \delta^D(\vec{x} - \vec{x}_i) \delta^D(\vec{x} - \vec{x}_j) n_i n_j f_i f_j \rangle \frac{1}{V^2} \qquad (A2)$$

$$= \sum_{i,j} \delta^D(\vec{x}_i - \vec{x}_j) \langle n_i^2 \rangle \langle f_i^2 \rangle \frac{1}{V^2} \qquad (A3)$$

$$= \sum_{i} \delta^{D}(0) \langle n_{i} \rangle \langle f_{i}^{2} \rangle \frac{1}{V^{2}} \qquad (A4)$$

$$= \sum_{i} \langle n_{i} \rangle \frac{\langle f_{i}^{2} \rangle}{V} \qquad (A5)$$

where we used the fact that $n_i$ and $f_i$ are uncorrelated, $\langle n_i^2 \rangle = \langle n_i \rangle$ due to Poisson sampling and we assumed that only equal cells are correlated. In the last step we used that $\delta^D(0) = V = 4\pi$. In our case, the random field we consider is the intrinsic ellipticity of the galaxies $\epsilon_S$. This is because, as already mentioned in Sec. II, we assume the observed ellipticity $\epsilon$ to be the sum of the shear $\gamma$ and the intrinsic ellipticity $\epsilon_S$, and neglect correlations between $\gamma$ and $\epsilon_S$ as given by intrinsic alignments. We denote the intrinsic ellipticity dispersion as $\sigma_\epsilon$ (with a typical value $\sigma_\epsilon \approx 0.3$). Expressing the field $f$ in a spherical basis as

$$f_{\epsilon m}(k) = \sqrt{\frac{2}{\pi}} \sum_{i} n_i f_i j_{\epsilon}(k\chi_i) Y_{\epsilon m}(\hat{n}_i), \qquad (A6)$$

and taking into account the redshift distribution of galaxies, one arrives at Eq. (12).
[1] M. Bartelmann and P. Schneider, Weak gravitational lensing, Phys. Rep. **340**, 291 (2001).

[2] H. Hoekstra and B. Jain, Weak gravitational lensing and its cosmological applications, Annu. Rev. Nucl. Part. Sci. **58**, 99 (2008).

[3] M. Kilbinger, Cosmology with cosmic shear observations: A review, Rep. Prog. Phys. **78**, 086901 (2015).

[4] J. Fluri, T. Kacprzak, R. Sgier, A. Réfrégier, and A. Amara, Weak lensing peak statistics in the era of large scale cosmological surveys, 2018.

[5] C.-A. Lin and M. Kilbinger, A new model to predict weak-lensing peak counts. I. Comparison with N-body simulations, Astron. Astrophys. **576**, A24 (2015).

[6] A. Peel, C.-A. Lin, F. Lanusse, A. Leonard, J.-L. Starck, and M. Kilbinger, Cosmological constraints with weak-lensing peak counts and second-order statistics in a large-field survey, Astron. Astrophys. **599**, A79 (2017).

[7] Z. M. Vallis, C. G. R. Wallis, and T. D. Kitching, On the effect of projections on convergence peak counts and Minkowski functionals, Astron. Comput. **24**, 84 (2018).

[8] A. Peel, V. Pettorino, C. Giocoli, J.-L. Starck, and M. Baldi, Breaking degeneracies in modified gravity with higher (than 2nd) order weak-lensing statistics, arXiv:1805.05146.

[9] D. J. Bacon, A. R. Refregier, and R. S. Ellis, Detection of weak gravitational lensing by large-scale structure, Mon. Not. R. Astron. Soc. **318**, 625 (2000).

[10] M. L. Brown, A. N. Taylor, D. J. Bacon, M. E. Gray, S. Dye, K. Meisenheimer, and C. Wolf, The shear power spectrum from the COMBO-17 survey, Mon. Not. R. Astron. Soc. **341**, 100 (2003).

[11] L. Van Waerbeke, Y. Mellier, T. Erben, J. C. Cuillandre, F. Bernardeau, R. Maoli, E. Bertin, H. J. McCracken, O. Le Fèvre, B. Fort, M. Dantel-Fort, B. Jain, and P. Schneider, Detection of correlated galaxy ellipticities from CFHT data: First evidence for gravitational lensing by large-scale structures, Astron. Astrophys. **358**, 30 (2000).

[12] C. Heymans *et al.*, CFHTLenS tomographic weak lensing cosmological parameter constraints: Mitigating the impact of intrinsic galaxy alignments, Mon. Not. R. Astron. Soc. **432**, 2433 (2013).

[13] S. Joudaki, C. Blake, C. Heymans, A. Choi, J. Harnois-Deraps, H. Hildebrandt, B. Joachimi, A. Johnson, A. Mead, D. Parkinson, M. Viola, and L. van Waerbeke, CFHTLenS revisited: Assessing concordance with Planck including astrophysical systematics, Mon. Not. R. Astron. Soc. **465**, 2033 (2017).

[14] T. D. Kitching *et al.*, 3D cosmic shear: Cosmology from CFHTLenS, Mon. Not. R. Astron. Soc. **442**, 1326 (2014).

[15] H. Hildebrandt *et al.*, KiDS-450: Cosmological parameter constraints from tomographic weak gravitational lensing, Mon. Not. R. Astron. Soc. **465**, 1454 (2017).

[16] F. Köhlinger, M. Viola, B. Joachimi, H. Hoekstra, E. van Uitert, H. Hildebrandt, A. Choi, T. Erben, C. Heymans, S. Joudaki, D. Klaes, K. Kuijken, J. Merten, L. Miller, P. Schneider, and E. A. Valentijn, KiDS-450: The tomographic weak lensing power spectrum and constraints on cosmological parameters, Mon. Not. R. Astron. Soc. **471**, 4412 (2017).

[17] E. van Uitert *et al.*, KiDS+GAMA: Cosmology constraints from a joint analysis of cosmic shear, galaxy-galaxy lensing, and angular clustering, Mon. Not. R. Astron. Soc. **476**, 4662 (2018).

[18] R. Laureijs, J. Amiaux, S. Arduini, J. Auguères, J. Brinchmann, R. Cole, M. Cropper, C. Dabin, L. Duvet, A. Ealet *et al.*, Euclid definition study report, arXiv:1110.3193.

[19] P. A. Abell, J. Allison, S. F. Anderson, J. R. Andrew, J. R. P. Angel, L. Armus, D. Arnett, S. J. Asztalos, T. S. Axelrod *et al.*, LSST Science Book, Version 2.0 (2009).

[20] W. Hu, Power spectrum tomography with weak lensing, Astrophys. J. Lett. **522**, L21 (1999).

[21] A. Heavens, 3D weak lensing, Mon. Not. R. Astron. Soc. **343**, 1327 (2003).

[22] A. Spurio Mancini, R. Reischke, V. Pettorino, B. M. Scháfer, and M. Zumalacárregui, Testing (modified) gravity with 3D and tomographic cosmic shear, Mon. Not. R. Astron. Soc. **480**, 3725 (2018).

[23] G. Pratten, D. Munshi, P. Valageas, and P. Brax, 3D weak lensing: Modified theories of gravity, Phys. Rev. D **93**, 103524 (2016).

[24] P. L. Taylor, T. D. Kitching, and J. D. McEwen, Preparing for the cosmic shear data flood: Optimal data extraction and simulation requirements for stage IV dark energy experiments, Phys. Rev. D **98**, 043532 (2018).

[25] P. L. Taylor, T. D. Kitching, J. D. McEwen, and T. Tram, Testing the cosmic shear spatially-flat universe approximation with generalized lensing and shear spectra, Phys. Rev. D **98**, 023522 (2018).

[26] F. Lanusse, A. Rassat, and J.-L. Starck, 3D galaxy clustering with future wide-field surveys: Advantages of a spherical Fourier-Bessel analysis, Astron. Astrophys. **578**, A10 (2015).

[27] 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**, 4452 (2016).

[28] 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**, 3272 (2017).

[29] B. Zieser and P. M. Merkel, The cross-correlation between 3D cosmic shear and the integrated Sachs-Wolfe effect, Mon. Not. R. Astron. Soc. **459**, 1586 (2016).

[30] Y. Ayaita, B. M. Schäfer, and M. Weber, Investigating clustering dark energy with 3D weak cosmic shear, Mon. Not. R. Astron. Soc. **422**, 3056 (2012).

[31] A. Grassi and B. M. Schäfer, Detecting baryon acoustic oscillations by 3D weak lensing, Mon. Not. R. Astron. Soc. **437**, 2632 (2014).

[32] A. F. Heavens, T. D. Kitching, and A. N. Taylor, Measuring dark energy properties with 3D cosmic shear, Mon. Not. R. Astron. Soc. **373**, 105 (2006).

[33] T. D. Kitching, A. F. Heavens, A. N. Taylor, M. L. Brown, K. Meisenheimer, C. Wolf, M. E. Gray, and D. J. Bacon, Cosmological constraints from COMBO-17 using 3D weak lensing, Mon. Not. R. Astron. Soc. **376**, 771 (2007).

[34] T. D. Kitching, A. F. Heavens, and S. Das, 3D weak gravitational lensing of the CMB and galaxies, Mon. Not. R. Astron. Soc. **449**, 2205 (2015).

[35] J. N. Goldberg, A. J. Macfarlane, E. T. Newman, F. Rohrlich, and E. C. G. Sudarshan, Spin-spherical harmonics and δ, J. Math. Phys. (N.Y.) **8**, 2155 (1967).
[36] E. Newman and R. Penrose, An approach to gravitational radiation by a method of spin coefficients, *J. Math. Phys. (N.Y.)* **3**, 566 (1962).

[37] P. G. Castro, A. F. Heavens, and T. D. Kitching, Weak lensing analysis in three dimensions, *Phys. Rev. D* **72**, 023516 (2005).

[38] M. Abramowitz, I. A. Stegun, and R. H. Romer, Handbook of mathematical functions with formulas, graphs, and mathematical tables, *Am. J. Phys.* **56**, 958 (1988).

[39] B. Leistedt, J. D. McEwen, T. D. Kitching, and H. V. Peiris, 3D weak lensing with spin wavelets on the ball, *Phys. Rev. D* **92**, 123010 (2015).

[40] T. D. Kitching and A. F. Heavens, Unequal-time correlators for cosmology, *Phys. Rev. D* **95**, 063522 (2017).

[41] P. M. Merkel and B. M. Schäfer, Intrinsic alignments and 3D weak gravitational lensing, *Mon. Not. R. Astron. Soc.* **434**, 1808 (2013).

[42] D. Levin, Fast integration of rapidly oscillatory functions, *J. Comput. Appl. Math.* **67**, 95 (1996).

[43] D. Levin, Analysis of a collocation method for integrating rapidly oscillatory functions, *J. Comput. Appl. Math.* **78**, 131 (1997).

[44] N. Kaiser, Weak gravitational lensing of distant galaxies, *Astrophys. J.* **388**, 272 (1992).

[45] N. Kaiser, Weak lensing and cosmology, *Astrophys. J.* **498**, 26 (1998).

[46] M. Loverde and N. Afshordi, Extended limber approximation, *Phys. Rev. D* **78**, 123506 (2008).

[47] J. Zuntz, M. Paterno, E. Jennings, D. Rudd, A. Manzotti, S. Dodelson, S. Bridle, S. Sehrish, and J. Kowalkowski, Cosmosis: Modular cosmological parameter estimation, *Astron. Comput.* **12**, 45 (2015).

[48] A. Kosowsky, Efficient computation of hyperspherical bessel functions, arXiv:astro-ph/9805173.

[49] U. Seljak and M. Zaldarriaga, A line-of-sight integration approach to cosmic microwave background anisotropies, *Astrophys. J.* **469**, 437 (1996).

[50] M. LoVerde and N. Afshordi, Extended limber approximation, *Phys. Rev. D* **78**, 123506 (2008).

[51] T. D. Kitching, J. Alsing, A. F. Heavens, R. Jimenez, J. D. McEwen, and L. Verde, The limits of cosmic shear, *Mon. Not. R. Astron. Soc.* **469**, 2737 (2017).

[52] L. Amendola et al., Cosmology and fundamental physics with the Euclid satellite, 2016.

[53] D. Blas, J. Lesgourgues, and T. Tram, The cosmic linear anisotropy solving system (CLASS). Part II: Approximation schemes, *J. Cosmol. Astropart. Phys.* 07 (2011) 034.

[54] A. Lewis, A. Challinor, and A. Lasenby, Efficient computation of CMB anisotropies in closed FRW models, *Astrophys. J.* **538**, 473 (2000).

[55] J. Lesgourgues, The cosmic linear anisotropy solving system (CLASS) III: Comparison with CAMB for LambdaCDM, arXiv:1104.2934.

[56] K. M. Górski, E. Hivon, A. J. Banday, B. D. Wandelt, F. K. Hansen, M. Reinecke, and M. Bartelmann, HEALPix: A framework for high-resolution discretization and fast analysis of data distributed on the sphere, *Astrophys. J.* **622**, 759 (2005).

[57] H. S. Xavier, F. B. Abdalla, and B. Joachimi, Improving lognormal models for cosmological fields, *Mon. Not. R. Astron. Soc.* **459**, 3693 (2016).

[58] E. A. Lim and D. Simon, Can we detect hot/cold spots in the CMB with Minkowski functionals?, *J. Cosmol. Astropart. Phys.* 01 (2012) 048.

[59] K. R. Mecke, T. Buchert, and H. Wagner, Robust morphological measures for large-scale structure in the Universe, *Astron. Astrophys.* **288**, 697 (1994).

[60] J. Schmalzing and K. M. Gorski, Minkowski functionals used in the morphological analysis of cosmic microwave background anisotropy maps, *Mon. Not. R. Astron. Soc.* **297**, 355 (1998).

[61] H. Hadwiger, *Vorlesungen über Inhalt, Oberfläche und Isoperimetrie* (Springer, New York, 1957).

[62] J. Schmalzing and T. Buchert, Beyond genus statistics: A unifying approach to the morphology of cosmic structure, *Astrophys. J. Lett.* **482**, L1 (1997).

[63] S. Winitzki and A. Kosowsky, Minkowski functional description of microwave background Gaussianity, *New Astron.* **3**, 75 (1998).

[64] C. Hikage, T. Matsubara, P. Coles, M. Liguori, F. K. Hansen, and S. Matarrese, Limits on primordial non-Gaussianity from Minkowski functionals of the WMAP temperature anisotropies, *Mon. Not. R. Astron. Soc.* **389**, 1439 (2008).

[65] T. Buchert, M. J. France, and F. Steiner, Model-independent analyses of non-Gaussianity in Planck CMB maps using Minkowski functionals, *Classical Quantum Gravity* **34**, 094002 (2017).

[66] A. Ducout, F. R. Bouchet, S. Colombi, D. Pogosyan, and S. Prunet, Non-Gaussianity and Minkowski functionals: Forecasts for Planck, *Mon. Not. R. Astron. Soc.* **429**, 2104 (2013).

[67] C. P. Novaes, A. Bernui, G. A. Marques, and I. S. Ferreira, Local analyses of Planck maps with Minkowski functionals, *Mon. Not. R. Astron. Soc.* **461**, 1363 (2016).

[68] P. A. R. Ade et al. (Planck Collaboration), Planck 2013 results. XXIV. Constraints on primordial non-Gaussianity, *Astron. Astrophys.* **571**, A24 (2014).

[69] P. A. R. Ade et al. (Planck Collaboration), Planck 2015 results. XVII. Constraints on primordial non-Gaussianity, *Astron. Astrophys.* **594**, A17 (2016).

[70] C. Hikage, E. Komatsu, and T. Matsubara, Primordial non-Gaussianity and analytical formula for Minkowski functionals of the cosmic microwave background and large-scale structure, *Astrophys. J.* **653**, 11 (2006).

[71] C. Park, Y.-Y. Choi, M. S. Vogeley, J. Richard Gott, III, J. Kim, C. Hikage, T. Matsubara, M.-G. Park, Y. Suto, and D. H. Weinberg (SDSS Collaboration), Topology analysis of the sloan digital sky survey. I. scale and luminosity dependence, *Astrophys. J.* **633**, 11 (2005).

[72] D. Babich, Optimal estimation of non-Gaussianity, *Phys. Rev. D* **72**, 043003 (2005).

[73] E. Komatsu and D. N. Spergel, Acoustic signatures in the primary microwave background bispectrum, *Phys. Rev. D* **63**, 063002 (2001).

[74] T. Okamoto and W. Hu, Angular trispectra of CMB temperature and polarization, *Phys. Rev. D* **66**, 063008 (2002).

[75] A. Dresden, The fourteenth western meeting of the american mathematical society, *Bull. Am. Math. Soc.* **26**, 385 (1920).

[76] R. Penrose, A generalized inverse for matrices, *Math. Proc. Cambridge Philos. Soc.* **51**, 406 (1955).
[77] S. Appleby, P. Chingangbam, C. Park, K. P. Yogendran, and P. K. Joby, Minkowski tensors in three dimensions: Probing the anisotropy generated by redshift space distortion, Astrophys. J. **863**, 200 (2018).

[78] L. Gleser, A. Nusser, B. Ciardi, and V. Desjacques, The morphology of cosmological reionization by means of Minkowski functionals, Mon. Not. R. Astron. Soc. **370**, 1329 (2006).

[79] C. Hikage, J. Schmalzing, T. Buchert, Y. Suto, I. Kayo, A. Taruya, M. S. Vogeley, F. Hoyle, J. R. Gott, III, and

J. Brinkmann, Minkowski functionals of SDSS galaxies I: Analysis of excursion sets, Publ. Astron. Soc. Jpn. **55**, 911 (2003).

[80] S. Yoshiura, H. Shimabukuro, K. Takahashi, and T. Matsubara, Studying topological structure of 21-cm line fluctuations with 3D Minkowski functionals before reionization, Mon. Not. R. Astron. Soc. **465**, 394 (2017).

[81] P. J. E. Peebles, *The large-scale structure of the universe* (Princeton University Press, Princeton, 1980).