{'fraction_non_alphanumeric': 0.04291261156774746, 'fraction_numerical': 0.011917368429899315, 'mean_word_length': 4.530724179603979, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 10, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
3D weak lensing with spin wavelets on the ball

We construct the spin flaglet transform, a wavelet transform to analyze spin signals in three dimensions. Spin flaglets can probe signal content localized simultaneously in space and frequency and, moreover, are separable so that their angular and radial properties can be controlled independently. They are particularly suited to analyzing of cosmological observations such as the weak gravitational lensing of galaxies. Such observations have a unique 3D geometrical setting since they are natively made on the sky, have spin angular symmetries, and are extended in the radial direction by additional distance or redshift information. Flaglets are constructed in the harmonic space defined by the Fourier-Laguerre transform, previously defined for scalar functions and extended here to signals with spin symmetries. Thanks to various sampling theorems, both the Fourier-Laguerre and flaglet transforms are theoretically exact when applied to bandlimited signals. In other words, in numerical computations the only loss of information is due to the finite representation of floating point numbers. We develop a 3D framework relating the weak lensing power spectrum to covariances of flaglet coefficients. We suggest that the resulting novel flaglet weak lensing estimator offers a powerful alternative to common 2D and 3D approaches to accurately capture cosmological information. While standard weak lensing analyses focus on either real or harmonic space representations (i.e., correlation functions or Fourier-Bessel power spectra, respectively), a wavelet approach inherits the advantages of both techniques, where both complicated sky coverage and uncertainties associated with the physical modeling of small scales can be handled effectively. Our codes to compute the Fourier-Laguerre and flaglet transforms are made publicly available.


I. INTRODUCTION
Weak gravitational lensing by large-scale structure, for example cosmic shear, is an observational probe that has the potential to constrain both the geometry and growth of structure of the Universe. It is a sensitive probe of dark energy physics, extensions to general relativity, and neutrino mass and hierarchy (see Refs. [1][2][3][4][5] and references therein for more details). The basic measurements of weak gravitational lensing from galaxy surveys are the thirdflattening, or third-eccentricity, of galaxy images (colloquially refereed to as 'ellipticity'), and galaxy sizes. These contain information about the intrinsic unlensed shape of the galaxies and the additional ellipticity (called 'shear'), as well as size changes caused by the weak gravitational lensing effect along the line of sight. This angular information can be supplemented with the redshift of the galaxies (either photometric or spectroscopic), yielding a catalogue of galaxy positions, shapes, sizes and redshifts. These can then be compared to predictions from cosmological models in order to constrain their parameters.
Due to the cosmological principle, i.e., assumptions of isotropy and homogeneity, the mean of the cosmic shear averaged over all galaxies in a sufficiently large survey is expected to be zero. However, the power spectrum of the shear (pair correlations of modes in 3D Fourier space) is expected to be nonzero and depends on both the geometry of the Universe and the power spectrum of matter perturbations in the 3D volume within which the galaxy sample lies. The 3D power spectrum of the galaxy shear estimates themselves is known as '3D cosmic shear'. There exist various strategies and approximations to efficiently compute the 3D cosmic shear power spectrum, including linking angular and radial scales, and binning in redshift-known as 'tomography'-that are more widely used than the full 3D case because of computational ease.
From a theoretical perspective, 3D cosmic shear has been developed in a series of papers [6][7][8][9] that introduced the relation between the underlying matter power spectrum, the geometric lensing kernel, and the on-sky measurements of galaxies shear and size changes. Because the gravitational lensing signal is caused by a tidal effect around massive on the sky and along the line of sight (e.g., the ability to distinguish galaxies from stars and to estimate their redshifts).
A number of 2D and 3D methods, including 3D wavelets [38][39][40], have been developed to analyze projected weak lensing observables as well as the polarization of the CMB, which is also a spin signal (e.g., Refs. [41][42][43][44][45]). However, these approaches do not natively deal with 3D spherical geometry and spin symmetries simultaneously. Wavelets in 3D spherical coordinates were constructed in Ref. [38] using B-spline kernels and the Fourier-Bessel basis. While these are directly connected to the Fourier-Bessel formalism, they can probe only 3D isotropic features in scalar signals. Hence, they do not account for the angular spin symmetries or the unique 2+1D geometry of weak lensing observables. By contrast, Ref. [39] introduced 'flaglets' to probe the angular and radial scales separately. This approach was based on the Fourier-Laguerre transform, a novel separable 3D harmonic transform with a sampling theorem and an analytical connection to the Fourier-Bessel basis. More recently, Ref. [40] proposed a construction of separable 3D needlets based on the Fourier transform in the radial direction.
In this paper, we extend the framework of Ref. [39] and present novel Fourier-Laguerre and flaglet transforms supporting spin symmetries and also directional features in the angular direction. We detail sampling theorems and efficient algorithms to compute these transforms exactly, exploiting the recently constructed spin directional wavelet transform on the sphere and sampling theorem on the rotation group, and their corresponding fast algorithms, developed in Refs. [46][47][48]. Thanks to the separability of the radial and angular components, the novel spin Fourier-Laguerre and flaglet transforms natively support the 2+1D spherical geometry of the weak lensing data while accounting for angular spin symmetries. They can also be related to the standard 3D weak lensing framework in the Fourier-Bessel basis thanks to the analytical connection between the Fourier-Laguerre and Fourier-Bessel transforms. Finally, we discuss the advantages of working in flaglet space (instead of the standard correlation functions, Fourier, or Fourier-Bessel approaches) that arise from the dual localization properties of flaglets in pixel and frequency space, for example in dealing with complicated angular masks or small scales where physical modeling is less certain.
This paper is structured as follows. In Sec. II we give a short review of the 3D weak lensing formalism. We adopt the language of observational cosmology and focus on the observables that can be measured from galaxy survey data and connected to cosmological models. In Secs. III and IV we step back from this formalism and give a formal mathematical presentation of the new analysis techniques (transforms and operators) for radial, spherical and 3D fields. This leads to the definition of the novel spin Fourier-Laguerre and spin flaglet transforms (in Secs. III and IV, respectively). In Sec. V, we combine the two viewpoints, apply the spin flaglet transform to 3D weak lensing, and study the properties of the shear flaglet coefficients. We present conclusions in Sec. VI. Further technical details on the special functions and useful approximations used in this paper are presented in appendixes.

II. 3D WEAK LENSING
In this section we give a brief introduction to the 3D weak lensing formalism, which has been developed in a series of recent papers [6][7][8][9]. Here we focus on the basic principles of 3D weak lensing and discuss the lensing potential, observable quantities, and how observables and theory can be compared to constrain cosmological models and parameters. For further background we refer the reader to the excellent exposition of 3D weak lensing given in Ref. [7]. Details of practical aspects and data-related complications are addressed in Ref. [12] (e.g., pseudo-C methods, real and imaginary covariance structures, E and B-mode decomposition, and weighting due to shape measurement biases). We touch these issues briefly in this section, and in Sec. V, but leave detailed investigations on data and simulations to future work.

A. Lensing potential
The weak gravitational lensing effect is commonly expressed in terms of the lensing potential φ, which depends on the integrated deflection angle along the line of sight, sourced by the local Newtonian potential Φ, where r is the comoving distance, n the angular position on the sky, and c is the speed of light in a vacuum. The geometrical factor reads for cosmologies with positive, flat and negative global curvatures, with curvature K = 1, 0, and −1, respectively. This expression assumes the linear regime and also the Born approximation, i.e., that the path of photons is unperturbed by lenses. The Newtonian potential Φ can be related to perturbations in the underlying matter density δ via Poisson's equation, where Ω M is the dimensionless matter density, H 0 is the current value of the Hubble parameter, and a(r) is the dimensionless scale factor. The 3D gradient is defined relative to comoving coordinates. This relation can be expressed conveniently in terms of the Fourier-Bessel basis since the Fourier-Bessel basis functions are eigenfunctions of the Laplacian operator ∇ 2 , with eigenvalues −k 2 . Consequently, the harmonic representation of Eq. (3) reads where and m label angular harmonic modes and k labels the radial wavenumber (in units of h −1 Mpc). The Fourier-Bessel transform is defined explicitly in Sec. III D by Eqs. (37) and (38), and discussed extensively in the same section. The dependency on r must be retained in Eq. (4) to account for the time evolution of the field: we apply the Fourier-Bessel transform to the homogenous field existing everywhere at the cosmological time corresponding to comoving distance r.
The harmonic representation of the lensing potential then reads where j are the spherical Bessel functions and R + denotes the positive real half-line, i.e., [0, ∞). The difference between Eq. (5) and the expression shown in Ref. [7] is due only to the different conventions used for the spherical Bessel transform, where we adopt the convention of Refs. [6,39] (see Appendix A and Sec. III D).

B. Observables
Weak lensing generates distortions in the observations of a background field, which may be characterised by spin quantities of several orders (see, e.g., Ref. [49]), with spin number s. In the weak lensing regime (i.e., away from the critical curve of lensing masses, where there are not multiple images of sources) four distortions can be produced [49]: the size magnification 0 κ [50]; the shear 2 γ; and the flexion, the combined effect of the first-flexion 1 F (a centroid shift) and third-flexion 3 G (a trefoil distortion). We introduce the generic notation s χ to denote these spin quantities, which are related to the lensing potential φ by (see, e.g., Ref. [7]) where ð andð are spin raising and spin lowering operators [51][52][53][54]. Further details about these operators and spin spherical harmonics in general are given in Sec. III and Appendix A. We now focus on the spin-2 quantity of shear their flaglet covariances). The shear field ±2 γ is typically decomposed into its real and imaginary parts, ±2 γ = γ 1 ± iγ 2 , and is related to the lensing potential by [7] 2 γ(n, r) = γ 1 (n, r) + iγ 2 (n, r) = ð 2 φ E (n, r) + iφ B (n, r) /2, 2 γ * (n, r) = −2 γ(n, r) = γ 1 (n, r) − iγ 2 (n, r) =ð 2 φ E (n, r) − iφ B (n, r) /2, where * denotes complex conjugation and we have represented the lensing potential by its parity even and odd components given by the E-and B-mode signals, respectively, denoted by superscripts E and B. Alternatively, the E-and B-mode signals may be related to the shear bỹ whereφ E andφ B are normalized versions of φ E and φ B and, when expressed in their Fourier-Bessel representations (see Sec. III D), are related byφ where The shear induced by gravitational lensing produces an E-mode signal only since density (scalar) perturbations cannot induce a parity odd B-mode component [7]. In the absence of systematic effects (see, e.g., Refs. [55,56]), φ E = φ and φ B = 0. The 3D cosmic shear formalism nevertheless includes B-modes since the B-mode signal-and EB-mode cross correlations-computed from data can be used as a null-test to search for residual systematics.
The spin ±2 cosmic shear signal ±2 γ can also be represented in Fourier-Bessel space, which leads to a simple harmonic connection between the shear and the lensing potential. Since shear is a spin quantity, it is decomposed into its spin Fourier-Bessel coefficients (see Sec. III D), denoted ±2 γ m (k). Scalar harmonic E-and B-mode signals may be computed from ±2 γ m (k) by where yet another normalization is assumed (as is standard practice). It also follows that By Eqs. (13) and (14) and either Eqs. (7) and (8) or Eqs. (9) and (10), it follows that the lensing potential and the shear are related in harmonic space indirectly by Consequently, the shear can be related directly to the lensing potential in Fourier-Bessel space by In fact, Eq. (18) can be recovered directly from the harmonic representation of Eq. (6). Nevertheless, we expose the E-and B-mode decomposition since it is the standard approach and is of additional practical use in studying residual systematics, as outlined. Finally, note that galaxy distances are not directly measurable; only their redshifts can be estimated from galaxy colours or spectra. Thus, the 3D shear signal ±2 γ(r) in real space is typically constructed from redshift space observations via the relation r(z), which assumes a fiducial reference cosmological model [12].

C. Cosmology
Eqs. (4) and (5) show that the lensing potential contains information about the matter fluctuations and geometry of the Universe, via δ and F K , respectively. In the former, most of the information arises from the first nonzero moment of the field, the 2-point covariance of the the fluctuations, described by the matter power spectrum P (k; r): where the angle-brackets (without a bar) indicate an ensemble average over several realizations of the field given the same cosmological model and power spectrum P (k; r), δ D denotes the 1D Dirac delta function and δ K the Kronecker delta symbol. Note that we have assumed homogeneity and isotropy and that the field is evaluated at a fixed comoving distance [7]. Higher order statistics of the density field are of less interest because they are more difficult to model and measure in data. For this reason we focus on the 2-point power spectrum of the cosmic shear and other weak lensing observables. The 3D lensing power spectrum C φφ (k, k ) is defined by and is dependent on P (k; r) and F K . Note that the 3D lensing potential is not homogenous and isotropic since it is given by a 2D projection of the Newtonian potential between an observer and source at radius r. Consequently, the lensing potential is homogenous and isotropic in the angular direction but not in the radial direction [7]. By Eq. (18), and equivalently by Eq. (6) for s = 2, the covariance of the 3D cosmic shear 2 γ in Fourier-Bessel space is related to the lensing power spectrum by Note that since the shear transform coefficients are complex quantities, their joint probability distribution must include the correlation between the real and imaginary parts [57,58]. This may be accounted for by creating an 'Affix' covariance, as shown in Ref. [12], or by adopting real spherical harmonics.
Further observational aspects such as the inclusion of photometric redshift information or a varying number-density of sources can also be included as shown in Ref. [9]. A further complicating aspect is that the observed galaxy ellipticity is not a direct measure of the shear but is a combination of the shear and the intrinsic (unlensed) ellipticity of the galaxy: e obs (n, r) = e I (n, r) + γ(n, r), to linear order. These effects can be taken into account by modeling auto-correlations between the intrinsic ellipticities e I (n, r) and cross-correlations between the intrinsic ellipticities and the shear as described in Ref. [59] (correlations with CMB lensing can also be included, as also shown in Ref. [59]). In this paper we are concerned with linking observable power spectra to the model power spectrum that should be assumed to include all of these effects. However, for simplicity and clarity we will refer to only the 'shear' power spectrum and keep the notation C φφ (k, k ). For more details on the additional modeling aspects we refer the reader to the references provided.
We do not consider higher-order statistics, such as the bi-spectrum or the tri-spectrum [59]. These can be expressed as generalizations of the power spectrum method presented. Note that in the Fourier-Bessel transforms we have assumed a flat geometry, and that the more general case requires the use of ultra spherical Bessel functions, where the geometry is now a generalized manifold. However, as pointed out in Ref. [11], the use of the spherical Bessel functions is well motivated in the limit of 1 and k (curvature scale) −1 . The standard approach of 3D cosmic shear, as performed in Ref. [12], is to measure the Fourier-Bessel coefficients of the shear 2 γ m (k) from a galaxy survey and then use them to constrain C φφ (k, k ) (hence, P (k; r) and F K ) through Eq. (21). This harmonic-space approach can be challenging due to both observational and technical complications. As mentioned previously, it is difficult to simultaneously deal with complicated masks, mode coupling and scale mixing of the Fourier-Bessel coefficients. We will come back to these issues in Sec. V when we detail our new method to constrain C φφ from flaglet coefficients.
We now step back from weak lensing and adopt a more general signal processing viewpoint to define the novel 3D spin Fourier-Laguerre and flaglet transforms (which are suitable for the analysis of arbitrary 3D spin fields). In Sec. V we apply these transforms specifically to weak lensing observables and construct an estimator related to the 3D lensing power spectrum.

III. SPIN FOURIER-LAGUERRE TRANSFORM
After concisely reviewing signals and transforms defined on the radial line R + , sphere S 2 , and rotation group SO(3), we construct an exact, separable Fourier-Laguerre transform to analyze signals of arbitrary spin defined on S 2 × R + , i.e., the 3D space formed by the product of the sphere S 2 and the radial line R + . To subsequently construct a directional wavelet transform in this space (presented in Sec. IV), it is necessary to also consider signals and transforms defined on SO(3)×R + , i.e., the space formed by the production of the rotation group SO(3) and the radial line R + . We thus also construct a Wigner-Laguerre harmonic transform on SO(3)×R + . Sampling theorems and fast algorithms to compute these transform exactly for bandlimited signals are presented. In addition, the translation and convolution operators needed to construct a 3D spin wavelet transform are outlined. Finally, we derive the relation between the Fourier-Laguerre and Fourier-Bessel transforms, allowing one to exactly compute the Fourier-Bessel coefficients of signals bandlimited in Fourier-Laguerre space. This property is of use in Sec. V to connect the Fourier-Laguerre and flaglet transforms with the 3D weak lensing formalism. Note that in what follows brackets with a separating bar denote inner products (and not ensemble averages). In other words, f |g = S dυf g * , where the integral runs over the space S of interest, i.e., where f and g are defined, with measure dυ.
A. Transforms on R + , S 2 , and SO (3) We first consider a square integrable complex signal on the radial line f ∈ L 2 (R + ), which can be expanded using the spherical Laguerre transform introduced in Ref. [39], with forward and inverse, respectively. The Laguerre basis functions are specified by where L (α) p is the p-th generalized Laguerre polynomial of order α (see Appendix A for further details about the basis functions K p ). In this paper we consider the case α = 0. Note that τ ∈ R + is a scaling factor to map the R + sampling theorem to any finite interval [0, R] of interest. We refer the reader to [39] for more details on the spherical Laguerre transform. We take r 2 dr as the natural measure on R + since we aim to construct 3D transforms in spherical coordinates, where the volume-invariant measure involves r 2 dr (in any case, an alternative measure could be adopted if it were desired).
We now consider a square integrable complex function on the sphere s f (n) ∈ L 2 (S 2 ) with n = (θ, φ), where θ ∈ [0, π] denotes colatitude and φ ∈ [0, 2π) longitude. A spin function s f with spin number s ∈ Z transforms under a local rotation by ϑ ∈ [0, 2π) in the tangent plane centered on n as s f (n) = e −isϑ s f (n), where the prime denotes the rotated function [51][52][53][54]. For s = ±2, such as the cosmic shear field, s f is invariant under local rotations of ±π. Note that the sign convention adopted for the argument of the complex exponential differs from the original definition [52] but is identical to the convention used in the context of the polarization of the CMB [53,54]. The natural harmonic transform on the sphere that accounts for spin symmetry is the spin spherical harmonics transform, with forward and inverse, respectively, where s Y m denotes the spin spherical harmonics (using the Condon-Shortley phase convention; see Appendix A for more details). The usual invariant measure of the sphere is dΩ(n) = sin θdθdφ. Note that spin signals have s f m = 0, ∀ < |s|.
The spherical Laguerre and spherical harmonics transforms can be combined naturally into a 3D separable trans- shown on the 3D Fourier-Laguerre sampling scheme with L = P = 80, plotted up to r = R/2 with R the radius of the final sample of the radial sampling scheme [39]. This plot shows that the 3D basis functions successfully probe harmonic modes in the angular and radial directions separately. Furthermore, notice that the basis functions are not well localized in space, similar to standard Fourier bases.
form, as shown in Ref. [39] for the spin-0 case. Before we extend this construction to higher spin numbers in the next section, we turn to the rotation group, SO (3), and its standard harmonic transform: the Wigner transform [48,60]. The rotation group is the natural manifold on which to construct a spherical wavelet transform probing directional features. This is due to the rotation properties of spherical harmonics and the natural convolution operator on the sphere (recalled in Sec. III E). Thus, we consider a square integrable complex signal defined on the rotation group are the Euler angles. Its forward and inverse Wigner transforms are given by respectively, where D mn are Wigner functions and the invariant measure on the rotation group reads dµ(ρ) = sin β dβdαdγ (see Appendix A for further details).
B. Fourier-Laguerre transform on S 2 × R + and Wigner-Laguerre transform on SO(3)×R + We now combine the previous transforms and extend the Fourier-Laguerre tranform presented in Ref. [39] so that its angular part can support spin signals on S 2 , in order to match the nature of weak lensing observations. More precisely, we now consider 3D complex functions in spherical coordinates s f (n, r) ∈ L 2 (S 2 ×R + ). We create a separable transform by combining the spherical Laguerre transform with the spin spherical harmonics, leading to the forward and inverse transforms respectively. The spin property implies that s f m,p = 0, ∀ < |s|.
For functions f (ρ, r) ∈ L 2 (SO(3)×R + ) we define the Wigner-Laguerre forward and inverse transforms as respectively. In these expressions, the mn indices refer to the angular modes (in the spherical harmonics or Wigner transforms), while p corresponds to the radial mode (from the spherical Laguerre transform). An example of a basis function is shown in Fig. 1. Finally, as in the individual transforms, the summations for the harmonic modes run over In what follows, we will omit these bounds for conciseness and use the terminology 'angular' for both the S 2 and SO(3) parts of the signals and transforms.

C. Sampling theorems and numerical implementation
In practice, evaluating the integrals in Eqs. (29) and (31) numerically requires quadrature rules. While generic numerical integration methods can be adopted, this problem can be tackled more accurately by appealing to sampling theorems. When considering bandlimited signals, sampling theorems lead to exact quadrature rules, which allow one to discretize the space such that the previous integrals can be computed exactly. In practice, signals are transformed with accuracy close to the level of machine precision, with errors due to the (finite, approximate) representation of floating point numbers only. This approach is particularly desirable when performing numerous transforms successively or when high numerical accuracy is required (which is typically the case in modern observational cosmology).
The radial and angular parts of the Fourier-Laguerre and Wigner-Laguerre transforms are fully separable. This implies that one can simply combine sampling theorems on the sphere, the rotation group and the radial line to obtain sampling theorems on S 2 × R + and SO(3)×R + . We now recall these sampling theorems and the notation for the sampling points and quadrature weights.
On R + , we use the spherical Laguerre quadrature rule presented in Ref. [39]. bandlimited signals have f p = 0, ∀p > P , with P the radial band-limit, while the sampling nodes and quadrature weights are denoted by r k and w k , respectively. This exact quadrature, with corresponding nodes and weights, is used to compute Eq. (22).
On the sphere, we consider the sampling theorem developed in Ref. [61]: a function s f with band limit L satisfies s f m = 0, ∀ ≥ L and can be completely described by the values taken on the equiangular, separable pixelization: n ij = (θ i , φ j ) with i = 0, . . . , L − 1 and j = 0, . . . , 2L − 1. The spherical harmonics coefficients of Eq. (25) can be calculated with a discrete sum of s f (n ij ) with quadrature weights w ij . In principle, other sampling theorems could be used (e.g., Refs. [62,63]), but note that at fixed band limit this sampling scheme requires the fewest samples [61]. pixelization schemes not based on sampling theorems could also be adopted (e.g., HEALPix [64]), in which case larger errors due to approximate numerical quadrature are expected. On SO (3), a function f with angular band limit L and azimuthal band limit N satisfies f mn = 0, ∀ ≥ L, |n| > N . The sampling theorem of Ref. [61] was extended to SO(3) in Ref. [48]: the sampling points, also equiangular and separable, are denoted by ρ hij = (α i , β j , γ h ) with i = 0, . . . , L − 1, j = 0, . . . , 2L − 1, and h = 0, . . . , 2N − 1. The quadrature weights are denoted by w hij and can be used to compute the Wigner coefficients of Eq. (27) exactly. The theoretically exact spin Fourier-Laguerre transform for bandlimited signals on which are the discrete versions of Eqs. (29) and (30), respectively. The theoretically exact Wigner-Laguerre transform for bandlimited signals on SO(3)×R + reads which are the discrete versions of Eqs. (31) and (32), respectively. Although we express the forward transforms of Eqs. (33) and (35) using discrete quadrature above, in practice we do not compute these expressions explicitly but rather compute them using fast algorithms, drawing on Refs. [48,61], as discussed below. The sums in these equations are finite and adjusted to the band limits (for the mnp indices) and to the sampling nodes (for the hijk indices). In what follows we will omit the pixel indices and bounds of the summations for conciseness. We will also use the integral (continuous) forms of the transforms for clarity, since the sampling theorem guarantees that these can be evaluated exactly for bandlimited signals (which is typically the case in practical applications).
As described above, any bandlimited signal with radial and angular band limits P and L, respectively, can be represented exactly by ∼ 2P L 2 samples on S 2 × R + thanks to the combination of the sampling theorem for the pixelization of Ref. [61] and the spherical Laguerre sampling theorem [39]. Similarly, signals on the rotation group, with the same band limits P and L, and an additional azimuthal band limit N , are captured by ∼ 4P N L 2 samples on SO(3)×R + , thanks also to the corresponding sampling theorem on the rotation group [48].
In terms of computational complexity, the separability of the angular and radial components in the Fourier-Laguerre and Wigner-Laguerre transforms is an essential property. Fast algorithms exist to compute transforms on the sphere and the rotation group [47,48,61,[65][66][67]. In particular, we use the implementations detailed in Refs. [48,61], which sets the complexity of the spin spherical harmonics and Wigner transforms to O(L 3 ) and O(N L 3 ), respectively, by exploiting fast Fourier transforms on the ring torus. The spherical Laguerre transform simply scales as O(P ), and the quadrature weights and basis functions are computed through recurrence relations [39]. Therefore, the spin Fourier-Laguerre and Wigner-Laguerre transforms scale as O(P L 3 ) and O(N P L 3 ), respectively. Note that the spin number s is a parameter that does not change the complexity or the accuracy of any of the transforms, which is often not the case for alternative approaches.

D. Connection to the Fourier-Bessel transform
In contrast to the Fourier-Laguerre construction, the standard nonseparable basis for 3D spin functions s f (n, r) ∈ L 2 (S 2 ×R + ) (such as weak lensing observables) is the combination of spherical harmonics with spherical Bessel functions j . These define the Fourier-Bessel transform, already used in Sec. II, with forward and inverse Unlike in the Fourier-Laguerre case, the Fourier-Bessel basis functions and transform are not separable: the angular modes are coupled with the radial modes k through the spherical Bessel function j (kr). This is because the basis functions are solutions of the (isotropic) Laplacian operator in spherical coordinates, with eigenvalues −k 2 . It is important to note that the Fourier-Bessel transform is notoriously difficult to evaluate for generic functions because it does not admit a sampling theorem [39,68]. Therefore, one must appeal to an approximate quadrature rule for the radial integrals (see e.g., Refs. [6,38,69]), which is challenging because the spherical Bessel functions have infinite support and oscillate rapidly.
Exact analytical formula only exist for simple forms of signals, but the first exact formula to compute the Fourier-Bessel transform of a wide class of signals was presented in Ref. [39]. The starting point of this approach is to decompose the spherical Bessel functions in the spherical Laguerre basis, i.e., The Fourier-Bessel coefficients of s f can then be expressed in terms of its Fourier-Laguerre coefficients: This sum is finite if s f is bandlimited in spherical Laguerre space, and s f m,p is also evaluated exactly thanks to the exact quadrature rules of the Fourier-Laguerre transform. Thus, Eq. (41) provides a way to compute the Fourier-Bessel transform of signals described in Fourier-Laguerre space exactly. It was also shown in Ref. [39] that J ,p (k) admits an analytical form, involving the moments of j (kr). Although J ,p (k) can still be challenging to evaluate numerically for high , it does not depend on the signal s f under consideration and can be tabulated. We will see that J ,p (k) appears in the treatment of weak lensing observables, natively expressed in the Fourier-Bessel basis.

E. Rotation, translation, and convolutions
We now construct rotation and convolution operators on S 2 × R + and SO(3)×R + , generalizing the operators of the scalar Fourier-Laguerre transform [39,68]. In general, such operators are essential for constructing a meaningful wavelet transform, where the wavelet coefficients of a signal are its convolution with the wavelets. Since wavelets have compact support in real and frequency space, this approach allows the transform to extract well defined scales and directions.
We first recall that the rotation of a spin function on the sphere s f ∈ L 2 (S 2 ) is naturally defined by where x is the Cartesian vector corresponding to n and R −1 ρ is the 3D rotation matrix corresponding to the rotation operator R ρ . In harmonic space, the rotation conveniently reduces to [46] ( Furthermore, we define the translation of a radial function f ∈ R + in spherical Laguerre space as [39,68] ( This is analogous to the harmonic representation of the translation operator for functions on the infinite (Cartesian) line R (but not analogous to the real space expression of the translation operator). In addition, such a translation operator can be expressed in terms of convolution with a shifted Dirac delta function, again in analog to translation on R, as detailed in Ref. [68]. We now consider functions on S 2 × R + , and introduce a directional convolution operator such that its action on f, g yields a function on SO(3)×R + defined by We also introduce an axisymmetric convolution operator (a special case of ) such that its action between f and an axisymmetric function h (satisfying R (0,0,γ) h = h, ∀γ) yields a function on S 2 × R + defined by where R n is the rotation operator for axisymmetric functions, i.e., R n = R (φ ,θ ,0) . In what follows, R ρ and R n refer to directional and axisymmetric rotations, respectively, since ρ and n denote angles on SO(3) and S 2 , respectively. More detailed discussions about rotation, translation, and convolution operators on the sphere and the radial line can be found in Refs. [46] and [39,68], respectively. All operators can be evaluated exactly via their harmonic representations and all transforms from pixel to harmonic space and conversely are exact thanks to the sampling theorems adopted.

IV. SPIN FLAGLET TRANSFORM
We now extend the scalar, axisymmetric flaglet transform of Ref. [39] to analyze signals of arbitrary spin and to probe their features in the angular direction. This new separable flaglet transform is defined in Sec. IV A, where we present the core equations and rely on the previous translation and convolution operators. However, to highlight the generic properties of the transform we do not specify the details of the flaglet kernels s Ψ ij and s Φ until Sec. IV B, where we also present our efficient implementation using sampling theorems and multiresolution algorithms.

A. Transform
The flaglet coefficients of a spin signal s f ∈ L 2 (S 2 × R + ), denoted by W sΨ ij s f , are defined on SO(3)×R + as the directional convolution of s f with the flaglets s Ψ ij ∈ L 2 (S 2 × R + ) where i = I 0 , . . . , I and j = J 0 , . . . , J denote the angular and radial scales, respectively. Flaglets are localized in scale, position and angular orientation, and designed to capture the directional, high-frequency content of the signal of interest. A scaling function s Φ ∈ L 2 (S 2 × R + ) is introduced to capture the low frequency content of the signal in the scaling coefficients W s Φ s f . In this work we analyze the low frequency part of the signal using axisymmetric flaglets because its directional content is typically of lower interest. In the context of weak-lensing, largest scales (low k in C φφ ) have far fewer modes, contain less information (due to the reduced amplitude of C φφ ), and are more difficult to deal with due to observational masks and other considerations mentioned above. However, extending the formalism presented here to support directional scaling functions is possible and straightforward.
We define W s Φ sf ∈ L 2 (S 2 × R + ) as the axisymmetric convolution of s f with the axisymmetric scaling function s Φ Since s Φ is axisymmetric, the Fourier-Laguerre coefficients s Φ m,p are nonzero for m = 0 only. In what follows we will use a compressed notation for all sums over scales i and j, assumed to run over I 0 , . . . , I and J 0 , . . . , J, respectively. These bounds, as well as the flaglets and scaling functions s Ψ ij and s Φ, are consistently defined in the next section in order to achieve exact reconstruction of the input (bandlimited) signal.
The wavelet transform makes direct use of the rotation, translation and convolution operators defined previously. Hence, the forward flaglet transform (or analysis step, equivalent to Eq. (47)) can also be expressed in Wigner-Laguerre space (on SO(3)×R + ) in a straightforward manner by The scaling coefficients (on S 2 × R + ) in Fourier-Laguerre space read The reconstruction of the original signal is achieved through the inverse flaglet transform (synthesis step) or in Fourier-Laguerre space In order for the transform to be exact, i.e., for the wavelet coefficients to capture all the information content of s f , the flaglets and scaling function must be defined to satisfy Eqs. (49), (50), and (52). In other words, they must be chosen such that their Fourier-Laguerre coefficients satisfy the admissibility condition  Figs. 3 and 4. Note that we do not show the hybrid kernel η λν , included in the scaling function to satisfy the admissibility condition at low p and , see Eqs. (53) and (61). We also do not include the directional component ζ m in this figure, which is nevertheless shown in real space in Fig. 3.

B. Wavelet construction
So far we have constructed a generic flaglet transform relying on the properties of the Fourier-Laguerre space. To uniquely characterise this transform, we need to specify the scaling function s Φ, the flaglets s Ψ ij , and the bounds I 0 , I, J 0 , J to satisfy Eq. (53). We follow the construction of scalar flaglets [39] and scale-discretized wavelets [46,70]. Other types of wavelets could be used; we opted for scale-discretized wavelets because they exhibit good localization properties [71] and can be easily extended to probe directional features [46,70]. They exhibit a compact representation in harmonic space that also allows several computational improvements, which are essential to make a 3D spin directional transform tractable. We use tiling functions κ and η defined as where which is unity for t < 1/λ, zero for t > 1, and is smoothly decreasing from unity to zero for t ∈ [1/λ, 1]. In these equations, λ is the wavelet dilation parameter and characterises the density of the tiling. In what follows, λ will refer to the wavelet tiling in the angular direction, while a second parameter ν is introduced to tile the radial direction (with functions κ ν and η ν constructed in the same manner). Finally, the remaining free function s λ (t) (also defined for ν) must have compact support in [ 1 λ , 1] and is defined as The tiling of the radial and angular harmonic spaces constructed with these functions is shown in Fig. 2 for λ = ν = 3, I 0 = J 0 = 2, L = P = 243.
In addition, as for the scalar axisymmetric flaglet transform [39], a hybrid tiling function is needed in order to construct a suitable scaling function Φ and satisfy Eq. (53), With the tiling functions κ and η defined in Eq. (54) and Eq. (55), we construct the flaglets in Fourier-Laguerre space as The extra function ζ controls the directionality component [47,48,70] and is parameterised by an azimuthal band limit N , such that ζ m = 0, ∀ , m with |m| ≥ N . It is defined in harmonic space as to satisfy various azimuthal symmetries and normalized such that m ζ m = 1, for all values of for which ζ m are nonzero for at least one value of m (see, e.g., , [71]). In this expression, a = 1 for N − 1 odd and a = i for N − 1 even, Finally, the axisymmetric scaling function is needed to capture the information from the regions ≤ λ I0 and p ≤ ν J0 , which are not probed by the flaglets defined above. Thus, it is constructed to satisfy Eq. (53) and reads These definitions entirely characterise the spin, directional flaglet transform. The free parameters of the transform are λ, ν, N, I 0 , and J 0 . λ and ν define the support of the flaglets in Fourier-Laguerre space and therefore the scales and features the flaglet coefficients extract, as shown in Fig. 2. I and J are the maximum angular and radial flaglet scales captured, respectively, and are fixed by the radial and angular band limits of the signal: I = log λ (L − 1) and J = log ν (P − 1) , to satisfy Eq. (53). N is used to specify the number of directions to probe. I 0 and J 0 are the first angular and radial scales of interest (larger scales being captured by the scaling function) and can be freely chosen provided they satisfy 0 ≤ I 0 ≤ I and 0 ≤ J 0 ≤ J.
The flaglets (and scaling function) tile the angular and radial frequency domains and p, as illustrated in Fig. 2, while satisfying the admissibility condition of Eq. (53). Wavelets are well localized simultaneously in the spatial domain, both in position and orientation, and the harmonic domain. It is shown in Ref. [71] that scale-discretized wavelets on the sphere exhibit excellent concentration properties, both in the scalar setting and in the spin setting. They are also steerable [47,48,70]. Flaglets naturally inherit these properties. In particular, as shown in Ref. [68] for the scalar setting, the flaglet transform forms a tight Parseval frame, i.e., the norm of the input signal is conserved.
Figs. 3 and 4 show spin-0 and spin-2 flaglets resulting from our construction, with parameters λ = ν = 3 and I 0 = J 0 = 2. In contrast to the Fourier-Laguerre basis functions shown in Fig. 1, which are not localized in real space (but are delta functions in harmonic space), all flaglets have good spatial localization properties in both radial and angular dimensions. They are also localized in Fourier-Laguerre space by construction, as shown in Fig. 2 and Eq. (59). Note that the particular i = 2, 3 and j = 2, 3 scales shown in Figs. 3 and 4 are highlighted as thicker lines in Fig. 2. Fig. 3 shows spin-0 flaglets (real part only), comparing the N = 2, 3 flaglets to the axisymmetric N = 1 case (first shown in Ref. [39]). The angular components have even and odd symmetry for N = 2, 3, respectively, by the construction of ζ [46]. Fig. 4 shows the real and imaginary parts and the complex modulus of spin-2 axisymmetric (N = 1) flaglets.

C. Multiresolution algorithm and implementation
As described previously, all the integrals and convolutions related to the Fourier-Laguerre and Wigner-Laguerre transforms can be computed exactly via their harmonic space representations and by appealing to sampling theorems to analyze or reconstruct signals in pixel space. Let us consider an input signal s f ∈ L 2 (S 2 × R + ) with radial and angular band limits P and L, respectively. As before, such a signal can be represented exactly by ∼ 2P L 2 samples and P L 2 spherical harmonic coefficients. The full algorithm to perform the flaglet transform is described by the pipeline below. The first and the third operations are spin Fourier-Laguerre transforms, while the second is the Spin s = 0 flaglets with λ = ν = 3, I 0 = J 0 = 2 and N = 3 We have implemented the exact and efficient algorithm described previously to compute the spin directional flaglet transform in the existing flag and flaglet codes (http://www.flaglets.org). The latter relies on the ssht (http://www.spinsht.org), s2let (http://www.s2let.org) and so3 (http://www.sothree.org) codes for the angular transforms and sampling theorems. The fftw (http://www.fftw.org) code is used to compute Fourier transforms. The core algorithms of the Fourier-Laguerre, Wigner-Laguerre and flaglet transforms are implemented in C and are able to handle large band limits and billions of samples on S 2 × R + and SO(3)×R + . Interfaces to the core functions as well as convenient data manipulation and plotting routines are provided in Matlab and Python.

V. APPLICATION TO 3D WEAK LENSING
In previous sections we developed the spin flaglet formalism to extract content localized in both space and frequency from 3D signals of arbitrary spin, which has the added benefit that radial and angular modes are separable. This is a general framework that can be applied to arbitrary 3D spin signals. We now apply this approach to represent 3D cosmic shear and other weak lensing observables, and relate the covariance of their flaglet decompositions to the lensing power spectrum. In particular, we highlight how one can exploit the properties of flaglets (separability and simultaneous localization in frequency and pixel space) to deal with partial sky coverage and small-scale modeling uncertainties. We also examine the approximations commonly adopted in cosmic shear analyzes and apply them to the flaglet covariance of weak lensing observables. Note that we focus on using axisymmetric wavelets to increase simplicity and readability of the final estimator. However, the latter is straightforwardly extended to directional flaglets, and this does not affect any of the advantages and properties discussed below.

A. 3D cosmic shear
We relate the covariance of the flaglet representation of the cosmic shear field to the lensing power spectrum. As shown in Sec. II, the 3D cosmic shear signal 2 γ is characterised by the lensing power spectrum C φφ (k, k ) defined in Fourier-Bessel space by where all relevant physical effects are assumed to be modelled (as discussed in Sec. II C). As in Sec. II, the brackets in this section denote ensemble average over cosmological realisations, rather than inner products (which are distinguished by a vertical bar between the functions of the inner product). We now consider the wavelet transform of 2 γ, using axisymmetric spin flaglets (i.e., N = 1). In this case all directional convolutions become axisymmetric convolutions . The flaglet coefficients live on S 2 × R + and read We consider the covariance between flaglet scales, Injecting the Fourier-Bessel decomposition of 2 γ into the integral form of the convolution of W 2Ψ ij 2γ (n, r), one obtains where 2 F ij m (k, n, r) is the wavelet transform of the Fourier-Bessel basis functions: Note that this formalism can be straightforwardly extended to directional wavelets: directional convolutions then replace axisymmetric convolutions, and the angles n, n ∈ S 2 become ρ, ρ ∈ SO(3) since the flaglet coefficients then live on SO(3)×R + . The previous expression can be simplified by appealing to the harmonic representation of the translation and rotation operators, which when applied to axisymmetric flaglets 2 Ψ ij (e.g., Refs. [46,72]) yield The functions 2 F ij m (k, n, r) then take the simple form As shown previously, the functions J ,p admit an analytical form [39] and can be tabulated. Furthermore, 2 H ij (k, r), and consequently 2 F ij m (k, n, r), is straightforward to compute since only a one dimensional summation over a finite number of p samples is required.
Further simplifications can be made by exploiting the spherical harmonic addition theorem: where P are the Legendre polynomials. By the additional theorem one can infer that the wavelet covariance only depends on the angle ∆θ, where n · n = cos(∆θ). Thus, C ij,i j (n, n , r, r ) becomes The angular dependence of the flaglet covariance depends only on the angular separation between pixels on the sky. This is expected since we are analyzing the 2-point fluctuations of the signal with axisymmetric flaglets and have assumed statistical homogeneity and isotropy. Eq. (73) is one of the key results of this paper. One can build a full likelihood analysis of 3D weak lensing observables and constrain cosmological models and parameters using spin flaglets (similar to that of Ref. [12]). On the left hand side of Eq. (73) is a quantity that can be computed from data using a wavelet decomposition (here using axisymmetric flaglets, but a similar expression can be obtained with directional flaglets). On the right hand side is the theoretical expectation of this quantity related to cosmological parameters through the 3D lensing power spectrum C φφ (k, k ). Even though C φφ (k, k ) mixes radial and angular modes, these are isolated into scales (i, j, i , j ) by the kernels 2 H ij (k, r). In practice, the flaglet covariance can be calculated for a pair of flaglet coefficients ij and i j and a pair of radii r and r by averaging over pairs of pixels separated by ∆θ. This approach naturally accounts for the 2+1D nature of the shear observables and takes advantage of the convenient connection between the Fourier-Laguerre and Fourier-Bessel transforms. It takes full advantage of the separability of the Fourier-Laguerre and flaglet transforms. Going through flaglet space allows one to exploit the excellent localization properties of flaglets in both pixel and frequency space, a property absent in all alternative approaches. Simultaneous localization in both pixel and frequency space allows one to cut or filter regions of the sky due to unobserved, unreliable or contaminated data, at the same time as cutting or filtering harmonic modes. For example, one may wish to cut or filter small-scale harmonic modes where physical modeling is less certain, e.g., high-k (and k ) modes from C φφ (k, k ). This is possible by studying the kernels 2 H ij (k, r) and only considering the flaglet scales ij, i j that probe the physical scales of interest in C φφ (k, k ). Thus, the spin flaglet 3D weak lensing approach is a flexible way to avoid the complications of correlation functions and Fourier-Bessel representations, which are typical of standard approaches.

B. Generalization to other weak lensing quantities
As shown in Eq. (6), spin-2 shear distortions are not the only effect of matter fluctuations on the observed properties of galaxies [49]. Other observables are expected to be measured at much lower signal-to-noise ratio than the shear distortion but they nevertheless contain important cosmological information. Thus, we generalize the flaglet covariance derived previously to other lensing distortions of different spin. For this purpose we introduce a general scaling factor so that we can write s χ m (k) = M ,s φ m (k), with s χ defined in Eq. (6). In this unified approach, we have C ij,i j (n, n , r, r ) = 2 π m (M ,s ) 2 with s F ij m (k, n, r) generalizing Eq. (70) to other values of spin. As before, the extension to directional wavelets is straightforward. Exploiting the spherical harmonic addition property, again, leads to the simplification C ij,i j (∆θ, r, r ) = 2 π (M ,s ) 2 with s H ij (k, r) defined as previously. This expression generalizes the 3D cosmic shear power spectrum to the case of lensing distortions of other spin. As before, this quantity can be calculated from data by averaging over pairs of pixels separated by ∆θ in bins of ij, i j , r, r . The theoretical prediction through the lensing potential, the Legendre polynomials and the s H kernels is readily computable. Finally, this formalism can also be extended to support cross-terms between the various spin components, which would be required for a combined analysis.

C. Approximations
We have related the covariance of weak lensing observables in flaglet space to the underlying 3D lensing power spectrum in the general setting. Several approximations are used in the literature to make the computation of shear and other lensing quantities more approachable. Some of these are lossless-or nearly so-while others result in a loss of information with respect to the general case. The most commonly used approximation is to use cosmic shear 'tomography', in which flat-sky and Limber approximations are typically also made. We apply these approximations below to the flaglet covariance of spin lensing quantities, where appropriate, and show that they result in representations that are more readily computable but are a lossy representation of the 3D lensing power spectrum.

Flat-sky approximation
Weak gravitational lensing of galaxy images has the strongest signal-to-noise on the scales of groups or clusters (e.g., Ref. [55]) and therefore there is relatively little signal in the power spectrum on large angular scales. In addition, for surveys with relatively small sky coverage and/or surveys with several small observational fields, effects due to the spherical geometry of the setting are expected to have a relatively small impact. However, for future surveys (e.g., Euclid [73,74] and LSST [75]) with observations over increasingly greater coverages of the sky, flat-sky approximations will become increasingly less accurate.
The flat-sky approximation assumes that the angular extent of the observational field is small and hence the geometry of the angular component is assumed to be planar (i.e., Euclidean). In this setting the spherical harmonics are approximated by the product of Bessel functions (of the first kind) and complex exponentials [7,76]. In the case of Eqs. (73) and (76), however, the sum over the product of spherical harmonic functions reduces to a Legendre polynomial (due to the spherical harmonic addition theorem), which in the flat-sky case would be approximated by a zeroth order Bessel function. However, given the computational ease of computing Legendre polynomials the full-sky case is readily computable and so a flat-sky approximation is not required. This differs to the computation of C φφ (k, k ), where a flat-sky approximation can reduce computational time considerably and has so-far been assumed in current applications to existing data [11,12], where sky coverage is relatively small.

Limber approximation
The Limber approximation [77] assumes that the evolution of radial modes is small over the survey volume under consideration. Under this approximation the spherical Bessel functions may be approximated by Dirac delta functions: as shown in Refs. [8,59], and rederived in Appendix B for the spherical Bessel convention adopted herein. The use of this approximation simplifies matters considerably: an 8-nested integral equation for C φφ (k, k ) becomes a single integral (see, e.g., Refs. [9,12]). However, even though this approximation primarily affects large scales 100, it is severe and can cause errors in the inferred cosmological parameters of tens of percent [9]. Exploiting , r 2 H i j * + 1/2 r , r , where the H ij kernels are those given in the Limber-approximated case of Eq. (79). The final approximation for tomography with binning in redshift is C ijj k (∆θ, r bin (z a ), r bin (z a )) R + dr R + dr C ij,i j (∆θ, r, r ) W (r, r , r bin (z a ), r bin (z a )), where the weight function W describes if a galaxy is in the bin-pair combination (r bin (z a ), r bin (z a )) or not. Bin weight functions are usually designed as top-hat functions in spectroscopic redshift z a and corrections are applied for 'leakage' due to photometric redshifts (see, e.g., Ref. [9]). The methods presented here can also be extended to deal with redshift uncertainties and more optimal weight functions.

VI. CONCLUSIONS
We have developed a novel approach to analyze weak gravitational lensing observations using 3D wavelets (flaglets) and the Fourier-Laguerre transform. We have shown that the covariance of the flaglet coefficients of weak lensing observables (such as cosmic shear, size distortions and flexion) can be directly related to the 3D lensing power spectra, exploiting the analytical connection between the Fourier-Bessel and Fourier-Laguerre transforms. Thanks to the separability of their radial and angular components, the Fourier-Laguerre and flaglet transforms can take advantage of existing fast algorithms and sampling theorems developed for the radial line, the sphere, and the rotation group. Furthermore, this approach is suited to dealing with the 2+1D nature of weak lensing observations and cosmological data sets in general. For example, typical data complications and systematic uncertainties are given on the sky and in the radial direction separately. The excellent localization properties of flaglets in both pixel and frequency space can be exploited to deal with these effects and alleviate the complications of weak lensing analyzes (e.g., complicated sky masks and uncertainties of the small-scale modeling) without complicating the estimation of the flaglet coefficients and covariance. Therefore, the method introduced here is a powerful alternative to common real and harmonic-space approaches where these effects are difficult to deal with [11,12]. Lensing observables and covariances are more easily modelled in the Fourier-Bessel basis, where unreliable or unmodelled small scales can also be filtered out. However, the mode-mixing due to partial sky coverage is difficult to handle in Fourier-Bessel space. By contrast, real space correlation function measurements can naturally deal with complicated survey window functions, but filtering scales and modeling covariances is then significantly more difficult. The dual localization of wavelets in pixel and harmonic space gives them the advantages of both approaches. In future work we will demonstrate the ability of this new approach to robustly extract cosmological information by applying it to simulations and real observations of weak lensing signals.
We have updated the publicly available flag and flaglet (http://www.flaglets.org) codes to compute the spin Fourier-Laguerre and flaglet transforms. These rely on the fast algorithms implemented in the following codes, also publicly available: s2let (http://www.s2let.org) for the spin directional spherical wavelets, ssht (http:// www.spinsht.org) for the spherical harmonics transform, so3 (http://www.sothree.org) for the Wigner transform, and fftw (http://www.fftw.org) for the Fourier transform. In all these codes, the core routines are implemented in C, and we have provided numerous wrappers and interfaces in Python, Idl, and Matlab, to call the core routines, and to manipulate and visualize data sets. As highlighted in Ref. [39] for the scalar setting, these codes can deal with large band limits and large 2D and 3D data sets, which can be pixelized into billions of samples and manipulated with no loss of information.
The orthogonality and completeness of the spin spherical harmonic functions s Y m read respectively. The spin spherical harmonics are the canonical (complete) basis for the representation of functions defined on the sphere S 2 , while accounting for the symmetry of spin-s signals. Spin raising and lowering operators, ð andð respectively, increment and decrement the spin order of a spin-s function and read [51][52][53][54] ð = − sin s θ ∂ ∂θ respectively. Using these operators, spin-s spherical harmonics s Y m can be constructed from scalar (spin-0) harmonics by We assumed the Condon-Shortley phase convention and define the factor respectively. The Wigner functions provide a complete basis for the representation of functions defined on the rotation group SO (3). Note that we adopt D * mn as basis functions since this simplifies connections to wavelet transforms.