© The Authors 2024

1 Introduction

Scattering transforms are a recently developed class of summary statistics for the study of non-Gaussian processes (Mallat 2012; Bruna & Mallat 2013). These statistics, which are built from successive wavelet convolutions and pointwise non-linearities such as a modulus, are inspired by neural networks but do not require any training steps. Introduced recently in astrophysics (Allys et al. 2019, 2020), scattering transforms have since demonstrated their ability to characterise highly nonGaussian processes, for instance for parameter estimation and classification tasks in fields as varied as the interstellar medium (Regaldo-Saint Blancard et al. 2020; Saydjari et al. 2021; Lei & Clark 2023), the large-scale structures (LSSs) of the Universe (Allys et al. 2020; Cheng et al. 2020; Cheng & Ménard 2021; Valogiannis & Dvorkin 2022a,b), or the epoch of reionisation (Greig et al. 2022; Hothi et al. 2024).

Another feature of scattering transforms is that they allow one to build very efficient generative models of physical fields, in a maximum entropy framework (Bruna & Mallat 2019). This allows one to sample new approximate realisations of a given process relying only on its scattering transform statistics, which can be estimated from a small amount of data, sometime even a single example image (Allys et al. 2020; Régaldo-Saint Blancard et al. 2023; Price et al. 2023; Cheng et al. 2024b). One application of such generative models is for training data augmentation for machine learning applications. For instance, it has been shown in Jeffrey et al. (2022), with simulated data, that such scattering transform models constructed from a single polarised microwave dust foreground patch can be sufficient to separate primordial B-modes in the cosmic microwave background (CMB) polarisation from this dust emission, even in an artificially challenging mono-frequency approach. Moreover, the framework on which these generative models have been constructed has led to the development of new statistical component separation algorithms, which have for instance been successfully applied to astrophysical data (Regaldo-Saint Blancard et al. 2021; Delouis et al. 2022; Auclair et al. 2024) and seismic signals (Siahkoohi et al. 2023a,b).

While these promising scattering transform generative models have mainly been developed for 2D planar data, the adaptation of these tools to spherical data is necessary for cosmological analysis, especially for the next generation of full sky surveys, such as LiteBIRD, the Lite (Light) satellite for the studies of B-mode polarisation and Inflation from cosmic background Radiation Detection (LiteBIRD Collaboration 2023), the Euclid space telescope (Laureijs et al. 2011) or the Vera C. Rubin Observatory (LSST Science Collaboration 2009). The adaptation of a first generation of scattering transforms to spherical signals was already introduced in McEwen et al. (2022) and has been used as a form of dimensionality reduction for other machine learning purposes. In this paper, we extend state-of-the-art scattering transforms (Morel et al. 2023; Cheng et al. 2024b) named the scattering covariances (SCs) to spherical fields.

The extension of scattering transforms to spherical data raises certain difficulties, namely, the definition of a directional spherical convolution with oriented filters such as wavelets (McEwen et al. 2015b, 2018), and the transposition of the planar translations, which appear in certain scattering transform representations. As a first step, we restricted ourselves to fields with statistically homogeneous fields with properties that do not depend on the position on the sphere. The generalisation beyond statistical homogeneity will be presented in a forthcoming paper. This naturally led us to cosmological fields, such as a weak lensing field from the LSSs of the Universe, and maps of the thermal Sunyaev-Zeldovitch (tSZ) effect of the CMB. We also consider a map of the Venus surface. In this paper, for all these spherical data, scattering transform generative models were constructed and validated from one single full-sky image.

In Sect. 2, we present the extension of the SC statistics to spherical fields. THEN, in Sect. 3 we present SC-based generative models and discuss their numerical implementation. Finally, in Sect. 4, we present the results obtained with these models for the four non-Gaussian spherical fields studied. Our conclusions are presented in Sect. 5.

2 Scattering covariance on the sphere

The SCs, or scattering spectra, were previously introduced in Morel et al. (2023) for 1D data and in Cheng et al. (2024b) for 2D planar maps. In this paper, we extend these statistics to spherical maps. This section introduces sampling schemes, directional convolutions, and wavelet transforms on the sphere, after which we define the SC statistics.

2.1 Sampling of spherical maps

A spherical field can be represented by its spherical harmonic transform, which is the spherical equivalent of the Fourier transform for planar maps. In the following, we worked with the usual spherical coordinates ω = (θ, φ), with co-latitude θ and longitude φ. With these coordinates, the spherical harmonic coefficients I_ℓm of a spherical field I(ω) defined over the sphere 𝕊² correspond to the projection onto the spherical harmonics Y_ℓm (ω): $$I_{\ell m}={\displaystyle {\int }_{{\mathbb{S}}^2}I}(\omega )Y_{\ell m}^{*}(\omega )\text{d}\Omega (\mathit{\omega }),$$(1)

where dΩ(ω) = sin θdθdφ is the solid angle element. The field can then be represented by its harmonic expansion, given by $$I(\omega )={\displaystyle \sum _{\ell =0}^{L-1}{\displaystyle \sum _{m=-\ell }^{\ell }{I}_{\ell m}}}{Y}_{\ell m}(\mathit{\omega }).$$(2)

The ℓ index is called the multipole and is inversely proportional to angular scales on the sky, while the order m at a given ℓ goes from −ℓ to ℓ and captures the anisotropic component of I(ω). The maximum value of ℓ considered, ℓ_max = L − 1, determines the smallest scale in the transform. For a real field, the coefficients satisfy the relation $$I_{\ell m}={(-1)}^m{I}_{\ell -m}^{*}.$$(3)

The numerical computation of the forward and inverse spherical harmonic transform depends on the sampling in pixel space of the spherical map, for which different choices can be made. For example, in cosmology, the community often adopts Healpix sampling (Górski et al. 2005), in which all pixels have the same area, which can be an advantage in practice. With this sampling, the map resolution is given by the nside parameter, the number of pixels being equal to 12 × nside². However, with this sampling, the spherical harmonic transform (as well as the Wigner transform defined below) is not accurate and needs to be refined iteratively. An alternative is instead to use an equiangular sampling, such as that defined in McEwen & Wiaux (2011), for which these transforms can be computed exactly (to machine precision). With this sampling, abbreviated by MW, the angular dimensions of all pixels are the same, and maps are stored as 2D arrays of shape (N_θ, N_φ) = (L, 2L − 1). In this paper, the SC statistics computations support various sampling schemes, including Healpix, MW, and others, while internal calculations typically adopt sampling schemes, such as MW, that afford exact spherical transforms for improved accuracy.

Another operation required, which is also sampling dependent, is the average on the sphere, defined as: $$\langle I(\omega )\rangle =\frac{1}{4\pi }{\displaystyle {\int }_{{\mathbb{S}}^2}I}(\omega )\text{d}\Omega (\omega ).$$(4)

In pixel space, this computation corresponds to $$\langle I(\omega )\rangle =\frac{1}{4\pi }{\displaystyle \sum _{p}I}({\omega }_p)\delta {\Omega }_p,$$(5)

where the sum is done over all pixels p, whose angular positions are noted ω_p. For approximate quadrature, δΩ_p can simply represent solid angle at pixel p. Alternatively, some sampling schemes exhibit exact quadrature (McEwen & Wiaux 2011), in which case δΩ_p denotes quadrature weights. When I_ℓm has been computed, one directly has $$\langle I(\omega )\rangle =\frac{1}{2\sqrt{\pi }}{I}_{00}.$$(6)

2.2 Wavelet transform on the sphere

The SC statistics are computed from wavelet transforms, which are obtained by convolving an initial map with a set of wavelet filters, where each filter extracts the local information at a particular scale and direction. Wavelet filters need to be localised both in pixel and harmonic space. In this work, we followed McEwen et al. (2015b, 2018) and defined the wavelets in harmonic space in separated form as $${\Psi }_{\ell m}^j=\sqrt{\frac{2\ell +1}{8{\pi }^2}}{\kappa }_{\ell }^j{\zeta }_{\ell m},$$(7)

in order to control their angular and directional localisation properties separately, respectively, by kernel ${\kappa }_{\ell }^j$, with filter scale j, and directional component ζ_ℓm . For the explicit definition of these two functions, one can refer to, for instance, McEwen et al. (2015b). The wavelets were designed to satisfy excellent spatial localisation and asymptotic non-correlation properties (McEwen et al. 2018). Moreover, their directional structure in m was set to zero for even/odd m to enforce odd/even symmetry in φ, resulting in odd/even symmetry in φ for N − 1 odd/even (McEwen et al. 2018).

In Figure 1, the left panel shows the Ψ_ℓm coefficients of one filter at a specific scale j, and the right panel, a cut at m = 0 of the full filter set. We note that with our convention, the j scales are ordered with ℓ multipoles, meaning that when j increases, the corresponding angular scale decreases (i.e. ℓ increases). Filters are maximum at ℓ ≃ η^j, with support within ℓ є [η^(j−1),η^{(j +1)}], where η defines the wavelet dilation parameter. In this paper we use dyadic scaling, corresponding to η = 2. In the following, when performing a convolution with a filter set, the range of scales probed by the wavelets is given by J_min ≤ j ≤ J_max where J_min ≥ 0 and $J_{\max }=\text{ceil}\left(\frac{\log (L-1)}{\log \eta }\right)$. The number of scales is given by J = J_max − J_min + 1. The angular resolution of the wavelets is parameterised by an integer N, allowing for the sample of 2N − 1 directions (see below)¹.

Wavelet transforms were computed from convolutions between the field under study and this set of wavelets. Various convolutions on the sphere can be considered (Roddy & McEwen 2021). In this work, we followed the standard directional convolution formalism presented in, for instance, McEwen et al. (2015b). These convolutions produce a set of spherical maps filtered at different scales (labelled by j) and orientations (labelled by γ).

The directional convolution I ★ Ψ^j of a field I with a wavelet Ψ^j consists in applying a rotation by a set ρ = (α,β,γ) of Euler angles of the wavelet Ψ^j initially located at the north pole, before computing an inner product between the wavelet and the field I: $$(I\star {\Psi }^j)(\mathit{\rho })={\displaystyle {\int }_{\Omega }I}(\mathit{\omega }){[R_{\rho }{\Psi }^j(\mathit{\omega })]}^{*}\text{d}\Omega ,$$(8)

where R_ρ is the rotation by Euler angles ρ, and ∗ stands for complex conjugation. From (I ★ Ψ^j)(ρ), we can identify (β, α) with the spherical coordinates ω = (θ, φ) and γ to the orientation that is probed in the convolution. In this way, we obtain oriented wavelet coefficients, with the shorthand notation $$(I\star {\Psi }^{j,\gamma })(\mathit{\omega })\equiv (I\star {\Psi }^j)(\alpha =\phi ,\beta =\theta ,\gamma ).$$(9)

While (I ★ Ψ^j,γ) is a convenient notation, which also matches with previous work, we emphasise that there is no Ψ^j,γ oriented wavelet by itself.

In practice, the directional convolution can be computed accurately and efficiently in Wigner space, which is the Fourier space associated with 3D rotations described by Euler angles. In this space, the directional convolution between a field I and a wavelet Ψ^j yields (McEwen et al. 2013, 2015b) $${(I\star {\Psi }^j)}_{mn}^{\ell }=\frac{8{\pi }^2}{2\ell +1}{I}_{\ell m}{\Psi }_{\ell n}^{j*},$$(10)

where I_ℓm and ${\Psi }_{\ell n}^j$ are the spherical harmonic coefficients of I and Ψ^j, respectively, and ${({\Psi }^j\star I)}_{mn}^{\ell }$ are the Wigner coefficients of the convolved field, that is, the Fourier representation of the directional wavelet coefficients defined over Euler angles ρ = (α, β, γ). To return to the spatial domain, we computed an inverse Wigner transform, defined as $$(I\star {\Psi }^{j,\gamma })(\mathit{\omega })\equiv (I\star {\Psi }^j)(\mathit{\rho })={\displaystyle \sum _{\ell =0}^L\frac{2\ell +1}{8{\pi }^2}}{\displaystyle \sum _{m,n=-\ell }^{\ell }{(I\star {\Psi }^j)}_{mn}^{\ell }}{D}_{mn}^{\ell *}(\mathit{\rho }),$$(11)

where $D_{mn}^{\ell }(\mathit{\rho })$ are the Wigner-D matrices (Varshalovich et al. 1988). Fast (inverse) Wigner transform algorithms can then be leveraged for efficient computation (McEwen et al. 2015a; Price & McEwen 2024). By computing the wavelet transform through harmonic space as described, pixelisation artefacts are avoided. Although the wavelets satisfy an admissibility condition such that the field can be recovered exactly from its wavelet coefficients (McEwen et al. 2015b), we are only concerned with the forward wavelet transform in this work. In the following, we also grouped the scale and orientation under a single index λ = (j, γ), writing I ★ Ψ^λ for the wavelet transform of the field I at a given oriented scale λ.

While computing convolutions through harmonic representations is highly accurate, it involves moderate computational cost since generalised Fourier transforms on the sphere and space of rotations must be computed (albeit using fast algorithms). An alternative would be to compute the convolutions in pixel space as done in Delouis et al. (2022). However this can introduce pixelisation artefacts. A future avenue to consider is hybrid discrete-continuous approaches, as they have been shown to be highly computationally efficient while also avoiding discretisation artefacts (Ocampo et al. 2023).

In order to optimise our numerical implementation, in particular the memory usage, all convolutions were computed in a multi-scale framework, where the map resolution is tuned to the scale at which the convolution is made. See Leistedt et al. (2013) for more details.

Figure 2 illustrates orthographic projections of the directional spherical wavelets for two scales and three angles, viewed looking down from the North pole.

Fig. 1 Harmonic representation of the wavelet filters. Left: real part of the ${\Psi }_{\ell m}^{j=4}$ filter. Right: cut at m = 0 of the full filter set for j spanning from Jmin = 1 to Jmax = 8.

Fig. 2 Spatial representation of the wavelet filters. The circles contain orthographic projections of the directional spherical wavelets for two scales and three angles, viewed looking down from the North pole. In this case, we set N = 3 so that the angle γ takes 2N − 1 = 5 values between 0 and 360 deg, but only three are shown here.

2.3 Scattering covariance coefficients

Scattering transforms cover several types of summary statistics (see, e.g. Allys et al. (2019, 2020); Cheng et al. (2024b)). In this work, we consider the SCs, or scattering spectra, previously introduced in Morel et al. (2023) and Cheng et al. (2024b). We chose the SCs because they only rely on convolutions and not on translations as the wavelet phase harmonic statistics (Allys et al. 2020), which are difficult to univocally define on spherical maps.

The SC statistics characterise the power and sparsity at given scales, as well as interaction between different scales. They are built from successive applications of wavelet transforms and modulus operators, followed by average and covariance computations (assuming that we work with homogeneous processes). For a detailed introduction to these statistics, the reader is invited to refer to Cheng et al. (2024b). We considered two coefficients at a single scale j₁ and a single angle γ₁, that is, λ₁ = (j₁,γ₁): $$S_1^{{\lambda }_1}=\langle |I\star {\Psi }^{{\lambda }_1}|\rangle ,$$(12) $$S_2^{{\lambda }_1}=\langle |I\star {\Psi }^{{\lambda }_1}{|}^2\rangle ,$$(13)

and two coefficients that characterise the couplings between two and three oriented scales²: $$S_3^{{\lambda }_1,{\lambda }_2}=\text{Cov}\left[I\star {\Psi }^{{\lambda }_1},|I\star {\Psi }^{{\lambda }_2}|\star {\Psi }^{{\lambda }_1}\right],$$(14) $$S_4^{{\lambda }_1,{\lambda }_2,{\lambda }_3}=\text{Cov}\left[|I\star {\Psi }^{{\lambda }_3}|\star {\Psi }^{{\lambda }_1},|I\star {\Psi }^{{\lambda }_2}|\star {\Psi }^{{\lambda }_1}\right],$$(15)

where ⟨·⟩ corresponds to the mean over the sphere, defined in Sect. 2.1, and where the covariances are defined as Cov[XY] = ⟨XY^∗⟩ − ⟨X⟩⟨Y^∗⟩ for two complex fields X and Y. Note that in our case, the wavelet transforms are always zero mean, since the wavelets are mean-free. Since taking the modulus of a field mainly displaces its frequency support toward lower frequency (Zhang & Mallat 2021; McEwen et al. 2022), it is sufficient to consider terms for which j₁ ≤ j₂ ≤ j₃ . Moreover, we introduce the additional parameter δ_j, which corresponds to the maximum distance between pairs of scales whose interactions are characterised: that is, the calculation of S₃ and S₄ is restricted to pairs of scales such that j₂ − j₁ ≤ δ_j and j₃ − j₁ ≤ δ_j.

The dimension of S₁ and S₂ coefficients is JΘ, with J the number of scales and Θ = 2N − 1 the number of orientations. Regarding S₃ and S₄, their dimensions are approximately J²Θ² and J³Θ³ or JδjΘ² and JδjΘ³ when considering a maximum distance between scales δ_j. The exact number of coefficients are given in Table 1.

The power spectrum of the field is mainly characterised by the $S_2^{{\lambda }_1}$ coefficients defined in Eq. (13). These terms correspond to the average of the power spectrum over the ℓ-wavelet bandpasses (Cheng et al. 2024b). However, they only constrain the power spectrum over a small number of bands and this is usually not precise enough for modelling purpose. Increasing the number of scales j that we probe can be done by decreasing the wavelet dilation parameter η. However, this leads to a large increase of the total number of SC coefficients. For this reason, we considered additional $S_2^{{\lambda }_1}$ coefficients, built with a second filter set with η′ <η (η = 2 and η′ ≃ 1.58) and a single orientation N′ = 1 (isotropic filters). These coefficients are called $S_2^{{\lambda }_1^{\prime }}$ and they allow us to constrain the power spectrum over thinner ℓ-bins.

For physics fields the power spectra can typically be modelled by a power law, at least over certain scales (Cheng et al. 2024b). This leads to SC coefficients, which can vary over several orders of magnitude, since their amplitude is controlled by the I ★ Ψ^λ terms, which filters the initial I field over the j frequency band of the wavelet. This amplitude discrepancy can lead to ill-conditioned optimisations. Following previous works (see, for instance Cheng et al. 2024b), we avoided this issue by normalising the SC statistics from the S₂ coefficients of a reference field that we note S _2,ref. We thus define $${\overline{S}}_1^{{\lambda }_1}=\frac{S_1^{{\lambda }_1}}{\sqrt{S_{2,\text{ref}}^{{\lambda }_1}}},{\overline{S}}_2^{{\lambda }_1}=\frac{S_2^{{\lambda }_1}}{S_{2,\text{ref}}^{{\lambda }_1}},{\overline{S}}_2^{{\lambda }_1^{\prime }}=\frac{S_2^{{\lambda }_1^{\prime }}}{S_{2,\text{ref}}^{{\lambda }_1^{\prime }}},$$(16) $${\overline{S}}_3^{{\lambda }_1,{\lambda }_2}=\frac{S_3^{{\lambda }_1,{\lambda }_2}}{\sqrt{S_{2,\text{ref}}^{{\lambda }_1}{S}_{2,\text{ref}}^{{\lambda }_2}}},{\overline{S}}_4^{{\lambda }_1,{\lambda }_2,{\lambda }_3}=\frac{S_4^{{\lambda }_1,{\lambda }_2,{\lambda }_3}}{\sqrt{S_{2,\text{ref}}^{{\lambda }_2}{S}_{2,\text{ref}}^{{\lambda }_3}}}.$$(17)

When constructing generative models below, we will take the target field as the reference field, which will allows us to deal with coefficients whose values will be at most of order unity.

3 Generative modelling

In this section, we describe how to build generative models from the SC statistics of a given field. We also give some details on the numerical implementation of the generative models and the associated computational cost.

3.1 Maximum entropy generative model

We built generative models under SC constraints. These are maximum entropy micro-canonical models, which are approximately sampled by gradient descent. We refer the readers to Bruna & Mallat (2019) for more details.

These models were constructed from statistics Φ estimated from a target field x_t; in this paper, the target field is a single full-sky map. The associated micro-canonical set Ω_ε of width ε is $${\mathbf{\Omega }}_{\epsilon }=\left\{x:\Vert \mathbf{\Phi }(x)-\mathbf{\Phi }(x_t){\Vert }^2<\epsilon \right\},$$(18)

where ||.|| is the Euclidean norm. The micro-canonical maximum entropy model is the model of maximal entropy defined over Ω_ε, which has a uniform distribution over this set.

In this paper, we approximated this sampling with a gradient descent approach, which consisted in transporting a higher entropy white Gaussian distribution into a distribution supported in Ω_ε. In practice, each new sample was obtained by first drawing a white noise realisation, and then performing a gradient descent in pixel space, or in harmonic space, using a loss function $$\mathbf{ℒ}(x)=\Vert \mathbf{\Phi }(x)-\mathbf{\Phi }(x_t){\Vert }^2.$$(19)

The typical width ε of the micro-canonical ensemble was then fixed by the number of iterations used in the gradient descent. The numerical details of this implementation are given in Sect. 3.2.

In our case, the summary statistics Φ(x) that we considered are the mean over pixels 〈x〉, its variance Var(x), and the normalised SC statistics, defined in Sect. 2.3. Thus, we have $$\mathbf{\Phi }(x)=\left\{\langle x\rangle ,\mathrm{Var}(x),{\overline{S}}_1^{{\lambda }_1},{\overline{S}}_2^{{\lambda }_1},{\overline{S}}_2^{{\lambda }_{1^{\prime }}},{\overline{S}}_3^{{\lambda }_1,{\lambda }_2},{\overline{S}}_4^{{\lambda }_1,{\lambda }_2,{\lambda }_3}\right\}.$$(20)

We note that the target statistics Φ(x_t) were evaluated from a single full-sky image, and that the SC generative models were then built from this single set of constraints. In this respect, our approach differs from machine learning-based approaches, which generally require training on a large, and potentially very expensive, dataset.

3.2 Details on the numerical implementation

In this work, we considered the SC defined on the sphere, constructed using spherical wavelet transforms which in turn depend on efficient spherical harmonic and Wigner transforms (see Sect. 2.2). As outlined in Sect. 3.1, new images were drawn from a SC micro-canonical model by minimising the loss defined in Eq. (19). A plethora of algorithms have been developed to solve such optimisation problems; however, they typically require gradient information, which in turn requires that each component of the loss be differentiable. Consequently, we required the spherical SC and thus the spherical wavelet, spherical harmonic, and Wigner transforms all to be differentiable. Recently, open-source JAX software that is differentiable and graphic processing unit (GPU)-accelerated has been developed for all of these transforms, including s2fft³ for spherical harmonic and Wigner transforms (Price & McEwen 2024) and s2wav⁴ for spherical wavelet transforms (Price et al. 2024). As part of the current work, we have developed a new open-source software implementing the spherical SC called s2scat⁵, which builds on top of s2fft and s2wav.

For a given target field, we computed a generated field by minimising the loss defined in Eq. (19) through a gradient descent in harmonic space with different initial conditions. These initial conditions are Gaussian white noises sampled in the spherical harmonic domain, that is, all of their I_ℓm real and imaginary parts were drawn from the same Gaussian distribution such that the total variance of the target field was reproduced. In this way, the starting angular power spectrum, as defined in Eq. (21), was flat.

To avoid a repeated spherical harmonic transform as the first step at each iteration in the computations, we chose to perform the gradient descent in the spherical harmonic domain rather than in pixel space. The variables we iterated on during the loss minimisation were thus the I_ℓm coefficients. Because the maps are real and thanks to relation (3), we could iterate on the I_ℓm with m ≥ 0 only. The loss minimisation was done through a gradient descent with the L-BFGS algorithm described in Byrd et al. (1995), using the JAX auto-differentiable framework (Bradbury et al. 2018) and the jaxopt package (Blondel et al. 2022). We stopped the optimisation after ~400 iterations, which in our experiment was the typical time for the loss function to decrease by about four orders of magnitude and to reach a plateau at values around 0.1 (meaning, since the loss was not normalised by the number of coefficients, that coefficients are on average constrained at sub-percent accuracy).

3.3 Computational cost

As outlined in Sect. 2, computation of the SC statistics requires repeated spherical convolutions with subsequent non-linear activation functions, in this case modulus operators. Although directional spherical convolutions can be naively computed in pixel space with complexity 𝓞(L⁵) (McEwen et al. 2007), they are more efficiently evaluated in harmonic space with complexity 𝓞(NL³) (McEwen et al. 2007, 2013, 2015b). Furthermore, excellent accuracy can be achieved by computing convolutions in harmonic space since pixelisation artefacts are avoided.

We must repeatedly map to and from spherical harmonic space within our generative model using s2fft (Price et al. 2023). Two operating modes are provided by s2scat: one computes and caches the reduced Wigner d-functions, which are then used at runtime (pre-compute mode); and the other computes these functions on the fly through recursive algorithms (on- the-fly mode). Conceptually, the pre-compute mode is fast but requires 𝓞(L³) memory, whereas the on-the-fly mode is slower but requires at most 𝓞(L²) memory. When running on GPUs for harmonic bandlimits L ≤ 512, we recommend that one adopt the pre-compute mode, deferring to the on-the-fly algorithms at higher resolutions. Although with GPU memory increasing rapidly with hardware developments, it is likely that the precompute mode will be able to be run at higher bandlimits on the latest and upcoming GPUs.

High-level benchmarking results are presented in Table 2. In each case, we consider an azimuthal bandlimit of N = 3, which corresponds to five directions on the sphere, and the full set of anisotropic SC. Our benchmarking was performed on a single NVIDIA A100 GPU with 40GB of on-board memory, although in practice s2scat can be distributed across a large number of GPUs. For completeness, we recorded the time for both a forward and gradient evaluation, in addition to the time required for just-in-time (JIT) compilation. One can also utilise the jax.vmap API, which allows one to batch calls to the maximum entropy model presented in Sect. 3, resulting in more optimal GPU utilisation. For example, suppose we sample from our micro-canonical model with 100 iterations of a first-order optimiser (e.g. ADAM; Kingma & Ba 2014). Generating a single new image at L = 256 takes ~4s, whereas a batched call to generate 20 such images takes ~12s which is ~0.5s per new sample. Furthermore, the jax.pmap API allows one to batch calls across GPU devices, therefore accelerating generation linearly with the number of available GPUs.

4 Validation of the generative models

In this section, we constructed SC generative models from four astrophysical fields showing different types of structures. We then compared the generated fields with the target one. We first did a visual comparison of the maps to assess the quality of the spatial texture reproduction. We then computed various statistics in order to quantitatively evaluate our generative model. For each type of data, we drew 50 samples of the micro-canonical ensemble and computed the mean and the standard deviation over those 50 realisations. For the LSS and the tSZ fields, while the optimisation was performed on the logarithm of the maps, we compared the statistics on the raw image, after taking the inverse transform.

Before any comparison, all maps were filtered in harmonic space, keeping only the I_ℓm coefficients such that ℓ_min ≤ ℓ ≤ L with ${\ell }_{\min }={\eta }^{J_{\min }}$ the central frequency of the wavelet associated with J_min . In this way, we considered only the scales that were constrained during the optimisation. We note however that, if necessary, it is possible to constrain in a similar scheme scales up to ℓ_min = 0.

We also propose a comparison with samples from a Gaussian model built from the power spectrum of the target field. We produced 50 Gaussian realisations using the synfast method from the Healpix package (Górski et al. 2005), which can be used to construct such a Gaussian model from the angular power spectrum of a target field. This allowed us to quantify the contribution of our models built from SC statistics compared to purely Gaussian statistics.

In Appendix A, we directly compared the SC statistics of the target and the generated fields. This was an additional check for the quality of the generative model. Indeed, we expected them to match as they are part of the coefficients constrained during the optimisation.

4.1 Description of the set of maps

We first present the astrophysical and cosmological fields from which we constructed SC generative models. They were expected to have homogeneous statistical properties on the sphere. This property is essential since this is assumed when computing statistics through spatial averages. This requirement, as well as a possible way to avoid this constraint, is discussed in Sect. 4.4. The four different fields, denoted as follows, are:

• LSS, a LSS simulation of weak lensing, from the CosmoGrid dataset (Kacprzak et al. 2023; Fluri et al. 2022);

• tSZ, a tSZ effect simulated map from Simons observatory Galactic foreground simulations (Ade et al. 2019), which were produced using the Sehgal et al. (2010) model with modifications to better match the recent measurements;

• Venus, a map of the Venus planet from Science On a Sphere database⁶;

• CMB, a CMB temperature map produced using the Python Sky Model software (PySM) from Thorne et al. (2017).

We refer the readers to the references given above for more details on these fields. As a Gaussian field, the CMB map is a good null test for our method, as presented in Appendix B. The other fields all originate from non-linear physical processes, and thus have highly non-Gaussian structures, as can be seen in the left column of Fig. 3. The diversity of these fields illustrates the generality and the versatility of our method, which could be used for various physical datasets.

While all of these simulated maps are available in Healpix format, the computation of the SC was done directly from the harmonic space. The conversion to this space was done using L − 1 = ℓ_max = 2nside, and thus acts as a low-pass filter operation. This implies that spatial frequencies at ℓ ≥ L are filtered out in this operation. We note that calculating the SC and performing the optimisation directly in spherical harmonic space means that there is no particular constraint on the sampling on the target map, even if the internal SC calculation steps are based here on MW or alternative sampling schemes to improve accuracy.

During the optimisation process, all maps were normalised such that their mean is zero and their standard deviation is one. In addition, the LSS and tSZ fields are highly non-Gaussian, making them difficult to model directly even with SC statistics. This is why we instead chose to model the logarithm of these maps⁷. This logarithmic transform brings the distribution closer to a Gaussian one, and reduces in particular the weight of the high amplitude tail of the probability density function (PDF), allowing for better SC generative models. At the end of the optimisation, however, we took the inverse transform for these maps, and we assessed the quality of the generative model on the raw images.

The generative models were run on MW maps with a resolution of L = 256, which corresponds to L(2L − 1) = 130 816 pixels. These real signals have L² = 65 536 complex harmonic coefficients. For the directional wavelets used to build the SC, we considered a dyadic scaling η = 2 and N = 3. This led to J_max = 8 and 2N − 1 = 5 orientations. As shown in Table 1, the value of J_min was tuned to each field in order to take into account the largest spatial scales that compose the maps. The number J of scales that we probed is also given in the Table. The maximum distance δ_j between scales was fixed to five, which allows us to divide the total number of SC coefficients by approximately two without degrading the quality of the generative model. Concerning the additional $S_2^{{\lambda }^{\prime }}$^′ coefficients, we chose an axisymmetric filter set with N′ = 1 and a scaling given by η′ ≃ 1.58. The exact total number of terms making up the summary statistics Φ(x) is given in Table 1.

Fig. 3 Visual validation of the generative model. From top to bottom, we show the maps for the LSS, tSZ, and Venus fields. The left column is the original target field. The central column shows one sample of the generated maps. The right column shows a Gaussian field with the same power spectrum as the target. For LSS and tSZ, we plotted the logarithm of the fields in order to better see the texture details. The colour bars are identical within each field.

4.2 Visual validation

As a first test, we can visually compare the target field and the generated ones, as presented in Fig. 3 for the LSS, tSZ, and Venus fields. They appear to be visually very similar to the original maps, which clearly shows that the SC statistics capture an important part of the non-Gaussian texture of the field. On the contrary, the structures are not reproduced in the Gaussian realisations shown on the right column. For LSS and tSZ, we plot the logarithm of the fields, which allows us to better see the textures. In addition, in Fig. 4, we show a zoom on a smaller region to better visualise the details in the spatial structures. In Appendix C, we also show four realisations of the fields starting from different initial conditions. This shows the ability of our generative models to sample independent realisations while capturing the overall texture of the fields.

4.3 Statistical validation with standard summary statistics

Following a similar approach to Cheng et al. (2024b) and Price et al. (2023), we compared summary statistics between the target field and the generated field. As previously, we show the mean and the standard deviation computed over 50 realisations. The summary statistics we chose to compare are:

• the PDF of the map;

• the angular power spectrum;

• the three Minkowski functionals.

The PDF of the maps and the three Minkowski functionals were performed on Healpix maps. This was done by projecting the output I_ℓm from the loss minimisation onto the Healpix map by an inverse spherical harmonic transform at the end of the generative process.

4.3.1 Probability density function

The PDF for the LSS, tSZ, and Venus fields are shown in Fig. 5, computed on the Healpix maps. On the first row, we show the PDF with a linear y-axis scaling, while on the second row, we show them with a logarithmic y-axis scaling in order to better exhibit the tails of the distributions on several orders of magnitude. The target fields are shown in blue and the generated ones in red. In yellow, we also show the comparison with the Gaussian realisations. By definition, the PDF of these realisations presents a Gaussian profile in linear scale and a parabolic profile in logarithmic scale.

While the Venus field has a PDF that only slightly differs from the Gaussian case, the two other fields clearly have nonGaussian features with large tails. The comparison of the target and generated PDFs with the Gaussian PDF also allows us to better see their non-symmetric shape, which is characteristic of non-Gaussian features. As we can see, the PDFs for SC models are well reproduced on at least three orders of magnitude. The results obtained for the LSS fields, which are very good up to five orders of magnitude, are especially striking. On the other hand, the results for the tSZ field begin to push the expressive limit of our generative models. For the Venus maps, we have identified the abrupt jump in the histogram as a flaw in the data used, which does not particularly illustrate a limitation of our maximum entropy SC model.

Fig. 4 Zoom-in on texture details. We show a zoom-in on a region for the LSS (top), tSZ (middle), and Venus (bottom) fields. Similarly to the previous figure, for LSS and tSZ, we plotted the logarithm of the fields in order to better see the texture details. The colour bars are identical within each field.

4.3.2 Angular power spectrum

We calculated the power spectrum in the usual way as $$C_{\ell }=\frac{1}{2\ell +1}{\displaystyle \sum _{m=-\ell }^{m=\ell }|}{I}_{\ell m}{|}^2,$$(21)

where the normalisation factor 1/(2ℓ + 1) yields a flat power spectrum in case of white Gaussian noise.

As discussed in Sect. 2.3, the power spectrum was constrained during the optimisation through the S₂ and the additional $S_2^{\prime }$ coefficients. We note however that these coefficients did not constrain the full power spectrum, as each term constrains only a weighted power spectrum over the frequency support of the associated wavelet.

The third row of Fig. 5 shows the results for the generation of the LSS, tSZ, and Venus fields, from left to right. Power spectra are well reproduced over all scales, even when they vary over up to four orders of magnitude. However, small oscillations around the target can be seen in the generated power spectra. These are residual features related to the frequency bands of the wavelets, which illustrate the trade-off between the quality of reproduction we want to achieve and the number of filters we use, that is to say, the computational efficiency of our generative model.

We note that in this paper, we include 11 of these additional $S_2^{\prime }$ coefficients. This number, and the precise shape of the wavelets used, could however be tuned to better reproduce the power spectrum of the target. However, care must be taken not to over-constrain these terms, as all of the samples generated would then have a power spectrum very close to that of the target, which does not necessarily correspond to a good generative model. The introduction of $S_2^{\prime }$ terms is an improvement with respect to previous work, since it allows us to have better power spectrum constraints without significantly increasing the overall number of SC statistics.

4.3.3 Minkowski functionals

Finally, we computed the Minkowski functionals on the Healpix maps. These standard non-Gaussian statistics characterise the topology of the level sets of the field. In two dimensions, there are three Minkowski functionals, V₀(u), V₁(u), and V₂(u), that depend on a pixel value threshold u. We computed them using Pynkowski software (Carones et al. 2024). We refer the reader to this publication for the complete definition of those statistics. The result is shown in Fig. 6 only for the LSS field, while the others are shown in Appendix D. In yellow, as a comparison, we show the case of the Gaussian realisations. For the generated fields and the Gaussian realisations, we plotted the mean (solid) and the standard deviation (shadow envelope) computed over 50 realisations. Thus, the SC models encompass very well these non-Gaussian statistical features.

4.4 Limitations and discussions

This work is a first implementation of generative models from state-of-the-art spherical SCs. As a first proof of concept, we constructed and validated these models on various cosmological fields, most of them simulated. However, building such models from real data, or using these tools to perform statistical component separation, may require some additional work to deal with their own specificity. In this section, we comment on some of these limitations and how to overcome them.

A first limitation of our current work is the fact that our models assume the statistical homogeneity of the fields studied. However, the ability to deal with non-homogeneous physical processes is usually required when modelling astrophysical fields, such as the Galactic emissions, whose properties typically vary strongly with latitude on the sky. An efficient way to deal with this issue is to rely on different masks in pixel space, for which statistical properties can be constrained independently (Delouis et al. 2022). However, this requires a trade-off between the size and number of masks: using a larger number of masks gives a better description of large-scale variations in statistical properties but increases the variance of the SC statistics estimates on each mask due to the smaller number of pixels used, in addition to increasing the total number of SC coefficients and the computational and memory cost.

A second limitation is the map resolution that we can achieve. For now, the generation of a new map at L = 256 and N = 3 takes ≤ 1s on a single GPU. This is thanks in part to a large number of pre-computed matrices necessary for the Wigner transform, which are stored in memory (several Gigabytes). Specifically this memory is cubic with L, which is prohibitive at high L. When increasing the resolution beyond L ∼ 1024, we usually reach the GPU memory limit and the coefficients have to be computed on the fly, although this of course depends on available GPU specifications. Critically, the on-the-fly approach dramatically reduces memory requirements so that generations at high L are at least feasible but at the cost of a significantly increased computation time. In future work, we will explore further optimisations for high L. A key avenue we will explore is the introduction of hybrid wavelet convolutions, which operate efficiently in pixel and harmonic space at high and low resolutions, respectively (see e.g. Delouis et al. 2022; Ocampo et al. 2023).

Fig. 5 Statistical validation. The PDF and angular power spectra for the LSS, tSZ, and Venus fields (left to right). The first row shows the PDF with a linear y-axis scaling, while the second row shows the same PDF with a logarithmic y-axis. The third row shows the angular power spectra. The target is shown in blue, the generated fields in red, and the Gaussian realisations in yellow. We plotted the mean (solid line) and the standard deviation (shadow envelope) over 50 realisations.

Fig. 6 Minkowski functionals. We plotted the three Minkowski functionals V0, V1 , and V2 for the LSS field. Blue is the target, red the generated fields, and yellow the Gaussian fields. For the generated fields and the Gaussian realisations, we plotted the mean (solid) and the standard deviation (shadow envelope) computed over 50 realisations.

5 Conclusions

The main result of this paper is the extension of state-of-the-art scattering transforms to spherical fields. We have worked with the last generation of scattering transform statistics, named SCs (Cheng et al. 2024b), which were previously introduced for 1D and 2D planar fields. They have the advantage of relying only on successive wavelet transforms and modulus, as well as on covariances, and do not require any translations. We have also used state-of-the-art directional wavelets on the sphere, computed in spherical harmonic space (McEwen et al. 2015b). The numerical implementation of this work, s2scat, is open-source and publicly available. Furthermore, it is fully auto-differentiable, using the JAX Python framework (Bradbury et al. 2018) and building on the s2fft/ s2wav packages (Price & McEwen 2024; Price et al. 2024).

These developments allow us to build generative models of full sky spherical fields without the need for large training datasets. In fact, our method holds in the limit case of a single data realisation. The performance of those generative models was validated quantitatively on different fields: a LSS weak lensing field, maps of the tSZ effect and of the CMB, and a map of the Venus surface, for which they performed extremely well. The diversity in terms of structures between the maps shows the impressive ability of SCs to comprehensively characterise very different non-Gaussian textures.

This work introduces a new powerful innovative approach for spherical data, and it opens interesting perspectives for astrophysical applications. In particular, we plan to use it for the study and the modelling of CMB astrophysical foregrounds. The first goal will be to have a tool to produce multiple realisations of the different astrophysical components, for example, using the AGORA simulations (Omori 2024). Then, scattering transforms could play a role in component separation, relying on both recently developed scattering transform-based statistical component separation approaches, as well as investigating how classical component separation methods could benefit from scattering transforms, using the non-Gaussianities as an additional lever arm for disentangling different components.

Finally, we also point out that SC statistics provide highly informative sets of statistics, which could be very useful in tasks such as parameter inference, such as simulation-based inference (SBI; Cranmer et al. 2020), from large cosmological surveys (see, for instance Régaldo-Saint Blancard et al. 2024; Gatti et al. 2024; Cheng et al. 2024a). This could be all the more useful as the compression factor they enable, compared with a direct description in pixel space, becomes extremely large at high resolution, due to logarithmic scale binning. This property could be further enhanced by using compression schemes such as the one presented in Cheng et al. (2024b), which can make SCs a very informative and versatile compressed set of statistics.

Data availability

We make our code available to the community so that this work can be easily reproduced and developed further; https://github.com/astro-informatics/s2scat

Acknowledgements

We thank Sixin Zhang and Anthony Banday for very helpful discussions all along the project. Most of the simulations were produced on the MesoPSL calculation centre, for which we thank the administrators. This work was supported in part by a CNES postdoctoral fellowship. MAP and JDM are supported in part by EPSRC (grant number EP/W007673/1) and STFC (grant number ST/W001136/1).

Appendix A SC statistics

In Figure A.1 we show the normalised SC coefficients ${\overline{S}}_1,{\overline{S}}_2,\text{\hspace{0.17em} and \hspace{0.17em}}{\overline{S}}_3$ of the LSS field. We plot the coefficients of the target field in blue, the generated ones in red, as well as the Gaussian realisations in yellow. Coefficients are plot following the lexicographic order. We chose not to show ${\overline{S}}_4$ for readability because of the large number of coefficients.

Regarding the Gaussian realisations, shown in yellow, we expect the ${\overline{S}}_3$ coefficients to be equal to zero up to the correlations induced by the overlapping between wavelet bands. As we can see, for ${\overline{S}}_3$, the mean is centred on zero.

By construction, ${\overline{S}}_2$ for the target field is equal to one. This is because we have considered the S₂ coefficients of the target as the reference to normalise all the coefficients, as described in Sect. 2.3. In this way, all the normalised coefficients are of the order of the unit. As we can see, SC statistics are well constrained by the optimisation, the generated coefficients in red well overlap the target coefficients in blue. ${\overline{S}}_3$ coefficients strongly differ from the Gaussian field, showing clear non-Gaussian signatures.

Fig. A.1 Normalised SC coefficients. We plot ${\overline{S}}_1,{\overline{S}}_2,\text{\hspace{0.17em} and \hspace{0.17em}}{\overline{S}}_3$ for the logarithm of the LSS field. We show the coefficients from the target field (blue), the generated fields (red), and the equivalent Gaussian realisations (yellow). The mean and the standard deviation over 50 realisations are shown as a solid line with a shadow envelope.

Appendix B CMB map as a null test

It is important to check that the generative model behaves as we expect for a Gaussian field. This is an important validation for a maximum entropy generative model. This is why we tested it on the CMB map. The result is shown in Fig. B.1. The upper parts shows the target map, a generated field and a Gaussian realisation. As expected, the three maps look similar. We also plot the PDF and the power spectrum which match very well.

Fig. B.1 Generative model for a CMB map. The upper part shows (left to right): the target map, a generated field, and a Gaussian realisation. The second row shows the PDF (linear and logarithmic scales) and the angular power spectrum.

Appendix C Multiple realisations

Figure C.1 shows multiple realisations obtained from the generative model, changing the initial Gaussian random noise. For each field we show four maps out of the 50 that we computed. This is to illustrate the visual similarity between the realisations.

Fig. C.1 Multiple realisations. Four generative models of the LSS, tSZ, and Venus fields (from top to bottom), obtained by changing the initial Gaussian random noise. In total, we ran 50 realisations for each field. The colour scales are identical within each field.

Appendix D Minkowski functionals for the three others fields

Figure D.1 shows the Minkowski functionals for the tSZ, Venus, and CMB maps.

Fig. D.1 Additional Minkowski functionals. The three Minkowski functionals V0, V1, and V2 for the tSZ (upper left), Venus (upper right), and CMB (bottom) fields. Blue is the target, red the generated fields, and yellow the Gaussian fields. For the generated fields and the Gaussian realisations, we plot the mean (solid) and the standard deviation (shadow envelope) computed over 50 realisations.

References

All Tables

All Figures

	Fig. 1 Harmonic representation of the wavelet filters. Left: real part of the ${\Psi }_{\ell m}^{j=4}$ filter. Right: cut at m = 0 of the full filter set for j spanning from Jmin = 1 to Jmax = 8.
In the text

	Fig. 2 Spatial representation of the wavelet filters. The circles contain orthographic projections of the directional spherical wavelets for two scales and three angles, viewed looking down from the North pole. In this case, we set N = 3 so that the angle γ takes 2N − 1 = 5 values between 0 and 360 deg, but only three are shown here.
In the text

Fig. 3 Visual validation of the generative model. From top to bottom, we show the maps for the LSS, tSZ, and Venus fields. The left column is the original target field. The central column shows one sample of the generated maps. The right column shows a Gaussian field with the same power spectrum as the target. For LSS and tSZ, we plotted the logarithm of the fields in order to better see the texture details. The colour bars are identical within each field.

In the text

	Fig. 4 Zoom-in on texture details. We show a zoom-in on a region for the LSS (top), tSZ (middle), and Venus (bottom) fields. Similarly to the previous figure, for LSS and tSZ, we plotted the logarithm of the fields in order to better see the texture details. The colour bars are identical within each field.
In the text

Fig. 5 Statistical validation. The PDF and angular power spectra for the LSS, tSZ, and Venus fields (left to right). The first row shows the PDF with a linear y-axis scaling, while the second row shows the same PDF with a logarithmic y-axis. The third row shows the angular power spectra. The target is shown in blue, the generated fields in red, and the Gaussian realisations in yellow. We plotted the mean (solid line) and the standard deviation (shadow envelope) over 50 realisations.

In the text

	Fig. 6 Minkowski functionals. We plotted the three Minkowski functionals V0, V1 , and V2 for the LSS field. Blue is the target, red the generated fields, and yellow the Gaussian fields. For the generated fields and the Gaussian realisations, we plotted the mean (solid) and the standard deviation (shadow envelope) computed over 50 realisations.
In the text

Fig. A.1 Normalised SC coefficients. We plot ${\overline{S}}_1,{\overline{S}}_2,\text{\hspace{0.17em} and \hspace{0.17em}}{\overline{S}}_3$ for the logarithm of the LSS field. We show the coefficients from the target field (blue), the generated fields (red), and the equivalent Gaussian realisations (yellow). The mean and the standard deviation over 50 realisations are shown as a solid line with a shadow envelope.

In the text

	Fig. B.1 Generative model for a CMB map. The upper part shows (left to right): the target map, a generated field, and a Gaussian realisation. The second row shows the PDF (linear and logarithmic scales) and the angular power spectrum.
In the text

	Fig. C.1 Multiple realisations. Four generative models of the LSS, tSZ, and Venus fields (from top to bottom), obtained by changing the initial Gaussian random noise. In total, we ran 50 realisations for each field. The colour scales are identical within each field.
In the text

	Fig. D.1 Additional Minkowski functionals. The three Minkowski functionals V0, V1, and V2 for the tSZ (upper left), Venus (upper right), and CMB (bottom) fields. Blue is the target, red the generated fields, and yellow the Gaussian fields. For the generated fields and the Gaussian realisations, we plot the mean (solid) and the standard deviation (shadow envelope) computed over 50 realisations.
In the text